Quantum localization and protein-assisted vibrational energy flow in cofactors

. Quantum effects inﬂuence vibrational dynamics and energy ﬂow in biomolecules, which play a central role in biomolecule function, including control of reaction kinetics. Lifetimes of many vibrational modes of proteins and their temperature dependence, as determined by quantum golden-rule-based calculations, exhibit trends consistent with experimental observation and distinct from estimates based on classical modeling. Particularly notable are quantum coherence effects that give rise to localization of vibrational states of sizable organic molecules in the gas phase. Even when such a molecule, for instance a cofactor, is embedded in a protein, remnants of quantum localization survive that inﬂuence vibrational energy ﬂow and its dependence on temperature. We discuss these effects on the mode-damping rates of a cofactor embedded in a protein, using the green ﬂuorescent protein chromophore as a speciﬁc example. We ﬁnd that for cofactors of this size embedded in their protein and solvent environment at room temperature a golden-rule calculation often overestimates the mode-damping rate.


Introduction
Vibrational energy flow in molecules mediates rates of chemical reactions, influencing the kinetics of conformational isomerization [1]- [6], bond breaking [7,8] and charge transfer [9,10]. There is currently considerable interest in the role of vibrational dynamics and energy flow in protein reaction kinetics [11]- [14], including enzyme kinetics [15], allostery [16]- [18], photochemical isomerization [19] and excited state proton transfer (ESPT) [20,21]; methods to search for and identify energy transport pathways through biomolecules are rapidly being developed [22]- [26]. A broad understanding of the nature of energy flow through different parts of a protein can be combined with ultrafast spectroscopic techniques to reveal changes in protein structure on fast time scales [27]- [33]. Ultrafast photochemical reactions in proteins such as rhodopsin [34], bacteriorhodopsin [35] and photoactive yellow protein [19,36] (PYP) occur from vibrationally unrelaxed states where impulsively excited vibrations appear to be coupled to the reaction coordinate. Coherent vibrations are observed during the primary events of vision [37] and the PYP photocycle [19,38,39]. Recent ultrafast Raman studies [40] of green fluorescent protein (GFP) have revealed that impulsive excitation of a single low-frequency wagging mode of the chromophore assists ESPT. Excess vibrational energy appears to generally facilitate these photochemical reactions in proteins rather than immediately dissipate in a random fashion.
Quantum mechanical effects on vibrational dynamics and energy flow in molecules are critical. Most vibrational modes of organic molecules (even as large as proteins) are unexcited at room temperature. The dynamics of many low-frequency vibrations that extend over large biomolecules are strongly influenced by the molecular geometry [41]- [46] and can be described reasonably well by beginning with a classical perspective. However, quantum effects play an important role in determining the rate of energy flow among the large majority of protein vibrations. Quantum coherence can give rise to localization of vibrational states of sizable organic molecules [47]- [52], including vibrationally cold biomolecules [53,54], which in some cases exhibit mode selective chemistry [55,56]. Although quantum localization of vibrational states does not strictly occur in large biomolecules such as proteins at physiological temperatures, as discussed below, coherence giving rise to localization can still persist in parts of the molecule and in cofactors in the presence of coupling to the protein and solvent environment and influence vibrational mode lifetimes and energy flow. 3 In this paper we examine quantum localization and the effects of coupling to the environment on the vibrational lifetimes of a cofactor embedded in a protein molecule. For a modest-sized organic molecule, such as a cofactor in the absence of the protein and solvent environment, energy flow among its vibrational states may be restricted, specifically when anharmonic coupling is smaller than the order of the mean spacing between resonantly coupled zero-order vibrational state energies. That picture can change when coupling to the environment is accounted for. Cofactor modes that couple directly to the environment may decay rapidly, with a lifetime that can be estimated using the golden rule. For a cofactor mode that couples only weakly to the environment, however, effects of localization found for the isolated cofactor can persist. In this case, dephasing effects give rise to mode damping, but at a rate typically below that predicted by the golden rule. We shall see that the lifetimes of these vibrational modes can be enhanced manyfold over a golden rule estimate by coherence effects, long enough to contribute to vibration-assisted conformational changes and charge transfer. As an illustrative example, we consider the lifetimes of the GFP chromophore modes, where we find, based on analogies to similar organic molecules, substantial enhancement due to coherence even in the presence of the protein environment at room temperature. We focus, in particular, on how the cofactor size affects localization of vibrational energy at physiological temperatures and on the temperature dependence of the coherence effects, which have not previously been addressed.
In the following section, we discuss quantum localization of vibrational states in molecules, summarizing a theoretical approach developed to characterize the vibrational states of large molecules [51,52]. In section 3 we apply this approach to the vibrational states of a cofactor. In section 4 we turn to the calculation of lifetimes of vibrational modes of proteins based on golden rule calculations. There we see that a quantum mechanical calculation of the vibrational lifetime can capture the temperature dependence correctly for many protein vibrations, whereas calculations based on classical modeling do not. Nevertheless, these calculations can overestimate the mode-damping rate in a cofactor if coherence effects that give rise to localization in the isolated cofactor are not taken into account. In section 5 we discuss the role of coupling to the protein and solvent environment on the mode-damping rate in a cofactor. Concluding remarks are presented in section 6.

