The red blood cell's primary physiological function is to deliver oxygen that enters the body through the lungs, to tissues and to return carbon dioxide back to the lungs, to be expelled from the body. In order to carry out these respiratory functions, basic but critical metabolic functionality in the red cell must be maintained; including maintaining adequate ATP levels to maintain cell shape and the preservation of reduced glutathione to counteract oxidative stresses. These capabilities are directly affected to varying degrees in patients with metabolic enzyme deficiencies. The two most common of such deficiencies in the human red cell are in G6PD and PK [17, 18]. We considered dual perturbations in silico where genetic variations in the properties of these two enzymes are crossed with changing environmental conditions.
1. Baseline responses of the 'normal' red cell: elucidating the role of pooled variables and transporters
The 'normal' or 'nominal' set of parameters in the dynamic model of red cell metabolism (abbreviated here as the nRBC) formed a reference case from which parameters can be perturbed. The dynamic responses nRBC (Figure 1, top) were analyzed under conditions of increased redox and energy loads to provide a baseline, or the 'normal' response to such stresses. The nRBCunder increased redox vs. increased energy loads resulted in different flux states, notably shifting between reliance on the glycolytic pathway versus the pentose phosphate pathway. The computed steady state concentrations and fluxes are provided in Additional file 1.
The evaluation of the computed nRBC metabolic network responses and the resulting internal changes from the imposition of redox or energy loads led to the observation that changes in transport fluxes reflected changes in internal functional states (Additional file 2, Figure 1):
-
As the redox load increased, flux through glycolysis decreased and flux through the pentose phosphate pathway increased. At very high redox loads, the flux through phosphoglucoisomerase (PGI) even reversed direction to achieve maximal total flux through the oxidative branch of the pentose phosphate pathway. These expected shifts in pathway uses were accompanied by some unexpected changes in transport fluxes. There was a dramatic decrease in adenosine (ADO) uptake and a moderate increase in the adenosine deaminase (ADA) flux. Hypoxanthine (HX) excretion increased and inosine (INO) uptake decreased. Internally, the fluxes through adenine phosphoribosyltransferase (ADPRT) and phosphoribosylpyrophosphate synthetase (PRPPSYN) increased significantly.
-
As the energy load is increased, flux through the pentose phosphate pathway decreased and flux increased through glycolysis. The adenosine kinase (AK) flux decreased, with a relatively small decrease in AMP deaminase (AMPDA), adenine phosphoribosyltransferase (ADPRT), IMP nucleotidase (IMPase) and purine nucleoside phosphorylase (PNPase). The phosphoribosylpyrophosphate synthetase (PRPPSYN) flux decreased. These changes in the flux state change the state of transporters. HX excretion decreased and INO uptake increased slightly. There was a reduction in adenine (ADE) uptake and a parallel increase in adenosine (ADO) uptake.
Taken together, these in silico results showed that the ratio of the adenosine to adenine uptake reflected whether the cell was responding to an increased redox load or an energy load.
Interpretation
Red cell metabolic functions have previously been interpreted [18, 19] in terms of pooled concentration variables and ratios thereof; such as redox and phosphate capacities and potentials (Figure 1, middle). This approach led to a functional (or physiological) description of the state of red cell metabolism, rather than a detailed biochemical description. The computations of changes in fluxes and concentrations at varying redox and energy loads in the nRBC can be interpreted within this functional network description.
The alterations in the computed internal function state of the nRBC in response to environmental parameters can be reduced to changes in the states of transport fluxes (Figure 1, bottom). Increases in energy demands will shift the internal pathways to maintain the phosphate potential and ADO uptake will increase, while ADE uptake decreases and HX secretion decreases. Conversely, when increased oxidative stresses are experienced, the fluxes shift towards the redox potential with decrease in ADO uptake and an increase in ADE uptake and HX secretion.
These simulated responses provide the baseline case for the nRBC and they can be characterized in terms of a functional description of the state of nRBC metabolism (Figure 1). The key result here was that the environmental loads shift the use of pathways and the redox/energy pools. The state of these pools (their charge = occupancy/capacity [18, 19]) was adjusted by alterations in key transport fluxes that are accessible to measurement.
Dual perturbations in silico
Since the simulated flux uptake patterns reflected changes in the internal states of the nRBC to environmental perturbations, it suggested the possibility to characterize the responses different genetic variants (nRBC versus vRBC; 'v' for 'variant') in the same fashion. The responses of the nRBC were again used as a reference. This led to the notion of a dual perturbation experiment in silico as a way to analyze and design experiments to determine the effects of genetic variation. Differences in functional capabilities were assessed by comparisons between different states in two different environmental conditions for the nRBC vs. a vRBC (i.e. any one of the 6 G6PD variants), see Figure 2. Calculating the differences between the rows in Figure 2 and then comparing the nRBC to the vRBC resulted in a dual perturbation comparison, in which the one perturbation reflected an environmental change and the other results from a genetic mutation. Such in silico studies were carried out for G6PD variants and PK variants whose altered kinetic parameters were measured from patients [17–20].
The steady state flux distributions for the nRBC along with the six G6PD variants (vRBCi, i = 1, 6) were calculated at a high and a low redox load (see METHODS). The comprehensive set of steady state fluxes and concentrations calculated for the cell under different conditions led to the definition of a sensitivity parameter relating the intracellular redox state of the cell and extracellular transporters,
where Rhigh and Rlow reflect the redox state of the cell (R = NADPH/(NADP + NADPH)) under high and low redox loads, respectively, and ν is a corresponding transport flux. A logarithmic sensitivity measure was used since this proved to be a sensitive relationship. The computational results illustrated how the change in the redox state of the cell (summarized as the relative size of the reduced NADP pool), an intracellular quantity, was tracked by transport fluxes, an extracellular quantity (see Figure 3). Hence, changes in the internal states of cells were detected by measuring the exchange rates.
To further elaborate the relationship between the cofactor pools and transport fluxes in the cell, the differences between the high and low exchange fluxes were calculated for each variant and then normalized to the corresponding flux difference in the nRBC. A value of 1 indicates no difference compared to the nRBC, whereas values greater than or less than one reflect differences between the vRBCs from the nRBC. The computational results showed that the transport fluxes alone are enough to distinguish between the less severe (A+ and A-, vRBC1 and vRBC2, respectively) and more severe (Iwate, Niigata, Yamaguchi, and Portici, vRBC3 through vRBC6, respectively) SNP variants (Figure 4). This finding was interesting because it suggested that changes in transport fluxes could identify changes in internal functional states, differentiating between normal and pathophysiological situations, and that changes were identified using relative changes in the fluxes.
PK dual perturbation studies were carried out by applying different energy loads instead of redox loads (Figure 5). All of the vRBCs considered differed from the nRBC, but to varying degrees. Decreased lactate production by all of the variants reflected a decreased glycolytic flux due to the decreased activity of PK. The increased transport of the other metabolites was a result of the need to balance all carbon uptake and secretion. The Mantova variant was predicted to be less severe than the other vRBCs. In these cases, the altered transport fluxes reflected a decreased total pool size of ATP internally.
2. Dynamic response of vRBCs to redox loads: elucidation of altered time-scale hierarchy
In order to further characterize the differences that occur between the nRBC and the vRBCs, the dynamic states of the vRBCs were analyzed by comparing changes in the dynamic sequence of events that occur when they respond to a perturbation. Different biochemical inter-conversions, illustrated in Figure 1B, occured on different time scales [21] when the nRBC was exposed to an environmental perturbation. This sequence of events can be different in a vRBC and how much they differ in silico from the nRBC could be a measure of the pathological severity of the genetic variation.
Dynamic responses of the nRBC
The dynamic response of a network can be characterized through the comprehensive identification of the time-sequence of pool formation between the metabolites in a network [22]. This approach for determining the temporal structure in metabolic network responses was used as a basis for the analysis how it responds to perturbations. As the nRBC responded to an increased redox load, the shift in flux from glycolysis to the pentose pathway was also reflected in changes in the time scales for pooling of metabolites (see Additional file 2, Figure 1). The shifts in pooling time scales were observed throughout the network, however the most pronounced shifts involved the interactions between GL6P and the hexose-phosphates, the pentose phosphates and glutathione/NADPH, and the adenosine phosphates (AMP, ADP, ATP) with the nucleotide precursor and degradation products.
Metabolite pooling over progressive time scales involves a complex series of interactions. Pool formation on fast time scales generally reflect fast achievement of equilibrium between isomers and other metabolites whose steady state concentration ratios are close to their equilibrium ratios [25, 26]. Pooling on slower time scales involves interactions between aggregate pools between different pathways in the network. These interactions are illustrated in Figure 6, in order to highlight the distinction between pooling within pathways and interactions between pathways.
Responses of Glucose-6-Phosphate Dehydrogenase Deficient RBCs
The G6PD variants due to causal SNP mutations were analyzed in terms of their effect on the time progression of the pooling process. Modal decomposition [23] was carried out for each of these cases under low and high redox loads. The pooling between metabolites over time was then determined for the normal cell (Additional file 2, Figure 1) and compared to each of the genetic variants under high redox loads (see Additional file 2, Figures 2, 3, 4, 5, 6, 7) [22].
Changes in metabolite pooling over progressive time scales for the non-chronic hemolytic variants (A+ and A-; these patients exhibit a relatively mild form of hemolytic anemia [17]) were very similar to the normal cell (Additional file 2, Figures 1, 2, 3), with the significant differences in pooling between metabolites occurring between members of the pentose phosphate pathway and some of the nucleotide precursors. The chronic hemolytic anemia variants in contrast, exhibited more pronounced changes in the dynamic sequence of interactions. The largest differences involved; 1) interactions between the oxidative branch of the pentose phosphate pathway and the first half of glycolysis, and 2) interactions among members within the non-oxidative branch of the pentose phosphate pathway and the nucleotide salvage pathway metabolites. The alterations in dynamics of the oxidative branch in the pentose phosphate pathway resulted in disruptions of the interactions between the different pathways on slower time scales (Figure 6).
Considering network functions in the context of the functional pools described in the middle panel in Figure 1, a summary of the alterations in metabolic states under increased redox loads of the six genetic variants considered are shown in Figure 7. The areas of the boxes reflect the size of the respective metabolite pools at the steady state. It is immediately apparent that the chronic hemolytic variants have substantially reduced abilities to counteract oxidative stresses.
Additionally there were changes in the dynamics of the chronic hemolytic variants, primarily among the metabolites that contribute to the adenosine phosphate potential and the pentose phosphate potential. The last four variants exhibited altered dynamic interactions among members of the nucleotide salvage pathways and the non-oxidative branch of the pentose phosphate pathway. These changes are denoted by slanted arrows in Figure 7.
3. Dynamic response of the nRBC to energy loads
Similar computations of interactions among metabolites across time scales were applied to characterize the dynamic response of the nRBC to changes in energy loads (see Additional file 2, Figure 8 and METHODS). As discussed above, when the nRBC responds to an increased energy load, flux through glycolysis is increased at the cost of a reduced flux through the pentose pathway. When the dynamic pooling among metabolites for the nRBC at a normal versus an increased load were compared, there was broad shifting of pooling between metabolites, with significant changes occurring between glycolytic intermediates and adenosine metabolites, and between nucleotide precursors and pentose phosphates.
Pyruvate Kinase Deficient Cases
In order to characterize the time hierarchical functionality of PK genetic variants, the dynamic pooling correlations were calculated for each vRBC under high-energy load conditions. Each vRBC was then compared to the nRBC under high-energy load conditions (see Additional file 2, Figures 9–15). The salient observations from these computations were: 1) the vRBCs all exhibited comparatively similar patterns of variation from the responses of nRBC (unlike the G6PD variants considered above), and 2) the changes between high and low energy load in the variants were all similar to the changes observed for the nRBC. Since the energy charge remained fairly constant for all of the variants as well, the changes resulting from PK deficiency were reflected primarily as changes in the pool size of available ATP.