Dual-topology/dual-coordinate free-energy simulation using QM/MM force field

We have developed a dual-topology/dual-coordinate free-energy simulation method for use with a QM/MM force field. By combining two parallel processes into one alchemical process, we are able to compute the double free-energy difference (delta deltaF) within a single simulation, which eliminates half of the expensive quantum-mechanical simulation in general. The method has been tested in computing the solvation free-energy differences of several molecular pairs and shows close agreement with experimental results.     

A reversible minimum-to-minimum mapping method for the calculation of free-energy differences

A general method is introduced for the calculation of the free-energy difference between two systems, 0 and 1, with configuration spaces omega(0), omega(1) of the same dimensionality. The method relies upon establishing a objective mapping between disjoint subsets gamma(i)(0) of omega(0) and corresponding disjoint subsets gamma(i)(1) of omega(1), and averaging a function of the ratio of configurational integrals over gamma(i)(0) and gamma(i)(1) with respect to the probability densities of the two systems. The mapped subsets gamma(i)(0) and gamma(i)(1) need not span the entire configuration spaces omega(0) and omega(1). The method is applied for the calculation of the excess chemical potential mu(ex) in a Lennard-Jones (LJ) fluid. In this case, omega(0) is the configuration space of a (N-1) real molecule plus one ideal-gas molecule system, while omega(1) is the configuration space of a N real molecule system occupying the same volume. Gamma(i)(0) and gamma(i)(1) are constructed from hyperspheres of the same radius centered at minimum-energy configurations of a set of "active" molecules lying within distance a from the ideal-gas molecule and the last real molecule, respectively. An algorithm is described for sampling gamma(i)(0) and gamma(i)(1) given a point in omega(0) or in omega(1). The algorithm encompasses three steps: "quenching" (minimization with respect to the active-molecule degrees of freedom), "mutation" (gradual conversion of the ideal-gas molecule into a real molecule, with simultaneous minimization of the energy with respect to the active-molecule degrees of freedom), and "excitation" (generation of points on a hypersphere centered at the active-molecule energy minimum). These steps are also carried out in reverse, as required by the bijective nature of the mapping. The mutation step, which establishes a reversible mapping between energy minima with respect to the active degrees of freedom of systems 0 and 1, ensures that excluded volume interactions emerging in the process of converting the ideal-gas molecule into a real molecule are relieved through appropriate rearrangement of the surrounding active molecules. Thus, the insertion problem plaguing traditional methods for the calculation of chemical potential at high densities is alleviated. Results are presented at two state points of the LJ system for a variety of radii a of the active domain. It is shown that the estimated values of mu(ex) are correct in all cases and subject to an order of magnitude lower statistical uncertainty than values based on the same number of Widom [J. Chem. Phys. 39, 2808 (1963)] insertions at high fluid densities. Optimal settings for the new algorithm are identified and distributions of the quantities involved in it [number of active molecules, energy at the sampled minima of systems 0 and 1, and free-energy differences between subsets gamma(i)(0) and gamma(i)(1) that are mapped onto each other] are explored.     

Robust and accurate method for free-energy calculation of charged molecular systems

A new approach is presented to eliminate the problem of creation and/or annihilation of atoms in free-energy calculations of charged molecular systems. The method employs a damping potential in the Ewald summation scheme, which is an exact solution of the electrostatics for three-dimensional periodic systems. The proposed method enables entire molecules to be mutated from a noninteracting (ideal) state in an efficient and robust way, thus providing a means by which accurate absolute free energies of structurally complex molecules can be determined. This methodology will enable chemical and phase equilibria to be determined for large molecular species with significant charge distributions, e.g., biomolecules and drugs.     

Mapping the backbone dihedral free-energy surfaces in small peptides in solution using adiabatic free-energy dynamics

