New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Vibrational energy relaxation in proteins

Edited by Bruce J. Berne, Columbia University, New York, NY, and approved January 12, 2005 (received for review December 6, 2004)
Abstract
An overview of theories related to vibrational energy relaxation (VER) in proteins is presented. VER of a selected mode in cytochrome c is studied by using two theoretical approaches. One approach is the equilibrium simulation approach with quantum correction factors, and the other is the reduced model approach, which describes the protein as an ensemble of normal modes interacting through nonlinear coupling elements. Both methods result in similar estimates of the VER time (subpicoseconds) for a CD stretching mode in the protein at room temperature. The theoretical predictions are in accord with previous experimental data. A perspective on directions for the detailed study of time scales and mechanisms of VER in proteins is presented.
When a protein is excited by ligand binding, ATP attachment, or laser pulses, vibrational energy relaxation (VER) occurs. Energy initially “injected” into a localized region flows to the rest of the protein and surrounding solvent. VER in large molecules (including proteins) is an important problem for chemical physics (1, 2). Even more significant is the challenge to relate VER to fundamental reaction processes, such as a conformational change or electron transfer of a protein, associated with protein functions. The development of an accurate understanding of VER in proteins is an essential step toward the goal of controlling protein dynamics (3).
Because of the advance of laser technology, there have been many experimental studies of VER in proteins (4–17). These experimental works are impressive, but it is difficult to derive detailed information from the experimental data alone. Theoretical approaches, including atomicscale simulations, can provide more detailed information. In turn, experimental data can be used to refine simulation methods and empirical force fields. This combination of experimental and theoretical studies of protein structures and dynamics has begun to blossom. As experimental methods develop further and theoretical approaches grow in accuracy, the relationship will become fruitful.
There have been many theoretical tools (see Theories) developed to analyze VER in proteins. Some aspects of VER in proteins can be explained by perturbative formulas based on the equilibrium condition of the bath (see Cyt c), but the use of the perturbative formulas may be too restrictive to describe protein dynamics generally at room temperature. In this article, we not only discuss the success of such established methods but also present a perspective on the study of VER in proteins.
Theories
In this section, we present a selective overview of theories appropriate for the study of VER in proteins. For the most part, these theories have been developed to deal with VER in liquids, solids, or glasses. For recent reviews, see refs. 18–20. We refer to two distinct categories; one is based on equilibrium dynamics and Fermi's golden rule, whereas the other is based on nonequilibrium dynamical models.
Fermi's Golden Rule. If (i) there is a clear separation between the system and bath, (ii) the coupling between them is weak enough, and (iii) the bath is assumed to be at thermal equilibrium, we can use quantum mechanical perturbation theory to derive a VER rate 1/T _{1} through Fermi's golden rule (19, 20), where the force–force correlation function ζ_{qm}(t) is defined as follows: where is the quantum mechanical force applied to the relaxing bond (system) considered, m_{S} is the system mass, ω _{S} is the system frequency, β is an inverse temperature, and the bracket in Eq. 2 indicates a quantum mechanical average.
However, this time correlation function is very hard to calculate numerically. As a result, many approximate schemes have been proposed to address this limitation. Some of the most successful approaches are mentioned below.
Landau–Teller–Zwanzig (LTZ) formula. The most simple approximation is to take the classical limit (ℏ → 0) of Eq. 1, Here, the bracket denotes a classical ensemble average. This equation is called the LTZ formula, which has been applied to the study of VER in liquids (21). This strategy was used by Sagnella and J.E.S. to discuss the VER of CO in photolyzed CO myoglobin (Mb*CO) (22). This approximation should be good for lowfrequency modes, but it becomes questionable for highfrequency modes because of quantum effects. Thus, advanced methods have been proposed to address this deficiency of the LTZ formula.
Quantum correction factor (QCF). The first alternative to the LTZ formula is the QCF method. The basic idea of the QCF method is to relate a quantum mechanical correlation function with its classical analog (23). When this procedure is done for the force autocorrelation function in Eq. 1, the final expression for the VER rate is as follows: where Q(ω _{S} ) is the QCF for the considered VER process and Q_{H} (ω _{S} ) is the QCF for a onephonon relaxation process (harmonic QCF), In refs. 19 and 20, this result was expressed as , which is correct in the limit βℏω _{S} » 1, as was appropriate for those studies.
If the relaxation process is the linear resonance (1:1 Fermi resonance), then Q(ω _{S} ) = Q_{H} (ω _{S} ); i.e., (24). Skinner and coworkers provided a theoretical framework for organizing and expanding on various QCFs appropriate for specific dynamical processes, depending on the underlying mechanism of VER. Although this strategy has been criticized (25–27), it is known that the QCF method works well for specific problems (19, 20, 28).
Reduced model approach. An alternative approach to address the shortcomings of the LTZ formula is to use the reduced model approach (18–20), which exploits a normal mode picture of a protein. By representing the Hamiltonian in terms of system, bath, and interaction terms, the residual interaction term may be expanded perturbatively as follows: Calculating the force from this Hamiltonian and substituting it into Fermi's golden rule in Eq. 1, we can derive a lowestorder VER rate as (19, 20), where In refs. 19 and 20, m_{S} in the perturbative formulas should read m_{S} = 1 because massweighted coordinates were employed.
The original formula contains delta functions, and we have included a width parameter γ to broaden the delta functions for numerical calculations. Another well known formula to describe a VER rate exists, the Maradudin–Fein formula (18, 29). with a width parameter γ. These two formulas are numerically similar for small γ, and equivalent for the limit of γ → 0 (30). This treatment is quantum mechanically exact given the approximate truncated form of the interaction Hamiltonian. We have found that the truncation error (the contribution from higherorder terms) can be a serious problem, especially for proteins. For a more accurate treatment of VER, we must appeal to more advanced methods, as described below.
Other (advanced) approaches. Methods that complement the three methods described above involve calculating the force autocorrelation function ζ(t) appearing in Fermi's golden rule using different levels of approximations. Shi and Geva (27) used a semiclassical approximation (31) for ζ(t) and showed that even the slow relaxation of neat liquid oxygen (at 77 K) can be well reproduced by their method. From their study, it was shown that the short time dynamics of ζ(t) is important to predict the correct VER rate. This result implies that the short time approximation may be adequate for an accurate description of VER. Various timedependent selfconsistent field methods (32) or path integral methods (33) should be applicable to calculate ζ(t). For other methods, see refs. 34–36.
To derive Fermi's golden rule, we have used the Bader–Berne correction (24), which holds only for harmonic systems. Bader et al. extended this to an anharmonic system within a classical framework (37) and found that the VER of such a system can be nonexponential in time and is significantly affected by the character of the bath. This consideration is important when one studies the VER of CO in MbCO, especially for the VER of a highly excited CO bond.
Nonequilibrium Simulation. The equilibrium simulation methods described above are based on Fermi's golden rule and invoke several assumptions as described above. These assumptions might be invalid in some cases. Because VER is a nonequilibrium phenomenon, the appeal of nonequilibrium approaches is quite natural.
Classical approaches. Classical nonequilibrium simulations to investigate VER in proteins were first conducted by Henry et al. (38). In conjunction with their experimental studies, they employed classical molecular dynamics simulations of heme cooling in Mb and cytochrome c (Cyt c) in vacuum and found that heme cooling occurred on two time scales: short (1–4 ps) and long (20 ps for Mb and 40 ps for Cyt c). Nagaoka and coworkers carried out the similar simulations for Mb in vacuum and obtained similar time scales (39). Importantly, they found that the normal mode frequencies localized in the proprionate side chains of the heme are resonant with the water vibrational frequencies.
The J.E.S. group executed several numerical simulations for Mb and Cyt c in water. Sagnella and J.E.S. showed that the VER for Mb in water can be described by a single exponential with a VER time of a few picoseconds (40). Furthermore, they suggested that the main doorway of VER is due to the coupling between the proprionate side chains and water, which is in accord with Nagaoka's and Hochstrasser's observations. Bu and J.E.S. supported this view through simulations of mutant Mbs and Mb variants having structurally modified heme groups (41). They also investigated VER of Cyt c in water and found that the VER presents a biphasic exponential decay with two VER times: fast (a few picoseconds) and slow (tens of picoseconds) (42).
Kidera's group studied VER in proteins from a different perspective (43). They excited a single normal mode in Mb and examined the vibrational energy transfer (VET) between normal modes. As is well known, VET is caused by (nonlinear) Fermi resonance: if the frequency matching is good, and the coupling between normal modes is strong enough, there will be VET. This picture is very useful to characterize VET at low temperatures. However, at high temperature, nonresonant VET occurs. They numerically found that the amount of VET is proportional to a reduced model energy including up to thirdorder coupling elements (see also Reduced model approach).
Quantum approaches. For all but the simplest systems, quantum approaches for nonequilibrium simulations are approximate and timeconsuming. Nevertheless, these methods can overcome problems inherent to classical simulations. There are two categories: vibrationally quantum methods and electronically quantum methods.
Hahn and Stock used a reduced model (consisting of the retinal rotation and other environmental degrees of freedom) to describe the pump–probe spectroscopy for the retinal chromophore in rhodopsin (44). Flores and Batista, employing the same model, suggested the possibility of controlling the retinal rotation by two (chirped) laser pulses (45). To solve the quantum dynamics for the large system, they employed timedependent selfconsistent field methods (32). Notably, vibrational SCF methods have been used to calculate vibrational energy levels for a small protein [bovine pancreatic trypsin inhibitor (BPTI)] (46).
The combination of classical simulations for vibrational motions and quantum calculations for electronic structure, in some portion of the molecule, has been widely used for the calculations of up to moderatesized molecules. One cuttingedge application to a large system is the calculation of the photoisomerization of bacteriorhodopsin in the excited chromophore state by Hayashi et al. (47). In their treatment, a portion of the retinal chromophore including three double bonds was treated as the quantum mechanical region, and the complement, including the protein and water, was treated as the molecular mechanical region. During the simulations, nonadiabatic transitions occur between two electronic states (S_{0} and S_{1}), which was treated semiclassically. They showed numerically that only one bond (C_{13} C_{14}) rotates unidirectionally because of the coupling with the protein, and they found that several other bonds can twist in any direction if there is no protein.
Cyt c
In this section, we focus on one protein, Cyt c and review the recent theoretical studies about this protein. There are several reasons to select this protein as a prototypical one. (i) Cyt c is a relatively small protein with 1,745 atoms. Other proteins of a similar scale are Mb, bovine pancreatic trypsin inhibitor (BPTI), human lysozyme, etc. (ii) The detailed xray structure is known for Cyt c. (iii) nCyt c has a function of electron transfer. The basic theoretical and computational works on Cyt c were summarized by Wolynes and coworkers (48). Wang et al. studied VER in Cyt c using their hydrodynamical method (49). Garcia and Hummer found anomalous diffusion for some principal components of Cyt c in water (50). Here, we describe the results of our studies on the VER of Cyt c by using two different methods, the QCF method described in Quantum correction factor and the reduced model approach described in Reduced model approach, and we compare them with the experimental results of Romesberg and coworkers (12, 13, 17).
Quantum Correction Factor Approach for Cyt c. Bu and J.E.S. (52) employed the QCF approach (see Quantum correction factor) to estimate the VER rate of a CD bond in the terminal methyl group of Met80 in Cyt c (see Fig. 1). Their calculations were carried out by using the program charmm (53), and Cyt c was surrounded by water molecules at 300 K. In Fig. 2, we show the force autocorrelation function and its power spectrum. With the CD bond frequency ω _{S} = 2,133 cm^{–1}, we find , so that the classical VER time is 1.0 ∼ 2.5 ps.
Because the CD bond frequency is located in a transparent region of the vibrational density of states, with no other state overlapping with this frequency (52), we concluded that there is no linear resonance (1:1 resonance). Thus, to use the QCF method, we need to assume nonlinear resonances corresponding to multiphonon VER processes. If the VER process assumes that two lowerfrequency bath modes, having frequencies ω _{A} and ω _{S} – ω _{A} , are each excited by one quantum of vibrational energy, the appropriate QCF (harmonic–harmonic QCF) is (23), Alternatively, if the VER process is one that leads to the excitation of one bath vibrational mode of frequency ω _{A} , with the remaining energy ℏ(ω _{S} – ω _{A} ) being transferred to lower frequency bath rotational and translational modes, the appropriate QCF (harmonic–harmonic–Schofield QCF) is (23), We needed to determine the value of ω _{A} to use these formulas. From the normal mode and anharmonic coefficient calculations carried out in Reduced model approach for Cyt c, we found that the CD mode is strongly resonant with two lowerfrequency modes, ω_{1655} = 685.48 cm^{–1} and ω_{3823} = 1,443.54 cm^{–1}, where ω _{S} – ω_{1655} – ω_{3823} = 0.03 cm^{–1} for the original parameters of charmm. We might be able to choose ω _{A} = 1,443.54 cm^{–1}1 or 685.48 cm^{–}. If we choose ω _{A} = 1,443.54 cm^{–} at 300 K, for the harmonic–harmonic QCF and for the harmonic–harmonic–Schofield QCF. Thus, we find .
Reduced Model Approach for Cyt c. H.F. and colleagues (19, 20) took the reduced model approach (see Reduced model approach) to investigate the VER for the same CD bond stretching in Cyt c. However, in their calculation, all modes represent normal modes, so the CD “bond” turned out to be the CD “mode.” Using the formulas given in Reduced model approach, they calculated the VER rate for the CD mode (ω _{CD} = 2,129.1 cm^{–1}) and other lowfrequency modes (ω_{3,330} = 1,330.9 cm^{–1}, ω_{1,996} = 829.9 cm^{–1}, ω_{1,655} = 685.5 cm^{–1}) as a function of the width parameter γ (Fig. 3). To this end, they needed to calculate anharmonic coupling coefficients according to the following formula: where U_{ik} is an orthogonal matrix that diagonalizes the (massweighted) hessian matrix at the mechanically stable structure K_{ij} , and K_{ij} (±Δq_{S} ) is a hessian matrix calculated at a shifted structure along the direction of a selected mode with a shift ±Δq_{S} .
If we take γ ≃ 3 cm^{–1} (20), we have T _{1} ≃ 0.1 ps, which agrees with the subpicosecond time scale for VER predicted using the QCF method (). We also see that the lowfrequency modes have longer VER time, a few picoseconds, which agrees with the similar calculations by Leitner's group (18). In Fig. 3B , we show the temperature dependence of the VER rate. At low temperatures, the VER rate becomes flat as a function of temperature. At these lower temperatures, the VER is caused by the remaining quantum fluctuation associated with zero point energy.
Related Experiment. Here, we discuss the related experiment by Romeberg and coworkers (12, 13, 17). They measured the shifts and widths of the absorption spectra for different forms of Cyt c; the widths of the spectra [full width at half maximum (FWHM)] were found to be Δω_{FWHM} ≃ 6.0 ∼ 13.0 cm^{–1}. If we can neglect inhomogeneous effects, the estimate of the VER time becomes, which corresponds to T _{1} ≃ 0.4 ∼ 0.9 ps. This estimate is similar to the QCF prediction by using Eq. 4 (≃ 0.3 ∼ 1.0 ps) and larger than the estimate by the reduced model approach by using Eqs. 11 or 15 (≃0.1 ps). This conclusion might be due to the strong resonance between the three modes (4,357, 3,823, and 1,655), which forms a peak near γ ≃ 0.03 cm^{–1} (Fig. 3A ). This resonance causes an increase in the VER rate, so we can say that this estimate of the VER rate is too large. However, there is no peak for the lowfrequency modes for γ < 10 cm^{–1}; the estimate of the VER rate does not seem to be affected by the resonances.
Note also that Romesberg and coworkers studied Met803D, methionine with three deuteriums on the terminal methyl group, whereas we examined Met801D, with one deuterium. It is known that the original charmm force field does not give an accurate frequency corresponding to such an absorption peak. However, the density functional theory calculation for the methionine leads to much better results (Matt Cremeens, personal communication). Clearly, we must improve our force field parameters according to density functional theory or other ab initio methods and examine how further optimization of the parameters affects the resonance structures and the VER rate of the protein.
Concluding Remarks
In this article, we have described theoretical (Theories) approaches to the study of VER in proteins. We have examined VER of a CD stretching bond (mode) in Cyt c from the QCF approach (Quantum correction factor) and the reduced model approach (Reduced model approach). For the CD mode in Cyt c (in vacuum) at room temperature, both approaches yield similar results for the VER rate, which is also very similar to an estimate derived from an experiment by Romesberg and coworkers. Our work demonstrates both the feasibility and accuracy of a number of theoretical approaches to estimate VER rates of selected modes in proteins.
There are advantages and disadvantages of the QCF approach and reduced model approach to the prediction of VER rates in proteins. The QCF method is simple and applicable even for a large molecule like a protein. However, the VER mechanism may not be known apriori, and it must be supplemented by other methods such as anharmonic coefficient calculations. Furthermore, the method relies on the local mode picture, which is easily applicable for highfrequency (localized) modes but not for lowfrequency (delocalized) modes. The reduced model approach is quantum mechanically exact and easily applicable for VER of lowfrequency modes. However, the anharmonic coefficient calculation is cumbersome, even for the thirdorder coupling terms in Cyt c. Moreover, such a Talyor series expansion has not been shown to converge at low order for general systems (26). Our preliminary calculations show that the classical VER dynamics using an isolated methionine does not seem to be affected by including the fourthorder coupling elements (Fig. 4A ), but we need to examine this issue further with quantum mechanical approaches, such as timedependent selfconsistent field methods. There is also an unsolved problem of the width parameter. Actually, this problem is not peculiar to the reduced model approach. The introduction of the width corresponds to coarsegraining, which also appears in the QCF approach when one averages the power spectrum of the force autocorrelation function. The most “ab initio” approach to solve this problem is a rigorous quantum mechanical treatment of the tier structure of energy levels in the protein (54). The other appealing way is to regard γ as a hopping rate between conformational substates (55, 56), or as the inverse of a frequency correlation time, that may be derived from estimates of the frequency fluctuation (57–59).
Because both the QCF and reduced model approaches are based on Fermi's golden rule, there is a limitation for the strength of the interaction between the system and bath. As noted in ref. 20, it might be critical to apply the perturbation method to this situation of VER of the CD mode in Cyt c. Other methods without this deficiency are needed. Promising approaches include nonequilibrium molecular dynamics methods (60), timedependent selfconsistent field methods (32), mixed quantum–classical methods (61), and semiclassical methods (27). Another important issue is to calculate not only VER rates but the physical observables related to the experiment data, such as absorption spectra or 2DIR signals (62, 63). In this case, we also need to deal with the effects of dephasing (decoherence) as well as those of VER.
The accuracy of the force field parameters is the most annoying problem. Our preliminary calculations show that the VER rate in Cyt c can vary by two orders of magnitude when we change the bond force constant by 10% (Fig. 4B ). This situation is similar to that of a reaction rate calculation, in which one must determine the activation energy accurately. Any inaccuracy in the activation energy causes an exponentially large deviation in the rate constant. This problem would be solved by means of ab initio quantum dynamics (Quantum approaches) or the reparametrization of the force field using experimental data or accurate ab initio calculations. Given sufficient accuracy in the force field, we will be in a position to discuss the relation between VER and functions in a protein. As is well known, the dynamics of proteins related to functions are well described by large amplitude (and low frequency) principal components (64, 65). The connection between principal components and VER should be investigated. The ergodic measure (66) would be a good device to examine this issue. As suggested by refs. 4 and 6, collective motions in proteins can be important for the fast VER in proteins. The collective motions near the protein surface including solvation dynamics of water (67–71) might be relevant for the VER and functions. Cyt c and Mb remain excellent target proteins to investigate these fundamental issues of protein dynamics and its relation to functions.
Acknowledgments
We thank Dr. Lintao Bu for collaboration and Dr. Matt Cremeens for density functional theory calculations. This work was supported by National Science Foundation Grant CHE0316551 and the Boston University Center for Computational Science.
Footnotes

