Computing free energies using nested sampling-based approaches

Abstract

Nested sampling (NS) has emerged as a powerful statistical mechanical sampling technique to compute the partition function of atomic and molecular systems. From the partition function all thermodynamic quantities can be calculated in absolute terms, including absolute free energies and entropies. In this article, we provide a brief overview of NS within a Bayesian context, as well as overviews of how NS is used to compute the partition functions and thermodynamic quantities in the canonical and isothermal-isobaric ensembles. Then we introduce a new scheme, Coupling Parameter Path Nested Sampling, to estimate the free energy difference between two systems with different potential energy functions. The method uses a NS simulation to traverse the same path through phase space as would be covered in traditional coupling parameter-based methods such as thermodynamic integration and perturbation approaches. We demonstrate the new method with two case studies and confirm its accuracy by comparison to conventional methods, including Widom test particle insertion and thermodynamic integration. The proposed method provides a powerful alternative to traditional coupling parameter-based free energy simulation methods.     

From the sequence to the superstructural properties of DNAs

A theoretical model for predicting intrinsic and induced DNA superstructures as well as their thermodynamic properties is presented. Intrinsic sequence-dependent superstructures are evaluated by integrating local deviations from the canonical B-DNA of the different dinucleotide steps. Induced superstructures are obtained by adopting the principle of minimum deformation free energy, evaluated in the Fourier space, in the framework of first-order elasticity. Finally dinucleotide stacking energies and melting temperatures are considered to account for local flexibility. In fact the two scales are strongly correlated. The model works very satisfactorily in predicting the sequence-dependent effects on the DNA experimental behavior, such as the gel electrophoresis retardation, the writhe transitions in topologically constrained domains, the thermodynamic constants of circularization reactions as well as the nucleosome thermodynamic stability constants.     

A Statistical Thermodynamic Model for Investigating the Stability of DNA Sequences from Oligonucleotides to Genomes

We describe the development and testing of a simple statistical mechanics methodology for duplex DNA applicable to sequences of any composition and extensible to genomes. The microstates of a DNA sequence are modeled in terms of blocks of basepairs that are assumed to be fully closed (paired) or open. This approach generates an ensemble of bubblelike microstates that are used to calculate the corresponding partition function. The energies of the microstates are calculated as additive contributions from hydrogen bonding, basepair stacking, and solvation terms parameterized from a comprehensive series of molecular dynamics simulations including solvent and ions. Thermodynamic properties and nucleotide stability constants for DNA sequences follow directly from the partition function. The methodology was tested by comparing computed free energies per basepair with the experimental melting temperatures of 60 oligonucleotides, yielding a correlation coefficient of −0.96. The thermodynamic stability of genic/nongenic regions was tested in terms of nucleotide stability constants versus sequence for the Escherichia coli K-12 genome. It showed clear differentiation of the genes from promoters and captures genic regions with a sensitivity of 0.94. The statistical thermodynamic model presented here provides a seemingly new handle on the challenging problem of interpreting genomic sequences.     

A Statistical Thermodynamic Model for Investigating the Stability of DNA Sequences from Oligonucleotides to Genomes

We describe the development and testing of a simple statistical mechanics methodology for duplex DNA applicable to sequences of any composition and extensible to genomes. The microstates of a DNA sequence are modeled in terms of blocks of basepairs that are assumed to be fully closed (paired) or open. This approach generates an ensemble of bubblelike microstates that are used to calculate the corresponding partition function. The energies of the microstates are calculated as additive contributions from hydrogen bonding, basepair stacking, and solvation terms parameterized from a comprehensive series of molecular dynamics simulations including solvent and ions. Thermodynamic properties and nucleotide stability constants for DNA sequences follow directly from the partition function. The methodology was tested by comparing computed free energies per basepair with the experimental melting temperatures of 60 oligonucleotides, yielding a correlation coefficient of −0.96. The thermodynamic stability of genic/nongenic regions was tested in terms of nucleotide stability constants versus sequence for the Escherichia coli K-12 genome. It showed clear differentiation of the genes from promoters and captures genic regions with a sensitivity of 0.94. The statistical thermodynamic model presented here provides a seemingly new handle on the challenging problem of interpreting genomic sequences.     

Membrane hydrophobicity determines the activation free energy of passive lipid transport