Ramachandran surfaces for the alanine di- and tripeptides in gas phase and solution are mapped out using the recently introduced adiabatic free-energy dynamics (AFED) approach introduced by Rosso et al. (J. Chem. Phys. 2002, 116, 4389) as applied to the CHARMM22 force field. It is shown that complete surfaces can be mapped out with an order of magnitude of greater efficiency with the AFED approach than they can using the popular umbrella sampling method. In the alanine dipeptide, it is found, in agreement with numerous other studies using the CHARMM22 force field, that the lowest free-energy structure is the extended beta conformation, (phi, psi) = (-81, 81), while in solution, the extended beta, (phi, psi) = (-81, 153) and right-handed alpha-helical, (phi, psi) = (-81, 63) conformations are nearly isoenergetic. In solution, a secondary minimum at (phi, psi) = (63, -81), corresponding to a C(7)ax conformation, occurs approximately 2.3 kcal/mol above the global free-energy minimum. The alanine tripeptide, a system that has received considerably less attention in the literature, is found to exhibit a similar structure to the alanine dipeptide with the extended beta conformation being the free-energy minimum in the gas phase and the beta and right-handed alpha-helical conformations being isoenergetic in solution. These studies indicate that the AFED method can be a powerful tool for studying multidimensional free-energy surfaces in complex systems.     

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.     

Scalable free energy calculation of proteins via multiscale essential sampling

A multiscale simulation method, "multiscale essential sampling (MSES)," is proposed for calculating free energy surface of proteins in a sizable dimensional space with good scalability. In MSES, the configurational sampling of a full-dimensional model is enhanced by coupling with the accelerated dynamics of the essential degrees of freedom. Applying the Hamiltonian exchange method to MSES can remove the biasing potential from the coupling term, deriving the free energy surface of the essential degrees of freedom. The form of the coupling term ensures good scalability in the Hamiltonian exchange. As a test application, the free energy surface of the folding process of a miniprotein, chignolin, was calculated in the continuum solvent model. Results agreed with the free energy surface derived from the multicanonical simulation. Significantly improved scalability with the MSES method was clearly shown in the free energy calculation of chignolin in explicit solvent, which was achieved without increasing the number of replicas in the Hamiltonian exchange.     

Peptide folding using multiscale coarse-grained models

The multiscale coarse-graining (MS-CG) method has been previously used to describe the equilibrium properties of peptides. The present study reveals that MS-CG models of alpha-helical polyalanine and the beta-hairpin V 5PGV 5 possess the capacity to efficiently refold in simulations initiated from unfolded configurations. The MS-CG peptides exhibit free energy landscapes that are funneled toward folded configurations and two-state folding behavior, consistent with the known characteristics of small, rapidly folding peptides. Moreover, the models demonstrate enhanced sampling capabilities when compared to systems with full atomic detail. The significance of these observations with respect to the theoretical basis of the MS-CG approach is discussed. The MS-CG peptides were used to reconstruct atomically detailed configurations in order to evaluate the extent to which MS-CG ensembles embody all-atom peptide free energy landscapes. Ensembles obtained from these reconstructed configurations display good agreement with the all-atom simulation data used to generate the MS-CG models and also corroborate the presence of features observed in the MS-CG peptide free energy landscapes. These findings suggest that MS-CG models may be of significant utility in the study of peptide folding.     

An algorithm to find minimum free-energy paths using umbrella integration

The calculation of free-energy barriers by umbrella sampling and many other methods is hampered by the necessity for an a priori choice of the reaction coordinate along which to sample. We avoid this problem by providing a method to search for saddle points on the free-energy surface in many coordinates. The necessary gradients and Hessians of the free energy are obtained by multidimensional umbrella integration. We construct the minimum free-energy path by following the gradient down to minima on the free-energy surface. The change of free energy along the path is obtained by integrating out all coordinates orthogonal to the path. While we expect the method to be applicable to large systems, we test it on the alanine dipeptide in vacuum. The minima, transition states, and free-energy barriers agree well with those obtained previously with other methods.     

Absolute free-energy calculations of liquids using a harmonic reference state

