Development of thermodynamic optimum searching (TOS) to improve the prediction accuracy of flux balance analysis

Flux balance analysis (FBA) has been widely used in calculating steady‐state flux distributions that provide important information for metabolic engineering. Several thermodynamics‐based methods, for example, quantitative assignment of reaction directionality and energy balance analysis have been developed to improve the prediction accuracy of FBA. However, these methods can only generate a thermodynamically feasible range, rather than the most thermodynamically favorable solution. We therefore developed a novel optimization method termed as thermodynamic optimum searching (TOS) to calculate the thermodynamically optimal solution, based on the second law of thermodynamics, the minimum magnitude of the Gibbs free energy change and the maximum entropy production principle (MEPP). Then, TOS was applied to five physiological conditions of Escherichia coli to evaluate its effectiveness. The resulting prediction accuracy was found significantly improved (10.7–48.5%) by comparing with the ¹³C‐fluxome data, indicating that TOS can be considered an advanced calculation and prediction tool in metabolic engineering. Biotechnol. Bioeng. 2013; 110: 914–923. © 2012 Wiley Periodicals, Inc.     

Consistent Estimation of Gibbs Energy Using Component Contributions

Standard Gibbs energies of reactions are increasingly being used in metabolic modeling for applying thermodynamic constraints on reaction rates, metabolite concentrations and kinetic parameters. The increasing scope and diversity of metabolic models has led scientists to look for genome-scale solutions that can estimate the standard Gibbs energy of all the reactions in metabolism. Group contribution methods greatly increase coverage, albeit at the price of decreased precision. We present here a way to combine the estimations of group contribution with the more accurate reactant contributions by decomposing each reaction into two parts and applying one of the methods on each of them. This method gives priority to the reactant contributions over group contributions while guaranteeing that all estimations will be consistent, i.e. will not violate the first law of thermodynamics. We show that there is a significant increase in the accuracy of our estimations compared to standard group contribution. Specifically, our cross-validation results show an 80% reduction in the median absolute residual for reactions that can be derived by reactant contributions only. We provide the full framework and source code for deriving estimates of standard reaction Gibbs energy, as well as confidence intervals, and believe this will facilitate the wide use of thermodynamic data for a better understanding of metabolism.     

Thermodynamic Analysis of Energy Transfer in Acidogenic Cultures

A global thermodynamic analysis, normally used for pure cultures, has been performed for steady‐state data sets from acidogenic mixed cultures. This analysis is a combination of two different thermodynamic approaches, based on tabulated standard Gibbs energy of formation, global stoichiometry and medium compositions. It takes into account the energy transfer efficiency, ?, together with the Gibbs free energy dissipation, ΔG_o, analysis of the different data. The objective is to describe these systems thermodynamically without any heat measurement. The results show that ? is influenced by environmental conditions, where increasing hydraulic retention time increases its value all cases. The pH effect on ? is related to metabolic shifts and osmoregulation. Within the environmental conditions analyzed, ? ranges from 0.23 for a hydraulic retention time of 20 h and pH 4, to 0.42 for a hydraulic retention time of 8 h and a pH ranging from 7–8.5. The estimated values of ΔG_o are comparable to standard Gibbs energy of dissipation reported in the literature. For the data sets analyzed, ΔG_o ranges from –1210 kJ/mol_x, corresponding to a stirring velocity of 300 rpm, pH 6 and a hydraulic retention time of 6 h, to –20744 kJ/mol_x for pH 4 and a hydraulic retention time of 20 h. For average conclusions, the combined approach based on standard Gibbs energy of formation and global stoichiometry, used in this thermodynamic analysis, allows for the estimation of Gibbs energy dissipation values from the extracellular medium compositions in acidogenic mixed cultures. Such estimated values are comparable to the standard Gibbs energy dissipation values reported in the literature. It is demonstrated that ? is affected by the environmental conditions, i.e., stirring velocity, hydraulic retention time and pH. However, a relationship that relates this parameter to environmental conditions was not found and will be the focus of further research.     