The collective behavior of lipids with diverse chemical and physical features determines a membrane’s thermodynamic properties. Yet, the influence of lipid physicochemical properties on lipid dynamics, in particular interbilayer transport, remains underexplored. Here, we systematically investigate how the activation free energy of passive lipid transport depends on lipid chemistry and membrane phase. Through all-atom molecular dynamics simulations of 11 chemically distinct glycerophospholipids, we determine how lipid acyl chain length, unsaturation, and headgroup influence the free energy barriers for two elementary steps of lipid transport: lipid desorption, which is rate limiting, and lipid insertion into a membrane. Consistent with previous experimental measurements, we find that lipids with longer, saturated acyl chains have increased activation free energies compared to lipids with shorter, unsaturated chains. Lipids with different headgroups exhibit a range of activation free energies; however, no clear trend based solely on chemical structure can be identified, mirroring difficulties in the interpretation of previous experimental results. Compared to liquid-crystalline phase membranes, gel phase membranes exhibit substantially increased free energy barriers. Overall, we find that the activation free energy depends on a lipid’s local hydrophobic environment in a membrane and that the free energy barrier for lipid insertion depends on a membrane’s interfacial hydrophobicity. Both of these properties can be altered through changes in lipid acyl chain length, lipid headgroup, and membrane phase. Thus, the rate of lipid transport can be tuned through subtle changes in local membrane composition and order, suggesting an unappreciated role for nanoscale membrane domains in regulating cellular lipid dynamics.     

Transition modes in Ising networks: an approximate theory for macromolecular recognition.

For a statistical lattice, or Ising network, composed of N identical units existing in two possible states, 0 and 1, and interacting according to a given geometry, a set of values can be found for the mean free energy of the 0-->1 transition of a single unit. Each value defines a transition mode in an ensemble of nu N = 3N - 2N possible values and reflects the role played by intermediate states in shaping the energetics of the system as a whole. The distribution of transition modes has a number of intriguing properties. Some of them apply quite generally to any Ising network, regardless of its dimension, while others are specific for each interaction geometry and dimensional embedding and bear on fundamental aspects of analytical number theory. The landscape of transition modes encapsulates all of the important thermodynamic properties of the network. The free energy terms defining the partition function of the system can be derived from the modes by simple transformations. Classical mean-field expressions can be obtained from consideration of the properties of transition modes in a rather straightforward way. The results obtained in the analysis of the transition mode distributions have been used to develop an approximate treatment of the problem of macromolecular recognition. This phenomenon is modeled as a cooperative process that involves a number of recognition subsites across an interface generated by the binding of two macromolecular components. The distribution of allowed binding free energies for the system is shown to be a superposition of Gaussian terms with mean and variance determined a priori by the theory. Application to the analysis of the biologically interaction of thrombin with hirudin has provided some useful information on basic aspects of the interaction, such as the number of recognition subsites involved and the energy balance for binding and cooperative coupling among them. Our results agree quite well with information derived independently from analysis of the crystal structure of the thrombin-hirudin complex.     

Solvation effects and driving forces for protein thermodynamic and kinetic cooperativity: how adequate is native-centric topological modeling?

What energetic and solvation effects underlie the remarkable two-state thermodynamics and folding/unfolding kinetics of small single-domain proteins? To address this question, we investigate the folding and unfolding of a hierarchy of continuum Langevin dynamics models of chymotrypsin inhibitor 2. We find that residue-based additive Gō-like contact energies, although native-centric, are by themselves insufficient for protein-like calorimetric two-state cooperativity. Further native biases by local conformational preferences are necessary for protein-like thermodynamics. Kinetically, however, even models with both contact and local native-centric energies do not produce simple two-state chevron plots. Thus a model protein's thermodynamic cooperativity is not sufficient for simple two-state kinetics. The models tested appear to have increasing internal friction with increasing native stability, leading to chevron rollovers that typify kinetics that are commonly referred to as non-two-state. The free energy profiles of these models are found to be sensitive to the choice of native contacts and the presumed spatial ranges of the contact interactions. Motivated by explicit-water considerations, we explore recent treatments of solvent granularity that incorporate desolvation free energy barriers into effective implicit-solvent intraprotein interactions. This additional feature reduces both folding and unfolding rates vis-à-vis that of the corresponding models without desolvation barriers, but the kinetics remain non-two-state. Taken together, our observations suggest that interaction mechanisms more intricate than simple Gō-like constructs and pairwise additive solvation-like contributions are needed to rationalize some of the most basic generic protein properties. Therefore, as experimental constraints on protein chain models, requiring a consistent account of protein-like thermodynamic and kinetic cooperativity can be more stringent and productive for some applications than simply requiring a model heteropolymer to fold to a target structure.     

Linked-function origins of cooperativity in a symmetrical dimer