↵ † To whom correspondence should be addressed. Email: straub{at}bu.edu.

Author contributions: J.E.S. designed research; H.F. performed research; and H.F. and J.E.S. wrote the paper.

This paper was submitted directly (Track II) to the PNAS office.

Abbreviations: VER, vibrational energy relaxation; LTZ, Landau–Teller–Zwanzig; Mb, myoglobin; Cyt c, cytochrome c; QCF, quantum correction factor.
 Copyright © 2005, The National Academy of Sciences
References
 ↵
 ↵

↵
Dian, B. C., Longarte, A. & Zwier, T. S. (2002) Science 296 , 2369–2373. pmid:12029064
 ↵

↵
Mizutani, Y. & Kitagawa, T. (1997) Science 278 , 443–446. pmid:9334299

Sagnella, D. E., Straub, J. E., Jackson, T. A., Lim, M. & Anfinrud, P. A. (1999) Proc. Natl. Acad. Sci. USA 96 , 14324–14329. pmid:10588704
 ↵
 ↵

Xie, A., van der Meer, A. F. G. & Austin, R. H. (2002) Phys. Rev. Lett. 88 , 01810214. pmid:11800992
 ↵

↵
Leitner, D. M. (2005) Adv. Chem. Phys., in press.

↵
Fujisaki, H., Bu, L. & Straub, J. E. (2005) Adv. Chem. Phys., in press.

