Virtual‐system‐coupled adaptive umbrella sampling to compute free‐energy landscape for flexible molecular docking

A novel enhanced conformational sampling method, virtual‐system‐coupled adaptive umbrella sampling (V‐AUS), was proposed to compute 300‐K free‐energy landscape for flexible molecular docking, where a virtual degrees of freedom was introduced to control the sampling. This degree of freedom interacts with the biomolecular system. V‐AUS was applied to complex formation of two disordered amyloid‐β (Aβ_30–35) peptides in a periodic box filled by an explicit solvent. An interpeptide distance was defined as the reaction coordinate, along which sampling was enhanced. A uniform conformational distribution was obtained covering a wide interpeptide distance ranging from the bound to unbound states. The 300‐K free‐energy landscape was characterized by thermodynamically stable basins of antiparallel and parallel β‐sheet complexes and some other complex forms. Helices were frequently observed, when the two peptides contacted loosely or fluctuated freely without interpeptide contacts. We observed that V‐AUS converged to uniform distribution more effectively than conventional AUS sampling did. © 2015 Wiley Periodicals, Inc.     

Two‐dimensional replica‐exchange method for predicting protein–ligand binding structures

We have developed a two‐dimensional replica‐exchange method for the prediction of protein–ligand binding structures. The first dimension is the umbrella sampling along the reaction coordinate, which is the distance between a protein binding pocket and a ligand. The second dimension is the solute tempering, in which the interaction between a ligand and a protein and water is weakened. The second dimension is introduced to make a ligand follow the umbrella potential more easily and enhance the binding events, which should improve the sampling efficiency. As test cases, we applied our method to two protein‐ligand complex systems (MDM2 and HSP 90‐alpha). Starting from the configuration in which the protein and the ligand are far away from each other in each system, our method predicted the ligand binding structures in excellent agreement with the experimental data from Protein Data Bank much faster with the improved sampling efficiency than the replica‐exchange umbrella sampling method that we have previously developed. © 2013 Wiley Periodicals, Inc.     

Free‐energy differences between states with different conformational ensembles

Multiple conformations separated by high‐energy barriers represent a challenging problem in free‐energy calculations due to the difficulties in achieving adequate sampling. We present an application of thermodynamic integration (TI) in conjunction with the local elevation umbrella sampling (LE/US) method to improve convergence in alchemical free‐energy calculations. TI‐LE/US was applied to the guanosine triphosphate (GTP) to 8‐Br‐GTP perturbation, molecules that present high‐energy barriers between the anti and syn states and that have inverted preferences for those states. The convergence and reliability of TI‐LE/US was assessed by comparing with previous results using the enhanced‐sampling one‐step perturbation (OSP) method. A linear interpolation of the end‐state biasing potentials was sufficient to dramatically improve sampling along the chosen reaction coordinate. Conformational free‐energy differences were also computed for the syn and anti states and compared to experimental and theoretical results. Additionally, a coupled OSP with LE/US was carried out, allowing the calculation of conformational and alchemical free energies of GTP and 8‐substituted GTP analogs. © 2013 Wiley Periodicals, Inc.     

Multidimensional virtual-system coupled canonical molecular dynamics to compute free-energy landscapes of peptide multimer assembly

An enhanced-sampling method termed multidimensional virtual-system coupled canonical molecular dynamics (mD-VcMD) method is developed. In many cases, generalized-ensemble methods realizing enhanced sampling, for example, adaptive umbrella sampling, apply an effective potential, which is derived from temporarily assumed canonical distribution as a function of one or more arbitrarily defined reaction coordinates. However, it is not straightforward to estimate the appropriate canonical distribution, especially for cases applying multiple reaction coordinates. The current method, mD-VcMD, does not rely on the form of the canonical distribution. Therefore, it is practically useful to explore a high-dimensional reaction-coordinate space. In this article, formulation of mD-VcMD and its evaluation with the simple molecular models consisting of three or four alanine peptides are presented. We confirmed that mD-VcMD efficiently searched 2D and 3D reaction-coordinate spaces defined as interpeptide distances. Direct comparisons with results of long-term canonical MD simulations revealed that mD-VcMD produces correct canonical ensembles. © 2019 Wiley Periodicals, Inc.     