Localization of vibrational states
Molecular vibrations, such as the vibrations of a biomolecule, can be described by a quantum mechanical system of coupled nonlinear oscillators, which in turn are coupled to a solvent environment. We begin by explicitly addressing the oscillator system. The corresponding classical system ultimately explores the entire phase space subject to specific constraints, such as constant energy, although there are conditions where transport in the phase space, via Arnold diffusion, may be very slow in many-dimensional systems [57,58]. The quantum mechanical oscillator system, on the other hand, can exhibit a localization transition to a relatively small number of states on the energy shell [52]. We discuss the origin of this transition and how to locate it for a specific molecule in this section. While quantum localization strictly exists only in the absence of coupling to the environment [52,59,60], remnants of the transition in this case still influence the rate of energy flow among vibrational states of the molecule and therefore the damping rate of its vibrational modes. We discuss the effects of coupling to the environment in the next section.
The low-order terms in V couple states close to one another in the vibrational quantum number space. We can thus picture [52] the topology of the state space as an N-dimensional lattice where each lattice site is identified by the number of vibrational quanta in each mode. For example, one zero-order state for a three-mode molecule could be identified by the lattice site (304), where one mode has three vibrational quanta, the second mode is unexcited, and the third mode has four vibrational quanta. These zero-order states, or sites on the lattice, are coupled locally to nearby sites by matrix elements arising from low-order anharmonic terms in the potential. We are interested in knowing whether the coupling between sites is sufficiently large for the molecule to explore all states on the energy shell or whether it is small enough so that the molecule remains localized to one or a relatively small number of sites.
While the energy at each site, which is given by ε α = α n ah ω α , has a definite value that depends on the frequencies of all the oscillators, it is useful in exploring the general properties of quantum energy flow to take these energies to be random from a distribution determined by the distribution of oscillator frequencies. Matrix elements arising from the anharmonic coupling in equation (1), which couple states locally in the quantum number space, also have specific values, but to explore the general properties of energy flow in the vibrational state space we can take these too as random from an appropriate distribution. For this reason we refer to this approach as local random matrix theory (LRMT) [63]. It is based on an analysis of a random matrix ensemble, where the particular molecule of interest is one member of an ensemble of similar molecules, all of which have zero-order energies and coupling terms falling within a similar range. An important constraint on members of the ensemble is that they obey selection rules, as introduced to a variety of random matrix ensembles [64]- [66], in this case due to the local anharmonic coupling between states.
The problem of vibrational energy flow in molecules, where zero-order states are sites on a lattice with random energies and coupled locally to one another, resembles the problem of single-particle quantum transport on a many-dimensional disordered lattice. Theoretical approaches to address the condensed phase problem and the emergence of Anderson localization [67] can thus be brought to bear on describing vibrational energy flow in molecules. Exploiting this connection and building on the approach of Abou-Chacra, Anderson and Thouless (AAT) to address the condensed phase problem [68], Logan and Wolynes [52] identified criteria for the localization of vibrational states, which occurs at a critical value, of the 5 order of one, of the product of the anharmonic matrix elements and local density of resonantly coupled states. The criterion for energy to flow over the energy shell in the vibrational state space is akin to a percolation threshold; if the average coupling between sites on the energy shell, or the density of resonantly coupled states, is too small, then transitions between states are confined to a local cluster on the energy shell, while sufficiently large coupling between sites, or sufficiently high density of resonantly coupled states, permits transitions among all states on the energy shell.
More formally, localization of vibrational states of a molecule depends on the self-energy of the states. The imaginary part is finite for extended states and proportional to the energy transfer rate from that state, while it is infinitesimally small for a localized state. To calculate the real and imaginary parts of the self-energy, following the approach of AAT, it is convenient to approximate the topology of the vibrational state space as a Cayley tree. This approximation becomes better justified as the number of oscillators in the system becomes large, the regime of interest when addressing the coupled vibrations of biomolecules. Starting with the Feenberg renormalized perturbation series (RPS), only the lowest order terms for the real and imaginary parts of the site self-energy are needed for a state space with a Cayley tree topology. The lowest order terms in the Feenberg RPS for the real, E j , and imaginary, j , parts of the site self-energy are given by where the absence of coupling to the environment corresponds to the limit η → 0, which we consider in this section; we address finite η in section 5. Q represents the distance in quantum number space between two states. It corresponds to the total number of quanta lost and gained in all oscillators in making a transition from the state labeled by j to states labeled by k; for example, cubic coupling can give rise to at most an exchange of three quanta among all the oscillators, in which case Q = 3, whereas quartic coupling can give rise to at most an exchange of four quanta among the oscillators, in which case Q = 4. |V Q, jk | is a matrix element coupling state j to state k, which lies a distance Q away; and K Q is the number of such states. To establish criteria for energy flow in the vibrational state space, the probability distribution of j must be solved self-consistently in the limit η → 0. In one approach to the self-consistent analysis, which we follow here, only the most probable value for the imaginary part of the self-energy, mp , is found self-consistently [51,52]. In a second approximate approach, the average inverse participation ratio is solved self-consistently [69]. In both approaches, the eigenstates of H are assumed localized, so that mp must be infinitesimally small and the average value for the inverse participation ratio must be finite. A range of molecular parameters is then identified for which mp is in fact infinitesimally small, or, similarly, where the average inverse participation ratio is finite. The results for both the selfconsistent approaches are nearly the same, differing only by a constant of order one. Details of the analysis, including incorporation of higher-order anharmonic terms and contributions of higher-order terms in the Feenberg RPS, can be found in [51]. Solving for mp we find that 6 vibrational energy flow is unrestricted when [51,52] T while quantum transport is localized in the vibrational state space at energy E when T (E) is less than 1. In equation (3), ρ Q is the local density of states that lies a distance Q away in quantum number space, and |V Q | is the average coupling matrix element to such states. If T (E) < 1, transitions only occur among a relatively small and local subset of states on the energy shell, the number of which can be estimated with the theory [69]. Above the transition, T (E) > 1, energy flows among all states of the energy shell. Schofield and Wolynes [70,71] have argued that energy flow in the vibrational state space both just above and well beyond the transition can be described by a random walk in quantum number space, a picture that has been supported by numerical calculations over a wide range of time scales [72]- [76]. The state-to-state energy transfer rate can be estimated by LRMT from the most probable value of the imaginary part of the site self-energy, which is given by [51] and is related to the most probable energy transfer rate in the vibrational state space by (4) goes over to a golden rule-like expression when T (E) 1 through a crossover region just above the transition. The golden rule expression, i.e. average value of imaginary part of the site self-energy, is The average value is finite at any energy, but when T (E) < 1 it is a poor estimate for the transition rate between states on the energy shell, since in the averaging there are many states for which is 0, and a very small number of states on the energy shell for which is very large. (In the limit N → ∞, the latter are of measure 0.) Indeed, when T (E) < 1, the most probable value of vanishes, i.e. mp = 0 in the absence of coupling to the environment.

