Effects of number of parallel runs and frequency of bias-strength replacement in generalized ensemble molecular dynamics simulations

Multicanonical MD

The McMD method efficiently explores the conformational space of a molecular system, by applying a biasing energy term. The Hamiltonian H of the system is (1) $$H=K+E_{\text{mc}},$$ where K and E_mc denote the kinetic energy and multicanonical energy, respectively. E_mc is defined as follows: (2) $$E_{\text{mc}}=E+RT\ln {\text{P}}_{\text{c}}\left(E,T\right),$$ where E is the potential energy, and the second term corresponds to the bias potential. R is the gas constant, and P_c(E, T) denotes the canonical distribution at the temperature T: (3) $${\text{P}}_{\text{c}}\left(E,T\right)=\frac{\text{n}\left(E\right)\exp \left(-\frac{E}{RT}\right)}{{\text{Z}}_{\text{c}}\left(T\right)},$$ where n(E) denotes the density of states, and Z_c(T) is the partition function of the canonical distribution at the temperature T. With this definition, the potential energy distribution of an ensemble obtained from the McMD, or the multicanonical distribution, P_mc(E), becomes uniform: (4) $$\begin{array}{c}{\text{P}}_{\text{mc}}\left(E\right)=\frac{\text{n}\left(E\right)\exp \left(-\frac{E_{\text{mc}}}{RT}\right)}{{\text{Z}}_{\text{mc}}\left(T\right)}\\ =\frac{\text{n}\left(E\right)\exp \left(-\frac{E}{RT}\right)}{{\text{P}}_{\text{c}}\left(E,T\right){\text{Z}}_{\text{mc}}\left(T\right)}=\frac{{\text{Z}}_{\text{c}}\left(T\right)}{{\text{Z}}_{\text{mc}}\left(T\right)}=const.\end{array}$$

As a result, the McMD method performs a random walk in the potential energy space and generates a uniform distribution of potential energy in a resultant ensemble. After a multicanonical ensemble is obtained, a canonical ensemble at any temperature in a sampled energy range can be generated by reweighting the probability of existence of each snapshot.

Equations (3) and (4) include an analytical form of n(E), which is usually unknown a priori. Therefore, n(E) is approximated as a parametric function, e.g., the polynomial function, and its parameters are estimated by iterations of McMD simulations to make P_mc(E) near-uniform. In the ith iteration, the bias potential is calculated using Eq. (2) with the canonical distribution obtained from the (i–1)th iteration, i.e., ${\mathrm{P}}_{\mathrm{c}}^{i-1}\left(E,T\right)$. As the result of the ith iteration, we obtain ${\mathrm{P}}_{\mathrm{m}\mathrm{c}}^i\left(E\right)$. ${\mathrm{P}}_{\mathrm{c}}^i\left(E,T\right)$ can be calculated as (5) $${\mathrm{P}}_{\mathrm{c}}^i\left(E,T\right)={\mathrm{P}}_{\mathrm{m}\mathrm{c}}^i\left(E\right){\mathrm{P}}_{\mathrm{c}}^{i-1}\left(E,T\right).$$

See Higo et al. (2012) for details.

Virtual-system coupled McMD

Virtual-system coupled McMD (V-McMD) introduces a virtual system, which interacts with the molecular system, and the multicanonical ensemble is calculated for the entire system consisting of these two subsystems (Higo, Umezawa & Nakamura, 2013). In practice, this method can be roughly interpreted as a combination of McMD and the simulated tempering method. The simulated tempering method replaces the system temperature with the Metropolis criterion and performs a canonical simulation until the next replacement trial. On the other hand, in V-McMD, the potential energy space is split into several regions (Fig. S1), and the molecular system is trapped in one of these regions. With a certain time interval (t_VST), the molecular system replaces the region to be trapped. The state variable governing which region traps the molecular system is called the “virtual state,” and the system defined by the virtual state is called the “virtual system.” The energy range of each virtual state is defined to be overlapped with the adjacent virtual states. When the molecular system has the potential energy E_k in the overlapped region of the ith and (i + 1)th virtual states, the state transition between these two virtual states can occur. Because this transition does not change the atomic coordinates or potential energy, the Metropolis criterion of this state transition is always satisfied. The time interval of virtual-state transitions (t_VST) should be determined arbitrarily by users. See Higo, Umezawa & Nakamura (2013) for details.

Trivial-trajectory parallelization

