Multiscale implementation of infinite-swap replica exchange molecular dynamics

Abstract

Replica exchange molecular dynamics (REMD) is a popular method to accelerate conformational sampling of complex molecular systems. The idea is to run several replicas of the system in parallel at different temperatures that are swapped periodically. These swaps are typically attempted every few MD steps and accepted or rejected according to a Metropolis-Hastings criterion. This guarantees that the joint distribution of the composite system of replicas is the normalized sum of the symmetrized product of the canonical distributions of these replicas at the different temperatures. Here we propose a different implementation of REMD in which (i) the swaps obey a continuous-time Markov jump process implemented via Gillespie's stochastic simulation algorithm (SSA), which also samples exactly the aforementioned joint distribution and has the advantage of being rejection free, and (ii) this REMD-SSA is combined with the heterogeneous multiscale method to accelerate the rate of the swaps and reach the so-called infinite-swap limit that is known to optimize sampling efficiency. The method is easy to implement and can be trivially parallelized. Here we illustrate its accuracy and efficiency on the examples of alanine dipeptide in vacuum and C-terminal β-hairpin of protein G in explicit solvent. In this latter example, our results indicate that the landscape of the protein is a triple funnel with two folded structures and one misfolded structure that are stabilized by H-bonds.

Keywords: HMM; REMD; SSA; importance sampling; protein-folding.

Fig. S1.

Reconstructed energy profiles [$-k_{B}T\log \rho \left(x\right)$] at $k_{B}T=0.08$ for the 1D model. From left to right: ISREMD, REMD-SSA at $\nu =\mathrm{150,000}/\mathrm{\Delta }t,\mathrm{15,000}/\mathrm{\Delta }t,\mathrm{1,500}/\mathrm{\Delta }t,150/\mathrm{\Delta }t$, REMD-SSA using 40 replicas with $\nu =\mathrm{1,500}$. The average number of jumps for these ν are $\sim {10}^3$, ${10}^4$, ${10}^5$, and ${10}^6$. The blue line is the analytical results and the orange circle is the numerical results. All results are from ${10}^5$ MD steps.

Fig. S2.

Running estimates for $p_L$ in the 1D example from REMD-SSA (orange), and standard REMD at rate 0.05 per step (yellow), 0.01 per step (purple), and 0.001 per step (green). The blue curve shows the analytical value of $p_L$. The insets show the energy profiles reconstructed by the various methods after ${10}^5$ steps. For REMD-SSA, $\nu =\mathrm{1,500}$ per step.

Fig. S3.

First row: reconstructed energy profiles [$-k_{B}T\log P\left(X\right)$] at $k_{B}T=0.04$ for 1D toy model using the temperature set $\left\{0.04,0.10,0.22,0.53,1.25\right\}$. From left to right: ISREMD and REMD-SSA at $\nu =\mathrm{15,000}/\mathrm{\Delta }t,\mathrm{1,500}/\mathrm{\Delta }t,150/\mathrm{\Delta }t$. The average number of jumps for these ν are $\sim {10}^5$, ${10}^4$, and ${10}^3$. The simulation is of ${10}^5$ MD steps. Second row: reconstructed energy profiles from standard REMD at $k_{B}T=0.08$ for 1D toy model using the tempearature set $\left\{0.08,0.16,0.32,0.63,1.25\right\}$. Swap rate from left to right: $0.02/\mathrm{\Delta }t$, $0.05/\mathrm{\Delta }t$, $0.1/\mathrm{\Delta }t$, and $0.5/\mathrm{\Delta }t$. Simulation length is ${10}^6$ steps. The blue line is the analytical results and the orange circle is the numerical results.

Fig. 1.

