Methyl Group Dynamics and the Onset of Anharmonicity in Myoglobin

Apr 10, 2008 - The role of methyl groups in the onset of low-temperature anharmonic dynamics in a crystalline protein at low temperature is investigat...
10 downloads 16 Views 508KB Size
5522

J. Phys. Chem. B 2008, 112, 5522-5533

Methyl Group Dynamics and the Onset of Anharmonicity in Myoglobin M. Krishnan,†,‡ V. Kurkal-Siebert,†,§ and Jeremy C. Smith*,†,‡ Interdisciplinary Center for Scientific Computing (IWR), UniVersity of Heidelberg, Im Neuenheimer Feld 368, D-69120, Heidelberg, Germany, and Center for Molecular Biophysics, Oak Ridge National Laboratory, 1 Bethel Valley Road, Oak Ridge, Tennessee ReceiVed: August 17, 2007; In Final Form: February 11, 2008

The role of methyl groups in the onset of low-temperature anharmonic dynamics in a crystalline protein at low temperature is investigated using atomistic molecular dynamics (MD) simulation. Anharmonicity appears at ∼150 K, far below the much-studied solvent-activated dynamical transition at ∼220 K. A significant fraction of methyl groups exhibit nanosecond time scale rotational jump diffusion at 150 K. The splitting and shift in peak position of both the librational band (around 100 cm-1) and the torsional band (around 270-300 cm-1) also differ significantly among methyl groups, depending on the local environment. The simulation results provide no evidence for a correlation between methyl dynamics and solvent exposure, consistent with the hydration-independence of the low-temperature anharmonic dynamics observed in neutron scattering experiments. The calculated proton mean-square fluctuation and methyl NMR order parameters show a systematic nonlinear dependence on the rotational barrier which can be described using model functions. The methyl groups that exhibit many rotational excitations are located near xenon cavities, suggesting that cavities in proteins act as activation centers of anharmonic dynamics. The dynamic heterogeneity and the environmental sensitivity of motional parameters and low-frequency spectral bands of CH3 groups found here suggest that methyl dynamics may be used as a probe to investigate the relation between low-energy structural fluctuations and packing defects in proteins.

1. Introduction The characterization of the internal dynamics of proteins is essential to understanding the mechanisms of their biological functions. Heterogeneous competing protein-protein, proteinsolvent, and intraprotein interactions give rise to a complex energy landscape, leading to a rich and wide spectrum of dynamical relaxational processes.1-4 The atomic motions in proteins can vary from localized vibrational motion, with periods of a few femtoseconds, to collective relaxations on timescales of 10-12 to 1 s in which a significant fraction of the atoms in the protein take part.1 A wide range of experimental techniques has been deployed to elucidate the nature of these dynamic processes.2-16 It is essential to resolve the different types of motions present and investigate their roles in biological function. It is known that proteins exhibit a temperature-dependent dynamical transition around 180-220 K.8,9,17-20 Below this temperature, the dynamics is similar to that of a glassy material, while at temperatures above ∼220 K, protein atoms exhibit liquid-like dynamics. It has been suggested that many proteins function only when the temperature is above the dynamical transition.21-24 The existence of dynamical transitions has been unambiguously shown in many biological systems.2-5,8-25 It was believed that the dynamics of a protein changes from harmonic to anharmonic across the transition.8,9 However, neutron scattering experiments and molecular dynamics simula* To whom correspondence should be addressed. E-mail: smithjc@ ornl.gov. † University of Heidelberg. ‡ Oak Ridge National Laboratory. § Current address: Molecular Modelling, Polymer Physics, BASF, Aktiengesellschaft, 67056 Ludwigshafen, Germany.

tion have shown signatures of anharmonic dynamics even below the ∼220 K dynamical transition for many proteins.15,26-28 The dynamic processes associated with this low-temperature anharmonicity, and how these motions may be related to global dynamical changes at the dynamical transition, are yet to be understood. Furthermore, such investigations can shed light on the relaxation processes that involve conversion between lowlying taxonomic conformational sub-states of the rugged potential energy surface of proteins.29,30 A recent neutron scattering study on hen egg-white lysozyme showed the existence of a low-temperature onset of anharmonicity at around 100 K, the origin of which was suggested to be methyl group rotation.27,28 In neutron scattering experiments, the main contribution to the scattered protein intensity arises from the nonexchangeable hydrogen atoms. A significant fraction of nonexchangeable hydrogens in proteins reside on CH3 groups: 26% in lysozyme, for example. Thus, the CH3 groups contribute significantly to the scattered intensity. Also, it has been suggested that a dominant contribution of the relaxation observed in dry myoglobin neutron scattering is due to methyl dynamics.10 1H NMR relaxation studies have also investigated the reorientational dynamics of C-H bond vectors of methyl groups,6 and 1H NMR experiments on dry lysozyme have shown that 70% of the total proton relaxation is due to methyl dynamics.31 Spectroscopic experiments on a variety of methyl-containing condensed-phase systems have shown the existence of packingsensitive dynamical modes of CH3 groups in the low-frequency domain (below 300 cm-1).32-36 These modes have been used as probes of solid-solid-phase transitions in methyl halides,37,38 hydrocarbon crystals,39,40 polymers41 and for understanding glassy dynamics.42 For proteins, the region below 300 cm-1

10.1021/jp076641z CCC: $40.75 © 2008 American Chemical Society Published on Web 04/10/2008

Onset of Anharmonicity in Myoglobin has a large contribution from collective dynamics in which a number of protein atoms move cooperatively. Some of these collective motions are essential for a protein to function. Unlike highly collective modes, the atomic displacements in CH3-active modes below 300 cm-1 are highly localized on the methyl groups present in the protein side chains. In proteins, the methyl groups tend to occur in the hydrophobic core. Unlike hydrocarbon crystals, in which methyl groups often experience almost identical environments, each methyl group in the hydrophobic region of a protein experiences a different environment. These environmental differences give rise to heterogeneity in the methyl dynamics.43-45 Thus, it may be, in principle, possible to use methyl dynamics as a probe to explore the intrinsic environmental heterogeneity in structure and dynamics existing in protein cores. The variation in the microenvironment around methyl groups gives rise to a wide distribution of their correlation times. NMR experiments have provided information on the site-specific motions of methyl-bearing side chains.46,47 For example, in ref 46, methyl dynamics was used as a probe to explore the changes in the conformational entropy of calmodulin upon complexation with the calmodulin-binding domain of smooth-muscle myosin. Significant changes were found in the dynamics of side chains while the backbone of calmodulin was apparently unaffected by the binding process. The temperature dependence of the fast dynamics of methyl-bearing side chains in this complex has also been investigated.47 Also, the dynamics of methyl groups in proteins has been used as a probe to study the dynamical transition,48 hydrophobic core fluidity,49,50 side-chain flexibility,6,51-53 and protein-ligand binding6,46,54 and to probe the effects of mutations on side-chains in proteins.55 In the present study, we investigate site-specific dynamics of methyl groups in hydrated myoglobin at the temperature corresponding roughly to the onset of anharmonic dynamics using MD simulation and normal-mode analysis (NMA). This investigation provides insights into the anharmonic relaxation processes that involve interconversion between the low-lying conformational sub-states of the potential energy surface of proteins and provides information that can complement sitespecific NMR experiments. 2. Simulation Details Simulations were performed of crystalline carbonmonoxymyoglobin (PDB 1A6G) in a monoclinic unit cell with the following lattice parameters: a ) 63.80 Å, b ) 30.63 Å, c ) 43.42 Å, R ) β ) γ ) 90°. The unit cell consists of two carbonmonoxy-myoglobin molecules. The X-ray crystal structure of a single unit cell, solved at 1.15 Å resolution,56 was generated and used as the starting structure of the MD simulation, performed with a primary box replicated with periodic boundary conditions. Hydration was carried out by initially filling up all of the spaces in the crystal with water, energy minimizing, heating to 300 K, equilibrating for 10 ps with velocity scaling in the NVE ensemble and performing constant pressure (NPT) MD simulation at 300 K for 100 ps. A total of 26 chloride ions were added per unit cell leading to an electrically neutral system. Thus, the simulated unit cell consisted of 8188 atoms. Calculations were performed with the CHARMM program57 and interaction parameters of the CHARMM22 set.58 Water molecules were represented by the TIP3P model.59 Electrostatic interactions were computed using the particle mesh Ewald method60 for which the direct sum cutoff was 1 Å and the

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5523

