Weighted Histogram Analysis Research Articles

Efficient Reconstruction of Complex Free Energy Landscapes by Multiple Walkers Metadynamics

Recently, we have introduced a new method, metadynamics, which is able to sample rarely occurring transitions and to reconstruct the free energy as a function of several variables with a controlled accuracy. This method has been successfully applied in many different fields, ranging from chemistry to biophysics and ligand docking and from material science to crystal structure prediction. We present an important development that speeds up metadynamics calculations by orders of magnitude and renders the algorithm much more robust. We use multiple interacting simulations, walkers, for exploring and reconstructing the same free energy surface. Each walker contributes to the history-dependent potential that, in metadynamics, is an estimate of the free energy. We show that the error on the reconstructed free energy does not depend on the number of walkers, leading to a fully linear scaling algorithm even on inexpensive loosely coupled clusters of PCs. In addition, we show that the accuracy and stability of the method are much improved by combining it with a weighted histogram analysis. We check the validity of our new method on a realistic application.

Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method: “Umbrella integration”

We present a method to analyze biased molecular-dynamics and Monte Carlo simulations, also known as umbrella sampling. In the limiting case of a strong bias, this method is equivalent to thermodynamic integration. It employs only quantities with easily controllable equilibration and greatly reduces the statistical errors compared to the standard weighted histogram analysis method. We show the success of our approach for two examples, one analytic function, and one biological system.

Through the channel and around the channel: Validating and comparing microscopic approaches for the evaluation of free energy profiles for ion penetration through ion channels.

Microscopic calculations of free energy profiles for ion transport through biological ion channels present a very serious challenge to modern simulation approaches. The main problem is due to the major convergence problems associated with the heterogeneous landscape of the electrostatic environment in ion channels and with the need to evaluate the profile associated with the transfer of the ion from bulk water to the channel environment. This problem is compounded by the lack of reliable and relevant benchmarks that can discriminate between alternative approaches. The present study is aimed at reducing the above problems by defining benchmarks that are directly relevant to ion channels and can also give converging results. This is done by constructing a series of models of a truncated gramicidin channel with different numbers of water molecules and by comparing the profiles for going around the channel and through the channel. These discriminating models are then used to validate and compare the adiabatic charging free energy perturbation (FEP) approach combined with an umbrella sampling approach (Warshel, A. J. Phys. Chem. 1982, 86, 2218) and the potential of mean force (PMF) approach used frequently in studies of ion channels. It is found that both approaches work quite well until one moves to the case of the fully solvated channel. In this limit, the PMF approach may give different results for the overall work of going through the channel and around the channel, while the FEP approach gives physically consistent results. The present benchmark also indicates that the weighted histogram analysis method (WHAM) approach does not offer a significant advantage over earlier approaches at least as much as studies of ion channels are concerned. Finally, it is concluded that the FEP approach may be more useful in evaluating the overall barrier for moving ions from water to ion channels and that in some cases it might be beneficial to use the FEP approach for selective points along the channel and then to connect these points by PMF calculations.

Conformational Free Energies of 1,2-Dichloroethane in Nanoconfined Methanol

Monte Carlo simulations have been used to construct free energy surfaces of 1,2-dichloroethane dissolved in methanol confined in hydrophobic spherical cavities of varying size (10-15 A) and solution density (0.6-0.79 g/cm3). The free energy surfaces are functions of two variables: the (center-of-mass) distance from the cavity wall of 1,2-dichloroethane and the Cl-C-C-Cl dihedral angle. Umbrella sampling and the weighted histogram analysis method were used to obtain accurate results for the free energy in these two degrees of freedom. Our results indicate that the conformational equilibrium and the barrier to internal rotation of the 1,2-dichloroethane depend on the position in the cavity. The results are discussed in the context of the solvent density, orientational distributions, and packing effects.

Conformational Equilibria and Free Energy Profiles for the Allosteric Transition of the Ribose-binding Protein

