Influence of various parameters in the replica-exchange molecular dynamics method: Number of replicas, replica-exchange frequency, and thermostat coupling time constant.

The replica-exchange molecular dynamics (REMD) method has been used for conformational sampling of various biomolecular systems. To maximize sampling efficiency, some adjustable parameters must be optimized. Although it is agreed that shorter intervals between the replica-exchange attempts enhance traversals in the temperature space, details regarding the artifacts caused by these short intervals are controversial. In this study, we revisit this problem by performing REMD simulations on an alanine octapeptide in an implicit solvent. Fifty different sets of conditions, which are a combination of five replica-exchange periods, five different numbers of replicas, and two thermostat coupling time constants, were investigated. As a result, although short replica-exchange intervals enhanced the traversals in the temperature space, they led to artifacts in the ensemble average of the temperature, potential energy, and helix content. With extremely short replica-exchange intervals, i.e., attempted at every time step, the ensemble average of the temperature deviated from the thermostat temperature by ca. 7 K. Differences in the ensembles were observed even for larger replica-exchange intervals (between 100 and 1,000 steps). In addition, the shorter thermostat coupling time constant reduced the artifacts found when short replica-exchange intervals were used, implying that these artifacts are caused by insufficient thermal relaxation between the replica-exchange events. Our results will be useful to reduce the artifacts found in REMD simulations by adjusting some key parameters.

The replica-exchange molecular dynamics (REMD) method is a powerful approach to investigate conformational ensembles of molecular systems. Although this method requires some adjustable parameters, there is no golden standard for setting of the parameters. Effects of the parameters on sampling efficiency and resultant ensemble have been argued. Our evaluation found that shorter intervals for replica-exchange attempts enhance conformational sampling but cause the artifacts due to insufficiency of thermal relaxation. In particular, the ensemble average of temperature deviates from the target temperature of the thermostat. The results of this evaluation will be helpful to adjust the parameters for REMD simulations.

Molecular dynamics (MD) simulations are a promising computational approach to investigate the microscopic behavior of molecular systems and have been successfully applied to various problems, e.g., drug discovery [1,2]. Although recent advances in computational technology have extended the scope of MD simulations [3,4], insufficiency of the simulation time scale is still considered a major limitation of this approach. Equilibration of a molecular system with a high degree of freedom requires a long time scale, which cannot be attained by conventional MD methods with current general-purpose computers [5-7].

To overcome this limitation, a variety of extended ensemble approaches have been extensively developed [8,9]. The replica-exchange MD (REMD) method [10], also known as the parallel tempering MD, is one of the generalized ensemble MD methods. The REMD method consists of performing several constant-temperature MD simulations on the same molecular system at different temperatures and exchanging the temperature of the thermostat between pairs of simulations based on the Metropolis criterion. Each simulation, or replica, traverses a wide range of temperatures, and thus, various molecular conformations can be efficiently sampled. Applying a re-weighting method to the multiple trajectories obtained at different temperatures generates a conformational ensemble in an equilibrium state at a target temperature [11,12]. In the past decades, many variant methods of REMD have been developed, e.g., Hamiltonian REMD [13], surface tension REMD [14], solute tempering [15,16], and mass scaling REMD [17]. REMD-based enhanced conformational sampling methods have been applied to investigate the conformational ensembles of a variety of molecular systems in an equilibrium state [18].

Since a major bottleneck of the REMD method is its high computational cost induced by the parallel computations of many replicas, optimizing the sampling efficiency is essential. In particular, optimizing the values of several adjustable parameters, such as the number of replicas (referred to as Nrep, hereafter) and the time interval between the replica-exchange attempts (referred to as tatt, hereafter), is not straightforward, and many evaluations of these parameters have been reported. Periole & Mark reported that for the conformational sampling of a heptapeptide with explicit water molecules, a short interval (tatt=0.1 ps) resulted in a high frequency of back exchange and high autocorrelation of the potential energy [19]. Abraham & Gready also suggested utilizing tatt ≥1.0 ps to avoid the autocorrelation of the potential energy [20]. Jani et al. simulated the folding–unfolding equilibrium of an engrailed homeodomain with three different combinations of tatt and simulation temperature and found that these settings impact the resultant ensembles [21]. On the other hand, Sindhikara et al. reported that short tatt enhance traversals in the temperature space without side effects based on simulations of short peptides with both implicit and explicit solvents [22,23], and Rosta & Hummer supported this conclusion [24]. Although it is clear that the shorter intervals enhance traversals in the temperature space, the recommended lower limit and artifacts generated are still debated.