Figure 1. Time-averaged mean square displacement as a function of temperature. The straight lines are fits to the data for different temperature ranges (solid line, 0 to 150 K; dashed line, 150 to 220 K; dotted line, above 220 K) and are shown as a guide to the eye.

reciprocal space interaction energies were computed on a 64 × 32 × 32 grid using sixth-degree B splines. The system was energy minimized to a root-mean square (rms) force gradient of 10-3 kcal/mol/Å. Subsequently, the system was uniformly heated to 150 K during 15 ps and equilibrated for 100 ps with velocity scaling in the NVE ensemble with P ) 1 bar and T ) 150 K. Equilibration was continued for an additional 200 ps at constant temperature and pressure conditions without velocity rescaling. The temperature and pressure coupling were enforced with the No´se-Hoover algorithm61-63 using the temperature and pressure piston masses of 2000 kcal ps2 and 500 au, respectively. Subsequently, the NPT production runs were performed for 10 ns. The equations of motion were integrated with a time step of 1 fs. Coordinates were written out every 50 fs. Before analysis, all coordinate sets were superposed on a primary-box reference structure so as to remove overall unit cell translation and rotation. Additional production runs of 1 ns duration were also performed sequentially at temperatures starting from 10 to 300 K in steps of 10 K. Normal-mode analysis (NMA) of a single unit cell of myoglobin crystal, replicated with periodic boundary conditions, was performed using the CHARMM program. To do this, the potential energy surface around the energy minimized structure of this crystal was approximated as harmonic and the massweighted Hessian matrix was diagonalized. The eigenvalues of the mass-weighted Hessian matrix provide the vibrational frequencies while the eigenvectors give the direction of the atomic displacements. The eigenvectors were used for the assignment of methyl-specific spectral bands. Additionally, gas-phase simulation of individual methylbearing residues were performed using standard Verlet leapfrog molecular dynamics. The gas-phase dynamics of methyl groups was compared with that in the crystalline environment to understand the role of collisions due to packing interactions in facilitating the rotational motion. 3. Results and Discussion 3.1. Onset of Anharmonicity. Since the aim is to investigate the dynamics of methyl groups at the onset of anharmonicity, the temperature at which the protein begins to exhibit anharmonic dynamics must be identified. For this, the mean square displacement (MSD) of atoms as a function of temperature was calculated (Figure 1). The MSD increases linearly at low temperatures then exhibits two slope changes: one at ∼150 K and the other at ∼220 K. The significant change in gradient at ∼220 K is the solvent-driven dynamical transition as observed in many biological systems.2-12,15,16 Here, we concentrate on

5524 J. Phys. Chem. B, Vol. 112, No. 17, 2008

Krishnan et al.

Figure 2. Schematic representation of a methyl group with atom labels; H1, H2, and H3 are the methyl hydrogens while C1 denotes the methyl carbon. The atom label C2 stands for the atom which is bonded with C1 while C3 represents the heavy atom which forms a bond with C2. The symmetry axis of the methyl group is also shown.

the activation of dynamic processes giving rise to the increase in displacements at 150 K. This low-temperature onset of anharmonicity has been observed in other proteins,15,26 and the dynamical processes associated with this transition were investigated recently in lysozyme using neutron scattering and MD simulations.27,28 In lysozyme, the low-temperature anharmonicity was observed at 100 K and was attributed to the onset of methyl dynamics. It was also demonstrated that the anharmonic dynamics observed at ∼100 K is independent of hydration level,27,28 while the dynamical transition at ∼200-220 K is observed only at hydration levels greater than ∼0.2 g of water/g of protein. 3.2. Methyl Dynamics at 150 K. Carboxy-myoglobin consists of 153 amino-acid residues that are distributed among 8 R helices and the loops that connect these helices. The methylbearing residues are leucine (LEU), iso-leucine (ILE), alanine (ALA), valine (VAL), methionine (MET), and threonine (THR). These residues constitute ∼39% of the entire sequence of myoglobin. Each myoglobin has 94 methyl groups thus amounting to a total of 188 methyl groups in the crystalline unit cell. In the following analysis, we restrict ourselves to addressing a few specific questions: Are methyl groups dynamically active far below the ∼220 K dynamical transition of proteins (particularly near low-temperature anharmonicity)? If yes, what is the nature of their dynamics? Do methyl groups exhibit anharmonic and heterogeneous dynamics at such low temperatures? If their dynamics is anharmonic and heterogeneous, what are the microscopic reasons for such a behavior? Also, what can we learn about the protein interior from the delocalized and heterogeneous dynamics of the methyl groups? In order to answer these questions, we have carried out investigations of site-specific methyl dynamics in the myoglobin crystal at 150 K. Also, we have characterized the role of the environment of each methyl group in its dynamics. In the following discussion, we present the answers to the above questions sequentially. 3.2.1. Dynamical ActiVity of Methyl Groups Far Below Dynamical Transition. In the following discussion, we assign labels to the atoms of a methyl group (Figure 2). Label C1 represents the methyl carbon while the C1-C2 and C2-C3 are the bonds medial to the methyl group. The molecular unit shown in Figure 2 forms a fragment of a protein side chain. For instance, C2 and C3 for the side chains of ILE, VAL, and THR

Figure 3. 10 ns probability density map of a tagged proton of a methyl group around its carbon is shown for representative CH3 groups in crystalline myoglobin at 150 K. The gray sphere represents the methyl carbon. As the symmetry axes of the methyl groups are oriented randomly in the protein structure, comparison of the proton density maps of different methyl groups requires the use of a fixed frame of reference for all CH3 groups. This is accomplished using the following procedure: for each CH3 group, we perform two consecutive rotations such that the normal to the plane formed by the atoms C1, C2, and C2 (see Figure 2) orients along the x axis, and the bond vector connecting atom C1 and C2 is along z axis.

correspond to Cβ and CR, and for ALA, these correspond to CR and N while for MET they correspond to S and Cγ, respectively. NMR experiments aimed at understanding protein side-chain motions study the relaxation data of the C-H bond vectors of methyl groups. These experiments have shown that methyl relaxations are due mainly to either rotation of C-H vectors about the symmetry axis or the reorientational dynamics of the symmetry axis.6 We first investigate the rotation of C-H vectors about the symmetry axis. This type of rotational methyl-group dynamics is sensitive to the local packing density and thus has been used earlier to characterize the solid-solid transition in methyl halides37,38 and the high pressure solid phases of n-alkanes.39,40 The direct way to investigate this dynamics is by calculating the probability density of a tagged methyl hydrogen around the methyl carbon. In Figure 3, the calculated probability density map for representative methyl groups is shown. Three lobes are observed around C1 for LEU107δ1 and ILE112δ1, while THR51γ2 and ILE30δ1 exhibit two lobes, and LEU61δ2 and LEU29δ1 show a single lobe. The presence of nonzero probabilities at more than one site is due to mobility of the tagged hydrogens of the methyl groups. Moreover, the absence of density between the lobes indicates that these tagged hydrogens exhibit jump-like motions from one site to the other. At 150 K, it is found that 11% and 17% of the methyl groups in myoglobin exhibit trimodal and bimodal distributions, respectively, over the 10 ns simulation period. These delocalization processes of methyl hydrogens can

