Molecular dynamics (MD) trajectories based on classical equations of motion can be used to sample the configurational space of complex molecular systems. However, brute-force MD often converges slowly due to the ruggedness of the underlying potential energy surface. Several schemes have been proposed to address this problem by effectively smoothing the potential energy surface. However, in order to recover the proper Boltzmann equilibrium probability distribution, these approaches must then rely on statistical reweighting techniques or generate the simulations within a Hamiltonian tempering replica-exchange scheme. The present work puts forth a novel hybrid sampling propagator combining Metropolis-Hastings Monte Carlo (MC) with proposed moves generated by non-equilibrium MD (neMD). This hybrid neMD-MC propagator comprises three elementary elements: (i) an atomic system is dynamically propagated for some period of time using standard equilibrium MD on the correct potential energy surface; (ii) the system is then propagated for a brief period of time during what is referred to as a “boosting phase,” via a time-dependent Hamiltonian that is evolved toward the perturbed potential energy surface and then back to the correct potential energy surface; (iii) the resulting configuration at the end of the neMD trajectory is then accepted or rejected according to a Metropolis criterion before returning to step 1. A symmetric two-end momentum reversal prescription is used at the end of the neMD trajectories to guarantee that the hybrid neMD-MC sampling propagator obeys microscopic detailed balance and rigorously yields the equilibrium Boltzmann distribution. The hybrid neMD-MC sampling propagator is designed and implemented to enhance the sampling by relying on the accelerated MD and solute tempering schemes. It is also combined with the adaptive biased force sampling algorithm to examine. Illustrative tests with specific biomolecular systems indicate that the method can yield a significant speedup.

Classical molecular dynamics (MD) and Metropolis-Hasting Monte Carlo (MC) simulations based on detailed atomic models are powerful tools to study the properties of complex biomolecular systems.1–3 While simulations based on realistic all-atom (AA) models arguably offer the most detailed information, such models evolve on a complex and rugged energy surface and their dynamics are often burdened by a host of slow processes. For this reason, achieving an adequate sampling of all the relevant configurations of a system from straight MD or MC simulation is often challenging.

A number of schemes have been proposed to enhance sampling and accelerate the exploration of configurational space by smoothing the underlying potential energy surface of a system.4–19 Among those schemes, two simple and attractive ideas are the hyperdynamics accelerated MD (aMD)9,10 and the replica-exchange with solute tempering (REST2) algorithms.16 However, to recover the Boltzmann equilibrium distribution of the proper Hamiltonian, one must either carry out a post-hoc reweighting analysis or generate the simulation within the context of a Hamiltonian tempering replica-exchange scheme involving multiple copies of the system.20 While the former suffers from most configurations being statistically meaningless after reweighting, the latter substantially increases the computational cost of the simulation.

One avenue to address these issues is to carry out the propagation in such a way that the resulting configurational sampling reflects the proper Hamiltonian. For example, it is possible to reformulate the self-guided Langevin dynamics (SGLD) in such a way as to restore microscopic detailed balance.19 In the present effort, we wish to explore hybrid dynamical propagation schemes that combine the strength of non-equilibrium molecular dynamics (neMD) and Metropolis Monte Carlo (MC) to achieve enhanced sampling.21–26 In hybrid neMD-MC, the value of some chosen variable or coupling parameter is altered gradually in a time-dependent and controlled fashion, while the remaining degrees of freedom are allowed to evolve freely according to the normal equations of motion. The configuration generated by such a non-equilibrium “switching” trajectory is then treated as a candidate that must be either accepted or rejected via a Metropolis criterion to generate the equilibrium Boltzmann distribution.

The hybrid neMD-MC propagator designed here comprises an “equilibrium phase,” a “boosting phase,” and a “Metropolis MC step”:

  • Equilibrium phase: An atomic system is dynamically propagated for some period of time using standard equilibrium MD on the correct potential energy surface.

  • Boosting phase: The system is then propagated for a brief period of time τ via a time-dependent Hamiltonian that is evolved toward the perturbed potential energy surface and then back to the correct potential energy surface.

  • Metropolis MC step: The resulting configuration at the end of the neMD trajectory is then accepted or rejected according to a Metropolis criterion before returning to step 1. Using a symmetric switching schedule for ramping the Hamiltonian up and down, as well as keeping detailed balance with symmetric two-end momentum reversal, ensures that the algorithm strictly produces a Boltzmann equilibrium distribution.27 

In contrast to the hybrid neMD-MC simulations guided by a coarse-grained (CG) model introduced previously,25 the algorithm described above relies on a perturbed Hamiltonian, but does not require the construction of a CG model to generate the target configuration. The hybrid neMD-MC sampling propagator implemented here rests on two schemes during the boosting phase, the so-called hyperdynamics accelerated MD (aMD)9,10 and the replica-exchange with solute tempering (REST2).16 Nevertheless, the strategy allows virtually any number of variations. Furthermore, the hybrid propagator may be naturally combined with a number of enhanced sampling strategies and free energy techniques.17 For example, preliminary results are also shown using a time-dependent bias along a collective variable determined via the adaptive biasing force (ABF) approach.28 In Sec. II, we formulate the theoretical basis of the hybrid neMD-MC sampling propagator. The performance of the method is then illustrated with specific biomolecular systems. Our results indicate that the method can yield a significant speedup for biomolecular systems.

