A database of thermodynamic properties of the reactions of glycolysis, the tricarboxylic acid cycle, and the pentose phosphate pathway

Abstract

A database of thermodynamic properties is developed, which extends a previous database of glycolysis and tricarboxylic acid cycle by adding the reactions of the pentose phosphate pathway. The raw data and documented estimations of solution properties are made electronically available. The database is determined by estimation of a set of parameters representing species-level free energies of formation. The resulting calculations provide thermodynamic and network-based estimates of thermodynamic properties for six reactions of the pentose phosphate pathway for which estimates are not available in the preexisting literature. Optimized results are made available in ThermoML format. Because calculations depend on estimated hydrogen and metal cation dissociation constants, an uncertainty and sensitivity analysis is performed, revealing 23 critical dissociation constants to which the computed thermodynamic properties are particularly sensitive.

Database URL:http://www.biocoda.org/thermo

Introduction

Reliable and self-consistent databases of thermodynamic properties for biochemical reactions are necessary for accurate analysis of biochemical systems (1–6). A recently developed database of thermodynamic properties for the reactions of glycolysis and the tricarboxylic acid cycle that was constituted from measured equilibrium data (7) represents a refinement to the Alberty database (8) in that it accounts for the ionic strength and interactions of biochemical reactants and metal cations (Mg²⁺, Ca²⁺, Na⁺ and K⁺) in estimating the derived properties from the raw data. The database of Li et al. (7) is a framework that can be extended and refined by adding the underlying raw experimental data in the database and/or refining the underlying model assumptions. Here an updated database is developed by adding the reactions of the pentose phosphate pathway into the original database.

As in Li et al. (7), thermodynamic properties (reference and ⁠) values are estimated by minimizing the difference between model predictions and experimental data. Apparent equilibrium constants for biochemical reactions are estimated by these derived thermodynamic properties and compared to experimental data measured under non-standard conditions. The basic formulae described in Li et al. (7) are used here, which account for temperature, ionic interactions and ion binding (8–13) effects to convert standard-state reference quantities to experimental state quantities.

Input data

Database of measured equilibrium constants

Raw experimental data are obtained from original reports (14–19). As in Li et al. (7), we preferentially select studies rated as A- and B-quality by Goldberg et al. (20–22). Free cation concentrations and ionic strength associated with these original studies are estimated based on the conditions reported in the original sources. All data and calculations are explicitly documented in the database. As we have previously established (7), each measurement entry in the raw-data database provides the following information: (i) enzyme name (EC number); (ii) experimental temperature, pH, ionic strength, apparent equilibrium constant, free metal cation concentrations, buffer and experimental method; (iii) quality rating from Goldberg et al. (20–22); (iv) notes on experiments and strategies of estimations and approximations applied in calculations; and (v) reference information. Up-to-date versions of the experimental database are made available at the URL http://www.biocoda.org/thermo, or by contacting the authors.

Database of reactions and estimated standard reaction enthalpies

There are eight reactions in the pentose phosphate pathway: glucose 6-phosphate dehydrogenase (EC 1.1.1.49, G6PD), 6-phosphogluconolactonase (EC 3.1.1.31, PGL), 6-phosphogluconate dehydrogenase (EC 1.1.1.44, PGD), ribose-5-phosphate isomerase (EC 5.3.1.6, R5PI), ribulose-phosphate 3-epimerase (EC 5.1.3.1, RUPE), transketolase (EC 2.2.1.1, TKL), transaldolase (EC 2.2.1.2, TAL) and transketolase 2 (EC 2.2.1.1, TKL2). The first three reactions are the oxidative pentose phosphate pathway, and the rest are the reductive pentose phosphate pathway (23). After adding these reactions into our original database, there are 33 reactions in total, which are shown in the Table 1 (first 33 reactions). Table 1 lists EC numbers, reaction names and the abbreviations employed, the reference reaction stoichiometries, and estimated standard reaction enthalpies at T = 298.15 K, I = 0 M. In the reference reaction, a superscript is used to indicate the charge of chemical species. For example, the abbreviation for references species for glucose (GLC⁰) is distinguished from the abbreviation for the biochemical reactant GLC.

Reaction enthalpies (⁠⁠) for two reactions in the pentose phosphate pathway (PGD and R5PI) can be estimated using van't Hoff equation because data on apparent equilibrium constants at different temperatures are available. Neither equilibrium data at different temperatures nor prior values of are available for the other six reactions of the pentose phosphate pathway; hence the symbol ‘#’ is used to denote the absence of data. For these cases, the value is set to zero in further calculations.

Database of reactant and dissociation constants

