Calibration of force-field dependency in free energy landscapes of peptide conformations by quantum chemical calculations

The free energy landscapes of peptide conformations were calibrated by ab initio quantum chemical calculations, after the enhanced conformational diversity search using the multicanonical molecular dynamics simulations. Three different potentials of mean force for an isolated dipeptide were individually obtained by the multicanonical molecular dynamics simulations using the conventional force fields, AMBER parm94, AMBER parm96, and CHARMm22. Each potential of mean force was then calibrated based upon the umbrella sampling algorithm from the adiabatic energy map that was calculated separately by the ab initio molecular orbital method, and all of the calibrated potentials of mean force coincided well. The calibration method was also applied to the simulations of a peptide dimer in explicit water models, and it was shown that the calibrated free energy landscapes did not depend on the force field used in the classical simulations, as far as the conformational space was sampled well. The current calibration method fuses the classical free energy calculation with the quantum chemical calculation, and it should generally make simulations for biomolecular systems much more reliable when combining with enhanced conformational sampling.     

Multidimensional free-energy calculations using the weighted histogram analysis method

The recently formulated weighted histogram analysis method (WHAM)¹ is an extension of Ferrenberg and Swendsen's multiple histogram technique for free-energy and potential of mean force calculations. As an illustration of the method, we have calculated the two-dimensional potential of mean force surface of the dihedrals gamma and chi in deoxyadenosine with Monte Carlo simulations using the all-atom and united-atom representation of the AMBER force fields. This also demonstrates one of the major advantages of WHAM over umbrella sampling techniques. The method also provides an analysis of the statistical accuracy of the potential of mean force as well as a guide to the most efficient use of additional simulations to minimize errors. © 1995 John Wiley & Sons, Inc.     

A united-residue force field for off-lattice protein-structure simulations. I. Functional forms and parameters of long-range side-chain interaction potentials from protein crystal data

A two-stage procedure for the determination of a united-residue potential designed for protein simulations is outlined. In the first stage, the long-range and local-interaction energy terms of the total energy of a polypeptide chain are determined by analyzing protein-crystal data and averaging the all-atom energy surfaces. In the second stage (described in the accompanying article), the relative weights of the energy terms are optimized so as to locate the native structures of selected test proteins as the lowest energy structures. The goal of the work in the present study is to parameterize physically reasonable functional forms of the potentials of mean force for side-chain interactions. The potentials are of both radial and anisotropic type. Radial potentials include the Lennard-Jones and the shifted Lennard-Jones potential (with the shift parameter independent of orientation). To treat the angular dependence of side-chain interactions, three functional forms of the potential that were designed previously to describe anisotropic systems are evaluated: Berne-Pechukas (dilated Lennard-Jones); Gay-Berne (shifted Lennard-Jones with orientation-dependent shift parameters); and Gay-Berne-Vorobjev (the same as the preceding one, but with one more set of variable parameters). These functional forms were used to parameterize, within a short-distance range, the potentials of mean force for side-chain pair interactions that are related by the Boltzmann principle to the pair correlation functions determined from protein-crystal data. Parameter determination was formulated as a generalized nonlinear least-squares problem with the target function being the weighted sum of squares of the differences between calculated and “experimental” (i.e., estimated from protein-crystal data) angular, radial-angular, and radial pair correlation functions, as well as contact free energies. A set of 195 high-resolution nonhomologous structures from the Protein Data Bank was used to calculate the “experimental” values. The contact free energies were scaled by the slope of the correlation line between side-chain hydrophobicities, calculated from the contact free energies, and those determined by Fauchere and Pliška from the partition coefficients of amino acids between water and n-octanol. The methylene group served to define the reference contact free energy corresponding to that between the glycine methylene groups of backbone residues. Statistical analysis of the goodness of fit revealed that the Gay-Berne-Vorobjev anisotropic potential fits best to the experimental radial and angular correlation functions and contact free energies and therefore represents the free-energy surface of side-chain-side-chain interactions most accurately. Thus, its choice for simulations of protein structure is probably the most appropriate. However, the use of simpler functional forms is recommended, if the speed of computations is an issue. © 1997 by John Wiley & Sons, Inc. J Comput Chem 18: 849–873, 1997     

