Replica exchange with dynamical scaling

Calculations of Absolute Free Energies, Enthalpies, and Entropies for Drug Binding

Computer simulation methods can aid in the rational design of drugs aimed at a specific target, typically a protein. The affinity of a drug for its target is given by the free energy of binding. Binding can be further characterized by the enthalpy and entropy changes in the process. Methods exist to determine exact free energies, enthalpies, and entropies that are dependent only on the quality of the potential model and adequate sampling of conformational degrees of freedom. Entropy and enthalpy are roughly an order of magnitude more difficult to calculate than the free energy. This project combines a replica exchange method for enhanced sampling, designed to be efficient for protein-sized systems, with free energy calculations. This approach, replica exchange with dynamical scaling (REDS), uses two conventional simulations at different temperatures so that the entropy can be found from the temperature dependence of the free energy. A third replica is placed between them, with a modified Hamiltonian that allows it to span the temperature range of the conventional replicas. REDS provides temperature-dependent data and aids in sampling. It is applied to the bromodomain-containing protein 4 (BRD4) system. We find that for the force fields used, the free energies are accurate but the entropies and enthalpies are not, with the entropic contribution being too positive. Reproducing the entropy and enthalpy of binding appears to be a more stringent test of the force fields than reproducing the free energy.

The pH Response of a Peptoid Oligomer

Polypeptoids are N-substituted glycine polymers, which differ from peptides in the placement of the side chain on the amide nitrogen rather than the C_α carbon. A peptoid with a chiral side chain containing both an aromatic group and carboxylic acid has a structure that responds to pH changes. All-atom molecular dynamics simulations using a force field specifically tuned for peptoids were carried out with an advanced sampling method for the peptoid (S)-N-(1-carboxy-2-phenylethyl)glycine in the high and low pH limits. The simulations show that the structure changes from mostly cis amide bonds at low pH to mostly trans bonds at high pH. The structural changes are driven by side chain-backbone hydrogen bonds at low pH and side chain repulsions and increased water contact at high pH.

Simulation studies of polypeptoids using replica exchange with dynamical scaling and dihedral biasing

Polypeptoids differ from polypeptides in that the amide bond can more frequently adopt both cis and trans conformations. The transition between the two conformations requires overcoming a large energy barrier, making it difficult for conventional molecular simulations to adequately visit the cis and trans structures. A replica-exchange method is presented that allows for easy rotations of the amide bond and also an efficient linking to a high temperature replica. The method allows for just three replicas (one at the temperature and Hamiltonian of interest, a second high temperature replica with a biased dihedral potential, and a third connecting them) to overcome the amide bond sampling problem and also enhance sampling for other coordinates. The results indicate that for short peptoid oligomers, the conformations can range from all cis to all trans with an average cis/trans ratio that depends on side chain and potential model.

Molecular dynamics simulations of the flexibility and inhibition of SARS-CoV-2 NSP 13 helicase

The helicase protein of the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) is both a good potential drug target and very flexible. The flexibility, and therefore its function, could be reduced through knowledge of these motions and identification of allosteric pockets. Using molecular dynamics simulations with enhanced sampling, we determined key modes of motion and sites on the protein that are at the interface between flexible domains of the proteins. We developed an approach to map the principal components of motion onto the surface of a potential binding pocket to help in the identification of allosteric sites.

Atomistic Simulations of Thermal Unfolding

This tutorial will provide a practical overview of the use of atomistic simulations to study thermal unfolding of biomolecules, in particular small proteins and RNA oligomers. The tutorial focuses on the use of atomistic, all atom simulations of biomolecules in explicit solvent, to study (reversible) thermal unfolding. The simulation methods described here have also been applied to study biomolecules using implicit solvent and coarse-grained models. We do not intend to provide an up-to-date review of the vast literature of biomolecular dynamics, enhanced sampling methods, force field developments, and applications of these methods. The purpose of this tutorial is to provide basic guidelines into the use of these methods to the starting scientist.