According to the theory of TTP, trajectories of multiple independent McMD runs with the same molecular system and different initial conditions can be treated as a single trajectory of an McMD simulation by concatenating the trajectories in an arbitrary order. This theory requires the condition that the initial coordinates of each run are sampled from the multicanonical distribution. Because the initial coordinates of production runs can be obtained from the near-uniform potential energy distribution generated by iterative simulations, it is expected that this condition holds. The McMD method with the TTP theory, which is called the TTP–McMD method, can be considered as a hybrid Monte Carlo sampler, by assuming that the system transitions from the last snapshot of the ith run (the microscopic state m_il) to the first snapshot of the jth run (the microscopic state m_jf) via a Monte Carlo step (Fig. S2). See Ikebe et al. (2010) for details.

Simulation protocol

We studied the three molecular systems, which are Trp-cage, PGA8, and PGA20 in an explicitly solvated cubic periodic boundary cell, by using the TTP-V-McMD method. Random coil structures of Trp-cage, PGA8, and PGA20 were constructed using the Modeller software (Webb & Sali, 2016) without any template. The termini of the PGAs were capped with acetyl and N-methyl groups, and the termini of the Trp-cage were not capped. Each of these molecular models was plased into a cubic box filled by water molecules; the number of water molecules were 5,097, 2,879, and 3,800 for Trp-cage, PGA8, and PGA20, respectively. In addition, a Cl^– ion was added to the Trp-cage model to cancel the net charge of the system. The net charge of the PGA models was zero because all the Glu residues were protonated.

The system was relaxed by using the GROMACS software (Pronk et al., 2013). Energy minimizations were successively applied using the steepest descent and conjugate gradient methods. Then, an MD simulation under a constant-pressure ensemble with the Berendsen barostat was performed for 1 ns. In the first half of the simulation, gradual heating from 10 to 300 K was applied. In the simulation, the positions of the heavy atoms of the Trp-cage were restrained, the bond lengths were not constrained, and the integration time step (Δt) was 0.5 fs. Subsequently, an additional constant-pressure relaxation was applied for 1 ns with Δt = 2.0 fs, and the covalent bonds to hydrogen atoms were constrained using the LINCS method (Hess et al., 1997; Hess, 2008). The final configuration of each model was used for the TTP-V-McMD simulations. The cell dimensions of these configurations were 54.0378, 44.6116, and 49.1174 Å for Trp-cage, PGA8, and PGA20, respectively.

For each model, the following steps were performed using our MD simulation program, which is called myPresto/omegagene and is tailored for GE simulations (Kasahara et al., 2016). The protein conformation was randomized with a constant-temperature simulation at 800 K. By using 30 snapshots taken from a trajectory with an interval of 300 ps, 30 independent runs were simulated with a gradual decrease in the temperature from 629 to 296 K to estimate the density of states. Successively, the TTP-V-McMD simulations were iteratively performed while updating the estimation of the density of states (Higo, Umezawa & Nakamura, 2013). A total of 84 production runs were performed (N_run = 84) for each of three different interval times for the virtual-state transitions (t_VST) meaning the interval times for bias-potential replacement: t_VST = 0.002, t_VST = 0.2, and t_VST = 20 ps. The simulation length of each run (t_run) was 50 ns except for the Trp-cage model with t_VST = 0.2 ps, t_run, of which the simulation length was 200 ns. In total, 50.4 μs of trajectories were simulated as production runs. The virtual system was divided into seven states that cover the energy range corresponding to the canonical distribution from 296 to 629 K. The velocity scaling method (Berendsen et al., 1984) was applied to maintain the system temperature.

For the potential parameters, the AMBER ff99SB-ILDN force field (Lindorff-Larsen et al., 2010), the ion parameter presented by Joung & Cheatham (2008), and the TIP3P water model (Jorgensen et al., 1983) were applied. The electrostatic potential was calculated using the zero-multipole summation method, which is a non-Ewald scheme (Fukuda, 2013; Fukuda, Kamiya & Nakamura, 2014). The zero-dipole condition with the damping factor α = 0 was used (Fukuda, Yonezawa & Nakamura, 2011; Fukuda et al., 2012).

Comparison of simulated ensembles among different settings

On the basis of the trajectories obtained from of the TTP-V-McMD production runs, the effects of the simulation conditions, i.e., the time interval for bias-strength replacement (t_VST), the number of independent runs (N_run), and the simulation time of each run (t_run), were assessed.

