INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS REFERENCES

The temperature dependence of the internal motion of proteins and its relation to function have been much studied using a variety of experimental and theoretical techniques (e.g., Knapp et al., 1982; Parak and Knapp, 1984; Doster et al., 1989; Smith et al., 1990; Smith, 1991, 2000; Rasmussen et al., 1992; Demmel et al., 1997; Fitter et al., 1997; Cordone et al., 1998; Ostermann et al., 2000). These studies have revealed a smooth change in the slope of the temperature dependence of the average mean square displacement of atoms in proteins, u², in the temperature range ~170-240 K. At temperatures below this transition, the dynamics appears harmonic and vibrational. Above the transition, anharmonic motions are also present. This change in dynamics resembles that seen in the liquid-glass transition (Angell, 1995). The anharmonic dynamics involve diffusive motion and allow sampling of different conformational substates (Elber and Karplus, 1987; Karplus and Petsko, 1990; Frauenfelder et al., 1991; Kneller and Smith, 1994). Some studies have found a correlation between the onset of the anharmonic motion and protein activity (Parak et al., 1980; Rasmussen et al., 1992; Ferrand et al., 1993; Ding et al., 1994; Ostermann et al., 2000).

In the present paper, we use computer simulation to examine two aspects of this temperature dependence. The first concerns the role of solvent. Experimental and computational work has indicated that the activation of anharmonic dynamics is strongly influenced by solvent (Ferrand et al., 1993; Demmel et al., 1997; Fitter et al., 1998a,b; Fitter, 1999; Vitkup et al., 2000; Réat et al., 2000). However, the question as to whether solvent is required for the activation of anharmonic dynamics is as yet unanswered, i.e., it is not known whether an isolated protein in the absence of solvent would have any transition from harmonic to anharmonic dynamics, or whether it remains harmonic and vibrational up to room temperature. Neutron-scattering experiments on dry protein samples show no clear evidence of dynamical transition behavior (Ferrand et al., 1993; Fitter et al., 1997, 1998b; Lehnert et al., 1998; Fitter, 1999; Réat et al., 2000). However, a definitive answer to this question requires the obtention of a harmonic reference state. Here, molecular dynamics (MD) simulations are performed on a model system, an isolated bovine pancreatic trypsin inhibitor (BPTI) molecule, as a function of temperature. To set a harmonic reference, normal mode analysis (NMA) is performed. In NMA, the harmonic approximation to the potential energy function is made, and the resulting atomic trajectories are a superposition of vibrational normal modes. In contrast, in MD simulations, the full anharmonic potential energy function is used. The temperature dependence of u² is calculated here from both NMA and MD. Deviation of the MD displacements from the NMA results indicates the presence of anharmonic dynamics.

The second question addressed here concerns the relationship between temperature-dependent dynamical quantities in proteins and neutron-scattering measurements. Neutron scattering is a technique of major importance for examining dynamical transition properties and is complementary to MD/NMA, probing similar length and timescales (Lovesey, 1984; Smith, 1991). Much of the dynamical transition work has involved examination of elastic neutron scattering (e.g., Ferrand et al., 1993; Fitter et al., 1997; Daniel et al., 1998, 1999; Lehnert et al., 1998; Cordone et al., 1999; Dunn et al., 2000; Réat et al., 2000; Zaccai, 2000). Atomic motions lead to a reduction in elastic scattering intensity in a manner closely analogous to that seen in crystallographic Bragg scattering. However, the interpretation of the elastic neutron-scattering data in terms of u² involves approximations in the description of the dynamics of the atoms in a protein and its relation to scattering amplitudes. In view of the volume of elastic neutron scattering work that has recently appeared, it is now timely to examine the approximations made. Therefore, here we use the MD simulation as a test "experimental data set." The measurable quantity, the dynamic structure factor, is calculated from the simulation trajectories and then subjected to data analysis similar to that performed experimentally. The effect of the approximations on the derived u² can then be determined by comparison with u² calculated exactly from the MD simulations.

Two approximations in particular are examined. The first of these is the so-called "Gaussian approximation" for the q-dependence of the dynamic structure factor (where q is the scattering vector). Second, the effect of the existence of a distribution of mean square displacements in the protein, rather than one uniform value, is explored. This "motional heterogeneity" leads to an additional non-Gaussian contribution to the dynamic structure factor. Finally, the effect of the instrumental energy resolution on the u² obtained is also probed. It is found that each of these effects reduces the u² obtainedthe present calculations provide estimates of how large these effects are and under which experimental conditions they are likely to exert an influence.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS REFERENCES