Herein, we revisit the evaluation of different REMD parameters using a poly-alanine peptide with an implicit solvent model following the protocol presented by Sindhikara et al. [22]. To evaluate the artifacts generated by the short interval between replica-exchange attempts, we tested a wide range of tatt, including extremely frequent attempts (tatt=0.001, 0.005, 0.01, 0.1, and 1 ps). In addition, we considered five different numbers of replicas (Nrep=8, 12, 16, 20, and 32) as well as the thermostat coupling time constant τ, whose effects have not yet been reported. Since the previous study implies that the artifacts caused by short tatt arise because of insufficient thermal relaxation between the replica-exchanges [19], τ should have a key role. In total, 5×5×2=50 sets of REMD simulations were performed to assess the effects of these three adjustable parameters on the resultant ensemble and on the sampling efficiency. The obtained ensembles were characterized by the helix content, potential energy distribution, and temperature distribution. The simulation time required for a round trip in the temperature space was assessed to evaluate the sampling efficiency. As a result, though shorter tatt enhance traversals in the temperature space, which is consistent with the previous works, they also tend to generate artifacts in the ensembles. In addition, longer τ caused more significant artifacts for shorter tatt. This result encourages the use of longer tatt to minimize artifacts, and there is therefore a trade-off between the accuracy of the ensemble and the sampling efficiency.

Methods

Simulation Settings

The REMD simulations were performed on a poly-alanine octapeptide with an implicit solvent model under various conditions (tatt, Nrep, and τ). As the initial structure for all simulations, an ideally extended structure of the poly-alanine octapeptide was built using the tLEaP software included in the AMBER suite. The N- and C-termini were capped with acetyl and N-methyl groups, respectively. For the potential calculations, the AMBER parm99SB force field [25] and the generalized Born (GB) implicit solvent model with the OBC(II) parameter [26] were utilized. The cut-off distance was large enough to cover the entire system. Covalent bonds involving hydrogen atoms were constrained by using the LINCS method [27,28]. The integration time step (Δt) was 1.0 fs. The Nosé–Hoover thermostat [29,30] was applied to generate NVT ensembles. The target temperature of each replica was defined as a geometric series ranging from 300 to 800 K.

We performed REMD simulations with 50 sets of conditions; i.e., the combination of five time intervals between replica-exchange attempts (tatt=0.001, 0.005, 0.01, 0.1, and 1 ps), five total number of replicas (Nrep=8, 12, 16, 20, and 32), and two coupling time constants for the Nosé-Hoover thermostat (τ=0.2 and 2 ps). The simulation time was 10 ns for each replica following the previous evaluation by Sindhikara et al. [22]. In total, 8.8 μs of REMD simulations were examined. All the simulations were performed using the Gromacs (ver. 4.6.1) software [31] with a modification for the velocity-rescaling of the thermostat [32]. In order to analyze the raw distribution at each thermostatting temperature, we did not apply re-weighting methods, e.g., the weighted histogram analysis [11] and multistate Bennett acceptance ratio methods [12].

Definitions

Herein, several parameters are used to describe the structural ensembles and simulation conditions. This sub-section summarizes all the terms used in this paper.