Let us consider a classical system with time-dependent Hamiltonian Hx;λ(t), where x represents all coordinates and momenta and λ(t) is a time-dependent coupling parameter. Here, λ = 0 denotes the unperturbed Hamiltonian, H0, while λ = 1 corresponds to the perturbed Hamiltonian, H1, offering the maximum boost. In the conventional hybrid neMD-MC scheme for alchemical changes to the Hamiltonian, one varies λ from 0 to 1 in a neMD trajectory in order to generate a proposed move from x to x′. The final configuration at the end of the non-equilibrium switching period corresponds to a different Hamiltonian. For example, in constant-pH simulations based on a hybrid neMD-MC algorithm, this could be the unprotonated or protonated state of a residue.21,22 In contrast, the boosting phase of the hybrid neMD-MC propagator starts and ends with the same unperturbed Hamiltonian, H0. In other words, the system is first brought to a perturbed Hamiltonian, H1, associated with a potential energy surface that is smoother and easier to sample (λ = 0 → λ = 1) and then brought back to the unperturbed Hamiltonian H0 (λ = 1 → λ = 0). A random walk in phase space for a system at thermodynamic equilibrium must obey microscopic detailed balance,

π(x)Tp(xx)Ta(xx)=π(x)Tp(xx)Ta(xx),
(1)

where π(x) is the equilibrium probability of being in state x, Tp(xx′) is the transition probability for generating a proposed move from x to x′, and Ta(xx′) is the transition probability for accepting a proposed move from x to x′. Here, π(x) = Q−1 exp[−H(x)/kBT], where Q is the canonical partition function. After the boosting phase is completed, the energy difference ΔH=HxB,λ=0HxA,λ=0 is calculated, and a standard Metropolis acceptance criterion

Ta=min1,eβΔH
(2)

is used to enforce microscopic detailed balance (β ≡ 1/kBT, where kB is Boltzmann’s constant and T is the absolute temperature). The hybrid propagator with the different stages is depicted schematically in Fig. 1.

FIG. 1.

Hybrid neMD-MC propagator scheme comprising an equilibrium MD phase, a non-equilibrium MD boosting phase, and an acceptance or rejection via a Metropolis MC step.

FIG. 1.

Hybrid neMD-MC propagator scheme comprising an equilibrium MD phase, a non-equilibrium MD boosting phase, and an acceptance or rejection via a Metropolis MC step.

Close modal

A “symmetric two-ends momentum reversal” prescription is used in association with the non-equilibrium switches,22,25,27 by which one randomly chooses with equal probability to carry out these trajectories as a “forward” or “backward” MD propagation. Forward propagation means that the non-equilibrium switch trajectory is simply continued using current positions and momenta, whereas backward propagation means that the momenta of the system are flipped at the beginning and at the end of the non-equilibrium switch trajectory (i.e., it is equivalent to a propagation going back in time). It was shown previously that the symmetric two-ends momentum reversal prescription increases the sampling rate by reducing the likelihood that different regions of configurational space remain isolated from one another.25,27 The energy difference ΔH in Eq. (2) may be further decomposed as a sum of two terms, q and w, where q is the heat exchange between the system and an external heat bath during propagation and w is the non-equilibrium work done during the perturbation. Since it is not trivial to keep track of heat-exchange during stochastic integration (such as Langevin dynamics), we instead opt to use a deterministic integrator such that the heat exchange is formally zero. The energy difference between states A and B is equal to the non-equilibrium work applied during the trajectory augmented by the shadow work done. The shadow work corresponds to the non-equilibrium work caused by the integrator error associated with the finite time step.29 During the boosting phase, we used a schedule symmetric in time for perturbing the Hamiltonian by the coupling parameter, λ(t),

λ(t;τ,τ1)=t/τ1(0t<τ1),1(τ1t<ττ1),(τt)/τ1(ττ1t<τ),
(3)

where τ is the total switching time and τ1 is the ramping time (see Fig. 1). The middle region, during which λ is a constant equal to 1, is the boosting phase that lasts for τ − 2τ1.

In practice, the time-dependent coupling parameter λ(t) was implemented via stepwise variations with a sequence of ns discrete steps,

λ(t;τ,τ1)=0,0t<τ1/ns,1/ns,τ1/nst<2τ1/ns,(ns1)/ns,τ1(ns1)/nst<τ1,1,τ1t<ττ1,11/ns,ττ1+τ1/nst<ττ1+2τ1/ns,1/ns,τ2τ1/nst<ττ1/ns,0,ττ1/nst<τ
(4)

as depicted in Fig. 1 (see also the Tcl pseudocode in the  Appendix).

The hybrid neMD-MC propagator algorithm was implemented according to the following steps:

  • The system is dynamically propagated for some period of time using unbiased equilibrium MD with the unperturbed Hamiltonian.

  • The forward/backward direction for the switch is randomly chosen with equal probability for the neMD switching trajectory; the momenta are flipped if a backward propagation is chosen.

    • The system is propagated via the time-dependent Hamiltonian H[x, λ(t)] toward the perturbed potential energy surface (λ = 0 → λ = 1) in a time τ1.

    • The system is propagated with the perturbed Hamiltonian for some time to enhance barrier crossings (λ = 1) for a time τ − 2τ1.

    • The system is propagated via a time-dependent Hamiltonian H(t) evolved back toward the unperturbed Hamiltonian (λ = 1 → λ = 0) in a time τ1.

    • The momenta are flipped again if the switch involved a backward propagation.

  • The resulting configuration is accepted or rejected according to the Metropolis-Hasting criterion: If the move is accepted, we repeat the cycle; if the move is rejected, we return to the conformation at the end of step (i).

The performance of the hybrid neMD-MC propagator depends on our ability to generate a proposition likely to help overcome the barriers in a rugged potential energy surface. For this purpose, the parameters controlling the schedule, τ, τ1, and the choice of the perturbed Hamiltonian, H1, are critical. Two popular perturbation schemes designed to enhance sampling by deforming the potential energy surface were considered here for the boosting phase. The first one is the accelerated MD (aMD),9,10