Parameterization of the Hamiltonian Dielectric Solvent (HADES) Reaction‐Field Method for the Solvation Free Energies of Amino Acid Side‐Chain Analogs

Optimization of the Hamiltonian dielectric solvent (HADES) method for biomolecular simulations in a dielectric continuum is presented with the goal of calculating accurate absolute solvation free energies while retaining the model’s accuracy in predicting conformational free‐energy differences. The solvation free energies of neutral and polar amino acid side‐chain analogs calculated by using HADES, which may optionally include nonpolar contributions, were optimized against experimental data to reach a chemical accuracy of about 0.5 kcal mol^?1. The new parameters were evaluated for charged side‐chain analogs. The HADES results were compared with explicit‐solvent, generalized Born, Poisson–Boltzmann, and QM‐based methods. The potentials of mean force (PMFs) between pairs of side‐chain analogs obtained by using HADES and explicit‐solvent simulations were used to evaluate the effects of the improved parameters optimized for solvation free energies on intermolecular potentials.     

Computing accurate potentials of mean force in electrolyte solutions with the generalized gradient-augmented harmonic Fourier beads method

We establish the accuracy of the novel generalized gradient-augmented harmonic Fourier beads (ggaHFB) method in computing free-energy profiles or potentials of mean force (PMFs) through comparison with two independent conventional techniques. In particular, we employ umbrella sampling with one dimensional weighted histogram analysis method (WHAM) and free molecular dynamics simulation of radial distribution functions to compute the PMF for the Na(+)-Cl(-) ion-pair separation to 16 A in 1.0M NaCl solution in water. The corresponding ggaHFB free-energy profile in six dimensional Cartesian space is in excellent agreement with the conventional benchmarks. We then explore changes in the PMF in response to lowering the NaCl concentration to physiological 0.3 and 0.1M, and dilute 0.0M concentrations. Finally, to expand the scope of the ggaHFB method, we formally develop the free-energy gradient approximation in arbitrary nonlinear coordinates. This formal development underscores the importance of the logarithmic Jacobian correction to reconstruct true PMFs from umbrella sampling simulations with either WHAM or ggaHFB techniques when nonlinear coordinate restraints are used with Cartesian propagators. The ability to employ nonlinear coordinates and high accuracy of the computed free-energy profiles further advocate the use of the ggaHFB method in studies of rare events in complex systems.     

Free-energy based pair-additive potentials for bulk Ni-Al systems: Application to study Ni-Al reactive alloying

We present new numerical pair-additive Al, Ni, and Al-Ni potentials by force-matching (FM) ionic force and virial data from single (bulk liquid) phase ab initio molecular dynamics (MD) simulations using the Born-Oppenheimer method. The potentials are represented by piece-wise functions (splines) and, therefore, are not constrained to a particular choice of analytical functional form. The FM method with virial constraint naturally yields a potential which maps out the ionic free-energy surface of the reference ensemble. To further improve the free energetics of the FM ensemble, the FM procedure is modified to bias the potentials to reproduce the experimental melting temperatures of the reference (FCC-Al, FCC-Ni, B2-NiAl) phases, the only macroscopic data included in the fitting set. The performance of the resultant potentials in simulating bulk metallic phases is then evaluated. The new model is applied to perform MD simulations of self-propagating exothermic reaction in Ni-Al bilayers at P = 0-5 GPa initiated at T = 1300 K. Consistent with experimental observations, the new model describes realistically a sequence of peritectic phase transformations throughout the reaction and at a realistic rate. The reaction proceeds through interlayer diffusion of Al and Ni atoms at the interface with formation of B2-NiAl in the Al melt. Such material responses have, in the past, been proven to be difficult to observe with then-existing potentials.     

Free energy perturbation study of water dimer dissociation kinetics

