Feb 052008
Blogging on Peer-Reviewed ResearchBecause even the simplest case of enzymatic catalysis involves multiple steps (i.e. association, chemistry, and dissociation), semantic issues have obscured the question of what role protein structural dynamics play. For instance, the opening of the lids in Adenylate kinase appears to be rate-limiting, but there is as yet no evidence that dynamics play any role in the phosphotranfer. Thus, one can argue that dynamics are critical (because they are rate-limiting) or irrelevant (because they contribute nothing to the chemical step). Some authors, Arieh Warshel for example, have argued that the latter finding is general—that the highly successful transition-state stabilization model leaves no place for dynamics in chemical catalysis. Some researchers, however, believe that dynamics may play a significant role in chemistry when catalysis proceeds by other routes, especially hydrogen tunneling. In this vein, Matthew Meyer, Diana Tomchick, and Judith Klinman presented evidence in last week’s Proceedings of the National Academy of Sciences that protein dynamics play a significant role in the chemical step of catalysis for the protein soybean lipoxygenase (SLO-1).

Hydrogen tunneling refers to the idea that a hydrogen may in some instances “tunnel through” an energy barrier, going directly from substrate to product without wasting any time or kT in actually surmounting that barrier. A reaction that primarily reflects tunneling should have three properties. The reaction rate should vary only weakly with temperature, because tunneling mostly divorces the rate from kT. Using an Arrhenius plot, this translates to a low calculated energy of activation. Replacing the reactive hydrogen with deuterium should enormously inhibit the reaction—a large kinetic isotope effect (KIE)—because the doubling of mass makes tunneling less likely. Nonetheless, the difference in the energies of activation (ΔEa) for hydrogen and deuterium calculated from an Arrhenius plot should also be small, if tunneling is the main mechanism. Conceivably an enzyme could enhance the rate of such a reaction by using dynamics to promote favorable vibrational modes for the tunneling to occur, or to temporarily adopt a disfavored structure that puts reactive groups at a better distance for tunneling.

SLO-1 exhibits all three features, and so Klinman’s group has used it as model system to try and understand whether and how enzymes promote tunneling reactions. In this case they have performed a series of mutations of isoleucine 553. In addition to existing data on WT and I553A, they analyzed the KIEs, crystal structures, and activation energies of I553L, I553V, and I553G mutants. They find that the protein as a whole is not much distorted by any of these mutations, nor do any of them appear to have significant effects on the binding pocket or substrate dissociation constant—the Kd for each mutant is around 3 μM, while WT is around 10 μM. Yet as the bulk of the mutated residue decreases, the ΔEa increases significantly. In addition, the kinetic isotope effects decrease much more sharply with temperature than for WT. Again, the magnitude of this increased temperature dependence appears to vary inversely with the bulk of the side chain at I553. From these features, and a decline in the magnitude of the pre-exponential factor, the authors argue that dynamics play an important role in encouraging the tunneling reaction.

Well, you know how I love long-range dynamic effects in proteins. But I don’t quite buy it here. An argument based on negatives is never very satisfying in any case, and here we have a lot of questions. For one thing, to say that the structure hasn’t changed significantly seems overly simplistic to me. The lowest-energy structure may not have been seriously deformed, but absent crystal packing and with a little extra kT around, the average structure might be different. Additionally, these structures were obtained in the absence of ligand. Even though the energetics of binding do not appear to differ significantly for these mutants, the structure of the enzyme or ligand may have changed in the bound state. Even if I accept that the existing structures argue for a completely identical binding site in the free state, and this could realistically be debated, that’s no guarantee that the same is true of the bound state.

I’m not denying that the evidence is strongly suggestive, and getting direct data may be difficult given the size of the protein. However, almost all of these side chains are methylated. That means that perdeuteration of the protein in concert with specific side-chain labeling could be especially fruitful in directly observing the dynamics of the pocket during catalysis. It should also in principle be possible to label other side chains in the region to determine how mutations affect the whole area. It is likely that the substrate can also be labeled, possibly with a nucleus or tag that is poorly relaxed. T2 is likely to be a major challenge here, but not necessarily an insurmountable one, especially since the enzyme appears to be folded and active up to around 50 °C. The presence of the iron is certainly a complicating factor; however, the dynamics should be the same whether catalysis is occurring or not, so possibly it could be substituted with a non-paramagnetic metal. It wouldn’t be the easiest project in the world, but it should be doable.

I think these results are intriguing, and I’d love to hear Warshel’s take on them. Nonetheless, I feel I have to reserve judgment at least until dynamic changes in the pocket that seem to correlate with the kinetic observations are directly demonstrated by NMR or some other method.

Meyer, M.P., Tomchick, D.R., Klinman, J.P. (2008). Enzyme structure and dynamics affect hydrogen tunneling: The impact of a remote side chain (I553) in soybean lipoxygenase-1. Proceedings of the National Academy of Sciences, 105(4), 1146-1151. DOI: 10.1073/pnas.0710643105OPEN ACCESS ARTICLE

Nov 082007

Blogging on Peer-Reviewed Research

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

Nov 082007

Blogging on Peer-Reviewed Research

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 S2. 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 CM(t) (for molecular motion) and the second one CE(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 CT(t), the total correlation function:

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
τ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 CT(t) translates into this model of the spectral density:
This should be familiar to NMR spectroscopists, though I have omitted scaling factors for clarity. S2, τ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 S2, τ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.