↵
Fujisaki, H., Bu, L. & Straub, J. E. (2005) in Normal Mode Analysis: Theory and Applications to Biological and Chemical Systems, eds. Cui, Q. & Bahar, I., in press.
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
Henry, E. R., Eaton, W. A. & Hochstrasser, R. M. (1986) Proc. Natl. Acad. Sci. USA 83 , 8982–8986. pmid:3024159
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
Brooks, B. R., Bruccoleri, R. E., Olafson, B. D., States, D. J., Swaminathan, S. & Karplus, M. (1983) J. Comp. Chem. 4 , 187–217.
 ↵
 ↵
 ↵
 ↵

Ma, J., Huo, S. & Straub, J. E. (1997) J. Am. Chem. Soc. 119 , 6454.
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

Higo, J., Sasai, M., Shirai, H., Nakamura, H. & Kugimiya, T. (2001) Proc. Natl. Acad. Sci. USA 98 , 5961–5964. pmid:11344268

Trek, M. & Tobias, D. J. (2002) Phys. Rev. Lett. 88 , 13810114. pmid:11955127

Fenimore, P. W., Frauenfelder, H., McMahon, B. H. & Parak, F. G. (2002) Proc. Natl. Acad. Sci. USA 99 , 16047–16051. pmid:12444262
 ↵
Citation Manager Formats
More Articles of This Classification
Biological Sciences
Related Content
 No related articles found.