Insights into the Thermal Response of a Poly(ethylene oxide)-poly(propylene oxide)-poly(ethylene oxide) Triblock Polymer in Water

A thermal responsive block copolymer made up of ethylene oxide (EO) and propylene oxide (PO) blocks was simulated with optimized atomistic potentials and enhanced sampling methods over a range of temperatures. The results for the L42 pluronic polymer (EO)₄(PO)₂₂(EO)₄, which is known to undergo a transition in this temperature range, and the similarly sized (EO)₃₀ polymer, which does not, are compared. The thermal responsive L42 polymers in a dilute solution tend to aggregate, and this tendency gets stronger as temperature increases. The poly(ethylene oxide) polymer shows no such tendency. The aggregation is stabilized by the hydrophobic contact of the propylene oxide methyl groups, which outweighs a small loss in hydrogen bonds between the ether oxygens and water.

An Implementation of Replica Exchange with Dynamical Scaling for Efficient Large-Scale Simulations

An implementation of the replica exchange with dynamical scaling (REDS) method in the commonly used molecular dynamics program GROMACS is presented. REDS is a replica exchange method that requires fewer replicas than conventional replica exchange while still providing data over a range of temperatures and can be used in either constant volume or constant pressure ensembles. Details for running REDS simulations are given, and an application to the human islet amyloid polypeptide (hIAPP) 11-25 fragment shows that the model efficiently samples conformational space.

A Hybrid Hamiltonian for the Accelerated Sampling along Experimental Restraints

In this article, we present an enhanced sampling method based on a hybrid Hamiltonian which combines experimental distance restraints with a bias dependent from multiple path-dependent variables. This simulation method determines the bias-coordinates on the fly and does not require a priori knowledge about reaction coordinates. The hybrid Hamiltonian accelerates the sampling of proteins, and, combined with experimental distance information, the technique considers the restraints adaptively and in dependency of the system's intrinsic dynamics. We validate the methodology on the dipole relaxation of two water models and the conformational landscape of dialanine. Using experimental NMR-restraint data, we explore the folding landscape of the TrpCage mini-protein and in a second example apply distance restraints from chemical crosslinking/mass spectrometry experiments for the sampling of the conformation space of the Killer Cell Lectin-like Receptor Subfamily B Member 1A (NKR-P1A). The new methodology has the potential to adaptively introduce experimental restraints without affecting the conformational space of the system along an ergodic trajectory. Since only a limited number of input- and no-order parameters are required for the setup of the simulation, the method is broadly applicable and has the potential to be combined with coarse-graining methods.

Enriched Conformational Sampling of DNA and Proteins with a Hybrid Hamiltonian Derived from the Protein Data Bank

In this article, we present a method for the enhanced molecular dynamics simulation of protein and DNA systems called potential of mean force (PMF)-enriched sampling. The method uses partitions derived from the potentials of mean force, which we determined from DNA and protein structures in the Protein Data Bank (PDB). We define a partition function from a set of PDB-derived PMFs, which efficiently compensates for the error introduced by the assumption of a homogeneous partition function from the PDB datasets. The bias based on the PDB-derived partitions is added in the form of a hybrid Hamiltonian using a renormalization method, which adds the PMF-enriched gradient to the system depending on a linear weighting factor and the underlying force field. We validated the method using simulations of dialanine, the folding of TrpCage, and the conformational sampling of the Dickerson⁻Drew DNA dodecamer. Our results show the potential for the PMF-enriched simulation technique to enrich the conformational space of biomolecules along their order parameters, while we also observe a considerable speed increase in the sampling by factors ranging from 13.1 to 82. The novel method can effectively be combined with enhanced sampling or coarse-graining methods to enrich conformational sampling with a partition derived from the PDB.

Replica-Based Protein Structure Sampling Methods: Compromising between Explicit and Implicit Solvents