The adjustable parameters for the REMD simulations are tatt, Nrep, Tk, and τ. tatt denotes the time interval between replica-exchange attempts. Nrep is the number of replicas. Tk is the target temperature of the k-th thermostat, where k=1, 2, ..., Nrep. τ represents the coupling time constant for the Nosé–Hoover thermostat. The REMD simulations were characterized by the parameters Pacc, tacc, tround, FH, FG, Rg, T, and V. Pacc denotes the acceptance ratio of the replica-exchange attempts. tacc is the expected time for accepting the replica-exchange attempts, calculated as tacc=tatt/Pacc. tround is the time required for a round trip in the temperature space; the round trip is defined as the traversals from a lowest/highest temperature (k=1 or k=Nrep) to a highest/lowest and back to the initial (lowest/highest) temperature. FH and FG indicate the ensemble average of fraction of α- and 310-helices in an ensemble at 300 K, respectively. These secondary structural elements were assessed using the DSSP software [33], and the fractions were calculated as the ratio of the number of residues with a given secondary structural element to the total number of residues. The radius of gyration (Rg) is measured with the g_gyrate module of the Gromacs package. T and V are the temperature and potential energy, respectively, of a sampled snapshot. The operators < > and σ( ) denote the average and standard deviation, respectively. In particular, the statistics over an ensemble at the k-th temperature are described as < >k and σ( )k.

Results and Discussion

The first three sub-sections describe the results for τ=2.0 ps. The first sub-section, Traversals in temperature space, discusses the sampling efficiency of the REMD simulations in terms of the time required for a round trip in the temperature space. The characteristics of the resultant ensembles are discussed in the second and third sub-sections, Helix content and radius of gyration and Potential energy and temperature distributions, respectively. The differences between the simulations for τ=2.0 and 0.2 ps are discussed in the fourth sub-section, Effects of the thermostat coupling constant.

Traversals in temperature space

For the REMD simulations with τ=2.0 ps and combinations of tatt and Nrep, the sampling efficiencies were assessed in terms of the average time for a round trip in the temperature space (tround; Table 1). In general, the shorter tatt exhibited faster traversals in the temperature space, consistent with previous reports [22,23]. In addition, for each time interval tatt, a distinct highest-efficiency value of Nrep was observed; the Nrep values that exhibit the fastest round trips for each tatt were 32, 16, 12, 12, and 12 for tatt=0.001, 0.005, 0.01, 0.1, and 1 ps, respectively. The most efficient Nrep value tends to decrease with an increase in tatt and reaches a plateau at tatt ≥0.01 ps. Although the larger Nrep showed higher acceptance ratios owing to larger overlaps of the energy distribution between neighboring replicas (Table 2), the larger Nrep require larger number of acceptances for a round trip. Thus, unless extremely frequent replica-exchange attempts were performed, the higher Nrep did not show better performance. An acceptance ratio of Pacc~0.5 seems to be the best condition for tatt ≥0.01 ps. Note that since Pacc depends on the degree of freedom of the system, the most efficient conditions can vary for different systems.

Table 1

The acceptance ratio for replica-exchange attempts, Pacc

	τ	2	2	2	2	2	0.2	0.2	0.2	0.2	0.2
Nrep	tatt	1	5	10	100	1000	1	5	10	100	1000
8		0.257	0.236	0.253	0.247	0.249	0.259	0.234	0.253	0.246	0.249
12		0.470	0.446	0.471	0.464	0.466	0.472	0.448	0.471	0.462	0.465
16		0.597	0.578	0.599	0.587	0.591	0.597	0.579	0.600	0.590	0.592
20		0.674	0.662	0.679	0.669	0.669	0.675	0.663	0.679	0.670	0.672
32		0.790	0.791	0.798	0.795	0.795	0.791	0.791	0.798	0.795	0.796

Table 2

Average time for a round trip in the temperature space, tround (ps)

	τ	2	2	2	2	2	0.2	0.2	0.2	0.2	0.2
Nrep	tatt	1	5	10	100	1000	1	5	10	100	1000
8		38.47	82.02	81.36	204.95	983.39	17.33	36.82	52.65	198.77	777.28
12		33.77	62.05	69.95	199.36	680.41	14.96	34.94	45.44	202.64	797.09
16		33.06	59.78	74.86	219.64	829.89	16.57	35.38	46.39	195.76	859.40
20		33.21	64.18	78.60	219.28	857.71	17.50	37.48	54.94	210.11	890.37
32		28.37	65.33	84.64	279.36	999.38	16.53	42.73	63.03	257.33	1102.94