Genome-scale thermodynamic analysis of Escherichia coli metabolism

Genome-scale metabolic models are an invaluable tool for analyzing metabolic systems as they provide a more complete picture of the processes of metabolism. We have constructed a genome-scale metabolic model of Escherichia coli based on the iJR904 model developed by the Palsson Laboratory at the University of California at San Diego. Group contribution methods were utilized to estimate the standard Gibbs free energy change of every reaction in the constructed model. Reactions in the model were classified based on the activity of the reactions during optimal growth on glucose in aerobic media. The most thermodynamically unfavorable reactions involved in the production of biomass in E. coli were identified as ATP phosphoribosyltransferase, ATP synthase, methylene-tetra-hydrofolate dehydrogenase, and tryptophanase. The effect of a knockout of these reactions on the production of biomass and the production of individual biomass precursors was analyzed. Changes in the distribution of fluxes in the cell after knockout of these unfavorable reactions were also studied. The methodologies and results discussed can be used to facilitate the refinement of the feasible ranges for cellular parameters such as species concentrations and reaction rate constants.     

A scalable algorithm to explore the Gibbs energy landscape of genome-scale metabolic networks

The integration of various types of genomic data into predictive models of biological networks is one of the main challenges currently faced by computational biology. Constraint-based models in particular play a key role in the attempt to obtain a quantitative understanding of cellular metabolism at genome scale. In essence, their goal is to frame the metabolic capabilities of an organism based on minimal assumptions that describe the steady states of the underlying reaction network via suitable stoichiometric constraints, specifically mass balance and energy balance (i.e. thermodynamic feasibility). The implementation of these requirements to generate viable configurations of reaction fluxes and/or to test given flux profiles for thermodynamic feasibility can however prove to be computationally intensive. We propose here a fast and scalable stoichiometry-based method to explore the Gibbs energy landscape of a biochemical network at steady state. The method is applied to the problem of reconstructing the Gibbs energy landscape underlying metabolic activity in the human red blood cell, and to that of identifying and removing thermodynamically infeasible reaction cycles in the Escherichia coli metabolic network (iAF1260). In the former case, we produce consistent predictions for chemical potentials (or log-concentrations) of intracellular metabolites; in the latter, we identify a restricted set of loops (23 in total) in the periplasmic and cytoplasmic core as the origin of thermodynamic infeasibility in a large sample (10(6)) of flux configurations generated randomly and compatibly with the prior information available on reaction reversibility.     

Thermodynamics-based metabolic flux analysis

A new form of metabolic flux analysis (MFA) called thermodynamics-based metabolic flux analysis (TMFA) is introduced with the capability of generating thermodynamically feasible flux and metabolite activity profiles on a genome scale. TMFA involves the use of a set of linear thermodynamic constraints in addition to the mass balance constraints typically used in MFA. TMFA produces flux distributions that do not contain any thermodynamically infeasible reactions or pathways, and it provides information about the free energy change of reactions and the range of metabolite activities in addition to reaction fluxes. TMFA is applied to study the thermodynamically feasible ranges for the fluxes and the Gibbs free energy change, Delta(r)G', of the reactions and the activities of the metabolites in the genome-scale metabolic model of Escherichia coli developed by Palsson and co-workers. In the TMFA of the genome scale model, the metabolite activities and reaction Delta(r)G' are able to achieve a wide range of values at optimal growth. The reaction dihydroorotase is identified as a possible thermodynamic bottleneck in E. coli metabolism with a Delta(r)G' constrained close to zero while numerous reactions are identified throughout metabolism for which Delta(r)G' is always highly negative regardless of metabolite concentrations. As it has been proposed previously, these reactions with exclusively negative Delta(r)G' might be candidates for cell regulation, and we find that a significant number of these reactions appear to be the first steps in the linear portion of numerous biosynthesis pathways. The thermodynamically feasible ranges for the concentration ratios ATP/ADP, NAD(P)/NAD(P)H, and H(extracellular)(+)/H(intracellular)(+) are also determined and found to encompass the values observed experimentally in every case. Further, we find that the NAD/NADH and NADP/NADPH ratios maintained in the cell are close to the minimum feasible ratio and maximum feasible ratio, respectively.     