Molecular dynamics simulations

The model system of BPTI consists of 892 atoms, 58 residues, and 4 internal water molecules. This model system is an approximation of a dry powder sample, as has been used in several neutron experiments (Ferrand et al., 1993; Fitter et al., 1997, 1998b; Lehnert et al., 1998; Fitter, 1999; Réat et al., 2000). To mimic typical experimental conditions, the exchangeable hydrogen atoms were replaced by deuterium, leaving 315 hydrogen atoms.

The system was simulated using CHARMM (Brooks et al., 1983) version 27 with all-atom parameter set 22 (Mackerell et al., 1998). The four internal, structural water molecules were represented by the TIP3P potential (Jorgensen et al., 1983). The simulations were performed in the NVE ensemble. A timestep of 0.002 ps was used with SHAKE (Ryckaert et al., 1977) applied to constrain bonds containing hydrogen or deuterium atoms. Nonbonded and electrostatic interactions were truncated using the shifting function (Steinbach and Brooks, 1994) at 13.0 Å and the dielectric constant was set to ₀.

Simulations were performed at 19 different temperatures: 80, 100, 120, 140 K, then in steps of 10 K to 280 K, and, finally 300 K. The starting structure for the first temperature simulated, 80 K, was the energy-minimized BPTI crystal structure (Parkin et al., 1996) Protein Data Bank reference (1BPI) (Berman et al., 2000). The starting structures for the rest of the simulations at increasing temperatures were the final structures from the preceding temperature. The systems at each temperature were equilibrated for 50 ps and then data collected every 0.1 ps for 2000 ps, i.e., a total simulation time of 2050 ps per temperature. The total simulation time was 38.95 ns, requiring 209 CPU hours running in parallel on 64 processors on an IBM SP2 computer, i.e., 13,376 CPU hours. After the final simulation at 300 K the root mean square deviation of the backbone heavy atoms from the crystal structure was 2 Å, and that calculated using the secondary structural elements only was within 1 Å of the crystal structure.

Normal modes

Normal modes were computed using the VIBRAN module of CHARMM using the same model system and potential function as for the MD simulations. At low temperatures, proteins are trapped in multiple minima leading to structural inhomogeneity (Elber and Karplus, 1987; Frauenfelder et al., 1991; Ostermann et al., 2000). This leads to dynamical inhomogeneity reflected in broadening of theoretical neutron and infrared spectra (Lamy et al., 1996). To take this inhomogeneity into account, 20 structures were taken from the 300-K simulation trajectories at regular intervals and were energy minimized to a root mean square energy gradient of <10⁷ kcal mol¹ Ź. Normal modes were computed in each minimum.

The eigenvalues and eigenvectors from the 20 normal-mode analyses were used to calculate the mean square displacements. Both classical and quantum mechanical mean square displacements were calculated and compared. The classical mean square displacements were calculated and compared. The classical mean square displacement averaged over the atoms in the protein, u²_class, is given by

(1)

where N_atom is the number of atoms, is the eigenvector for atom  in mode k, k_B is the Boltzmann constant, T is the temperature, m is the atomic mass, and _k is the mode eigenvalue, or frequency. The equivalent quantum mechanical expression is

(2)

where  = h/2 and h is the Planck constant.

Neutron-scattering properties

The nMOLDYN package (Kneller et al., 1995) was used to calculate neutron-scattering properties from the MD trajectories.

Dynamic structure factor

(3)

(4)

(5)

(6)

(7)

(8)

,coh

,inc

24

(9)

(10)

Gaussian approximation

(11)

(12)

(13)

Motional heterogeneity

(14)

Timescale dependence

(15)

(16)

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS REFERENCES

Dynamical transition

In Fig. 1, u²(t) is shown at times, t, ranging from 10 to 1000 ps for each of the 19 temperatures sampled. The quantity is shown from the MD simulations, together with the quantum and classical u² from the harmonic analyses. The quantum-mechanical u² is ~30% larger than u²_class at 80 Kthis difference drops to below 10% at 300 K. u²_qm is formally nonlinear in T (see Eq. 2). However, the deviation of u²_qm from linearity in Fig. 1 is clearly small, and therefore u² T is a good approximation in the harmonic regime.