Onset of Anharmonicity in Myoglobin

Figure 4. (a) Time evolution of the angle, θ, between a C-H bond vector of a methyl group at time t and at time zero is shown for selected methyl groups of myoglobin at 150 K that exhibit different dynamical behavior. (b) The rotational trajectories of representative mobile methyl groups of myoglobin crystal are shown here.

contribute significantly to the residual entropy of proteins thereby playing a key role in determining the thermodynamic properties of proteins at low temperatures.47,64,65 3.2.2. Nature of Methyl Dynamics at Low Temperature. The probability density maps in Figure 3 indicate that the hydrogens of some methyl groups are delocalized and exhibit jump-like motion. However, Figure 3 does not convey the time scales associated with these site-specific dynamical events. For this, the time evolution of the C-H bond vector connecting C1 and the tagged proton of a methyl group must be calculated. In Figure 4a is shown the time dependence of the angle (θ) subtended by the component of C-H bond vector that is perpendicular to C1-C2 vector at time t with the same vector at time zero for the same methyl groups shown in the proton density map (see Figure 3). Although the time evolution of methyl dihedral angle would provide the same information, the initial value of θ is zero while it can be nonzero for the dihedral angle. The following observations are evident from this figure: (i) The C-H vectors of a significant fraction of methyl groups undergo sharp 120° transitions from one orientation to another. (ii) The frequency of occurrence of these sharp transitions is very different among the methyl groups. Indeed, some methyl groups did not show any transitions during the entire 10 ns trajectory. The sharp orientational change of C-H bond vectors is due to the anharmonic jump rotational transition of methyl groups between potential energy minima. Furthermore, the wide variation in the time required for the C-H vectors to reorient implies heterogeneity in their dynamics. The jump-like motions

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5525 of hydrogens will lead to quasi-elastic neutron scattering of the type that has been observed at 150 K in proteins.15,19,27,28 In Figure 4b is shown the time evolution of θ for a few mobile methyl groups from each methyl-bearing residue type in myoglobin. It is evident from this figure that the rates of rotational transitions in ILE methyl groups (particularly, ILE112δ1, ILE28δ1 (2), ILE75δ1 (2), ILE111δ1, and ILE107δ1) are higher than those for other residue types. The CH3 groups in LEU residues (LEU137δ2, LEU69δ2 (2), LEU104δ2, LEU149δ2, LEU72δ2) also exhibit rotational transitions but at lower rate. Correspondingly, 16 out of 36 ILE methyl groups exhibit an 120° orientational change at least once during the course of the simulation compared with 18 out of 72 for LEU. A total of 10 out of 32 of VAL, 3 out of 34 of ALA, and 2 out of 10 of THR residues exhibit 120° transitions while all of the MET methyl groups undergo many rotations at the onset of low-temperature anharmonicity. In the CHARMM force field, methyl torsional motion is strongly influenced by the “intrinsic” dihedral potential, which is of the form Eφ ) Kφ(1 + cos 3φ), where φ is the dihedral angle and Kφ is the force constant, although contributions can also arise from nonbonded and other bonded terms.66,67 The methyl groups in ILEδ and MET residues have relatively lower intrinsic potential energy dihedral barriers than other methyl groups. The value of the force constant associated with the torsional motion of δ-methyl groups in ILE is 0.16 kcal/mol (Kφ ) 0.2 kcal/mol for MET) thus summing to an intrinsic barrier of 2.88 kcal/mol (1.68 kcal/mol for MET) while the other methyl groups in the protein have intrinsic barriers of 3.6 kcal/ mol. Thus, the methyl groups in ILEδ and MET residues rotate more around their symmetry axis than other methyl groups. The wide variation in the frequency of rotational transitions observed among ILE δ-methyl groups indicates that crystal interactions play a key role in modulating the effective rotational barriers. 3.2.3. Methyl Rotational Time Correlation Function. The methyl dynamics can be further understood using rotational time correlation functions of the C-H bond vectors. Let rˆCH(t + τ) denote the unit vector pointing along the C-H vector of a methyl group at time t + τ. A useful time correlation function of this vector is calculated as follows:

CNMR(t) ) 〈P2[rˆCH(0) .rˆCH(t)]〉

(1)

where P2[x] is the second-order Legendre polynomial and is given as

1 P2[x] ) [3x2 - 1] 2

(2)

The subscript NMR indicates that this correlation function can in principle be used to interpret NMR relaxation experiments. The bracket symbols 〈 〉 indicate the ensemble average over all three C-H vectors of a methyl group and over all time origins. Dynamical NMR relaxation experiments can be analyzed using a “model-free” approach which relates CNMR(t) to a generalized order parameter, S2 and correlation times associated with backbone or side-chain motions.68 S2 is defined as

S2 ) lim CNMR(t) tf∞

(3)

and the relaxation time, τS, is given as

τS )

1 1 - S2

∫0∞ (CNMR(t) - S2) dt

(4)

5526 J. Phys. Chem. B, Vol. 112, No. 17, 2008

Krishnan et al.

Figure 6. Power spectrum obtained by the Fourier transformation of the autocorrelation function of the force acting on methyl group hydrogen is shown for representative methyl groups.

Figure 5. (a) Rotational time correlation function, CNMR(t), of a C-H bond vector of a methyl group is shown for the methyl groups from Figure 4. (b) The behavior of CNMR(t) for all methyl groups from one of the myoglobin molecules in the crystalline unit cell, whose C-H bond vectors exhibit an 120° orientational change at least once during the course of the 10 ns simulation.

In Figure 5b is shown CNMR(t) for all methyl groups that undergo an 120° orientational change at least once in the 10 ns. Variation in the relaxation behavior of CNMR(t) suggests that the motional parameters vary widely among the methyl groups, thus again indicating heterogeneity in their dynamics. 3.2.4. Site-Specific Spectral Analysis of Methyl Relaxations. A further useful way of characterizing the influence of sitespecific local molecular packing and atomic fluctuations on methyl dynamics is to perform spectral analysis of individual methyl groups. Here, we use NMA, MD-derived time correlation functions and the theoretical neutron scattering methods to perform the spectral analysis. NMA describes only the harmonic dynamics while the MD-derived spectra also include anharmonic contributions. We first calculate the forces acting on methyl group hydrogens due to their interaction with the environment. Let Fi,j(t) be the force that is felt by the jth hydrogen atom of the ith methyl group at time t. The force-force time correlation function, CF,i(t), of the ith methyl group is defined as follows,

〈∑ 〈∑ 3

Thus, S2 describes the amplitude of motion while τS gives the relaxation rate associated with the motion. These motional parameters are extensively used to understand biological processes such as protein-ligand binding46 and thermodynamic properties of proteins and crystalline amino acids.69 Figure 5a shows CNMR(t) for the same representative methyl groups for which the rotational trajectories were shown in Figure 4a. Initially, all of these functions relax to a value close to 0.9 in less than 1 ps, because of fluctuations of the C-H bond vector around its equilibrium orientation. Subsequently, the behavior of CNMR(t) differs among methyl groups. For instance, the CNMR(t) of ILE107δ1 and ILE112δ1 decay to a value close to zero at around 2.5 ns while those of THR51γ2 and ILE30δ1 relax to only 0.6 over the same time period. ILE61δ2 and LEU29δ1 CNMR(t) remain constant at 0.9 at all times. The relaxation of this function to zero for ILE107δ1 and ILE112δ1 indicates that these methyl groups are mobile, leading to decorrelation of CNMR(t) while the value close to 1 for ILE61δ2 and LEU29δ1 implies that these methyl groups are not mobile. A value between 0.3 and 0.7 means that the methyl groups concerned are sterically hindered and undergo rotational transitions infrequently. The methyl dynamics at 150 K is clearly not fully explored over the 10 ns time scale of the simulation. For a methyl group having ideal tetrahedral geometry, CNMR(t) decays to 1/9 (i.e., the ideal value of a converged S2 is 1/9), unless it is completely frozen.70,71 The present analysis serves mainly to pinpoint the short time (nanosecond and below) dynamics incipient at the glass transition.