An efficient approach is described for using accurate ab initio calculations to determine the rates of elementary condensation and evaporation processes that lead to nucleation of aqueous aerosols. The feasibility of the method is demonstrated in an application to evaporation rates of water dimer at 230 K. The method, known as ABC-FEP (ab initio/classical free energy perturbation), begins with a calculation of the potential of mean force for the dissociation (evaporation) of small water clusters using a molecular dynamics (MD) simulation with a model potential. The free energy perturbation is used to calculate how changing from the model potential to a potential calculated from ab initio methods would alter the potential of mean force. The difference in free energy is the Boltzmann-weighted average of the difference between the ab initio and classical potential energies, with the average taken over a sample of configurations from the MD simulation. In principle, the method does not require a highly accurate model potential, though more accurate potentials require fewer configurations to achieve a small sampling error in the free energy perturbation step. To test the feasibility of obtaining accurate potentials of mean force from ab initio calculations at a modest number of configurations, the free energy perturbation method has been used to correct the errors when some standard models for bulk water (SPC, TIP4P, and TIP4PFQ) are applied to water dimer. To allow a thorough exploration of sampling issues, a highly accurate fit to results of accurate ab initio calculations, known as SAPT-5s, as been used a proxy for the ab initio calculations. It is shown that accurate values for a point on the potential of mean force can be obtained from any of the water models using ab initio calculations at only 50 configurations. Thus, this method allows accurate simulations of small clusters without the need to develop water models specifically for clusters.     

Level-Set Variational Implicit-Solvent Modeling of Biomolecules with the Coulomb-Field Approximation

Central in the variational implicit-solvent model (VISM) [Dzubiella, Swanson, and McCammon Phys. Rev. Lett.2006, 96, 087802 and J. Chem. Phys.2006, 124, 084905] of molecular solvation is a mean-field free-energy functional of all possible solute-solvent interfaces or dielectric boundaries. Such a functional can be minimized numerically by a level-set method to determine stable equilibrium conformations and solvation free energies. Applications to nonpolar systems have shown that the level-set VISM is efficient and leads to qualitatively and often quantitatively correct results. In particular, it is capable of capturing capillary evaporation in hydrophobic confinement and corresponding multiple equilibrium states as found in molecular dynamics (MD) simulations. In this work, we introduce into the VISM the Coulomb-field approximation of the electrostatic free energy. Such an approximation is a volume integral over an arbitrary shaped solvent region, requiring no solutions to any partial differential equations. With this approximation, we obtain the effective boundary force and use it as the "normal velocity" in the level-set relaxation. We test the new approach by calculating solvation free energies and potentials of mean force for small and large molecules, including the two-domain protein BphC. Our results reveal the importance of coupling polar and nonpolar interactions in the underlying molecular systems. In particular, dehydration near the domain interface of BphC subunits is found to be highly sensitive to local electrostatic potentials as seen in previous MD simulations. This is a first step toward capturing the complex protein dehydration process by an implicit-solvent approach.     

Overcoming free energy barriers using unconstrained molecular dynamics simulations

Association of unconstrained molecular dynamics (MD) and the formalisms of thermodynamic integration and average force [Darve and Pohorille, J. Chem. Phys. 115, 9169 (2001)] have been employed to determine potentials of mean force. When implemented in a general MD code, the additional computational effort, compared to other standard, unconstrained simulations, is marginal. The force acting along a chosen reaction coordinate xi is estimated from the individual forces exerted on the chemical system and accumulated as the simulation progresses. The estimated free energy derivative computed for small intervals of xi is canceled by an adaptive bias to overcome the barriers of the free energy landscape. Evolution of the system along the reaction coordinate is, thus, limited by its sole self-diffusion properties. The illustrative examples of the reversible unfolding of deca-L-alanine, the association of acetate and guanidinium ions in water, the dimerization of methane in water, and its transfer across the water liquid-vapor interface are examined to probe the efficiency of the method.     

Determination of side‐chain‐rotamer and side‐chain and backbone virtual‐bond‐stretching potentials of mean force from AM1 energy surfaces of terminally‐blocked amino‐acid residues,for coarse‐grained simulations of protein structure and folding. I. The method