The thermodynamic origins of substrate binding cooperativity in a dimeric enzyme that can bind one substrate (A) and one allosteric ligand (X) to each of two identical subunits are discussed. It is assumed that maximal activity is not subject to allosteric modification and that the substrates and allosteric ligands achieve binding equilibrium in the steady state. Each uniquely ligated form is assumed to be capable of exhibiting unique binding properties, and only the principles of thermodynamic linkage are used to constrain the system further. The explicit relationship between the Hill coefficient, the concentration of X, and the magnitudes of the relevant coupling free energies and dissociation constants is derived. In the absence of X only the homotropic coupling between substrate sites contributes to a nonhyperbolic substrate saturation profile. An allosteric ligand, X, can alter the cooperativity in two distinct ways, one mechanism being manifested when X is saturating and the only only when X is present at saturating concentrations. By evaluating the concentration of substrate required to produce half-maximal velocity as a function of [X], as well as the Hill coefficients when X is absent and fully saturating, the dissociation and coupling constants most important for understanding the mechanisms of allosteric action in an enzyme of this type can be determined.     

Linkage between cooperative oxygenation and subunit assembly of cobaltous human hemoglobin.

The thermodynamic linkage between cooperative oxygenation and dimer-tetramer subunit assembly has been determined for cobaltous human hemoglobin in which iron(II) protoporphyrin IX is replaced by cobalt(II) protoporphyrin IX. The equilibrium parameters of the linkage system were determined by global nonlinear least-squares regression of oxygenation isotherms measured over a range of hemoglobin concentrations together with the deoxygenated dimer-tetramer assembly free energy determined independently from forward and reverse reaction rates. The total cooperative free energy of tetrameric cobalt hemoglobin (over all four binding steps) is found to be 1.84 (+/- 0.13) kcal, compared with the native ferrous hemoglobin value of 6.30 (+/- 0.14) kcal. Detailed investigation of stepwise cooperativity effects shows the following: (1) The largest change occurs at the first ligation step and is determined on model-independent grounds by knowledge of the intermediate subunit assembly free energies. (2) Cooperativity in the shape of the tetrameric isotherm occurs mainly during the middle two steps and is concomitant with the release of quaternary constraints. (3) Although evaluation of the pure tetrameric isotherm portrays identical binding affinity between the last two steps, this apparent noncooperativity is the result of a "hidden" oxygen affinity enhancement at the last step of 0.48 (+/- 0.12) kcal. This quaternary enhancement energy is revealed by the difference in subunit assembly free energies of the triply and fully ligated species and is manifested visually by the oxygenation isotherms at high versus low hemoglobin concentration. (4) Cobaltous hemoglobin dimers exhibit apparent anticooperativity of 0.49 (+/- 0.16) kcal (presumed to arise from heterogeneity of subunit affinities).(ABSTRACT TRUNCATED AT 250 WORDS)     

Thermodynamics of RNA-RNA binding

BACKGROUND: Reliable prediction of RNA-RNA binding energies is crucial, e.g. for the understanding on RNAi, microRNA-mRNA binding and antisense interactions. The thermodynamics of such RNA-RNA interactions can be understood as the sum of two energy contributions: (1) the energy necessary to 'open' the binding site and (2) the energy gained from hybridization. METHODS: We present an extension of the standard partition function approach to RNA secondary structures that computes the probabilities Pu[i, j] that a sequence interval [i, j] is unpaired. RESULTS: Comparison with experimental data shows that Pu[i, j] can be applied as a significant determinant of local target site accessibility for RNA interference (RNAi). Furthermore, these quantities can be used to rigorously determine binding free energies of short oligomers to large mRNA targets. The resource consumption is comparable with a single partition function computation for the large target molecule. We can show that RNAi efficiency correlates well with the binding energies of siRNAs to their respective mRNA target. AVAILABILITY: RNAup will be distributed as part of the Vienna RNA Package, www.tbi.univie.ac.at/~ivo/RNA/     

Electronic structure calculation study of metal complexes with a phytosiderophore mugineic acid

Structural and thermodynamic properties of biologically important metal-mugineic acid complexes have been studied from the theoretical side in order to understand the metal-chelating mechanism of phytosiderophore mugineic acid at an atomic level. Density-functional theory methods combined with the polarizable continuum model (PCM) have been employed to obtain free energies of complex formation and redox potentials for metal-mugineic acid complexes in solution. It has been found that the free energies of complex formation calculated at the B3LYP/PCM level of theory are in moderate agreement with available experimental results. The inclusion of explicit water molecules interacting with the carboxylic groups in deprotonated mugineic acid through strong hydrogen-bonds is found to further improve the calculated free energies of complex formation.     

Tuning hydrogen storage of carbon nanotubes by mechanical bending: theoretical study