CF,i(t) )

j)1 3

j)1

Fi,j (t)‚Fi,j (0)

〉 〉

(5)

Fi,j (0)‚Fi,j (0)

where the angle brackets denote averaging over all time origins. The function CF,i(t) exhibits a fast decay (t < 5 ps) together with characteristic long-time fluctuations (not shown). The associated power spectrum, CF,i(ω), is calculated by Fourier transforming CF,i(t):

CF,i (ω) )

∫0∞ CF,i (t) cos(ωt) dt

(6)

Since low-frequency modes are highly sensitive to the environmental effects, we consider here only spectral features below 300 cm-1. In Figure 6 is shown the behavior of CF,i(ω) for a few representative methyl groups that are dynamically different. Two spectral bands are found that are specific to methyl vibrations: one between 50 and 125 cm-1 and the other around 290 cm-1. The spectral feature around 290 cm-1 has earlier been attributed to methyl rotational dynamics.32,39 The frequency and the line width of the spectral features differ among the methyl groups because of differences in the environment. For example, the lowest frequency peak is centered around 70 cm-1 for LEU61δ2 and LEU29δ1 while for the other residues it is observed around 90 cm-1. Further, the density of states between the two peaks observed in the spectral region between

Onset of Anharmonicity in Myoglobin

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5527

50 and 125 cm-1 is close to zero for LEU61δ2 and LEU29δ1 while significant overlap of these peaks can be seen for the other methyl groups. Similarly, the spectral feature observed around 290 cm-1 is also different for these two methyl groups than the rest: sharp peaks are present for LEU61δ2 and LEU29δ1 while the other methyl groups exhibit broad peaks with spectral splittings. Examination of the rotational trajectories and the CNMR(t) show that LEU61δ2 and LEU29δ1 do not exhibit rotational excitations during the 10 ns MD trajectory, while the other CH3 groups undergo hindered rotation. Experiments have shown similar spectral bands in these spectral regions.32,33,37 For instance, incoherent inelastic neutron scattering experiments have observed torsional excitations of methyl groups around 270 cm-1 32 while Raman spectra of crystalline methyl halides have shown CH3 librational bands around 50-120 cm-1.37 Highresolution inelastic neutron scattering spectroscopy experiments combined with NMA for probing vibrations in the enzyme Staphylococcal nuclease have observed a sharp peak at 235 cm-1 (269 cm-1 in NMA) which was attributed to methyl torsion.33 3.2.5. Assignment of Spectral Bands. In order to understand the dynamic processes associated with the low-frequency methyl-specific spectral bands, we employ two different approaches: (1) spectral analysis of the torque acting on the C-H vectors of the CH3 groups and (2) NMA. In the former approach, we define ⊥ ri,j ) ri,j - nˆ (ri,j‚nˆ )

(7)

⊥ Fi,j ) Fi,j - nˆ (Fi,j‚nˆ )

(8)

where nˆ is the unit normal to the plane formed by hydrogen atoms of a methyl group, ri,j is the C-H bond vector (connecting the methyl carbon and ith hydrogen), and Fi,j is the force acting on the ith hydrogen of the jth methyl group. The component of the torque acting on a methyl proton is calculated as ⊥ ⊥ × Fi,j τi,j ) ri,j

(9)

Figure 7. (a) Power spectrum obtained by the Fourier transformation of the autocorrelation function of individual (solid line) and total (dashed line) torque acting on the methyl protons shown for the same methyl groups reported in Figure 6. (b) The atomic displacement vectors associated with the (A) librational motion at 114 cm-1 and (B) rotational motion at 240 cm-1 of a methyl group in ILE107δ1 residue.

3

Tj )

τi,j ∑ i)1

(10)

The autocorrelation functions of these torques are given as

〈 〈∑ 3

Cjτ(t)

)

τi,j(t)‚τi,j(0) ∑ i)1 3

(11)

τi,j(0)‚τi,j(0)

i)1

CjT(t) )

〉 〉

〈Tj(t)‚Tj(0)〉

〈Tj(0)‚Tj(0)〉

(12)

where the angle brackets represent averaging over all time origins. The CjT(t) characterizes the collective rotation of methyl protons (methyl rotation) while Cjτ(t) provides insight into the uncorrelated motion of individual methyl protons (for example, methyl libration). The power spectra obtained from these correlation functions are shown for selected methyl groups in Figure 7 for the same frequency range as that of Figure 6. Interestingly, the spectra obtained from Cjτ(t) reproduce the low-frequency spectral band (between 50 and 125 cm-1) observed in the spectra calculated from the force autocorrelation

function. Thus, it is evident from this analysis that the spectral band observed between 50 and 125 cm-1 arises from a librational-type motion of the methyl hydrogens that changes the instantaneous orientation of the methyl symmetry axis. The spectra obtained from CjT(t) (also shown in Figure 7) show spectral bands between 230 and 320 cm-1 indicating that methyl rotation is excited in this frequency range. This interpretation is confirmed by inspection of the eigenvectors obtained from the NMA approach. In Figure 7 (bottom panel), displacement eigenvectors for the methyl hydrogens in ILE107δ1 are shown. The nature of the atomic displacements shown for a mode at 240 cm-1 (similar methyl active modes were observed between 230 and 310 cm-1 where methyl groups exhibit large-amplitude harmonic displacements) indicates that the spectral peak observed around 290 cm-1 can be safely assigned to methyl rotation while the librational motion of methyl protons is evident from the mode at 114 cm-1 shown in Figure 7. 3.2.6. Comparison with Theoretical Neutron Scattering Spectra. The spectra obtained from the force-force auto correlation function can be compared with the incoherent neutron scattering spectra calculated from the MD trajectory. The incoherent dynamic structure factor is defined as follows:

Sinc(q b, ω) )

1 2π

∫-∞∞ Iinc(qb, t) e-iωt dt

(13)

5528 J. Phys. Chem. B, Vol. 112, No. 17, 2008

Krishnan et al.

Figure 8. Theoretical neutron scattering spectra (bottom profiles) of selected methyl groups compared with corresponding spectra obtained from force auto correlation function (top profiles). The theoretical neutron scattering spectra were calculated using “nMoldyn” software72

Here, pq b and pω are the momentum and energy transfer, respectively, and Iinc(q b, t) denotes the intermediate scattering function given by,

Iinc(q b, t) )

1

2 bi(inc) 〈e-iqb‚br (0) eiqb‚br (t)〉 ∑ N i i

i

(14)