There are seven reactants introduced to the reactant database by adding the pentose phosphate pathway into the thermodynamic database: erythrose 4-phosphate (E4P), 6-phosphoglucono-δ-lactone (PGLT), 6-phospho-d-gluconate (PGN), ribose 5-phosphate (R5P), ribulose 5-phosphate (RU5P), sedoheptulose 7-phosphate (S7P) and xylulose 5-phosphate (X5P). For the five sugar phosphates (E4P, R5P, RU5P, S7P and X5P), cation dissociation constants for only R5P can be found in NIST database. Since E4P, RU5P, S7P and X5P are structurally similar substances to R5P (all have similar near neighbors to the phosphate group), a pragmatic approach to estimate the necessary dissociation constants is to use the values for R5P (24). Specifically, E4P, RU5P, S7P and X5P are assumed to have dissociation properties equal to those for R5P in our calculations. For PGLT, we use the pK_H1 value of 5.99 from Alberty (10). For PGN, Casazza et al. (15) report that the hydrogen ion dissociation constant for the carboxylic acid is >1.0E-4, while the constant for the phosphate ester is the same as that of glucose 6-P [pK_H1 = 5.99 (10)]. Here we arbitrarily assign the value of 4.995 to the pK_H1 for PGN, which is the average of 4 and 5.99. The uncertainties of these assignments are considered in the uncertainty and sensitivity analysis below.

Basic thermodynamic and ion binding data for biochemical reactants and associated reference species are listed in Table 2. Each entry in the table contains the following information: (i) detailed name of reactant; (ii) reference species abbreviation; (iii) reactant abbreviation; (iv) number of protons in reference species; and (v) the dissociation constants (provided as pK) and the corresponding dissociation enthalpies Dissociation constants and enthalpies are tabulated at 298.15 K and 0.1 M ionic strength. The symbol ‘#’ is used to indicate absence of data. For these cases, the pK's are assumed to be infinite (no binding) with corresponding dissociation constants equal to zero. In the calculations, all the pK and values are adjusted to a common reference state of T = 298.15 K and I = 0.

Estimation of standard-state thermodynamic quantities

Given the compiled raw experimental data on the reactions and reactants, the thermodynamic model is used to estimate reference and values for the reference reactions and species. Note that in our thermodynamic model the values are adjustable parameters estimated to obtain the best fit to the biochemical equilibrium data. As in previous studies (8), values of for oxidized species of certain redox pairs are arbitrarily set to zero. Thus, these values are not true free energies of the reactions of formation for these chemical species; instead, they are parameters for which the thermodynamic model makes optimal predictions for these interdependent biochemical reactions.

Because the calculated are interrelated [e.g. in the reductive pentose phosphate pathway, five carbon sugars (X5P and R5P) are converted into three carbon (GAP) and six carbon (F6P) sugars which can then be utilized by the pathway of glycolysis], the database of may be rigorously extended only by recalculating the entire database using all of the raw data. Therefore, model fitting is based on total 686 data entries for the network of all 33 reactions (the first 33 reactions listed in Table 1). Standard free energies of formation for all the reference species are unable to be estimated independently because there are 29 stoichiometrically independent reactions and in total 40 reactants in our system. There are four reactions PDH, IDH2, AKGDH and SDH (last four reactions in Tables 1 and 4) in the TCA cycle for which direct measurements are not available in the literature. For these reactions, values are calculated using Alberty's database (8) and set as constraints to perform the optimization of the thermodynamic model. IDH2 is not an independent reaction in the overall network of 33 reactions. Therefore, values of for 32 references species may be estimated from data on the 29 reactions with the three constraints. The values of for eight species are set to fixed values which are obtained from Alberty (8), as shown in Table 3 (values of for CoQ⁰ and CoQH₂⁰ do not come into these calculations). A constrained nonlinear optimization procedure with the fmincon solver (Mathworks, Inc.) is used to analyze the whole data set. By weighting in inverse proportion to the number of data points available for a given reaction and minimizing the difference between model predictions and experimental data, a simultaneous solution of standard reaction Gibbs energies is obtained for the entire data set.

Results

Estimated Gibbs free energies of reaction and formation

Figure 1 illustrates model predictions versus experimental data for all data used in the analysis. The predicted values are obtained based on accounting for the biochemical state associated with a given experimental measurement in the raw-data database. Data points in the pentose phosphate pathway are shown as filled squares. Other data points are shown as open squares. The data span almost 14 orders of magnitude, and reveal good agreement between model predictions and measured data.