The effects of mechanical bending on tuning the hydrogen storage of titanium functionalised (4,0) carbon nanotube have been assessed using density functional theory calculations with reference to the ultimate targets of the US Department of Energy (DOE). The assessment has been carried out in terms of physisorption, gravimetric capacity, projected densities of states, statistical thermodynamic stability and reaction kinetics. The Ti atom binds at the hollow site of the hexagonal ring. The average adsorption energies (?0.54 eV) per hydrogen molecule meet the DOE target for physisorption (?0.20 to ?0.60 eV). The curvature attributed to the bending angle has no effect on the average adsorption energies per H₂ molecule. With no metal clustering, the system gravimetric capacities are expected to be as large as 9.0 wt%. The reactions of the deformed (bent) carbon nanotube have higher probabilities of occurring than those of the un-deformed carbon nanotube. The Gibbs free energies, enthalpies and entropies meet the ultimate targets of the DOE for all temperatures and pressures. The closest reactions to zero free energy occur at (378.15 K/2.961 atm.) and reverse at (340 and 360 K/1 atm.). The translational component is found to exact a dominant effect on the total entropy change with temperature. Favourable kinetics of the reactions at the temperatures targeted by DOE are reported regardless of the applied pressure. The more preferable thermodynamic properties assigned to the bending nanotube imply that hydrogen storage can be improved compared to the nonbending nanotube.     

Energy landscapes and properties of biomolecules

Thermodynamic and dynamic properties of biomolecules can be calculated using a coarse-grained approach based upon sampling stationary points of the underlying potential energy surface. The superposition approximation provides an overall partition function as a sum of contributions from the local minima, and hence functions such as internal energy, entropy, free energy and the heat capacity. To obtain rates we must also sample transition states that link the local minima, and the discrete path sampling method provides a systematic means to achieve this goal. A coarse-grained picture is also helpful in locating the global minimum using the basin-hopping approach. Here we can exploit a fictitious dynamics between the basins of attraction of local minima, since the objective is to find the lowest minimum, rather than to reproduce the thermodynamics or dynamics.     

Ionic control of immobilized enzymes. Kinetics of acid phosphatase bound to plant cell walls.

When an enzyme is bound to an insoluble polyelectrolyte it may acquire novel kinetic properties generated by Donnan effects. It the enzyme is homogeneously distributed within the matrix, a variation of the electrostatic partition coefficient, when substrate concentration is varied, mimics either positive or negative co-operativity. This type of non-hyperbolic behaviour may be distinguished from true co-operativity by an analysis of the Hill plots. If the enzyme is heterogeneously distributed within the polyelectrolyte matrix, an apparent negative co-operativity occurs, even if the electrostatic partition coefficient does not vary when substrate concentration is varied in the bulk phase. If the partition coefficient varies, mixed positive and negative co-operativities may occur. All these effects must be suppressed by raising the ionic strength in the bulk phase. Attraction of cations by fixed negative charges of the polyanionic matrix may be associated with a significant decrease of the local pH. The magnitude of this effect is controlled by the pK of the fixed charges groups of the Donnan phase. The local pH cannot be much lower than the value of this pK. This effect may be considered as a regulatory device of the local pH. Acid phosphatase of sycamore (Acer pseudoplatanus) cell walls is a monomeric enzyme that displays classical Michaelis-Menten kinetics in free solution. However, when bound to small cell-wall fragments or to intact cells, it has an apparent negative co-operativity at low ionic strength. Moreover a slight increase of ionic strength apparently activates the bound enzymes and tends to suppress the apparent co-operativity. At I0.1, or higher, the bound enzyme has a kinetic behavior indistinguishable from that of the purified enzyme in free solution. These results are interpreted in the light of the Donnan theory. Owing to the repulsion of the substrate by the negative charges of cell-wall polygalacturonates, the local substrate concentration in the vicinity of the bound enzyme is smaller than the corresponding concentration in bulk solution. The kinetic results obtained are consistent with the view that there exist at least three populations of bound enzyme with different ionic environments: a first population with enzyme molecules not submitted to electrostatic effects, and two other populations with molecules differently submitted to these effects. The theory allows one to estimate the proportions of enzyme belonging to these populations, as well as the local pH values and the partition coefficients within the cell walls.     

An evaluation of Poisson-Boltzmann electrostatic free energy calculations through comparison with experimental mutagenesis data