View larger version (21K): [in this window] [in a new window]   FIGURE 1   Temperature dependence of average mean square displacements, u2(t), from MD simulations for times, t, of () 10 ps, () 100 ps, () 500 ps, and (*) 1000 ps. Values are also shown from the normal mode analysis, in the classical (), and quantum mechanical (- - - -) cases. The error bars on the normal mode results are two times the standard deviation in the harmonic mean square displacements over the 20 normal mode analyses. The u2 values are weighted by the atomic incoherent scattering lengths.

At 80 K, the classical MD u²(t) at 10, 100, 500, and 1000 ps and the harmonic u²_class all coincide. The lowest frequency vibration in the harmonic analysis is 6.7 cm¹, corresponding to a period of 5 ps. Therefore, the modes are adequately sampled in timescales as low as 10 ps. Consequently, at 80 K, the dynamics is entirely described by fast vibrational motions in the protein.

The MD u²(t) at temperatures above 100 K is timescale dependent. Already at 150 K, one finds u²(10 ps) < u²(100 ps) < u²(500 ps) < u²(1000 ps), indicating the presence of motions over a range of timescales between 10 and 1000 ps. The timescale dependence is considerably larger at higher temperatures, and at 270-300 K u²(500 ps) and u²(1000 ps) are about twice as large as u²(10 ps). In Fig. 2 is shown the time evolution of u²(t) at selected temperatures. This confirms the presence of significant long-timescale motion at T  150 K. 

View larger version (23K): [in this window] [in a new window]   FIGURE 2   Time evolution of the mean square displacement, u2(t) from MD simulation at temperatures of 80, 150, 240, and 300 K.

The MD profiles in Fig. 1 for 100 ps and longer times deviate significantly from the classical harmonic results for T > 100 K, indicating the presence of anharmonic dynamics above 100 K. For these timescales, anharmonic motions contribute the major part of u²(t) at 300 K. The faster, 10-ps timescale dynamics starts to deviate from harmonic behavior at a higher temperature, T > ~180 K. That the onset of fast anharmonic dynamics is at a higher temperature than for slow dynamics is consistent with experimental results on enzyme solutions that also demonstrate a timescale dependence of u² (Daniel et al., 1999; Dunn et al., 2000). At timescales greater than 100 ps, the u²(t) temperature profiles are considerably rougher at high temperatures than at shorter timescales. This reflects the presence of long-timescale (10¹⁰ s) structural transitions in the protein that are incompletely sampled in the simulations. The presence of 10¹⁰-s timescale motions at high temperatures is confirmed in Fig. 2.

To our knowledge, in powdered dry protein systems, deviation from harmonic behavior has not been reported experimentally. The present results suggest that this may be due to the absence, in the experimental analyses, of a harmonic reference profile, such as is provided here by the NMA. Indeed, visual inspection suggests that the MD results in Fig. 1 should be able to be reasonably well fitted by straight lines up to 240 K, but with steeper gradients than the NMA profile. In the absence of the harmonic NMA baseline, then, the approximate linearity of the MD profiles would lead to the erroneous conclusion that little anharmonic behavior is present. However, as the present MD/NMA comparison shows, significant deviation from harmonic dynamics is present at temperatures as low as ~100 K.

At higher temperatures (>240 K) a deviation from linearity is seen, especially on the longer timescales, but which is at most 10% (at 300 K). This deviation is due to increased flexibility in the loops at the C-terminal and N-terminal ends of the protein.

Gaussian approximation

The validity of the Gaussian approximation, i.e., Eqs. 12-13, is now explored by examination of the scattering calculated from individual atoms. We examine the scattering from the 300-K trajectory. At lower temperatures, the scattering may be somewhat closer to Gaussian due to reduced anharmonicity of the dynamics. S(, ) was calculated for each individual atom from the 300-K MD trajectory using Eqs. 7 and 8 without invoking the Gaussian approximation. The S(, ) thus obtained was convoluted with an instrumental energy resolution function (see Sect. Dynamic Structure Factor above). Two different typical energy resolutions were used, corresponding to 100 and 20 µeV full width half maximum. Subsequently, the elastic scattering intensity, S(q,  = 0), was estimated by integrating S(, ) over ± where is the half width half maximum of the energy resolution function. Deviation from linearity in ln S(q,  = 0) vs. q² then indicates the presence of non-Gaussian scattering.