Regarding the expected time to accept a replica-exchange attempt, log(tacc) was proportional to log(tround) for each value of Nrep (Fig. 1). In addition, the larger Nrep tended to cause longer tround for the same tacc because of the larger acceptance number required for a round trip.

Figure 1

Relationship between the expected time for an acceptance (tacc; the horizontal axis) and the average time required for a round trip in the temperature space (tround; the vertical axis) with (A) τ=2 ps and (B) τ=0.2 ps. The conditions Nrep=8, 12, 16, 20, and 32 are shown in red, green, blue, cyan, and purple, respectively. The lines indicate least-square fitted functions.

Helix content and radius of gyration

The fractions of α-helix conformations (FH) in the ensembles at 300 K are summarized in Figure 2. The three longest intervals between replica-exchange attempts (tatt ≥0.01 ps) exhibited similar FH values regardless of the number of replicas (FH=0.838~0.854). On the contrary, when tatt=0.001 ps, larger Nrep resulted in higher FH. Although the FH values when tatt=0.005 ps also depended on Nrep, the trend was the opposite (smaller Nrep resulted in lower FH), and the deviations of the FH values with respect to those obtained with longer tatt were subtle. In particular, when tatt=0.005 ps and Nrep=32, the fraction of α-helix conformations (FH=0.870) was similar to those with longer tatt. While the three longest intervals (tatt ≥0.01 ps) exhibited similar FH values, subtle differences were observed; the Pearson correlation coefficient (PCC) of the FH values for the five Nrep conditions between tatt=0.01 and 1 ps was −0.187 and that between tatt=0.1 and 1 ps was 0.902. These results indicate that tatt has a significant impact on the conformational ensemble, and longer tatt produce more robust resultant ensembles in terms of FH with a wide range of Nrep. The effects of tatt on FH converged when tatt ≥0.1 ps. Note that the 310-helix content, FG, was negatively correlated with FH (the PCC between FH and FG was −0.99; Supplementary Fig. S1).

Figure 2

α-helix content in the resultant ensembles with (A) τ=2 ps and (B) τ=0.2 ps. The horizontal and vertical axes indicate Nrep and FH, respectively. The conditions tatt=0.001, 0.005, 0.01, 0.1, and 1 ps are shown in red, green, blue, cyan, and purple, respectively. The error bars indicate the standard errors among 20 sub-trajectories generated by uniformly dividing the last 9 ns of the original trajectory.

An experimental study of a poly-alanine heptapeptide suggested that the poly-proline II conformation was the dominant structure [34], which is inconsistent with our simulation results that indicate high helix tendencies (FH~0.85). However, a direct comparison with experimentally measured helix contents is not straightforward because of the limitation of the current force fields and implicit solvent models [35]. Even with explicit solvent models, the radius of gyration of denatured peptides tends to be underestimated by some standard force fields [36,37], and issues caused by finite-size effects arising from the periodic boundary conditions have also been raised [38-41]. In addition, experimental quantification of helix content of oligopeptides is not straightforward [42]. Nevertheless, this study focuses only on the relative differences among the REMD configurations, and consistency with experimental measurements is out of the scope of this work.

We also analyzed the peptide conformations in terms of Rg and found that Rg was negatively correlated with FH (PCC was −0.99; Supplementary Fig. S2). In addition, we confirmed the convergence of Rg by calculating its cumulative average as a function of the simulation time (Supplementary Fig. S3). Rg values converged in the simulation time, and there were no clear differences among simulations carried out using various REMD conditions.

Potential energy and temperature distributions