The structure of a protein is often not completely accessible by experiments. In silico, replica exchange molecular dynamics (REMD) is the standard sampling method for predicting the secondary and tertiary structures from the amino acid sequence, but it is computationally very expensive. Two recent adaptations from REMD, temperature intervals with global exchange of replicas (TIGER2) and TIGER2A, have been tested here in implicit and explicit solvents. Additionally, explicit, implicit, and hybrid solvent REMD are compared. On the basis of the hybrid REMD (REMDh) method, we present a new hybrid TIGER2h algorithm for faster structural sampling, while retaining good accuracy. The implementations of REMDh, TIGER2, TIGER2A, and TIGER2h are provided for nanoscale molecular dynamics (NAMD). All the methods were tested with two model peptides of known structure, (AAQAA)₃ and HP7, with helix and sheet motifs, respectively. The TIGER2 methods and REMDh were also applied to the unknown structure of the collagen type I telopeptides, which represent bigger proteins with some degree of disorder. We present simulations covering more than 180 μs and analyze the performance and convergence of the distributions of states between the particular methods by dihedral principal component and secondary structure analysis.

Toward structure prediction of cyclic peptides

Cyclic peptides are a promising class of molecules that can be used to target specific protein-protein interactions. A computational method to accurately predict their structures would substantially advance the development of cyclic peptides as modulators of protein-protein interactions. Here, we develop a computational method that integrates bias-exchange metadynamics simulations, a Boltzmann reweighting scheme, dihedral principal component analysis and a modified density peak-based cluster analysis to provide a converged structural description for cyclic peptides. Using this method, we evaluate the performance of a number of popular protein force fields on a model cyclic peptide. All the tested force fields seem to over-stabilize the α-helix and PPII/β regions in the Ramachandran plot, commonly populated by linear peptides and proteins. Our findings suggest that re-parameterization of a force field that well describes the full Ramachandran plot is necessary to accurately model cyclic peptides.

Coarse kMC-based replica exchange algorithms for the accelerated simulation of protein folding in explicit solvent

In this paper, we present a coarse replica exchange molecular dynamics (REMD) approach, based on kinetic Monte Carlo (kMC).

Superposition-Enhanced Estimation of Optimal Temperature Spacings for Parallel Tempering Simulations

Effective parallel tempering simulations rely crucially on a properly chosen sequence of temperatures. While it is desirable to achieve a uniform exchange acceptance rate across neighboring replicas, finding a set of temperatures that achieves this end is often a difficult task, in particular for systems undergoing phase transitions. Here we present a method for determination of optimal replica spacings, which is based upon knowledge of local minima in the potential energy landscape. Working within the harmonic superposition approximation, we derive an analytic expression for the parallel tempering acceptance rate as a function of the replica temperatures. For a particular system and a given database of minima, we show how this expression can be used to determine optimal temperatures that achieve a desired uniform acceptance rate. We test our strategy for two atomic clusters that exhibit broken ergodicity, demonstrating that our method achieves uniform acceptance as well as significant efficiency gains.

Optimization of Umbrella Sampling Replica Exchange Molecular Dynamics by Replica Positioning

The positioning of sampling windows in an umbrella sampling simulation has an effect on the rate of convergence and computational efficiency. When such simulation is coupled with a Hamiltonian replica exchange setup, we show that such positioning can be optimized for maximal convergence of the results. We present a method for estimating the exchange acceptance ratio (EAR) between two arbitrary positions on a reaction coordinate in umbrella sampling replica exchange (USRE) molecular dynamics (MD). We designed a scoring function to optimize the position of the set of replicas (windows). By maximizing the scoring function, we make EAR the same for all neighbor replica pairs, increasing the efficiency of the method. We tested our algorithm by sampling a torsion for butane in implicit solvent and by studying a salt bridge in explicit solvent. We found that the optimized set of replicas recovers the correct free energy profile much faster than for equally spaced umbrellas.