Plots of ln S(q,  = 0) vs. q² calculated with 20-µeV energy resolution are shown for two extreme cases in Fig. 3. Atom H from the backbone of ARG 39 has the lowest H-atom u, of 0.6 Å², and a linear plot, demonstrating Gaussian behavior. In contrast, an H atom from the side chain -CH₃ in MET 52 has the largest u, of 19 Å² and clearly exhibits non-Gaussian behavior over the q-range presented.

View larger version (17K): [in this window] [in a new window]   FIGURE 3   Normalized ln S(q,  = 0) vs. q2 profiles for (), the hydrogen atom with the lowest u = 0.6 Å2; and (- - - -), the hydrogen atom with the highest u = 19 Å2. The elastic scattering intensity was calculated using a 20-µeV energy resolution function at 300 K.

The Gaussian approximation is strictly valid for q  0. Therefore, the accuracy of the approximation will depend on the q-range examined. We now obtain an estimate of the proportion of atoms exhibiting non-Gaussian dynamics at 300 K in two different q ranges and at two different instrumental energy resolutions. The q²-ranges chosen are 0.01 < q² < 1.44 Å² (hereafter called "low-q") and 0.16 < q² < 23.04 Å² (hereafter called "high-q"). These are similar to the ranges commonly used in obtaining u² from neutron-scattering data. For example, at the Institut Laue-Langevin, Grenoble, the low-q range is available on the instruments IN6 and IN16 and the high-q range is available on IN13.

In each of the four cases, u² was obtained by fitting a straight line to ln S(q,  = 0) vs. q². In Fig. 4, the goodness-of-fit (R²) for each of the 315 hydrogen atoms are shown. R² = 1 indicates a perfect fit. In the low-q range (Fig. 4, a and b), R² is above 0.95 for 97% of the atoms, indicating that the Gaussian approximation is valid for these atoms over this q range. For low-q, there is no clear correlation between R² and u. In contrast, in the high-q range (Fig. 4, c and d), ~20% of the atoms have R² worse (smaller) than 0.95. There is also a correlation between R² and u in this case, the fit being poorer for atoms with large u. For both high- and low-q, there appears to be no relation between the effectiveness of the Gaussian approximation in the examination of individual atomic motions and the instrumental resolution, or, in other words, the timescale of the motions examined. In summary, the results in Fig. 4 indicate that the Gaussian approximation is reasonable over the low-q range examined but not the high-q range, and that variation of the instrumental energy resolution between 20 and 100 µeV has little effect on the accuracy of the approximation.

View larger version (31K): [in this window] [in a new window]   FIGURE 4   Goodness-of-fit (R2) versus long-time u for all hydrogen atoms at 300 K for different instrumental energy resolutions and in different q2 ranges: (a) 100 µeV, 0.01 < q2 < 1.44 Å2; (b) 20 µeV, 0.01 < q2 < 1.44 Å2; (c) 100 µeV, 0.16 < q2 < 23.04 Å2; and (d) 20 µeV, 0.16 < q2 < 23.04 Å2. The u are normalized to u at 80 K. R2 is given by 1   (yi  yi,fit)2/ (yi  )2 where yi refers to ln S(q,  = 0) and yi,fit are the fitted values of ln S(q,  = 0).

Motional heterogeneity

We now examine the effect on u² of motional heterogeneity. In Figs. 5 and 6, u²(t) is shown calculated from the MD simulations for t = 10, 100, and 1000 ps. Also plotted is u²(t)_G, obtained by inserting the individual atom u(t) into Eq. 13, summing to obtain S(q,  = 0) for the whole protein, and then taking the slope of ln S(q,  = 0) vs. q². ln S(q,  = 0) vs. q², thus obtained at 300 K, is shown in Fig. 7.

View larger version (20K): [in this window] [in a new window]   FIGURE 5   Temperature dependence of mean square displacement. The solid lines are u2(t) calculated directly from the MD simulation using Eqs. 15 and 16. The dot-dashed lines are u2G calculated from S(q, ) over the low-q range of 0.01 < q2 < 1.44 Å2. All values are weighted by incoherent scattering lengths and have been normalized to the value at 80 K.