Quantum localization in a cofactor
In the absence of coupling to the environment, vibrational states of molecules are localized and energy flow is severely restricted when T (E) < 1. The quantum ergodicity threshold of T (E) = 1 has been calculated using the theory summarized in the previous section for numerous organic molecules, including peptides [53,54] and photoactive molecules [77,78]. Moreover, the influence of quantum localization and sluggish energy flow on unimolecular reaction kinetics such as the rate of conformational isomerization has been examined [6], [79]- [82], revealing substantial deviations from transition state theory predictions resulting from these effects. Using ab initio values for the cubic anharmonic constants, or estimates based on values for other organic molecules of similar size, and direct count of the local density of vibrational states coupled to a given state on the energy shell, threshold energies ranging from about 1000 to 2000 cm −1 have been found for many organic molecules [47,53,54,77,78] with roughly 50-100 vibrational modes.  (3), as a function of temperature, where the average energy corresponding to a given temperature is used, as discussed in the text. The solid line corresponds to values of T calculated for the GFP chromophore by scaling the results for trans-stilbene, which also has 14 non-hydrogen atoms, as discussed in the text. Estimates for T of a chromophore with two more non-hydrogen atoms (dashed) and two fewer (dotted) are also shown.
Particularly relevant to exploring the possibility of quantum localization in photoactive cofactors are results for trans-stilbene, a 72-mode molecule that undergoes photoisomerization. For this molecule a localization threshold of about 1300 cm −1 was determined theoretically [77,78], consistent with experimental measurements on trans-stilbene [83] for the onset of irreversible vibrational relaxation in the S 1 state.
The value of T (E) calculated with equation (3) is sensitive to the energy, E, of the molecule. Since we will consider cofactors in contact with a thermal bath, to calculate T(E) we take the thermal energy to be E = α n α h ω α , where n α = (e¯h ω α /k B T − 1) −1 . (We use T for temperature, to distinguish from T for the transition parameter defined by equation (3).) Using the vibrational frequencies for trans-stilbene in the S 1 state, we then find E ≈ 1632 cm −1 at 300 K, an energy at which we have calculated [77] T ≈ 9, above the range of T where remnants of localization would significantly influence energy transport in the vibrational state space.
Using stilbene as a guide, we attempt here to roughly estimate the range of values of T (E) for photoactive cofactors in proteins. For instance, the GFP chromophore has the same number, 14, of non-hydrogen atoms as stilbene, but the low-frequency vibrations of the chromophore embedded in the protein are apparently higher. The 125 cm −1 mode observed experimentally appears to be the lowest frequency mode that can be characterized as a chromophore mode [40]. This frequency is 97 cm −1 higher than the lowest frequency mode of S 1 trans-stilbene in the gas phase. We can bring the low-frequency modes for the GFP chromophore roughly into line with the range used for the calculation of S 1 trans-stilbene by adding 100 cm −1 to the trans-stilbene vibrational mode frequencies. We plot values of T for this 14-atom (not counting hydrogen atoms) chromophore in figure 1. The values of T are plotted as a function of E corresponding to the temperature indicated in the figure.
We observe that, for molecules about the size of the GFP chromophore, values of T appear to be smaller than 1 even at about 300 K. We have also calculated estimates for T for somewhat 8 smaller and larger cofactors, adding or removing mode frequencies randomly in the range of vibrational frequencies for the GFP chromophore, and we plot the results in figure 1. The plot suggests that the size of these cofactors might be about the limit where vibrational lifetimes are influenced by quantum localization, at least around room temperature. The figure also suggests that for a very large organic molecule, such as the protein as a whole, the value of T at room temperature is very much larger than 1, so that transition rates between vibrational states would typically be well described by the golden rule formula. It is only when one part of the molecule is relatively weakly coupled to the rest, such as a cofactor embedded in a protein, that we turn first to addressing energy transfer within the cofactor itself, and then consider the effects of coupling to the protein/solvent environment.
The values of T that are plotted in figure 1 are seen to be exponentially dependent on ambient temperature. We note that due to the strong sensitivity of T to the ambient temperature the lifetimes of low-frequency cofactor modes may increase dramatically when lowering the temperature, since, as we shall see in section 5, mode lifetimes depend sensitively on T when T is less than or not far above 1.