The ribose-binding protein (RBP) is a sugar-binding bacterial periplasmic protein whose function is associated with a large allosteric conformational change from an open to a closed conformation upon binding to ribose. The crystal structures of RBP in open and closed conformations have been solved. It has been hypothesized that the open and closed conformations exist in a dynamic equilibrium in solution, and that sugar binding shifts the population from open conformations to closed conformations. Here, we study by computer simulations the thermodynamic changes that accompany this conformational change, and model the structural changes that accompany the allosteric transition, using umbrella sampling molecular dynamics and the weighted histogram analysis method. The open state is comprised of a diverse ensemble of conformations; the open ribose-free X-ray crystal conformations being representative of this ensemble. The unligated open form of RBP is stabilized by conformational entropy. The simulations predict detectable populations of closed ribose-free conformations in solution. Additional interdomain hydrogen bonds stabilize this state. The predicted shift in equilibrium from the open to the closed state on binding to ribose is in agreement with experiments. This is driven by the energetic stabilization of the closed conformation due to ribose-protein interactions. We also observe a significant population of a hitherto unobserved ribose-bound partially open state. We believe that this state is the one that has been suggested to play a role in the transfer of ribose to the membrane-bound permease complex.

Absolute free energies of binding of peptide analogs to the HIV‐1 protease from molecular dynamics simulations

The constants of binding of five peptide analogs to the active site of the HIV-1 aspartic-protease are calculated based on a novel sampling scheme that is efficient and does not introduce any approximations in addition to the energy function used to describe the system. The results agree with experiments. The squared correlation coefficient of the calculated vs. the measured values is 0.79. The sampling scheme consists of a series of molecular dynamics integrations with biases. The biases are selected based on an estimate of the probability density function of the system in a way to explore the conformational space and to reduce the statistical error in the calculated binding constants. The molecular dynamics integrations are done with a vacuum potential using a short cutoff scheme. To estimate the probability density of the simulated system, the results of the molecular dynamics integrations are combined using an extension of the weighted histogram analysis method (C. Bartels, Chem. Phys. Letters 331 (2000) 446-454). The probability density of the solvated ligand-protein system is obtained by applying a correction for the use of the short cutoffs in the simulations and by taking into account solvation with an electrostatic term and a hydrophobic term. The electrostatic part of the solvation is determined by finite difference Poisson-Boltzmann calculations; the hydrophobic part of the solvation is set proportional to the solvent accessible surface area. Setting the hydrophobic surface tension parameter equal to 8 mol(-1) K(-1) A(-2), absolute binding constants are in the muM to nM range. This is in agreement with experiments. The standard errors determined from eight repeated binding constant determinations are a factor of 14 to 411. A single determination of a binding constant is done with 499700 steps of molecular dynamics integration and 4500 finite difference Poisson-Boltzmann calculations. The simulations can be analyzed with respect to conformational changes of the active site of the HIV-1 protease or the ligands upon binding and provide information that complements experiments and can be used in the drug development process.

Thermodynamics of folding and association of lattice-model proteins

Closely related to the "protein folding problem" is the issue of protein misfolding and aggregation. Protein aggregation has been associated with the pathologies of nearly 20 human diseases and presents serious difficulties during the manufacture of pharmaceutical proteins. Computational studies of multiprotein systems have recently emerged as a powerful complement to experimental efforts aimed at understanding the mechanisms of protein aggregation. We describe the thermodynamics of systems containing two lattice-model 64-mers. A parallel tempering algorithm abates problems associated with glassy systems and the weighted histogram analysis method improves statistical quality. The presence of a second chain has a substantial effect on single-chain conformational preferences. The melting temperature is substantially reduced, and the increase in the population of unfolded states is correlated with an increase in interactions between chains. The transition from two native chains to a non-native aggregate is entropically favorable. Non-native aggregates receive approximately 25% of their stabilizing energy from intraprotein contacts not found in the lowest-energy structure. Contact maps show that for non-native dimers, nearly 50% of the most probable interprotein contacts involve pairs of residues that form native contacts, suggesting that a domain-swapping mechanism is involved in self-association.

