Rafael Brüschweiller’s lab has for several years been engaged, with some success, in an effort to improve MD simulations to the point where they can predict S2 values that match NMR results for side chains. While a cold-eyed analysis of the correlations they’ve obtained so far (r values near 0.6) might not be very favorable, the comparisons aren’t that bad. Their most recent communication to JACS (citation at the end of the post), appearing online last week, displays a marked improvement in the correlations between simulation and experiment. The r values are still not that close to 1, but the current results appear to be a significant step forward.
So, what was done differently? Showalter et al. use an altered version of the AMBER99 forcefield that has a modified dihedral angle potential. This had good results in simulating the dynamics of backbone amide moieties, although MD has historically done a reasonably good job with these anyway. In this communication they simulated the side-chain motions of calbindin and compared them to experimental values. They calculated spectral densities J(ω) from their simulated correlation functions, though they are required in this case to make use of experimentally determined molecular correlation times. They do a strikingly good job of predicting the J(ω) at the Larmor frequency of deuterium and also at twice that frequency (figure filched from paper):
This really is an amazingly good job. Yet, as you can see from their figure 2 (a part of it is at right), the S2 values they obtain aren’t very close to those that are derived from NMR experiments. This is also reflected in the relatively poor agreement at J(0) (r=0.86). This seems to be very odd, because the magnitude of the spectral density at J(0) is very strongly dependent on τm, which they took from an NMR experiment. As is evident from the model-free expression for J(ω), the order parameter scales this term. Keep in mind that τm is typically on the order of 10-9 seconds while τe is on the order of 10-11 seconds—this means that the second term in the spectral density expression is negligible at J(0). Because they took their τm from experimental data, the decreased correlation at J(0) indicates that their correlation functions converged to inaccurate values.
Giving the data a once-over, it appears that their fitted order parameters were mostly high. It’s not clear to me why this should be so, except that over-constraint of the backbone may be affecting the side chains. From the supplementary information it appears that fits of backbone order parameters were also slightly higher in MD than experiment.
The most notable failure is not much help because it missed low. The significant outlier in the J(0) plot is threonine 45, shown in orange at right—this figure is made from PDB structure 3ICB, which was used in this simulation. The correlation function for this residue fails to converge, largely due to sampling of an alternate ψ angle. This behavior is consistent with the observations of low order parameters in that particular loop of the protein. From the crystal structure it appears that the hydroxyl moiety of Thr 45 is capping a helix (blue). It’s possible that the misbehavior of this particular residue is due to some miscalibration of the force-field that doesn’t accurately capture this capping interaction. The altered backbone dihedral angle potential may be overwhelming the hydrogen bonding interaction, resulting in the aberrantly low J(0) fit for this methyl group.
In general, the simulations for threonines and valines were not as accurate as those for other types of residues, which seems a little strange. These were also unusual in that they missed low, while, as I mentioned, on average residues tended to miss high. The branched nature of these amino acids causes some steric interactions with the backbone, so one would expect an improved potential to help these residues the most. However, if the altered potential is causing unwarranted excursions from the equilibrium structure, as seems to be the case with Thr 45, then valines and threonines, whose motions are at least partially controlled by steric interactions with the backbone, might be the most strongly affected.
I should stress that it’s not necessary for the simulation to produce too much backbone motion to get this result. If the backbone dynamics are the wrong kind of motion that could have this effect even if the model-free parameters for the backbone derived from the simulation appear to be accurate.
It isn’t terribly clear why an improved backbone potential should increase the correlation of side-chain order parameters between MD and NMR. Showalter et al. venture no explanation, and my own research and that of others hasn’t shown any particular linkage between backbone dynamics and side-chain dynamics, except in the case of alanine residues. It may be that a more accurate depiction of backbone motions contributes to a more accurate dynamic environment generally. Or, the motions of the backbone and side chains could be related in unexpected ways. An in-depth analysis of the simulation probing for these correlations could be very instructive. Regardless, these results are a significant, encouraging step towards using MD simulations to interpret the findings of NMR dynamics experiments.
Showalter, S. A.; Johnson, E.; Rance, M.; Bruschweiler, R. “Toward Quantitative Interpretation of Methyl Side-Chain Dynamics from NMR by Molecular Dynamics Simulations” J. Am. Chem. Soc. (Communication); 2007;ASAP Article.