For systems involving highly and oppositely charged proteins, electrostatic forces dominate association and contribute to biomolecular complex stability. Using experimental or theoretical alanine-scanning mutagenesis, it is possible to elucidate the contribution of individual ionizable amino acids to protein association. We evaluated our electrostatic free energy calculations by comparing calculated and experimental data for alanine mutants of five protein complexes. We calculated Poisson-Boltzmann electrostatic free energies based on a thermodynamic cycle, which incorporates association in a reference (Coulombic) and solvated (solution) state, as well as solvation effects. We observe that Coulombic and solvation free energy values correlate with experimental data in highly and oppositely charged systems, but not in systems comprised of similarly charged proteins. We also observe that correlation between solution and experimental free energies is dependent on dielectric coefficient selection for the protein interior. Free energy correlations improve as protein dielectric coefficient increases, suggesting that the protein interior experiences moderate dielectric screening, despite being shielded from solvent. We propose that higher dielectric coefficients may be necessary to more accurately predict protein-protein association. Additionally, our data suggest that Coulombic potential calculations alone may be sufficient to predict relative binding of protein mutants.     

A computer simulation of functional group contributions to free energy in water and a DPPC lipid bilayer

A series of all-atom molecular dynamics simulations has been performed to evaluate the contributions of various functional groups to the free energy of solvation in water and a dipalmitoylphospatidylcholine lipid bilayer membrane and to the free energies of solute transfer (Delta(DeltaG(o))X) from water into the ordered-chain interior of the bilayer. Free energies for mutations of the alpha-H atom in p-toluic acid to six different substituents (-CH3, -Cl, -OCH3, -CN, -OH, -COOH) were calculated by a combined thermodynamic integration and perturbation method and compared to literature results from vapor pressure measurements, partition coefficients, and membrane transport experiments. Convergence of the calculated free energies was indicated by substantial declines in standard deviations for the calculated free energies with increased simulation length, by the independence of the ensemble-averaged Boltzmann factors to simulation length, and the weak dependence of hysteresis effects on simulation length over two different simulation lengths and starting from different initial configurations. Calculated values of Delta(DeltaG(o))X correlate linearly with corresponding values obtained from lipid bilayer transport experiments with a slope of 1.1 and from measurements of partition coefficients between water and hexadecane or decadiene, with slopes of 1.1 and 0.9, respectively. Van der Waals interactions between the functional group of interest and the acyl chains in the ordered chain region account for more than 95% of the overall potential energy of interaction. These results support the view that the ordered chain region within the bilayer interior is the barrier domain for transport and that solvation interactions within this region resemble those occurring in a nonpolar hydrocarbon.     

Solvent denaturation and stabilization of globular proteins

Statistical thermodynamic theory has recently been developed to account for the stabilities of globular proteins. Here we extend that work to predict the effects of solvents on protein stability. Folding is assumed to be driven by solvophobic interactions and opposed by conformational entropy. The solvent dependence of the solvophobic interactions is taken from transfer experiments of Nozaki and Tanford on amino acids into aqueous solutions of urea or guanidine hydrochloride (GuHCl). On the basis of the assumption of two pathways involving collapse and formation of a core, the theory predicts that increasing denaturant should lead to a two-state denaturation transition (i.e., there is a stable state along each path separated by a free energy barrier). The denaturation midpoint is predicted to occur at higher concentrations of urea than of GuHCl. At neutral pH, the radius of the solvent-denatured state should be much smaller than for a random-flight chain and increase with either denaturant concentration or number of polar residues in the chain. A question of interest is whether free energies of folding should depend linearly on denaturant, as is often assumed. The free energy is predicted to be linear for urea but to have some small curvature for GuHCl. Predicted slopes and exposed areas of the unfolded states are found to be in generally good agreement with experiments. We also discuss stabilizing solvents and compare thermal with solvent denaturation.     

High and low density intracellular water.

The proposal that liquid water consists of microdomains of rapidly-exchanging polymorphs of high and low density is examined for its impact upon roles of water in biology. It is assumed that the two polymorphs persist in solution and adjacent to surfaces and that solutes partition asymmetrically between them. It transpires that chaotropes are solutes which partition preferentially into low density water and displace the water equilibrium toward the high density polymorph. Kosmotropes. both ionic and non-polar, partition into high density water and induce low density water. Displacement of the water equilibrium at constant temperature and pressure has a thermodynamic cost which can be high. This appears to be a dominant factor in folding of proteins and DNA, aggregation of biopolymers and insolubility of non-polar kosmotropes. Cells control both the concentration of proteins and the selection of small solutes to produce an intracellular environment most conducive to co-ordinated enzyme function. Intracellular water has similar microdomains to bulk water, but surfaces and solutes redistribute them. Average properties, as measured by NMR are similar, but local properties on a nm scale may differ widely. Enzymes apparently use these local differences to activate cations for transport, induce movement and for synthesis.