Comparison of efficiency and bias of free energies computed by exponential averaging, the Bennett acceptance ratio, and thermodynamic integration

Recent work has demonstrated the Bennett acceptance ratio method is the best asymptotically unbiased method for determining the equilibrium free energy between two end states given work distributions collected from either equilibrium and nonequilibrium data. However, it is still not clear what the practical advantage of this acceptance ratio method is over other common methods in atomistic simulations. In this study, we first review theoretical estimates of the bias and variance of exponential averaging (EXP), thermodynamic integration (TI), and the Bennett acceptance ratio (BAR). In the process, we present a new simple scheme for computing the variance and bias of many estimators, and demonstrate the connections between BAR and the weighted histogram analysis method. Next, a series of analytically solvable toy problems is examined to shed more light on the relative performance in terms of the bias and efficiency of these three methods. Interestingly, it is impossible to conclusively identify a "best" method for calculating the free energy, as each of the three methods performs more efficiently than the others in at least one situation examined in these toy problems. Finally, sample problems of the insertion/deletion of both a Lennard-Jones particle and a much larger molecule in TIP3P water are examined by these three methods. In all tests of atomistic systems, free energies obtained with BAR have significantly lower bias and smaller variance than when using EXP or TI, especially when the overlap in phase space between end states is small. For example, BAR can extract as much information from multiple fast, far-from-equilibrium simulations as from fewer simulations near equilibrium, which EXP cannot. Although TI and sometimes even EXP can be somewhat more efficient in idealized toy problems, in the realistic atomistic situations tested in this paper, BAR is significantly more efficient than all other methods.

Temperature Weighted Histogram Analysis Method, Replica Exchange, and Transition Paths

We analyzed the data from a replica exchange molecular dynamics simulation using the weighted histogram analysis method to combine data from all of the temperature replicas (T-WHAM) to obtain the room-temperature potential of mean force of the G-peptide (the C-terminal beta-hairpin of the B1 domain of protein G) in regions of conformational space not sampled at room temperature. We were able to determine the potential of mean force in the transition region between a minor alpha-helical population and the major beta-hairpin population and identify a possible transition path between them along which the peptide retains a significant amount of secondary structure. This observation provides new insights into a possible mechanism of formation of beta-sheet secondary structures in proteins. We developed a novel Bayesian statistical uncertainty estimation method for any quantity derived from WHAM and used it to validate the calculated potential of mean force. The feasibility of estimating regions of the potential of mean force with unfavorable free energy at room temperature by T-WHAM analysis of replica exchange simulations was further tested on a system that can be solved analytically and presented some of the same challenges found in more complex chemical systems.

Theoretical calculations of homoconjugation equilibrium constants in systems modeling acid–base interactions in side chains of biomolecules using the potential of mean force

The potentials of mean force (PMFs) were determined for systems forming cationic and anionic homocomplexes composed of acetic acid, phenol, isopropylamine, n-butylamine, imidazole, and 4(5)-methylimidazole, and their conjugated bases or acids, respectively, in three solvents with different polarity and hydrogen-bonding propensity: acetonitrile (AN), dimethyl sulfoxide (DMSO), and water (H(2)O). For each pair and each solvent a series of umbrella-sampling molecular dynamics simulations with the AMBER force field, explicit solvent, and counterions added to maintain a zero net charge of a system were carried out and the PMF was calculated by using the Weighted Histogram Analysis Method (WHAM). Subsequently, homoconjugation-equilibrium constants were calculated by numerical integration of the respective PMF profiles. In all cases but imidazole stable homocomplexes were found to form in solution, which was manifested as the presence of contact minima corresponding to hydrogen-bonded species in the PMF curves. The calculated homoconjugation constants were found to be greater for complexes with the OHO bridge (acetic acid and phenol) than with the NHN bridge and they were found to decrease with increasing polarity and hydrogen-bonding propensity of the solvent (i.e., in the series AN > DMSO > H(2)O), both facts being in agreement with the available experimental data. It was also found that interactions with counterions are manifested as the broadening of the contact minimum or appearance of additional minima in the PMF profiles of the acetic acid-acetate, phenol/phenolate system in acetonitrile, and the 4(5)-methylimidazole/4(5)-methylimidzole cation conjugated base system in dimethyl sulfoxide.