Generating reservoir conformations for replica exchange through the use of the conformational space annealing method

Temperature replica exchange molecular dynamics (T-REM) has been successfully used to improve the conformational search for model peptides and small proteins. However, for larger and more complicated systems, the use of T-REM is computationally intensive since the complexity of the free energy landscape and number of replicas required increase with system size. Achieving convergence of systems with slow transition kinetics is often difficult. Several methods have been proposed to overcome the size and convergence speed issues of standard T-REM. One of these is the Reservoir Replica Exchange Method (R-REM), in which the conformational search and temperature equilibration are separated by exchanging with a pre-existing reservoir of structures. This approach allows the integration of computationally efficient search algorithms with replica exchange. The Conformational Space Annealing (CSA) method has been shown to be able to determine the global energy minimum of proteins efficiently and has been used in structure prediction successfully. CSA uses a genetic algorithm to generate a diverse set of conformations to determine the minimum energy structure. We combine these methods by using conformations generated by the CSA method to build a reservoir. R-REM is then used to seed the top replica with the structures from the reservoir; fast convergence at every temperature is observed. The efficiency of this method is then demonstrated with model peptides and small proteins, and significant improvement of efficiency is observed while maintaining the overall shape of the free energy landscape.

Replica exchanging self-guided Langevin dynamics for efficient and accurate conformational sampling

This work presents a replica exchanging self-guided Langevin dynamics (RXSGLD) simulation method for efficient conformational searching and sampling. Unlike temperature-based replica exchanging simulations, which use high temperatures to accelerate conformational motion, this method uses self-guided Langevin dynamics (SGLD) to enhance conformational searching without the need to elevate temperatures. A RXSGLD simulation includes a series of SGLD simulations, with simulation conditions differing in the guiding effect and/or temperature. These simulation conditions are called stages and the base stage is one with no guiding effect. Replicas of a simulation system are simulated at the stages and are exchanged according to the replica exchanging probability derived from the SGLD partition function. Because SGLD causes less perturbation on conformational distribution than high temperatures, exchanges between SGLD stages have much higher probabilities than those between different temperatures. Therefore, RXSGLD simulations have higher conformational searching ability than temperature based replica exchange simulations. Through three example systems, we demonstrate that RXSGLD can generate target canonical ensemble distribution at the base stage and achieve accelerated conformational searching. Especially for large systems, RXSGLD has remarkable advantages in terms of replica exchange efficiency, conformational searching ability, and system size extensiveness.

MDMS: Molecular Dynamics Meta-Simulator for evaluating exchange type sampling methods

Replica exchange methods have become popular tools to explore conformational space for small proteins. For larger biological systems, even with enhanced sampling methods, exploring the free energy landscape remains computationally challenging. This problem has led to the development of many improved replica exchange methods. Unfortunately, testing these methods remains expensive. We propose a Molecular Dynamics Meta-Simulator (MDMS) based on transition state theory to simulate a replica exchange simulation, eliminating the need to run explicit dynamics between exchange attempts. MDMS simulations allow for rapid testing of new replica exchange based methods, greatly reducing the amount of time needed for new method development.

Replica exchange statistical temperature molecular dynamics algorithm

The replica exchange statistical temperature molecular dynamics (RESTMD) algorithm is presented, designed to alleviate an extensive increase of the number of replicas required as system size increases in the conventional temperature replica exchange method (tREM), and to obtain improved sampling in individual replicas. RESTMD optimally integrates multiple STMD (Phys. Rev. Lett. 2006, 97, 050601) runs with replica exchanges, giving rise to a flat energy sampling in each replica with a self-adjusting weight determination. The expanded flat energy dynamic sampling range allows the use of significantly fewer STMD replicas while maintaining the desired acceptance probability for replica exchanges. The computational advantages of RESTMD over conventional REM and single-replica STMD are explicitly demonstrated with an application to a coarse-grained protein model. The effect of two different kinetic temperature control schemes on the sampling efficiency is explored for diverse simulation conditions.

