- Research article
- Open Access
Novel Plasmodium falciparum metabolic network reconstruction identifies shifts associated with clinical antimalarial resistance
BMC Genomicsvolume 18, Article number: 543 (2017)
Malaria remains a major public health burden and resistance has emerged to every antimalarial on the market, including the frontline drug, artemisinin. Our limited understanding of Plasmodium biology hinders the elucidation of resistance mechanisms. In this regard, systems biology approaches can facilitate the integration of existing experimental knowledge and further understanding of these mechanisms.
Here, we developed a novel genome-scale metabolic network reconstruction, iPfal17, of the asexual blood-stage P. falciparum parasite to expand our understanding of metabolic changes that support resistance. We identified 11 metabolic tasks to evaluate iPfal17 performance. Flux balance analysis and simulation of gene knockouts and enzyme inhibition predict candidate drug targets unique to resistant parasites. Moreover, integration of clinical parasite transcriptomes into the iPfal17 reconstruction reveals patterns associated with antimalarial resistance. These results predict that artemisinin sensitive and resistant parasites differentially utilize scavenging and biosynthetic pathways for multiple essential metabolites, including folate and polyamines. Our findings are consistent with experimental literature, while generating novel hypotheses about artemisinin resistance and parasite biology. We detect evidence that resistant parasites maintain greater metabolic flexibility, perhaps representing an incomplete transition to the metabolic state most appropriate for nutrient-rich blood.
Using this systems biology approach, we identify metabolic shifts that arise with or in support of the resistant phenotype. This perspective allows us to more productively analyze and interpret clinical expression data for the identification of candidate drug targets for the treatment of resistant parasites.
Three billion people are at risk for malaria infection globally and treatment approaches are failing. Malaria is caused by Plasmodium parasites, and most deaths are associated with human-infective P. falciparum. Without an efficacious vaccine, antimalarials are essential to combat the severity and spread of disease. Combination therapies are implemented to preserve antimalarial efficacy and slow resistance development [1,2,3]; despite this approach, this eukaryotic pathogen has developed resistance to every antimalarial on the market [4,5,6].
Typically, resistance is conferred by genomic changes that lead to drug export or impaired drug binding (for example ); however, non-genetic mechanisms have also been implicated in Plasmodium resistance development [8,9,10] and other pathogenic organisms, such as Pseudomonas aeruginosa  (reviewed in ). These laboratory-based studies provide insight into metabolic flexibility but the presence of relatively few examples limit our understanding of this method of adaptation, especially in malaria. Here, we aim to look beyond genetic mechanisms of resistance to identify resistance-associated metabolic adaptation. We hypothesize that metabolic changes must occur to support the resistance phenotype and resistance-conferring mutations. Ultimately, these changes, or ‘shifts,’ are required to increase the fitness of resistant parasites, or support the development of additional genetic changes that affect fitness. Metabolic or phenotypic ‘background’ could be as important as genetic background in the development of resistance.
In clinical malaria infections, artemisinin resistance is established in Southeast Asia [13,14,15]. This phenotype is correlated with mutations in the P. falciparum Kelch13 gene [13, 14, 16, 17] and changes in both signalling pathways [18,19,20,21] and organellar function [22,23,24,25,26,27,28,29]. Overall, due to the complexity of artemisinin’s mechanism of killing (see citations above and [30,31,32,33,34,35,36]), it has been challenging to separate the causes and effects of resistance. For this reason, there are few novel solutions to antimalarial resistance beyond altering the components of combination therapies to regain efficacy (e.g. artemisinin-atovaquone-proguanil ). We aim to gain a new perspective on resistance by viewing it through a ‘metabolic lens’. By characterizing the metabolic shifts that occur during or after resistance acquisition, we can begin to understand more about what it takes to support new functions, such as novel signalling (e.g. PI3K signalling is affected by PfKelch13 mutations [14, 15, 20, 37, 38]), drug detoxification (e.g. regulating ROS stress associated with artemisinin treatment [24, 25, 30, 33]), or stage alterations (e.g. dormancy of early ring stages [18, 39,40,41,42]) in resistant parasites. Once we identify these compensatory changes, we can potentially target them. Plasmodium metabolic genes are better characterized than signalling pathways, as (for example) PlasmoDB identifies 43 3D7 genes associated with the term ‘signalling’ as opposed to 1112 3D7 genes associated with the term ‘metabolism’ , and many antimalarials target metabolic functions [44,45,46,47]. Moreover, metabolism has been described as the best-understood cellular process , making interpreting metabolic analyses more tractable. Ultimately, if we can identify targetable conserved metabolic differences that arise with or in support of resistance, we can develop more robust antimalarial combination therapies aimed at preventing resistance.
Here, we use a systems biology approach to analyze the metabolic profile associated with resistant and sensitive parasites. First, to maximize the accuracy of our predictions, we curated an existing genome-scale network reconstruction of asexual blood-stage P. falciparum metabolism. Using constraint-based metabolic modeling, we integrated transcriptomic data from over 300 clinical isolates from Cambodia and Vietnam with varying levels of artemisinin sensitivity. This approach identified innate metabolic differences that arise with or in support of the resistant phenotype, despite large clinical variability, over multiple genetic backgrounds. Additionally, we were able to explore the functional consequences of expression changes by predicting essential enzymes within these distinct metabolic contexts; these enzymes are candidate drug targets for the prevention of drug resistance.
Analysis of artemisinin sensitive and resistant transcriptomes
In order to investigate the presence of a distinct metabolic phenotype in artemisinin resistant parasites, we analysed a previously published expression dataset of clinical isolates from Southeast Asia (NCBI Gene Expression Omnibus accession: GSE59097). Patient blood samples were collected immediately prior to beginning artemisinin combination therapy, and their relative expression was evaluated via microarray . This dataset profiles (1) in vivo artemisinin naïve parasites, providing a view of the innate differences between sensitive and resistance parasites, and (2) a diverse population of parasites collected from multiple collection sites across two countries, allowing us to summarize variable resistant phenotypes that laboratory adapted parasites and in vitro assays cannot practically encompass.
We confined our analysis of this previously published expression data to ring-stage parasites from Cambodia and Vietnam, two countries that had clear resistant and sensitive parasite populations as defined by parasite clearance half-life, an in vivo phenotypic measure of resistance, and PfKelch13 mutations, a commonly-used genetic marker of resistance (Fig. 1a & b). There were 97 and 24 ring-stage resistant parasite expression profiles from Cambodia and Vietnam, respectively; resistant parasites are defined by both the presence of PfKelch13 mutations and a parasite clearance half-life of more than 5 h. There were 141 and 43 ring-stage sensitive parasite expression profiles from Cambodia and Vietnam, respectively, as defined by wild-type PfKelch13 alleles and clearance half-life of less than 5 h. Despite obvious genotypic and phenotypic separation (Fig. 1a & b), artemisinin sensitive and resistant parasites do not separate well by hierarchical clustering of expression data (Fig. 1c).
Additionally, when comparing sensitive parasites to resistant parasites in either country, the fold change of transcript expression is moderate; no genes exhibited notable differential expression across both analyses (fold change >2 or <0.5 for both Cambodia and Vietnam sample sets, data not shown). Among metabolic genes specifically, expression differences are small (maximum fold change 0.6 and 1.6) and few are both significant and conserved between data sets (11 in common from 174 in Cambodia and 37 in Vietnam; Additional file 1: Figure S1A & B). Large amounts of transcriptional variation (due to stage-dependent expression, genotypic variability, and host-pathogen interactions) across the population of clinical parasites may hide differences in the data sets. Moreover, we built a Random Forest classifier with expression data to predict resistance outcomes; the classifier predicted resistance poorly, with only 30.77% sensitivity (indicating only 30.77% of resistant samples were correctly identified) and 97.96% specificity (indicating 97.96% of sensitive samples were correctly identified) (Additional file 2: Figure S2A).
Although the expression data classifier performed poorly, a similar classifier built from metadata associated with each sample (patient and parasite characteristics) was highly predictive of resistance status with 85.71% sensitivity and 88.91% specificity (Additional file 2: Figure S2B). In our analysis, two specific mutations and collection site were the most predictive of resistance status; removing any of these three variables decreased classifier accuracy by over 20%. If Kelch13 mutations are used to predict resistance (rather than used to define resistance), Kelch13 mutations are most predicative of resistance (data not shown). Thus, metadata better predicts resistance than expression data. In order to deconvolve this innate variability and identify functional cellular changes associated with varying levels of artemisinin sensitivity, we integrated metabolic expression data into a genome-scale metabolic model of blood-stage P. falciparum.
Manual metabolic network curation
To maximize the predictive ability of the metabolic network model, we curated an existing, well-validated reconstruction of asexual blood-stage P. falciparum (iTH366, ) to improve its scope, and species- and stage-specificities. Our curated reconstruction, iPfal17, includes all metabolic reactions encoded by characterized genes in the parasite’s genome, summarizing metabolic behavior during the asexual blood-stage parasite. It is larger in scope from the previously published version due to the addition of 268 reactions (Table 1 , Additional file 3: Table S1 & S2), with 9.6% more enzymatic reactions and 2.3% more reactions with gene annotations. We also added 124 genes to the network (Table 1 & Additional file 3: Table S1). It is larger in scope and gene coverage than a recent de novo reconstruction (Table 1). iPfal17 has gene annotations for 80.0% of enzymatic reactions, and 20.5% of transport and exchange reactions (Fig. 2). iPfal17 includes 25.4% of the 1178 EC annotations in the P. falciparum genome, adding 14 EC numbers  (Additional file 3: Table S1). We evaluated enzyme complex or isozyme status and replaced 7 gene-protein-reaction relationships (Additional file 3: Table S1).
Following curation, the species- and stage-specificity of the model was also improved. Gene annotations were evaluated against PlasmoDB resources , resulting in 124 additional gene annotations (Additional file 3: Table S1). Importantly, we removed cellular import of pyrimidines from the host erythrocyte, as P. falciparum relies on de novo synthesis (Additional file 3: Table S2) [47, 51]. Blood-stage specificity was improved by removing genes only used in other life stages (specifically, the gene encoding chitinase ). Additionally, 77 functionally unnecessary reactions were removed due to a lack of genetic and biochemical support (Additional file 3: Table S2). Reactions necessary for growth were added manually (Additional file 3: Table S1). Reactions were individually curated, changing metabolite utilization and stoichiometry (Additional file 3: Table S1).
The iPfal17 reconstruction contains five compartments: extracellular space and four intracellular compartments (cytoplasmic, mitochondrial, apicoplast, and food vacuole, Additional file 3: Table S1). Few studies since the Plata, et al. reconstruction (iTH366) investigated protein localization and therefore, few changes were made to compartmental assignments; the food vacuole compartment, containing two reactions, was added in this version of the reconstruction (Additional file 3: Table S1). As in iTH366, reactions with unknown localization were placed within the cytoplasm (Additional file 3: Table S1) . Again, similar to iTH366, a mitochondrial inner matrix was not added, as there is no evidence that the blood-stage parasite requires a proton gradient for energy production [51, 54, 55]. Nonpolar metabolites generated in one compartment and utilized in another were transported as needed for network functionality by assuming passive diffusion .
We also included annotations that will accelerate future curation efforts. First, we did not remove blocked reactions (those that do not carry flux due to their lack of connectivity to other components of the network) because further research may add connectivity to these network components. iPfal17 contains 303 blocked reactions and 78 dead-end metabolites (specifically, 32 metabolites are not consumed and 46 are not produced). For example, 4-pyridoxate (a byproduct of vitamin B6 biosynthesis) is included; production is supported by bioinformatic analyses of the parasite genome, but the metabolite function or excretion pathway is not known. Second, citations are included within iPfal17 to identify the date of discovery and degree of literature support for each reaction (Additional file 3: Table S1 & S2). Literature support was only added to modified reactions, resulting in 231 citations (Table 1 & Additional file 3: Table S1).
Metabolomics curation of biomass reaction
For the newly curated iPfal17 model, we modified the Plasmodium biomass reaction to better represent in vitro data (Table 2). We added tRNA-ligated amino acids to the amino acid requirements to force protein production, rather than only demanding free amino acids. Additionally, lipid classes were added based on recently published metabolomics findings; phosphatidylinositol, phosphatidylglycerol, sphingomyelin, diacylglycerides, and triglycerides were added due to their observed increase in abundance between uninfected and infected erythrocytes . Phosphatidylcholine ethers, acyl phosphatidylgycerol, lysophosphatidylinositol, bis(monoacyl-glyceryl)phosphate, and monosialodihexosylganglioside were excluded from the biomass reaction, as there is no known Plasmodium catabolism or import pathways for these lipids . Analysis of metabolomics data enabled further curation of the biomass reaction with the addition of malate, alpha-ketoglutarate, and glutathione (both reduced and oxidized) [28, 57, 58]. Importantly, we included the requirement for cellular export of lactate and hemozoin. Lactate is measured in extracellular in vitro metabolomics and in vivo via blood acidosis; it is the terminal product of glycolysis, the sole energy production pathway used by the blood-stage parasite [59,60,61,62]. By requiring lactate export, we force the model to utilize glycolytic energy metabolism. Similarly, hemoglobin degradation is essential for the blood-stage parasite to produce free amino acids. Parasites can also import and synthesize some amino acids, but the breakdown of hemoglobin (and subsequent production of its byproduct, hemozoin) is necessary for growth [23, 63, 64]. Thus, by requiring hemozoin export, we force the in silico parasite to degrade hemoglobin as the primary pathway for amino acid production.
iPfal17 validation and functional requirements
To validate the model against experimental results, essential metabolic tasks of blood-stage growth were identified and evaluated (Table 3). These tasks simulate experimental manipulations of the parasite or culturing environment, or clinical observations. For example, the parasite is able to grow with glucose as the sole carbon source and hypoxanthine as the purine source, and the parasite’s induction of blood acidosis via lactate [65,66,67]. Additional tasks include the parasite’s failure to grow in the presence of antimetabolites for riboflavin, nicotinamide, thiamine, and pyridoxine . We defined this set of tasks to provide a framework for curation and validation efforts of future network reconstructions. Although iPfal17 fails to pass all metabolic task simulations, we believe this is the most comprehensive and accurate model to date due to the curation efforts and results from tests of the metabolic tasks. Failures generally exist in pathways that currently contain many reversible reactions (i.e. tasks 5a–b for glycolysis) or if the experimental evidence is not mechanistic (i.e. tasks 1a–d) or fully characterized (i.e. task 4; Table 3).
We also evaluated predictions of the effects of gene knockouts and enzyme inhibitors using previously published experimental results (Table 4; Additional file 3: Table S7). Our updated model had improved accuracy of gene and reaction essentiality predictions, compared to previous models (Table 4). We predict that there are 159 essential reactions, and 107 lethal single gene knockouts (Additional file 3: Table S3 & S4). Of experimentally validated knockouts, iPfal17 accurately predicts essentiality of 79.5% of genes and enzymes tested in P. falciparum and 61.4% for those tested in P. berghei (Tables 3 & 4 , Additional file 3: Table S6); predictions are also more accurate for gene knockouts and are less accurate in predicting enzyme inhibition (Tables 3 & 4).
Integration of expression data into the metabolic model
With this high-quality metabolic network reconstruction, we integrated expression data from sensitive and resistant parasites collected in Cambodia and Vietnam into iPfal17 using the Metabolic Adjustment for Differential Expression algorithm (MADE ). MADE constrains gene utilization in the model to maximally account for statistically significant changes in expression while maintaining network functionality requirements (e.g. parasite viability). MADE integrates differential expression by minimizing the difference between significant expression changes (up/down) and model constraints (gene usage); this avoids arbitrary thresholding and ensures the condition-specific model is consistent with experimental data. Essential genes and genes supported by expression data (by having no change in expression or being upregulated in that condition) remain in the model that represents that condition. Conversely, if a gene is significantly down regulated in a condition and not functionally necessary for metabolism, the reactions catalyzed by the encoded enzyme will be removed from the model. Therefore, condition-specific models contain a reduced network with a subset of reactions annotated in the original curated reconstruction that are either necessary for network functionality (as defined by the objective function) and/or are supported by expression data . Thus, MADE generates functional condition-specific models representing the cell’s metabolic capability given the condition-specific expression.
MADE integration of sensitive and resistant expression data from both countries generated four condition-specific models (Fig. 3). By comparing these models, we identified differences in gene and pathway utilization between resistant and sensitive parasites that are consistent between the isolates from the two countries (Additional file 4: Fig. S3). First, we conducted an enrichment analysis on genes that remain in (i.e. can be utilized by) each constrained model by comparing to the unconstrained curated model. As expected, all four models were enriched with genes involved in pathways with many essential reactions or little redundancy, such as transport reactions, tRNA synthesis, purine metabolism, and others (Additional file 5: Fig. S4, see model). Sensitive (wild type) models corresponding to isolates from both Cambodia and Vietnam are uniquely enriched with the utilization of genes involved in the metabolism of nicotinate/nicotinamide (p-value = 1.47*10−2), glutamate (p-value = 1.28*10−13), and selenocysteine (p-value = 5.85*10−4). Thus, sensitive models contain more reactions in these pathways than the unconstrained model, resulting from increased expression of these pathways in sensitive parasites (Additional file 5: Fig. S4). Resistant models from both countries are uniquely enriched with the utilization of genes involved in pyrimidine (p-value = 2.18*10−7), polyamine (p-value = 4.39*10−4), redox reactions (p-value = 5.13*10−5), and central carbon metabolism (glycolysis [p-value = 4.39*10−4] and the pentose phosphate pathway [p-value = 6.06*10−3]). Thus, resistant models have a larger proportion of their total reactions associated with these pathways than the original unconstrained model, whereas sensitive models do not have this enrichment. This indicates that these pathways are upregulated in resistant parasites and may remain important for metabolism in the resistant state (Additional file 5: Fig. S4).
Identification of conserved and uniquely essential pathways
Beyond general differences in pathway utilization, which encompasses both essentiality and pathway-level differences in expression, artemisinin sensitive and resistant parasites have unique essential genes and reactions. To identify these essential reactions and provide insight on targetable metabolic enzymes in the clinical isolates, we performed in silico single gene and reaction deletions with each of the four condition-specific models. Datasets from the parasites from each country were initially analyzed separately and then lists were compared to ensure resistance-associated trends are reproducible and observed in independent analyses. As expected, we identified many essential functions conserved in all models (Additional file 3: Table S5), which is consistent with an active core metabolism required for basic parasite survival. Importantly, 21 reactions were essential in only resistant models, but not sensitive models (Table 5). Theoretically, drugs targeting these reactions would kill resistant parasites and have no effect on sensitive parasites; thus, there would be no selective pressure within the sensitive parasite population to develop resistance to these drugs. This list included serine hydroxymethyltransferase (PFL1720w in folate metabolism), the glycine cleavage system (PFL1550w and others in folate metabolism), thiamine diphosphokinase (PFI1195c in cofactor metabolism, specifically thiamine diphosphate), fumarate hydratase and malate dehydrogenase (PFI1340W and PFF0895w, respectively, in the mitochondrial electron transport chain and TCA cycle), and fructose hexokinase (PFF1155W in glycolysis; Table 5). We also identified 12 reactions that were essential only in artemisinin sensitive parasites (Table 6). Drugs targeting these reactions should not be combined with artemisinin, as they would not kill (and may select for) resistant parasites. Fortunately, no existing drug targets were found in this list of essential genes and reactions (Table 6). Among those identified were sphingomyelin synthase 2 (PFF1215w) and several transport reactions, which furthers our understanding of condition-specific intra-organellar function (Tables 5 & 6). Overall, our systems biology-based approach reveals unique metabolic phenotypes associated with artemisinin sensitivity; these differences were not detected in the original analysis of the expression dataset or by separately analyzing Cambodian or Vietnamese isolates ( and data not shown).
Systems biology approaches enable unbiased analyses of antimalarial resistance phenotypes. Here, we describe a newly curated metabolic network reconstruction of the malaria parasite that can serve as a platform for the analysis of gene expression and other ‘omics data, and as a tool to generate testable hypotheses regarding essential genes and metabolic phenotypes. In particular, we used this network reconstruction to characterize key metabolic dependencies in resistant and sensitive parasites. We revealed emergent patterns in pathway activity, differential utilization of organelles, metabolic flexibility, and targetable weakness of resistant parasites.
Data-driven model curation improves predictive capability
Several P. falciparum reconstructions have been generated since the publication of iTH366, including those highlighting unique developmental stages within the blood-stage asexual cycle by integrating stage-specific expression , de novo reconstructions to implement novel modeling approaches , integrated host and pathogen networks , and those exploring the other life stages of the parasite [73, 74]. iPfal17 represents the most comprehensive and validated metabolic reconstruction of the asexual blood-stage malaria parasite, P. falciparum, to date. With iPfal17, we can simulate growth and predict gene and reaction essentiality and integrate datasets to probe targeted phenotypes, like resistance. It is larger in scope than previous models, includes more gene annotations, and documents literature citations associated with its components (Table 1 & Additional file 3: Table S1, Fig. 2). Moreover, invalid reactions have been removed, improving accuracy (Additional file 3: Table S2). These curation efforts improve the model validity by better recapitulating experimental results, removing functions known to not occur in the asexual blood-stage parasite, and adding functions for which there is experimental evidence. Thus, gene and reaction knockout predictions generated with this model are more accurate. Moreover, iPfal17 has greater interpretability as reaction citations are included and accessible to users.
iPfal17 is similar in functional distribution and scope to other high quality models of apicomplexans, despite its reduced genome size. The P. falciparum genome is 23.3 MB and contains 5423 genes (excluding the antigenic var. genes) [43, 75, 76]; iPfal17 accounts for the function of 987 metabolites, 730 enzymatic reactions, 1195 total reactions, and 488 genes (Table 1). For reference, the network reconstruction for Toxoplasma gondii, with a genome of 80 MB with 8000 genes [77, 78], accounts for 1019 metabolites, 1089 enzymatic reactions, 3387 total reactions, and 527 genes ; with a genome of 32.8 MB and 8272 genes , the Leishmania major network reconstruction accounts for 1101 metabolites, 1047 enzymatic reactions, 1112 total reactions, and 560 genes . These parasites all have notably poor genome annotation (40–60% of the genes are unknown) [43, 75] and, thus, have fewer associated genes than many other reconstructions (e.g. the E. coli and S. cerevisiae reconstructions account for 1366 and 910 genes , respectively [82, 83]).
Intracellular parasites, like Plasmodium, require more exchange and transport reactions as they obtain many nutrients from the host environment [57, 84,85,86,87]. This reliance on the host for metabolic function permits the parasite to increase fitness by reducing its genome and hijacking host function. P. falciparum does just that: the parasite remodels the host erythrocyte, generating a vesicular network for protein translocation and increasing host cell permeability for nutrient acquisition from the host serum [88,89,90,91,92]. Thus, the apicomplexan network reconstructions include more transport reactions, many of which are not genetically mapped. Additionally, we chose to exclude an erythrocytic host compartment from the extracellular environment, despite the parasite’s intra-host growth [57, 66, 68, 93]. Other recent reconstructions [72, 94] have added this compartment, but the erythrocytic compartment is unlikely to improve model function due to the gross disruption of the host membrane as a barrier [57, 66, 68, 93].
We generated gene and reaction essentiality predictions with our curated network model, prior to integration of expression data, and found results largely consistent with previous models  (Table 4). We identified 159 essential reactions and 107 essential metabolic genes (Additional file 3: Table S3 & S4); 24 of these have been empirically tested in cultured P. falciparum parasites (Table 4 , and in P. berghei- Additional file 3: Table S6). iPfal17 better predicts experimentally determined essential reactions than previous models, across a broad set of metabolic pathways (Table 4 and data not shown). iPfal17 predictions fail when essential genes or reactions are involved closely with spontaneous reactions (i.e. lactoylglutathione lyase is downstream of a spontaneous reaction and upstream of nonmetabolic redox products), are in pathways with uncharacterized mechanisms (i.e. plasmepsin II in hemoglobin degradation) or if experimental evidence is contradictory (i.e. heme biosynthesis pathway; Table 4).
Because pharmacological enzyme inhibition can be quite noisy and genetic modification has been challenging in Plasmodium, the development of CRISPR-Cas9 and other technologies will make it possible to integrate new experimental observations into the model with increasing accuracy [95,96,97,98]. Until then, the model can be used to identify enzyme inhibitors with off-target effects. For example, within the heme biosynthesis pathway, pharmacological inhibition of aminolevulinic acid dehydrogenase and protoporphyrinogen oxidase kills blood-stage parasites ; however, disrupting the genes encoding the first (aminolevulinic acid dehydrogenase) and last (ferrochetalase) genes is not lethal in blood-stage parasites [100, 101]. iPfal17 predictions are consistent with the gene knockout experiments in P. falciparum of Ke, et al., suggesting that the enzyme inhibitors used by Ramya, et al. have off target effects (Table 4 and Additional file 3: Table S7). iPfal17 also fails to predict the lethal nature of adenosine deaminase in purine-free conditions . Adenosine deaminase converts adenosine to hypoxanthine; as 38 reactions produce AMP, which then generate hypoxanthine products, we propose adenosine deaminase may be essential for nonmetabolic functions or the inhibitor of adenosine deaminase has off - target effects. Furthermore, these results generate hypotheses about the differential metabolic capabilities of P. falciparum and P. berghei, as experimental results in the rodent parasite conflict with some P. falciparum predictions (Tables 3 & 4 , Additional file 3: Table S6).
Data integration reveals distinct metabolic patterns
The integration of expression data from clinical parasites into our network reconstruction highlights the differential utilization of metabolic genes and reveals metabolic shifts associated with variation in innate artemisinin sensitivity (Additional file 4: Figure S3 & Additional file 5: Figure S4). Enriched metabolic pathways detected in sensitive and resistant models are consistent with previous experimental observations. For example, resistant models are uniquely enriched with genes involved in pyrimidine biosynthesis and mitochondrial redox reactions. This finding is consistent with the importance of mitochondrial function in surviving artemisinin stress [26, 29] and the physical interactions between artemisinin and proteins involved in glycolysis, nucleotide synthesis, and the mitochondria in mammalian cells and P. falciparum [103,104,105]. Additionally, the metabolic disruption of the redox reactions in the electron transport chain upon artemisinin treatment (via decreased production of orotate and fumarate, presumably via dihydroorotate dehydrogenase and succinate dehydrogenase enzymes [22, 28, 106]) suggests that changes in these pathways may be important for survival in the presence of the drug. Thus, this metabolic network analysis approach allows us to filter out noise from diverse clinical isolates to identify alternative utilization of pathways associated with artemisinin resistance. However, due to the nature of this type of analysis, these enrichment results do not implicate specific reactions that are uniquely active in artemisinin sensitive or resistant parasites.
Condition-specific models have unique metabolic requirements
Upon integration of expression data and the identification of differentially utilized pathways above, we next used these models to predict targetable differences in sensitive and resistant parasites by identifying reactions that are essential within the context of the metabolic network (Additional file 3: Table S5, Tables 5 & 6). We identified (1) differences in intra-organellar function, (2) metabolic flexibility of scavenging and biosynthesis pathways, and (3) targetable weakness of resistant parasites. These metabolic shifts primarily reside in mitochondrial metabolism, as well as folate and polyamine metabolism. Together, these results highlight the overall plasticity of P. falciparum metabolism and opportunities for further development of potential drug targets.
Interestingly, several transport reactions are found to be differentially essential in our constrained models (Tables 5 & 6). Many transport reactions (79.5%) have no associated gene due to the incomplete characterization of the P. falciparum genome (Fig. 2). They are included in the model due to biochemical evidence or functional necessity (i.e. a metabolite is produced in one compartment but it is a substrate for an enzyme in another). Transcriptomic data integration does not constrain their behavior explicitly: expression integration reduces the total number of reactions in a model, forcing transport of metabolites among organelles if within-compartment biosynthesis is non-functional. Function within organelles requires transport and loss of function reduces transport needs. Specifically, several mitochondrial and apicoplast transport reactions are uniquely essential in the sensitive and resistant parasite populations (Fig. 4). In resistant models, this includes the mitochondrial transport of metabolites associated with the TCA cycle and electron transport chain (fumarate, oxaloacetate, and NADPH) and those involved in generation of folates (tetrahydrofolate, glycine, CO2, and NH4+) (Fig. 4a). In sensitive models, apicoplast transport of ADP, ATP, and phosphate is essential (Fig. 4a). Overall, these results indicate that sensitive and resistant parasites are differentially utilizing pathways within these organelles, and have unique requirements for transport of essential substrates. This observation is consistent with previous studies and our enrichment results highlighting the influence of mitochondrial metabolism on survival in the presence of artemisinin [26, 29]. Moreover, oxygen transport into the cell and then into the mitochondria is only essential in sensitive parasites, further predicting differential use of the mitochondria in these parasites as oxygen serves as the terminal step in the electron transport chain. Resistant parasites are predicted to generate oxygen within the mitochondria via superoxide dismutase as opposed to transport (Fig. 4a).
We also identify differential utilization of transport pathways from the extracellular environment into the parasite. Plasmodium metabolism contains redundancies; for many essential metabolites, the parasite’s genome encodes one or more biosynthetic pathways, while there is also evidence for a parallel host-scavenging pathway  (e.g. lipid  and amino acid [56, 63] scavenging). Upon model integration, we find that artemisinin resistant and sensitive parasites utilize some of these metabolic pathways in alternative ways (Fig. 4a). Bioinformatic analyses indicate Plasmodium can either scavenge or synthesize putrescine  and adenosylmethionine  (two essential polyamines and precursors to spermidine [50, 107]). Similar redundancy has been identified for the acquisition of p-aminobenzoate, a folate precursor generated by branch of glycolysis necessary for nucleotide synthesis ([43, 108]; Fig. 4a&b, Tables 5 & 6). These metabolites are measurable via blood sample metabolomics [108, 109]; therefore, host scavenging is a viable option for blood-stage parasites. We predict that sensitive parasites rely on the import of putrescine, adenosylmethionine, and p-aminobenzoate. Resistant parasite expression supports either host scavenging or direct biosynthesis due to parasite survival upon reaction knockout in silico. Thus, we expect that resistant parasites are more metabolically flexible for these metabolites; perhaps resistant parasites have failed to appropriately modulate their transition to the nutrient-rich blood-stage environment, and this unexpected flexibility is evolutionarily beneficial once confronted with artemisinin treatment.
Interestingly, recent metabolomics studies demonstrate that intra-parasitic putrescine levels are decreased upon artemisinin treatment . Furthermore, protein interaction studies indicate artemisinin covalently binds with spermidine synthase and adenosylmethionine synthetase . Activity in both the biosynthetic and scavenging pathway of putrescine and adenosylmethionine may allow resistant parasites to compensate for artemisinin’s effect on polyamines. The essential role of polyamines is well established in Plasmodium [111, 112]. In other organisms, these compounds stabilize DNA and RNA  and signal a pause in the cell cycle . In the presence of artemisinin, perhaps polyamines act to stabilize the genome from oxidative stress [24, 30, 33] and trigger dormancy [18, 19]. As resistant parasites are more likely to survive dormancy, flexibility in polyamine metabolism could provide more routes for artemisinin survival [29, 115].
Our systems biology approach also identifies metabolic weaknesses of resistant parasites; these weaknesses can be used to identify drug targets for combination therapies (Fig. 5). For example, we identified the mitochondrial import of fumarate and subsequent conversion to oxaloacetate (via fumarate hydratase, PFI1340W, and malate dehydrogenase, PFF0895W) to be uniquely essential in resistant parasites (Fig. 5a, Table 5). Expression data from sensitive parasites supports mitochondrial import of malate and utilization of malate:quinone oxidoreductase (PFF0815W) to generate oxaloacetate from malate, bypassing the need for fumarate and the associated enzymes, fumarate hydratase and malate dehydrogenase. We predict that inhibitors of fumarate transport or fumarate hydratase and malate dehydrogenase would specifically kill artemisinin resistant parasites, offering an example of enhanced metabolic flexibility of sensitive parasites and a potential artemisinin-combination therapy target. The TCA cycle is essential during the mosquito-stage of parasite development [61, 116], but not the blood-stage [60, 61]; this once again highlights the possibility that resistant parasites exhibit incomplete transition to the metabolic state most appropriate for nutrient-rich blood.
Additionally, we identified serine hydroxymethyltransferase (SHMT) and thiamine diphosphokinase as potential drug targets of resistant parasites (Table 5 , Fig. 5b & c); see below for discussion of SHMT. Both the import of thiamine and thiamine diphosphokinase are essential only in resistant parasites (Fig. 5c), and we predict inhibition of import or enzyme activity would specifically target resistant parasites. These reactions are relatively uncharacterized as the parasite can likely synthesize thiamine diphosphate (vitamin B1) de novo . Thus, this approach can generate novel hypotheses and be utilized for the identification of novel drug targets, and, importantly, targets to help prevent the development of resistance.
Data-driven model implementation highlights knowledge gaps
Although iPfal17 represents our best understanding of intra-erythrocytic P. falciparum biochemistry as the most comprehensive reconstruction to date, predictions occasionally contradict published experimental results. These results illuminate experimental complexities and incompletely characterized pathways. For example, our model predicted that cytosolic SHMT is only essential in resistant parasites (Fig. 5b left). In sensitive parasites, the essential metabolites can be generated by SHMT or the mitochondrial glycine cleavage system, given the reversible nature of these enzymes [118, 119]. Therefore, in our sensitive models, neither SHMT nor the glycine cleavage system is essential when knocked out individually. This observation conflicts with the literature, as SHMT is essential in cultured parasites [118, 120, 121]. Thus, iPfal17 is unable to predict this intricacy of parasite metabolism, revealing interesting regulatory effects, an uncharacterized location dependency for metabolite generation, or in vivo/in vitro differences in enzyme reversibility (Additional file 6).
Similarly, model integration reveals that protein localization influences essentiality predictions. We predicted that the cyclical oxidization and reduction of glutathione, a key regulator of oxidative stress [122,123,124,125], and supporting reactions were essential only in resistant parasites when the glutathione redox system was located within the mitochondria (data not shown). This is consistent with artemisinin’s induction of reactive oxygen species, the parasite’s obvious need to survive this stress [24, 33,34,35], data showing artemisinin sensitivity is correlated with glutathione levels in rodent Plasmodium , and artemisinin’s inhibition of mammalian glutathione s-transferases . However, upon moving these reactions to the cytosolic and apicoplast compartments (as supported by ), these reactions were no longer essential. Thus, model analysis challenges the integration of previously incomparable datasets by demonstrating that this localization and role of glutathione yield different predictions. Future studies will be required to clarify these findings.
Here, we have presented both a novel blood-stage-specific P. falciparum metabolic network reconstruction, iPfal17, and investigation of the metabolic differences between artemisinin sensitive and resistant parasites. Antimalarial resistance is a major public health problem and we demonstrate that constraint-based modeling can be used to reveal metabolic shifts that arise with or in support of the resistant phenotype and discrepancies between otherwise incomparable datasets. We find inherent differences in artemisinin resistant and sensitive parasite metabolism, even before artemisinin treatment. Artemisinin resistant parasites have major metabolic shifts in the mitochondria and in the synthesis of folates and polyamines, indicating incomplete transition to the metabolic state most appropriate for the blood-stage environment. These findings generate areas of future research to elucidate Plasmodium biochemistry, understand the evolution of artemisinin resistant parasites, and tackle antimalarial resistance.
Normalized preprocessed data was obtained from GEO (GSE59097) . Probes on the microarray platform GPL18893 were annotated using NCBI’s stand-alone BLAST correcting the gene labels for 647 probes. Only top hits were used; specifically, hits with greater than 95% identity, no gaps, and a score of over 100 were used (Additional file 3: Table S8). The R package limma was used to compare artemisinin sensitive and resistant samples collected from Cambodia and Vietnam . Samples with predominantly ring-stage parasites with no detectable gametocytes were used. Resistant parasites were defined as both having at least one mutant Kelch13 allele and a parasite clearance half-life of greater than 5 h (Fig. 1) [49, 128]. Sensitive parasites were defined by having at least no mutant Kelch13 alleles and a parasite clearance half-life of less than 5 h. Random Forest classifiers were built using the R package randomForest, using all ring-stage samples . The metadata classifier used the variables listed in Additional file 2: Figure S2, as outlined in the original study . Cambodian and Vietnamese ring-stage transcriptomes were compared separately to ensure patterns associated with resistance status were reproducible across phylogenies. These countries were chosen for large number of isolates and prevalence of resistance. Microarray probes were screened to remove non-metabolic genes and to keep only one probe per gene (consistent with standard practice). Multiple testing correction was conducted using a false discovery rate [130, 131].
Gene expression data with calculations of fold changes and associated adjusted p-value were incorporated into our curated model using the Metabolic Adjustment for Differential Expression (MADE) algorithm. MADE utilizes statistical significance of gene expression changes along with network context to assign binary gene states (‘on’/‘off’) to each metabolic gene. This constrains the network by limiting flux through reactions mapped to ‘off’ genes while maintaining growth, or a similar objective. An 80% growth threshold was used given that there is no reported evidence that resistant and sensitive parasites produce variable biomass as measured by the size of ring-stage parasites; while varying this threshold affects sensitive parasite biomass yield, it does not affect essentiality predictions (data not shown). Essential genes were predicted for the resultant condition-specific models (Fig. 3) by conducting single gene and reaction deletions with established algorithms . Consensus lethal gene and reaction deletions from the Cambodian and Vietnamese parasite models were used.
Flux analysis and metabolic tasks
Flux balance analysis (FBA) is an approach to explore metabolic phenotypes in silico . FBA simulates steady-state flux values for each of the network’s reactions that maximize subsequent flux through an objective function given a set of constraints. We chose biomass production as the objective reaction, consistent with previous studies interrogating gene essentiality [50, 53, 79, 134], and permitted flux through all transport reactions. Constraints on the system include conservation of mass, reversibility of reactions, and reaction localization. Flux variability analysis (FVA) uses a related approach to find the range of fluxes permissible given system constraints .
We simulated in vitro experiments and in vivo data to evaluate the model; these are our metabolic ‘tasks’ that the reconstruction should pass. We simulate in vitro growth requirements by modifying media components or access to particular metabolites. Metabolite import or production was eliminated from the reconstruction, and subsequent biomass production was observed. Effects of enzyme inhibition, gene knockouts, and metabolite production were also used to evaluate the model. Lethal modifications were defined as changes that resulted in no production of biomass; growth-reducing modifications were defined as producing less than 90% of unconstrained flux value [81, 134].
The COBRA Toolbox 2015, Tiger Toolbox (version 1.3.1), and MATLAB R2013b were used for model generation and flux simulations.
Manual curation of an existing P. falciparum metabolic network reconstruction  was conducted by a literature review and reference to generic and Plasmodium-specific databases (KEGG, Expasy, and PlasmoDB, MPMP) [43, 136,137,138]. Data obtained from these sources were used to evaluate the inclusion of reactions as well as their stoichiometry, reversibility, localization, and gene annotations. Genetically and biochemically supported reactions were kept and new reactions were added. Reactions were removed if (1) explicitly determined to be false or (2) were nonfunctional and not supported biochemically or genetically. Spontaneous reactions (reactions that occur without enzymes) are noted to differentiate from orphan reactions (reactions with unknown enzyme catalysts).
In order to assess gene essentiality, we used a biomass reaction as the modeling objective function. Thus, flux through this reaction, simulating cellular growth, was maximized for all in silico experimental procedures. We used the biomass reaction from a previous study  with modifications. Curation of the biomass reaction was informed by metabolites detected in metabolomics studies [28, 56,57,58]; if possible, metabolite ratios were predicted from metabolomics data. We curated the biomass reaction with consideration of published essentiality data; metabolites detected in metabolomics experiments with no known catabolism or import pathways were excluded from the biomass reaction.
We predicted essentiality by performing single deletion studies with both genes and reactions and double gene deletion studies in our curated model and each expression-constrained sensitive and resistant models. All simulations were performed in an in silico red blood cell environment (Additional file 3: Table S9). Gene deletions were simulated by removing the gene of interest from the model. This change results in the inhibition of flux through all reactions that require that gene to function. If the model could not produce biomass with these constraints, the gene was deemed essential. Growth reducing phenotypes were also observed and noted. For reaction deletion studies, we removed reactions sequentially. Subsequent growth effects were used to determine reaction essentially. Consensus results for resistant or sensitive models are discussed.
Flux balance analysis
Flux variability analysis
Metabolic Adjustment for Differential Expression algorithm
tricarboxylic acid cycle
Schwartz, E. and T. Lachish, Artemisinin-based combination therapy (ACT) versus atovaquone-proguanil: do not choose between but, rather, combine them. Evidence Based Medicine, 2016. 21(2): p. 64–64.
Olliaro PL, Taylor WR. Developing artemisinin based drug combinations for the treatment of drug resistant falciparum malaria: a review. J Postgrad Med. 2004;50(1):40.
Eastman RT, Fidock DA. Artemisinin-based combination therapies: a vital tool in efforts to eliminate malaria. Nat Rev Microbiol. 2009;7(12):864–74.
Chakraborty A. Emerging drug resistance in plasmodium falciparum: a review of well-characterized drug targets for novel antimalarial chemotherapy. Asian Pac J of Tro Dis. 2016;6(7):581–8.
Cowman AF, et al. Malaria: biology and disease. Cell. 2016;167(3):610–24.
Plowe CV, et al. World antimalarial resistance network (WARN) III: molecular markers for drug resistant malaria. Malar J. 2007;6(1):1–10.
Sidhu AB, Verdier-Pinard D, Fidock DA. Chloroquine resistance in plasmodium falciparum malaria parasites conferred by pfcrt mutations. Science. 2002;298
Guler JL, et al. Atovaquone tolerance in plasmodium falciparum parasites selected for high-level resistance to a dihydroorotate dehydrogenase inhibitor. Antimicrob Agents Chemother. 2015;59(1):686–9.
Herman JD, et al. A genomic and evolutionary approach reveals non-genetic drug resistance in malaria. Genome Biol. 2014;15(11):511.
Gabryszewski, S.J., et al., Evolution of Fitness Cost-Neutral Mutant PfCRT Conferring P. falciparum 4-Aminoquinoline Drug Resistance Is Accompanied by Altered Parasite Metabolism and Digestive Vacuole Physiology. PLOS Pathog, 2016. 12(11): p. e1005976.
Meylan S, et al. Carbon sources tune antibiotic susceptibility in Pseudomonas Aeruginosa via tricarboxylic acid cycle control. Cell Chem Biol. 2017; 24:195-206.
El-Halfawy OM, Valvano MA. Non-genetic mechanisms communicating antibiotic resistance: rethinking strategies for antimicrobial drug design. Expert Opin Drug Discov. 2012;7(10):923–33.
Ashley EA, et al. Spread of artemisinin resistance in plasmodium falciparum malaria. N Engl J Med. 2014;371(5):411–23.
Miotto O, et al. Genetic architecture of artemisinin-resistant plasmodium falciparum. Nat Genet. 2015;47(3):226–34.
Straimer J, et al. K13-propeller mutations confer artemisinin resistance in Plasmodium falciparum clinical isolates. Science. 2015;347(6220):428–31.
Ariey F, et al. A molecular marker of artemisinin-resistant plasmodium falciparum malaria. Nature. 2014;505(7481):50–5.
Brown TS, et al. Plasmodium falciparum field isolates from areas of repeated emergence of drug resistant malaria show no evidence of hypermutator phenotype. Infection, Genetics and Evolution: J Mol Epidemiol Evol Genet Infect Dis. 2015;30:318–22.
Cheng Q, Kyle DE, Gatton ML. Artemisinin resistance in Plasmodium falciparum: A process linked to dormancy? International Journal for parasitology. Drugs and Drug Resist. 2012;2:249–55.
Codd A, et al. Artemisinin-induced parasite dormancy: a plausible mechanism for treatment failure. Malar J. 2011;10
Straimer J, et al. Drug resistance. K13-propeller mutations confer artemisinin resistance in Plasmodium falciparum clinical isolates. Sci. 2015;347(6220):428–31.
Mbengue A, et al. A molecular mechanism of artemisinin resistance in plasmodium falciparum malaria. Nat. 2015;520(7549):683–7.
Ying-Zi Y, Little B, Meshnick SR. Alkylation of proteins by artemisinin: effects of heme, pH, and drug structure. Biochem Pharmacol. 1994;48(3):569–73.
Dalal S, Klemba M. Amino acid efflux by asexual blood-stage plasmodium falciparum and its utility in interrogating the kinetics of hemoglobin endocytosis and catabolism in vivo. Mol Biochem Parasitol. 2015;201(2):116–22.
Klonis N, et al. Artemisinin activity against plasmodium falciparum requires hemoglobin uptake and digestion. Proc Natl Acad Sci U S A. 2011;108(28):11405–10.
Wang J, et al. Artemisinin Directly Targets Malarial Mitochondria through Its Specific Mitochondrial Activation. PLoS One. 2010;5(3):e9582.
Chen N, et al. Fatty acid synthesis and pyruvate metabolism pathways remain active in Dihydroartemisinin-induced dormant ring stages of plasmodium falciparum. Antimicrob Agents Chemother. 2014;58(8):4773–81.
Vega-Rodríguez J, et al. Implications of glutathione levels in the plasmodium berghei response to chloroquine and artemisinin. PLoS One. 2015;10(5):e0128212.
Cobbold SA, et al. Metabolic dysregulation induced in plasmodium falciparum by Dihydroartemisinin and other front-line antimalarial drugs. J Infect Dis. 2016;213(2):276–86.
Peatey CL, et al. Mitochondrial membrane potential in a small subset of artemisinin-induced dormant plasmodium falciparum parasites in vitro. J Infect Dis. 2015;212(3):426–34.
Meshnick SR. Artemisinin: mechanisms of action, resistance and toxicity. Int J Parasitol. 2002;32(13):1655–60.
Eckstein-Ludwig U, et al. Artemisinins target the SERCA of plasmodium falciparum. Nature. 2003;424(6951):957–61.
Golenser J, et al. Current perspectives on the mechanism of action of artemisinins. Int J Parasitol. 2006;36(14):1427–41.
Efferth T, Oesch F. Oxidative stress response of tumor cells: microarray-based comparison between artemisinins and anthracyclines. Biochem Pharmacol. 2004;68(1):3–10.
Antoine T, et al. Rapid kill of malaria parasites by artemisinin and semi-synthetic endoperoxides involves ROS-dependent depolarization of the membrane potential. J Antimicrob Chemother. 2014;69(4):1005–16.
Sun C, et al. Two distinct and competitive pathways confer the cellcidal actions of artemisinins. Microbial Cell. 2015;2(1):14–25.
Li, W., et al., Yeast model uncovers dual roles of mitochondria in the action of artemisinin. PLoS Genet, 2005. 1(3): p. e36.
Kamau E, et al. K13-propeller polymorphisms in plasmodium falciparum parasites from sub-Saharan Africa. J Infect Dis. 2015;211(8):1352–5.
Isozumi R, et al. Novel mutations in K13 propeller Gene of artemisinin-resistant <i>plasmodium falciparum</i>. Emerg Infect Dis. 2015;21(3):490–2.
Klonis N, et al. Altered temporal response of malaria parasites determines differential sensitivity to artemisinin. Proc Natl Acad Sci U S A. 2013;110(13):5157–62.
Mok S, et al. Artemisinin resistance in Plasmodium falciparum is associated with an altered temporal pattern of transcription. BMC Genomics. 2011:12(1).
Teuscher F, et al. Artemisinin-induced dormancy in plasmodium falciparum: duration, recovery rates, and implications in treatment failure. J Infect Dis. 2010;202(9):1362–8.
Witkowski B, et al. Increased tolerance to artemisinin in plasmodium falciparum is mediated by a quiescence mechanism. Antimicrob Agents Chemother. 2010;54(5):1872–7.
Collaborativea TPGD, PlasmoDB: An integrative database of the Plasmodium falciparum genome. Tools for accessing and analyzing finished and unfinished sequence data. Nucleic Acids Res. 2001;29(1):66–9.
Fidock DA, et al. Mutations in the P. falciparum digestive vacuole transmembrane protein PfCRT and evidence for their role in chloroquine resistance. Mol Cell. 2000:6.
Peterson DS, Walliker D, Wellems TE. Evidence that a point mutation in dihydrofolate reductase-thymidylate synthase confers resistance to pyrimethamine in falciparum malaria. Proc Natl Acad Sci U S A. 1988;85(23):9114-9118.
Siregar JE, et al. Direct evidence for the atovaquone action on the plasmodium cytochrome bc1 complex. Parasitol Int. 2015;64(3):295–300.
Phillips MA, Rathod PK. Plasmodium dihydroorotate dehydrogenase: a promising target for novel anti-malarial chemotherapy. Infect Disord Drug Targets. 2010;10(3):226–39.
Fuhrer T, et al. Genomewide landscape of gene–metabolome associations in Escherichia coli. Mol Syst Biol. 2017;13(1):907.
Mok S, et al. Population transcriptomics of human malaria parasites reveals the mechanism of artemisinin resistance. Science. 2015;347(6220):431–5.
Plata G, et al. Reconstruction and flux-balance analysis of the plasmodium falciparum metabolic network. Mol Syst Biol. 2010;6
Painter HJ, et al. Specific role of mitochondrial electron transport in blood-stage plasmodium falciparum. Nature. 2007;446(7131):88–91.
Langer RC, Vinetz JM. Plasmodium ookinete-secreted chitinase and parasite penetration of the mosquito peritrophic matrix. Trends Parasitol. 2001;17(6):269–72.
Thiele I, Palsson BO. A protocol for generating a high-quality genome-scale metabolic reconstruction. Nat Protocols. 2010;5(1):93–121.
Sturm A, et al. Mitochondrial ATP synthase is dispensable in blood-stage plasmodium berghei rodent malaria but essential in the mosquito phase. Proc Natl Acad Sci. 2015;112(33):10216–23.
Ginsburg H. Abundant proton pumping in plasmodium, but why? Trends Parasitol. 2002;18(11):483–6.
Gulati S, et al. Profiling the essential nature of lipid metabolism in asexual blood and gametocyte stages of plasmodium falciparum. Cell Host Microbe. 2015;18(3):371–81.
Olszewski KL, et al. Host-parasite interactions revealed by plasmodium falciparum metabolomics. Cell Host Microbe. 2009;5(2):191–9.
Sana TR, et al. Global Mass Spectrometry Based Metabolomics Profiling of Erythrocytes Infected with Plasmodium falciparum. PLoS One. 2013;8(4):e60840.
Biddau M, Müller S. Carbon Metabolism of Plasmodium falciparum. Compr Anal Parasite Biol: Metab to Drug Discov. 2016:371.
Ke H, et al. Genetic investigation of tricarboxylic acid metabolism during the plasmodium falciparum life cycle. Cell Rep. 2015;11(1):164–74.
MacRae JI, et al. Mitochondrial metabolism of sexual and asexual blood stages of the malaria parasite plasmodium falciparum. BMC Biol. 2013;11(1):67.
Yeh I. Computational analysis of plasmodium falciparum metabolism: organizing genomic information to facilitate drug discovery. Genome Res. 2004;14(5):917–24.
Liu J, et al. Plasmodium falciparum ensures its amino acid supply with multiple acquisition pathways and redundant proteolytic enzyme systems. Proc Natl Acad Sci U S A. 2006;103(23):8840–5.
Krugliak M, Zhang J, Ginsburg H. Intraerythrocytic plasmodium falciparum utilizes only a fraction of the amino acids derived from the digestion of host cell cytosol for the biosynthesis of its proteins. Mol Biochem Parasitol. 2002;119(2):249–56.
Asahi H, et al. Hypoxanthine: a low molecular weight factor essential for growth of erythrocytic plasmodium falciparum in a serum-free medium. Parasitol. 1996;113(01):19–23.
Geary TG, et al. Nutritional requirements of plasmodium falciparum in culture. III. Further observations on essential nutrients and antimetabolites. The J Protozool. 1985;32(4):608–13.
Miller LH, et al. The pathogenic basis of malaria. Nature. 2002;415(6872):673–9.
Geary TG, Divo AA, Jensen JB. Nutritional requirements of plasmodium falciparum in culture. II. Effects of antimetabolites in a semi defined medium. The J Protozool. 1985;32(1):65–9.
Jensen PA, Papin JA. Functional integration of a metabolic network model and expression data without arbitrary thresholding. Bioinform. 2011;27(4):541–7.
Fang X, Reifman J, Wallqvist A. Modeling metabolism and stage-specific growth of plasmodium falciparum HB3 during the intraerythrocytic developmental cycle. Mol BioSyst. 2014;10(10):2526–37.
Chiappino-Pepe A, et al. Bioenergetics-based modeling of plasmodium falciparum metabolism reveals its essential genes, nutritional requirements, and thermodynamic bottlenecks. PLoS Comput Biol. 2017;13(3):e1005397.
Wallqvist A, et al. Metabolic host responses to malarial infection during the intraerythrocytic developmental cycle. BMC Syst Biol. 2016;10(1):58.
Tymoshenko S, et al. Functional genomics of plasmodium falciparum using metabolic modelling and analysis. Briefings in Functional Genomics. 2013;12(4):316–27.
Bazzani S, Hoppe A, Holzhütter H-G. Network-based assessment of the selectivity of metabolic drug targets in Plasmodium falciparum with respect to human liver metabolism. BMC Syst Biol. 2012:6(1).
Aurrecoechea, C., et al., EuPathDB: the eukaryotic pathogen genomics database resource. Nucleic Acids Res, 2017: 45:D581-D591.
Gardner MJ, et al. Genome sequence of the human malaria parasite plasmodium falciparum. Nat. 2002;419(6906):498–511.
Gajria B, et al. ToxoDB: an integrated Toxoplasma gondii database resource. Nucleic Acids Res. 2008;36(Database issue):D553–6.
Xia D, et al. The proteome of Toxoplasma gondii: integration with the genome provides novel insights into gene expression and annotation. Genome Biol. 2008;9(7):R116.
Tymoshenko S, et al. Metabolic needs and capabilities of <italic>Toxoplasma gondii </italic> through combined Computational and experimental analysis. PLoS Comput Biol. 2015;11(5):e1004261.
Ivens AC, et al. The genome of the Kinetoplastid parasite, Leishmania major. Sci. 2005;309(5733):436.
Chavali AK, et al. Systems analysis of metabolism in the pathogenic trypanosomatid Leishmania major. Mol Syst Biol. 2008;4(1):177.
Orth, J.D., et al., A comprehensive genome-scale reconstruction of Escherichia coli metabolism—2011. Mol Syst Biol, 2011. 7: p. 535–535.
Heavner BD, Price ND. Comparative analysis of yeast metabolic network models highlights progress, opportunities for metabolic reconstruction. PLoS Comput Biol. 2015;11(11):e1004530.
McConville M. Open questions: microbes, metabolism and host-pathogen interactions. BMC Biol. 2014;12(1):18.
Wellems TE, Fairhurst RM. An evolving picture of the interactions between malaria parasites and their host erythrocytes. Cell Res. 2012;22(3):453–6.
Imlay L, Odom AR. Isoprenoid metabolism in apicomplexan parasites. Current Clinical Microbiology Reports. 2014;1(3–4):37–50.
Mazumdar J, Striepen B. Make it or take it: fatty acid metabolism of apicomplexan parasites. Eukaryot Cell. 2007;6(10):1727–35.
Lanzer M, et al. Maurer's clefts: a novel multi-functional organelle in the cytoplasm of plasmodium falciparum-infected erythrocytes. Int J Parasitol. 2006;36(1):23–36.
Bannister LH, et al. A brief illustrated guide to the ultrastructure of plasmodium falciparum asexual blood stages. Parasitol Today. 2000;16(10):427–33.
Baumeister S, et al. Evidence for the involvement of plasmodium falciparum proteins in the formation of new permeability pathways in the erythrocyte membrane. Mol Microbiol. 2006;60(2):493–504.
Ginsburg H, et al. New permeability pathways induced in membranes of plasmodium falciparum infected erythrocytes. Mol Biochem Parasitol. 1983;8(2):177–90.
Staines HM, et al. Solute transport via the new permeability pathways in plasmodium falciparum–infected human red blood cells is not consistent with a simple single-channel model. Blood. 2006;108(9):3187–94.
Divo AA, et al. Nutritional requirements of plasmodium falciparum in culture. I. Exogenously supplied dialyzable components necessary for continuous growth. J Protozool. 1985;32(1):59–64.
Phaiphinit S, et al. In silico multiple-targets identification for heme detoxification in the human malaria parasite plasmodium falciparum. Infect Genet Evol. 2016;37:237–44.
Ghorbal M, et al. Genome editing in the human malaria parasite plasmodium falciparum using the CRISPR-Cas9 system. Nat Biotech. 2014;32(8):819–21.
Lee MCS, Fidock DA. CRISPR-mediated genome editing of plasmodium falciparum malaria parasites. Genome Med. 2014;6(8):63.
Wagner JC, et al. Efficient CRISPR-Cas9-mediated genome editing in plasmodium falciparum. Nat Meth. 2014;11(9):915–8.
Lu J, et al. A redesigned CRISPR/Cas9 system for marker-free genome editing in plasmodium falciparum. Parasit Vectors. 2016;9(1):1.
Ramya TNC, et al. Inhibitors of nonhousekeeping functions of the apicoplast defy delayed death in plasmodium falciparum. Antimicrob Agents Chemother. 2007;51(1):307–16.
Ke H, et al. The heme biosynthesis pathway is essential for plasmodium falciparum development in mosquito stage but not in blood stages. J Biol Chem. 2014;289(50):34827–37.
Nagaraj VA, et al. Malaria parasite-synthesized Heme is essential in the mosquito and liver stages and complements host Heme in the blood stages of infection. PLoS Pathog. 2013;9(8):e1003522.
Ho M-C, et al. Structural and metabolic specificity of methylthiocoformycin for malarial adenosine deaminases. Biochemistry. 2009;48(40):9618–26.
Zhou Y, Li W, Xiao Y. Profiling of multiple targets of artemisinin activated by hemin in cancer cell proteome. ACS Chem Biol. 2016;11(4):882–8.
Ismail HM, et al. Artemisinin activity-based probes identify multiple molecular targets within the asexual stage of the malaria parasites plasmodium falciparum 3D7. Proc Natl Acad Sci. 2016;113(8):2080–5.
Prieto JH, et al. Large-Scale Differential Proteome Analysis in Plasmodium falciparum under Drug Treatment. PLoS One. 2008;3(12):e4098.
Creek DJ, et al. Metabolomics-based screening of the malaria box reveals both novel and established mechanisms of action. Antimicrob Agents Chemother. 2016;60(11):6650–63.
Pretzel J, et al. Characterization and redox regulation of plasmodium falciparum methionine adenosyltransferase. J Biochem. 2016;160(6):355–67.
Salcedo-Sora JE, et al. The molecular basis of folate salvage in plasmodium falciparum characterization of two folate transporters. J Biol Chem. 2011;286(52):44659–68.
Wishart DS, et al. HMDB: the human metabolome database. Nucleic Acids Res. 2007;35(suppl 1):D521–6.
Wang J, et al. Haem-activated promiscuous targeting of artemisinin in plasmodium falciparum. Nat Commun. 2015;6:10111.
Singh S, et al. Characterization of Simian malarial parasite (plasmodium knowlesi)-induced putrescine transport in rhesus monkey erythrocytes a NOVEL PUTRESCINE CONJUGATE ARRESTS IN VITRO GROWTH OF SIMIAN MALARIAL PARASITE (PLASMODIUM KNOWLESI) AND CURES MULTIDRUG RESISTANT MURINE MALARIA (Plasmodium YOELII) INFECTION IN VIVO. J Biol Chem. 1997;272(21):13506–11.
Assaraf YG, et al. Effect of polyamine depletion on macromolecular synthesis of the malarial parasite, plasmodium falciparum, cultured in human erythrocytes. Biochem J. 1987;242(1):221–6.
Stevens L. The biochemical role of naturally occurring polyamines in nucleic acid synthesis. Biol Rev. 1970;45(1):1–25.
Mandal S, et al. Depletion of cellular polyamines, spermidine and spermine, causes a total arrest in translation and growth in mammalian cells. Proc Natl Acad Sci. 2013;110(6):2169–74.
Teuscher F, et al. Phenotypic changes in artemisinin-resistant plasmodium falciparum lines in vitro: evidence for decreased sensitivity to dormancy and growth inhibition. Antimicrob Agents Chemother. 2012;56(1):428–31.
Srivastava A, et al. Stage-specific changes in plasmodium metabolism required for differentiation and adaptation to different host and vector environments. PLoS Pathog. 2016;12(12):e1006094.
Müller S, Kappes B. Vitamin and cofactor biosynthesis pathways in plasmodium and other apicomplexan parasites. Trends Parasitol. 2007;23(3):112–21.
Maenpuen S, et al. Characterization of plasmodium falciparum serine hydroxymethyltransferase—a potential antimalarial target. Mol Biochem Parasitol. 2009;168(1):63–73.
Salcedo E, Sims PFG, Hyde JE. A glycine-cleavage complex as part of the folate one-carbon metabolism of plasmodium falciparum. Trends Parasitol. 2005;21(9):406–11.
Pornthanakasem W, et al. Plasmodium serine hydroxymethyltransferase: indispensability and display of distinct localization. Malar J. 2012;11(1):387.
França TCC, Pascutti PG, Ramalho TC. A three-dimensional structure of plasmodium falciparum serine hydroxymethyltransferase in complex with glycine and 5-formyl-tetrahydrofolate. Homology modeling and molecular dynamics. Biophys Chem. 2005;115(1):1–10.
Schulz JB, et al. Glutathione, oxidative stress and neurodegeneration. Eur J Biochem. 2000;267(16):4904–11.
Mittler R. Oxidative stress, antioxidants and stress tolerance. Trends Plant Sci. 2002;7(9):405–10.
Becker K, et al. Oxidative stress in malaria parasite-infected erythrocytes: host–parasite interactions. Int J Parasitol. 2004;34(2):163–89.
Färber PM, et al. Molecular cloning and characterization of a putative glutathione reductase gene, the PfGR2 gene, from plasmodium falciparum. Eur J Biochem. 1996;239(3):655–61.
Kehr S, et al. Compartmentation of redox metabolism in malaria parasites. PLoS Pathog. 2010;6(12):e1001242.
Ritchie, M.E., et al., Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res, 2015. 43(7): p. e47-e47.
White LJ, et al. Defining the In Vivo Phenotype of Artemisinin-Resistant Falciparum Malaria: A Modelling Approach. PLoS Med. 2015;12(4):e1001823.
Liaw A, Wiener M. Classification and regression by randomForest. R News. 2002;2(3):18–22.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc. Series B (Methodological). 1995:289–300.
Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. Ann Stat. 2001:1165–88.
Becker SA, et al. Quantitative prediction of cellular metabolism with constraint-based models: the COBRA toolbox. Nat Protoc. 2007;2
Orth JD, Thiele I, Palsson BO. What is flux balance analysis? Nat Biotech. 2010;28(3):245–8.
Oberhardt MA, et al. Metabolic network analysis of Pseudomonas Aeruginosa during chronic cystic fibrosis lung infection. J Bacteriol. 2010;192
Gudmundsson S, Thiele I. Computationally efficient flux variability analysis. BMC Bioinform. 2010;11(1):489.
Kanehisa M, et al. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 2016;44(D1):D457–62.
Gasteiger E, et al. ExPASy: the proteomics server for in-depth protein knowledge and analysis. Nucleic Acids Res. 2003;31(13):3784–8.
Ginsburg H. Progress in in silico functional genomics: the malaria metabolic pathways database. Trends Parasitol. 2006;22(6):238–40.
Yeh E, DeRisi JL. Chemical Rescue of Malaria Parasites Lacking an Apicoplast Defines Organelle Function in Blood-Stage Plasmodium falciparum. PLoS Biol. 2011;9(8):e1001138.
Woodrow CJ, Burchmore RJ, Krishna S. Hexose permeation pathways in plasmodium falciparum-infected erythrocytes. Proc Natl Acad Sci U S A. 2000;97(18):9931–6.
Anfinsen CB, et al. STUDIES ON MALARIAL PARASITES : VIII. FACTORS AFFECTING THE GROWTH OF PLASMODIUM KNOWLESI IN VITRO. J Exp Med. 1946;84(6):607–21.
Jiang L, et al. Potent and selective activity of a combination of thymidine and 1843U89, a folate-based thymidylate synthase inhibitor, against plasmodium falciparum. Antimicrob Agents Chemother. 2000;44(4):1047–50.
Yu M, et al. The fatty acid biosynthesis enzyme FabI plays a key role in the development of liver-stage malarial parasites. Cell Host Microbe. 2008;4(6):567–78.
Vaughan AM, et al. Type II fatty acid synthesis is essential only for malaria parasite late liver stage development. Cell Microbiol. 2009;11(3):506–20.
Deng X, et al. Structural plasticity of malaria dihydroorotate dehydrogenase allows selective binding of diverse chemical scaffolds. J Biol Chem. 2009;284(39):26999–7009.
McRobert L, McConkey GA. RNA interference (RNAi) inhibits growth of plasmodium falciparum. Mol Biochem Parasitol. 2002;119(2):273–8.
Nguyen, C., et al., Deoxyuridine triphosphate nucleotidohydrolase as a potential antiparasitic drug target. (0022–2623 (Print)).
Thornalley, P.J., R.J. Strath M Fau - Wilson, and R.J. Wilson, Antimalarial activity in vitro of the glyoxalase I inhibitor diester, S-p-bromobenzylglutathione diethyl ester. (0006–2952 (Print)).
Hanada K, et al. Plasmodium falciparum phospholipase C hydrolyzing sphingomyelin and lysocholinephospholipids is a possible target for malaria chemotherapy. J Exp Med. 2002;195(1):23–34.
Silva, A.M., et al., Structure and inhibition of plasmepsin II, a hemoglobin-degrading enzyme from Plasmodium falciparum. (0027–8424 (Print)).
Hoepfner D, et al. Selective and specific inhibition of the plasmodium falciparum lysyl-tRNA synthetase by the fungal secondary metabolite cladosporin. Cell Host Microbe. 2012;11(6):654–63.
Patzewitz EM, Wong EH, Muller S. Dissecting the role of glutathione biosynthesis in plasmodium falciparum. Mol Microbiol. 2012;83(2):304–18.
Sturm N, et al. Identification of proteins targeted by the thioredoxin superfamily in plasmodium falciparum. PLoS Pathog. 2009;5(4):e1000383.
Belorgey D, Lanfranchi DA, Davioud-Charvet E. 1,4-naphthoquinones and other NADPH-dependent glutathione reductase-catalyzed redox cyclers as antimalarial agents. Curr Pharm Des. 2013;19(14):2512–28.
Muller S. Role and regulation of glutathione metabolism in plasmodium falciparum. Mol. 2015;20(6):10511–34.
Pastrana-Mena R, et al. Glutathione reductase-null malaria parasites have normal blood stage growth but arrest during development in the mosquito. J Biol Chem. 2010;285(35):27045–56.
We acknowledge the members of the Guler, Papin, and Petri labs for their thoughtful conversations and insight. We would also like to thank Dr. Paul Jensen, Dr. Edik Blais, and Gregory Medlock for project feedback.
The study was financed by institutional funding from the University of Virginia (JG) and by the National Institute of Allergy and Infectious Disease R21AI119881 (JG and JP). MC is supported by an institutional training grant (T32GM008136).
Availability of data and materials
The supporting data and materials of this article are included within the article and additional files; additionally, our model and code is available on https://github.com/gulermalaria/iPfal17.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Distribution of genome-wide expression data demonstrates moderate differential expression between sensitive and resistant parasites. Fold change values from differential expression between sensitive and resistant parasites from Cambodia (A) and Vietnam (B) with significantly differentially expressed genes in red. Fold change is the ratio of mean expression in resistant parasites to mean expression in sensitive parasites, for each respective country. (PDF 820 kb)
Artemisinin resistance is better predicted by metadata classifier than expression classifier. Using full expression profile (A) or metadata (B), including patient and parasite features (see Methods), we can classify samples as artemisinin sensitive or resistant by Random forest analysis. Of the top 25 most important variables (gene probes) in the expression classifier, 12 encoded exported proteins, four genes of complete unknown function, three encoded a putative kinase and putative phosphatases, one encoded a component of dynein, four were uncharacterized genes though to be involved in protein folding or trafficking, and one encoded a transcription factor. Abbreviation key: (all from  unless noted otherwise). aprs_mutation = apicoplast ribosomal protein S10 (PF3D7_1460900.1) mutation ; fd_mutation = ferredoxin (PF3D7_1318100) mutation ; Field_site = location at which blood sample was collected; mdr_mutation = multidrug resistance protein 2 (PF3D7_1447900) mutation; partner_drug = Partner drug (Artemisinin based combination treatment) administered from day 3 onwards; crt_mutation2 = second CRT (PF3D7_0709000)  mutation measured; crt_mutation1 = first CRT (PF3D7_0709000)  mutation measured; Patient_age_yr = patient age in years; pRBC_sampling_vol_uL = Volume of packed RBC collected (uL); RNA_yield_ug = Amount of Total RNA isolated for each sample (μg); Patient_temp_c = patient temperature at time of admission in Celsius; ART_drug = Type and dosage of artemisinin drug given once a day on days 0, 1 and 2; asexual_parasite_count = Total asexual parasite densities per uL on admission; total_parasite_1000 = total number of parasites in whole sample of infected RBC collected (pRBC collection vol. * total parasite count per uL) divided by 1000; SRCC_asexual_stage = Spearman rank correlation coefficient of the gene expression for the isolate sample to the projected hpi; Kmeans_Grp = expression group (see ); Asexual_stage_hpi = Projected hours post invasion (hpi) of the parasite asexual stage; Gender = Patient gender; gam_count = Total gametocyte parasite densities per uL on admission; Hct_percent = patient hematocrit (%) on admission; Sampling_Time_24_hr = Time of sample collection in 24 h format. (PDF 38 kb)
Modifications and reaction additions in iPfal17 curation. Table S2. Reactions deleted from Plata et al. model (iTH366) in generating iPfal17. Table S3. Predicted lethal reactions in wild-type blood-stage Plasmodium falciparum. Table S4: Predicted lethal genes in wild-type blood-stage Plasmodium falciparum. Table S5. Consensus predicted lethal reactions across 4 expression-constrained models. Table S6. PlasmoGem orthologous results. Table S7: Extended Table 4. Table S8. Microarray platform blast results. Table S9: Metabolites in the in silico extracellular environment. (XLSX 1104 kb)
Functional differences in data-driven sensitive and resistant models. Gene states from four condition-specific models, the results of MADE integration, cluster by sensitivity not by location. Active genes in red/blue, with genes removed from expression-constrained models in white. (TIFF 1795 kb)
Artemisinin sensitive and resistant parasites utilize different metabolic genes and pathways. Enrichment analysis of gene utilization in sensitive and resistant parasite models demonstrates functional differences in expression data integration. Consensus gene utilization from resistant and sensitive models (both Cambodian and Vietnamese datasets) were used and compared to unconstrained model. Black = p < 0.001, grey = p < 0.01, light grey = p < 0.05, white = non significant. Abbreviation key: Aminosugars = amino sugar metabolism; ArgPro = arginine and proline metabolism; AsnAsp = asparagine and aspartate metabolism; BiosynthesisCytochrome = biosynthesis of cytochromes; CoABiosynthesis = coenzyme-A biosynthesis; Dolichol = dolichol metabolism; Exchange = exchange reactions; FattyAcidSynthesis = fatty acid synthesis; FolateBiosynthesis = folate biosynthesis; Glu = glutamate metabolism; Glycolysis = glycolysis; GlySer = glycine and serine metabolism; GPIAnchorBiosynthesis = GPI anchor biosynthesis; Hemoglobin = hemoglobin degradation (including hemozoin formation); InositolPhosphate = inositol phosphate metabolism; Isoprenoids = isoprenoid metabolism; LeuIleVal = leucine, isoleucine, and valine metabolism; Lys = lysine metabolism; MannoseFructose = mannose and fructose metabolism; MetPolyamine = methionine and polyamine metabolism; MitochondrialElectronFlow = mitochondrial electron transport chain; MitochondrialTCACycle = mitochondrial tricarboxylic acid cycle; NicotinateNicotinamide = nicotinate and nicotinamide metabolism; Nitrogen = nitrogen metabolism; PentosePhosphateCycle = pentose phosphate cycle; PheTyr = phenylalanine and tyrosine metabolism; Phosphatidylcholine = phosphatidylcholine metabolism; PhosphatidyletanolaminePhosphatidylserine = phosphatidyletanolamine and phosphatidylserine metabolism; Porphyrin = porphyrin metabolism; Propionate = propionate metabolism; Purine = purine metabolism; Pyrimidine = pyrimidine metabolism; Pyruvate = pyruvate metabolism; Redox = redox metabolism; RedoxMitochondrialAntioxidantSystem = mitochondrial redox metabolism; Riboflavin = riboflavin (vitamin B2) metabolism; Selenocysteine = selenocysteine metabolism; ShikimateBiosynthesis = shikimate biosynthesis; SphingomyelinCeramide = sphingomyelin and ceramide metabolism; Terpenoid = terpenoid metabolism; Thiamine = thiamine biosynthesis; Transport = transport reactions; tRNA = tRNA and protein synthesis; Trp = tryptophan metabolism; Ubiquinone = ubiquinone metabolism; UtilizationPhospholipids = utilization of phospholipids; VitB6 = pyridoxal (vitamin B6) metabolism. (PDF 6 kb)
iPfal17 in SBML format. (XML 2000 kb)