Long-Time Conformational Transitions of Alanine Dipeptide in Aqueous Solution: Continuous and Discrete-State Kinetic Models

We present an analysis of the thermodynamics, conformational dynamics, and kinetics of the solvated alanine dipeptide molecule. Solvation was treated in the framework of the OPLS/analytic generalized born with nonpolar interactions effective potential model. The effective free energy map was generated in a series of multiwindow umbrella sampling all-atom simulations using the weighted histogram analysis method. A Brownian dynamics approach was used to examine the room-temperature dynamics of the dipeptide. To emulate the realistic dynamics of conformational interconversions of the alanine dipeptide in water, the parameters that govern the time evolution of the reduced model were chosen to mimic the intermediate-time dynamics of a model which treats both the solute and the solvent explicitly. The characteristic time (mean first passage time) for the “fast” C7eq → αR transition was found to be around 249 ps in good agreement with the results of previously reported studies. From an ensemble of microsecond Brownian dynamics trajectories we were able to numerically estimate with high statistical precision the mean transition time for the “slow” transition αR → C7ax, which was found to be approximately 11 ns. Conformational kinetics of the dipeptide was further analyzed by introducing a discrete-state kinetic model consisting of four states (αR, β/C5/C7eq, αL, and C7ax) and 12 rate constants, which reproduces the long-time behavior of the Brownian dynamics simulations. The simulations described here complement and help to motivate single-molecule spectroscopic studies of alanine dipeptide and other peptides in aqueous solution. The close correspondence established between the continuous dynamics and a discrete-state kinetic model provides a possible route for studying the long-time kinetic behavior of larger polypeptides using detailed all-atom effective potentials.

Theoretical Calculations of Heteroconjugation Equilibrium Constants in Systems Modeling Acid−Base Interactions in Side Chains of Biomolecules Using the Potential of Mean Force

The potentials of mean force (PMFs) were determined for systems forming molecular and ionic heterocomplexes (in which the proton is transferred from a neutral acid to a neutral base molecule) composed of acetic acid and phenol with amines in three solvents with different polarities and hydrogen-bonding propensities: acetonitrile (AN), dimethyl sulfoxide (DMSO), and water. For each pair and each solvent, a series of umbrella-sampling molecular dynamics simulations with the AMBER force field and explicit solvent molecules were carried out and the PMF was calculated by using the weighted histogram analysis method (WHAM). Subsequently, heteroconjugation equilibrium constants were calculated by numerical integration of the respective PMF profiles. All complexes except those involving acetic acid as a proton donor in water were found to form in solution, which was manifested by the presence of contact minima corresponding to hydrogen-bonded species in the PMF curves. The calculated heteroconjugation constants were found to be greater for complexes with proton transfer than for those without proton transfer for the acetic acid-n-butylamine system in all solvents and for the phenol-imidazole and phenol-n-butylamine systems in DMSO when the same reference state (molecular acid and molecular base) was considered. For the acetic acid-n-butylamine system in acetonitrile and dimethyl sulfoxide, the calculated constants of ionic heterocomplex formation are greater than the values determined by potentiometric titration, while the calculated constants of molecular heterocomplex formation, as well as the calculated constants of heterocomplex formation from acetic acid and 4(5)-methylimidazole, are in good agreement with the experimental data, which suggests that the force field applied overestimates the energy of the interactions of oppositely charged ions or the energy of the interactions of hydrogen-bonding species with solvent molecules. The constants were found to decrease with an increasing polarity and hydrogen-bonding propensity of the solvent (i.e., in the series AN > DMSO > H 2 O), this being in agreement with the available experimental data.