Lifetimes of protein vibrational modes
For larger biomolecules under physiological conditions it is reasonable to initially neglect the effects of quantum localization on energy flow and to estimate the lifetimes of vibrations with quantum first-order perturbation theory. Unless some part of a large biomolecule, such as a cofactor, is relatively weakly coupled to it, energy transfer among vibrational states of the biomolecule should be well described by the golden rule. We thus calculate the lifetimes of most vibrational modes of proteins as for many condensed phase systems, adopting the golden rule for the damping rate of a mode with excess energy in a system where the vibrational modes are otherwise thermally populated. One relevant general feature of vibrational energy flow, as discussed in the previous sections, is its local nature, i.e. energy is typically transferred via low-order anharmonic coupling that gives rise to Fermi resonances. It is therefore reasonable to start with the lowest order anharmonic interactions, i.e. the cubic terms, to calculate the mode damping rate in a protein.
The damping rate of a vibrational mode α, W α , can be written as the sum of terms that can be characterized as 'decay' and 'collision', [84]- [86], where n α is the occupation number of mode α, which at temperature T we take to be n α = (exp(hω α /k B T) − 1) −1 , and αβγ are the coefficients of the cubic terms in the expansion of the interatomic potential in normal coordinates. The delta functions ensure energy conservation. In a finite system such as the protein considered here, calculations of the rate are carried out by replacing the delta functions with rectangles centered about the frequency ω α such that the width is larger than the spacing between the modes; we have used a width typically containing at least 3-4 modes. For most modes in the main band of protein vibrations we have found that the rate does not vary significantly upon variation of this width. However, in some cases this substitution leads to artifacts of the chosen width, particularly if energy decays from vibrational modes that are strongly localized in space and there is substantial heterogeneity in the coupling terms.
In figure 2(a) we plot the results of a calculation of the damping rate for vibrational modes of GFP, which includes the chromophore and water molecules contained inside the β-barrel, as a function of frequency to 2000 cm −1 at temperatures of 15, 135 and 300 K; the force field used for the calculation is from MOIL [87]. We observe substantial fluctuations in the decay rate and also plot the running average over 100 modes. There is significant dependence of the decay rate on temperature at lower frequency, commensurate with the factor of 20 rise in temperature, but less so at higher frequency, and in some frequency regions almost none at all. As discussed below, the classical vibrational mode damping rate is proportional to the temperature. The weak temperature dependence of the lifetime observed in many frequency regions, such as the amide I region of 1600-1700 cm −1 , as seen in the calculations plotted in figure 2(a) and in pump-probe experiments on other proteins [88], emerges from quantum mechanical perturbation theory.
This weak temperature dependence indicates anisotropic flow of energy through the vibrational states, with energy transfer from higher-frequency modes disproportionately to modes of intermediate frequency rather than low-frequency modes. Time-resolved studies of other relatively high-frequency vibrations, in addition to those of the amide I region reveal a similar trend, including pump-probe studies of myoglobin-CO, which indicate that decay of the CO stretch, about 1950 cm −1 , is essentially independent of temperature from about 6 to 310 K [89,90]. Low-frequency modes are abundant in proteins, and typically extend over many more atoms of the protein than do higher-frequency modes, as illustrated in figure 2(b), where we plot the participation ratio for each mode, α, of GFP, defined by [91] I α = ( n,r ∈x,y,z |c α,n,r | 4 ) −1 , where n is an atom and r is the x, y or z coordinate, and c is a coefficient of the normalized eigenmode. The participation ratio is a measure of the number of degrees of freedom that contribute to a vibrational mode, and we observe an overall trend of sizable I α to about 100 or 150 cm −1 , then a steep decline and plateau by about 300 cm −1 to, on average, a relatively small number of degrees of freedom contributing to vibrations.
Flow of energy to low-frequency modes of proteins is mediated by the anharmonic matrix elements that couple them to the higher-frequency modes often interrogated experimentally. Given the high density of low-frequency modes in a protein, matrix elements coupling a given high-frequency mode to a pair of other modes must be typically small if the frequency of one of the pair is small to explain the observed temperature dependence. For the matrix element coupling a triplet of modes to be appreciable, the three modes must overlap in space, i.e. there must be spatially overlapping modes in resonance for energy transfer [92,93]. Modes of higher frequency are typically spatially localized, i.e. displacements are large for a relatively small number of atoms of the protein, as illustrated in figure 2(b).
One consequence of spatial localization of normal modes in disordered systems is that frequencies of normal modes whose localization centers overlap in space are generally very different, which gives the appearance of 'repulsion' of mode frequencies between pairs of nearby localized modes [94]. We illustrate this effect in figure 3 for the backbone of GFP, which includes amide I vibrations, an effect we have illustrated earlier for myoglobin [95,96]. To reveal the repulsion between mode frequencies for localized modes close in space, we locate the atom on the backbone with the largest displacement in each normal mode, α, and calculate the probability, P( ω), that for another mode, β, for which the backbone atom with the largest displacement lies a certain distance away from that in mode α, the difference in frequency between them is ω = |ω α − ω β |. We consider only those localized modes α with frequency ω α ranging between 1000 and 2000 cm −1 . Probabilities are calculated for ω in intervals of 20 cm −1 to ω = 600 cm −1 . The dashed-line histogram is an average over all pairs of modes centered on the backbone of GFP, regardless of the distance between the localization centers of modes α and β. We observe that for any pair of modes α and β there is a greater chance that the difference between their vibrational frequencies is small than large. However, if we consider only pairs of modes localized to atoms within 2 Å of each other, about a chemical bond length, we obtain the solid-line histogram. There we see that if pairs of localized modes lie close in space, there is a propensity for their frequency difference to be large, in this case about 500 cm −1 , rather than small, say, below 200 cm −1 .
The lack of temperature dependence of the mode-damping rate observed in the amide I region and the weak temperature dependence seen in many other spectral regions are the result of first-order perturbation theory in quantum mechanics. Determination of the classical damping rate depends on the approach that is used, but the classical rate is generally smaller and depends on temperature. Stock [97] has shown that the ratio of the classical to quantum rate, when cubic coupling to modes β and γ gives rise to damping of mode α, is (1 + /2)(k B T/ω β + k B T/ω γ ) − (k B T) 2 /ω β ω γ , where is a constant between 0 and 1. The classical damping rate thus depends on temperature regardless of the frequency of the modes into which energy is transferred, in contrast to quantum relaxation.