UaMD[x,α,E]=U[x],U[x]E,U[x]+(EU[x])2α+(EU[x]),U[x]<E,
(5)

where α is a tuning parameter that determines the depth of the modified potential energy basins lying below the minimum threshold energy E (the aMD potential is flat when U[x] is equal or smaller than E if α is zero and increasingly unperturbed when α becomes very large).9,10 The aMD prescription may be applied to the total potential energy function or to various contributions such as the torsional potentials. In the present study, we have used Eq. (5) to reduce the energy barriers in the torsional potentials. To construct a time-dependent Hamiltonian for the boosting phase of the hybrid-aMD propagator, the threshold energy E was kept fixed while the tuning parameter α in Eq. (5) was effectively replaced by α/λ(t), with λ(t) following the symmetric schedule given in Eq. (3). The singularity at the first step is avoided by setting E equal to zero until λ(t) reaches its first non-zero value at time tτ1/ns following Eq. (4).

The second perturbation scheme considered here is the replica-exchange with solute tempering (REST2),16 

UREST2[x,γ]=Uv+γUu+γ1/2Uuv,
(6)

where Uv is the solvent potential energy, Uu is the solute potential energy, and Uuv is the solute-solvent interaction potential energy. It should be noted that we use the acronym REST2 here to indicate solute interaction tempering, even in the absence of multiple replicas. Let the coupling parameter γmax correspond to the maximum solute tempering allowed, traditionally expressed as T/Tu in terms of an effective temperature Tu ascribed to the solute.16 To construct a time-dependent Hamiltonian for the boosting phase of the hybrid-REST2 propagator, the coupling parameter γ in Eq. (6) was effectively replaced by the time-dependent form γ(t) = λ(t)(γmax − 1) + 1, with λ(t) following the symmetric schedule given in Eq. (3). As a result, the Hamiltonian is not perturbed when λ = 0, and the coupling parameter of the solute-tempering boost reaches its maximum of γmax when λ = 1.

The hybrid neMD-MC propagator was tested on several biologically relevant molecular systems. The first objective of these tests is to show that the propagator does produce the correct Boltzmann equilibrium distribution. The second objective of these tests is to show that the propagator achieves convergence of the equilibrium properties more rapidly than simple unbiased MD. The hybrid neMD-MC propagator was implemented as a general Tcl script for the simulation program NAMD.30 The Tcl script is given in the  Appendix.

The first test case is deca-alanine in vacuum with a dielectric constant ϵ = 20. Eight individual 200-ns trajectories were generated using (i) equilibrium MD, (ii) two hybrid propagators using either REST2 or aMD as a boost, (iii) two hybrid propagators randomly accepting proposed moves by means of mean-acceptance ratio from the corresponding propagators, and (iv) two accelerating schemes, namely, REST2 and aMD. The equations of motion were integrated with a time step of 1 fs. The temperature was maintained at 300 K using Langevin dynamics with a damping coefficient of 1.0 ps−1. Nonbonded short range interactions were truncated at 9 Å with a switching function effective from 8 Å. The parameters of the switching schedule were set to teqMD = 5 ps, τ = 5 ps, τ1 = 2 ps, number of switches ns = 39, maximum REST2 boost γ = 0.8, and maximum aMD boost (E, α) = (50, 10). The CHARMM 22 protein force field31 was employed to model deca-alanine, allowing a direct comparison with a previous theoretical study.28 A second test case considers again deca-alanine in vacuum, but this time with a dielectric constant, ϵ, of 1. This system is normally difficult to sample using simple unbiased equilibrium MD trajectories because of the strong internal backbone-backbone hydrogen bonds, making the usage of some form of importance sampling strategy a necessity. Here, the adaptive biasing force (ABF) method was considered.32,33 Results from ABF with standard MD and from ABF with hybrid neMD-MC boosted by aMD and REST2 were examined. Eight individual 40-ns trajectories were generated for each case, using α-helix and C5-extended conformation as initial coordinates. After optimization, the switching schedule was set to teqMD = 7 ps, τ = 2 ps, τ1 = 1 ps, number of switches ns = 18, maximum REST2 boost γ = 0.5, and maximum aMD boost (E, α) = (40, 1). Extensive computations based on the multiple-walker ABF method32,33 were performed to provide a reference potential of mean force (PMF) that can be used to assess the converge of the different methods. The computation was carried out with eight independent walkers, corresponding to an aggregated simulation time of 8 μs.

A third test case considers the folding/unfolding of the acetyl and NH2 terminally capped (AAQAA)3 peptide solvated in explicit water. The 15 residue peptide was simulated using the hybrid-REST2 neMD-MC propagator. As a basis of comparison, the system was also simulated with simple unbiased MD, as well as multiple-copy REST2 and temperature replica-exchange MD (TREMD) simulations.