where bi(inc) is the incoherent neutron scattering length of atom i while b ri(t) is the position vector of atom i at time t. Since the incoherent neutron scattering length of the hydrogen atom is much higher than that of other atoms, in proteins, the dynamic structure factor obtained from neutron scattering experiments is dominated by contributions from the nonexchangeable hydrogens. In Figure 8, Sinc(q b, ω) together with the corresponding forceforce spectrum for selected methyl groups of myoglobin is shown. The calculated neutron scattering spectra reproduce the spectral features obtained from the force-force correlation function. The advantage of using force-force spectra in interpreting neutron scattering is that one can isolate the contributions of specific interactions. For example, by obtaining the time series of the force acting on a methyl group due to the solvent molecules alone and performing a power spectral analysis, the effect of the solvent on the dynamics of methyl groups can be determined. Here, these spectra were found to be featureless and did not have methyl-specific spectral bands below 300 cm-1 thus indicating weak coupling between the solvent and the methyl dynamics (see next section). 3.2.7. Effect of SolVent on Methyl Dynamics. Does solvent exposure of a methyl group play a role in determining its dynamics? In order to answer this question, we first estimate the degree of solvent exposure by calculating the radial distribution function, gC1-WO(r), of the carbon atom of any given methyl group with the oxygen atoms (denoted as WO) of the water molecules. In Figure 9, spherically averaged gC1-WO(r) around selected methyls from different locations of the protein at 150 K are presented. The distribution functions for ILE112δ1, LEU2δ1, LEU137δ2, and ILE112γ2 exhibit well-defined first coordination shells followed by additional sub-peaks at larger distances, indicating that these methyl-bearing side chains are exposed to the solvent. The absence of intense peaks in the gC1-WO(r) of ILE107δ1, LEU32δ1, LEU104δ2(2), and ILE107γ2 (2) indicates a lack of a first solvation shell around these methyl groups, indicating that these side chains are buried. The integrated coordination number (CN) of the water molecules around the carbon atoms of the methyl groups, obtained by integrating gC1-WO(r), is provided in Table 1. The coordination number at 7 Å varies from 1 to 32, and provides

Figure 9. Radial distribution of water oxygen atoms around methyl carbons is shown for representative methyl groups of crystalline myoglobin. These functions were obtained by averaging over 10 ns MD trajectory. Here, MB2 represents the second myoglobin molecule in the crystalline unit cell.

TABLE 1: Fraction of Methyl Groups (Given as a Percentage) That Are Surrounded by CN Number of Water Molecules within a Given Distance, r, from the Methyl Carbon CN

P% (r ) 4.0 Å)

P% (r ) 5.0 Å)

P% (r ) 6.0 Å)

P% (r ) 7.0 Å)

0-4 4-8 8-12 12-16 16-20 20-24 24-28 28-32 32-36

84.0 16.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

48.9 26.6 22.3 2.1 0.0 0.0 0.0 0.0 0.0

25.5 28.2 18.1 17.6 9.0 1.6 0.0 0.0 0.0

13.3 13.8 17.0 16.0 13.3 10.6 9.0 4.8 2.1

a good measure of the relative burial of methyl groups: those with a coordination number less than 6 at 7 Å can be considered as buried while CN > 20 at 7 Å are exposed to solvent. We now examine the correlation between the solvent exposure and the dynamics of CH3 groups by comparing the amplitudes of atomic displacements of the hydrogen atoms and the solvent coordination number. The amplitudes of atomic displacements were calculated using both the NMA approach and the MD trajectory. In the NMA approach, the magnitude of fluctuations of methyl protons are related to the eigenvectors and the vibrational frequencies through the following expression: 3N

〈u 〉nma ) KBT 2

3

∑ ∑ i)1 j)1

b e ji‚e bji mjωi2

(15)

where b eji represents three components corresponding to the jth atom of the ith eigenvector, ωi is the vibrational frequency of the ith normal mode, T is the temperature, KB is the Boltzmann constant, mj is the mass of the jth atom, and N denotes the total number of atoms in the system. The atomic displacements calculated from NMA provide only the harmonic fluctuations. To evaluate the anharmonic part of methyl dynamics, 〈u2〉 of hydrogen atoms of methyl groups were calculated from the MD trajectory. In Figure 10, the 〈u2〉 of methyl hydrogens calculated from NMA and MD plotted against CN are shown. It is evident from the MD 〈u2〉 that many methyl groups exhibit small amplitude motions with 〈u2〉 around 0.79 Å2, while a significant fraction of CH3 groups show largeamplitude displacements with an average 〈u2〉 greater than 1.5 Å2. This large-amplitude anharmonic excitation of methyl groups

Onset of Anharmonicity in Myoglobin

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5529

Figure 10. (a) Mean-square methyl hydrogen atom fluctuation calculated from NMA approach plotted versus solvent coordination number (CN). (b) Average mean squared fluctuations of each methyl group in myoglobin, evaluated from the MD trajectory, is plotted against CN.

is seen both at low and high values of CN, again indicating a negligible role of solvent exposure in methyl dynamics. To further examine the solvent effect, the solvent density distribution was calculated around each CH3 group and the NMR relaxation functions were compared. In Figure 11a, the arrangement of water molecules around selected “mobile” and “rigid” methyl groups is shown. For ILE107δ1, LEU104δ2(2), LEU32δ1, and ILE107γ2(2), there are less than 6 water molecules within a distance of 7 Å, while for LEU137δ2, ILE112δ1, LEU2δ1, and ILE112γ2, there are more than 25 water molecules within this distance. For comparison, in Figure 11b, CNMR(t) for the mobile methyl groups is shown. Irrespective of the degree of solvent exposure, the CNMR(t) relaxes to close to zero for the mobile methyl groups, while for the rigid methyl groups, the CNMR(t) fluctuates around 0.9 (not shown). The above results indicate little or weak coupling between the solvent and the methyl dynamics thus providing a possible explanation for the hydration independence of activation of low-temperature anharmonic dynamics observed in recent neutron scattering experiments.27,28 3.2.8. Rotational Barrier and Motional Parameters. We now perform a quantitative examination of the effect on the rotational barrier of the microenvironment of methyl groups. This is done by comparing the rotational barriers with the dynamical properties (〈u2〉 and NMR order parameter) of the methyl groups. The rotational barriers were calculated from the rotational trajectories of the methyl groups (see Figure 4). The probability distribution functions, P(θ), of θ were obtained from the rotational trajectories, and the barrier heights (∆E) were calculated from the potential of mean force (PMF) profiles that were obtained from P(θ) using the following equation:

U (θ) ) -KBT log P(θ)

(16)

Since we study methyl dynamics at a low temperature, the methyl rotational transitions are infrequent on the 10 ns simulation time scale. The low energy conformations dominate the Boltzmann distribution and there is incomplete sampling of the conformations close to the barrier. This incomplete sampling prevents the determination of the exact PMF profiles

Figure 11. (a) Atomic density map of water oxygen around a given methyl group for selected “mobile” (ILE107δ1, LEU137δ2, LEU104δ2(2), ILE112δ1) and “rigid” (LEU32δ1, LEU2δ1, ILE107γ2(2), ILE112γ2) methyl groups and (b) their corresponding rotational time correlation functions, CNMR(t). The probability maps are calculated in the same fixed frame of reference as in Figure 3.

and the barrier heights from brute force MD simulations. Thus, in our calculations, the effective barrier height is taken as the highest statistically significant energy (difference in the PMF between the highest probable methyl conformation and the lowest probable conformation with a least nonzero probability) sampled during methyl rotational dynamics. In Figure 12 is shown the 〈u2〉 of the methyl protons plotted against the rotational barrier. Most of the CH3 groups have barrier heights around 2.8 kcal/mol while the rest of the methyl groups have rotational barriers between 2.4 and 2.8 kcal/mol. Recent NMR relaxation experiments have determined the barrier to be 2.8 kcal/mol and suggested that a typical methyl rotational barrier in proteins ranges from 2.5 to 3.0 kcal/mol.73 It is evident from Figure 12 that the 〈u2〉 increases systematically with decreasing barrier height. This result suggests that a relatively subtle barrier reduction gives rise to large amplitude motions of methyl protons. The average 〈u2〉 calculated around ∆E ) 2.8 kcal/mol is 0.8 Å2 while that evaluated around ∆E ) 2.45 kcal/mol is 1.94 Å2. Thus, 10% reduction in the barrier height gives rise to more than a 2-fold increase in the 〈u2〉. This barrier dependence of the 〈u2〉 can be fitted to the following functional form:

〈u2〉 ) a0 + a1 ea2(∆E- 2.81-a3)

(17)

The solid line in Figure 12a shows the best fit to the data using eq 17 and the best fitting parameters are a0 ) 2.10 Å2, a1 ) -4.66 Å2, a2 ) 7.10 (kcal/mol)-1 and a3 ) 0.15 kcal/mol. The

5530 J. Phys. Chem. B, Vol. 112, No. 17, 2008

Krishnan et al.

Figure 12. (a) Atomic fluctuations (〈u2〉) of methyl protons versus rotational barrier (∆E). (b) NMR order parameter (S2 defined in eq 3) of methyl groups versus ∆E. The data corresponding to representative mobile methyl groups shown in Figures 4a and 5a are also shown here with different symbols.

five data points (corresponding to LEU86δ1, LEU86δ2, ILE142δ1, LEU2δ1(2), LEU2δ2(2)) that are scattered away from the solid lines were not considered while fitting the data to the model function. We have also investigated the dependence of NMR order parameter (S2 defined in eq 3) on the barrier height. In Figure 12b is shown the NMR order parameter versus the rotational barrier. Here, again, a systematic dependence is observed of S2 on ∆E. The values of S2 for CH3 groups with ∆E around 2.8 kcal/mol are around 0.9 suggesting that these methyl groups are motionally restricted while a decrease in ∆E leads to a significant decrease in S2: A reduction in barrier height of ∼0.3 kcal/mol decreases S2 to close to zero. Figure 12b indicates that the NMR order parameter (thus methyl dynamics) is indeed highly sensitive to ∆E. This dependence can be fitted to a sigmoidal type function defined as follows:

S2 ) a0 +

1 1+e

-a1(∆E - a2)

(18)

The solid line in Figure 12b shows the best fit to the data using eq 18, with parameters a0 ) 0.01, a1 ) 24.34 (kcal/mol)-1 and a2 ) 2.68 kcal/mol. Equation 18 has a very simple physical meaning: the higher the rotational barrier, the more restricted the motion of the methyl group. The parameter a2 is the value of the barrier height at which the value of the order parameter is equal to 0.5 (i.e., S2 ) 0.5; this would mean that the methyl group is partially mobile and partially restricted). Methyl groups that have barrier heights larger than a2 are considered to be more sterically hindered while the dynamics of methyl groups with ∆E less than a2 are highly mobile because of less steric hindrance. The rate of decay of S2 with respect to ∆E is determined by the parameter a1. Equations 17 and 18 relate methyl motional parameters to rotational barriers. Although these relations can be used to estimate methyl barriers from the 〈u2〉 and S2, the latter two quantities will be time and temperature dependent; that is, increasing the temperature or the MD trajectory length will lead

Figure 13. (a,b) Time evolution of CNMR(t) (left panel) and the rotational trajectories (right panel) of “catalyzed” methyl groups are compared with their gas phase (G) dynamics.

to changes in 〈u2〉 and S2 for those methyl groups that are not converged within 10 ns. Similar incomplete convergence was reported in earlier atomistic MD simulations on other proteins at room temperature.43,44 However, a systematic dependence of the motional parameters on ∆E is evident (eqs 18 and 17). 3.2.9. Barrier Reduction by Crystal EnVironment. It has been shown that in crystals of small peptides the microenvironment around methyl groups facilitate rotational dynamics by reducing the effective barrier to rotation.74,75 To understand this “catalytic” effect, the dynamics of methyl groups in myoglobin crystal is compared with the results obtained from simulation of individual methyl-bearing residues in the gas phase, and the results are shown in Figure 13. The C-H bond vectors of the LEUδ2 methyl group do not exhibit rotational excitations in the gas phase (indicated as G) while LEU137δ2, LEU69δ2(2), and LEU104δ2 undergo many rotational excitations in 10 ns. Correspondingly, for these methyl groups, the CNMR(t) relax faster than LEUδ2 (G). Thus, for these groups, there is significant barrier reduction in the crystal. Similar effects are also seen for a few CH3 groups in valine residues, as shown in Figure 13b. Barrier reduction was also observed for CH3 groups in some isoleucine residues. However, rotational excitations of the C-H bond vectors (in ILEδ1) were also observed in the gas phase. In total, it was found that around 10% of the methyl groups in myoglobin crystal exhibit barrier reduction relative to the gas phase. 3.2.10. Xenon CaVities Act as ActiVation Centers of Anharmonic Dynamics. A recent investigation of the role of crystal environment in catalyzing methyl rotation has suggested that van der Waals interactions due to the nonbonded environment of the side chain can reduce the barrier to methyl rotation.75 The essential role of van der Waals interactions in such local motions41,75,76 suggests that variation of protein packing density around methyl groups could play a dominant role in the dynamic

Onset of Anharmonicity in Myoglobin

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5531

Figure 14. Distribution of distances between carbon atoms of methyl groups and closest xenon cavities in the crystal structure of carboxymyoglobin is shown as a histogram. The fraction of methyl groups (given in percentage) that are located within a given distance, r, away from the closest xenon cavities is also shown (solid line with a circle symbol).

TABLE 2: Distances between the Carbon Atoms of Methyl Groups That Exhibit Many Rotational Excitations and the Nearest Xenon Cavities Where XEn Represents the nth Xenon Cavity in Myoglobin index

methyl group

nearest xenon cavity

distance [Å]

1 2 3 4 5 6 7

LEU137 (Cδ2) LEU69 (Cδ2) LEU104 (Cδ2) ILE111 (Cδ1) ILE107 (Cδ1) ILE75 (Cδ1) ILE112 (Cδ1)

XE3 XE4 XE1 XE2 XE4 XE3 XE4

0.0 0.0 0.0 0.0 0.0 3.1 6.1

heterogeneity observed here. The local packing density influences the compressibility and other thermodynamic and kinetic properties of proteins.77 Most proteins have “packing defects”, that is, cavities, with volumes ranging from 30 to 100 Å3.78,79 Acting as source of free volume, these packing defects reduce the thermodynamic stability of proteins79 and may act as channels for the transport of ligands to and from the active site.78,79 For metmyoglobin, X-ray crystallographic studies on the protein equilibrated with xenon gas led to the observation of four xenon-binding sites in interior spaces of this protein, indicating the presence of cavities large enough to accommodate xenon atoms.78,79 These xenon cavities form parts of pathways for migration of ligands to and from the heme group. To understand whether the xenon cavities play a role in the activation of anharmonic dynamics in proteins, we examined at the location of the methyl groups that exhibited many rotational excitations during the course of the 150 K simulation. In Figure 14 is shown the distribution of distances between carbon atoms of methyl groups and closest xenon cavities in the crystal structure of carboxy-myoglobin. It can be seen that ∼35% of methyl groups in myoglobin is located within 4 Å of the xenon cavities. The methyl groups of LEU and ILE residues that exhibit many rotational excitations and faster relaxation relative to the gas phase are located near xenon cavities (see Table 2). In Figure 15 are shown these mobile methyl groups together with xenon cavities. The methyl groups of LEU104δ2, LEU69δ2, LEU137δ2, ILE111δ1, and ILE107δ1 are part of xenon cavities while ILE75δ1 was located at 3.1 Å away from the xenon cavity (XE3). The δ-methyl group in ILE112 is located at 6.1 Å away from XE4. These distances are much less than half of the radius of gyration of the protein (Rg ) 16.29 Å calculated from the last 2.5 ns of the MD trajectory). The above result suggests

