Adaptive lambda square dynamics simulation: An efficient conformational sampling method for biomolecules

A novel, efficient sampling method for biomolecules is proposed. The partial multicanonical molecular dynamics (McMD) was recently developed as a method that improved generalized ensemble (GE) methods to focus sampling only on a part of a system (GEPS); however, it was not tested well. We found that partial McMD did not work well for polylysine decapeptide and gave significantly worse sampling efficiency than a conventional GE. Herein, we elucidate the fundamental reason for this and propose a novel GEPS, adaptive lambda square dynamics (ALSD), which can resolve the problem faced when using partial McMD. We demonstrate that ALSD greatly increases the sampling efficiency over a conventional GE. We believe that ALSD is an effective method and is applicable to the conformational sampling of larger and more complicated biomolecule systems. © 2013 Wiley Periodicals, Inc.     

Theory for trivial trajectory parallelization of multicanonical molecular dynamics and application to a polypeptide in water

Trivial trajectory parallelization of multicanonical molecular dynamics (TTP-McMD) explores the conformational space of a biological system with multiple short runs of McMD starting from various initial structures. This method simply connects (i.e., trivially parallelizes) the short trajectories and generates a long trajectory. First, we theoretically prove that the simple trajectory connection satisfies a detailed balance automatically. Thus, the resultant long trajectory is regarded as a single multicanonical trajectory. Second, we applied TTP-McMD to an alanine decapeptide with an all-atom model in explicit water to compute a free-energy landscape. The theory imposes two requirements on the multiple trajectories. We have demonstrated that TTP-McMD naturally satisfies the requirements. The TTP-McMD produces the free-energy landscape considerably faster than a single-run McMD does. We quantitatively showed that the accuracy of the computed landscape increases with increasing the number of multiple runs. Generally, the free-energy landscape of a large biological system is unknown a priori. The current method is suitable for conformational sampling of such a large system to reduce the waiting time to obtain a canonical ensemble statistically reliable.     

Exhaustive search of the configurational space of heat-shock protein 90 with its inhibitor by multicanonical molecular dynamics based dynamic docking

Multicanonical molecular dynamics based dynamic docking was used to exhaustively search the configurational space of an inhibitor binding to the N-terminal domain of heat-shock protein 90 (Hsp90). The obtained structures at 300 K cover a wide structural ensemble, with the top two clusters ranked by their free energy coinciding with the native binding site. The representative structure of the most stable cluster reproduced the experimental binding configuration, but an interesting conformational change in Hsp90 could be observed. The combined effects of solvation and ligand binding shift the equilibrium from a preferred loop-in conformation in the unbound state to an α-helical one in the bound state for the flexible lid region of Hsp90. Thus, our dynamic docking method is effective at predicting the native binding site while exhaustively sampling a wide configurational space, modulating the protein structure upon binding.     

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.     

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.     

Free Energy Calculations using a Swarm‐Enhanced Sampling Molecular Dynamics Approach

Free energy simulations are an established computational tool in modelling chemical change in the condensed phase. However, sampling of kinetically distinct substates remains a challenge to these approaches. As a route to addressing this, we link the methods of thermodynamic integration (TI) and swarm‐enhanced sampling molecular dynamics (sesMD), where simulation replicas interact cooperatively to aid transitions over energy barriers. We illustrate the approach by using alchemical alkane transformations in solution, comparing them with the multiple independent trajectory TI (IT‐TI) method. Free energy changes for transitions computed by using IT‐TI grew increasingly inaccurate as the intramolecular barrier was heightened. By contrast, swarm‐enhanced sampling TI (sesTI) calculations showed clear improvements in sampling efficiency, leading to more accurate computed free energy differences, even in the case of the highest barrier height. The sesTI approach, therefore, has potential in addressing chemical change in systems where conformations exist in slow exchange.     

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     

Variation of free‐energy landscape of the p53 C‐terminal domain induced by acetylation: Enhanced conformational sampling

The C‐terminal domain (CTD) of tumor suppressor protein p53 is an intrinsically disordered region that binds to various partner proteins, where lysine of CTD is acetylated/nonacetylated and histidine neutralized/non‐neutralized. Because of the flexibility of the unbound CTD, a free‐energy landscape (FEL) is a useful quantity for determining its statistical properties. We conducted enhanced conformational sampling of CTD in the unbound state via virtual system coupled multicanonical molecular dynamics, in which the lysine was acetylated or nonacetylated and histidine was charged or neutralized. The fragments were expressed by an all‐atom model and were immersed in an explicit solvent. The acetylation and charge‐neutralization varied FEL greatly, which might be convenient to exert a hub property. The acetylation slightly enhanced alpha‐helix structures that are more compact than sheet/loop conformations. The charge‐neutralization produced hairpins. Additionally, circular dichroism experiments confirmed the computational results. We propose possible binding mechanisms of CTD to partners by investigating FEL. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