For the optimization, different switching schedules were tested while varying τ, τ1, the maximum boost (REST2 or aMD), and number of switches. After optimization, the parameters of the switching schedule were set to teqMD = 2 ps, τ = 24 ps, τ1 = 8 ps, number of switches ns = 790, and maximum REST2 boost γ = 0.7. The system was equilibrated at 300 K with initial structure of complete α-helical conformation solvated with 3964 water molecules in a cubic cell of initial dimensions equal to 55 × 55 × 55 Å3. The equations of motion were integrated with a time step of 2 fs. The temperature was maintained at 300 K using Langevin dynamics with a damping coefficient of 1 ps−1. Nonbonded short-range interactions were truncated at 12 Å with a switching function effective from 10 Å. The CHARMM 36 force field was used to describe the peptide, with the TIP3P water model. The REST2 calculation was done with 10 replicas ranging up to the same maximum boost of γ = 0.7 as in the hybrid neMD-MC simulations and exchange attempts every 0.2 ps. The TREMD simulation was done with 32 replicas ranging from 278 K to 416 K and exchange attempts every 10 ps, matching exactly the protocol of Best et al.34 However, the total length of the present TREMD simulation is much shorter (1.875 ns per replica and 60 ns in total) than that of Best et al.34 (150 ns per replica and 4.8 μs in total). To provide an objective comparison of the actual computational cost of these single- and multiple-copy approaches, the results are reported in terms of MD steps per replica times the number of replica involved.

We first examine the results for deca-alanine in vacuum with a dielectric constant ϵ = 20. The main purpose of this test is to validate the hybrid neMD-MC propagator. Because of the high dielectric constant, this is a convenient toy model, easy to sample directly by brute-force MD. The PMF as a function of the end-to-end distance extracted from unbiased MD is shown in Fig. 2. It is a double-well potential, with basins at 6 Å and 12 Å, separated by a small barrier of about 0.5 kcal/mol. The PMF calculated using the hybrid propagator is also shown in Fig. 2. The comparison shows that the hybrid method correctly reproduces the PMF extracted from the unbiased MD trajectory. The small error bars corresponding to a 99% confidence interval are indicative of the high convergence. As expected, the PMFs generated by pure REST2 and aMD without any reweighting do not reproduce the correct result. The biased PMFs nonetheless serve to display the overall free energy surface when the system is propagated at the highest boosting level using these two perturbation schemes designed to enhance sampling by deforming the potential energy surface. Finally, the random acceptance simulations with the same mean acceptance ratio from the hybrid-aMD (0.1522 ± 0.002) and hybrid-REST2 simulations (0.716 ± 0.003) deviate from the correct PMF. This shows that the acceptance criterion with the symmetric two-ends momentum reversal is critical to obey detailed balance and obtain correct results.

FIG. 2.

Validation test using deca-alanine with dielectric constant of 20. SA-Random stands for hybrid scheme that randomly accepts proposed moves using the same mean acceptance ratio achieved from hybrid simulations.

FIG. 2.

Validation test using deca-alanine with dielectric constant of 20. SA-Random stands for hybrid scheme that randomly accepts proposed moves using the same mean acceptance ratio achieved from hybrid simulations.

Close modal

Generally, the efficiency of neMD-MC depends on maximizing the production of uncorrelated configurations while trying to minimize the effort spent to generate these.26 The main drawback of neMD-MC simulations is that each new attempted MC move requires a non-equilibrium MD simulation; a low acceptance rate necessarily implies that a large fraction of computer time is discarded by the algorithm. Whereas long switches might yield a high acceptance probability but are computationally prohibitive, short switches are computationally inexpensive, but expected to yield vanishingly low acceptance probabilities. The most efficient algorithm is obviously a compromise balancing between these two opposing factors. It is generally unclear as to how to systematically achieve this balance.

In conventional neMD-MC schemes, the switch is meant to allow transitions between two different discrete states, thus longer switch times lead to higher acceptance rates as the transformation approaches the adiabatic limit.26 In a previous study of a constant-pH algorithm, an optimal switching time for changing the protonation state of propionic acid in an explicit solvent was determined to be about 15 ps.26 However, the hybrid propagator described here is not designed to produce transitions between discrete states. One implication is that a high acceptance probability, i.e., the number of accepted moves over the total number of attempted moves, does not guarantee that one is using the most effective perturbation scheme and boosting schedule. For example, the accepted candidate configuration may be too close to the starting configuration if the perturbation is too mild. Within the constraint of a fixed budget of computer time, increasing indefinitely the switching time also implies a concomitant decrease in the time available to perform equilibrium MD. Under these premises, simply maximizing the acceptance probability is not a useful criterion for optimizing efficiency.

Several factors could affect the efficiency of the hybrid neMD-MC scheme described here. As a simple measure of efficiency of the hybrid neMD-MC propagation, we define a “conformational evolution speed” given by the net root-mean-square deviation (RMSD) of the peptide at a time t relative to the peptide configuration at time t − ΔT. Monitoring the conformational evolution speed, while including the total cost to generate the non-equilibrium switches, allows us to assess the net efficiency of the hybrid neMD-MC propagator in terms of the perturbation scheme and boosting schedule. As illustrated in Fig. 3, it is indeed possible to identify a maximum efficiency in the case of the deca-alanine system. Based on this metric, a switching time of about 1.8 ps yields the most efficient hybrid neMD-MC propagator for the deca-alanine system.

FIG. 3.

Conformational evolution speed for deca-alanine in vacuum with a dielectric constant of 5 as a function of the boosting schedule. The conformational evolution speed is root-mean-square deviation (RMSD) of the peptide (calculated from all heavy atoms) at a time t relative to its configuration at a previous time t −ΔT. For each case, the time interval ΔT was set to its switching time. The blue line is a fitted function for visual guide.

FIG. 3.

Conformational evolution speed for deca-alanine in vacuum with a dielectric constant of 5 as a function of the boosting schedule. The conformational evolution speed is root-mean-square deviation (RMSD) of the peptide (calculated from all heavy atoms) at a time t relative to its configuration at a previous time t −ΔT. For each case, the time interval ΔT was set to its switching time. The blue line is a fitted function for visual guide.

Close modal