Figure 15. (a) Surface representation of cavity forming atoms shown together with a cartoon representation of myoglobin. (b) Locations of methyl groups that exhibit many rotational excitations shown together with the xenon cavities in myoglobin. The spheres represent the positions of the methyl carbons while the xenon cavities are shown by surfaces formed by atoms that form these cavities. These cavities were determined using CASTp software with a probe radius of 1.4 Å.80 The distances between the methyl carbons and the nearest cavity forming atoms are given in Table 2. The atom labels of the methyl carbons are also given in Table 2.

that methyl groups situated close to xenon cavities are activated at the onset of low-temperature anharmonic dynamics due to the availability of free volumes near packing defects. Thus, internal cavities may play a key role in triggering anharmonic dynamics in proteins. 4. Conclusions The present all-atom molecular dynamics simulations of crystalline myoglobin investigate the role of methyl dynamics in low-temperature anharmonicity of proteins. The simulated system begins to exhibit anharmonic dynamics at around 150 K. A significant fraction of the nonexchangeable hydrogens in the protein (i.e., those dominating neutron scattering) are present in methyl groups, and these may play a key role in the onset of low-temperature anharmonic dynamics (150 K). The calculated atomic proton probability density maps and the rotational trajectories of a significant fraction of methyl groups exhibit delocalized, jump-like motion of methyl protons at 150 K within the 10 ns time scale of the simulation. This jump diffusion will produce quasi-elastic neutron scattering, as has been observed in proteins at low temperatures. The calculated rotational time correlation functions of C-H bond vectors of individual methyl groups suggest a wide distribution of rotational relaxation rates among methyl groups, indicating heterogeneity in their dynamics.

5532 J. Phys. Chem. B, Vol. 112, No. 17, 2008 Using spectral analysis of individual methyl groups, methylspecific dynamical modes are found in the low-frequency region (below 300 cm-1) of the spectrum. The power spectra obtained by Fourier transformation of the autocorrelation function of forces acting on the methyl protons in the MD match well with the calculated incoherent neutron scattering spectra. Two spectral bands are observed: one around 290 cm-1 and the other between 50 and 125 cm-1. NMA and analysis of the autocorrelation functions of the torques acting on the C-H bond vectors of methyl groups show that the band between 50 and 125 cm-1 arises from librational motions of the methyl groups while the feature present around 290 cm-1 is the methyl rotational mode. The spectral splittings, the peak positions, and the line widths of these low-frequency methyl-active modes differ among the CH3 groups in proteins depending on the local packing densities. The role of solvent exposure in methyl dynamics was also investigated by comparing solvent distribution functions and the rotational trajectories of methyl groups. The methyl dynamics is uncorrelated with the degree of solvent exposure, consistent with the hydration independence of the activation of low-temperature anharmonic dynamics observed in recent simulations and neutron scattering experiments.19,27,28 This contrasts with the strong hydration dependence observed for the higher-temperature (∼220 K) transition.17-20,24,25 Calculated methyl motional parameters, that is, the meansquare fluctuation and the NMR order parameter, show a systematic nonlinear dependence on the rotational barrier. The methyl groups that exhibit many rotational excitations (at the onset of low-temperature anharmonicity) during the course of the simulation are found to be located near xenon cavities. Thus, the present results suggest that the xenon cavities in proteins may act as activation centers of anharmonic dynamics at low temperatures. The individual methyl groups in a protein have different steric environments that give rise to a wide distribution of associated relaxation times. The environments of some of the methyl groups may facilitate rotational dynamics, reducing the effective barrier to rotation. This “catalytic” effect of the microenvironment around methyl groups has been seen in crystals of small peptides and other proteins74,75 and in other, nonbiological systems.76 The correlation between the location of the mobile methyl groups and the positions of xenon cavities observed in the present investigation indicates the essential role of internal cavities in protein anharmonic dynamics. This correlation is also consistent with the reduction in the flexibility of protein side chains during protein-ligand binding processes as observed in site-specific NMR experiments that use methyl groups as a probe of side chain dynamics.46 Upon binding to a protein, the ligand molecule induces conformational changes in protein that can result in a decrease of the free volumes of these cavities thereby reducing the mobility of methyl groups that were mobile in the uncomplexed protein. Further theoretical studies on the energetics of methyl groups in proteins and their temperature dependence would be very useful to provide further insight into the role of such local relaxations in understanding complex biological processes such as protein-ligand binding. Detailed examination of the contribution of methyl group dynamics to neutron scattering can also be envisaged, as has been performed for the alanine dipeptide67 and L-alanine.66 The advent of next-generation neutron sources (such as the Spallation Neutron Source at Oak Ridge National Laboratory) together with advanced deuteration facilities suggests a line of neutron scattering experimentation, involving specific hydrogenation of methyl-containing residues, such as ILE and MET.

