Metabolic host responses to malarial infection during the intraerythrocytic developmental cycle

Background The malarial parasite Plasmodium falciparum undergoes a complex life cycle, including an intraerythrocytic developmental cycle, during which it is metabolically dependent on the infected human red blood cell (RBC). To describe whole cell metabolic activity within both P. falciparum and RBCs during the asexual reproduction phase of the intraerythrocytic developmental cycle, we developed an integrated host-parasite metabolic modeling framework driven by time-dependent gene expression data. Results We validated the model by reproducing the experimentally determined 1) stage-specific production of biomass components and their precursors in the parasite and 2) metabolite concentration changes in the medium of P. falciparum-infected RBC cultures. The model allowed us to explore time- and strain-dependent P. falciparum metabolism and hypothesize how host cell metabolism alters in response to malarial infection. Specifically, the metabolic analysis showed that uninfected RBCs that coexist with infected cells in the same culture decrease their production of 2,3-bisphosphoglycerate, an oxygen-carrying regulator, reducing the ability of hemoglobin in these cells to release oxygen. Furthermore, in response to parasite-induced oxidative stress, infected RBCs downgraded their glycolytic flux by using the pentose phosphate pathway and secreting ribulose-5-phosphate. This mechanism links individually observed experimental phenomena, such as glycolytic inhibition and ribulose-5-phosphate secretion, to the oxidative stress response. Conclusions Although the metabolic model does not incorporate regulatory mechanisms per se, alterations in gene expression levels caused by regulatory mechanisms are manifested in the model as altered metabolic states. This provides the model the capability to capture complex multicellular host-pathogen metabolic interactions of the infected RBC culture. The system-level analysis revealed complex relationships such as how the parasite can reduce oxygen release in uninfected cells in the presence of infected RBCs as well as the role of different metabolic pathways involved in the oxidative stress response of infected RBCs. Electronic supplementary material The online version of this article (doi:10.1186/s12918-016-0291-2) contains supplementary material, which is available to authorized users.

Text S1. Calculation of metabolic fluxes in uninfected human red blood cells Our objective was to obtained metabolic fluxes in human red blood cells (RBCs) that were neither infected nor cocultured with the malaria parasite Plasmodium falciparum. We did this by calculating a set of fluxes that satisfied the constraints in a human erythrocyte metabolic network [1] and had the best fit to increasing/decreasing rates of extracellular metabolite concentrations, as estimated from experimental metabolomics data of uninfected RBC cultures [2].

METHODS
We calculated the fluxes of uninfected RBCs through the following three steps.

Step I: Estimating experimental rates of extracellular concentration changes
We first estimated the experimental means and 95% confidence intervals (CIs) of the relative change rates (in h -1 ) of the extracellular concentrations for each of the 51 metabolites that were present in the metabolic network [1] and whose extracellular metabolomic data were experimentally determined [2]. The data included the means and standard deviations derived from triplicate measurements of each metabolite's extracellular concentration at 8, 16, 24, 32, 40, and 48 h relative to the level at time t = 0 [2].
For each metabolite, we calculated the mean and 95% CI of its relative change rate generated by 10,000 random simulations. For each simulation, we sampled the metabolite's relative concentration at each time point from a normal distribution with the corresponding mean and standard deviation [2] and performed a linear regression of the sampled time-series concentrations to determine the relative change rate.
We further calculated the experimental means and 95% CI for the extracellular concentration change rate (in mmol/[h•10 12 RBC]) in the medium for 24 metabolites. Thus, we multiplied the corresponding values for the relative change rates (as obtained above) and the initial concentrations (in mmol/l) of the metabolites in the medium [3] and divided the products by the RBC concentration of 0.11 trillion/l (10 12 RBC/l) as calculated from a hematocrit level of 1% [2]. Although the lack of initial conditions for the other 27 metabolites prevented us from estimating their concentration change rates, for each of these metabolites, we qualitatively labeled the experimental concentration change direction as an increase (or a decrease) if the 95% CI of the metabolite's relative change rates included only positive (or negative) values. We were unable to determine the direction for a metabolite for which the corresponding 95% CI included both positive and negative values.
Step II: Determining the computational exchange rates Given that the metabolites' extracellular concentration changes stemmed from their uptake or secretion between RBCs and the medium, we calculated the computational exchange rates by 1) minimizing the differences between the computational rates and experimental means of the concentration changes for metabolites with available experimental means (the aforementioned 24 metabolites) and 2) minimizing the discrepancy in directions between the computational exchanges and experimental concentration changes for metabolites with available experimental observations (part of the 27 metabolites mentioned above). We did this by solving the following optimization problem: where L j comprises non-negative slack variables, RE denotes the set of the exchange reactions of the 24 metabolites for which we estimated the experimental means of concentration change rates in step I, v exp,j indicates the estimated mean rate for the metabolite with exchange reaction j and the corresponding v j denotes the computational exchange rate for which positive and negative values indicate secretion and uptake of the metabolite, RI and RD represent the set of exchange reactions of metabolites for which we only determined the experimental concentration change directions as increases and decreases, respectively, S RBC , lb, and ub indicate the stoichiometric matrix, lower bounds, and upper bounds, respectively, of the fluxes derived from the RBC metabolic network [1], and S' comprises the coefficients for additional constraints, which set the flux through the pentose phosphate pathway and that through the Rapoport-Luebering shunt to be 3.4% and 145%, respectively, of the glucose uptake rate [4].
Step III: Minimizing uptake rates and overall fluxes Due to the non-unique nature for the solution of the optimization problem in step II, we selected the metabolic flux distribution satisfying minimal uptake rates and minimal overall fluxes by sequentially solving two optimization problems. The first problem was as follows: where RU represents the set of uptake reactions of the metabolites for which we determined neither the rates nor directions of the experimental concentration changes and L j * was derived from the solution of the optimization in step II. α j comprises the concentrations in Roswell Park Memorial Institute (RPMI) 1640 medium [3] of the metabolites for uptake reaction j in RU. If the concentration for a metabolite was unavailable, we arbitrarily set it to a small value (10 -6 mmol/l).
Given the solution of the above equation, we finally calculated the fluxes of uninfected RBCs by minimizing the overall fluxes, as follows: where v j * comprises the fluxes through uptake reaction j in RU in the solution of the first optimization problem in step III.