The optimal estimates of listed in Table 4 (⁠⁠), compared to the values from Li et al. (7) (⁠⁠). of first 25 reactions are almost the same as the results obtained previously (7). For the remaining eight reactions of the pentose phosphate pathway, since the of the reference species for E4P, PGLT, PGN and S7P are not available in the Alberty (8) or Goldberg's database (25), the symbol ‘#’ is used to denote the lack of a value. For reactions R5PI and RUPE, the absolute difference between our current model predictions and the values reported by Goldberg et al. (25) is within a reasonable margin, 0.21 and 0.07 kJ/mol, respectively. Predicted versus experimental measures are plotted in Figure 2 for these two reactions (R5PI and RUPE). The current database and the Goldberg database yield similar results. Recall that the dissociation properties of R5P, RU5P and X5P are assumed to be the same; the binding polynomials for the reactants on the left- and right-hand sides of these reactions are identical. Thus, the computed apparent equilibrium constants for these reactions do not depend on the experimental ionic composition. The validity of this assumption will be considered in the uncertainty and sensitivity analysis below. In addition, the predicted thermodynamic properties are not affected by the ionic strength for these reactions.

Optimal predicted associated with the predictions are listed in Table 5, compared to the optimal values of our previous version of the database (7). After the seven reactants of the pentose phosphate pathway and three constraints (⁠ of reaction PDH, AKGDH and SDH) are introduced, most of the estimated values of are shifted substantially. However, from Table 4 we can see that these shifts do not change the predicted compared to previous predictions. This is because the shifts in do not change the optimized results (apparent or ⁠) of our model. Recall that here the estimated values represent parameters in a thermodynamic model for this set of interdependent biochemical reactions. Since the number of independent reactions for which data exist (29 reactions) is smaller than the number of values to represent the system (40 reactants), values for several reference species are set to either existing values or values reported elsewhere (Table 3). As a result, these values are not physical constants. Rather, they are parameters in a thermodynamic network model that together form a self-consistent picture of the thermodynamics of reactions of the set of reactants studied here.

Predicted apparent Gibbs free energies under physiology conditions

The contributions of non-standard elements to the differences between the standard free energy for the reference chemical reaction (⁠⁠) and the physiological energy values (⁠⁠) are listed in the last four columns of Table 4. The results can be divided to three cases: (i) if the stoichiometric coefficient of H⁺ is non-zero in the reference chemical reaction (v_H ≠ 0), the differences are substantial and pH contributes significantly, accounting for at least 73% of the difference, e.g. reactions GLK, PFK and GAP; (ii) if v_H = 0 (⁠⁠) and (⁠⁠), the differences tend to be relatively small compared with case (i) and temperature and/or binding polynomials contribute most to the differences, e.g. reactions PGI, TPI and ACON; (iii) if v_H = 0 (⁠⁠) and ionic strength and/or binding polynomials contribute most to the differences, e.g. reactions FBA, NPTH and PCL.

For the reactions of the pentose phosphate pathway, only the first three reactions G6PD, PGL and PGD show substantial difference between free energy for the reference chemical reaction (⁠⁠) and for biochemical reaction (⁠⁠) under physiological conditions. These three reactions have v_H ≠ 0. They are the oxidative portion of the pentose phosphate pathway (15, 23), which produce NADP_red and are essentially ‘irreversible’ in vivo (23).

Dissociation constants uncertainty and sensitivity analysis

The pK values listed in Table 2 are taken as the average value when there are several values (≥2) available in NIST database (27). For these pK values, the average value may not represent the best choice to be used in the model, i.e. some value among those available values may be more accurate than others. For some pK values, there exists only one estimate or no direct estimates. In order to predict the impact of uncertainty of these values on the model output, an uncertainty and sensitivity analysis is performed.

The product (US) combining uncertainty U and sensitivity S can be used to check the overlapping effect of uncertainty and sensitivity. For example, recall that we arbitrarily assign the value of 4.995 to the pK_H1 for PGN. Since the value is not available in NIST database, its uncertainty U is set to the average number 0.0609. If we consider the theoretical range of 4–5.99 discussed above, then the calculated uncertainty U is ∼0.4. For this case, because the computed pK_H1 sensitivity is 0.0225, the US product is <0.01, which is small enough that the value of pK_H1 for PGN has no substantial effect on the model output. Similarly, recall that E4P, RU5P, S7P and X5P are assumed to have dissociation properties equivalent to R5P in our calculations. Although the computed pK sensitivities of E4P, RU5P, S7P and X5P are 3.81E-6, 0.1007, 0.1433 and 0.1597, respectively, all the US products of the pK_H1 for these reactants are <0.01. Therefore, the assumption of equivalent dissociation properties for E4P, R5P, RU5P, S7P and X5P does not substantially effect our calculations.

