On easy implementation of a variant of the replica exchange with solute tempering in GROMACS

To reduce the number of replicas required in the conventional replica exchange method for huge systems, recently the replica exchange with solute tempering (REST) method was proposed. Here we showed that a variant of REST realized by rescaling the force-field parameters can be performed with GROMACS 4 without changing the code. We tested the variant REST for alanine dipeptide and an N-terminal peptide from p53 confirming its performance nearly equal to the original REST.     

Multiple scaling replica exchange for the conformational sampling of biomolecules in explicit water

A multiple scaling replica exchange method for the efficient conformational sampling of biomolecular systems in explicit solvent is presented. The method is a combination of the replica exchange with solute tempering (REST) technique and a Tsallis biasing potential. The Tsallis biasing increases the sampling efficiency, while the REST minimizes the number of replicas needed. Unbiased statistics can be obtained by reweighting of the data using a weighted histogram analysis technique. The method is illustrated by its application to a ten residue peptide in explicit water.     

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.     

Serial replica exchange

Parallel tempering (or the replica exchange method (REM)) is a powerful method for speeding up the sampling of conformational states of systems with rough energy landscapes, like proteins, where stable conformational states can be separated by large energy barriers. The usual implementation of the REM is performed on local computer clusters (or parallel processors) where the different replicas must be run synchronously. Here, we present serial replica exchange (SREM), a method that is equivalent to the standard REM in terms of efficiency yet runs asynchronously on a distributed network of computers. A second advantage is the method's greatly enhanced fault tolerance, which enables the study of biological systems on worldwide distributed computing environments, such as Folding@Home. For proof of concept, we apply the SREM to a single alanine dipeptide molecule in explicit water. We show that the SREM reproduces the thermodynamic and structural properties determined by the REM.     

Modified replica exchange simulation methods for local structure refinement

Parallel tempering, also known as replica exchange molecular dynamics (REMD), has recently been successfully used to study the structure and thermodynamic properties of biomolecules such as peptides and small proteins. For large systems, however, applying REMD can be costly since the number of replicas needed increases as the square root of the number of degrees of freedom in the system. Often, enhanced sampling is only needed for a subset of atoms, such as a loop region of a large protein or a small ligand binding to a receptor. In such applications, it is often reasonable to assume a weak dependence of the structure of the larger region on the instantaneous conformation of the smaller region of interest. For these cases, we derived two variant replica exchange methods, partial replica exchange molecular dynamics (PREMD) and local replica exchange molecular dynamics (LREMD). The Hamiltonian for the system is separated, with replica exchange carried out only for terms involving the subsystem of interest while the remainder of the system is maintained at a single temperature. The number of replicas required for efficient exchange thus depends on the number of degrees of freedom in the fragment needing refinement rather than on the size of the full system. The method can be applied to much larger systems than was previously practical. This also provides a means to preserve the integrity of the structure outside the refinement region without introduction of restraints. LREMD takes this weak coupling approximation a step further, employing only a single representation of the large fragment that simultaneously interacts with all of the replicas of the subsystem of interest. This is obtained by combining replica exchange with the locally enhanced sampling approximation (LES), reducing the computational expense of replica exchange simulations to near that of a single standard molecular dynamics (MD) simulation. Use of LREMD also permits the use of LES without requiring the specification of a single temperature, a known difficulty for standard LES simulations. We tested these two methods on the loop region of an RNA hairpin model system and find significant advantages over standard MD and REMD simulations.     

Catalytic mechanism of RNA backbone cleavage by ribonuclease H from quantum mechanics/molecular mechanics simulations