Hybrid quantum/classical path integral approach for simulation of hydrogen transfer reactions in enzymes

A hybrid quantum/classical path integral Monte Carlo (QC-PIMC) method for calculating the quantum free energy barrier for hydrogen transfer reactions in condensed phases is presented. In this approach, the classical potential of mean force along a collective reaction coordinate is calculated using umbrella sampling techniques in conjunction with molecular dynamics trajectories propagated according to a mapping potential. The quantum contribution is determined for each configuration along the classical trajectory with path integral Monte Carlo calculations in which the beads move according to an effective mapping potential. This type of path integral calculation does not utilize the centroid constraint and can lead to more efficient sampling of the relevant region of conformational space than free-particle path integral sampling. The QC-PIMC method is computationally practical for large systems because the path integral sampling for the quantum nuclei is performed separately from the classical molecular dynamics sampling of the entire system. The utility of the QC-PIMC method is illustrated by an application to hydride transfer in the enzyme dihydrofolate reductase. A comparison of this method to the quantized classical path and grid-based methods for this system is presented.     

Enhanced sampling simulations of DNA step parameters

A novel approach for the selection of step parameters as reaction coordinates in enhanced sampling simulations of DNA is presented. The method uses three atoms per base and does not require coordinate overlays or idealized base pairs. This allowed for a highly efficient implementation of the calculation of all step parameters and their Cartesian derivatives in molecular dynamics simulations. Good correlation between the calculated and actual twist, roll, tilt, shift, and slide parameters is obtained, while the correlation with rise is modest. The method is illustrated by its application to the methylated and unmethylated 5′‐CATGTGACGTCACATG‐3′ double stranded DNA sequence. One‐dimensional umbrella simulations indicate that the flexibility of the central CG step is only marginally affected by methylation. © 2014 Wiley Periodicals, Inc.     

Rapid sampling of all‐atom peptides using a library‐based polymer‐growth approach

We adapted existing polymer growth strategies for equilibrium sampling of peptides described by modern atomistic forcefields with a simple uniform dielectric solvent. The main novel feature of our approach is the use of precalculated statistical libraries of molecular fragments. A molecule is sampled by combining fragment configurations—of single residues in this study—which are stored in the libraries. Ensembles generated from the independent libraries are reweighted to conform with the Boltzmann‐factor distribution of the forcefield describing the full molecule. In this way, high‐quality equilibrium sampling of small peptides (4–8 residues) typically requires less than one hour of single‐processor wallclock time and can be significantly faster than Langevin simulations. Furthermore, approximate, clash‐free ensembles can be generated for larger peptides (up to 32 residues in this study) in less than a minute of single‐processor computing. We discuss possible applications of our growth procedure to free energy calculation, fragment assembly protein‐structure prediction protocols, and to “multi‐resolution” sampling. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011     

Repeated‐annealing sampling combined with multicanonical algorithm for conformational sampling of bio‐molecules

A novel conformational sampling method (repeated‐annealing sampling method) is proposed to execute an efficient conformational sampling at a reasonable computational cost. In the method, a molecular dynamics simulation is done with repeating an elemental process. An elemental process consists of four subprocesses: high‐temperature run, annealing, room‐temperature run, and fast heating. The sampling is done automatically according to a temperature‐control schedule. The room‐temperature run is treated with the multicanonical algorithm, and the other subprocesses are done with the conventional molecular dynamics algorithm. The method, differing from the generalized ensemble methods recently developed, is not warrantable to give the canonical ensemble because of the nonphysical process in the annealing. However, we observed that the slower the annealing and the longer the high‐temperature run, the closer the sampled conformations to those of the canonical ensemble. A test was performed with tri‐N‐acetyl‐D ‐glucosamine in vacuo, and the results were compared with those from the conventional multicanonical simulation. Not only the reweighted canonical distribution function but also the energy landscape were in good agreement with those from the conventional multicanonical simulation. The potential of mean force also showed a fairly good agreement with that from the conventional multicanonical simulation in the room‐temperature region. © 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1098–1106, 2001     

Hybrid quantum/classical molecular dynamics for a proton transfer reaction coupled to a dissipative bath