RESULTS
First, we attempted to use the original iAB-RBC-283 metabolic network of RBCs [1] to reproduce the experimental concentration changes of the metabolites in the medium. We calculated the metabolic fluxes of uninfected RBCs, including the secretion or uptake fluxes of the metabolites present in the medium, and compared these computational exchange fluxes with the experimental concentration change rates for 24 metabolites for which the experimental rates were derived from metabolomic data [2] and initial concentration [3]. Table 1 shows that the computational rates were 0 for 18 metabolites, for 11 of which, such as arginine, the rates were within the experimental 95% CI, indicating quantitative agreements between computational results and experimental data for these metabolites. However, the disagreement of the remaining seven metabolites indicated that the RBC network was unable to capture the decrease in the concentration of R-pantothenate and the increases in thiamin, methionine, phenylalanine, (iso)leucine, tryptophan, and valine.
To address the above discrepancies, we further modified the iAB-RBC-283 metabolic network [1]. Due to the lack of the R-pantothenate consumption pathway in iAB-RBC-283, we added the pantothenate kinase reaction, which converts R-pantothenate into coenzyme A, in the global human metabolic network Recon1 [5]. Given the inability of human to synthesize thiamin (vitamin B 1 ), we added a hypothetical reaction of thiamin accumulation, by assuming that this metabolite's concentration increase in the medium was due to the secretion of molecules already contained in RBCs before the initial time point. We also added a degradation reaction of Albumax II (lipid-rich bovine serum albumin protein) and set the related coefficients of amino acids and lipids based on the protein's amino acid sequence [6] and experimentally measured lipid compositions [7], respectively. This was done based on the assumption that methionine, phenylalanine, (iso)leucine, tryptophan, and valine could be generated through the degradation of Albumax II protein in the medium [2] by proteolytic enzymes of RBCs [8,9], given the lack of de novo syntheses of these amino acids in humans.
Finally, given the modified metabolic network, we recalculated the metabolic fluxes and compared the computational and experimental changes in metabolite concentrations in the medium. Table 2 shows the quantitative agreements for most metabolites (20 of 24 metabolites).    Thus, instead of deriving the biomass production levels from only one simulation based on the given expression data, we created an ensemble of simulations using the experimental expression levels as a baseline, but with a small amount of added noise, and averaged the calculated results.
We added normally distributed random noise of zero mean and 10% standard deviation to the expression values of a gene g at time t, i.e., ĝ t = g t + f(μ,σ), where f(μ,σ) is a normal probability distribution with mean µ = 0 and standard deviation σ = 0.1·g t , g t is noise-free gene expression data, and ĝ t is the gene expression data after the noise addition. We repeated this procedure 20 times, calculated the metabolic fluxes during the IDC for each instance, and calculated a final averaged flux results for each metabolite. To ensure reproducibility of the results, we generated all random numbers using the same random generator algorithm [10] and recorded the seeds.
Additional file 1: Figure S2 (page 13 below) shows the results for the biomass production levels using this procedure. The effect is to minimize the hourly spikes, providing an overall smoothing that diminishes the exaggerated sensitivity of the model to gene expression variaions. This is noticeable in the reduction in variations associated with highly variable DNA and ubiquinol-8 production levels during the trophozoite stage and an overall smoothing out of the co-factor production levels. Thus, we can attribute the bulk of the observed high-frequency variability seen in these metabolite levels to an artificial high sensitive of the model to gene expression variations. The overall observations based on the averaged data discussed in the manuscript of general, low-frequency stage-specific trends and their differences are robust with respect to this model behavior. data for DNA and RNA [11] and for phospholipids [12]. The predictions of the new modeling framework closely correspond to our earlier developed model that did not explicitly treat host red blood cell metabolism [13], however, the results are not quantitatively the same. The horizontal bars indicate the length of the time intervals. The colors of these bars represent the simulation results (blue) and experimental data (green).