We use quantum mechanics/molecular mechanics simulations to study the cleavage of the ribonucleic acid (RNA) backbone catalyzed by ribonuclease H. This protein is a prototypical member of a large family of enzymes that use two-metal catalysis to process nucleic acids. By combining Hamiltonian replica exchange with a finite-temperature string method, we calculate the free energy surface underlying the RNA-cleavage reaction and characterize its mechanism. We find that the reaction proceeds in two steps. In a first step, catalyzed primarily by magnesium ion A and its ligands, a water molecule attacks the scissile phosphate. Consistent with thiol-substitution experiments, a water proton is transferred to the downstream phosphate group. The transient phosphorane formed as a result of this nucleophilic attack decays by breaking the bond between the phosphate and the ribose oxygen. In the resulting intermediate, the dissociated but unprotonated leaving group forms an alkoxide coordinated to magnesium ion B. In a second step, the reaction is completed by protonation of the leaving group, with a neutral Asp132 as a likely proton donor. The overall reaction barrier of ～15 kcal mol(-1), encountered in the first step, together with the cost of protonating Asp132, is consistent with the slow measured rate of ～1-100/min. The two-step mechanism is also consistent with the bell-shaped pH dependence of the reaction rate. The nonmonotonic relative motion of the magnesium ions along the reaction pathway agrees with X-ray crystal structures. Proton-transfer reactions and changes in the metal ion coordination emerge as central factors in the RNA-cleavage reaction.     

Hybrid Hamiltonian replica exchange molecular dynamics simulation method employing the Poisson-Boltzmann model

A hybrid Hamiltonian replica exchange molecular dynamics simulation scheme based on explicit water model hybrided with Poisson-Boltzmann model is brought out. In this method the motions of atoms are governed by potential energy obtained from explicit water model. However, the exchanges between different replicas under different temperatures are controlled by the solvation energies of the solute calculated using the Poisson-Boltzmann model. In order to get the correct canonical ensembles, the van der Waals radii, which are used to define the dielectric boundary, have to be optimized. The conformational spaces of three distinct pentapeptides, Met-enkephalin, alanine 5, and glycine 5, are explored. We find that with the optimized radii the structural ensembles are nearly identical to those obtained by standard replica exchange simulations while the number of replica needed is reduced greatly.     

Effect of Cholesterol Molecules on Aβ1-42 Wild-Type and Mutants Trimers

Alzheimer’s disease displays aggregates of the amyloid-beta (Aβ) peptide in the brain, and there is increasing evidence that cholesterol may contribute to the pathogenesis of the disease. Though many experimental and theoretical studies have focused on the interactions of Aβ oligomers with membrane models containing cholesterol, an understanding of the effect of free cholesterol on small Aβ42 oligomers is not fully established. To address this question, we report on replica exchange with a solute tempering simulation of an Aβ42 trimer with cholesterol and compare it with a previous replica exchange molecular dynamics simulation. We show that the binding hot spots of cholesterol are rather complex, involving hydrophobic residues L17–F20 and L30–M35 with a non-negligible contribution of loop residues D22–K28 and N-terminus residues. We also examine the effects of cholesterol on the trimers of the disease-causing A21G and disease-protective A2T mutations by molecular dynamics simulations. We show that these two mutations moderately impact cholesterol-binding modes. In our REST2 simulations, we find that cholesterol is rarely inserted into aggregates but rather attached as dimers and trimers at the surface of Aβ42 oligomers. We propose that cholesterol acts as a glue to speed up the formation of larger aggregates; this provides a mechanistic link between cholesterol and Alzheimer’s disease.     

Replica exchange with dynamical scaling

A replica exchange method is presented which requires fewer replicas and is designed to be used for large systems. In this method, dynamically scaled replicas are placed between conventional replicas at broadly spaced temperatures. The potential of the scaled replicas is linearly scaled by a dynamical variable which varies between 0 and 1. When the variable is near either end point the replica can undergo exchanges with one of its neighboring replicas. Two different versions of the method are presented for a model system of a small peptide in water. The scaled replica can replace many replicas and the method can be up to ten times more efficient than conventional replica exchange.     

Simulation of the pressure and temperature folding/unfolding equilibrium of a small RNA hairpin