For the (AAQAA)3 system in an explicit solvent, we monitored the fraction of secondary-structure elements in the peptide chain as a function of switching time. As illustrated in Fig. 4, the efficiency is fairly poor when the switching time is too short. The simulation started with the peptide in an α-helical conformation, which persisted for an extended period of time. In contrast, the efficiency is considerably improved by employing a switching time of 24 ps. While it is possible that this could lead to further improvement, the behavior of the system was not examined for switching times longer than 24 ps. In practice, the efficiency of the (AAQAA)3 system in an explicit solvent is already very good with 24 ps according to Fig. 4. As shown in Table I, the acceptance probability Pa is about 0.92 with a switching time of 24 ps. This is approaching the limit afforded by the shadow work, which would be present with normal equilibrium MD propagation.

FIG. 4.

Optimization test for the (AAQAA)3 system in an explicit solvent by monitoring time cumulative average for the α-helicity. 4 independent trajectories were simulated for all cases, and 4 additional independent simulations were operated for MD and 24 ps switch. The peptide is defined to be in the α-helical conformation when the (ϕ, ψ) backbone dihedral angles of three consecutive residues fulfill the conditions |−65° − ϕ| < 35° and |−37° − ψ| < 30°.

FIG. 4.

Optimization test for the (AAQAA)3 system in an explicit solvent by monitoring time cumulative average for the α-helicity. 4 independent trajectories were simulated for all cases, and 4 additional independent simulations were operated for MD and 24 ps switch. The peptide is defined to be in the α-helical conformation when the (ϕ, ψ) backbone dihedral angles of three consecutive residues fulfill the conditions |−65° − ϕ| < 35° and |−37° − ψ| < 30°.

Close modal
TABLE I.

Acceptance probability for the (AAQAA)3 peptide.a

τ (ps)τ1 (ps)Pa
10 0.42 
12 0.60 
14 0.72 
16 0.79 
24 0.92 
τ (ps)τ1 (ps)Pa
10 0.42 
12 0.60 
14 0.72 
16 0.79 
24 0.92 
a

τ is the total switching time, τ1 is the ramping time, and the boosting phase lasts τ − 2τ1; 4 neMD-MC simulations were carried out for a switching time of 10, 12, 14, and 16 ps, and 8 neMD-MC simulations were carried out for a switching time of 24 ps.

Deca-alanine is an extremely well-characterized benchmark system for testing the ability of new sampling schemes to capture subtle conformational equilibria.28,35 Its global free-energy minimum corresponds to a fully α-helical conformation resulting from the formation of robust ii + 4 intra-molecular hydrogen bonds. When generating a PMF along its end-to-end distance, the metastable states that exist along the orthogonal degrees of freedom in the “rugged” region between 4 and 12 Å can dramatically slow the convergence of conformational sampling. Furthermore, there is a very large free energy difference of about 30 kcal/mol between the completely folded form (14 Å) and the family of unfolded structures (32 Å). As a result, biased simulations must be used to overcome the considerable imbalance in the equilibrium probability along the end-to-end distance. Here, we use the adaptive biasing force (ABF) method.32,33 It should be emphasized that the necessity to adopt some form of importance sampling method to explore the end-to-end distance of the deca-alanine system remains true, whether a hybrid neMD-MC propagator is used or not.

The main results are displayed in Fig. 5. Comparison of the various methods reveals that the hybrid propagator with ABF clearly converged faster toward the reference PMF. While the standard deviations of the PMF based on eight independent simulations, either MD or the hybrid, are comparable in regions of the reaction pathway that are easily sampled (>12 Å), those for the hybrid method are certainly lower in the rugged region (<12 Å). This is indicative that the hybrid method facilitates the sampling of the numerous metastable states lying along orthogonal degrees of freedom in this region. Interestingly, when simulating multiple independent trajectories with pure ABF, we observed that some trajectories occasionally remained “trapped” in metastable regions for an extended fraction of the simulation time. This occurred for ABF with both starting conformations, leading to an over-stabilization of the metastable states and resulting in slow convergence. If one consider these trapped simulations as outliers and discard them, then ABF-MD performs overall better than ABF-hybrid. Importantly, this highly undesirable trapping phenomenon was not observed for the ABF-hybrid simulations. In assessing the overall efficiency of the ABF-hybrid simulations, it is important to note that no data are collected about the gradient during the neMD switching trajectories. Thus, the total amount of data accumulated to estimate the mean gradient is necessarily smaller than with ABF-MD. The hybrid scheme with ABF having only 70% of equivalent MD data points, yet performing better, underscores its power.

FIG. 5.

PMF of the end-to-end distance of tests on deca-alanine in vacuum computed using ABF and ABF with hybrid-(REST2, aMD). Results using α-helix (top) and C5-extended (bottom) conformation as the starting structure are shown. 8 independent trajectories of 40 ns were simulated using ABF and ABF-hybrid.

FIG. 5.

PMF of the end-to-end distance of tests on deca-alanine in vacuum computed using ABF and ABF with hybrid-(REST2, aMD). Results using α-helix (top) and C5-extended (bottom) conformation as the starting structure are shown. 8 independent trajectories of 40 ns were simulated using ABF and ABF-hybrid.

Close modal