Absolute free-energy methods provide a potential solution to the overlap problem in free-energy calculations. In this paper, we report an extension of the previously published confinement method (J. Phys. Chem. B 2006, 110, 17212-20) to fluid simulations. Absolute free energies of liquid argon and liquid water are obtained accurately and compared with results from thermodynamic integration. The method works by transforming the liquid state into a harmonic, solid reference state. This is achieved using a special restraint potential that allows molecules to change their restraint position during the simulation, which circumvents the need for the molecules to sample the full extent of their translational freedom. The absolute free energy of the completely restrained reference state is obtained from a normal mode calculation. Because of the generic reference state used, the method is applicable to nonhomogeneous, diffusive systems and could provide an alternative method in situations in which solute annihilation fails due to the size of the solute. Potential applications include calculation of solvation energies of large molecules and free energies of peptide conformational changes in explicit solvent.     

Order parameters for macromolecules: application to multiscale simulation

Order parameters (OPs) characterizing the nanoscale features of macromolecules are presented. They are generated in a general fashion so that they do not need to be redesigned with each new application. They evolve on time scales much longer than 10(-14) s typical for individual atomic collisions/vibrations. The list of OPs can be automatically increased, and completeness can be determined via a correlation analysis. They serve as the basis of a multiscale analysis that starts with the N-atom Liouville equation and yields rigorous Smoluchowski/Langevin equations of stochastic OP dynamics. Such OPs and the multiscale analysis imply computational algorithms that we demonstrate in an application to ribonucleic acid structural dynamics for 50 ns.     

A multiscale simulation system for the prediction of drug-induced cardiotoxicity

The preclinical assessment of drug-induced ventricular arrhythmia, a major concern for regulators, is typically based on experimental or computational models focused on the potassium channel hERG (human ether-a-go-go-related gene, K(v)11.1). Even if the role of this ion channel in the ventricular repolarization is of critical importance, the complexity of the events involved make the cardiac safety assessment based only on hERG has a high risk of producing either false positive or negative results. We introduce a multiscale simulation system aiming to produce a better cardiotoxicity assessment. At the molecular scale, the proposed system uses a combination of docking simulations on two potassium channels, hERG and KCNQ1, plus three-dimensional quantitative structure-activity relationship modeling for predicting how the tested compound will block the potassium currents IK(r) and IK(s). The obtained results have been introduced in electrophysiological models of the cardiomyocytes and the ventricular tissue, allowing the direct prediction of the drug effects on electrocardiogram simulations. The usefulness of the whole method is illustrated by predicting the cardiotoxic effect of several compounds, including some examples in which classic hERG-based models produce false positive or negative results, yielding correct predictions for all of them. These results can be considered a proof of concept, suggesting that multiscale prediction systems can be suitable for being used for preliminary screening in lead discovery, before the compound is physically available, or in early preclinical development when they can be fed with experimentally obtained data.     

A multiscale simulation model for poly(ethylene oxide)

An overview of the recent development of a multiscale simulation of amorphous polymeric materials at the bulk density is presented. Poly(ethylene oxide), (PEO), (CH₃O-[CH₂-CH₂-O]nCH₃) was selected to illustrate the method. The model starts from an ab initio quantum chemistry to obtain the statistical weights of polymer conformation based on the rotational isomeric state (RIS) theory. PEO chains were then mapped to a coarse-grained model using the modified RIS model onto the second nearest neighbor diamond (2nnd) lattice. The average non-bonded interactions were treated by the discretized Lennard-Jones (LJ) potential. Bulk PEO melts with molecular weight up to 8000 g/mol was generated and equilibrated. The on-lattice properties such as molecular size and conformational statistics agree well with the theory. Fully atomistic amorphous PEO models can be obtained by the reverse-mapping procedure to recover the missing atoms. After an energy minimization step, properties including torsional angle distribution, solubility parameter and static neutron scattering structure factor are in good agreement with experimental results.     

