Glassy Dynamics of Protein Methyl Groups Revealed by Deuteron NMR
We investigated site-specific dynamics of key methyl groups in the hydrophobic core of chicken villin headpiece subdomain (HP36) over the temperature range between 298 and 140K using deuteron solid-state NMR longitudinal relaxation measurements. The relaxation of the longitudinal magnetization is weakly non-exponential (glassy) at high temperatures and exhibits a stronger degree of non-exponentiality below about 175K. In addition, the characteristic relaxation times deviate from the simple Arrhenius law. We interpret this behavior via the existence of distribution of activation energy barriers for the three-site methyl jumps, which originates from somewhat different methyl environments within the local energy landscape. The width of the distribution of the activation barriers for methyl jumps is rather significant, about 1.4 kJ/mol. Our experimental results and modeling allow for the description of the apparent change at about 175K without invoking a specific transition temperature. For most residues in the core, the relaxation behavior at high temperatures points to the existence of conformational exchange between the sub-states of the landscape, and our model takes into account the kinetics of this process. The observed dynamics are the same for dry and hydrated protein. We also looked at the effect of F58L mutation inside the hydrophobic core on the dynamics of one of the residues and observed a significant increase in its conformational exchange rate constant at high temperatures.
Introduction<!>Sample preparation<!>NMR spectroscopy<!>Main features of longitudinal relaxation<!>Theoretical background and implementation of the model<!>Distribution of methyl group activation energies and existence of conformational exchange<!>Results for F58L mutant labeled at L69 site<!>FMOC-Valine amino acid shows a much narrower distribution of activation energies<!>Relaxation data are not dependent on the hydration state<!>Conclusion<!>
<p>Hydrophobic core of proteins plays an important role in overall dynamic processes of the protein, which contribute to its biological functioning.1 A variety of approaches have been undertaken to study the dynamics of the hydrophobic cores, in which NMR spectroscopy plays a unique role due to its ability to provide site-specific resolution.2–13</p><p>Methyl groups are very sensitive probes of internal environment5,12–20 and the results on the dynamics that one obtains report not only on the rotation of the methyl groups themselves, but also on the dynamical changes experienced by the local environments of these groups. As a result, this approach is especially useful for observing kinetic and thermodynamic features of the free energy landscape evolving as a function of temperature. Various degrees of freedom inside the protein, as well as coupling with the solvent, make the landscape very complicated with many possible interactions between the energy states and a broad range of time scales for the transitions.14,20–22 Each spectroscopic technique which one employs to probe the dynamics is sensitive to a particular window of time scales, and, thus, different measurements are likely to detect flexibility of different degrees of freedom. The current challenge in the description of protein dynamics is to establish a relation between specific motional modes and protein function. For example, slow ms-μs modes have been recognized to be important in many events such as ligand binding, molecular recognition, enzyme catalysis, protein folding, dynamical transition, and protein-nucleic acid interactions.5,14,15,20,23–27</p><p>In this work, we continue to study the dynamics of the hydrophobic core using chicken villin headpiece subdomain (HP36)28 as a model system. This small globular model protein has been extensively used in various experimental and computational protein folding and dynamics studies.28–54 We combine solid-state deuteron NMR, known for its sensitivity to the details of molecular motions,3,14,16,19,20,55,56 with computational modeling techniques to look at the dynamics of methyl-bearing side chains. Due to its small size, HP36 protein can be synthesized chemically, which permits an easy introduction of selective isotope labels. Thus, the spectroscopic information that one obtains from deuteron NMR data of these samples is site-specific and has minimal effects from spin diffusion, and, as such, provides high quality data for the development of motional models. Earlier, we developed motional models for slow ms-μs time scale motions describing the dynamics of the whole methyl-bearing side chain in the 298 to 233K temperature range.53,57,58 These discrete jump models were based on both line shape and relaxation rates measurements and provided a detailed kinetic and thermodynamic description of the free energy landscape in the core More specifically, the line shapes were interpreted in terms of two modes of motions: large angle rotameric jumps between major and minor conformers of the side chains and restricted diffusion of the methyl axis on an arc.</p><p>Below 233K the line shapes are no longer sensitive to the slow dynamics of methyl side-chains and, thus, we need to interpret the relaxation data alone. The models that are self-consistent within a wide temperature range are likely to provide the best description of the dynamics in this situation. In a previous study we conducted the measurements of the longitudinal relaxation rates all the way down to 4K for a single methyl site of L69.53,54 The relaxation techniques are most sensitive to pico-second time scale dynamics.16 The motions on these time scales are governed by activation energy barriers for methyl group three-site jumps. The magnitude of the barriers reflect details of local free energy landscape involving small groups of atoms, such as individual side chains or groups of side chains within non-solvent exposed protein interior. The methyl group of L69 has demonstrated an abrupt change in the apparent activation energy value at 172K, which we previously interpreted as a change in the local environment of the methyl group.53 A second low-temperature event54 was detected at around 95K and corresponded to a drastic increase in the longitudinal relaxation times.</p><p>In this work, we investigate the main features of the activation barrier distribution for the three-site jumps of individual methyl groups in the hydrophobic core. We show that the states which give rise to this distribution are undergoing conformational (dynamical) exchange. Our key new observation is the presence of non-exponentiality in the 2H-longitudinal magnetization decay curves, the degree of which increases below ~175K. We refer to the non-exponential relaxation and the distribution of the activation energy barriers that it reflects by the term "glassy dynamics," as these are in some ways similar to the dynamical features often observed for structural glasses and polymers.59,60 The term "glassy dynamics" was also used to describe internal degrees of freedom of proteins in the context of protein folding.61,62 The dynamics of a system appear as glassy only when the conformational exchange process is on the time scale slower or on the order of the characteristic time of the experiment. Non-exponential 2H-NMR has been used to investigate glassy systems in a number of works.63–66 Roggatz et al.64 have also noted that this approach permits the probing of very slow conformational exchange processes, for which the characteristic time of the exchange is on the order of the longitudinal relaxation time. Based on our data, we developed a model which describes all of the main features of longitudinal relaxation and includes the effects of the conformational exchange. The model describes the entire temperature range of 298–140K without a need to invoke a specific transition temperature.</p><p>For this study we i) conducted extensive measurements and analysis of the magnetization decay curves for the methyl-bearing side chains of V50, L61, L69, and L75; ii) probed the dependence on the protein hydration; iii) investigated the relaxation behavior in 9-fluorenylmethyloxycarbonyl (FMOC)-Valine amino acid; and iv) looked at the effect of hydrophobic core mutation on the relaxation behavior of L69.</p><!><p>Powdered 5,5,5-d3 -fluorenylmethyloxycarbonyl (FMOC)-Leucine and d8-FMOC-Valine were purchased from Cambridge Isotopes Laboratories (Andover, MA). All leucine side chains have a single CD3 label, with a 50%/50% ratio for the two methyl groups. Commercially available Valine-d8 has all α-, β-, and γ-deuterons labeled. All protein samples were synthesized commercially by solid-state peptide synthesis with incorporation of the deuterated leucine or valine at selected sites. The samples were purified by Reverse-Phase HPLC. The identity and purity of the samples was confirmed by mass-spectroscopy and Reverse-Phase HPLC. Lyophilized powders were dissolved in water and pH was adjusted to about 6 using NaOH/HCl. All wild-type and mutant samples are expected to be folded in this range. The samples were then re-lyophilized. We have previously confirmed the refolding procedure for HP36 by measuring several chemical shift values, which agreed with the values reported for the folded protein in solution.54 The sample hydration was performed by exposing re-lyophilized powder to vapor diffusion in a dessicator until the water content reached about 35–40% by weight.</p><!><p>Data were collected on a 17.6 T spectrometer equipped with a static deuteron probe with the lowest operational temperature point of 140K. T1Z (Zeeman) measurements under static conditions were performed by the inversion recovery sequence using the multiple echo acquisition detection scheme (QCPMG).67 Briefly, QCPMG detection breaks the powder pattern spectrum into a series of spikes that follow the shape of the powder pattern. Unlike magic angle spinning, QCPMG detection does not suppress relaxation anisotropy.68 Ten to fifteen QCPMG echoes were collected with 104μs pulse spacing, corresponding to QCPMG spikelets (sidebands) spaced at 10 kHz intervals. The number of acquisitions varied between 512 and 4096 depending on the signal to noise ratio in each sample, as well as the precision of the data needed to define non-exponential decays. Seven to nine relaxation delays were collected. Relaxation data are reported for the spikelets which contain the main contribution originating from major singularities of the methyl powder patterns. This corresponds to ±10 kHz spikelets for all cases with the exception of the L75 site in the hydrated state at high temperature, for which the major singularity is at zero frequency (demonstrated in Figure S2).</p><p>As we demonstrated previously,57 the contribution of β- and α-deuterons of valine-labeled samples to the quadrupole echo line shapes can be filtered-out by employing short recycle delays, as the longitudinal relaxation times of γ-deuterons in our range of temperatures are significantly shorter than that of either β- or α-deuterons.</p><!><p>Longitudinal relaxation time measurements of methyl deuterons probe primarily pico-second time scale motions arising from three-site jumps around the C3 symmetry axis of the methyl group. For a single rate constant, there is an analytical expression which relates the value of the rate constant to the value of relaxation time T1Z.69 With complex systems such as biological molecules and polymers, one often encounters a distribution of rate constants which contribute to the apparent values of relaxation times T1Z. In this case, simulations are necessary to relate the parameters of the motional processes to observed relaxation times.</p><p>Note that the build-up of magnetization in the inversion-recovery experiment may or may not be single-exponential depending on time scales of the underlying motional processes. For example, if different sub-states within the conformational ensemble have somewhat different values of activation energies for methyl group rotations, a distribution of activation energies and rate constants is expected. Now, if the sub-states are in an exchange regime such that the rate of exchange between sub-states is faster than characteristic T1Z values corresponding to the 3-site jumps within each sub-state (referred to as the fast regime), then the apparent build-up of magnetization will be mono-exponential. The relaxation rate in this case is the population-weighted average of the relaxation rates for each sub-state. If, on the other hand, the exchange rate is slow compared to characteristic individual T1Z values, a multiple-exponential behavior is expected.</p><p>In order to investigate dynamics of main hydrophobic core methyl groups in HP36, we conducted longitudinal relaxation time measurements over a wide temperature range between 298 and 140K for the following methyl bearing side chains: V50, L61, L69, and L75. These side-chains have a similar extent of solvent-exposure47 and their location in the hydrophobic core is illustrated in Figure S1. Note that each sample had deuteron labels on only one of the residues. Side chains of leucine and valine amino acids have two methyl groups, and each of these groups contribute equally toward the NMR signal. The details of the labeling pattern are specified in the Sample preparation section.</p><p>Magnetization versus time curves for all of these methyl groups display weak non-exponential behavior for temperatures between 298K and about 175K. Below this temperature, the non-exponentiality becomes more pronounced and increases as the temperature is further lowered (Figure 1). This is the main new experimental result of this work. To our knowledge, it has not been previously characterized in detail for methyl groups in proteins and has not been explored in the earlier study on L69 residue.53</p><p>In order to characterize this behavior, the decays were fitted to a stretched-exponential function adapted for the inversion-recovery experiment: 
(1)I(t)-I(∞)=(I(0)-I(∞))e-(t/T1eff)β, in which I(t) is the signal intensity, 
T1eff is the effective relaxation time, and β is the parameter which reflects the degree of non-exponentiality, 0 < β ≤ 1. β less than 1 corresponds to a non-exponential behavior. Eq. (1) is one of the typical empirical functions that are used in descriptions of various dynamic effects in glasses.59,63 In the absence of conformational exchange, it can be viewed as an integral of ordinary mono-exponential relaxation functions over a distribution of the relaxation rates.70 The improvement in the quality of the fits with this function, as compared to the mono-exponential function, is illustrated by the plots of the residuals (bottom panels, Figure 1).</p><p>The values of 
T1eff and β are shown in Figure 2 for samples in the dry state. At high temperatures, the values of β range between 0.9 and 1, while below about 175K there is a gradual decrease in the values of β. The temperature dependence of 
T1eff shows an asymmetrical behavior and displays a higher value of the minimum compared to a single motional process (see also Figure S1 of Vugmeyster et al.53 for details). The non-Arrhenius temperature dependence coupled with non-exponential relaxation appears to be typical for systems undergoing glassy dynamics.59,63,71</p><p>The measurements shown in Figures 1 and 2 were done with the multiple-echo acquisition scheme,67 referred to as QCPMG (quadrupole Carr–Purcell–Meiboom–Gill). This scheme, described in more detail in the Materials and Methods section, was used for the enhancement of signal-to-noise ratio. The values shown correspond to the QCPMG spikelets around the position of the major singularities of the powder patterns. Typical QCPMG spectra are shown in Figure S2 of the Supporting Material and quadrupole echo powder patterns can be found in Vugmeyster et al.57 Each QCPMG spikelet echo includes signals from a range of crystallite orientations,67 which should lead in itself to non-exponentiality. This effect, however, is negligible, with β in the range of 0.98 <β <1 for correlation times in the pico-second to micro-second range. QCPMG detection, coupled with long relaxation delays, has allowed for detection and detailed characterization of the effects of the non-exponentiality in the magnetization decay.</p><!><p>The results of the relaxation behavior depicted in Figure 1 and 2 can be most naturally explained in terms of a model that assumes a distribution of conformational states. Each of the states has a different activation energy barrier for the three-site jumps, stemming from slight differences in the methyl environment. A methyl group samples the distribution of activation energies in the process of conformational exchange.</p><p>Note that in our model the states are distinguished by their values of the activation energies Ea, but not the potential energies. Assuming also a continuous distribution of conformations, we introduce time-dependent magnetization density m(Ea,t) such that ∫m(Ea,t)dEa =M (t) is the overall longitudinal magnetization. At equilibrium meq(Ea)= Meq f(Ea), where f(Ea) is the probability density for the system to be in the state with the activation energy Ea.</p><p>The decay of the longitudinal magnetization in the presence of conformational exchange can be described by Bloch-McConnell72 equations adapted for a continuous distribution of conformational states</p><p>
(2)dm(Ea,t)dt=-(m(Ea,t)-meq(Ea,t))T1Z(Ea)-∫dEa′(W(Ea,Ea′)m(Ea,t)-W(Ea′,Ea)m(Ea′,t)), in which T1Z (Ea) is the relaxation time determined by the activation energy Ea, meq (Ea) is the equilibrium magnetization density and 
W(Ea,Ea′) is the probability density of the transition from the state with activation energy Ea to that of 
Ea′. The initial conditions for m(Ea,0), depend on a particular experimental scheme and are described below.</p><p>The transition probability density 
W(Ea,Ea′) obeys the detailed balance condition given by 
W(Ea,Ea′)f(Ea)=W(Ea′,Ea)f(Ea′). As usual, the detailed balance condition can be viewed as following from the microscopic reversibility of the exchange process, which is described by the symmetric exchange rate constant 
K(Ea,Ea′)=K(Ea′,Ea). Then 
W(Ea,Ea′) can be expressed as 
W=K(Ea,Ea′)f(Ea′), in which K (E, E′) is independent of the density of states. Substituting into Eq. (2) gives</p><p>A simple choice for the probability density is a Gaussian form 
f(Ea)=12πσ2exp(-(Ea-〈Ea〉)22σ2). It is important to note that we assume f (Ea) to be temperature independent, and thus our model does not take into account the temperature dependence of the probability density originating from different potential energies of the states. We distinguish the states only by their values of the activation energies for the three-site jump process.</p><p>To perform numerical calculations and also to demonstrate a more direct connection with the usual form of the McConnell equations, we want to rewrite Eq.(3) in a discrete form. We introduce an equally spaced grid of the values of activation energy Ea,i (i = 1..N) and define mi =m(Ea,i,t), T1Z,i = T1Z (Ea,i), and Kij =K (Ea,j,Ea,i). Probability density is changed to probability distribution according to 
fi=f(Ea,i)/∑j=1Nf(Ea,j)≈f(Ea,i)ΔEa, where ΔEa is the grid step. Integration can be exchanged for summation. This requires that the variation in the integrand in Eq.(3) is small, or at least not large, on the ΔEa scale. Thus we have</p><p>Two limiting cases of Eq.(4) can be considered. In the so-called "strong collision limit," the rate of transition between the states does not depend on either initial or final values of the activation energies and, thus, we have Kij =kex. In this case Eq.(4) becomes</p><p>This limit has been considered in the description of relaxation of guest molecules in glass formers by Roggatz et al.64</p><p>It seems more likely, however, that short-range stochastic jumps in the activation energy space (diffusion limit) is a more appropriate description leading to the nearest neighbor interaction model. This model assumes that Kij =kD for |i−j| = 1 and Kij =0 otherwise. We further assumed that Kij is the same for all of the nearest neighbor pairs. Eq. (4) then can be rewritten as</p><p>
(6)dmidt=-mi-mieqT1Z,i-kDmi(fi+1+fi-1)+kDfi(mi+1+mi-1) where we set f0 = fN+1 =0 and m0 =mN+1 =0.</p><p>In terms of the continuous distribution, nearest neighbor jumps imply that 
K(Ea,Ea′) is different from zero only for 
∣Ea-Ea′∣~ΔEa, which ensures that the discretization procedure is valid in this limit. We can treat this case as diffusion in the activation energy space. The diffusion coefficient D(Ea), which depends on the activation energy, can be obtained by the standard procedure reducing Eq.(6) to a diffusion type equation. This will lead to 
D(Ea)=kDf(Ea)ΔEa3, in which we used fi ≈f (Ea,i)ΔEa. We introduce, then, a parameter characterizing the magnitude of the exchange rate, kex, as the inverse of the time for the diffusion process with the average diffusion coefficient 〈D(Ea)〉 to reach the characteristic width of the distribution σ</p><p>Defined in this way, kex can be used as a model parameter which does not depend on the number of grid points in the simulations.</p><p>The initial values of mi (t) are set by assuming an ideal inversion in the inversion recovery experiment: 
mi(0)=-mieq∝fi. We have selected a grid of activation energies which covers −3σ to +3 σ range with equal spacing, 
Ea,i=〈Ea〉+3σ(2i-1N-1-1). We have selected N = 11, for which case kex =0.06kD, as follows from Eq.(7).</p><p>The relaxation time in each conformer T1Z,i is defined by fast three-site jumps of the methyl group around the methyl axis, T1Z,i =T1Z (ki) .69 Since we have used the QCPMG detection scheme, strictly speaking we cannot directly apply the analytical formula for the 3-site hops because each point in the QCPMG spectrum has a different mix of frequencies compared to the regular quadrupolar echo detected powder pattern. Relaxation time anisotropy will then lead to somewhat different effective relaxation times compared to the quadrupolar echo detection scheme.68 We, thus, simulated the inversion recovery relaxation experiment with QCPMG detection using EXPRESS simulation program73 in order to establish a correspondence between the values of ki and T1Z,i. The value of quadrupolar coupling constant Cq was taken as 160 kHz53,54 We also note that the effect of the motions on the order of Cq determined by the line shape measurements is relatively minor, under about 7%,53 and was not included in the simulations. The individual three-site jump rate constants ki are assumed to follow the Arrhenius law k(T) =k0e−EaT, in which we take k0 to be the same for all states.</p><p>The simulated functions mi(t) are summed to obtain total simulated magnetization. 
M(t)=∑i=1Nmi(t)ΔEa. M(t) was fitted by the stretch-exponential function defined in Eq. (1) to yield simulated values of 
T1eff and β. We have also confirmed that the modeled values of 
T1eff and β do not change when a different value of N (N = 21) is taken.</p><p>We distinguish three motional regimes defined on the basis of the comparison between the values of kex and 
1/T1eff. In the fast regime 
kex≫1/T1eff, in the slow regime 
kex≫1/T1eff such that the distribution becomes static, and in the intermediate regime 
kex~1/T1eff. The effect of the dynamical averaging is to lower the value of 
T1eff and raise the value of β, as demonstrated in Figure 3. For the fast regime β ≈ 1 and 
1/T1eff=〈1/T1Z(Ea)〉, where 〈…〉 denote ensemble average. If the system remains in the slow regime for all temperatures, the temperature dependence of 
T1eff and β is driven by the temperature dependence of the width of the distribution of the three-site jump rate constants, given by σ/T according to the Arrhenius law. We note that around the minimum of 
T1eff (at 1000/T ≈ 5 for the parameters shown in Figure 3) there is an elevation in the values of β, which reflects the fact that around this point the relaxation times are weakly dependent on the three-site jump rate constants. In general, the change between the three regimes as a function of temperature can be somewhat complicated, because, while kex has monotonous dependence on temperature, 
T1eff does not. For the purpose of demonstrating the relaxation behavior in Figure 3, we have chosen the Arrhenius temperature dependence for kex, kex (T) kex,0e−Eex/T. If we limit our consideration to the case for which Eex is larger than 〈Ea〉 and also σ ≫ 〈Ea〉, then there is only one temperature at which 
kex=1/T1eff. In this case, for significantly higher temperatures, the system is in the fast regime while for lower temperatures it is in the slow regime. There is also some interval of temperatures for which the system is in the intermediate regime, as reflected in the cross-over between the two limiting behaviors of 
T1eff and β. In Figure 3 we give two examples of such a cross-over. The dependence of 
T1eff and β on the exact value of kex is most pronounced in the intermediate regime; for the other two regimes it is only important that it is either much larger (fast regime) or much smaller (slow regime) than 
T1eff.</p><!><p>In principle, it is possible to perform a direct fit of the simulated magnetization decay curves to the experimental data. However, we have chosen to fit the simulated magnetization decay curves to the stretched-exponential function, and then compare the obtained simulated values of 
T1eff and β to the experimental ones. This approach permits the capturing of the most essential features of the relaxation behavior.</p><p>The static case (kex = 0) of our model (depicted in Figure 3) appears to capture many of the features of the experimental relaxation behavior, in particular the increase of the non-exponentiality at low temperatures. Therefore, we first consider the fits to the model with kex set to zero for all temperatures. Our fitting parameters are the mean value of the activation energy <Ea>, the width of the distribution of activation energies σ, and the asymptotic value of the three-site hop rate constants at high temperatures k0.</p><p>The fits to this case, obtained via χ2 minimization, are shown as a solid line in Figure 2. This limit is sufficient to describe the relaxation behavior below the temperature of the minimum of 
T1eff and the apparent increase in the non-exponentiality below 175K for all residues. We are thus capable to describe this feature without the need to postulate a specific transition temperature, as was suggested in the earlier work.53</p><p>The central values of the activation energies and the widths of the Gaussian distributions are similar for all residues (Table 1). The values of the activation energies <Ea> in the range of 11–12 kJ/mol are typical for methyl groups in proteins.17 The values of the width of the distribution σ are ~1.4 kJ/mol. They can be compared to a typical distribution across multiple residues in proteins, which was determined to be around 3.8 kJ/mol based on 33 residues in α-spectrin SH3 domain.17 The obtained values of σ are also comparable to the difference between <Ea> for the four hydrophobic core residues in HP36.</p><p>For all residues except V50, the presence of dynamical averaging due to conformational exchange is apparent for higher temperatures, for which the data deviates from the static limit of the model. As elaborated in the previous section, the main fine effect of the dynamical averaging is to induce a steeper rise in the values of β when the exchange regime shifts from the static limit to the intermediate or fast. Thus, the data suggest the full model with the inclusion of the nonzero conformational exchange rate is necessary to describe the relaxation at high temperatures.</p><p>The relaxation behavior is most sensitive to exact values of kex only in the intermediate regime, at which is comparable to 
T1eff. For other regimes one can establish either the lower limit of kex (for the fast regime) or the upper limit (for the slow regime). For high temperatures at which the relaxation data deviates significantly from the static limit, we have employed the full model with the inclusion of the exchange effect. The intermediate regime is already apparent at the highest temperature of 298K for L61 and L69, while for L75 it occurs at a lower temperature of about 283K. At 298K L75 is already in the fast regime. The fits are demonstrated by the squares in Figure 2 and the values of kex for the intermediate regime temperatures are given in Table 1. For L61 and L69, the values of kex were determined at three temperatures. It is not possible to infer how kex varies with temperature within our precision.</p><p>L75 residue is the only methyl group in HP36 that is explicitly associated with its biological function of F-actin binding. It is interesting to note that this residue displays the highest rate constant at room temperature for both the rotameric motional mode and the mode detected in this work; thus, it is possible that the role of this residue in F-actin binding process is achieved by kinetic control.</p><p>The question remains regarding the physical basis of the observed conformational states, which we cannot answer based on the present work alone. Our model describes the main features of the relaxation behavior by introduction of multiple states with different activation energy values. However, this model does not connect the values of the activation energies to a structural description of these states. In principle, the states can arise either due to somewhat different orientations of the methyl axis or could include a scenario in which the activation energy is modified not by the changes in the orientation of the methyl axis itself, but rather by subtle reorientations of the methyl group surrounding. Glassy dynamics in polymers are often caused by large amplitude intramolecular motions of different segments within the same molecule that "freeze" with the change in temperature or other conditions.74 As a result, one obtains multiple states which differ in relative orientations of large molecular segments. Unlike polymers, globular folded proteins have a well-defined structure, especially in the core region. Therefore, the origin of multiple states in a well-structured hydrophobic core is expected to be rather different from the scenario in polymers.</p><p>As mentioned in the introduction, the same sites that are probed in this study were investigated by deuteron line shape analysis.57 However, we do not think that the distribution of states observed in this work is due to the motional modes detected from the line shape data. The primary reason is the difference in the time scales. The line shape technique is sensitive to time scales on the order of the quadrupolar coupling constant. Its effective value averaged over the three-site jumps is 53 kHz, which is much larger than 
1/T1eff and, therefore, any processes that contribute to the change in the line shape would effectively average the relaxation rates without giving rise to non-exponentiality. In Vugmeyster et al. 57 the line shapes were modeled by two modes, representing rotameric jumps and restricted diffusion of the methyl axis on an arc. An additional reason against interpreting the rotameric motion as the basis for the distribution of the three-site jump rates is that the populations of the alternative rotameric conformers were found to be small at low temperatures, such that they can not give rise to the observed non-exponentiality. The conclusion that the populations of alternative rotamers are minor at low temperatures is also consistent with the 13C magic angle spinning spectra of HP35 in frozen solutions at 25K by Hu et al.29 The fact that the conformational exchange process studied in this work occurs on a relatively slow time scale also means that this process has a negligible effect on the T1Z anisotropy, unlike the motional modes on the order of the quadrupolar coupling constant, which make a slight but measurable reduction in the T1Z anisotropy.53</p><p>As the width of the distributions σ are similar for all of the four residues located in different parts of the hydrophobic core, it is possible that these widths, which reflect different local conformational sub-states, are influenced by a common factor, such as the packing of the core. We also observe that the exchange rates constants kex at high temperatures are different for the four residues, and thus the kinetics is less influenced by the common environment of the core and is more sensitive to local conditions. The shorter side-chain of V50 compared to the longer leucine side-chains may be a factor why this residue has the lowest exchange rate constant, assuming that the mobility of the side-chain as a whole facilitates the exchange. While we cannot quantify the changes in the torsional angles and distances corresponding to the various substates, it is likely that these are of moderate amplitudes due to the possibility of the common packing factor as well as the time scales on the order of tens of milliseconds. Future advances in accelerated molecular dynamics techniques may pinpoint the structural basis for this phenomenon.</p><!><p>Three phenylalanine residues form the basis of the hydrophobic core in HP36 and are flanked by several methyl groups.30,47,48 In order to probe the effect of the removal of a neighboring phenylalanine ring on the observed dynamics of HP36, we compared the dynamics of L69 in the wild-type protein to the F58L mutant, in which the ring closest to L69 has been replaced by another hydrophobic side chain (Figure 4A). Note that previous work by Raleigh and McKnight groups41,42 identified that the F58L mutation preserves the main fold of HP36, but lowers the thermostability of the protein. Our earlier deuteron line shape measurements in the 298–233K temperature range57 indicated that this mutation increases the mobility of L69 side chain by increasing the rate constant and the activation energy for the rotameric jumps, as well as increasing the length of the arc for the restricted diffusion mode.</p><p>The results for the relaxation measurements (Figure 4B) indicate that the main difference between the wild-type and the mutant sample lies in the overall increase in the values of β at high temperatures, which immediately suggests faster conformational exchange. Indeed, if we fit the high temperature data using the full model presented above, we find very high conformational exchange rates at high temperatures. The exchange is in the fast regime (with respect to the relaxation times) until about 273K, and then slows down to the intermediate regime for 254 and 243K. The lower range value of kex at room temperature is 100 s−1 for F58L, as compared to the best-fit value of 21 s−1 for the wild-type protein.</p><p>These results, combined with the previous line shape measurements, indicate that the removal of the phenylalanine ring in the proximity of L69 significantly alters the kinetics and thermodynamics of the local landscape toward increased flexibility of the side chain at high temperatures. At lower temperature this flexibility is quenched due to increased values of activation energies for both the conformational exchange processes detected here and the rotameric jumps as assessed from the line shapes.</p><!><p>To investigate whether the observed glassy behavior in the protein requires the complexity of the protein environment, we also look at the NMR relaxation in amino acids that are modified with the FMOC group at the amino side. This bulky group creates a hydrophobic environment in the proximity of methyl groups, yet it does not contain all the complex elements of the protein cores.</p><p>We previously have found that FMOC-Leucine lacks the apparent change in the slope in the 
T1eff curve at ~175K.58 As the analysis of the data for FMOC-Leucine is somewhat complicated by the presence of two equally populated rotamers, we decided to investigate FMOC-Valine amino acid.</p><p>As shown in Figure 5, neither the change in the slope in the 
T1eff curve nor the onset of strong non-exponentiality is present in this sample. In addition, there is no elevation of the minimum compared to the one calculated for a single three-site jump rate. This immediately suggests a much narrower distribution of the activation energy barriers. The static limit is sufficient to fit the data with σ = 0.6 kJ/mol, which is more than two times smaller than the values found for the protein sites. This suggests that the complexity of the protein environment allows for a wide distribution of the activation energies, and, as a result, a more complex relaxation behavior.</p><!><p>The behavior depicted in Figures 1 and 2 is independent of whether the protein is in the hydrated or dry state for all sites probed, as illustrated in Figure 6 for L69. The hydrated state is 35–40% water content by weight, which is a typical amount corresponding to one layer of water molecules and thus defines the threshold at which most globular proteins are fully hydrated.75,76</p><p>Hydration has been found to have a strong effect on deuteron line shapes of these residues and the removal of water strongly reduced the dynamics of the rotameric jumps and restricted diffusion, probed in the 298–253K temperature range. We note that there is no contradiction between the line shape and the relaxation data. The direct contribution to the relaxation is under ~7 % for motions on the time scale on the order of Cq. Any distribution of the activation energy barriers that can be potentially created by these modes will not affect the apparent relaxation time because of the fast motional regime with respect to the relaxation time.</p><!><p>Deuteron Zeeman order relaxation behavior of several hydrophobic core methyl-bearing side chains in HP36 protein exhibited an intriguing temperature dependence. At high temperatures there is a tendency toward weak non-exponentiality for most sites, while below ~175K the degree of the non-exponentiality starts to increase when the temperature is lowered. We have modeled this behavior assuming the existence of a distribution of the activation energies for the three-site jumps of the methyl groups which corresponds to various sub-states with slightly different methyl environments. This model naturally predicts the apparent change of slope in the relaxation times curves below 175K, noted in the previous work,53 without postulating a specific transition temperature. The parameters of the distribution are similar for all side-chains. We note that the width of the distribution of the three-site jumps activation energy is relatively large, about 1.4kJ/mol, and is comparable to the difference in the central values of the activation energies among the residues. The width of the distribution is much narrower in the FMOC-Valine amino acid, which leads to very different relaxation features compared to the protein.</p><p>For most residues, with the exception of V50, the relaxation behavior at high temperatures suggests the existence of conformational exchange between the sub-states. The conformational exchange kinetics between the sub-states is taken into account in a quantitative manner with the use of McConnell equations, which are tailored for the description of conformers distinguished only by their values of the activation energies. At room temperature, the dynamics of L69 and L61 are in the intermediate regime, in which the exchange rate constants are comparable with the inverse of the longitudinal relaxation times (approximately 50ms at 298K). L75 residue, which is involved in F-actin binding, displays the highest value of the conformational exchange rate constant at room temperature and is in the fast regime with respect to the relaxation time.</p><p>Interestingly, the observed relaxation features are not dependent on the hydration state of the protein, indicating that interactions with the hydration layer are not necessary for the presence of the glassy behavior of methyl groups.</p><p>We also looked at the effects of F58L mutation inside the hydrophobic core on the behavior of leucine-69 side chain and found that the mutation has a profound effect on the kinetics at this site by significantly increasing the exchange rate constant at higher temperatures. Taken together with the previous line shape data, this result leads to the conclusion that the mutation significantly affects both the kinetics and thermodynamics of slow motions inside the hydrophobic core. Site-specificity and dependence on mutations of slow motions inside the core of HP36, as well as clear differences in comparison to FMOC-Valine and FMOC-Leucine amino acids, may indicate that the mobility of the core is optimized for the biological function of the protein.</p><!><p>Figure S1. Ribbon diagram of HP36 with key hydrophobic core residues labeled. Figure S2. QCPMG spectra for the four core methyl-bearing residues. This information is available free of charge via the Internet at http://pubs.acs.org</p>