As a final illustrative test, we consider the folding/unfolding of the 15 residue peptide (AAQAA)3 in an explicit solvent. This system has been the object of extensive simulation studies and is very well characterized.25,34 To highlight the conformational sampling challenge, the simulations were all started in the long-lived α-helix metastable state. Let us first compare the results from simple brute-force equilibrium MD and a hybrid neMD-MC propagator with REST2 boosting. Figure 6 shows the efficiency of the hybrid-REST2 propagator through time-cumulative averages of secondary structures (α-helix and β-strand). It is observed that the simulation converges towards the reference equilibrium average much faster with the hybrid propagator. This trend is systematic, as indicated by the 99% confidence interval estimated from the 8 independent trajectories. Rapid conformational fluctuations of the (AAQAA)3 peptide are visibly occurring with the hybrid-REST2 propagator, as shown by the time-series of the individual simulations (Figs. S1–S4 of the supplementary material). Conversely, the systems simulated via simple brute-force equilibrium MD remain essentially trapped in the starting α-helical conformation. While the brute-force equilibrium MD generated 13 times more data points than the hybrid propagator, the conformational space of the (AAQAA)3 peptide was explored much more efficiently with the hybrid propagator. Remarkably, REST2 (10 replicas) and TREMD (32 replicas) result in a weaker performance than the single-copy hybrid-REST2 neMD-MC propagator. The conformational sampling of (AAQAA)3 in solution remains inefficient with these multiple-copy approaches, despite the fact that there were frequent exchanges between the replicas (Figs. S5 and S6 of the supplementary material). The difference in performance is particularly glaring when comparing in terms of total computational cost of all approaches (Fig. 6). It is of interest to compare these results with the ABF calculations for the deca-alanine toy model discussed above. In the case of the ABF-hybrid simulations, nearly 70% of the data produced correspond to equilibrium sampling and were accumulated to estimate the mean gradient. Increasing the amount of non-equilibrium switches is not advantageous because it reduces the amount of data used to estimate the mean gradient needed in ABF. This points to the limitation associated with the inherent loss of data during the non-equilibrium switches in the context of ABF. When the data collected from simple brute-force equilibrium MD are sufficiently uncorrelated, as in the case of deca-alanine, then increasing the proportion of the sampling carried out via the hybrid propagator becomes statistically less efficient. Even if the hybrid propagator actually moves the system around conformational space rapidly, enough data points must still be accumulated to accurately determine the average gradient that is needed for the ABF method. This stands in contrast with the (AAQAA)3 simulation, for which the fraction of helicity could be determined accurately even though only 7.7% of the data produced corresponds to equilibrium sampling. The hybrid propagator can efficiently explore the conformational space accessible to the peptide in solution and the result without relying on accumulated equilibrium data as in the case of ABF. Even though the brute-force equilibrium MD simulations include about 13 times more data points than the hybrid simulation, it is extremely inefficient to explore the conformational space accessible to the (AAQAA)3 peptide in solution.

FIG. 6.

Cumulative time-average for the secondary structure elements of a (AAQAA)3 peptide in an explicit solvent using unbiased MD (single-copy), hybrid-REST2 neMD-MC (single-copy), REST2 with 10 replicas, and TREMD with 32 replicas. 8 independent trajectories of 60 ns were simulated using unbiased MD and the hybrid neMD-MC propagator, 6 ns each of 10 replicas were simulated using REST2, and 1.875 ns each of 32 replicas were simulated using TREMD. The total number of MD steps on the x-axis represents the actual computational costs of each approach as the number of steps per replica times the number of replica. The peptide is defined to be in the α-helical conformation when the (ϕ, ψ) backbone dihedral angles of three consecutive residues fulfill the conditions |−65° − ϕ| < 35° and |−37° − ψ| < 30°. The peptide is in the β-strand conformation when (ϕ, ψ) fulfill the conditions |−140° − ϕ| < 40° and |150° − ψ| < 30°.

FIG. 6.

Cumulative time-average for the secondary structure elements of a (AAQAA)3 peptide in an explicit solvent using unbiased MD (single-copy), hybrid-REST2 neMD-MC (single-copy), REST2 with 10 replicas, and TREMD with 32 replicas. 8 independent trajectories of 60 ns were simulated using unbiased MD and the hybrid neMD-MC propagator, 6 ns each of 10 replicas were simulated using REST2, and 1.875 ns each of 32 replicas were simulated using TREMD. The total number of MD steps on the x-axis represents the actual computational costs of each approach as the number of steps per replica times the number of replica. The peptide is defined to be in the α-helical conformation when the (ϕ, ψ) backbone dihedral angles of three consecutive residues fulfill the conditions |−65° − ϕ| < 35° and |−37° − ψ| < 30°. The peptide is in the β-strand conformation when (ϕ, ψ) fulfill the conditions |−140° − ϕ| < 40° and |150° − ψ| < 30°.

Close modal

Powerful hybrid neMD-MC algorithms can be designed to boost and accelerate the conformational sampling of complex molecular systems. At the heart of the hybrid neMD-MC algorithms is a transient non-equilibrium perturbation of the system, which is aimed at overcoming the barriers in the rugged energy landscapes. Importantly, the neMD-MC algorithms robustly generate the correct Boltzmann equilibrium distribution. A vast range of time-dependent perturbations are possible to construct a family of hybrid neMD-MC propagators adapted to different situations. In the present test, we boosted the conformational sampling by relying on Hamiltonian perturbations based on aMD and REST2, but the strategy allows virtually any accelerated methods to be introduced during the neMD steps to efficiently lower the free-energy barriers. The benchmark systems tested demonstrate the correctness and effectiveness of the hybrid propagator, emphasizing its faster convergence compared to equilibrium MD. Although the switching schedule was optimized heuristically in the present study, it is possible to infer it analytically when the Hamiltonian is perturbed linearly.26 Depending on the problem at hand, hybrid neMD-MC propagators may advantageously be combined with a number of established strategies and free energy methods,17 including umbrella sampling (US),36 adaptive biasing force (ABF),28 Hamiltonian tempering replica-exchange,37 and alchemical free energy perturbation.12 Future work will explore how the present hybrid neMD-MC propagators behave in the case of very complex biological objects and can be possibly tailored to address the rugged free-energy landscape underlying intricate processes therein.

