Can MD match model-free?

As I mentioned in my last post, a major challenge in the interpretation of protein motion has been the poor correlation between dynamics information arising from computer molecular dynamics simulations and NMR relaxation experiments. MD simulations complement the order parameters of the model-free formalism by providing detailed descriptions of motions that model-free parameters describe in a general way. However, the output of an MD simulation is only to be trusted if its specific model matches the experimental observations from NMR. Historically this has not been the case, especially when it comes to the description of side-chain motions.

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 S² 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 S² 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.

Read the rest...

The model-free dynamics formalism of Lipari and Szabo

I've seen a lot of really interesting articles in the past several days, though as always I admit that many things I find fascinating are sleep aids to others. I was going to write about the most esoteric article first, but as I started writing out an explanation I realized I needed to split things into two posts. The next post will be about an article I read last week in JACS ASAP pre-publication. But this post is about an article that appeared in JACS when I was about four years old. In it, Giovanni Lipari and Attila Szabo introduced a way of describing the dynamic information obtained from NMR so useful and successful that it has not been supplanted even after a quarter of a century. In addition, the article is really wonderfully written, and if you are a professional then I highly recommend reading it (the citation is at the bottom of this post). If you are not a professional, then read this post instead.

Regrettably, it is impossible to talk about this without using mathematics, but I promise you will not be forced to endure anything taxing, just exponents and multiplication. The Lipari-Szabo model-free formalism derives much of its power from its relative simplicity, so I promise this will not be too painful, and for those of my family who are interested, you will learn what all that junk I talk about means. Practicing NMR spectroscopists will notice a few differences in terminology and some (hopefully slight) conceptual fudging in what follows, but hey, you should be reading the original article anyway.

So. Let's imagine we have a really simple chemical system, like a single bond between two atoms, which we will call A and X (see right). Now, unless our AX compound is very, very cold it will be moving around all the time. Exactly how fast it moves will depend on what is around it (AX moves quicker through water than through honey) and what temperature it is (AX moves faster when warm), but unless we are at absolute zero, it will move. We can imagine two kinds of motion. AX can move forward or sideways: this is called translational motion, and having mentioned it, we won't worry about it anymore. In fact, we will temporarily pretend that it does not even happen. Also, AX can tumble—for instance, it could rotate around some point in the center of the bond, or swing around like a clock hand, or wobble at the ends (as shown). This is the motion we will consider.

Suppose you had a system like this and you wanted to come up with a way to show how it was moving. Well, one thing you could do is draw a picture or make a movie showing how quickly it tumbled. However, this would be a pretty awkward way to convey the data and might not be very informative. It might be better to try to come up with a single number that tells us how much AX has tumbled. Well, one way to do that would be to monitor AX over time, starting at some random time i, and compare the bond vector at any time i+t to the bond vector at time i. We could assign to any time t a number representing how likely it is to find AX in the same orientation at t that it had at i. Obviously, early on, this probability would be very high, and the larger t got, the smaller the probability would get. An example graph is over there on the left. This is called a correlation function, because it represents the correlation between the AX orientation at time i and the orientation at time i+t.

Obviously, for any real system the correlation function will be bumpy and difficult to represent mathematically. However, it turns out that we can approximate a correlation function reasonably well using an exponential decay, that is, a function based on a negative exponent of Euler's number. Thus we have:
where C(t) is the correlation function, t is the time (our x axis), and τ_m is a time constant, the value of which will depend on the speed with which AX tumbles. And now we have a single number that tells us about the tumbling. This number, τ_m, expresses the speed at which AX tumbles, and is called the rotational correlation time of the molecule. This particular number has been given an ungodly different number of subscripts, but I will use the original notation of Lipari and Szabo and call it τ_m—this means it is the molecular correlation time.

Well, so far so good, but a piddly system that just consists of two atoms isn't really that interesting. Your typical masochistic biomolecular NMR researcher is after molecules (proteins) that have hundreds or thousands of atoms, and even more bond vectors. Now, the picture would still be very simple if these very large molecules were perfectly rigid, but in fact proteins are very flexible and each bond vector is capable of moving on its own, independent of the overall molecular motion. So if we imagine that AX is part of a protein (so that X is bonded to other atoms and A is free) we have to figure out how to represent two motions. This might seem like a major problem, but we still have some hope.