FSATOOL: A useful tool to do the conformational sampling and trajectory analysis work for biomolecules

Reliable conformational sampling and trajectory analysis are always important to the study of the folding or binding mechanisms of biomolecules. Generally, one has to prepare many complicated parameters and follow a lot of steps to obtain the final data. The whole process is too complicated to new users. In this article, we provide a convenient and user-friendly tool that is compatible to AMBER, called fast sampling and analysis tool (FSATOOL). FSATOOL has some useful features. First and the most important, the whole work is extremely simplified into two steps, one is the fast sampling procedure and the other is the trajectory analysis procedure. Second, it contains several powerful sampling methods for the simulation on graphics process unit, including our previous mixing replica exchange molecular dynamics method. The method combines the advantages of the biased and unbiased simulations. Finally, it extracts the dominant transition pathways automatically from the folding network by Markov state model. Users do not need to do the tedious intermediate steps by hand. To illustrate the usage of FSATOOL in practice, we perform one simulation for a RNA hairpin in explicit solvent. All the results are presented. © 2019 Wiley Periodicals, Inc.     

Assessment of enveloping distribution sampling to calculate relative free enthalpies of binding for eight netropsin-DNA duplex complexes in aqueous solution

The performance of enveloping distribution sampling (EDS) simulations to estimate free enthalpy differences associated with seven alchemical transformations of A-T into G-C base pairs at the netropsin binding site in the minor groove of a 13-base pair DNA duplex in aqueous solution is evaluated. It is demonstrated that sufficient sampling can be achieved with a two-state EDS Hamiltonian even for large perturbations such as the simultaneous transformation of up to three A-T into three G-C base pairs. The two parameters required to define the EDS reference state Hamiltonian are obtained automatically using a modified version of a scheme presented in earlier work. The sensitivity of the configurational sampling to a variation of these parameters is investigated in detail. Although for relatively small perturbations, that is, one base pair, the free enthalpy estimate depends only weakly on the EDS parameters, the sensitivity is stronger for the largest perturbation. Yet, EDS offers various convenient measures to evaluate the degree of sampling and thus the reliability of the free enthalpy estimate and appears to be an efficient alternative to the conventional thermodynamic integration methodology to obtain free energy differences for molecular systems.     

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.     

Branch migration of Holliday junction in RuvA tetramer complex studied by umbrella sampling simulation using a path‐search algorithm

Branch migration of the Holliday junction takes place at the center of the RuvA tetramer. To elucidate how branch migration occurs, umbrella sampling simulations were performed for complexes of the RuvA tetramer and Holliday junction DNA. Although conventional umbrella sampling simulations set sampling points a priori, the umbrella sampling simulation in this study set the sampling points one by one in order to search for a realistic path of the branch migration during the simulations. Starting from the X‐ray structure of the complex, in which the hydrogen bonds between two base‐pairs were unformed, the hydrogen bonds between the next base‐pairs of the shrinking stems were observed to start to disconnect. At the intermediate stage, three or four of the eight unpaired bases interacted closely with the acidic pins from RuvA. During the final stage, these bases moved away from the pins and formed the hydrogen bonds of the new base‐pairs of the growing stems. The free‐energy profile along this reaction path showed that the intermediate stage was a meta‐stable state between two free‐energy barriers of about 10 to 15 kcal/mol. These results imply that the pins play an important role in stabilizing the interactions between the pins and the unpaired base‐pairs. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Replica state exchange metadynamics for improving the convergence of free energy estimates

Metadynamics (MTD) is a powerful enhanced sampling method for systems with rugged energy landscapes. It constructs a bias potential in a predefined collective variable (CV) space to overcome barriers between metastable states. In bias‐exchange MTD (BE‐MTD), multiple replicas approximate the CV space by exchanging bias potentials (replica conditions) with the Metropolis–Hastings (MH) algorithm. We demonstrate that the replica‐exchange rates and the convergence of free energy estimates of BE‐MTD are improved by introducing the infinite swapping (IS) or the Suwa‐Todo (ST) algorithms. Conceptually, IS and ST perform transitions in a replica state space rather than exchanges in a replica condition space. To emphasize this, the proposed scheme is called the replica state exchange MTD (RSE‐MTD). Benchmarks were performed with alanine polypeptides in vacuum and water. For the systems tested in this work, there is no significant performance difference between IS and ST. © 2015 Wiley Periodicals, Inc.     