See supplementary material for details (time-series, time-cumulative average, and accepted exchanges) about the simulations of the (AAQAA)3 peptide in an explicit solvent (MD, hybrid neMD-MC, REST2, and TREMD).

This work was supported by the National Science Foundation (NSF) through Grant No. MCB-1517221. Additional support was provided by the France and Chicago Collaborating in The Sciences (FACCTS) program (to B.R. and C.C.). An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program; this research used resources of the Argonne Leadership Computing Facility (ALCF), which is a DOE Office of Science User Facility supported under Contract No. DE-AC02-06CH11357. Additional computational resources were provided by the University of Chicago Research Computing Center and the Grand Equipement National de Calcul Intensif-Centre Informatique National de lEnseignement Superieur (GENCI-CINES) at Montpelier, France (to C.C.).

Schematic Tcl pseudocode for executing the hybrid neMD-MC propagator with the NAMD30 program:

1.
M.
Karplus
and
J. A.
McCammon
, “
Molecular dynamics simulations of biomolecules
,”
Nat. Struct. Mol. Biol.
9
,
646
652
(
2002
).
2.
N.
Metropolis
,
A. W.
Rosenbluth
,
M. N.
Rosenbluth
, and
A. H.
Teller
, “
Equation of state calculations by fast computing machines
,”
J. Chem. Phys.
21
,
1087
1092
(
1953
).
3.
W. K.
Hastings
, “
Monte Carlo sampling methods using Markov chains and their applications
,”
Biometrika
57
,
97
109
(
1970
).
4.
T.
Huber
,
A. E.
Torda
, and
W. F.
van Gunsteren
, “
Local elevation: A method for improving the searching properties of molecular dynamics simulation
,”
J. Comput.-Aided Mol. Des.
8
,
695
708
(
1994
).
5.
V.
Leone
,
F.
Marinelli
,
P.
Carloni
, and
M.
Parrinello
, “
Targeting biomolecular flexibility with metadynamics
,”
Curr. Opin. Struct. Biol.
20
,
148
154
(
2010
).
6.
X.
Wu
and
B. R.
Brooks
, “
Self-guided Langevin dynamics simulation method
,”
Chem. Phys. Lett.
381
,
512
518
(
2003
).
7.
A.
Damjanović
,
X.
Wu
,
B. R.
Brooks
 et al, “
