Scrutinizing molecular mechanics force fields on the submicrosecond timescale with NMR data

Abstract

Protein dynamics on the atomic level and on the microsecond timescale has recently become accessible from both computation and experiment. To validate molecular dynamics (MD) at the submicrosecond timescale against experiment we present microsecond MD simulations in 10 different force-field configurations for two globular proteins, ubiquitin and the gb3 domain of protein G, for which extensive NMR data is available. We find that the reproduction of the measured NMR data strongly depends on the chosen force field and electrostatics treatment. Generally, particle-mesh Ewald outperforms cut-off and reaction-field approaches. A comparison to measured J-couplings across hydrogen bonds suggests that there is room for improvement in the force-field description of hydrogen bonds in most modern force fields. Our results show that with current force fields, simulations beyond hundreds of nanoseconds run an increased risk of undergoing transitions to nonnative conformational states or will persist within states of high free energy for too long, thus skewing the obtained population frequencies. Only for the AMBER99sb force field have such transitions not been observed. Thus, our results have significance for the interpretation of data obtained with long MD simulations, for the selection of force fields for MD studies and for force-field development. We hope that this comprehensive benchmark based on NMR data applied to many popular MD force fields will serve as a useful resource to the MD community. Finally, we find that for gb3, the force-field AMBER99sb reaches comparable accuracy in back-calculated residual dipolar couplings and J-couplings across hydrogen bonds to ensembles obtained by refinement against NMR data.

Figure 1

Comparison of R_RDC computed from MD ensembles generated with various force fields and two different electrostatic treatments (see Methods). No error bars are given in this figure, as only a single trajectory of 1 μs was available; refer to Fig. 3 to judge the variation in the respective R_RDC.

Figure 2

Fit of RDC data (R_RDC^x) computed from different-sized windows between x = 1 ns and x = 250 ns. (a and b) R_RDC^x resolved by starting position within the 1000 ns trajectories gb3_CHARMM22_PME and UBI_AMBER03, respectively. (c) For each trajectory (for numbering on x axis, use key given in legends of Fig. 3), R_RDC^x has been averaged over all possible windows. (d) For each trajectory (for index, see Fig. 3), the ratio R_RDC^x/R_RDC²⁵⁰ has been averaged over all possible windows. For most simulations, the ratio is close to 1 with x = 50 ns and x = 100 ns. This demonstrates that no significant improvement in R_RDC is gained by going to the longer averaging window of x = 250 ns.

Figure 3

Time-resolved R_RDC for MD simulations with various force fields. The residual dipolar couplings were computed from the respective MD trajectory in overlapping windows of 50 ns length. The figure key shows the simulation index used on the x axis of Fig. 2d in parentheses. The results shown here are computed for NH RDCs. Results for CH and NC RDCs are shown in the Supporting Material.

Figure 4

Comparison of ^3hJ_exp couplings across hydrogen bonds for proteins ubiquitin (left) and gb3 (right) with couplings back-calculated from MD simulations carried out with force-fields AMBER03 (a and d), AMBER99sb (b and e), and OPLS/AA (c and f). The second panel compares ^3hJ back-calculated from MD simulations that were restarted from the crystal structure every 50 ns. (Dashed lines) Mean ± 0.25 Hz.

Figure 5

PCA of ubiquitin ensembles. For each panel, PCA is carried out over the combined set of structures taken from 2k39 and the respective force-field ensemble. The projection to principal components 1 and 2 (pc #1 and pc #2) are plotted on x and y axis, respectively. The coloring depicts the 50 ns window average of R_RDC for a given snapshot (compare to Fig. 3); all panels use the same color-scale with blue for the lowest R_RDC values and red for the highest R_RDC values. The first principal component (pc #1, x axis) contributes >30% to the overall motion for all trajectories.

Figure 6

Snapshots from ubiquitin ensembles (green) and x-ray structure (1ubi, gray). Arrows are shown between residue pairs whose hydrogen bond is an outlier in Fig. 4. (a) Ubiquitin structure taken from the high R_RDC region of the OPLS/AA-PME ensemble. The projection of the selected structure to the PCA coordinates shown in Fig. 5 is (1.03; 0.45); see red/orange cluster in upper-right corner in that figure. (b) Ubiquitin structure taken from the high R_RDC region of the AMBER03 ensemble. The two selected frames project into the green cluster in Fig. 5, with (0.38; 1.24) and (1.23; 0.23), for the dark-green and light-green structures, respectively.