Figure 3A illustrates that US products span eight orders of magnitude. Figure 3B illustrates the detailed distribution of the US products >0.01. All US products are <0.11. There are 23 cases for which US >0.01. These 23 US values belong to 15 reactants and four pKs as listed in Table 8. They are a subset of the pKs with sensitivity >0.1. They demonstrate the most important pKs which can make obvious impact on the model output. The largest four US values are as indicated in Figure 3B.

Database dissemination

ThermoML is an extensible markup language (XML)-based approach, which is an IUPAC standard for storage and exchange of thermodynamic property data (29–32). Our optimized results are stored in the standard ThermoML format with two small extensions to the current ThermoML schema (32): (i) adding ‘pseudo-Gibbs free energy of formation, kJ/mol’ in the list of ePropName in BioProperties of PureOrMixtureData; and (ii) adding ‘biochemical network calculation’ in ePredictionType of Prediction. We use the term ‘pseudo-Gibbs free energy of formation’ because in our thermodynamic model the estimated values represent adjustable parameters for the given set of interdependent biochemical reactions. In our ThermoML data files, the abbreviation name of reactant is also added in sCommonName in Compound. In the ThermoML reaction database, if the value of is not available, it is specified in sPredictionMethodDescription, and both nPropValue and nPropDigits are set to 0 in PropertyValue of NumValues. The ThermoML schema and our ThermoML data files are provided in Supplementary Data.

Discussion

We have updated our biochemical thermodynamic database by adding the reactions of the pentose phosphate pathway. To build the new database, each original publication has been studied on a case-by-case basis, case-specific assumptions and approximations determined and documented, and case-specific calculations performed and documented. Raw data and documented estimations of solution properties are made electronically available so that the updated database remains transparent and extensible.

The developed database is optimally self-consistent and consistent with the data available for reactions considered and the constraints. Theoretical predictions of apparent equilibrium constants optimally match experimental data on equilibrium constants. These apparent equilibrium constants are predicted based on the estimated species-level Gibbs free energies of formation and accounting for the effects of temperature, ionic interactions and hydrogen and metal cation binding. The new database provides thermodynamic and network-based estimates of thermodynamic properties for six reactions of the pentose phosphate pathway for which estimates are not available in the pre-existing literature.

These calculations demonstrate how network thermodynamic calculations are effectively extended. Adding raw experimental data on reactions and reactants of the pentose phosphate pathway into corresponding raw-data database, reaction database and reactant database, respectively, a previous thermodynamic network model was extended to include these elements.

Although the new optimal estimates of for reactants of glycolysis and tricarboxylic acid cycle are equal to previous estimates (7), most of the optimal predicted values which are associated with the predictions are shifted compared to the previous version. This result demonstrates that reliable and self-consistent extensions require the recalculation of entire set of species-level parameters. Furthermore, sensitivity analysis reveals that some formation energies can vary substantially without changing the optimized objective value significantly. While this set of optimal values represents a self-consistent set, combining these estimated formation energies with independently estimated parameters from other studies would require re-optimizing a combined database.

The uncertainty and sensitivity analysis of dissociation constants reveals 23 pKs most important to the model output. Additional experimental measurements of these parameter values are desirable.

Supplementary Data

Supplementary data are available at Database Online.

Acknowledgements

The authors are grateful to Robert Goldberg for advice and critical comments.

Funding

National Institutes of Health Heart Lung and Blood Institute (grant number HL072011). Funding for open access charge: National Institutes of Health Heart Lung and Blood Institute.

Conflict of interest. None declared.

References

Appendix A

Calculating the physiological free energy

Example

Nomenclature

B: an empirical constant taken to be 1.6 M^−1/2

E: the minimum squared difference between model simulations and experimental data

⁠: standard Gibbs free energy of reaction

⁠: standard Gibbs free energy of formation of  species i

⁠: apparent Gibbs free energy of reaction

⁠: enthalpy of proton/metal cation dissociation

⁠: standard enthalpy of reaction

I: ionic strength

K: chemical (reference reaction) equilibrium constant

⁠: apparent equilibrium constant

K_di: the dissociation constant of cation ion d

M: the number of reactions

N: the number of pKs which have several values available  in NIST database

N_j: the number of experimental data in each reaction

P_j: binding polynomial associated with reactant j

pK: negative logarithm of the dissociation constant

R: gas constant, 8.3145 J K⁻¹mol⁻¹

S: sensitivity

T: temperature

U: uncertainty

v_i: stoichiometric coefficient for species i in a given  reaction

z_i: valence of species i

⁠: the ionic strength contribution

⁠: the pH contribution

⁠: the binding polynomial contribution

⁠: the temperature contribution