Equilibrium sampling in biomolecular simulations

Equilibrium sampling of biomolecules remains an unmet challenge after more than 30 years of atomistic simulation. Efforts to enhance sampling capability, which are reviewed here, range from the development of new algorithms to parallelization to novel uses of hardware. Special focus is placed on classifying algorithms--most of which are underpinned by a few key ideas--in order to understand their fundamental strengths and limitations. Although algorithms have proliferated, progress resulting from novel hardware use appears to be more clear-cut than from algorithms alone, due partly to the lack of widely used sampling measures.

How hot? Systematic convergence of the replica exchange method using multiple reservoirs

We have devised a systematic approach to converge a replica exchange molecular dynamics simulation by dividing the full temperature range into a series of higher temperature reservoirs and a finite number of lower temperature subreplicas. A defined highest temperature reservoir of equilibrium conformations is used to help converge a lower but still hot temperature subreplica, which in turn serves as the high-temperature reservoir for the next set of lower temperature subreplicas. The process is continued until an optimal temperature reservoir is reached to converge the simulation at the target temperature. This gradual convergence of subreplicas allows for better and faster convergence at the temperature of interest and all intermediate temperatures for thermodynamic analysis, as well as optimizing the use of multiple processors. We illustrate the overall effectiveness of our multiple reservoir replica exchange strategy by comparing sampling and computational efficiency with respect to replica exchange, as well as comparing methods when converging the structural ensemble of the disordered Abeta(21-30) peptide simulated with explicit water by comparing calculated Rotating Overhauser Effect Spectroscopy intensities to experimentally measured values.

Generalized replica exchange method

We present a powerful replica exchange method, particularly suited to first-order phase transitions associated with the backbending in the statistical temperature, by merging an optimally designed generalized ensemble sampling with replica exchanges. The key ingredients of our method are parametrized effective sampling weights, smoothly joining ordered and disordered phases with a succession of unimodal energy distributions by transforming unstable or metastable energy states of canonical ensembles into stable ones. The inverse mapping between the sampling weight and the effective temperature provides a systematic way to design the effective sampling weights and determine a dynamic range of relevant parameters. Illustrative simulations on Potts spins with varying system size and simulation conditions demonstrate a comprehensive sampling for phase-coexistent states with a dramatic acceleration of tunneling transitions. A significant improvement over the power-law slowing down of mean tunneling times with increasing system size is obtained, and the underlying mechanism for accelerated tunneling is discussed.

Replica exchange statistical temperature Monte Carlo

The replica exchange statistical temperature Monte Carlo algorithm (RESTMC) is presented, extending the single-replica STMC algorithm [J. Kim, J. E. Straub, and T. Keyes, Phys. Rev. Lett. 97, 050601 (2006)] to alleviate the slow convergence of the conventional temperature replica exchange method (t-REM) with increasing system size. In contrast to the Gibbs-Boltzmann sampling at a specific temperature characteristic of the standard t-REM, RESTMC samples a range of temperatures in each replica and achieves a flat energy sampling employing the generalized sampling weight, which is automatically determined via the dynamic modification of the replica-dependent statistical temperature. Faster weight determination, through the dynamic update of the statistical temperature, and the flat energy sampling, maximizing energy overlaps between neighboring replicas, lead to a considerable acceleration in the convergence of simulations even while employing significantly fewer replicas. The performance of RESTMC is demonstrated and quantitatively compared with that of the conventional t-REM under varying simulation conditions for Lennard-Jones 19, 31, and 55 atomic clusters, exhibiting single- and double-funneled energy landscapes.

Replica exchange with nonequilibrium switches