We evaluated the potential energy and temperature distributions of the simulated ensembles at the lowest temperature (T1=300 K). For tatt=0.001 ps, the potential energy distribution is shifted compared with the other four tatt values; smaller and larger Nrep showed higher and lower potential energies, respectively (Fig. 3A, B and Supplementary Fig. S4A, B, C). The difference in the ensemble averages of the potential energy 1 with the different parameters showed a trend similar to the α-helix content (FH); the PCC between 1 and FH was −0.99 (Figs. 2A and 3C).

Figure 3

Potential energy distributions in ensembles at 300 K with (A, B, C) τ=2 ps and (D, E, F) τ=0.2 ps. The conditions tatt=0.001, 0.005, 0.01, 0.1, and 1 ps are shown in red, green, blue, cyan, and purple, respectively. (A, B, D, E) The distribution of the potential energy V in each ensemble are shown for (A, D) Nrep=8 and (B, E) Nrep=32. (C, F) Ensemble averages of the potential energy 1 for each Nrep and tatt conditions are summarized.

For the temperature distribution, the most significant shift was observed for tatt=0.001 ps with Nrep=32 (Fig. 4B). Its ensemble average of temperature 1 was 293.13 K despite that the target temperature of the thermostat was 300 K. Contrarily to the potential energy distribution, no clear shift of the temperature distribution was observed for Nrep=8 with tatt=0.001 ps (Fig. 4A); the distributions obtained for the other Nrep values are shown in Supplementary Figure S5A, B, C. The ensemble averages of temperatures 1 for various conditions are summarized in Figure 4C. This clearly shows that there is a problem if the frequency of the replica-exchange attempts is very high (tatt=0.001 ps); at this value, the ensemble average of temperatures exhibited a distinct behavior and largely deviated from the target temperature of the thermostat. For longer tatt, the deviations from the target temperature were smaller; the maximum errors from the target temperature were 7.07, 1.28, 0.91, 0.28, and 0.06 K for tatt=0.001, 0.005, 0.01, 0.1, and 1 ps, respectively. It is noteworthy that there were detectable differences even between the two largest tatt values.

Figure 4

Temperature distributions in ensembles at 300 K with (A, B, C) τ=2 ps and (D, E, F) τ=0.2 ps. Notations are the same with Figure 3. (A, B, D, E) The distribution of the temperature T in each ensemble are shown for (A, D) Nrep=8 and (B, E) Nrep=32. (C, F) Ensemble averages of the temperature 1 for each Nrep and tatt conditions are summarized.

Thermostatting errors were also observed in replicas with higher temperatures. The ratios of the ensemble average of the temperature over the target temperature k/Tk for various conditions are summarized in Figure 5 and Supplementary Figure S6. Interestingly, these ratios oscillated along the temperature index, k. In general, the ratios of the (k+2i)-th temperatures, where k and i denote positive integers, were similar. This zigzag profile was common, irrespective of Nrep, and the wavelength of the oscillation (of the temperature) was always equal to two replicas. This artificial oscillation tended to decay with an increase in tatt. We speculate that this oscillation is due to the REMD protocol, which performs two types of replica-exchange attempts consecutively: (i) exchanging the 1st–2nd, 3rd–4th, ... replica pairs and (ii) exchanging the 2nd–3rd, 4th–5th, ... replica pairs. If the system temperature is not properly relaxed in the time interval between the replica-exchange attempts, such a systematic trial scheme can generate the oscillation.

Figure 5

Ratios of the ensemble average of temperatures to the target temperature at k-th thermostat k/Tk with (A, B, C) τ=2 ps and (D, E, F) τ=0.2 ps for each Nrep condition: (A, D) Nrep=8 , (B, E) Nrep=16, and (C, F) Nrep=32. The results with for tatt=0.001, 0.005, and 0.1 ps are shown in red, green and cyan, and those with tatt=0.01 and 1 ps are not shown for clarity.

Effects of the thermostat coupling constant