For the Trp-cage model, we analyzed the FEL for various conformational ensembles on the basis of the two structural parameters: the root-mean-square deviation (RMSD) of Cα atoms from the native conformation (PDB ID: 1L2Y, model 1), which is denoted as RMSD_native, and the radius of gyration (R_g). The FEL is visualized as the map of the potential of mean forces (PMF) on the plane defined by these two parameters. We defined the reference ensemble as the ensemble calculated for the conditions of t_run = 200 ns, N_run = 84, and t_VST = 0.2 ps, because it is expected to have the highest reliability owing to its abundance of samples (it comprises a total of 16.8 μs of simulations). The FELs analyzed in various conditions were compared with the reference FEL with regard to the Pearson correlation coefficient of the PMF (PCC_PMF). To calculate the PCC_PMF for a pair of FELs, bins without samples in one of the two FELs were ignored. In addition, the probability of the existence of the native conformations in each ensemble (P_native) was measured to characterize each ensemble. The native conformations are defined as the conformations with RMSD_native ≤ 2.0 Å.

For the PGA models, the FELs were analyzed using principal component analysis (PCA) based on the Cα–Cα distances (28 and 190 dimensions for PGA8 and PGA20, respectively). The PCAs were performed using aggregations of trajectories with all the three t_VST conditions for each model. For each t_VST condition, the ensemble calculated from the entire trajectory (t_run = 50 ns and N_run = 84) was considered as the reference ensemble. The FELs were compared with regard to PCC_PMF, similar to the Trp-cage case.

To assess the effects of N_run and t_run, PCC_PMF (and P_native for the Trp-cage model) were calculated for ensembles with subsets of the reference trajectories. Because there are many possibilities to pick N_run runs from 84 runs and t_run-length trajectories from the entire set of trajectories, we analyzed them by using the bootstrap approach. We constructed an ensemble by taking a random sampling of N_run runs from 84 runs with replacement and repeated it 100 times. The statistics over the 100 ensembles were analyzed via simulation with this N_run setting. This process was repeatedly performed for N_run = 1, 2, …, 84. For the case of t_run, the trajectories were split into 5-ns bins, and an ensemble was constructed by taking a random sampling of t_run/5 bins with replacement. We confirmed that the results of the bootstrap analyses with 100 and 200 samples were consistent (Fig. S3).

The sampling efficiency was measured in terms of the frequency of traversals between low- and high-energy regions, which were defined as the ranges (E_min, E_low) and (E_high, E_max), respectively. Here, E_min and E_max denote the minimum and maximum potential energies in all the trajectories, respectively, and E_low and E_high are defined as follows.

X is an arbitrary parameter in the range of 0–0.5. We assessed X = 0.2 and 0.3. The traversal frequency F_travers^E was calculated as the number of traversals between the two energy regions during 1.0 ns. The traversal frequencies of RMSD_native and R_g (F_travers^RMSD and F_travers^Rg, respectively) were also analyzed.

FEL of folding–unfolding equilibrium of Trp-cage

For the Trp-cage model, we performed 34 iterations of TTP-V-McMD simulations while updating the estimation of the density of states, n(E), and obtained a near-uniform energy distribution (Fig. S4). On the basis of this estimation, we performed production runs with N_run = 84, t_run = 200 ns, and t_VST = 0.2 ps. This is called the reference setting hereinafter. The resultant canonical ensemble reweighted at 300 K is referred to as the reference ensemble.

The FEL of the reference ensemble projected on the RMSD_native–R_g plane is shown in Fig. 1A. The most stable basin corresponds to the native structure consisting of an α-helix at the N-terminus, a 3₁₀-helix at the middle, and a loop region at the C-terminus (the secondary structural elements were recognized by using the DSSP software) (Kabsch & Sander, 1983). For example, the RMSD_native of one of the most probable structures in this basin was 0.994 Å (Fig. 1B). The energy barrier (approximately 3.3 kcal/mol) was observed at RMSD_native ≈ 3 Å in a low-R_g regime. Around this barrier, the 3₁₀-helix at the middle of the peptide chain was partially deformed; this deformation can be the first step of an unfolding process (Fig. 1C). The details of the unfolding pathway are not discussed in this paper. The second basin was widely spread around RMSD_native = 4–7 Å and R_g = 7–9 Å. This corresponds to the unfolded state, and examples of the unfolded structures taken from this basin are shown in Figs. 1E and 1F. The difference in the PMF between the bottoms of the first and second stable basins was 1.014 kcal/mol, and the population of the native conformations (P_native) was 22.37%. The landscape is qualitatively similar to that calculated using the REMD method reported by another group (Day, Paschek & Garcia, 2010). Our TTP-V-McMD simulation successfully identified the native structure as the most stable basin in the energy landscape, by using the reference setting.