Hybrid QM/MM study of thio effects in transphosphorylation reactions: the role of solvation.

Transphosphorylation thio effects in solution are studied using hybrid QM/MM calculations with a d-orbital semiempirical Hamiltonian. Activated dynamics simulations were performed for a 3' ribose-phosphate model in an explicit 20 A sphere of TIP3P water surrounded by a solvent boundary potential, and free energy analysis was performed using the weighted histogram analysis method. Single thio-substitutions at all of the phosphoryl oxygen positions and a double thio-substitution at the nonbridging positions were considered. The reaction free energy profiles are compared with available experimental data, and the role of solvation on the barrier heights and reaction coordinate is discussed. These results provide an important step in the characterization of thio effects in reactions of biological phosphates that may aid in the interpretation of kinetic data and ultimately help to unravel the catalytic mechanisms of ribozymes.

Carbohydrate conformation in aqueous solution: calculation of a QM/MM potential of mean force

The aqueous solution conformation of model disaccharide, 4- O-α- d-xylopyranosyl-α- d-xylopyranose, was investigated using a combined quantum mechanical (QM)/molecular mechanical (MM) potential. In vacuo and aqueous solution conformational free energy surfaces were constructed from potential of mean force calculations, using weighted histogram analysis of 8.5 and 13.5 ns combined QM/MM molecular dynamics simulations, respectively. Calculations suggest a proportion of direct (5%) and water-bridged (34%) intersaccharide hydrogen bonds in aqueous solution, the latter consistent with a broad range of ϕψ space. Mulliken analysis indicates non-uniform polarisation of carbohydrate electron density, contrasting with electrostatic descriptions employed in force fields based on uniform over-polarisation at the HF/6-31G* level of theory.

The Development of Replica-Exchange-Based Free-Energy Methods

The calculation of relative free energies that involve large reorganizations of the environment is one of the great challenges of condensed-phase simulation. Such calculations are of particular importance in protein−ligand free-energy calculations. To meet this challenge, we have developed new free-energy techniques that combine the advantages of the replica-exchange method with free-energy perturbation (FEP) and finite-difference thermodynamic integration (FDTI). These new techniques are tested and compared with FEP, FDTI, and the adaptive umbrella weighted histogram analysis method (AdUmWHAM) on the challenging calculation of the relative hydration free energy of methane and water. This calculation involves a large solvent configurational change. Through the use of replica-exchange moves along the λ-coordinate, the configurations sampled along λ are allowed to mix, which leads to dramatic improvements in solvent configurational sampling, an efficient reduction of random sampling error, and a reduction o...

Exploring Flory's isolated-pair hypothesis: statistical mechanics of helix-coil transitions in polyalanine and the C-peptide from RNase A.

To evaluate Flory's isolated-pair hypothesis in the context of helical peptides, we explore equilibrium conformations of alpha-helix-forming polypeptides as a function of temperature by means of replica exchange molecular dynamics in conjunction with the CHARMM/GB implicit solvent force field and the weighted histogram analysis method. From these simulations, Zimm-Bragg parameters, s and sigma, of Ac-Alan-NMe are computed as a function of temperature. The values obtained for s(T) and sigma(T) remain unchanged along the length of the polypeptide except for very short chains and yield results consistent with measurements based on short helix-forming peptides but suggest larger s values than anticipated from polymer-based measurements. From direct estimates of the density of states for Ac-Alan-NMe (n = 3-20) and peptide constructs based on the C peptide from RNase A, the conformational entropy is calculated versus temperature. The calculated S(T) shows a clear proportionality to the chain length over a wide range of temperature. This is observed in polypeptides with both significantly branched and simple methyl (alanine) side chains. These results provide evidence for the validity of Flory's isolated pair hypothesis, at least in the context of helical peptides and helix-to-coil transitions in these peptides.

Molecular Simulation Study of the Potentials of Mean Force for the Interactions between Models of Like-Charged and between Charged and Nonpolar Amino Acid Side Chains in Water