Backbone relaxation coupled to the ionization of internal groups in proteins: A self-guided Langevin dynamics study
,”
Biophys. J.
95
,
4091
4101
(
2008
).
8.
X.
Wu
and
B. R.
Brooks
, “
Toward canonical ensemble distribution from self-guided Langevin dynamics simulation
,”
J. Chem. Phys.
134
(
13
),
134108
(
2011
).
9.
A. F.
Voter
, “
Hyperdynamics: Accelerated molecular dynamics of infrequent events
,”
Phys. Rev. Lett.
78
,
3908
(
1997
).
10.
D.
Hamelberg
,
J.
Mongan
, and
J. A.
McCammon
, “
Accelerated molecular dynamics: A promising and efficient simulation method for biomolecules
,”
J. Chem. Phys.
120
,
11919
(
2004
).
11.
M.
Fajer
,
D.
Hamelberg
, and
J. A.
McCammon
, “
Replica-exchange accelerated molecular dynamics (REXAMD) applied to thermodynamic integration
,”
J. Chem. Theory Comput.
4
,
1565
1569
(
2008
).
12.
W.
Jiang
and
B.
Roux
, “
Free energy perturbation Hamiltonian replica-exchange molecular dynamics (FEP/H-REMD) for absolute ligand binding free energy calculations
,”
J. Chem. Theory Comput.
6
,
2559
2565
(
2010
).
13.
L.
Maragliano
and
E.
Vanden-Eijnden
, “
A temperature accelerated method for sampling free energy and determining reaction pathways in rare events simulations
,”
Chem. Phys. Lett.
426
,
168
175
(
2006
).
14.
H.
Lei
and
Y.
Duan
, “
Improved sampling methods for molecular simulation
,”
Curr. Opin. Struct. Biol.
17
,
187
191
(
2007
).
15.
M.
Christen
and
W. F.
van Gunsteren
, “
On searching in, sampling of, and dynamically moving through conformational space of biomolecular systems: A review
,”
J. Comput. Chem.
29
,
157
166
(
2008
).
16.
L.
Wang
,
R.
Friesner
, and
B. J.
Berne
, “
Replica exchange with solute scaling: A more efficient version of replica exchange with solute tempering (REST2)
,”
J. Phys. Chem. B
115
,
9431
9438
(
2011
).
17.
A.
Mitsutake
,
Y.
Mori
, and
Y.
Okamoto
, “
Enhanced sampling algorithms
,”
Methods Mol. Biol.
924
,
153
195
(
2013
).
18.
R. C.
Bernardi
,
M. C.
Melo
, and
K.
Schulten
, “
Enhanced sampling techniques in molecular dynamics simulations of biological systems
,”
Biochim. Biophys. Acta
1850
,
872
877
(
2015
).
19.
X.
Wu
,
B. R.
Brooks
, and
E.
Vanden-Eijnden
, “
Self-guided Langevin dynamics via generalized Langevin equation
,”
J. Comput. Chem.
37
,
595
601
(
2016
).
20.
M.
Arrar
,
C. A. F.
de Oliveira
,
M.
Fajer
,
W.
Sinko
, and
J. A.
McCammon
, “
w-REXAMD: A Hamiltonian replica exchange approach to improve free energy calculations for systems with kinetically trapped conformations
,”
J. Chem. Theory Comput.
9
,
18
23
(
2013
).
21.
H. A.
Stern
, “
Molecular simulation with variable protonation states at constant pH
,”
J. Chem. Phys.
126
,
164112
(
2007
).
22.
H. A.
Stern
, “
Erratum: Molecular simulation with variable protonation states at constant pH
,”
J. Chem. Phys.
127
,
079901
(
2007
).
23.
A. J.
Ballard
and
C.
Jarzynski
, “
Replica exchange with nonequilibrium switches
,”
Proc. Natl. Acad. Sci. U. S. A.
106
,
12224
12229
(
2009
).
24.
J. P.
Nilmeier
,
G. E.
Crooks
,
D. D. L.
Minh
, and
J. D.
Chodera
, “
Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation
,”
Proc. Natl. Acad. Sci. U. S. A.
108
,
E1009
E1018
(
2011
).
25.
Y.
Chen
and
B.
Roux
, “
Enhanced sampling of an atomic model with hybrid nonequilibrium molecular dynamics–Monte Carlo simulations guided by a coarse-grained model
,”
J. Chem. Theory Comput.
11
,
3572
3583
(
2015
).
26.
B. K.
Radak
and
B.
Roux
, “
Efficiency in nonequilibrium molecular dynamics Monte Carlo simulations
,”
J. Chem. Phys.
145
,
134109
(
2016
).
27.
Y.
Chen
and
B.
Roux
, “
Efficient hybrid non-equilibrium molecular dynamics–Monte Carlo simulations with symmetric momentum reversal
,”
J. Chem. Phys.
141
,
114107
(
2014
).
28.
J.
Comer
,
J. C.
Gumbart
,
J.
Hénin
,
T.
Lelièvre
,
A.
Pohorille
, and
C.
Chipot
, “
The adaptive biasing force method: Everything you always wanted to know but were afraid to ask
,”
J. Phys. Chem. B
119
,
1129
1151
(
2015
).
29.
D. A.
Sivak
,
J. D.
Chodera
, and
G. E.
Crooks
, “
Using nonequilibrium fluctuation theorems to understand and correct errors in equilibrium and nonequilibrium simulations of discrete Langevin dynamics
,”
Phys. Rev. X
3
,
011007
(
2013
).
30.
J. C.
Phillips
,
R.
Braun
,
W.
Wang
,
J.
Gumbart
,
E.
Tajkhorshid
,
E.
Villa
,
C.
Chipot
,
R. D.
Skeel
,
L.
Kalé
, and
K.
Schulten
, “
Scalable molecular dynamics with NAMD
,”
J. Comput. Chem.
26
,
1781
1802
(
2005
).
31.
A. D.
MacKerell
, Jr.
,
D.
Bashford
,
M.
Bellott
,
R. L.
Dunbrack
, Jr.
,
J. D.
Evanseck
,
M. J.
Field
,
S.
Fischer
,
J.
Gao
,
H.
Guo
,
S.
Ha
,
D.
Joseph-McCarthy
,
L.
Kuchnir
,
K.
Kuczera
,
F. T. K.
Lau
,
C.
Mattos
,
S.
Michnick
,
T.
Ngo
,
D. T.
Nguyen
,
B.
Prodhom
,
W. E.
Reiher
 III
,
B.
Roux
,
M.
Schlenkrich
,
J. C.
Smith
,
R.
Stote
,
J.
Straub
,
M.
Watanabe
,
J.
Wiórkiewicz-Kuczera
,
D.
Yin
, and
M.
Karplus
, “
All-atom empirical potential for molecular modeling and dynamics studies of proteins
,”
J. Phys. Chem. B
102
,
3586
3616
(
1998
).
32.
J.
Comer
,
J.
Phillips
,
K.
Schulten
, and
C.
Chipot
, “
Multiple-replica strategies for free-energy calculations in NAMD: Multiple-walker adaptive biasing force and walker selection rules
,”
J. Chem. Theory Comput.
10
,
5276
5285
(
2014
).
33.
K.
Minoukadeh
,
T.
Lelièvre
, and
C.
Chipot
, “
Potential of mean force calculations: A multiple-walker adaptive biasing force approach
,”
J. Chem. Theory Comput.
6
,
1008
1017
(
2010
).
34.
R. B.
Best
,
X.
Zhu
,
J.
Shim
,
P. E. M.
Lopes
,
J.
Mittal
,
M.
Feig
, and
A. D.
MacKerell
, Jr.
, “
Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone ϕ, ψ and side-chain χ1 and χ2 dihedral angles
,”
J. Chem. Theory Comput.
8
,
3257
3273
(
2012
).
35.
S.
Park
,
F.
Khalili-Araghi
,
E.
Tajkhorshid
, and
K.
Schulten
, “
Free energy calculation from steered molecular dynamics simulations using Jarzynski’s equality
,”
J. Chem. Phys.
119
,
3559
3566
(
2003
).
36.
G. M.
Torrie
and
J. P.
Valleau
, “
Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling
,”
J. Comput. Phys.
23
,
187
199
(
1977
).
37.
W.
Jiang
,
Y.
Luo
,
L.
Maragliano
, and
B.
Roux
, “
Calculation of free energy landscape in multi-dimensions with Hamiltonian-exchange umbrella sampling on petascale supercomputer
,”
J. Chem. Theory Comput.
8
,
4672
4680
(
2012
).

Supplementary Material