Effects of simulation time for each run

The FELs of the Trp-cage model were drawn for a variety of t_run values under the condition of N_run = 84 and compared with the reference FEL. The FELs based on the trajectories of 0–25, 0–50, and 0–100 ns are shown in Figs. 2A–2C, respectively. The overall geometries of these FELs were qualitatively similar to the reference (Fig. 1A); their PCC_PMF values were 0.936, 0.936, and 0.994, respectively. The bootstrap statistics of PCC_PMF for each t_run value are summarized in Fig. 2D. For t_run = 200 ns, the bootstrap average and the standard deviation (SD) of PCC_PMF were 0.990 and 0.007, respectively. Even in the worst case among 100 randomly generated ensembles with t_run = 200 ns, PCC_PMF was 0.966. From this condition, a decrease in t_run yielded a slow decay of PCC_PMF, and PCC_PMF reached 0.9 at t_run ≈ 30 ns, which corresponds to 15% of the samples in the reference. Further decreasing t_run resulted in a steep decrease of PCC_PMF. Along with the decrease of the bootstrap average of PCC_PMF, the SD was increased. This means that an insufficient simulation time causes a loss of robustness of the results.

In contrast to the fact that the PCC_PMF decays in a shorter t_run than the reference, the balance between the folded and unfolded states (P_native) was almost constant regardless of t_run (Fig. 2E); the bootstrap average of P_native for t_run = 5–200 ns was in the range of 0.220 to 0.225. However, the SD of P_native was reduced with the increase of t_run; the SDs of P_native at t_run = 5, 50, and 200 ns were 0.07, 0.02, and 0.008, respectively. The loss of robustness due to the insufficiency of the simulation time is demonstrated in terms of not only the similarity of the entire FEL but also the stability of the native fold.

Effects of number of independent runs

As in the previous subsection, the effects of the reduction of N_run on the FELs were assessed under the condition of t_run = 200 ns. Examples of FELs with N_run = 10, 21, and 42 are shown in Figs. 3A–3C, respectively; the PCC_PMF values were 0.637, 0.939, and 0.993, respectively. Although the positions and wideness of the basins were similar to the reference, the FELs with a smaller N_run were smoother and lacked small bumps on the landscapes. The bootstrap statistics of PCC_PMF for various N_run values (Fig. 3D) were similar to those for t_run (Fig. 2D). The quantity of the samples required for PCC_PMF ≥ 0.9 was approximately one-fourth of the reference (N_run ≈ 21). The average (and the SD) of PCC_PMF at N_run = 21 and 42 were 0.906 (0.07) and 0.956 (0.04), respectively. Larger N_run values are needed to obtain robust results.

Regarding P_native, the influence of the reduction of N_run (Fig. 3E) differed from that of the reduction of t_run (Fig. 2E). A lower N_run resulted in the underestimation of the population of native conformations. P_native reached at plateau for N_run ≥ 21. A certain number of runs was needed to obtain robust results, and t_run = 200 ns was too short to reach equilibrium with a small number of trajectories for this system.

Balance between simulation time and number of runs

The evaluation for various t_run values with N_run = 84 runs (Fig. 2) and that for various N_run values with t_run = 200 ns (Fig. 3) indicate that reducing t_run produced better results than reducing N_run if the cumulative simulation time (N_run × t_run) was the same. Figure 4 shows direct comparisons of the results, indicating that high-N_run conditions resulted in a higher PCC_PMF and more similar values of P_native to the reference, with lower SDs, than long-t_run conditions. In particular, the qualitative difference between the two strategies is shown by the mean of P_native. Reducing N_run resulted in the significant underestimation of the fold stability, but reducing t_run did not.

In addition, we performed bootstrap analyses for all the combinations of 40-t_run settings (5, 10, 15, …, 200 ns) and 21-N_run settings (4, 8, 12, …, 84). The average values of PCC_PMF and P_native in all the conditions are presented in Fig. 5 and Fig. S5. The PCC_PMF was proportional to log(N_run × t_run). While the trend of P_native is ambiguous, the use of a larger number of samples resulted in a higher P_native. In the case where only small amount of data was available, a lower ratio of t_run/N_run (purple plots in Fig. 5) yielded better results.

Effects of frequency of bias-strength replacement

The parameter t_VST controls the frequency of the bias-strength switching in the TTP-V-McMD method. We investigated the effects of this parameter by comparing the TTP-V-McMD simulations of the Trp-cage model under the three conditions—t_VST = 0.002, 0.2, and 20 ps—with t_run = 50 ns for N_run = 84.