(A) Free energy surfaces along the two dihedral angles ϕ and ψ for AD in vacuum at 300 K. (Upper Left) Fifty-nanosecond REMD-MSSA simulation with four replicas. (Upper Right) Fifty-nanosecond ISREMD simulation with four replicas using the analytical weights calculated from Eq. 10. (Lower Left) Fifty-nanosecond REMD-MSSA simulation with 12 replicas. (Lower Right) Fifty-nanosecond standard REMD simulation with 12 replicas and a swapping rate of 0.5 ps⁻¹. All plots have $60\times 60$ bins and filtered by standard Gaussian kernel. Same level sets are used in the contour plots. (B) Trajectories $\sigma \left(i,t\right)$ of one replica. (C) Distributions of $\sigma \left(i,t\right)$ for the same replica after 0.8 ns of simulation. In B and C, the orange curve is from REMD-MSSA and the blue one is from standard REMD, both with 12 replicas.

Fig. S4.

AD: Temperature trajectories from REMD-SSA (orange) and standard REMD (blue) of all replicas in the case of using 12 temperatures. X, time ($\times$ 200 ps); Y, temperature index.

Fig. S5.

AD: Temperature distributions from REMD-SSA (orange) and standard REMD (blue) of all replicas in the case of 12 temperatures. X, temperature index; Y, population fraction. The results are from the 0.8-ns run.

Fig. 2.

Free energy surfaces of β-hairpin at 270 K along the two order parameters β-strand HB number and backbone radius of gyration RG (A) from 120 ns of standard REMD at swapping rate 1 ${\text{ps}}^{-1}$ and (B) from 100 ns of REMD-MSSA. (C) Trajectories of $\sigma \left(i,t\right)$ for a given replica i. (D) Distributions of $\sigma \left(i,t\right)$ for this replica, from a 36-ns-long run. (E) Relative entropy of the distributions of $\sigma \left(i,t\right)$ (in an ascending order). In C–E the orange curve is for REMD-MSSA and the blue one for standard REMD.

Fig. S6.

β-hairpin: Traces of β-strand H-bond number ${\text{N}}_{\text{H}}$ (black) and C_α rmsd (red) relative to the native conformation from representative replicas in the REMD-SSA simulation. β-sheet folded states are highlighted in cyan.

Fig. S7.

Free energy surface along β-strand number ($N_H$) and all HB number on backbone (all $N_H$) at 270 K constructed from 100-ns REMD-SSA. Shown around the FES are the typical conformations in basin a, b, and c with key hydrogen bonds and salt bridges on display.

Fig. S8.

β-hairpin: Temperature trajectories from REMD-SSA (orange) and standard REMD (blue) of 12 representative replicas. X: time (ns); Y: temperature index.

Fig. S9.

β-hairpin: Temperature distributions from REMD-SSA (orange) and standard REMD (blue) for the 12 replicas. X, temperature index; Y, population fraction. The results are from the 36-ns run.

Fig. 3.

(A) Native β-hairpin structure with H-bond registry indicated. (B) H–O contact map of for ${\beta }_N$, ${\beta }_2$, and M; in each map, the horizontal axis indicates the residue number of the amide H and the vertical axis the residue number of the carbonyl O. (C) FES at 270 K in the space of number of β-strand H-bonds ${\text{N}}_{\text{H}}$ and the backbone radius of gyration RG. (D) FES at 270 K in the space of two contact-map distances for conformations with more than two β-strand H-bonds. The insets in C and D are representative conformations of the native folded state, ${\beta }_N$; the folded state ${\beta }_2$; and a misfolded state, M. The FESs in C and D are built from a 100-ns-long run of REMD-MSSA. (E) Traces of ${\text{N}}_{\text{H}}$ and RG from a 220-ns-long bare MD simulation, where the blue is for ${\beta }_2$ and the green is for ${\beta }_N$.

Fig. S10.

Trajectories of potential energy, β-strand HB number, backbone radius of gyration, and $C_{\alpha }$ rmsd, from 220-ns MD runs. The blue is for state ${\beta }_2$ and the green is for state ${\beta }_N$.

Fig. S11.

Free energy barriers of the model in Eq. S23 at different α.

Fig. S12.

Probability distribution (Left) and free energy (Right) at $k_{B}T=0.025$. MD steps are ${10}^6$. The black line is analytical results and the red circles are from simulation.