Quantitative time-course metabolomics in human red blood cells reveal the temperature dependence of human metabolic networks

The temperature dependence of biological processes has been studied at the levels of individual biochemical reactions and organism physiology (e.g. basal metabolic rates) but has not been examined at the metabolic network level. Here, we used a systems biology approach to characterize the temperature dependence of the human red blood cell (RBC) metabolic network between 4 and 37 °C through absolutely quantified exo- and endometabolomics data. We used an Arrhenius-type model (Q10) to describe how the rate of a biochemical process changes with every 10 °C change in temperature. Multivariate statistical analysis of the metabolomics data revealed that the same metabolic network-level trends previously reported for RBCs at 4 °C were conserved but accelerated with increasing temperature. We calculated a median Q10 coefficient of 2.89 ± 1.03, within the expected range of 2–3 for biological processes, for 48 individual metabolite concentrations. We then integrated these metabolomics measurements into a cell-scale metabolic model to study pathway usage, calculating a median Q10 coefficient of 2.73 ± 0.75 for 35 reaction fluxes. The relative fluxes through glycolysis and nucleotide metabolism pathways were consistent across the studied temperature range despite the non-uniform distributions of Q10 coefficients of individual metabolites and reaction fluxes. Together, these results indicate that the rate of change of network-level responses to temperature differences in RBC metabolism is consistent between 4 and 37 °C. More broadly, we provide a baseline characterization of a biochemical network given no transcriptional or translational regulation that can be used to explore the temperature dependence of metabolism.