Table 1 summarizes the frequency of traversals between high- and low-potential energy regimes (F_trv^E), as defined in Eqs. (6) and (7) with X = 0.3 and 0.2, as well as the frequency of traversals between RMSD_native (F_trv^RMSD) and R_g (F_trv^Rg). The simulations with a shorter t_VST resulted in faster traversals in the potential energy space, indicating that with a shorter t_VST, a wider potential energy range can be sampled in a shorter time. However, faster traversal in the potential energy space does not ensure faster transition of the protein conformation. For both X = 0.2 and 0.3, although the setting of t_VST = 0.002 ps yielded the highest F_trv^E, this condition did not yield a higher F_trv^RMSD and F_trv^Rg compared to when a longer t_VST was used. This result indicates that the relaxation of the conformation requires a longer time than that of the potential energy. If a strong bias is applied and the system takes a high-potential energy state, it can return to low-energy states before conformational changes. Therefore, a moderate speed for traversals in the potential energy space is ideal for efficient conformational sampling. In the case of X = 0.2, t_VST = 0.2 ps exhibited the most frequent conformational changes.

In addition, the resultant ensembles were slightly affected by the setting of t_VST. We analyzed P_native for ensembles of various t_run values with N_run = 84 using the bootstrap method (Fig. S6). The results for all three t_VST values showed similar trends, i.e., near-constant average values and the gradual decay of the SD with the increase of t_run. While t_VST = 0.2 ps showed a smaller P_native than the other two t_VST settings, the difference was smaller than the SD. On the other hand, higher SD values were observed in the following order: t_VST = 0.2 > 20 > 0.002 ps. This is consistent with the order of F_trv^RMSD and F_trv^Rg (Table 1). The result indicates that more frequent traversals between high- and low-RMSD_native conformations make it possible to explore a wider region of the conformational space; thus, the population of the native conformation decreases, and the SD increases.

Regarding the PCC_PMF with the reference setting (t_VST = 0.2 ps, N_run = 84, and t_run = 200 ns), the average PCC_PMF values at t_run = 50 ns differed among different settings of t_VST (Fig. S6). This indicates that changing t_VST yields subtle differences in the resultant ensemble. Regarding the balance between t_run and N_run, the trends were similar for all the settings of t_VST (Fig. S7).

Effects of system complexity: comparison with the PGA models

We performed the same analyses for the molecular models of PGA8 and PGA20. In contrast to Trp-cage, these peptides did not exhibit a particular fold. The FELs of both PGA8 and PGA20 were unimodal distributions, the basins of which consisted of a variety of collapsed conformations (Fig. 6 for t_VST = 0.2 ps). The ensembles included short secondary structural elements but they were unstable. Although the shape of the small bumps in the basins differed depending on the simulation conditions, the overall geometries of the FELs were similar (Fig. S8 for t_VST = 0.002 ps and 20 ps).

Regarding the balance between t_run and N_run, Fig. 7 shows the bootstrap averages of PCC_PMF between the ensemble calculated by the full-length trajectory (t_run = 50 ns and N_run = 84) and those calculated by the reduced trajectories. No clear differences were found between the PCC_PMF curve with reduced t_run and that with reduced N_run for both the PGA8 and PGA20 (Fig. 7 for t_VST = 0.2 ps; Fig. S9 for the other conditions). A small number of long simulations exhibited the similar efficiency as that of many short simulations. In addition, no significant differences were found between the results of PGA8 and PGA20. It is noteworthy that the conformational space of PGA20 is considerably wider than that of PGA8 and similar to that of Trp-cage, because the conformational space volume of polypeptides is determined primarily by their length. Therefore, we concluded that the effects of balance between t_run and N_run are determined by the complexity of the FEL (e.g., existence of the free-energy barrier) rather than the conformational space volume. An increase in the number of runs is more effective for a system with more complex FEL.

For the PGA models, the frequencies of traversals in the potential energy and R_g spaces (F_trv^E and F_trv^Rg, respectively) are summarized in Table 2. Both the PGA8 and PGA20 models yielded similar trends as the Trp-cage model (Table 1). Although frequent replacements of bias-potential strength enhanced the traversals in the potential energy space, they did not enhance the conformational changes in terms of R_g. This implies that the conformational changes are much slower than the potential energy changes even if there is no free-energy barrier exists in the landscape. However, in contrast to the Trp-cage case, the drawback of the frequent replacement, that is, slow traversals in the conformational space, is unclear in the case of PGA20.