We introduce a replica exchange (parallel tempering) method in which attempted configuration swaps are generated using nonequilibrium work simulations. By effectively increasing phase space overlap, this approach mitigates the need for many replicas. We illustrate our method by using a model system and show that it is able to achieve the computational efficiency of ordinary replica exchange, using fewer replicas.

Bad Seeds Sprout Perilous Dynamics: Stochastic Thermostat Induced Trajectory Synchronization in Biomolecules

Molecular dynamics simulations starting from different initial conditions are commonly used to mimic the behavior of an experimental ensemble. We show in this article that when a Langevin thermostat is used to maintain constant temperature during such simulations, extreme care must be taken when choosing the random number seeds to prevent statistical correlation among the MD trajectories. While recent studies have shown that stochastically thermostatted trajectories evolving within a single potential basin with identical random number seeds tend to synchronize, we show that there is a synchronization effect even for complex, biologically relevant systems. We demonstrate this effect in simulations of alanine trimer and pentamer and in a simulation of a temperature-jump experiment for peptide folding of a 14-residue peptide. Even in replica-exchange simulations, in which the trajectories are at different temperatures, we find partial synchronization occurring when the same random number seed is employed. We explain this by extending the recent derivation of the synchronization effect for two trajectories in a harmonic well to the case in which the trajectories are at two different temperatures. Our results suggest several ways in which mishandling selection of a pseudorandom number generator initial seed can lead to corruption of simulation data. Simulators can fall into this trap in simple situations such as neglecting to specifically indicate different random seeds in either parallel or sequential restart simulations, utilizing a simulation package with a weak pseudorandom number generator, or using an advanced simulation algorithm that has not been programmed to distribute initial seeds.

Optimal replica exchange method combined with Tsallis weight sampling

A unified framework integrating the generalized ensemble sampling associated with the Tsallis weight [C. Tsallis, J. Stat. Phys. 52, 479 (1988)] and the replica exchange method (REM) has been proposed to accelerate the convergence of the conventional temperature REM (t-REM). Using the effective temperature formulation of the Tsallis weight sampling, it is shown that the average acceptance probability for configurational swaps between neighboring replicas in the combination of Tsallis weight sampling and REM (Tsallis-REM) is directly proportional to an overlap integral of the energy distributions of neighboring replicas as in the t-REM. Based on this observation, we suggest a robust method to select optimal Tsallis parameters in the conventional parametrization scheme and present new parametrization schemes for the Tsallis-REM, which significantly improves the acceptance of configurational swaps by systematically modulating energy overlaps between neighboring replicas. The distinguished feature of our method is that all relevant parameters in the Tsallis-REM are automatically determined from the equilibrium phase simulation using the t-REM. The overall performance of our method is explicitly demonstrated for various simulation conditions for the Lennard-Jones 31 atom clusters, exhibiting a double-funneled energy landscape.

Exchange frequency in replica exchange molecular dynamics

The effect of the exchange-attempt frequency on sampling efficiency is studied in replica exchange molecular dynamics (REMD). We show that sampling efficiency increases with increasing exchange-attempt frequency. This conclusion is contrary to a commonly expressed view in REMD. Five peptides (1-21 residues long) are studied with a spectrum of exchange-attempt rates. Convergence rates are gauged by comparing ensemble properties between fixed length test REMD simulations and longer reference simulations. To show the fundamental correlation between exchange frequency and convergence time, a simple model is designed and studied, displaying the same basic behavior of much more complex systems.

A temperature predictor for parallel tempering simulations

An algorithm is proposed that generates a set of temperatures for use in parallel tempering simulations (also known as temperature-replica exchange molecular dynamics simulations) of proteins to obtain a desired exchange probability Pdes. The input consists of the number of protein atoms and water molecules in the system, information about the use of constraints and virtual sites and the lower temperature limits. The temperatures generated yield probabilities which are very close to Pdes (correlation 97%), independent of force field and over a wide temperature range. To facilitate its use, the algorithm has been implemented as a web server at .