Krishnan et al. The investigation of low-temperature quasi-elastic scattering by these residues would afford an experimental counterpart to many of the present calculations. Acknowledgment. M.K. acknowledges Dr. K. Moritsugu and Roland Schulz for useful discussions. References and Notes (1) Brooks, C. L., III; Karplus, M.; Pettitt, B. M. Proteins : A Theoretical PerspectiVe of Dynamics, Structure, and Thermodynamics; Wiley: New York, 1988. (2) Knapp, E. W.; Fischer, S. F.; Parak, F. J. Chem. Phys. 1983, 78, 4701. (3) Diehl, M.; Doster, W.; Petry, W.; Schober, H. H. Biophys. J 1997, 73, 2726. (4) Re´ at, V.; Patzel, H.; Ferrand, M.; Pfister, C.; Oesterhelt, D.; Zaccai, G. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 4970. (5) Knapp, E. W.; Fischer, S. F.; Parak, F. J. Phys. Chem. 1982, 86, 5042. (6) Igumenova, T. I.; Frederick, K. K.; Wand, A. J. Chem. ReV. 2006, 106, 1672. (7) Ishima, R.; Torchia, D. A. Nat. Struct. Biol. 2000, 7, 740. (8) Doster, W.; Cusack, S.; Petry, W. Nature 1989, 337, 754. (9) Doster, W.; Cusack, S.; Petry, W. Phys. ReV. Lett. 1990, 65, 1080. (10) Doster, W.; Settles, M. Biochim. Biophys. Acta 2005, 1749, 173. (11) Re´at, V.; Dunn, R. V.; Ferrand, M.; Finney, J. L.; Daniel, R. M.; Smith, J. C. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 9961. (12) Daniel, R. M.; Dunn, R. V.; Finney, J. L.; Smith, J. C. Annu. ReV. Biophys. Biomol. Struct. 2003, 32, 69. (13) Daniel, R. M.; Smith, J. C.; Ferrand, M.; He´ry, S.; Dunn, R.; Finney, J. L. Biophys. J. 1998, 75, 2504. (14) Daniel, R. M.; Finney, J. L.; Re´at, V.; Dunn, R.; Ferrand, M.; Smith, J. C. Biophys. J. 1999, 77, 2184. (15) Kurkal, V.; Daniel, R. M.; Finney, J. L.; Tehei, M.; Dunn, R. V.; Smith, J. C. Biophys. J. 2005, 89, 1282. (16) Brown, K. G.; Erfurth, S. C.; Small, E. W.; Peticolas, W. L. Proc. Natl. Acad. Sci. U. S. A. 1972, 69, 1467. (17) Tournier, A. L.; Smith, J. C. Phys. ReV. Lett. 2003, 91, 208106-1. (18) Tournier, A. L.; Xu, J. C.; Smith, J. C. Biophys. J. 2003, 85, 1871. (19) Hayward, J. A.; Smith, J. C. Biophys. J. 2002, 82, 1216. (20) Hayward, J. A.; Finney, J. L.; Daniel, R. M.; Smith, J. C. Biophys. J. 2003, 85, 679. (21) Austin, R. H.; Beeson, K. W.; Eisenstein, L.; Frauenfelder, H.; Gunsalus, I. C. Biochemistry 1975, 14, 5355. (22) Rasmussen, B. F.; Stock, A. M.; Ringe, D.; Petsko, G. A. Nature 1992, 357, 423. (23) Ding, X.; Rasmussen, B. F.; Petsko, G. A.; Ringe, D. Biochemistry 1994, 33, 9285. (24) Ferrand, M.; Dianoux, A. J.; Petry, W.; Zaccai, G. Proc. Natl. Acad. Sci. U. S. A. 1990, 90, 9668. (25) Tarek, M.; Tobias, D. J. Phys. ReV. Lett. 2002, 88, 138101-1. (26) Receveur, V.; Calmettes, P.; Smith, J. C.; Desmadril, M.; Coddens, G.; Durand, D. Proteins: Struct. Funct. Genet. 1997, 28, 380. (27) Roh, J. H.; Novikov, V. N.; Gregory, R. B.; Curtis, J. E.; Chowdhuri, Z.; Sokolov, A. P. Phys. ReV. Lett. 2005, 95, 038101. (28) Roh, J. H.; Curtis, J. E.; Azzam, S.; Novikov, V. N.; Peral, I.; Chowdhuri, Z.; Gregory, R. B.; Sokolov, A. P. Biophys. J. 2006, 91, 2573. (29) Hartmann, H.; Parak, F.; Steigemann, W.; Petsko, G. A.; Ringe Ponzi, D.; Frauenfelder, H. Proc. Natl. Acad. Sci. U. S. A. 1982, 79, 4967. (30) Frauenfelder, H.; Sligar, S. G.; Wolynes, P. G. Science 1991, 254, 1598. (31) Andrew, E. R.; Bryant, D. J.; Cashell, E. M. Chem. Phys. Lett. 1980, 69, 551. (32) Be´e, M.; Renault, A.; Lajze´rowicz-Bonneteau, J.; Le Bars-Combe, M. J. Chem. Phys. 1992, 97, 7730. (33) Goupil-Lamy, A. V.; Smith, J. C.; Yunoki, J.; Parker, S. F.; Kataoka, M. J. Amer. Chem. Soc. 1997, 119, 9268. (34) Hayward, R. L.; Middendorf, H. D.; Wanderlingh, U.; Smith, J. C. J. Chem. Phys. 1995, 102, 5525. (35) Morelon, N. D.; Kneller, G. R.; Ferrand, M.; Grand, A.; Smith, J. C.; Be´e, M. J. Chem. Phys. 1998, 109, 2883. (36) Hartmann, C.; Joyeux, M.; Trommsdorff, H. P.; Vial, J. C.; von Borczyskowski, C. J. Chem. Phys. 1992, 96, 6335. (37) Binbrek, O. S.; Anderson, A.; Torrie, B. H. J. Chem. Phys. 1985, 82, 1468. (38) Prager, M.; Stanislawski, J.; Ha¨usler, W. J. Chem. Phys. 1987, 86, 2563. (39) Krishnan, M.; Balasubramanian, S. J. Phys. Chem. B 2005, 109, 1936. (40) Kavitha, G.; Narayana, C. J. Phys. Chem. B 2006, 110, 8777.

Onset of Anharmonicity in Myoglobin (41) Alvarez, F.; Alegria, A.; Colmenero, J.; Nicholson, T. M.; Davies, G. R. Macromolecules 2000, 33, 8077. (42) Qi, F.; Bo¨hmer, R.; Sillescu, H. Phys. Chem. Chem. Phys. 2001, 3, 4022. (43) Best, R. B.; Clarke, J.; Karplus, M. J. Am. Chem. Soc. 2004, 126, 7734. (44) Chatfield, D. C.; Augsten, A.; D’cunha, C.; Wong, S. E. J. Comput. Chem. 2003, 24, 1052. (45) Chatfield, D. C.; Szabo, A.; Brooks, B. R. J. Am. Chem. Soc. 1998, 120, 5301. (46) Lee, A. L.; Kinnear, S. A.; Wand, A. J. Nat. Struct. Biol. 2000, 7, 72. (47) Lee, A. L.; Wand, A. J. Nature 2001, 411, 501. (48) Curtis, J. E.; Tarek, M.; Tobias, D. J. J. Am. Chem. Soc. 2004, 126, 15928. (49) Best, R. B.; Clarke, J.; Karplus, M. J. Mol. Biol. 2005, 349, 185. (50) Best, R. B.; Rutherford, T. J.; Freund, S. M. V.; Clarke, J. Biochemistry 2004, 43, 1145. (51) Tugarinov, V.; Kay, L. E. J. Am. Chem. Soc. 2006, 128, 7299. (52) Chou, J. J.; Case, D. A.; Bax, A. J. Am. Chem. Soc. 2003, 125, 8959. (53) Tugarinov, V.; Kay, L. E. Chem. Biol. Chem. 2005, 6, 1567. (54) Johnson, E.; Chazin, W. J.; Rance, M. J. Mol. Biol. 2006, 357, 1237. (55) Millet, O.; Mittermaier, A.; Baker, D.; Kay, L. E. J. Mol. Biol. 2003, 329, 551. (56) Vojtechovsky, J.; Chu, K.; Brendzen, J.; Sweet, R. M.; Schlichting, I. Biophys. J. 1999, 77, 2153. (57) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (58) Mackerell, A. D., Jr.; Bashford, D.; Bellaott, R. L.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E., III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe,

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5533 M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586. (59) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (60) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (61) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384. (62) No´se, S.; Klein, M. L. Mol. Phys. 1983, 50, 1055. (63) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (64) Prabhu, N. V.; Lee, A. L.; Wand, A. J.; Sharp, K. A. Biochemistry 2003, 42, 562. (65) Li, Z.; Raychaudhuri, S.; Wand, A. J. Protein Sci. 1996, 5, 2647. (66) Micu, A. M.; Durand, D.; Quilichini, M.; Field, M. J.; Smith, J. C. J. Phys. Chem. 1995, 99, 5645. (67) Kneller, G. R.; Doster, W.; Settles, M.; Cusack, S.; Smith, J. C. J. Chem. Phys. 1992, 97, 8864. (68) Lipari, G.; Szabo, A. J. Am. Chem. Soc. 1982, 104, 4546. (69) LeMaster, D. M. J. Am. Chem. Soc. 1999, 121, 1726. (70) Wand, A. J.; Urbauer, J. L.; McEvoy, R. P.; Bieber, R. J. Biochemistry 1996, 35, 6116. (71) Chatfield, D. C.; Augsten, A.; D′Cunha, C. J. Biomol. NMR 2004, 29, 377. (72) Kneller, G. R.; Keiner, V.; Kneller, M.; Schiller, M. Comput. Phys. Commun. 1995, 91, 191. (73) Xue, Y.; Pavlova, M. S.; Ryabov, Y. E.; Reif, B.; Skrynnikov, N. R. J. Am. Chem. Soc. 2007, 129, 6827. (74) Kitson, D. H.; Hagler, A. T. Biochemistry 1988, 27, 7176. (75) Baudry, J.; Smith, J. C. J. Phys. Chem. B 2005, 109, 20572. (76) Baudry, J. J. Am. Chem. Soc. 2006, 128, 11088. (77) Halle, B. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 1274. (78) Tilton, R. F., Jr.; Kuntz, I. D.; Petsko, G. A. Biochemistry 1984, 23, 2849. (79) Brunori, M.; Gibson, Q. H. EMBO Reports 2001, 2, 674. (80) Binkowski, T. A.; Naghibzadeg, S.; Liang, J. Nuclei. Acids Res. 2003, 31, 3352.