The purpose of this study was to determine the potentials of mean force (PMF) of the interactions between models of like-charged and between models of charged and nonpolar amino acid side chains in water to design improved side chain−side chain interaction potentials for our united residue UNRES force field for protein-structure prediction. Restrained molecular dynamics with the AMBER force field, the TIP3P model of water, and the Ewald summation were used to carry out simulations, and the weighted histogram analysis method (WHAM) was used to calculate the PMFs as functions of solute−solute distances. The following types of systems were considered to model the interactions between like-charged side chains and the interactions between charged and nonpolar side chains: (i) a pair of positively charged ions (potassium, ammonium, and guanidine, respectively); (ii) a pair of negatively charged ions (chloride and acetate, respectively), and (iii) pairs of methane with potassium, ammonium, and chloride ions, re...

The Filling Potential Method: A Method for Estimating the Free Energy Surface for Protein−Ligand Docking

We developed a molecular simulation method suitable for estimation of binding free energy, called the filling potential method, based on the concept of the Taboo search, which is a type of self-avoiding random walk consisting of a cycle of local-minimum searches and transition state searches. The filling potential method is an umbrella potential sampling method, and enables the ligand molecule to drift from its local minima automatically. In the case of the filling potential method, the umbrella potential is a combination of Gaussian-type repulsive potentials, which are located on the trajectory of the ligand. Without setting the reaction coordinates a priori, this method searches for and determines the suitable reaction coordinates by successive generation of umbrella potentials based on its trajectory analysis. The weighted histogram analysis for these trajectories gives the binding free energy of the ligand to the receptor protein. It was applied to the complex of thermolysin and its inhibitors in expl...

Implementation of an adaptive umbrella sampling method for the calculation of multidimensional potential of mean force of chemical reactions in solution.

We describe the implementation of an adaptive umbrella sampling method, making use of the weighted histogram analysis method, for computing multidimensional potential of mean force for chemical reaction in solution. The approach is illustrated by investigating the effect of aqueous solution on the free energy surface for the proton transfer reaction of [H(3)N-H-NH(3)](+) using a combined quantum mechanical and molecular mechanical AM1/TIP3P potential.

Free Energy Surface, Reaction Paths, and Kinetic Isotope Effect of Short-Chain Acyl-CoA Dehydrogenase

We have applied molecular dynamics umbrella-sampling simulations and ensemble-averaged variational transition state theory with multidimensional tunneling (EA-VTST/MT) to explore the free energy surface, reaction paths, and dideuterium kinetic isotope effect (KIE) of the β-oxidation of butyryl-coenzyme A by short-chain acyl-CoA dehydrogenase. The potential energy surface is obtained by combined quantum mechanics−molecular mechanics (QM/MM) with specific reaction parameters and a simple valence bond term. The calculations include determination of the potential of mean force (PMF) in both two dimensions (2D) and one dimension (1D) by using the weighted histogram analysis method. The 2D PMF indicates that the hydride transfer is the rate-limiting step of the mechanism, and the 1D PMFs are used to calculate rate constants. We include full molecular dynamics of a 7824-atom reaction zone with stochastic boundary conditions in a first approximation of the quasiclassical rate constant. This first approximation is then corrected by transmission coefficients that account for dynamical recrossing and tunneling. For the calculation of transmission coefficients, the system is divided into a 43-atom primary zone, consisting of atoms directly involved in the enzyme reaction, and a 23 234-atom secondary zone, called the bath. In both approximations considered, the primary zone atoms are allowed to have kinetic energy in all steps. For calculation of the transmission coefficient in the static-secondary-zone (SSZ) approximation, the bath is frozen for each member of an ensemble of transition state configurations and associated reaction paths, which are then averaged, whereas in the equilibrium-secondary-zone (ESZ) approximation 7781 atoms of the bath are fully relaxed (with stochastic boundary conditions) along each reaction path. We found that the SSZ simulations underestimate the KIE, whereas the ESZ simulations overestimate the KIE.