First of all, let's pretend for a moment that the protein doesn't tumble. Well, in that case we could just use the same trick we did earlier, and draw up a correlation function to describe how AX moves inside the protein. Now, we have one difference from the earlier case. Because AX is now part of a much larger system, the correlation function does not decay to zero. Instead, it will decay to some static value, depending on how tightly constrained the bond motion is. If the bond is absolutely rigid, then that value is 1, meaning you have a 100% chance of finding the bond in the same position no matter when you look. If AX is completely flexible, then the ultimate value is 0. So, we'll call this number, which can take values from 0 to 1 and describes the rigidity of the bond, the generalized order parameter and denote it with S². Now, aside from the little wrinkle we just introduced, this is still an exponential decay. So, we need a correlation time, which we will denote with τ_e (an "extra" correlation time). So now we have a second correlation function:

We'll call our first function C_M(t) (for molecular motion) and the second one C_E(t) (for the extra motion). Obviously, these two correlation functions together will completely describe the motion of the bond. But how do we stick them together? If the two motions are correlated (i.e. the probabilities are not independent) this is really quite difficult to do in any general way. But, if the extra motion (the bond moving on its own) really is completely independent of the molecular motion (tumbling), then we can just multiply the functions to get C_T(t), the total correlation function:

where
So now we have an expression to describe the correlation function (i.e. the motion) of the bond, in terms of three parameters (τ_T is just a combination of two others so it doesn't count):

τ_m - the MOLECULAR correlation time
S² - the ORDER PARAMETER
τ_e - the EXTRA correlation time

That's great, but if we could actually measure the correlation function directly we'd hardly need all these parameters, right? We could just look at the correlation functions and compare them. But, it turns out that it is extremely difficult to measure the correlation function of a bond like this directly. Using NMR, however, we can measure something that is related, called the spectral density and denoted by J(ω). What the spectral density is and how NMR can be used to measure it is something for another time (this post is already really freakin' long). The upshot of all of it is, though, that our model of the correlation function C_T(t) translates into this model of the spectral density:
This should be familiar to NMR spectroscopists, though I have omitted scaling factors for clarity. S², τ_m, and τ_T mean exactly the same thing they did in equations 3 and 4 above. The ω represents a frequency, which should give you a little hint about what's going on here: the spectral density is the correlation function represented in terms of frequency, rather than time. This information can be extracted from various relaxation rates measured in NMR spectrometers.

So, now we have a lovely little model of bond motion and a way to get at the parameters experimentally. But it isn't really a model, because we haven't made any sort of assumptions about the nature of the movement involved. We've merely stipulated that there is a global motion, that there is an internal motion, and that these motions are independent. Hence this is called a model-free formalism for describing dynamics. The Lipari-Szabo model-free formalism is extremely powerful because it is very general—the main advantage of this approach is that it can be used to describe almost any motion of a bond. But, the main disadvantage of this approach is that it can be used to describe almost any motion of a bond. Using S², τ_m, and τ_e you can only find out how much a bond is moving, and how fast. The direction and distribution of motions is unknown. So if you want to know, for instance, whether two bonds move in a synchronous or asynchronous manner, the model-free approach can't help you.

Unless you have extremely good luck, you need computer simulations of protein motions to get those kind of answers. All-atom molecular dynamics simulations of proteins have the advantage of explicitly telling you exactly where every bond is moving at every moment. They have some disadvantages too, but the major one of importance for this discussion is that they have historically been very bad at reproducing the dynamics information that comes from NMR, especially when it comes to the motions of side chains. The resolution of that problem will be the subject of my next post.

Lipari, G. and Szabo, A. "Model-free approach to the interpretation of nuclear magnetic resonance relaxation in macromolecules. 1. Theory and range of validity" J. Am. Chem. Soc. 1982: 104(17); pp. 4546-4559.

Read the rest...