We report molecular dynamics simulations of the equilibrium folding/unfolding thermodynamics of the RNA tetraloop in explicit solvent. A replica exchange molecular dynamics study of the r(CGUUGCCG) oligomer that forms a hairpin is performed for 226 ns per replica, using 52 replicas. We are able to show the unbiased folding of all replicas starting from extended conformations. The equilibrium pressure-temperature free energy of folding, DeltaG(P,T), is calculated from the averaged energy, pressure, and specific volume change upon folding of the oligomer as a function of T at constant volume. We find that this oligomer is destabilized by increasing hydrostatic pressure, similar to the behavior of globular proteins.     

Hamiltonian replica exchange molecular dynamics using soft-core interactions

To overcome the problem of insufficient conformational sampling within biomolecular simulations, we have developed a novel Hamiltonian replica exchange molecular dynamics (H-REMD) scheme that uses soft-core interactions between those parts of the system that contribute most to high energy barriers. The advantage of this approach over other H-REMD schemes is the possibility to use a relatively small number of replicas with locally larger differences between the individual Hamiltonians. Because soft-core potentials are almost the same as regular ones at longer distances, most of the interactions between atoms of perturbed parts will only be slightly changed. Rather, the strong repulsion between atoms that are close in space, which in many cases results in high energy barriers, is weakened within higher replicas of our proposed scheme. In addition to the soft-core interactions, we proposed to include multiple replicas using the same Hamiltonian/level of softness. We have tested the new protocol on the GTP and 8-Br-GTP molecules, which are known to have high energy barriers between the anti and syn conformation of the base with respect to the sugar moiety. During two 25 ns MD simulations of both systems the transition from the more stable to the less stable (but still experimentally observed) conformation is not seen at all. Also temperature REMD over 50 replicas for 1 ns did not show any transition at room temperature. On the other hand, more than 20 of such transitions are observed in H-REMD using six replicas (at three different Hamiltonians) during 6.8 ns per replica for GTP and 12 replicas (at six different Hamiltonians) during 8.7 ns per replica for 8-Br-GTP. The large increase in sampling efficiency was obtained from an optimized H-REMD scheme involving soft-core potentials, with multiple simulations using the same level of softness. The optimization of the scheme was performed by fast mimicking [J. Hritz and C. Oostenbrink, J. Chem. Phys. 127, 204104 (2007)].     

Protein-ligand binding affinity by nonequilibrium free energy methods

Nonequilibrium (NE) free energy methods are embarrassingly parallel and may be very conveniently run on desktop computers using distributed computing software. In recent years there has been a proliferation of NE methods, but these approaches have barely, if at all, been used in the context of calculating protein-ligand binding free energies. In a recent study by these authors, different combinations of NE methods with various test systems were compared and protocols identified which yielded results as accurate as replica exchange thermodynamic integration (RETI). The NE approaches, however, lend themselves to extensive parallelization through the use of distributed computing. Here the best performing of those NE protocols, a replica exchange method using Bennett's acceptance ratio as the free energy estimator (RENE), is applied to two sets of congeneric inhibitors bound to neuraminidase and cyclooxygenase-2. These protein-ligand systems were originally studied with RETI, giving results to which NE and RENE simulations are compared. These NE calculations were carried out on a large, highly distributed group of low-performance desktop computers which are part of a Condor pool. RENE was found to produce results of a predictive quality at least as good as RETI in less than half the wall clock time. However, non-RE NE results were found to be far less predictive. In addition, the RENE method successfully identified a localized region of rapidly changing free energy gradients without the need for prior investigation. These results suggest that the RENE protocol is appropriate for use in the context of predicting protein-ligand binding free energies and that it can offer advantages over conventional, equilibrium approaches.     

Two-dimensional replica exchange approach for peptide-peptide interactions