A hybrid quantum/classical molecular dynamics approach is applied to a proton transfer reaction represented by a symmetric double well system coupled to a dissipative bath. In this approach, the proton is treated quantum mechanically and all bath modes are treated classically. The transition state theory rate constant is obtained from the potential of mean force, which is generated along a collective reaction coordinate with umbrella sampling techniques. The transmission coefficient, which accounts for dynamical recrossings of the dividing surface, is calculated with a reactive flux approach combined with the molecular dynamics with quantum transitions surface hopping method. The hybrid quantum/classical results agree well with numerically exact results in the spatial-diffusion-controlled regime, which is most relevant for proton transfer in proteins. This hybrid quantum/classical approach has already been shown to be computationally practical for studying proton transfer in large biological systems. These results have important implications for future applications to hydrogen transfer reactions in solution and proteins.     

Effects of system net charge and electrostatic truncation on all‐atom constant pH molecular dynamics

Constant pH molecular dynamics offers a means to rigorously study the effects of solution pH on dynamical processes. Here, we address two critical questions arising from the most recent developments of the all‐atom continuous constant pH molecular dynamics (CpHMD) method: (1) What is the effect of spatial electrostatic truncation on the sampling of protonation states? (2) Is the enforcement of electrical neutrality necessary for constant pH simulations? We first examined how the generalized reaction field and force‐shifting schemes modify the electrostatic forces on the titration coordinates. Free energy simulations of model compounds were then carried out to delineate the errors in the deprotonation free energy and salt‐bridge stability due to electrostatic truncation and system net charge. Finally, CpHMD titration of a mini‐protein HP36 was used to understand the manifestation of the two types of errors in the calculated pK_a values. The major finding is that enforcing charge neutrality under all pH conditions and at all time via cotitrating ions significantly improves the accuracy of protonation‐state sampling. We suggest that such finding is also relevant for simulations with particle mesh Ewald, considering the known artifacts due to charge‐compensating background plasma. © 2014 Wiley Periodicals, Inc.     

Efficient estimation of binding free energies between peptides and an MHC class II molecule using coarse‐grained molecular dynamics simulations with a weighted histogram analysis method

We estimate the binding free energy between peptides and an MHC class II molecule using molecular dynamics (MD) simulations with the weighted histogram analysis method (WHAM). We show that, owing to its more thorough sampling in the available computational time, the binding free energy obtained by pulling the whole peptide using a coarse‐grained (CG) force field (MARTINI) is less prone to significant error induced by inadequate‐sampling than using an atomistic force field (AMBER). We further demonstrate that using CG MD to pull 3–4 residue peptide segments while leaving the remaining peptide segments in the binding groove and adding up the binding free energies of all peptide segments gives robust binding free energy estimations, which are in good agreement with the experimentally measured binding affinities for the peptide sequences studied. Our approach thus provides a promising and computationally efficient way to rapidly and reliably estimate the binding free energy between an arbitrary peptide and an MHC class II molecule. © 2017 Wiley Periodicals, Inc.     

Prospects of transition interface sampling simulations for the theoretical study of zeolite synthesis

The transition interface sampling (TIS) technique allows large free energy barriers to be overcome within reasonable simulation time, which is impossible for straightforward molecular dynamics. Still, the method does not impose an artificial driving force, but it surmounts the timescale problem by an importance sampling of true dynamical pathways. Recently, it was shown that the efficiency of TIS when calculating reaction rates is less sensitive to the choice of reaction coordinate than those of the standard free energy based techniques. This could be an important advantage in complex systems for which a good reaction coordinate is usually very difficult to find. We explain the principles of this method and discuss some of the promising applications related to zeolite formation.     

Protein folding pathways extracted by OFLOOD: Outlier FLOODing method

The Outlier FLOODing method (OFLOOD) is proposed as an efficient conformational sampling method to extract biologically rare events such as protein folding. In OFLOOD, sparse distributions (outliers in the conformational space) were regarded as relevant states for the transitions. Then, the transitions were enhanced through conformational resampling from the outliers. This evidence indicates that the conformational resampling of the sparse distributions might increase chances for promoting the transitions from the outliers to other meta‐stable states, which resembles a conformational flooding from the outliers to the neighboring clusters. OFLOOD consists of (i) detections of outliers from conformational distributions and (ii) conformational resampling from the outliers by molecular dynamics (MD) simulations. Cycles of (i) and (ii) are simply repeated. As demonstrations, OFLOOD was applied to folding of Chignolin and HP35. In both cases, OFLOOD automatically extracted folding pathways from unfolded structures with ns‐order computational costs, although µs‐order canonical MD failed to extract them. © 2014 Wiley Periodicals, Inc.     