View larger version (20K): [in this window] [in a new window]   FIGURE 6   Temperature dependence of mean square displacement. The solid lines are u2(t) calculated directly from the MD simulation using Eqs. 15 and 16. The dot-dashed lines are u2G calculated from S(q, ) over the high-q range of 0.16 < q2 < 23.04 Å2. All values are weighted by incoherent scattering lengths and have been normalized to the value at 80 K.

View larger version (18K): [in this window] [in a new window]   FIGURE 7   Normalized ln S(q,  = 0) calculated at 300 K using the Gaussian approximation () and the fit to ln S(q,  = 0) in the low-q range of 0.01 < q2 < 1.44 Å2 (- - - -).

In Fig. 5-6, the low-q and high-q ranges are used to obtain u²(t)_G. It was shown in the previous section that the Gaussian approximation is valid for the low-q range, and its effect on u is invariant with the energy resolution (i.e., the timescale of the motions concerned). Therefore, in Fig. 5, in which the same low-q range was used, the use of the Gaussian approximation has little effect on the differences between u²(t)_G and the MD simulation u²(t). Rather, the differences arise from motional heterogeneity. Figure 7 shows that, at 300 K, motional heterogeneity leads to markedly non-Gaussian total scattering even when the individual atom scattering is assumed to be Gaussian. Figure 5 shows that, for the low-q range u²(t)_G underestimates u²(t) by ~10% at t = 10 ps, ~20% at 100 ps, and ~30% at 1000 ps. Consequently, motional heterogeneity has a significant effect on u² that is, in addition, timescale dependentas the timescale of the motions increase the error due to motional heterogeneity also increases.

The question now arises as to whether the timescale dependence of the effect of motional heterogeneity on the derived u² is also reflected in a timescale dependence of the distribution of u(t) in the MD simulations. The distribution of u at t = 10, 100, and 1000 ps, is shown in Fig. 8. The range of u values indeed increases with timescale, from 0.3-10 Ų at 10 ps, through 0.3-13 Ų at 100 ps, to 0.4-19 Ų at 1000 ps. Therefore, at longer times, u are more widely distributed. The distributions are approximately exponential in shape, as shown by the fitted curves in Fig. 8. The data in Fig. 8 are consistent with the conclusion drawn from Fig. 5 because the wider range of u at higher temperatures leads to the approximation made in Eq. 14 becoming less reliable.

View larger version (18K): [in this window] [in a new window]   FIGURE 8   Distribution N(u) of u for hydrogen atoms at 300 K (solid lines). A function of the form y = A exp(xB) was fitted to the distributions. The parameters A and B found from the fitting are: (a) A = 128.7, B = 1.0 for 10 ps; (b) A = 115.3, B = 0.9 for 100 ps; and (c) A = 93.3, B = 0.8 for 1000 ps. The fits to the distributions are shown as dashed lines. The vertical line in (a) indicates u2G = 1.836 Å2, calculated in the low-q range of 0.01 < q2 < 1.44 Å2 at 300 K.

In the high-q analysis (Fig. 6) the error in u² is much larger: ~50% at 10 ps, ~60% at 100 ps, and ~70% at 1000 ps, i.e., u²(t)_G are roughly half as large as those obtained for the low-q range. This increased error arises from the cumulative effects of making the Gaussian approximation and motional heterogeneity.

Experimental data treatment of MD-derived S(q, )

To further our practical understanding of the approximations introduced in elastic neutron-scattering data analysis, in Fig. 9, several methods for deriving the mean square displacements are compared. As reference values, the long-time (t = 1000 ps) MD u², calculated using Eq. 16, and the harmonic mean square displacements, u²_class and u²_qm are given. These are compared with mean square displacements obtained by analyzing S(q, ) from the MD trajectories as if it were experimentally derived, i.e., by subjecting the MD-derived S(q, ) to various data-analysis procedures adopted experimentally. In performing the experimental data treatments, we examine the effect of varying instrumental energy resolution and both effects that have been individually examined in the previous two sections, i.e., motional heterogeneity and the Gaussian approximation. The following experimental data treatments were performed:

1.	u2nores was obtained by calculating S(q, ) from Eqs. 7 and 8 without convolution with an instrumental resolution function, taking the elastic peak intensity S(q,  = 0), and fitting a straight line to ln S(q,  = 0) vs. q2 between 0.01 < q2 < 1.44 Å2. This analysis corresponds to a low-q analysis that would be made experimentally on an instrument of very high energy resolution, i.e., < 0.1 µeV. Because a low-q range is being used, the Gaussian approximation is valid, and motional heterogeneity is the only significant effect here.
2.	u2Elowq was calculated using the same procedure as in 1 above but S(q, ) is convoluted with a 20-µeV full width at half maximum resolution function, and S(q,  = 0) is consequently obtained by integrating S(q, ) over ± = 10 µeV. The low-q range (0.0 < q2 < 1.44 Å2) is again used. In this treatment, motional heterogeneity and the energy-resolution function (timescale dependence) will significantly affect the results, whereas the Gaussian approximation remains valid.
3.	u2Ehighq is calculated using the same procedure as in 2 above but using the high-q range (0.16 < q2 < 23.04 Å2). Motional heterogeneity, the Gaussian approximation, and the effect of limited instrumental energy resolution are all expected to significantly affect the results here.

View larger version (21K): [in this window] [in a new window]   FIGURE 9   Temperature dependence of the mean square displacement for: () long-time MD u2(1 ns); (×) u2nores; (*) u2Elowq; () u2Ehighq; () u2class; and (- - - -) u2qm. All mean square displacements are normalized with respect to the 80-K value.

Figure 9 shows that, above ~100 K the MD-derived long-time u² are significantly larger than u² obtained from any of the experimental data-treatment methods. u²_nores is ~30% lower than the long-time u². The only significant approximation made in the derivation of u²_nores is the neglect of motional heterogeneity. The 30% reduction in u² is in accord with that found in the previous section when comparing u² and u²_G over the same q range. It is also consistent with the observation that u²_G and u²_nores have similar values, although they were calculated using different methods.

In u²_Elowq, the effect of limited energy resolution is added to that of motional heterogeneity. The result is ~50% lower than the long-time u². Comparing u²_Elowq with u²_nores, shows that the 20-µeV energy resolution decreases u² by ~20%.

The q-range used in obtaining u² also has a large impact on the results. A high-q range introduces errors due to making the Gaussian approximation into u²_Ehighq, along with the effects already examined, motional heterogeneity and energy-resolution function. The combination of all three effects reduces u² by ~70%. Indeed, u²_Ehighq is lower than u² obtained with a harmonic model. Comparing u²_Ehighq with u²_Elowq indicates that the Gaussian approximation reduces u² by ~20%.

The various effects examined above all reduce the values of u² obtained from elastic neutron-scattering analyses. To a first approximation, reduction of u² can be considered equivalent to reduction of the timescale at which the motions are observed. Therefore, we determined the effective timescales of the various experimental u² by determining for which value of u²(t) from the simulation is closest to the u² obtained from each of the three data treatments. It was not possible to obtain a meaningful least-squares fit of u²(t) to u²_Ehighq. However, the fits for u²_nores and u²_Elowq were satisfactory (data not shown), leading to best agreement for t = 92 ps and 5 ps, respectively. This means that the effect of motional heterogeneity is equivalent to reducing the effective timescale of the motions examined by ~90% (92 ps compared to the MD simulation reference value of t = 1000 ps), and that of motional heterogeneity combined with the 20-µeV energy resolution function decreases the effective timescale of observable motions by ~99.5% (5 ps compared with 1000 ps).

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS REFERENCES

The present article examines aspects of the temperature dependence of protein dynamics and how information on this can be extracted from elastic neutron scattering. The theoretical model system examined is that of a single BPTI molecule. This is an approximation of a dry protein experimental system that neglects environmental effects such as the possible presence of protein-protein contacts in an experimental system. In the cases of dry proteins examined experimentally so far, to our knowledge, the presence of dynamical transition behavior (from harmonic to anharmonic dynamics with increasing temperature) has not yet been reported. However, the present simulations clearly show the activation of anharmonic motions on a range of timescales and at temperatures (~100-120 K) significantly lower than hitherto reported, also for solvated proteins. Therefore, although ample evidence exists that dynamical transition behavior is strongly coupled to solvent (see e.g., Ferrand et al., 1993; Demmel et al., 1997; Fitter et al., 1998a,b; Fitter, 1999; Vitkup et al., 2000; Réat et al.,