The replica exchange molecular dynamics (REMD) method has emerged as a standard approach for simulating proteins and peptides with rugged underlying free energy landscapes. We describe an extension to the original methodology--here termed umbrella-sampling REMD (UREMD)--that offers specific advantages in simulating peptide-peptide interactions. This method is based on the use of two dimensions in the replica cascade, one in temperature as in conventional REMD, and one in an umbrella sampling coordinate between the center of mass of the two peptides that aids explicit exploration of the complete association-dissociation reaction coordinate. To mitigate the increased number of replicas required, we pursue an approach in which the temperature and umbrella dimensions are linked at only fully associated and dissociated states. Coupled with the reweighting equations, the UREMD method aids accurate calculations of normalized free energy profiles and structural or energetic measures as a function of interpeptide separation distance. We test the approach on two families of peptides: a series of designed tetrapeptides that serve as minimal models for amyloid fibril formation, and a fragment of a classic leucine zipper peptide and its mutant. The results for these systems are compared to those from conventional REMD simulations, and demonstrate good convergence properties, low statistical errors, and, for the leucine zippers, an ability to sample near-native structures.     

A convective replica‐exchange method for sampling new energy basins

Replica‐exchange is a powerful simulation method for sampling the basins of a rugged energy landscape. The replica‐exchange method's sampling is efficient because it allows replicas to perform round trips in temperature space, thereby visiting both low and high temperatures in the same simulation. However, replicas have a diffusive walk in temperature space, and the round trip rate decreases significantly with the system size. These drawbacks make convergence of the simulation even more difficult than it already is when bigger systems are tackled. Here, we present a simple modification of the exchange method. In this method, one of the replicas steadily raises or lowers its temperature. We tested the convective replica‐exchange method on three systems of varying complexity: the alanine dipeptide in implicit solvent, the GB1 β‐hairpin in explicit solvent and the Aβ_25–35 homotrimer in a coarse grained representation. For the highly frustrated Aβ_25–35 homotrimer, the proposed “convective” replica‐exchange method is twice as fast as the standard method. It discovered 24 out of 27 free‐energy basins in less than 500 ns. It also prevented the formation of groups of replicas that usually form on either side of an exchange bottleneck, leading to a more efficient sampling of new energy basins than in the standard method. © 2012 Wiley Periodicals, Inc.     

Coulomb replica‐exchange method: Handling electrostatic attractive and repulsive forces for biomolecules

We propose a new type of the Hamiltonian replica‐exchange method (REM) for molecular dynamics (MD) and Monte Carlo simulations, which we refer to as the Coulomb REM (CREM). In this method, electrostatic charge parameters in the Coulomb interactions are exchanged among replicas while temperatures are exchanged in the usual REM. By varying the atom charges, the CREM overcomes free‐energy barriers and realizes more efficient sampling in the conformational space than the REM. Furthermore, this method requires only a smaller number of replicas because only the atom charges of solute molecules are used as exchanged parameters. We performed Coulomb replica‐exchange MD simulations of an alanine dipeptide in explicit water solvent and compared the results with those of the conventional canonical, replica exchange, and van der Waals REMs. Two force fields of AMBER parm99 and AMBER parm99SB were used. As a result, the CREM sampled all local‐minimum free‐energy states more frequently than the other methods for both force fields. Moreover, the Coulomb, van der Waals, and usual REMs were applied to a fragment of an amyloid‐β peptide (Aβ) in explicit water solvent to compare the sampling efficiency of these methods for a larger system. The CREM sampled structures of the Aβ fragment more efficiently than the other methods. We obtained β‐helix, α‐helix, 3₁₀‐helix, β‐hairpin, and β‐sheet structures as stable structures and deduced pathways of conformational transitions among these structures from a free‐energy landscape. © 2012 Wiley Periodicals, Inc.     

Temperature weighted histogram analysis method, replica exchange, and transition paths