Effects of protein and solvent environment on localization and mode damping in a cofactor
The most probable rate of energy transfer in the vibrational state space is proportional to mp (E), which we have until now expressed in the absence of interactions with the protein and solvent environment. We take this coupling into account by considering finite η. Contributions to η in the solid state problem usually include inelastic scattering due to electron-phonon interactions, which in the vibrational problem addressed here would correspond to population relaxation contributions to the dephasing rate that arise from the coupling of vibrational states of the cofactor to its environment; η can also represent rates of pure dephasing [52,98,99]. Here, to illustrate the effects of coupling of the cofactor to its environment, we use η to represent the rate of population relaxation of modes of the cofactor strongly coupled to its environment, for which we have previously calculated estimates. We then examine how such coupling affects energy relaxation from a specific mode of the cofactor that is weakly coupled to the protein and solvent environment.
Defining µ(E) = mp (E) + η, the statistical self-consistent analysis summarized in section 2 yields [51,52] which holds for any value of T(E). We take the average of the imaginary part of the site selfenergy, (E) , which is proportional to the mean energy transfer rate, as the golden rule result, which is given above by equation (5). While equation (7) is a result of the statistical selfconsistent analysis for the coupled many-oscillator model given by equation (1), where some of the oscillators are also coupled to the environment, this solution has the general form expected for the lifetime of a vibrational state that is coupled to a tier of other states, which in turn couples 6. An average mode damping rate of 2 ps −1 has been used, which, as seen in figure 2(a), is a representative value for modes above 100 cm −1 at 300 K.
to another tier of states [100], the latter mimicking the coupling of some of the vibrational states of the molecule to the environment. For this state space structure, we would expect that the lifetime, as exhibited by equation (7), generally decreases with increasing coupling of cofactor modes to the environment until such coupling becomes so large that the lifetime turns over, essentially a motional narrowing regime, at which point the lifetime increases with increasing coupling of the intermediate tier of states to the environment. These variations are illustrated in figure 4. For convenience, following Logan and Wolynes [52], we define dimensionless variables, mp (E) and η , as Equation (7) is quadratic in η with the solution [52] η (E) = Since both mp (E) and η are real and non-negative, when 0 T 2 and 0 η η max = [2 − T ]/2T 1/2 we choose the negative sign, whereas for η η max we choose the positive sign.
For T 2 only the positive sign is correct.
The dimensionless quantity mp corresponds to the ratio of the most probable rate of energy transfer in the vibrational state space to the golden rule estimate for that rate. Since relaxation of a vibrational state can be due to relaxation from any of the N modes of the molecule, the mode damping rate is ∼ mp /N . The relaxation rate of a vibrational state, mp , 13 is given by the product of mp and the golden rule rate, (E) . If we replace (E) with the golden rule result for the mode damping rate, W, calculated with equation (6), and use mp to refer to the ratio of the most probable mode damping rate to W, we can estimate the most probable mode damping rate as the product of W and mp . Similarly, we take the dimensionless dephasing rate, η , to be the ratio of the relaxation rate from modes of the cofactor to the protein and solvent environment to the golden rule result for the mode damping rate, W.
For W, we observe in figure 2(a) that the average mode damping rate for most modes of GFP is around 2-3 ps −1 at 300 K. We plot the most probable mode damping rate in figure 4, taking the average mode damping rate to be 2 ps −1 and solving equation (9) to find mp , as a function of the dephasing rate, which is the product of η and the average mode damping rate. We observe in figure 4 that the mode damping rate depends sensitively on η when T(E) is less than and even somewhat greater than 1. For a dephasing rate below 0.2 ps −1 , for example, the lifetimes of cofactor modes would be about an order of magnitude longer than predicted by equation (6) when T < 1. For a dephasing rate below 1 ps −1 and T < 0.5 the chromophore modes would be damped at a rate that is less than one-third the damping rate predicted by equation (6); and even if T is not far above 1 the lifetime can be about a factor of two longer than predicted by equation (6) for even a small dephasing rate.
We estimate η using the rate of population relaxation of those cofactor modes that strongly couple to the protein and solvent environment, using values that we have previously calculated for PYP [38,101]. Those values were calculated with equation (6) where α is a mode identifiable as a cofactor mode, and β or γ are protein or solvent modes. We use these values as estimates for η in the range 100-200 cm −1 here since, with the force field models used in the above GFP calculations we could not identify any modes in this frequency range as chromophore modes. We found for several modes of PYP between 100 and 200 cm −1 that could be characterized as chromophore modes that the lifetime due to coupling to the protein was about 600-700 fs. We assume that a specific cofactor mode of interest is only weakly coupled to the protein and that it is coupled to other modes of the cofactor, some of which relax in ≈700 fs.
Based on the data plotted in figure 1 we estimate T ≈ 0.5 for the chromophore of GFP around 300 K. For this T the mode damping rate is at least three times slower than that predicted with equation (6) for any rate of dephasing. Using as a crude estimate η corresponding to 700 fs, the lifetimes of the chromophore vibrations of GFP that are weakly coupled to the protein would be about a factor of three greater than the average lifetime predicted with equation (6), or about 1-2 ps. Our estimates for T plotted in figure 1 are of course very rough. If T were slightly smaller, say ≈0.2, the lifetime, as seen in figure 4, would be about a factor of ten greater than the average lifetime predicted with equation (6), or about 3-5 ps, comparable to the time for ESPT in GFP [40].