Orthogonal sampling in free‐energy calculations of residue mutations in a tripeptide: TI versus λ‐LEUS

In a recent article (Bieler et al., J. Chem. Theory Comput. 2014, 10, 3006), we introduced a combination of λ‐dynamics and local‐elevation umbrella‐sampling termed λ‐LEUS to calculate free‐energy changes associated with alchemical processes using molecular dynamics simulations. This method was suggested to be more efficient than thermodynamic integration (TI), because the dynamical variation of the alchemical variable λ opens up pathways to circumvent barriers in the orthogonal space (defined by the N – 1 degrees of freedom that are not subjected to the sampling enhancement), a feature λ‐LEUS shares with Hamiltonian replica‐exchange (HR) approaches. However, the mutation considered, hydroquinone to benzene in water, was no real challenge in terms of orthogonal‐space properties, which were restricted to solvent‐relaxation processes. In the present article, we revisit the comparison between TI and λ‐LEUS considering non‐trivial mutations of the central residue X of a KXK tripeptide in water (with X = G, E, K, S, F, or Y). Side‐chain interactions that may include salt bridges, hydrogen bonds or steric clashes lead to slow relaxation in the orthogonal space, mainly in the two‐dimensional subspace spanned by the central and ψ dihedral angles of the peptide. The efficiency enhancement afforded by λ‐LEUS is confirmed in this more complex test system and can be attributed explicitly to the improved sampling of the orthogonal space. The sensitivity of the results to the nontrivial choices of a mass parameter and of a thermostat coupling time for the alchemical variable is also investigated, resulting in recommended ranges of 50 to 100 u nm² and 0.2 to 0.5 ps, respectively. © 2015 Wiley Periodicals, Inc.     

Fragment Molecular Orbital method-based Molecular Dynamics (FMO-MD) as a simulator for chemical reactions in explicit solvation

Fragment Molecular Orbital based-Molecular Dynamics (FMO-MD, Komeiji et al., Chem Phys Lett 2003, 372, 342) is an ab initio MD method suitable for large molecular systems. Here, FMO-MD was implemented to conduct full quantum simulations of chemical reactions in explicit solvation. Several FMO-MD simulations were performed for a sphere of water to find a suitable simulation protocol. It was found that annealing of the initial configuration by a classical MD brought the subsequent FMO-MD trajectory to faster stabilization, and also that use of bond constraint in the FMO-MD heating stage effectively reduced the computation time. Then, the blue moon ensemble method (Sprik and Ciccotti, J Chem Phys 1998, 109, 7737) was implemented and was tested by calculating free energy profiles of the Menschutkin reaction (H3N + CH3Cl --> +H3NCH3 + Cl-) in the presence and absence of the solvent water via FMO-MD. The obtained free energy profiles were consistent with the Hammond postulate in that stabilization of the product by the solvent, namely hydration of Cl-, shifted the transition state to the reactant-side. Based on these FMO-MD results, plans for further improvement of the method are discussed.     

Detailed considerations for a balanced and broadly applicable force field: a study of substituted benzenes modeled with OPLS-AA

Modern classical force fields have been traditionally parameterized by attempting to maximize agreement to any number of experimental and/or quantum mechanical target properties. As these force fields are pushed towards obtaining quantitative estimates of often subtle energetic differences, stringent and consistent parameterization criteria, particularly in regard to charge distributions, are required to ensure that systematic errors cancel, that parameters are transferable between molecules, and that performance does not significantly deteriorate when using more approximate methods, such as with continuum solvent models. Relative free energies of hydration are presented here for 40 mono- and disubstituted benzenes modeled with the OPLS-AA force field; heats of vaporization and pure liquid densities at standard conditions are presented when experimental data is available. Overall agreement between OPLS-AA and experiment is remarkable (average error = 0.5 kcal/mol for DeltaDeltaG(hydration), 1.0 kcal/mol for DeltaH(vap) (0), 0.02 g/mL for densities), yet several functional groups are identified as having consistent and correctable errors (alkyl-, nitro-, and thiobenzenes). Relative free energies of hydration obtained with rigorous free energy perturbations using explicit solvent are also compared with energies from minimizations using a generalized Born model (GB). There is high correlation between these estimates (R = 0.99), and as demonstrated here, reparameterization of the aforementioned groups can be guided with rapid GB calculations.