We analyzed the data from a replica exchange molecular dynamics simulation using the weighted histogram analysis method to combine data from all of the temperature replicas (T-WHAM) to obtain the room-temperature potential of mean force of the G-peptide (the C-terminal beta-hairpin of the B1 domain of protein G) in regions of conformational space not sampled at room temperature. We were able to determine the potential of mean force in the transition region between a minor alpha-helical population and the major beta-hairpin population and identify a possible transition path between them along which the peptide retains a significant amount of secondary structure. This observation provides new insights into a possible mechanism of formation of beta-sheet secondary structures in proteins. We developed a novel Bayesian statistical uncertainty estimation method for any quantity derived from WHAM and used it to validate the calculated potential of mean force. The feasibility of estimating regions of the potential of mean force with unfavorable free energy at room temperature by T-WHAM analysis of replica exchange simulations was further tested on a system that can be solved analytically and presented some of the same challenges found in more complex chemical systems.     

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.     

On the efficiency of exchange in parallel tempering monte carlo simulations

We introduce the concept of effective fraction, defined as the expected probability that a configuration from the lowest index replica successfully reaches the highest index replica during a replica exchange Monte Carlo simulation. We then argue that the effective fraction represents an adequate measure of the quality of the sampling technique, as far as swapping is concerned. Under the hypothesis that the correlation between successive exchanges is negligible, we propose a technique for the computation of the effective fraction, a technique that relies solely on the values of the acceptance probabilities obtained at the end of the simulation. The effective fraction is then utilized for the study of the efficiency of a popular swapping scheme in the context of parallel tempering in the canonical ensemble. For large dimensional oscillators, we show that the swapping probability that minimizes the computational effort is 38.74%. By studying the parallel tempering swapping efficiency for a 13-atom Lennard-Jones cluster, we argue that the value of 38.74% remains roughly the optimal probability for most systems with continuous distributions that are likely to be encountered in practice.     

The dynamics of water exchange in gadolinium DOTA complexes studied by transition path sampling and potential of mean force methods

The mechanism of water exchange at the Gd centre of the two isomers of [Gd(iii)DOTA](-) (gadolinate(1-), [1,4,7,10-tetraazacyclododecane-1,4,7,10-tetraacetato(4-)-N1,N4,N7,N10,O1,O4,O7,O10]) has been explored using transition path sampling and potential of mean force methods to sample those regions of phase space inaccessible to standard molecular dynamics simulation. We find that there are definite differences in the details of the solvent rearrangement accompanying the exchange of the capping water molecule for the two isomers. We conclude that these solvent effects, rather than any differences in the binding energy of the capping water, are central in determining the exchange rate. We find that the potential of mean force studies yield absolute and relative rates of water exchange for the two isomers that are in good agreement with experiment.     

Adlayer morphologies and free energy landscapes of clusters of bis-fullerenes on model gold surfaces

There have been a few experimental reports of self-assembled adlayers of bis-fullerene molecules on solid substrates. Most of these studies suggest the adsorbate molecules are lying down on the surface, with the fullerene moieties almost close packed. However, very little theoretical work has been carried out on such systems, and little is known about the roles played by different parts of the potential energy in driving the self-assembly. We carry out a Temperature Replica Exchange Monte Carlo study here of two representative bis-fullerene molecules on a metal substrate. We use a coarse-grained model potential energy function, in which certain parameters can be varied within the range of their experimental uncertainty. The molecules investigated consist of two fullerene moieties bonded by a rigid bridging group. In particular, the effect of the strength of the fullerene interaction E(FG) with the substrate (nominally Au(111)) has been investigated in detail. To ensure efficient sampling of the rugged potential energy surfaces encountered in the simulations, we utilize replica exchange techniques. These enable us to construct free energy landscapes for the system. We find that for relatively low values of E(FG) the molecules form standing-up adlayers. By contrast, for higher values of E(FG), lying-down adlayers dominate. For one molecule, two different crystalline adlayer morphologies have been identified. The detailed structure of the lying-down layer is a function of the temperature and of the group used to bridge the fullerene moieties.