The coupling constant τ of the Nosé–Hoover thermostat is one of the key parameters of MD simulations. We examined the effects of τ by performing REMD simulations with two settings, τ=0.2 and 2 ps. First, the shortest τ enhanced traversals in the temperature space (Table 2), but τ did not affect the acceptance ratio of the replica-exchange attempts Pacc (Table 1). This behavior can be interpreted as follows. After a replica-exchange attempt is accepted, the two exchanged replicas should have similar temperatures. When successive replica-exchange attempts occur before thermal relaxation, the attempts to move these replicas to the previous temperature should have high acceptance probabilities. Namely, a longer τ, indicative of a slower thermal relaxation, increases the probability of repeating the round trip between two neighboring temperatures. Therefore, although the Pacc values were similar for both τ conditions, the time required for a round trip in the temperature space differed. Similarly, frequent back exchanges were also reported by Periole & Mark for short tatt [19]. Second, a linear relationship between the logarithm of tacc and that of tround was also observed for τ=0.2 ps (Fig. 1B). Third, the α- and 310-helix contents for the two τ conditions showed similar profiles (Fig. 2 and Supplementary Fig. S1). However, some differences can be found in their potential energy distributions. For τ=0.2 ps and Nrep=32, the distributions of V1 were similar for the various tatt settings (Fig. 3E), in contrast to the distributions when τ=2 ps (Fig. 3B). Fourth, although the thermostatting errors at Tk=300 K for the lowest tatt values were also observed for τ=0.2 ps, the errors were smaller than those for τ=2 ps (Fig. 4). However, an oscillation of the averaged temperatures was also observed (Fig. 5). In summary, the coupling time constant τ significantly affected the sampling efficiency and the ensemble generated by the REMD method, especially for short time intervals between the replica-exchange attempts. This suggests that artifacts arising from short tatt are related to insufficient thermal relaxation.

Conclusion

We evaluated the effects of parameters on REMD simulations, i.e., the number of replicas (Nrep), the time interval between replica-exchange attempts (tatt), and the coupling time constant for the thermostat (τ), in terms of the sampling efficiency and resultant ensembles of an alanine octapeptide. The results suggest that, although shorter tatt enhance traversals in the temperature space (Table 2), they also cause larger artifacts. Short tatt also shift the temperature distributions of the simulated ensembles from the target temperature of the thermostat (Fig. 4). In addition, the thermostatting errors artificially oscillate along with the thermostat-ID, k (Fig. 5). The potential energy distribution and helix content are also impacted by the value of tatt (Figs. 2 and 3). The other parameters, Nrep and τ, also significantly influence the generated ensembles whenever short time intervals are used.

Nrep is a key parameter that affects the sampling efficiency, which is measured by the time needed for a round trip in the temperature space (tround). Higher Nrep increase both the acceptance ratio (Pacc; Table 1) and the acceptance number required for a round trip; the former increases and the latter decreases the efficiency. In the molecular system studied here, Nrep=12 or 16 show the best performance in terms of tround unless tatt is too short (Table 2). Although higher Nrep require higher computational costs, those are not always effective to enhance the sampling. The logarithm of tround is proportional to that of tacc (Fig. 1), and higher Nrep result in larger tround for the same time interval for accepting the replica-exchange attempts (tacc).

Whereas various combinations of parameters were examined in this study, there are some other factors which can affect the results of REMD simulations, e.g., the system size, the degree of freedom of the system, the integration time step Δt, and the thermostatting algorithm. Note that Sindhikara et al. examined the effects of parameters with implicit and explicit solvents and got the same conclusions in both cases. The complexity of the method makes it difficult to determine optimal parameters, and it is not straightforward to find a universal rule for the parameter determination. Although Sindhikara et al. recommended 10–100 steps for the replica-exchange interval [23], our results imply that this setting is not always adequate. We demonstrated that employing a 10-steps interval (tatt=0.01 ps) generated artifacts. Whereas using a 100-steps interval (tatt=0.1 ps) clearly reduced the artifacts compared with shorter intervals, the ensemble still deviated from that obtained with a 1,000-steps interval (tatt=1 ps). Users should carefully choose this value to get a balance between sampling efficiency and tolerable artifacts.