Picosecond fluctuating protein energy landscape mapped by pressure–temperature molecular dynamics simulation
 ^{†}Physical Biology Center for Ultrafast Science and Technology, Arthur Amos Noyes Laboratory of Chemical Physics, California Institute of Technology, 1200 East California Boulevard M/C 127–72, Pasadena, CA 91125;
 ^{‡}Computational Molecular Biophysics, Interdisciplinary Center for Scientific Computing (IWR), University of Heidelberg, Im Neuenheimer Feld 368, D69120 Heidelberg, Germany;
 ^{§}University of Tennessee/Oak Ridge National Laboratory, Center for Molecular Biophysics, Oak Ridge National Laboratory, P.O. Box 2008, Oak Ridge, TN 378316164; and
 ^{¶}Laboratory of Molecular Design, Center for Bioinformatics, Institute of Molecular and Cellular Biosciences, University of Tokyo, Tokyo 1130032, Japan
See allHide authors and affiliations

Contributed by Ahmed H. Zewail, August 30, 2007 (received for review June 28, 2007)
Abstract
Microscopic statistical pressure fluctuations can, in principle, lead to corresponding fluctuations in the shape of a protein energy landscape. To examine this, nanosecond molecular dynamics simulations of lysozyme are performed covering a range of temperatures and pressures. The well known dynamical transition with temperature is found to be pressureindependent, indicating that the effective energy barriers separating conformational substates are not significantly influenced by pressure. In contrast, vibrations within substates stiffen with pressure, due to increased curvature of the local harmonic potential in which the atoms vibrate. The application of pressure is also shown to selectively increase the damping of the anharmonic, lowfrequency collective modes in the protein, leaving the more local modes relatively unaffected. The critical damping frequency, i.e., the frequency at which energy is most efficiently dissipated, increases linearly with pressure. The results suggest that an invariant description of protein energy landscapes should be subsumed by a fluctuating picture and that this may have repercussions in, for example, mechanisms of energy dissipation accompanying functional, structural, and chemical relaxation.
A fundamental concept in protein biophysics is that of an energy landscape. Originally introduced to explain the temperature dependence of the rebinding kinetics of ligands to myoglobin (1, 2), over time the energy landscape picture has been successfully used to describe such diverse phenomena as the temperatureinduced dynamical transition, protein folding, and enzyme function (3–5).
Valuable information on the physics of protein energy landscapes has been obtained through study of the temperature dependence of protein structure and dynamics (6). However, it is only recently that the experimental tools have become available to investigate the effect of pressure on the structure and dynamics of biomolecules (7, 8). Whereas temperature variations influence both the thermal energy and, due to thermal expansivity, the volume of a system, a pressure perturbation affects only the system volume, via the isothermal compressibility. For proteins, for example, pressure changes of a few kbar give rise to changes in free energy that are similar to the freeenergy difference between the native and unfolded state (8).
In Fig. 1, a schematic diagram of a protein energy landscape is shown. The region of minimal free energy, which accomodates the ensemble of native or folded states, is surrounded by regions of high free energy, where the protein is partly or wholly unfolded. Local minima in the energy landscape may give rise to functional or folding intermediate states. In contrast to temperature variations, the application of pressure alters the structure of the energy landscape and may, for example, result in nonnative configurations competing with the ensemble of lowest freeenergy states, thus leading to pressureinduced structural changes or unfolding. In this contribution, however, the focus is on pressureinduced changes in the nativestate region of the energy landscape.
On a microscopic length scale, i.e., the length scale of a single protein molecule, pressure is not constant but fluctuates with amplitudes in the kbarrange around the macroscopic mean value. These fluctuations give rise to instantaneous changes in the energy landscape topology on the characteristic time scale τ_{F} = d/v _{s}, where d is the molecular diameter and v _{s} the velocity of sound; for small globular proteins in solution, realistic values are in the order of d ≈ 1 nm and v _{s} ≈ 10^{3} m s^{−1}, hence τ_{F} ≈ 1 ps. Motions on this time scale are important for overall protein flexibility and, thus, may also be of relevance for functional dynamics.
In this contribution, the pressureinduced changes in the nativestate region of the energy landscape of the protein hen egg white lysozyme (HEWL) in solution are investigated by using molecular dynamics (MD) simulations in the pressure and temperature ranges from 0.1 to 1,000 MPa and from 20 to 320 K, respectively. HEWL has been shown to readily and reversibly denature above 350 MPa on the time scale of minutes (9). However, here the focus is on changes upon pressurization in the local curvature of, and barrier heights separating substates in the energy landscape of the native ensemble, as is illustrated in Fig. 1. Furthermore, the pressuretemperature dependence of dissipative forces acting on collective modes of protein motion are analyzed. The results suggest that the traditional, static view of a protein energy landscape should be subsumed by a fluctuating picture, and that this can potentially have functional consequences.
Results and Discussion
In Fig. 2, we show representative examples of picosecond fluctuations in pressure and atomic coordinates observed during the MD simulations. The pressure of the whole simulationbox rapidly fluctuates with an amplitude of ≈50 MPa around the reference value. The fluctuations in atomic position involve fast components on the subpicosecond time scale, which resemble vibrations in the underlying approximately harmonic potential, and slower components on the picosecond (or longer) time scale, which are due to variations in the position and force constant of that potential. Fig. 2 also illustrates how the underlying potential and force constant can be determined from the fluctuations in atomic position. The variation of the force constant on the picosecond time scale is also shown.
A measure for the overall motion present in a protein molecule is provided by the timedependent atomic meansquare displacement (MSD),
From the slope of 〈Δr ^{2}〉_{P}(T) in the linear regime (T ≤ 160 K) the effective force constant of the underlying harmonic potential can be derived by using k _{eff}(P) = 6k _{B}/(d〈Δr ^{2}〉_{P}/dT) (12). The results are shown in Fig. 4. k _{eff}(P) increases with increasing pressure corresponding to an overall stiffening of the protein. The slope of k _{eff}(P) is approximately linear in two distinct ranges of P, being 3.0 ± 0.2 pN/nm MPa and 1.6 ± 0.4 pN/nm MPa for P ≤ 400 MPa and P ≥ 700 MPa, respectively, with a transition around P* ≈ 480 MPa, indicating a qualitative change in the pressure response of the protein energy landscape.
The transition pressure P* is closely similar to the value of 400 MPa found previously for the pressureinduced dynamical transition in the atomic MSD for a crystalline protein, Staphylococcal nuclease (SNase) (13), between two approximately linear regions of 〈Δr ^{2}〉_{300K}(P), with d〈Δr ^{2}〉(P)/dP_{P<P*} > d〈Δr ^{2}〉(P)/dP_{P>P*}. To investigate whether a similar transition occurs in the present system of HEWL in solution, 〈Δr ^{2}〉_{T} is plotted against pressure in Fig. 5. For all temperatures, 〈Δr ^{2}〉_{T}(P) significantly decreases with increasing pressure but with the slope being nonlinear over the whole pressure range. Assuming this decrease in mobility with increasing P being due to the reduction in accessible volume, 〈Δr ^{2}〉/〈Δr _{0} _{2}〉 = (V/V _{0})^{2/3}, and using the relations β = −(∂(lnV)/∂P) and μ = ∂β^{−1}/∂P, where β is the compressibility and μ is the nonlinearity index, yields where c _{1,2} are integration constants. For lysozyme in solution the pressure dependence of 〈Δr ^{2}〉_{T}(P) can be well described by using Eq. 1 , as shown in Fig. 5. The fit results for μ are also shown in Fig. 5 and vary between 1.1 and 1.7. Experimental μ values for proteins range from 7 to 12 (14), and for an ideal harmonic system μ = 1 (15). This result indicates that lysozyme in the present MD simulations only partially relaxed at high pressure, but was rather elastically compressed and thus remained close to the native conformation.
The above analysis investigated the dynamics of individual atoms, averaged over the whole ensemble present in the protein, thus providing a meanfield measure for the local (cartesian realspace) energy landscape explored by a single atom. The curvature of this local energy landscape was found to increase with increasing pressure, whereas the effective barrier height separating substates remains unchanged. These results are pictorially described in Fig. 1 Upper Right.
In the following, the energy landscape is further examined but with the focus redirected from the above meanfield description to the characterization of collective dynamics. Collective coordinates were determined by using principal component (PC) analysis of the full MD trajectories. Apart from ≈10 modes at T > T _{g}, all PC modes were found to be approximately harmonic. Thus, Brownian dynamics in a harmonic potential was used to model the collective dynamics and, in particular, to determine the eigen and damping frequencies of each PC mode. Results for anharmonic and/or undersampled lowestfrequency modes (≈5–20 modes for each trajectory) were excluded in the following analysis. However, it is important to note that, due to the superposition of PC modes (while sampling mode m, the trajectory simultaneously samples all other modes), the effect of anharmonic (and also solvent) dynamics is implicitly included in all mode variables, i.e., in the time series of the projection σ_{m}(t). In the Brownian description, Eq. 2 , this effect is mediated by the random force, ξ.
In Fig. 6, we show the damping frequencies, Γ plotted against the eigenfrequencies, ω_{0} of the PC modes for selected pressure and temperature values. In general, the lowestfrequency PC modes are those with the largest amplitude. Fig. 6 Top shows the distribution of Γ_{P}(ω_{0}) for individual modes at 300 K. At all pressure values the distribution of damping frequencies has a spread of ≈0.3 THz (at given ω_{0}) and the average Γ increases with increasing ω_{0}. For the lowestfrequency modes Γ > ω_{0} and, therefore, the motion along these modes is overdamped. Γ = ω_{0} is critical damping. Critical damping establishes the fastest process for dissipating the energy contained in an oscillator and, therefore, the determination of the critical damping frequency Γ_{c} is of possible functional interest.
At 300 K, the application of pressure gives rise to an increase of Γ at any given eigenfrequency for modes in the range ω_{0} ≲ 3 THz, leading to a pressureinduced increase in the critical damping frequency, from 0.66 THz at 0.1 MPa to 1.0 THz at 1,000 MPa. The damping frequencies of PC modes with ω_{0} ≳ 3 THz show no significant dependence on pressure, indicating that the dynamics of the largestamplitude modes, i.e., those modes with the lowest eigenfrequencies, is more strongly affected by pressure than the more localized modes with higher frequencies.
The temperature dependence of Γ(ω_{0}) is shown for P = 0.1 MPa in Fig. 6 Bottom. In the regime of harmonic protein dynamics, i.e., for T < T _{g}, Γ(ω_{0}) shows little temperature dependence. However, in the presence of anharmonic dynamics, i.e., for T ≳ T _{g}, Γ(ω_{0}) significantly increases with increasing temperature. This suggests the existence of two distinct damping regimes: one that is associated with the dynamics within harmonic substates, and one of relatively larger strength that is associated with anharmonic and/or diffusive motions. This finding is pictorially described in Fig. 1 Lower Right, where an increased amplitude in the energy landscape roughness symbolizes relatively stronger damping.
Fig. 6 Middle shows a comparison of Γ_{P}(ω_{0}) for two temperatures, 100 K and 300 K. For PC modes with ω_{0} ≳ 1 THz, the main effect of the lower temperature is to reduce the damping frequency at any given ω_{0} by ≈0.4 THz. This is in accord with the above finding. However, for the lowestfrequency (ω_{0} ≲ 0.5 THz) modes, the temperatureinduced change in Γ is smaller. Because these modes describe collective displacements distributed over the whole protein, damping of these modes may be of particular importance. At temperatures below T _{g}, the surrounding solvent acts as a cage for protein dynamics (on the ns/ps time scale) that may provide a strong damping mechanism additional to proteininternal damping. With increasing frequency, the PC modes become more localized and, thus, the relative importance of protein/solvent interactions is reduced and proteininternal damping mechanisms become more dominant. This hypothesis is supported by the observation that, at 100 K, Γ(ω_{0}) initially decreases with increasing ω_{0}, as is also shown in Fig. 6. For T > T _{g}, with increasing temperature, this initial decrease in Γ(ω_{0}) becomes smaller and vanishes at ≈260 K, as can be seen for P = 0.1 MPa in Fig. 6 Bottom.
The pressure and temperature dependence of Γ(ω) of the lowfrequency PC modes also determines the critical damping frequency, Γ_{c}(P,T). Γ_{c}(P,T) was determined for all temperatures and pressures as the intersection of the Γ = ω_{0} line with a linear (T ≥ 260 K) or quadratic (T ≤ 240 K) fit to Γ(ω_{0}) for ω_{0} ≤ 1.5 THz. The results are shown in Fig. 7. At constant temperature, Γ_{c}(P) increases approximately linearly with increasing pressure. Because Γ(ω) is a slowly varying monotonic function for ω ≲ 1 THz, Γ_{c} ∝ P is equivalent to an approximately linear increase with pressure of Γ, which can also be seen in Fig. 6 Top and Middle. For T ≤ 220 K and T ≥ 240 K, the slope Γ′_{c} of Γ_{c}(P) is almost independent of temperature, the averages being (2.7 ± 0.1)10^{−4} THz MPa^{−1} and (3.56 ± 0.06)10^{−4} THz MPa^{−1}, respectively. This difference in slope can be related to the two damping regimes found above, which are associated with the harmonic and anharmonic dynamics below and above T _{g}, respectively.
The critical damping frequency at ambient pressure (i.e., at 0.1 MPa), Γ_{c} ^{0}, was also determined from the linear fits to Γ _{c} (P). For T ≲ 160 K, Γ_{c} ^{0}(T) fluctuates around a constant value of ≈0.43 ± 0.01 THz, whereas for T ≳ 180 K Γ_{c} ^{0}(T) increases linearly with increasing temperature. Again, this finding is consistent with the above observation that the damping is independent of temperature in the harmonic substates. However, when the temperature increases to or above T _{g}, the trajectories cumulatively sample more anharmonic regions of the energy landscape, in which Γ is larger and thus Γ_{c} increases.
Conclusions
Understanding the topology of and the dynamics on protein energy landscapes furnishes insight into the physical mechanisms governing protein folding and function. Here, the variation of pressure and temperature in MD simulations of lysozyme has revealed a pressureinduced increase in the curvature of the harmonic substates lining the ensemble of native states. The effective force constant associated with this curvature initially increases linearly with increasing pressure, but at ≈480 MPa exhibits a 2fold reduction in slope, indicating the presence of two pressure–response mechanisms. Because the application of pressure in the 100 MPa range does not significantly change the lengths of chemical bonds, both mechanisms are likely to originate from interactions governed by van der Waals forces. The question arises as to whether these mechanisms describe the different pressure responses of the protein and the solvent, or whether both can be attributed to either the protein or the solvent. For example, the application of pressure significantly decreases the tetrahedrality of the water structure (16). Because the protein dynamics is slaved to the solvent fluctuations (17), the comparatively more rigid tetrahedral water structure at lower pressure may then cause a stiffer protein response to pressure changes than the less rigid van der Waals liquidlike water structure at higher pressure.
Protein relaxation accompanying functional processes such as folding, ligand binding or chemical reaction often requires the dissipation of kinetic and/or potential energy contained within collective degrees of freedom (modes). This dissipation is fastest for modes that are critically damped and, because of the equipartition of energy among all modes and mode coupling (18), these modes contribute disproportionately high to energy dissipation and thus relaxation. Accordingly, the critical damping frequencies for collective modes of protein motions were determined. Γ_{c}(P) is found to depend approximately linearly on pressure indicating that the pressure dependence is governed by a single physical mechanism. In contrast, the temperature dependence of Γ_{c}(T) reveals the existence of two damping mechanisms with a transition temperature close to the protein dynamical transition. One temperatureindependent mechanism (present at all T) is associated with the motions within harmonic substates. The other mechanism (present for T ≳ T _{g}) originates from anharmonic and/or diffusive dynamics connecting the harmonic substates and its contribution increases with increasing temperature.
Finally, the above results suggest that an invariant description of the protein energy landscape is incomplete on the picosecond time scale. Fluctuations in volume, and thus in pressure, are inevitably present on the molecular length scale and lead to variations in both the topology and the dissipation characteristics of the energy landscape. For globular proteins (with an average size on the nanometer length scale), these fluctuations occur on the picosecond time scale, which is the predominant time scale for global vibrations (19). Therefore, pressurefluctuation induced variations in the protein energy landscape may be important for substate alterations of functional relevance and may expedite relaxation processes by broadening the range of the most effective, i.e., critical, damping.
Methods
Molecular Dynamics.
MD simulations were performed of hen egg white lysozyme over the temperature and pressure ranges 20–320 K (every 20 K) and 0.1–1,000 MPa (every 50 MPa). The initial structure 1GXV (20) (solved at atmospheric pressure using NMR) was taken from the Protein Data Bank (21) and centered in the rhombicdodecahedron primary simulation cell with initial box length of 70 Å. A total of 7,309 TIP4P water molecules (22), 15 sodium ions, and 24 chloride ions were added as a 0.1 M salinity solvent, resulting in an electrically neutral system comprising 31,236 atoms.
MD simulations were performed by using the Gromacs suite of programs (23, 24) with the allatom OPLSAA/L force field (25) and periodic boundary conditions. The OPLS force field was parameterized in conjunction with the TIP4P potential, which reproduces the anomalous diffusion coefficient of water in the present temperature and pressure ranges (16). Electrostatic interactions were computed by using the particle mesh Ewald method (26) with the directsum cutoff and Fourier grid spacing being 9 and 1.2 Å, respectively, and van der Waals interactions were cut off at 14 Å.
The system was energy minimized to a rootmeansquare (RMS) force gradient of 1.6 × 10^{−2} kJ mol^{−1}Å^{−1}, subsequently heated during 100 ps and then pressurized and equilibrated during 600 ps. Temperature and pressure coupling were enforced by using the extendedensemble Nosé–Hoover/Parrinello–Rahman algorithms (27–30) with a coupling time constant of 1 ps. The production phase at each temperature–pressure value was 1 ns long. Coordinates were saved with a sampling interval of 0.1 ps used in all subsequent analyses. For all simulations, the trajectoryaverage protein C _{α}atom RMSdeviation from the experimental structure was <1.7 Å, indicating that the protein structure remained close to the native state in all simulations. Also, the comparison with an additional 10 ns trajectory (at T = 300 K, P = 0.1 MPa) showed that the 1ns sampling of the native state is sufficient for the atomic RMS displacements to converge within 5% of the average value obtained during the 10ns trajectory.
Eigen and Damping Frequencies of Protein Collective Motions.
Collective protein motions present in each trajectory r _{T,P}(t) were determined by using principal component (PC) analysis as described in ref. 31, i.e., by diagonalizing the allproteinatom massweighted variance–covariance matrix yielding the eigenvectors v _{m}. The dynamics along the mth PC mode were investigated by analyzing the time series of the projection, σ_{m}(t) = (M ^{1/2}[r _{T, P}(t) − r _{T, P} ^{ave}]) · v _{m}, where M is the diagonal mass matrix and r _{T, P} ^{ave} denotes the average structure. The effective freeenergy landscape along the mode is given by G _{m}(σ) = −k _{B} T ln p _{σ} _{m}, where k _{B} is the Boltzmann constant and p_{σ} _{m} dσ is the probability of σ_{m}(t) ∈ [σ,σ + dσ). If the motion is harmonic, then G _{m}(σ) is also harmonic and the probability density p_{σm} is Gaussian. Anharmonic and/or quasiharmonic modes were found for T ≳ 200 K at all pressure values as has been reported for other systems (11, 13).
The projection σ_{m}(t) can be rationalized as a stochastic process sampling the underlying potential. Therefore, if G _{m}(σ) is approximately harmonic, σ_{m}(t) can be described as Brownian motion in a harmonic potential using the set of Langevin equations where Γ is the damping frequency, ω_{0} is the eigenfrequency, and ξ(t) is a random force that represent the effect of environmental degrees of freedom, i.e., ξ_{m}(t) implicitly represents the effect on σ_{m}(t) of the solvent dynamics but also of σ_{m′≠m}(t). In equilibrium, ξ and Γ are related by the fluctuation–dissipation theorem, 〈ξ(t)ξ(t′)〉 = 4Γk _{B} T δ(t − t′), where δ(t) is the Dirac Delta function.
The power spectral density (PSD) S_{υυ}(ω) can be calculated analytically (32)
with
Acknowledgments
This work was supported by the Physical Biology Center for Ultrafast Science and Technology. L.M. acknowledges support from a Deutsche Physikalische Gesellschaft/Deutsche Forschungsgemeinschaft fellowship to visit Japan before joining the group at The California Institute of Technology.
Footnotes
 ^{‖}To whom correspondence should be addressed. Email: zewail{at}caltech.edu

Author contributions: L.M., J.C.S., A.K., and A.H.Z. performed research; and L.M., J.C.S., A.K., and A.H.Z. wrote the paper.

The authors declare no conflict of interest.
 © 2007 by The National Academy of Sciences of the USA
References
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Gruner SM
 Katrusiak A ,
 McMillan P
 ↵
 ↵
 ↵

↵
 Zaccaï G
 ↵
 ↵

↵
 MoelwynHughes EA
 ↵

↵
 Fenimore PW ,
 Frauenfelder H ,
 McMahon BH ,
 Parak FG
 ↵
 ↵
 ↵

↵
 Berman HM ,
 Westbrook J ,
 Feng Z ,
 Gilliland G ,
 Bhat TN ,
 Weissig H ,
 Shindyalov IN ,
 Bourne PE
 ↵
 ↵

↵
 Lindahl E ,
 Hess B ,
 van derSpoel D
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Reichl LE
 ↵
 ↵