In this and the accompanying article, we report the development of new physics‐based side‐chain‐rotamer and virtual‐bond‐deformation potentials which now replace the respective statistical potentials used so far in our physics‐based united‐reside UNRES force field for large‐scale simulations of protein structure and dynamics. In this article, we describe the methodology for determining the corresponding potentials of mean force (PMF's) from the energy surfaces of terminally‐blocked amino‐acid residues calculated with the AM1 quantum‐mechanical semiempirical method. The approach is based on minimization of the AM1 energy for fixed values of the angles λ for rotation of the peptide groups about the C^α ··· C^α virtual bonds, and for fixed values of the side‐chain dihedral angles χ, which formed a multidimensional grid. A harmonic‐approximation approach was developed to extrapolate from the energy at a given grid point to other points of the conformational space to compute the respective contributions to the PMF. To test the applicability of the harmonic approximation, the rotamer PMF's of alanine and valine obtained with this approach have been compared with those obtained by using a Metropolis Monte Carlo method. The PMF surfaces computed with the harmonic approximation are more rugged and have more pronounced minima than the MC‐calculated surfaces but the harmonic‐approximation‐and MC‐calculated PMF values are linearly correlated. The potentials derived with the harmonic approximation are, therefore, appropriate for UNRES for which the weights (scaling factors) of the energy terms are determined by force‐field optimization for foldability. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Molecular cluster decay viewed as escape from a potential of mean force

We show that evaporation from a quasistable molecular cluster may be treated as a kinetic problem involving the stochastically driven escape of a molecule from a potential of mean force. We derive expressions for the decay rate, and a relationship between the depth of the potential and the change in system free energy upon loss of a molecule from the cluster. This establishes a connection between kinetic and thermodynamic treatments of evaporation, but also reveals differences in the prefactor in the rate expression. We perform constant energy molecular dynamics simulations of cluster dynamics to calculate potentials of mean force, friction coefficients and effective temperatures for use in the kinetic analysis, and to compare the results with the directly observed escape rates. We also use the simulations to estimate the escape rates by a probabilistic analysis. It is much more efficient to calculate the decay rate by the methods we have developed than it is to monitor escape directly, making these approaches potentially useful for the assessment of molecular cluster stability.     

Calculation of free energy surfaces using the methods of thermodynamic perturbation theory

A molecular dynamics method for determining the free energy difference between systems separated in configuration space has been developed. With this new approach, which is based on thermodynamic perturbation techniques, potentials of mean force for conformational changes may be calculated. As a test of the method, the potential of mean force and radial distribution function for liquid argon have been computed. The results are in good agreement with those obtained from an ordinary simulation.     

Toward efficient chemical potential calculations by expanded ensemble simulations; to make the free energy pathway fairly level

A scheme is suggested of how to construct good bias potentials ("balancing factors") to be used in expanded ensemble (EE) calculations of chemical potentials of solutions. A combination of two strategies are used: (i) to use a pathway for particle insertions that avoids large variations in free energy and (ii) to use calculated free energy derivatives to construct a bias potential that makes the pathway fairly level. Only a few very short simulations are needed to accomplish the latter, and then, a full EE simulation is done to obtain the chemical potential. By practical calculations of the chemical potential of benzene, cyclohexane, and benzylamine in water, it is shown that this method is at least equally efficient to the recent adaptive EE (AEE) method by Aberg et al. (J. Chem. Phys. 2004, 120, 3370). Furthermore, the new method provides an alternative strategy that complements existing EE methods.     

A general framework for non-Boltzmann Monte Carlo sampling

Non-Boltzmann sampling (NBS) methods have been extensively employed in recent years, mainly due to their ability to enhance ergodicity in simulations of complex systems. In addition, they make possible reliable computation of equilibrium properties (ensemble averages, free-energy differences, and potentials of mean force) over continuous ranges of thermodynamic conditions. In this work, we put forward a general and systematic framework for NBS methods that allows a single set of equations and procedures to be applied to diverse systems. Moreover, we show how to exploit simulation data most effectively by obtaining continuous profiles of any mechanical properties, including structural quantities not directly related to the ensemble parameters. Finally, we demonstrate the usefulness of the developed formulation by applying it to spin systems, Lennard-Jones fluids, and a model protein molecule (both in isolation and in the proximity of a flat wall).     

Development of a quantum mechanics-based free-energy perturbation method: use in the calculation of relative solvation free energies

Free-energy perturbation (FEP) is considered the most accurate computational method for calculating relative solvation and binding free-energy differences. Despite some success in applying FEP methods to both drug design and lead optimization, FEP calculations are rarely used in the pharmaceutical industry. One factor limiting the use of FEP is its low throughput, which is attributed in part to the dependence of conventional methods on the user's ability to develop accurate molecular mechanics (MM) force field parameters for individual drug candidates and the time required to complete the process. In an attempt to find an FEP method that could eventually be automated, we developed a method that uses quantum mechanics (QM) for treating the solute, MM for treating the solute surroundings, and the FEP method for computing free-energy differences. The thread technique was used in all transformations and proved to be essential for the successful completion of the calculations. Relative solvation free energies for 10 structurally diverse molecular pairs were calculated, and the results were in close agreement with both the calculated results generated by conventional FEP methods and the experimentally derived values. While considerably more CPU demanding than conventional FEP methods, this method (QM/MM-based FEP) alleviates the need for development of molecule-specific MM force field parameters and therefore may enable future automation of FEP-based calculations. Moreover, calculation accuracy should be improved over conventional methods, especially for calculations reliant on MM parameters derived in the absence of experimental data.     

Free energy surfaces for the interaction of D-glucose with planar aromatic groups in aqueous solution

Multidimensional potentials of mean force for the interactions in aqueous solution of both anomers of D-glucopyranose with two planar aromatic molecules, indole and para-methyl-phenol, have been calculated using molecular dynamics simulations with umbrella sampling and were subsequently used to estimate binding free energies. Indole and para-methyl-phenol serve as models for the side chains of the amino acids tryptophan and tyrosine, respectively. In all cases, a weak affinity between the glucose molecules and the flat aromatic surfaces was found. The global minimum for these interactions was found to be for the case when the pseudoplanar face of β-D-glucopyranose is stacked against the planar surfaces of the aromatic residues. The calculated binding free energies are in good agreement with both experiment and previous simulations. The multidimensional free energy maps suggest a mechanism that could lend kinetic stability to the complexes formed by sugars bound to sugar-binding proteins.     

The potential of mean force on halide ions near the Cu(100) surface

This paper reports an investigation of the phenomenon of specific adsorption of halide ions on a Cu(100) surface using Monte Carlo simulations. The system was modeled by considering each ion in a water lamina placed between two copper walls. The potentials used in simulations were constructed by fitting to results of quantum calculations. The solvent contribution to the potential of mean force (pmf) was calculated for each of the four halide ions using the free energy perturbation method. Given the difficulty of finding a reliable ion–metal potential, several alternatives, representing extremal models, were tested in combination with the solvent mean force on the ions, F⁻, Cl⁻, Br⁻ or I⁻. The results for the pmf on an ion near the metal surface are discussed in the light of the experimental data available. The sensitivity of the results to the type of ion–metal potential used in the simulations is stressed.     

Modification and optimization of the united-residue (UNRES) potential energy function for canonical simulations. I. Temperature dependence of the effective energy function and tests of the optimization method with single training proteins

We report the modification and parametrization of the united-residue (UNRES) force field for energy-based protein structure prediction and protein folding simulations. We tested the approach on three training proteins separately: 1E0L (beta), 1GAB (alpha), and 1E0G (alpha + beta). Heretofore, the UNRES force field had been designed and parametrized to locate native-like structures of proteins as global minima of their effective potential energy surfaces, which largely neglected the conformational entropy because decoys composed of only lowest-energy conformations were used to optimize the force field. Recently, we developed a mesoscopic dynamics procedure for UNRES and applied it with success to simulate protein folding pathways. However, the force field turned out to be largely biased toward -helical structures in canonical simulations because the conformational entropy had been neglected in the parametrization. We applied the hierarchical optimization method, developed in our earlier work, to optimize the force field; in this method, the conformational space of a training protein is divided into levels, each corresponding to a certain degree of native-likeness. The levels are ordered according to increasing native-likeness; level 0 corresponds to structures with no native-like elements, and the highest level corresponds to the fully native-like structures. The aim of optimization is to achieve the order of the free energies of levels, decreasing as their native-likeness increases. The procedure is iterative, and decoys of the training protein(s) generated with the energy function parameters of the preceding iteration are used to optimize the force field in a current iteration. We applied the multiplexing replica-exchange molecular dynamics (MREMD) method, recently implemented in UNRES, to generate decoys; with this modification, conformational entropy is taken into account. Moreover, we optimized the free-energy gaps between levels at temperatures corresponding to a predominance of folded or unfolded structures, as well as to structures at the putative folding-transition temperature, changing the sign of the gaps at the transition temperature. This enabled us to obtain force fields characterized by a single peak in the heat capacity at the transition temperature. Furthermore, we introduced temperature dependence to the UNRES force field; this is consistent with the fact that it is a free-energy and not a potential energy function. beta     

Sequence-dependent conformational energy of DNA derived from molecular dynamics simulations: toward understanding the indirect readout mechanism in protein-DNA recognition

Sequence dependence of DNA conformation plays a crucial role in its recognition by proteins and ligands. To clarify the relationship between sequence and conformation, it is necessary to quantify the conformational energy and specificity of DNA. Here, we make a systematic analysis of dodecamer DNA structures including all the 136 unique tetranucleotide sequences at the center by molecular dynamics simulations. Using a simplified conformational model with six parameters to describe the geometry of adjacent base pairs and harmonic potentials along these coordinates, we estimated the equilibrium conformational parameters and the harmonic potentials of mean force for the central base-pair steps from many trajectories of the simulations. This enabled us to estimate the conformational energy and the specificity for any given DNA sequence and structure. We tested our method by using sequence-structure threading to estimate the conformational energy and the Z-score as a measure of specificity for many B-DNA and A-DNA crystal structures. The average Z-scores were negative for both kinds of structures, indicating that the potential of mean force from the simulation is capable of predicting sequence specificity for the crystal structures and that it may be used to study the sequence specificity of both types of DNA. We also estimated the positional distribution of conformational energy and Z-score within DNA and showed that they are strongly position dependent. This analysis enabled us to identify particular conformations responsible for the specificity. The presented results will provide an insight into the mechanisms of DNA sequence recognition by proteins and ligands.     

Particle-based multiscale coarse graining with density-dependent potentials: application to molecular crystals (hexahydro-1,3,5-trinitro-s-triazine)

We describe the development of isotropic particle-based coarse-grain models for crystalline hexahydro-1,3,5-trinitro-s-triazine (RDX). The coarse graining employs the recently proposed multiscale coarse-graining (MS-CG) method, which is a particle-based force-matching approach for deriving free-energy effective interaction potentials. Though one-site and four-site coarse-grain (CG) models were parameterized from atomistic simulations of non-ordered (molten and ambient temperature amorphous) systems, the focus of the paper is a detailed study of the one-site model with a brief recourse to the four-site model. To improve the ability of the one-site model to be applied to crystalline phases at various pressures, it was found necessary to include explicit dependence on a particle density, and a new theory of local density-dependent MS-CG potentials is subsequently presented. The density-dependency is implemented through interpolation of MS-CG force fields derived at a preselected set of reference densities. The computationally economical procedure for obtaining the reference force fields starting from the interaction at ambient density is also described. The one-site MS-CG model adequately describes the atomistic lattice structure of α-RDX at ambient and high pressures, elastic and vibrational properties, pressure-volume curve up to P = 10 GPa, and the melting temperature. In the molten state, the model reproduces the correct pair structure at different pressures as well as higher order correlations. The potential of the MS-CG model is further evaluated in simulations of shocked crystalline RDX.