MetRxn: a knowledgebase of metabolites and reactions spanning metabolic models and databases

Background  

Increasingly, metabolite and reaction information is organized in the form of genome-scale metabolic reconstructions that describe the reaction stoichiometry, directionality, and gene to protein to reaction associations. A key bottleneck in the pace of reconstruction of new, high-quality metabolic models is the inability to directly make use of metabolite/reaction information from biological databases or other models due to incompatibilities in content representation (i.e., metabolites with multiple names across databases and models), stoichiometric errors such as elemental or charge imbalances, and incomplete atomistic detail (e.g., use of generic R-group or non-explicit specification of stereo-specificity).     

Standard thermodynamic formation properties for the adenosine 5'-triphosphate series.

The criterion for chemical equilibrium at specified temperature, pressure, pH, concentration of free magnesium ion, and ionic strength is the transformed Gibbs energy, which can be calculated from the Gibbs energy. The apparent equilibrium constant (written in terms of the total concentrations of reactants like adenosine 5'-triphosphate, rather than in terms of species) yields the standard transformed Gibbs energy of reaction, and the effect of temperature on the apparent equilibrium constant at specified pressure, pH, concentration of free magnesium ion, and ionic strength yields the standard transformed enthalpy of reaction. From the apparent equilibrium constants and standard transformed enthalpies of reaction that have been measured in the adenosine 5'-triphosphate series and the dissociation constants of the weak acids and magnesium complexes involved, it is possible to calculate standard Gibbs energies of formation and standard enthalpies of formation of the species involved at zero ionic strength. This requires the convention that the standard Gibbs energy of formation and standard enthalpy of formation for adenosine in dilute aqueous solutions be set equal to zero. On the basis of this convention, standard transformed Gibbs energies of formation and standard transformed enthalpies of formation of adenosine 5'-trisphosphate, adenosine 5'-diphosphate, adenosine 5'-monophosphate, and adenosine at 298.15 K, 1 bar, pH = 7, a concentration of free magnesium ions of 10(-3) M, and an ionic strength of 0.25 M have been calculated.     

Equilibrium calculations on systems of biochemical reactions at specified pH and pMg.

The use of G' in discussing the thermodynamics of biochemical reactions at a specified pH and pMg is justified by use of a Legendre transform of the Gibbs energy G. When several enzymatic reactions occur simultaneously in a system, the standard transformed Gibbs energies of reaction delta rG'0 can be used in a computer program to calculate the equilibrium composition that minimizes the transformed Gibbs energy at the specified pH and pMg. The calculation of standard transformed Gibbs energies of formation of reactants at pH 7 and pMg 3 is described. In addition a method for calculating the equilibrium concentrations of reactants is illustrated for a system with steady state concentrations of some reactants like ATP and NAD.     

dGPredictor: Automated fragmentation method for metabolic reaction free energy prediction and de novo pathway design

Group contribution (GC) methods are conventionally used in thermodynamics analysis of metabolic pathways to estimate the standard Gibbs energy change (Δ_rG′^o) of enzymatic reactions from limited experimental measurements. However, these methods are limited by their dependence on manually curated groups and inability to capture stereochemical information, leading to low reaction coverage. Herein, we introduce an automated molecular fingerprint-based thermodynamic analysis tool called dGPredictor that enables the consideration of stereochemistry within metabolite structures and thus increases reaction coverage. dGPredictor has comparable prediction accuracy compared to existing GC methods and can capture Gibbs energy changes for isomerase and transferase reactions, which exhibit no overall group changes. We also demonstrate dGPredictor’s ability to predict the Gibbs energy change for novel reactions and seamless integration within de novo metabolic pathway design tools such as novoStoic for safeguarding against the inclusion of reaction steps with infeasible directionalities. To facilitate easy access to dGPredictor, we developed a graphical user interface to predict the standard Gibbs energy change for reactions at various pH and ionic strengths. The tool allows customized user input of known metabolites as KEGG IDs and novel metabolites as InChI strings (https://github.com/maranasgroup/dGPredictor).     

Reduction of Fe(III) (Hydr)oxides with Known Thermodynamic Stability by Geobacter metallireducens

Sulfate-reducing and methanogenic microorganisms become inactive when the concentration of the electron donors drops below a threshold set by the minimum Gibbs free energy required for the bacterial metabolism to be maintained. Thus, their activity is thermodynamically controlled. In this paper we study if the activity of dissimilatory Fe(III) reducing bacteria is also limited by the thermodynamics of the reaction. We synthesized five Fe (III) (hydr)oxides (FHOs) of moderate stability and determined the solubility product (log K _SO (?39.1)-(?41.8)), in order to calculate their standard free energy of formation. K _SO values, estimated from the particle size did not correspond with experimentally determined ones. HCO₃ ^? and PIPES-buffered media, containing 45 mM FHO and either 1, 10, or 100 mM acetate were inoculated with Geobacter metallireducens. At the end of bacterial reduction, the Gibbs free energy of the reaction showed significant differences between the different FHOs. The termination of the bacterial activity was consequently not triggered thermodynamically. However, the non-dissolved Fe(II) (HCl-soluble minus soluble Fe(II)) showed an excellent correlation with the surface of the FHOs (15 μmol m^?2). It is therefore likely that the termination of the reaction was caused by blocking of the FHO surface with insoluble Fe(II), as has been previously reported in the literature. The ecological significance of both thermodynamic limitation and surface availability limitation is discussed for FHOs of different K _SO in environments with approximately neutral pH.     

Estimation of standard Gibbs energy changes of biotransformations

Contributions and corrections for the estimation of standard Gibbs energies are given. The group contribution method, applicable to both cyclic and acyclic compounds, permits the approximate estimation of the standard Gibbs energy of a biotransformation, given the stoichiometry and structures of the metabolites involved. Estimated standard Gibbs energies of formation for a number of acyclic biochemical compounds are provided.     

Determination of the free-energy change for repair of a DNA phosphodiester bond

The repair of phosphodiester bonds in nicked DNA is catalyzed by DNA ligases. Ligation is coupled to cleavage of a phosphoanhydride bond in a nucleotide cofactor resulting in a thermodynamically favorable process. A free energy value for phosphodiester bond formation was calculated using the reversibility of the T4 DNA ligase reaction. The relative number of DNA nicks to phosphodiester bonds in a circular plasmid DNA, formed during this reaction at fixed concentrations of ATP to AMP and PP(i), was quantified. At 25 degrees C, pH 7, the equilibrium constant (K(eq)) for the ligation reaction is 3.89 x 10(4) m. This value corresponds to a standard free energy (DeltaG degrees ') of -6.3 kcal mol(-1). By subtracting the known energy contribution due to hydrolysis of ATP to AMP and PP(i), DeltaG degrees ' for the hydrolysis of a DNA phosphodiester bond is -5.3 kcal mol(-1).     

Thermal versus guanidine-induced unfolding of ubiquitin. An analysis in terms of the contributions from charge-charge interactions to protein stability.

We have characterized the guanidine-induced unfolding of both yeast and bovine ubiquitin at 25 degrees C and in the acidic pH range on the basis of fluorescence and circular dichroism measurements. Unfolding Gibbs energy changes calculated by linear extrapolation from high guanidine unfolding data are found to depend very weakly on pH. A simple explanation for this result involves the two following assumptions: (1) charged atoms of ionizable groups are exposed to the solvent in native ubiquitin (as supported by accessible surface area calculations), and Gibbs energy contributions associated with charge desolvation upon folding (a source of pK shifts) are small; (2) charge-charge interactions (another source of pK shifts upon folding) are screened out in concentrated guanidinium chloride solutions. We have also characterized the thermal unfolding of both proteins using differential scanning calorimetry. Unfolding Gibbs energy changes calculated from the calorimetric data do depend strongly on pH, a result that we attribute to the pH dependence of charge-charge interactions (not eliminated in the absence of guanidine). In fact, we find good agreement between the difference between the two series of experimental unfolding Gibbs energy changes (determined from high guanidine unfolding data by linear extrapolation and from thermal denaturation data in the absence of guanidine) and the theoretical estimates of the contribution from charge-charge interactions to the Gibbs energy change for ubiquitin unfolding obtained by using the solvent-accessibility-corrected Tanford-Kirkwood model, together with the Bashford-Karplus (reduced-set-of-sites) approximation. This contribution is found to be stabilizing at neutral pH, because most charged groups on the native protein interact mainly with groups of the opposite charge, a fact that, together with the absence of large charge-desolvation contributions, may explain the high stability of ubiquitin at neutral pH. In general, our analysis suggests the possibility of enhancing protein thermal stability by adequately redesigning the distribution of solvent-exposed, charged residues on the native protein surface.     

Levels of thermodynamic treatment of biochemical reaction systems.

Equilibrium calculations on biochemical reaction systems can be made at three levels. Level 1 is the usual chemical calculation with species at specified temperature and pressure using standard Gibbs energies of formation of species or equilibrium constants K. Level 2 utilizes reactants such as ATP (a sum of species) at specified T, P, pH, and pMg with standard transformed Gibbs energies of formation of reactants or apparent equilibrium constants K'. Calculations at this level can also be made on the enzymatic mechanism for a biochemical reaction. Level 3 utilizes reactants at specified T, P, pH, and pMg, but the equilibrium concentrations of certain reactants are also specified. The fundamental equation of thermodynamics is derived here for Level 3. Equilibrium calculations at this level use standard transformed Gibbs energies of formation of reactants at specified concentrations of certain reactants or apparent equilibrium constants K". Level 3 is useful in calculating equilibrium concentrations of reactants that can be reached in a living cell when some of the reactants are available at steady-state concentrations. Calculations at all three levels are facilitated by the use of conservation matrices and stoichiometric number matrices for systems. Three cases involving glucokinase, glucose-6-phosphatase, and ATPase are discussed.     

Quantitative Assessment of Thermodynamic Constraints on the Solution Space of Genome-Scale Metabolic Models

Constraint-based methods provide powerful computational techniques to allow understanding and prediction of cellular behavior. These methods rely on physiochemical constraints to eliminate infeasible behaviors from the space of available behaviors. One such constraint is thermodynamic feasibility, the requirement that intracellular flux distributions obey the laws of thermodynamics. The past decade has seen several constraint-based methods that interpret this constraint in different ways, including those that are limited to small networks, rely on predefined reaction directions, and/or neglect the relationship between reaction free energies and metabolite concentrations. In this work, we utilize one such approach, thermodynamics-based metabolic flux analysis (TMFA), to make genome-scale, quantitative predictions about metabolite concentrations and reaction free energies in the absence of prior knowledge of reaction directions, while accounting for uncertainties in thermodynamic estimates. We applied TMFA to a genome-scale network reconstruction of Escherichia coli and examined the effect of thermodynamic constraints on the flux space. We also assessed the predictive performance of TMFA against gene essentiality and quantitative metabolomics data, under both aerobic and anaerobic, and optimal and suboptimal growth conditions. Based on these results, we propose that TMFA is a useful tool for validating phenotypes and generating hypotheses, and that additional types of data and constraints can improve predictions of metabolite concentrations.     

Thermodynamics-based Metabolite Sensitivity Analysis in metabolic networks

The increasing availability of large metabolomics datasets enhances the need for computational methodologies that can organize the data in a way that can lead to the inference of meaningful relationships. Knowledge of the metabolic state of a cell and how it responds to various stimuli and extracellular conditions can offer significant insight in the regulatory functions and how to manipulate them. Constraint based methods, such as Flux Balance Analysis (FBA) and Thermodynamics-based flux analysis (TFA), are commonly used to estimate the flow of metabolites through genome-wide metabolic networks, making it possible to identify the ranges of flux values that are consistent with the studied physiological and thermodynamic conditions. However, unless key intracellular fluxes and metabolite concentrations are known, constraint-based models lead to underdetermined problem formulations. This lack of information propagates as uncertainty in the estimation of fluxes and basic reaction properties such as the determination of reaction directionalities. Therefore, knowledge of which metabolites, if measured, would contribute the most to reducing this uncertainty can significantly improve our ability to define the internal state of the cell. In the present work we combine constraint based modeling, Design of Experiments (DoE) and Global Sensitivity Analysis (GSA) into the Thermodynamics-based Metabolite Sensitivity Analysis (TMSA) method. TMSA ranks metabolites comprising a metabolic network based on their ability to constrain the gamut of possible solutions to a limited, thermodynamically consistent set of internal states. TMSA is modular and can be applied to a single reaction, a metabolic pathway or an entire metabolic network. This is, to our knowledge, the first attempt to use metabolic modeling in order to provide a significance ranking of metabolites to guide experimental measurements.     

Systems of biochemical reactions from the point of view of a semigrand partition function

Semigrand partition functions contain all the thermodynamic information on reaction systems. When they are written for systems at specified pH, they yield the transformed Gibbs energy G' of the system and the thermodynamic properties that can be calculated from G'. When they are written for systems at specified pH and specified concentrations of coenzymes, they yield the further transformed Gibbs energy G" and properties that can be calculated from G". This is illustrated by considering: (1) a reactant that is a weak monoprotic acid at a specified pH; (2) a reaction between two pseudoisomer groups at a specified pH; and (3) the first five reactions of glycolysis. Equilibrium compositions in glycolysis are calculated at pH 7 and different steady-state concentrations of ATP and ADP.     

Reductive unfolding of serum albumins uncovered by Raman spectroscopy

The reductive unfolding of bovine serum albumin (BSA) and human serum albumin (HSA) induced by dithiothreitol (DTT) is investigated using Raman spectroscopy. The resolution of the S-S Raman band into both protein and oxidized DTT contributions provides a reliable basis for directly monitoring the S-S bridge exchange reaction. The related changes in the protein secondary structure are identified by analyzing the protein amide I Raman band. For the reduction of one S-S bridge of BSA, a mean Gibbs free energy of -7 kJ mol(-1) is derived by studying the reaction equilibrium. The corresponding value for the HSA S-S bridge reduction is -2 kJ mol(-1). The reaction kinetics observed via the S-S or amide I Raman bands are identical giving a reaction rate constant of (1.02 +/- 0.11) M(-1) s(-1) for BSA. The contribution of the conformational Gibbs free energy to the overall Gibbs free energy of reaction is further estimated by combining experimental data with ab initio calculations.     

Comparative thermodynamic analysis of DNA--protein interactions using surface plasmon resonance and fluorescence correlation spectroscopy

We report a kinetic and thermodynamic analysis of interactions between ssDNA and replication protein A (RPA) using surface plasmon resonance (SPR) and fluorescence correlation spectroscopy (FCS) at variable temperature. The two methods yield different values for the Gibbs free energy but nearly the same value for the reaction enthalpy of ssDNA-RPA complex formation. The Gibbs free energy was determined by SPR and FCS to be -62.6 and -54.7 kJ/mol, respectively. The values for the reaction enthalpy are -64.4 and -66.5 kJ/mol. It is concluded that the difference in Gibbs free energy measured by the two methods is due to different reaction entropies. The entropic contribution to the free energy at 25 degrees C is -1.8 kJ/mol for SPR and -11.8 kJ/mol for FCS. In SPR, the reaction is restricted to two dimensions because of immobilization of the DNA molecules to the sensor surface. In contrast, FCS is able to follow complex formation without spatial restrictions. In consequence, the reaction entropy determined from SPR experiments is lower than for FCS experiments.