Theoretical study of the neutral hydrolysis of methyl formate via a concerted and stepwise water-assisted mechanism using free-energy curves and molecular dynamics simulation

A procedure previously described by us is used for the theoretical study of chemical reactions in solution by means of molecular dynamics simulation, with solute–solvent interaction potentials LJ (12-6-1) derived from ab initio quantum calculations. We apply the procedure to the case of the neutral hydrolysis of methyl formate, HCOOCH₃ + 3H₂O → HCOOH + CH₃OH + 2H₂O in aqueous solution, via concerted and stepwise water-assisted mechanisms. We use the solvent as reaction coordinate, and the free-energy curves for the calculation of the activation energies. The theoretical calculation for the thermodynamics of this hydrolysis reaction in aqueous solution, assisted by three water molecules, is in agreement with the available experimental information. In particular our study gives values of ΔG ^≠ = 28.88 and 28.17 kcal/mol for the concerted and stepwise mechanisms, close to the experimental activation barrier of 28.8 kcal/mol, and a significant improvement over the values of 48.05 and 45.66 kcal/mol found in another similar study using the PCM model.     

Boundary problems in the calculation of light scattering intensities from liquids using computer simulation

Computer simulation of molecular dynamics in the Lennard-Jones model of liquid argon at the triple point density has been used to calculate depolarised scattering intensities based on a pairwise-additive dipole-induced-dipole mechanism. It is shown that large systematic errors arise from boundary effects as a result of the finite size of the model liquid samples. This severely limits the use of such calculations in evaluating light-scattering models.     

A constrained approach to multiscale stochastic simulation of chemically reacting systems

Stochastic simulation of coupled chemical reactions is often computationally intensive, especially if a chemical system contains reactions occurring on different time scales. In this paper, we introduce a multiscale methodology suitable to address this problem, assuming that the evolution of the slow species in the system is well approximated by a Langevin process. It is based on the conditional stochastic simulation algorithm (CSSA) which samples from the conditional distribution of the suitably defined fast variables, given values for the slow variables. In the constrained multiscale algorithm (CMA) a single realization of the CSSA is then used for each value of the slow variable to approximate the effective drift and diffusion terms, in a similar manner to the constrained mean-force computations in other applications such as molecular dynamics. We then show how using the ensuing Fokker-Planck equation approximation, we can in turn approximate average switching times in stochastic chemical systems.     

The Ising free-energy functional

We show that the Ising free-energy functional f(T, X) yields a second-order phase transition. The value of X, a generalized order parameter, which minimizes the functional is the real-order parameter, Ω = tanh[Ω(T_C/T)] where T_c = J/(2k) and J is the Ising coupling constant. The Ising theory is applicable to the second-order phase transitions in ferromagnets, ferroelectrics, and superconductors. The constant temperature slices of f(T, X) yields the Ising functional which is valid over the entire range of T and which is identical to the Landau functional for T/T_C close to unity. © 1996 John Wiley & Sons, Inc.     

Exploring the free-energy landscape of a short peptide using an average force

The reversible folding of deca-alanine is chosen as a test case for characterizing a method that uses an adaptive biasing force (ABF) to escape from the minima and overcome the barriers of the free-energy landscape. This approach relies on the continuous estimation of a biasing force that yields a Hamiltonian in which no average force is exerted along the ordering parameter xi. Optimizing the parameters that control how the ABF is applied, the method is shown to be extremely effective when a nonequivocal ordering parameter can be defined to explore the folding pathway of the peptide. Starting from a beta-turn motif and restraining xi to a region of the conformational space that extends from the alpha-helical state to an ensemble of extended structures, the ABF scheme is successful in folding the peptide chain into a compact alpha helix. Sampling of this conformation is, however, marginal when the range of xi values embraces arrangements of greater compactness, hence demonstrating the inherent limitations of free-energy methods when ambiguous ordering parameters are utilized.     

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.