Concluding remarks
The localization of vibrational states of sizable organic molecules arises from quantum coherence effects. While the vibrational states are strictly extended if there is interaction with the environment, quantum transport in the vibrational state space nevertheless remains influenced by the coherence effects that give rise to localization in the same molecule in the absence of such interactions. As discussed here, mode damping rates in a molecule embedded in a protein and solvent may be significantly different, often much slower, than rates predicted using the golden rule. This occurs if the product of the coupling between resonantly coupled states of the molecule and the density of resonantly coupled states is below a critical value of order one, a criterion for quantum ergodicity that bears some resemblance to percolation. If this criterion is not met, energy flow on the energy shell is localized to a cluster of nearby states rather than over the energy shell. When coupling to the environment is included, energy flows ergodically over the energy shell in the vibrational state space of the molecule, albeit often slower than would be predicted using a golden rule formula.
We have illustrated this for a cofactor coupled to its protein and solvent environment. The localization of vibrational energy in the cofactor and the influence of coupling to the protein and solvent environment on vibrational energy flow are analogous to the role of dephasing in the energy transport in photosynthetic proteins, where the interplay between quantum coherence effects and dephasing is critical to light energy transport [102]- [106]. Because of these effects, excess vibrational energy can facilitate fast reactions in a protein rather than simply dissipate more-or-less randomly, as assumed in many standard rate theories, e.g. the Rice-Ramsperger-Kassel-Marcus theory [34,107]. At the same time, reaction would not occur if all vibrational states of the cofactor were localized, i.e. in the absence of coupling to the protein. Energy flow is required to stabilize products, and coupling to the protein may allow energy transfer from the Franck-Condon active modes to other modes overlapping the reaction coordinate at a rate that assists reaction.
As an example, we have considered the chromophore of GFP. Calculations of the mode lifetimes of GFP at several temperatures using the golden rule with simple force field models reveal quantum effects on the mode-damping rate, specifically the weak temperature dependence in many spectral regions such as the amide I, consistent with experimental observation. However, the golden-rule-based calculations also indicate that all modes between 50 and 300 cm −1 decay in less than 1 ps at 300 K. Recent ultrafast stimulated Raman spectra indicate vibrational assistance of ESPT via impulsive excitation of a ≈125 cm −1 vibration of the chromophore that is coupled to the reaction coordinate [40], which can serve to facilitate the optimal structural environment for ESPT. As observed here, the actual lifetime of such a mode may be much longer than the golden rule estimate due to remnants of quantum localization. Assuming that a particular mode of interest, e.g. one that facilitates reaction, is weakly coupled to the protein and solvent environment, but interacts with other modes of the cofactor that are significantly coupled to the environment, we expect the lifetime of the mode(s) of interest to exceed that predicted by golden rule estimates by at least a factor of three or more, as discussed in section 5. This effect can be tested by varying the temperature. As illustrated in figure 1, the transition parameter, T, varies exponentially with temperature. The rate of mode damping is very sensitive to T when T is below or not far above 1, as seen in figure 4, and is therefore far more sensitive to variations in temperature than is the golden rule estimate for the vibrational mode damping rate, which over many spectral regions depends little on temperature.
Striking changes in vibrational energy flow with small variation of temperature have been observed in helical peptides in chloroform near 270 K [108]. These changes have been attributed to a dynamical transition of the peptide-solvent system, and the existence of a dynamical transition is corroborated by molecular simulations [108]. It would be interesting to examine whether the enhancement of energy flow over such a narrow range of temperature arises from a localization transition in the peptide, an enhancement of dephasing that is coupled to a dynamical transition of the solute-solvent system, or both. Hamm and co-workers have observed a rapid rise in the dephasing rate of the C O vibration with increasing the temperature [109] to about 270 K. Indications of increased energy transfer to the C O on some of the peptide bonds in the helix as the temperature is increased above 270 K are seen in classical molecular dynamics simulations [110], and it would be interesting to compare this with the results of a quantum mechanical treatment.