The temperature dependence of biological processes has been studied at the levels of individual biochemical reactions and organism physiology (e.g. basal metabolic rates) but has not been examined at the metabolic network level. Here, we used a systems biology approach to characterize the temperature dependence of the human red blood cell (RBC) metabolic network between 4 and 37°C through absolutely quantified exo-and endometabolomics data. We used an Arrhenius-type model (Q 10 ) to describe how the rate of a biochemical process changes with every 10°C change in temperature. Multivariate statistical analysis of the metabolomics data revealed that the same metabolic network-level trends previously reported for RBCs at 4°C were conserved but accelerated with increasing temperature. We calculated a median Q 10 coefficient of 2.89 ؎ 1.03, within the expected range of 2-3 for biological processes, for 48 individual metabolite concentrations. We then integrated these metabolomics measurements into a cell-scale metabolic model to study pathway usage, calculating a median Q 10 coefficient of 2.73 ؎ 0.75 for 35 reaction fluxes. The relative fluxes through glycolysis and nucleotide metabolism pathways were consistent across the studied temperature range despite the non-uniform distributions of Q 10 coefficients of individual metabolites and reaction fluxes. Together, these results indicate that the rate of change of network-level responses to temperature differences in RBC metabolism is consistent between 4 and 37°C. More broadly, we provide a baseline characterization of a biochemical network given no transcriptional or translational regulation that can be used to explore the temperature dependence of metabolism.
The rate of biological processes increases with increasing temperature. The dependence of biochemical rates on temperature has been studied since the late 19th century using an Arrhenius-type approach (1)(2)(3)(4)(5)(6). The metric used for evaluating such temperature dependence is the temperature coefficient (Q 10 coefficient). 2 If Q 10 ϭ 2 for a given process, then the rate of that process increases by a factor of 2 for every 10°C increase in temperature. The Q 10 can be calculated from the slope of a rate versus temperature plot, which is approximately linear over the biologically relevant temperature range of 0 -40°C (7). Individual enzymes have different Q 10 coefficients that are generally expected to be in the 2-3 range (3,4,6). Temperature dependence at the physiological level is determined using phenomenological measurements (such as growth rate) to study overall physiological changes (5,(7)(8)(9)(10). Changes at the physiological level depend on more than changes in the underlying individual biochemical reaction rates (11,12). For instance, various regulatory mechanisms (e.g. transcriptional, post-translational, and allosteric) determine how cells respond to temperature shock (13). The existence of extra layers of regulation complicates the effects of temperature change on the biochemical network. Biology is inherently multiscale, and the gap between observing the temperature dependence at the scale of an individual reaction and at the physiological level can be addressed through methods of systems biology (14).
Red blood cells (RBCs) represent an ideal cell type to study the temperature dependence of network-level metabolic biochemistry because of the absence of a nucleus and genetic material. This absence results in a lack of complicated transcriptional or translational regulation on metabolic enzyme activity. Allosteric and other regulation of enzymatic reaction rates is still present in RBCs, representing enzyme kinetic mechanisms and thus direct biochemical functions.
In this study, we investigated the temperature dependence of metabolism in the RBC at the network level by examining the United States Department of Energy Grant DE-SC0008701 and NHLBI, National Institutes of Health Grant R43HL123074. The authors declare that they have no conflicts of interest with the contents of this article. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. This article contains supplemental information, Data S1, Figs. S1-S16, and Tables S1 and S2. 1 To whom correspondence should be addressed. Tel.: 8585345668; Fax: 8588223120; E-mail: palsson@ucsd.edu.
rate of change of metabolite concentrations and metabolic reactions rates. We measured exo-and endometabolomics profiles in human RBCs stored at four different temperatures that span the range between the ex vivo (storage) and in vivo (body) temperatures: 4, 13, 22, and 37°C. We regressed the concentration profile of each metabolite across the measured temperature range to calculate its Q 10 value. We integrated these measurements with a cell-scale network reconstruction of RBC metabolism (15) that contains 216 metabolites (43% of which are measurable by quantitative metabolomic profiling) to calculate Q 10 coefficients for reaction fluxes and to observe pathway usage. By examining metabolite profiles in the context of a cell-scale metabolic model, we were able to assess temperature dependence on a network level, thus bridging the gap between studies at the reaction and physiological levels.

Measurement of the temperature dependence of RBC metabolism ex vivo
RBCs were collected using standard collection procedures and stored in SAGM medium (16) at 4, 13, 22, and 37°C. Start-ing 1 day after RBC collection (taken as time 0 in figures), metabolomic measurements were made in biological triplicate over time at each temperature (Fig. 1A). The data set was composed of 97 metabolites. In addition, standard blood bank quality control and assurance measurements (e.g. hemolysis and pH) were made at multiple time points over 21 (4, 13, and 22°C) and 7 days (37°C); all measured profiles are presented in supplemental Figs. S6 -S16 and Data S1.
As part of the baseline characterization, we observed the same metabolic changes that have previously been reported in the literature for RBC storage at 4°C. We observed the same previously documented accumulation of lactate, 5-oxoproline, and hypoxanthine (17)(18)(19) and the depletion of AMP and the phosphoglycerate pool (17,19). We measured the same high levels of intracellular malate previously reported in SAGM at 4°C (19). A complete discussion of the raw metabolomics data is presented in the supplemental information and Figs. S2 and S3.
To determine temperature dependence, we first needed to determine how relative storage time scaled across temperature (Food and Drug Administration regulations set the maximum storage time for RBCs at 42 days). To determine a network-level . Overlaying these plots on the same axes shows that the shape of the three-phase metabolic decay is conserved but accelerated with increasing temperature as evidenced by the location of the day 7 time point. See supplemental Fig. S1 for a more detailed PCA plot at each temperature. The numbers in parentheses represent the amount of variance explained by each component. Black arrows and roman numerals label the three metabolic shifts that occur over the storage period. C, the first principal component was plotted against the time vector at each temperature to determine the relative storage time at each temperature. Linear regression was used to estimate the rate of change, showing strong correlation between PC1 and time at each temperature. D, these rates of change were then used to estimate the change in metabolic rate for every 10°C (Q 10 ) from an Arrhenius-type log 2 (rate) versus temperature plot.
Temperature dependence of human metabolic networks Q 10 , we used multivariate statistical analysis on the metabolomics data to assess the impact of temperature changes on systemic metabolism. The network-level Q 10 was used to determine the time points that represented the same metabolic phase at each temperature. Once these time periods were defined, we calculated Q 10 coefficients for each measured metabolite using linear regression. The metabolomics data were then integrated into a mechanistic cell-scale model to calculate the rate of each reaction (i.e. the flux) in the network at each temperature; these calculated reaction rates were used to determine Q 10 coefficients for reactions and to assess pathway usage across temperature.

Network-level temperature dependence
Following published reports for RBC storage at 4°C (17,19,20), principal component analysis (PCA) was performed to obtain a global characterization of the data set. PCA is a multivariate statistical method that reduces the dimensionality of a complex data set by calculating the relative contribution of each measurement to the overall variability observed in the data. For each temperature, we performed PCA (Fig. 1B) on the timeseries concentration profiles of eight recently identified extracellular metabolites (adenine, glucose, hypoxanthine, lactate, malate, nicotinamide, 5-oxoproline, and xanthine) that robustly represent the RBC metabolome under storage conditions (20). These metabolites serve as qualitative biomarkers for the age of stored RBCs and have been shown to also be good quantitative predictors for other systemic metabolite concentrations (21,22). To make an accurate comparison across temperatures, the same loading coefficients were applied to the data at each temperature (see "Experimental procedures" for full details). PCA revealed the same metabolic "shifts" that have been previously observed at 4°C (17,19). These shifts separate three distinct metabolic states that can be reliably determined from the profiles of the biomarkers (20). During storage at 4°C, the two shifts in the PCA plot occur approximately at days 10 and 17. Here, these same metabolic states were observed to be conserved but notably accelerated with temperature as evidenced by the location of the day 7 time point at each temperature ( Fig.  1B and supplemental Fig. S1). We identified the duration of the first metabolic state at each temperature. We then used these time points to determine the starting and ending points for the linear regression that would be used to calculate individual metabolite and reaction Q 10 coefficients.
The first principal component at each temperature was highly correlated with time ( Fig. 1C), yielding a "network-level" Q 10 of 1.46 (R 2 ϭ 0.97) that describes how the system proceeds in storage time. The third metabolic state at 4°C is primarily characterized by a general loss of function as the RBC undergoes severe morphological changes (23,24), often leading to complications for transfusion patients (25,26). Thus, we only used measurements from the first two metabolic states to calculate the network-level Q 10 . These results show an overall three-state metabolic decay that is observed to accelerate with increased temperature.

Metabolite-level temperature dependence
To determine the temperature dependence of individual metabolites, we used the data from the first metabolic state at each temperature (identified from the PCA results in Fig. 1B). We made this choice because the data are believed to be the most accurate as the cells are still intact and metabolism is functioning the closest to its normal physiological state. We linearly regressed the concentration profile of each metabolite at each temperature and used these rates to calculate a Q 10 value (Fig. 1, C and D). Not all metabolite profiles could be accurately fit with a linear curve during the first state; to account for this, we did not include metabolites with an R 2 Ͻ 0.50. Q 10 coefficients for the 48 metabolites whose profiles could be estimated well with a linear fit are reported in Table 1. The calculated Q 10 coefficients span 1.28 (extracellular 5-oxoproline) to 5.89 (intracellular hypoxanthine).
The calculated metabolite Q 10 coefficients generally fall in the expected range of 2-3 for biochemical reactions (3, 4, 6) with a median of 2.89 (Fig. 2). The standard deviation of the  (20).

Temperature dependence of human metabolic networks
metabolite Q 10 coefficients was 1.03, indicating that although the distribution was centered in the expected range for biological measurements the temperature scaling was not uniform across the network. Interestingly, the biomarker pools that have been shown to robustly define the metabolic decay process (20) represented almost the full range of calculated Q 10 coefficients from 1.28 (5-Oxoproline) to 5.05 (malate); extracellular xanthine and hypoxanthine did not have calculated Q 10 coefficients due to the R 2 cutoff. Several metabolite Q 10 coefficients fell below 2.00 or above 3.00, including extracellular 5-oxoporline (1.28), ATP (3.95), ADP (4.48), and extracellular malate (5.05).

Reaction-level temperature dependence
Next, we investigated the temperature dependence of biochemical reactions in the metabolic network. Previous studies have investigated the temperature dependence of individual reactions (1-6), but our goal was to use systems biology approaches to determine the temperature dependence of all reactions in the network together. To this end, we used a mechanistic cell-scale model of the RBC (15) to calculate the flux state of the network (i.e. the flux through each reaction in the system). The flux through a reaction (a rate with units of mmol/h) was calculated at each temperature; these values were then used to calculate a Q 10 for each reaction using the same procedure shown in Fig. 1D. We tailored the model to the physiology at each temperature by integrating the metabolomics measurements for the first metabolic state into the model according to Bordbar et al. (27). See "Experimental procedures" for full details on flux modeling and metabolomics integration.
Model simulations yielded flux states for each temperature. We excluded transporters and reactions that carried flux at fewer than three temperatures to ensure the accuracy of the Q 10 calculations. The calculated Q 10 coefficients for 35 reactions are shown in Table 2 (only fits with R 2 Ն 0.50 were included in analysis). The distribution of reaction Q 10 coefficients (Fig. 2) was tighter than that of the metabolite Q 10 coefficients (standard deviation of 0.75 for reactions versus 1.03 for metabolites) but was still centered in the expected 2-3 range (median of 2.73). Several reactions in nucleotide metabolism and glutathione metabolism had Q 10 coefficients above 3.5.

Pathway usage across temperature
We then examined pathway usage across the 33°C temperature range studied. Each reaction in RBC reconstruction has previously been assigned to one of 18 different metabolic subsystems (15). We first investigated whether the reactions in any of these pathways shared a similar dependence on temperature (i.e. pathways in which reactions had similar Q 10 coefficients). We performed a pairwise comparison of the Euclidean distance between Q 10 coefficients of reactions from a given subsystem (n reactions) versus 100,000 random permutations of n reactions from all calculated Q 10 coefficients. The reactions in glycolysis (n ϭ 12) and nucleotide metabolism (n ϭ 7) both showed statistically significantly smaller distances between Q 10 coefficients than would be expected due to random chance (p Ͻ 7eϪ5 and p Ͻ 0.047, respectively). Notably, these two metabolic subsystems contain reactions that interact with several of the storage age biomarkers (glucose, lactate, hypoxanthine, and xanthine). In particular, the temperature dependence of the glycolytic reactions is dictated by glucose, the major input to that linear pathway (Fig. 3).
Although the enrichment of a given pathway for similar Q 10 coefficients depends on the reaction fluxes, this analysis does not directly measure whether the flux states between two temperatures were similar. To investigate the conservation of pathway usage across temperature, we normalized the data at each temperature to glucose uptake and calculated the percent distance between the flux states (i.e. the flux through each reaction) at 13, 22, and 37°C and the flux state at 4°C. Reactions that were unused at every temperature were excluded (n ϭ 61). Of the remaining 166 reactions, 33 (19.9%) had less than 50% difference at each of the three higher temperatures. All of the glycolytic reactions for which Q 10 coefficients were calculated (Fig. 3) were present in this subset of reactions.

Organizational structure of the network
To this point, we have investigated how the metabolomics measurements could be used to determine the temperature dependence of individual metabolites and reactions and of the network. When we integrated the metabolomics measurements into the cell-scale model, the network structure at each temperature changed due to the accumulation and/or depletion of metabolites. These results, however, do not provide any explanation as to why we observe certain Q 10 coefficients for certain reactions. To answer this question, we studied the organizational structure of the network in terms of the coupling of certain reaction fluxes. Two reaction fluxes are said to be "fluxcoupled" if the ratio of one to the other is constant (28); the flux

Temperature dependence of human metabolic networks
coupling of a network is a property inherent to its topology. A set of coupled reactions is simply the linking of coupled reactions into a single pathway.
To determine whether the coupling of the network changed with temperature, we defined the coupling characteristic of a network to be the mean of the coupled reaction sets. We compared the base network structure of the RBC metabolic network with that of each temperature, revealing that the coupling characteristic of the network decreased ϳ5-10% at each temperature. This result prompted us to ask whether this decoupling was more or less than would be expected due to random chance. We ran a permutation test, determining that the decou-pling of the networks observed with a change in temperature was significantly less than would be expected due to random chance (p Ͻ 6eϪ3; see "Experimental procedures" for details on the generation of random networks and the permutation test). Thus, the RBC metabolic network is robust against changes in temperature over the measured 33°C range.

Discussion
The temperature dependence of biological processes is of fundamental interest. Q 10 coefficients have been used to characterize the temperature dependence of individual biochemical reactions and of organism-level behavior. Although these studies have yielded valuable insights, there is a gap between studying temperature dependence at the biochemical and physiological levels (1-10). Systems biology principles can be used to move past the individual reaction level and assess temperature dependence on the network level, effectively bridging the gap between the previous work on temperature dependence at the reaction and physiological levels. In this study, we used deepcoverage metabolomics of human red blood cells in storage to investigate the temperature dependence of network-level biochemistry. We chose a range of temperatures that ranges from the storage temperature of RBCs (4°C) to body temperature (37°C). By studying temperatures in this range, we provide baseline data with no transcriptional or translational regulation that can be used to begin to understand the temperature dependence of metabolism in broader contexts. Specifically, these data address the "passive" control that temperature has over network fluxes and metabolites in a system that has most of central carbon metabolism but is not growing. The results obtained here have several primary implications.
First, we determined that although the rate of change for each metabolite and reaction increased with increasing temperature the network-level response was dampened (i.e. lower Q 10 ). This behavior is to be expected because of the difference in response times between individual reactions and a pathway. Notably, the scaling of metabolite and reaction temperature dependence was not observed to be uniform across the network (with high variability in Q 10 coefficients). The fact that these network-level behaviors are conserved is a surprising result because certain enzyme inactivations might be expected to be qualitatively disruptive of network behavior at various temperatures. The observed behavior indicates that individual metabolite and reaction Q 10 coefficients vary dramatically to preserve global network characteristics. This variability was particularly notable for the storage age biomarkers (20), an interesting result because of their ability to define the qualitative trend of the entire network. Overall, such variability can be expected due to the thermodynamics and kinetics associated with biochemical and enzymatic reactions. Approximately half of the calculated temperature coefficients fell in the 2-3 range (Fig. 2), which is generally accepted as the typical estimate for biological systems (3,4,6). Several of the metabolites (e.g. malate, hypoxanthine, and glutathione) and reactions (several glutathione synthesis and nucleotide metabolism reactions) that we calculated to have high Q 10 coefficients are relatively disconnected from the rest of the network and thus may not be as affected.

Temperature dependence of human metabolic networks
Second, the pathway analyses indicate that even without transcriptional or translational regulation RBCs maintain consistent flux through glycolysis and nucleotide metabolism across temperature. The storage age biomarkers in glycolysis (glucose and lactate) represent the primary exchanges for the RBC metabolic network. The preservation of these exchanges indicates that sustaining these reactions is inherent to the network topology and necessary to maintain physiological functions. The change in pathway usage across the rest of the network is dependent upon kinetics and thermodynamics, properties that are subject to change over a 33°C temperature range.
Third, it is important to be aware that the temperature dependences calculated here are based on ex vivo measurements, not single-enzyme in vitro assays. Such an assay would inherently be absent from any regulatory or network influences. Our results, although also generated in the absence of transcriptional or translational regulatory effects, suggest that the organizational structure of the network influences the temperature dependence of individual enzymes. Thus, we would not expect the systemic Q 10 coefficients calculated here to correlate with previously reported in vitro values. For example, previous studies have reported Q 10 coefficients for pyruvate kinase (PYK) in the range of 3.3-4.2 for fish (5), 1.4 -1.9 for bats (29), and 1.66 -1.69 for turtles (30); we calculated a systemic Q 10 coefficient of 2.59 for PYK. The use of a cell-scale model to calculate the flux coupling of various reactions in the metabolic network provides an unambiguous explanation for why sets of reactions in the network have similar Q 10 coefficients despite variation of individual Q 10 coefficients. Additionally, the flux coupling results suggest that the ex vivo Q 10 coefficients calculated here would be different from in vitro Q 10 coefficients and from each other because the network structure intrinsically constrains the temperature dependence of certain reactions. The networks at the measured temperatures displayed no significant loss of flux coupling compared with the base model, suggesting that the structural characteristics of the RBC network are robust despite the accumulation or depletion of intermediate metabolites.

Temperature dependence of human metabolic networks
We have described the temperature dependence of a human metabolic network over a temperature range of 33°C, from 4°C (the Food and Drug Administration-defined RBC storage temperature for transfusion) to 37°C (the in vivo temperature). Although the RBC represents a simple cellular system, the lack of complex regulatory motifs allows for a direct interrogation of the systems biochemistry underlying metabolism. The use of systems biology methods empowered us to assess the temperature dependence of an ex vivo metabolic system at the network level, helping to understand the relationship between the structure and function of the RBC metabolic network.
All RBC samples were analyzed three times: once in positive ionization mode using acidic chromatographic condition and twice in negative ionization mode using both acidic and basic chromatographic conditions. During acidic conditions, mobile phase A was 100% ACN, and mobile phase B was 100% H 2 O, both containing 0.1% formic acid. The following elution gradient was used during acidic condition: 0 min 99% A, 7 min 30% A, 7.1 min 99% A, 10 min 99% A. Basic conditions used ACN, 10 mM sodium bicarbonate (95:5) as mobile phase A and ACN, 10 mM sodium bicarbonate (5:95) as mobile phase B. During basic conditions, the following elution gradient was used: 0 min 99% A, 6 min 30% A, 6.5 min 99% A, 10 min 99% A.
In all conditions, the flow rate was set at 0.4 ml/min, column temperature was set at 45°C, and injection volume was 3.5 l. The mass spectrometer was operated using a 1.5-kV capillary voltage, 30-V sampling cone, and 5-V extraction cone. The cone and the desolvation gas flow were 50 and 800 liters/h, respectively. The source and desolvation gas temperatures were 120 and 500°C, respectively. MS spectra were acquired in centroid mode from m/z 50 to 1,000 using a scan time of 0.3 s. Leucine enkephalin (2 ng/l) was used as lock mass (m/z 556.2771 and 554.2615 in positive and negative experiments, respectively).
Identification of unexpected metabolites was achieved by integration, alignment, and conversion of MS data points into exact mass retention time pairs (MarkerLynx, v4.1, Waters). The identity of the unexpected metabolites was established by verifying peak retention time, accurate mass measurements, and tandem mass spectrometry against our in-house database and online databases, including HMDB (32) and METLIN (33). TargetLynx (v4.1, Waters) was used to integrate chromatograms of targeted metabolites. Extracted ion chromatograms were extracted using a 0.02-mDa window centered on the expected m/z for each targeted compound. Quantitation was performed by external calibration with reference standards. Details regarding the quantitative analysis (including the linear range and limit of detection) are reported in supplemental Tables S1 and S2. See the supplemental information for additional details

Multivariate statistical analysis
PCA has previously shown three distinct metabolic shifts occurring over the 42-day storage period for RBCs at 4°C: days 1-10, days 10 -17, and days 17-42 (17,19,20). In the data presented in this study, RBCs were stored at 4°C for 21 days; to ensure that all three shifts were captured in full, we used the metabolomics data from Bordbar et al. (19) to calculate the weightings for the principal components. PCA was performed on the Z-scores of the eight extracellular biomarkers (20): adenine, glucose, hypoxanthine, lactate, malate, nicotinamide, 5-oxoproline, and xanthine. These weights were then used to transform the data presented here so that a more representative comparison could be made across all temperatures. The components were rotated in the transformed space such that the day 0 measurement appears in the lower left corner of the plot.

Calculation of temperature coefficients
Each measurement was plotted against time to find the rate of change. The rate of change was calculated through simple linear regression according to Equation 1, where ŷ is the calculated response, ␤ 0 is the y-intercept, ␤ 1 is the regression coefficient (i.e. the slope), and ⑀ is the error (Fig. 1C).
To determine the goodness of fit, the coefficient of determination (R 2 ) was calculated using Equation 2, where ŷ is the calculated value of y and y is the mean of y. Once a rate was obtained for each measurement, Q 10 was calculated from the slope of the log 2 (rate) versus temperature plot (Fig.  1D). This procedure is outlined in the following text. The definition of Q 10 is on the Arrhenius equation (34), which states that the rate of a process (k) is exponentially related to the temperature of that process, where k represents the rate of the process, A is a constant factor that represents the frequency of molecular collision, E a is the activation energy, R is the gas constant, and T is the temperature of the process (35). The temperature coefficient is empirically defined as the ratio of the rates of two processes (6, 11, 34, 36 -38),

Temperature dependence of human metabolic networks
where k 1 and k 2 represent the rates of two processes and T 1 and T 2 represent the temperature at which these processes occur. This relationship is the slope of a plot of log 2 (rate) versus temperature (Fig. 1D). Using linear regression, the slope of the log 2 (rate) versus temperature plots was calculated for each measurement and used to calculate the Q 10 value from Equation 5, where m is defined as follows.
Metabolites whose R 2 values were less than 0.50 were excluded in the analysis. This cutoff was determined based on the distribution of R 2 values (supplemental Fig. S5); our goal was to maximize the amount of data captured while simultaneously minimizing the inclusion of noisy or poorly fit data. We used the same procedure for calculating the reaction Q 10 coefficients. We only included reactions that carried flux in at least three of the temperatures; transport reactions and reactions whose R 2 was less than 0.50 were excluded in the analysis. This cutoff was determined based on the distribution of R 2 values (supplemental Fig. S5); our goal was to maximize the amount of data captured while simultaneously minimizing the inclusion of noisy or poorly fit data. For reactions, we extended this cutoff to be based on the p value of the F statistic for the linear regression fit of the log 2 (rate) versus temperature plot; only those reactions whose p value was less than 0.05 were included. The R 2 for the reactions whose p values were less than 0.05 was greater than 0.80 for all reactions (i.e. no reactions were excluded based on the R 2 value from its linearly regressed fit).

Flux modeling
We used a modified version of the erythrocyte metabolic reconstruction iAB-RBC-283 (15), which was previously used for building personalized kinetic models (39). We integrated the metabolite concentrations into this model to predict the flux state of the network at each temperature using non-steadystate flux balance analysis (27). 2,3-DPG has previously been determined to be one of the more important metabolites in RBC physiology but was not measured here. To obtain the most accurate model of the entire network, we used existing metabolomics data (19) to predict the concentration profile of 2,3-DPG using the profiles of the eight biomarkers as input according to the workflow described in Yurkovich et al. (21).

Flux coupling
We used F2C2 (40) to calculate the flux coupling of the metabolic networks. When the metabolomics data are integrated into the metabolic model, the structure of the network is altered through the addition of source and sink reactions (27). To construct random networks, we needed to ensure that the models were feasible and mimicked the models created using the measured data. Therefore, we randomly chose nodes in the base RBC model that were "measured" as either accumulation (source) or depletion (sink); we used the measured data to determine 1) the number of randomly "measured" nodes, 2) the distribution of intracellular versus extracellular nodes, 3) the distribution of sources versus sinks, and 4) the bounds for the added sources and sinks.
The permutation test compared the difference in the flux coupling characteristic between the base network and the networks at the measured temperatures with the flux coupling characteristic between the base model and 1,000 randomly generated networks (supplemental Fig. S4). The p value was determined to be the fraction of random networks whose difference in the flux coupling characteristic was less than or equal to that of the measured temperature networks.