Free energy calculation using molecular dynamics simulation combined with the three-dimensional reference interaction site model theory. II. Thermodynamic integration along a spatial reaction coordinate

We propose the thermodynamic integration along a spatial reaction coordinate using the molecular dynamics simulation combined with the three-dimensional reference interaction site model theory. This method provides a free energy calculation in solution along the reaction coordinate defined by the Cartesian coordinates of the solute atoms. The proposed method is based on the blue moon algorithm which can, in principle, handle any reaction coordinate as far as it is defined by the solute atom positions. In this article, we apply the present method to the complex formation process of the crown ether 18-Crown-6 (18C6) with the potassium ion in an aqueous solution. The separation between the geometric centers of these two molecules is taken to be the reaction coordinate for this system. The potential of mean force (PMF) becomes the maximum at the separation between the molecular centers being ～4 A?, which can be identified as the free energy barrier in the process of the molecular recognition. In a separation further than the free energy barrier, the PMF is slightly reduced to exhibit a plateau. In the region closer than the free energy barrier, approach of the potassium ion to the center of 18C6 also decreases the PMF. When the potassium ion is accommodated at the center of 18C6, the free energy is lower by -5.7 ± 0.7 kcal/mol than that at the above mentioned plateau or converged state. By comparing the results with those from the free energy calculation along the coupling parameters obtained in our previous paper [T. Miyata, Y. Ikuta, and F. Hirata, J. Chem. Phys. 133, 044114 (2010)], it is found that the effective interaction in water between 18C6 and the potassium ion vanishes beyond the molecular-center-separation of 10 A?. Furthermore, the conformation of 18C6 is found to be significantly changed depending upon the 18C6-K(+) distance. A proper conformational sampling and an accurate solvent treatment are crucial for realizing the accurate PMF, and we believe that the proposed method is useful to evaluate the PMF in a solution. A discussion upon the PMF in terms of the three-dimensional distribution function for the solvent is also presented.     

Molecular dynamics of apo-adenylate kinase: a distance replica exchange method for the free energy of conformational fluctuations

A large domain motion in adenylate kinase from E. coli (AKE) is studied with molecular dynamics. AKE undergoes a large-scale rearrangement of its lid and AMP-binding domains when the open form closes over its substrates, AMP, and Mg2+-ATP, whereby the AMP-binding and lid domains come closer to the core. The third domain, the core, is relatively stable during this motion. A reaction coordinate that monitors the distance between the AMP-binding and core domains is selected to be able to compare with the results of energy transfer experiments. Sampling along this reaction coordinate is carried out by using a distance replica exchange method (DREM), where systems that differ by a restraint potential enforcing different reaction coordinate values are independently simulated with periodic attempts at exchange of these systems. Several methods are used to study the efficiency and convergence properties of the DREM simulation and compared with an analogous non-DREM simulation. The DREM greatly accelerates the rate and extent of configurational sampling and leads to equilibrium sampling as measured by monitoring collective modes obtained from a principal coordinate analysis. The potential of mean force along the reaction coordinate reveals a rather flat region for distances from the open to a relatively closed AKE conformation. The potential of mean force for smaller distances has a distinct minimum that is quite close to that found in the closed form X-ray structure. In concert with a decrease in the reaction coordinate distance (AMP-binding-to-core distance) the lid-to-core distance of AKE also decreases. Therefore, apo AKE can fluctuate from its open form to conformations that are quite similar to its closed form X-ray structure, even in the absence of its substrates.     

Ab initio prediction of protein-ligand binding structures by replica-exchange umbrella sampling simulations

We have developed a prediction method for the binding structures of ligands with proteins. Our method consists of three steps. First, replica-exchange umbrella sampling simulations are performed along the distance between a putative binding site of a protein and a ligand as the reaction coordinate. Second, we obtain the potential of mean force (PMF) of the unbiased system using the weighted histogram analysis method and determine the distance that corresponds to the global minimum of PMF. Third, structures that have this global-minimum distance and energy values around the average potential energy are collected and analyzed using the principal component analysis. We predict the binding structure as the global-minimum free energy state on the free energy landscapes along the two major principal component axes. As test cases, we applied our method to five protein-ligand complex systems. Starting from the configuration in which the protein and the ligand are far away from each other in each system, our method predicted the ligand binding structures in excellent agreement with the experimental data from Protein Data Bank.     

Quasi-Hamiltonian equations of motion for internal coordinate molecular dynamics of polymers

Conventional molecular dynamics simulations of macromolecules require long computational times because the most interesting motions are very slow compared to the fast oscillations of bond lengths and bond angles that limit the integration time step. Simulation of dynamics in the space of internal coordinates, that is, with bond lengths, bond angles, and torsions as independent variables, gives a theoretical possibility of eliminating all uninteresting fast degrees of freedom from the system. This article presents a new method for internal coordinate molecular dynamics simulations of macromolecules. Equations of motion are derived that are applicable to branched chain molecules with any number of internal degrees of freedom. Equations use the canonical variables and they are much simpler than existing analogs. In the numerical tests the internal coordinate dynamics are compared with the traditional Cartesian coordinate molecular dynamics in simulations of a 56 residue globular protein. For the first time it was possible to compare the two alternative methods on identical molecular models in conventional quality tests. It is shown that the traditional and internal coordinate dynamics require the same time step size for the same accuracy and that in the standard geometry approximation of amino acids, that is, with fixed bond lengths, bond angles, and rigid aromatic groups, the characteristic step size is 4 fs, which is 2 times higher than with fixed bond lengths only. The step size can be increased up to 11 fs when rotation of hydrogen atoms is suppressed. © 1997 by John Wiley & Sons, Inc. J Comput Chem 18 : 1354–1364, 1997     

Beta‐hairpin restraint potentials for calculations of potentials of mean force as a function of beta‐hairpin tilt,rotation, and distance

We have developed a set of restraint potentials for β‐hairpin tilt relative to the membrane normal, β‐hairpin rotation around the β‐hairpin axis, and hairpin–hairpin distance. Such restraint potentials enable us to characterize the molecular basis of specific β‐hairpin tilt and rotation in membranes and hairpin–hairpin interactions at the atomic level by sampling their conformational space along these degrees of freedom, i.e., reaction coordinates, during molecular dynamics simulations. We illustrate the efficacy of the β‐hairpin restraint potentials by calculating the potentials of mean force (PMFs) as a function of tilt and rotation angles of protegrin‐1 (PG‐1), a β‐hairpin antimicrobial peptide, in an implicit membrane model. The peptide association in the membrane is also examined by calculating the PMFs as a function of distance between two PG‐1 peptides in various dimer interfaces. These novel restraint potentials are found to perform well in each of these cases and are expected to be a useful means to study the microscopic driving forces of insertion, tilting, and rotation of β‐hairpin peptides in membranes as well as their association in aqueous solvent or membrane environments particularly when combined with explicit solvent models. © 2008 Wiley Periodicals, Inc. J Comput Chem, 2009     

An algorithm for the uniform sampling of iso-energy surfaces and for the calculation of microcanonical averages

In this article an algorithm is proposed to efficiently perform the uniform sampling of an iso-energy surface corresponding to a fixed potential energy U of a molecular system, and for calculating averages of certain quantities over microstates having this energy (microcanonical averages). The developed sampling technique is based upon the combination of a recently proposed method for performing constant potential energy molecular dynamics simulations [Rapallo, A. J Chem Phys 2004, 121, 4033] with well-established thermostatting techniques used in the framework of standard molecular dynamics simulations, such as the Andersen thermostat, and the Nose-Hoover chain thermostat. The proposed strategy leads to very accurate and drift-free potential energy conservation during the whole sampling process, and, very important, specially when dealing with high-dimensional or complicated potential functions, it does not require the calculation of the potential energy function hessian. The technique proved to be very reliable for sampling both low- and high-dimensional surfaces.