The Origin of the Time-Dependent Spectral Shift in the Fluorescence

J. Phys. Chem. B 2010, 114, 11323–11337

11323

Picosecond Protein Dynamics: The Origin of the Time-Dependent Spectral Shift in the Fluorescence of the Single Trp in the Protein GB1 Dmitri Toptygin,*,† Thomas B. Woolf,‡ and Ludwig Brand† Department of Biology, The Johns Hopkins UniVersity, Baltimore, Maryland 21218, and Department of Physiology, The Johns Hopkins UniVersity School of Medicine, Baltimore, Maryland 21205 ReceiVed: May 14, 2010; ReVised Manuscript ReceiVed: July 22, 2010

How a biological system responds to a charge shift is a challenging question directly relevant to biological function. Time-resolved fluorescence of a tryptophan residue reflects protein and solvent response to the difference in π-electron density between the excited and the ground state. In this study we use molecular dynamics to calculate the time-dependent spectral shift (TDSS) in the fluorescence of Trp-43 in GB1 protein. A new computational method for separating solvent, protein, and fluorophore contributions to TDSS is applied to 100 nonequilibrium trajectories for GB1 in TIP3P water. The results support several nontrivial conclusions. Both longitudinal and transverse relaxation modes of bulk solvent contribute to the TDSS in proteins. All relaxation components slower than the transverse relaxation of bulk solvent have significant contributions from both protein and solvent, with a negative correlation between them. Five exponential terms in the TDSS of GB1 are well separated by their relaxation times. A 0.036 ps term is due to both solvent (60%) and protein (40%). Two exponential terms represent longitudinal (τL ≈ 0.4 ps) and transverse (τD ≈ 5.6 ps) relaxation modes of TIP3P water. A 131 ps term is attributable to a small change in the tertiary structure, with the R-helix moving 0.2 Å away from the β-strand containing Trp-43. A 2580 ps term is due to the change in the conformation of the Glu-42 side chain that brings its carboxyl group close to the positively charged end of the excited fluorophore. Interestingly, water cancels 60% of the TDSS resulting from this conformational change. 1. Introduction A time-dependent spectral shift (TDSS)1 has been experimentally observed in the fluorescence of solvatochromic dyes in polar solvents2-10 and in the fluorescence of solvatochromic fluorophores incorporated in biological macromolecules, such as proteins11-35 and DNA.36,37 While there is a broad consensus that the TDSS in the emission from dyes in polar solvents reflects the relaxation dynamics of the solvent,3 there is still no agreement on the origins of TDSS in biological macromolecules. Some authors attributed experimentally observed TDSS in proteins to the relaxation of the protein matrix.11-14,16,19,24,27,33,34 Others argue that the relaxation of the solvent (water) is entirely responsible for TDSS in proteins.18,21,22,25,28,32,35 There are also several reports in which the authors found experimental evidence for contribution from both water and protein matrix to the TDSS.15,20,26,29-31 One approach that may help to separate the contributions of water and protein matrix to TDSS in proteins is based on the different time scales of these contributions. In bulk water, TDSS occurs on the subpicosecond time scale.4 TDSS in proteins covers a wide range of time scales, from femtoseconds to at least tens of nanoseconds.11-35 Abbyad et al.29 measured TDSS at seven different sites within the same protein. Both subpicosecond and slower (picosecond and nanosecond) relaxations were found to be contributing to TDSS at every site, however, the subpicosecond relaxations were found to be dominant at those sites where the fluorophore was in contact with water, * To whom correspondence should be addressed. Phone: (410) 516-7300. Fax: (410) 516-5213. E-mail: [email protected]. † Department of Biology. ‡ Department of Physiology.

and the slower relaxations were found to be dominant at the sites buried deep inside the protein. Based on these findings, Abbyad et al.29 attributed the subpicosecond TDSS to the relaxation of water, and the slower TDSS at t > 1 ps to the relaxation of the protein matrix. In accordance with this point of view, in those experiments where the time resolution was slower than 1 ps, the experimentally observed TDSS would be entirely due to protein relaxation, which is in agreement with the conclusions of most studies where the time resolution was insufficient to detect the relaxation of bulk water.11-13,19,24,27,34 A different point of view on the origin of slow TDSS in proteins was expressed by a group of authors.18,21,22,25,28,31,32,35 According to the latter point of view, the motion of water molecules in the vicinity of the protein is highly constrained, and this “biological water” layer is the sole origin of the slow TDSS in proteins.18,21,22,25 In later papers, it is suggested that the fluctuations of water and protein atoms are coupled and TDSS is attributed to the motion of water molecules in the protein hydration layer.28,31,32,35 Thus, there are still two diametrically opposite points of view on the origin of slow TDSS in proteins, and experimental evidence alone cannot help to reach a consensus regarding this subject. Molecular dynamics (MD) simulations offer a way to separate the contributions of water and protein atoms to TDSS. Numerous MD simulations of TDSS in proteins have been reported.21,30,31,38-44 In some of these reports the contributions of water and protein atoms to TDSS were separated,31,39,40,42,43 and here the points of view regarding the contributions of water and protein to slow TDSS diverged again. According to Nilsson and Halle,40 water cannot contribute to slow TDSS. According to Golosov and Karplus,43 depending on the location of the solvatochromic

10.1021/jp104425t  2010 American Chemical Society Published on Web 08/11/2010

11324

J. Phys. Chem. B, Vol. 114, No. 34, 2010

fluorophore, TDSS may be due to water, protein, or both water and protein. In addition, several groups of authors found independently that in those cases where both protein and solvent contribute to the spectral shift, the contributions from protein and water have a negative correlation, that is, when protein shifts the emission spectrum to the red, water shifts it to the blue and vice versa.39-42 Hassanali et al.42 explained this finding in terms of a competition between protein and water for the solvation of the fluorophore. Later, Halle and Nilsson45 published an eloquent explanation for the negative correlation based on a dielectric continuum model. It is likely that the differences in the methods of converting MD trajectories to spectral shifts contributed in part to the differences in the conclusions regarding the relative contributions of protein and water. Nonequilibrium hybrid quantum mechanics-molecular dynamics (QM-MD) simulations were used to calculate TDSS in all papers that came out of the group of Callis.38,39,46 A much faster approach, based on Coulomb’s equation, was employed by Hassanali et al.42 to estimate TDSS from nonequilibrium MD trajectories. A large group of authors used the linear-response method to estimate TDSS.30,40,41,43 This method makes it possible to estimate TDSS from correlations of random fluctuations in an equilibrium MD trajectory. Finally, Li et al.31,44 applied both the nonequilibrium MD and the linear-response method in their theoretical investigation of the TDSS in the emission of Trp-7 in apomyoglobin. In the present report we review the assumptions on which the non-QM method of calculating TDSS from nonequilibrium MD and the linear-response method are based and define the limits of applicability of these methods. Furthermore, we reexamine the methods of separating the contributions of water and protein to the TDSS. In previous studies the effects of translation and rotation of the fluorophore in the local electric field were not separated from that due to the changes in the electric field itself. This is equivalent to studying the motion of protein and solvent atoms in a reference frame that moves together with the fluorophore. Consider a hypothetical situation where all water molecules and most of the protein atoms (except for the fluorophore) are fixed in the laboratory frame, while the fluorophore is allowed to rotate. Based on the old method of separating the contributions, the TDSS observed in this hypothetical situation would be attributed in part to the protein and in part to water, even though water molecules did not move in this example. We correct this problem in our new method described below. Furthermore, we also split the protein contribution into the contributions of individual atoms, which makes it possible to visualize the physical mechanism of TDSS on the atomic level. We also emphasize the importance of comparing the results of MD simulations to an experimentally measured TDSS. Furthermore, for proteins and other macromolecules one must compare the absolute amplitudes of the spectral shifts rather than the normalized response functions S(t),40,43 also denoted as CF(t),3 CS(t),41 and C(t),30 because the absolute values at t ) 0 and t ) ∞ are not certain (t ) ∞ is never achieved in MD simulations of macromolecules in explicit solvent, and in experimental work neither t ) 0 nor t ) ∞ can be achieved due to the finite time resolution of the instrument and the limited excited-state lifetime of the fluorophore). For an experimental test of our new computational methods, we chose the fluorescence of the only native Trp residue W43 in the B1 domain of streptococcal protein G (GB1). In our previously published experimental study,27 the system was carefully examined for the presence of heterogeneity, and

Toptygin et al. heterogeneous population decay was ruled out based on the width conservation criterion.19 Moreover, out of all single-Trp protein variants examined in our laboratory, the native Trp in GB1 exhibits the largest-amplitude slow TDSS on the picosecond and nanosecond time scale.27 What is also important, the experimentally measured slow TDSS of Trp in GB1 is more than 90% over in just 2 ns after excitation, which makes it possible to limit nonequilibrium MD trajectories to just 2 ns in the excited state. For most other proteins studied in our laboratory, the TDSS does not even start to level off by t ) 20 ns,19,47 thus, at least a 10-fold longer MD simulation would be required for those proteins. GB1 protein was also experimentally studied by Cohen and Boxer and their co-workers;20,29 in both studies, the TDSS was observed in the fluorescence of a synthetic amino acid Aladan, which represents a solvatochromic fluorophore attached to alanine. The TDSS of Aladan incorporated in GB1 instead of the native Trp residue W4320 is qualitatively similar to the TDSS of the Trp-43 studied in our laboratory.27 All these considerations, as well as the known X-ray crystal structure48 and NMR solution structures49,50 of GB1 influenced our selection of Trp in GB1 as a model system for the present study. 2. Theory This section consists of four parts. In the first part we consider the relationship between TDSS and the local electric field. In the second part we introduce a new method of separating the contributions to TDSS resulting from the motions of different atoms, atomic groups, and molecules. In the third part we review the connection between TDSS in polar solvents and the theory of dielectric relaxation in homogeneous liquids. In the fourth part we consider the errors resulting from the use of the linearresponse method. 2.1. Relationship between TDSS and the Local Electric Field. Spectral shifts in both absorption and fluorescence emission of solvatochromic fluorophores have their origin in the Stark effect.51,52 The energies of the singlet electronic states S0 and S1 vary with the local electric field acting upon the fluorophore. The electric field is generated by charged and polar molecules and chemical groups in the environment of the fluorophore. Any motion of these molecules and groups results in a time variation of the local electric field. The Stark effect translates the variation of the local electric field into a TDSS. It is convenient to characterize the spectral shifts in fluorescence emission using either the peak frequency νpk or the center-ofgravity frequency νc of the instantaneous emission spectrum. In pump-dump-probe experiments like those described in ref 16, νpk and νc are the peak and the center-of-gravity of the instantaneous absorption spectrum. If the experimental timeresolved emission (or absorption) spectra are corrected as described in Appendix A (see Supporting Information), then the following equation holds true

hν(t) ) hν0 + ∆E(t)

(1)

Here h is the Planck constant; ν(t) represents the time variation of either νpk or νc; hν0 is a constant; for a strict definition of hν0, see Appendices A and B; and ∆E(t) is the time variation of the difference between the Stark shifts in the energies of the electronic states S0 and S1. Quantum mechanical (QM) calculations provide the most accurate way of calculating ∆E(t) for a fluorophore in an electric field. Hybrid QM-MD simulations of TDSS in proteins and

Origin of Time-Dependent Spectral Shift in GB1

J. Phys. Chem. B, Vol. 114, No. 34, 2010 11325

peptides have been reported by the group of Callis.38,39,53 Hybrid QM-MD simulations usually cover only a short period of time, from 238 to 30 ps39 for proteins and up to 60 ps for small peptides.53 In those cases where several nanoseconds of TDSS must be calculated, the use of QM-MD simulations is not desirable because of their slow speed. This has encouraged many authors to look for faster methods of calculating ∆E(t) from MD trajectories. In Appendix B (see Supporting Information) we derive the following approximate equation for calculating ∆E(t):

∆E(t) )

∑ ∆Qj(0)φ(rj, t)

(2)

j

Equation 2 has been derived from quantum-mechanical perturbation theory, which has long been used to describe the Stark effect in atoms and molecules.54 The summation is carried out over all fluorophore atoms, rj is the radius vector of atom j, and ∆Qj(0) is defined as the difference (0) (0) ∆Qj(0) ) Qj1 - Qj0

Qj1(0)

Qj0(0)

rlab ) Rrmol + u

(5)

rmol ) RT(rlab - u)

(6)

(3) 55

and represent Mulliken atomic charges on atom j of the fluorophore in the excited state S1 and in the ground state S0, respectively, when the fluorophore is acted upon by a constant and uniform electric field, see Appendix B in Supporting Information for details. The electric potential φ(r,t) is generated by the partial electric charges of all atoms, including solvent and protein atoms, except the atoms of the fluorophore,

φ(r, t) )

greater atomic displacements than those relative atomic motions that actually produce TDSS. It makes even less sense to use the reference frame attached to a small atomic group (e.g., one side chain) of the protein molecule or to a water molecule. The only reasonable choice is to attach the reference frame to the protein molecule. Because the molecule is flexible, the question how this can be practically accomplished is not trivial. This question is considered in Methods, and here we assume that the reference frame is attached to the protein in such manner that the rotation and translation of the protein as a whole does not change the coordinates of its atoms in this reference frame. Coordinate transformations between the laboratory reference frame and that attached to the protein molecule are described by the following equations:

1 4πε0

Q

∑ |r - ri i(t)|

(4)

i

Here ε0 is the dielectric permittivity of vacuum, index i counts all atoms except those that belong to the fluorophore, Qi is the partial electric charge on atom i, and ri is the radius vector of this atom. Equations 1-4 provide a fast way of calculating TDSS directly from a nonequilibrium MD trajectory. Nonequilibrium MD is triggered by switching the partial charges of the fluorophore atoms from Qj0(0) to Qj1(0) at t ) 0. TDSS calculated using eqs 1-4 from just one MD trajectory is overwhelmed by random noise, therefore, averaging over many trajectories is usually necessary. It must be emphasized that eq 2 is approximate. The errors resulting from the use of eq 2 and the limits of its applicability are considered in Appendix C (see Supporting Information). 2.2. Separation of Contributions from Different Motions. According to eqs 1, 2, the spectral shift depends on the radius vectors rj of the fluorophore atoms and also on the electric potential φ(r,t) generated by the partial charges of all atoms except the atoms of the fluorophore. Thus, the motion of the fluorophore as well as the motion of the surrounding atoms contribute to the spectral shift. Here we will attempt to separate the contributions resulting from the motions of the fluorophore, from the motions of other protein groups and from the motion of solvent molecules. Ultimately we will define the contribution to the TDSS from the motion of each individual atom. To be able to tell which atoms moved and how far, we will have to choose a reference frame. It makes little sense to consider atomic motions in the laboratory frame, because in this frame the rotation and translation of the protein as a whole results in much

Here rlab denotes the column vector of three Cartesian coordinates xlab, ylab, zlab in the laboratory reference frame; these are the coordinates used in MD simulations. Likewise, rmol denotes the column vector of three Cartesian coordinates xmol, ymol, zmol in the reference frame attached to the protein molecule. Throughout this paper equations containing radius vectors r without superscripts are valid regardless of the choice of the reference frame; for example, eqs 2 and 4 are valid both in the laboratory and in the molecular reference frame. If radius vectors from more than one reference frame appear in the same equation, then each radius vector is provided with a superscript to identify its reference frame, as in eqs 5 and 6. R denotes a 3 × 3 rotation matrix of the special orthogonal group SO(3); its transpose RT is equal to its inverse. u is the translation vector, which describes the position of the origin of the molecular reference frame in the laboratory reference frame. The rotation and translation of the protein molecule results in the time variation of the matrix R and vector u. For every protein atom n the vector of mean equilibrium coordinates 〈rnmol〉 is defined as the ensemble averages of rnmol at long times after excitation. In practice one can average rnmol over all MD trajectories and also over a time period that starts after the bulk of the TDSS is over and continues to the end of the trajectory. Note, that averaging is carried out on the coordinates of protein atoms only and in the reference frame attached to the protein; it makes no sense to average any Cartesian coordinates in the laboratory reference frame or to average the coordinates of solvent atoms in the protein reference frame. Using the mean coordinate vectors 〈rjmol〉 for the fluorophore atoms it is possible to split the energy ∆E(t) in eq 2 into the part ∆Ef(t) that varies due to the motion of the fluorophore itself and the part ∆Ee(t) that varies due to the changes in the electric field:

∆E(t) ) ∆Ef(t) + ∆Ee(t)

(7)

11326

J. Phys. Chem. B, Vol. 114, No. 34, 2010

∆Ef(t) )

Toptygin et al.

∑ ∆Qj(0)[φ(rjmol, t) - φ(〈rjmol〉, t)]

(8)

∑ ∆Qj(0)φ(〈rjmol〉, t)

(9)

j

∆Ee(t) )

j

where ∆Qj(0) is the charge difference defined in eq 3 and φ is the electric potential defined in eq 4. If the fluorophore does not move in the reference frame attached to the protein, then all the differences in square brackets on the right-hand side of eq 8 equal zero, and therefore, ∆Ef(t) also equals zero. On the other hand, in the hypothetical situation where all protein atoms except those of the fluorophore itself and all water atoms are fixed in the reference frame attached to the protein, the potentials φ(〈rjmol〉,t) are time-invariant, and therefore, ∆Ee(t) is also timeinvariant. In other words, without the motion of the fluorophore there is no time variation in ∆Ef(t), and without the motions of other atoms there is no time variation in ∆Ee(t). This justifies treating ∆Ef(t) as the contribution of the fluorophore atoms and ∆Ee(t) as the contribution of nonfluorophore protein atoms and solvent atoms. To split ∆Ee(t) further into the contributions of protein and solvent atoms, the electric potentials φ(〈rjmol〉,t) that appear in eq 9 must be expressed explicitly in terms of the coordinates of these atoms,

φ(〈rjmol〉, t) )

1 4πε0

Q

∑ |rmol(t) -i 〈rmol〉| i

i

(10)

1 4πε0

∆Qj(0)Qi

∑ |rmol(t) - 〈rmol〉| j

i

to the electric displacement D, then the relaxation of the electric field E does not contain the term exp(-t/τD). Instead, it contains a faster exponential term exp(-t/τL),59 as shown in the middle panel of Figure 1. The relationship between τL and τD is well-known,3,56-59

j

Here, index i counts all atoms except those that belong to the fluorophore, Qi is the partial charge on atom i, and rimol the coordinate vector for this atom. The time variation of φ(〈rjmol〉,t) results only from the time variation of rimol, because the partial atomic charges Qi and the mean coordinates 〈rjmol〉 are timeinvariant. The contribution of one nonfluorophore atom i to ∆Ee(t) and also to ∆E(t) can be obtained by taking just one term from the sum in eq 10 and substituting it in eq 9,

∆Ei(t) )

Figure 1. Graphical definitions of the transverse relaxation time τD (top), the longitudinal relaxation time τL (middle), and an example of a simple structure in which both exp(-t/τD) and exp(-t/τL) are present in the dielectric response (bottom). E(t) is the electric field, D(t) is the electric displacement, and θ(t) is Heaviside step function.

(11)

j

The contribution of one fluorophore atom j to ∆Ef(t) and also to ∆E(t) can be obtained by taking just one term from the sum in eq 8. The solvent contribution can be determined by summing up ∆Ei(t) from eq 11 over all water atoms. The contribution of the protein can be determined by summing up ∆Ei(t) from eq 11 over all nonfluorophore protein atoms and adding ∆Ef(t) from eq 8 to the sum. The protein contribution, water contribution, and contributions of individual atoms must be then averaged over all nonequilibrium trajectories (ensemble averaging). 2.3. Dielectric Relaxation of Bulk Solvent. Dielectric relaxation of the solvent usually has a significant contribution to the TDSS of fluorophores in proteins. Dielectric relaxation of polar solvents occurs on multiple time scales.56 The Debye model of dielectric relaxation in polar liquids contains only one characteristic time constant τD, which is commonly called the Debye relaxation time.3,56-59 An exponential term of the form exp(-t/τD) is directly observed in the time variation of the electric displacement D after a jump in the electric field E in the form of Heaviside step function θ(t),59 as shown in the top panel of Figure 1. However, if the step-function jump is applied

τL )

ε∞ τ εS D

(12)

In eq 12, εS is the static (low-frequency) dielectric permittivity of the solvent, ε∞ is the high-frequency dielectric permittivity of the solvent, τL is the longitudinal relaxation time, and τD is the transverse relaxation time and also the Debye relaxation time. The terms “transverse” and “longitudinal” are relevant in the case of a continuous homogeneous polar liquid without borders or foreign objects. In this case, the dielectric polarization density P can be separated into two independent parts: (i) the longitudinal part, for which ∇ · P * 0 and ∇ × P ) 0, and (ii) the transverse part, for which ∇ · P ) 0 and ∇ × P * 0.56 If the dielectric response is described by the Debye model, then both the longitudinal and the transverse component relax exponentially, but with different relaxation times τL and τD.56,59 The dielectric polarization near a spherical ion immersed in a continuous homogeneous polar solvent consists of the longitudinal part only, therefore a step-function jump in the charge of the ion would induce an exponential relaxation process with the characteristic time τL. Such a charge jump contradicts the principle of electric charge conservation; therefore, it cannot be achieved in practice. A step-function jump in the electric dipole moment is practically achievable, and it occurs when a solvatochromic fluorophore jumps from the ground state to the excited state. Maroncelli and Fleming57 have shown that for a point-dipole centered inside a spherical cavity a step-function jump in the dipole moment induces a relaxation process containing a single exponential term exp(-t/τF), where57,58

τF )

2ε∞ + εC τ 2εS + εC D

(13)

Origin of Time-Dependent Spectral Shift in GB1 In eq 13 εC is the dielectric permittivity of the spherical cavity, which is intended to represent fluorophore polarizability. From eqs 12 and 13, it follows that εC ) 0 results in τF ) τL, while εC ) ∞ results in τF ) τD. According to realistic estimates,3,57,58 for spherical fluorophores in solvents of high polarity τF is only slightly greater than τL, while τD is much greater than both τL and τF. The above example with a dipole jump in a spherical cavity shows that the observed dielectric relaxation cannot be always separated in two parts with the characteristic relaxation times τL and τD. Complete separation of the transverse and longitudinal polarization components can be accomplished only in the absence of borders and foreign objects (such as the spherical cavity in the above example). In the case of a parallel-plate capacitor, the concepts “transverse” and “longitudinal” become completely irrelevant. The electric field E between the plates of such a capacitor is uniform and so is the dielectric polarization density P, therefore, ∇ · P ) 0 and ∇ × P ) 0, which fits neither the definition of longitudinal polarization nor the definition of transverse polarization.56 Yet, with a parallel-plate capacitor one can observe both relaxation times τL and τD. If a step-function jump is applied to the voltage across the capacitor, then the charge on each plate of the capacitor will change with time as D(t) in the top panel of Figure 1. However, if a stepfunction jump is applied to the charge on the plates of the same capacitor, then the voltage across the capacitor will change with time as E(t) in the middle panel of Figure 1. In all previous examples, the relaxation curve contained only one exponential term, which could be exp(-t/τL), exp(-t/τF), or exp(-t/τD). Now the question is whether two or more exponential terms can be observed simultaneously in a relaxation curve if the solvent is described by the Debye relaxation model. A physical structure in which both exp(-t/τL) and exp(-t/τD) can be observed simultaneously is depicted in the bottom panel of Figure 1. A layer of a nonpolar solid dielectric (such as polyethylene or diamond), shown by the yellow color in Figure 1, is immersed in a polar liquid (such as water), shown by the blue color. A narrow cylindrical channel whose diameter is much smaller than the thickness of the layer, is drilled through the solid dielectric and is also filled with the polar liquid. The whole structure is placed between the plates of a parallel-plate capacitor and a step-function jump is applied to the charge on the plates of this capacitor. This creates a step-function jump in D1 (the electric displacement in the polar liquid). Because D1 and D2 are both normal to the interface between the two dielectrics, D1 ) D2 and, therefore, there is also a step-function jump in D2 (the electric displacement in the nonpolar dielectric). For the nonpolar dielectric the static dielectric constant and the high-frequency dielectric constant is the same quantity, which will be denoted εNP. The relationship between the electric field and the electric displacement in the nonpolar dielectric is instantaneous, that is, D2 ) εNPE2, which results in a stepfunction jump in E2 (the electric field in the nonpolar dielectric). Because E2 and E3 are both parallel to the walls of the cylindrical channel, E2 ) E3, and therefore, there is also a stepfunction jump in E3 (the electric field in the polar liquid inside the channel). A step-function jump in the electric field in the polar liquid results in a relaxation process that contains the exponential term exp(-t/τD). The relaxation of the polar liquid outside the channel will result in the term exp(-t/τL). Therefore, both τL and τD can be observed at the same time in an electrostatic experiment. It should also be possible to observe two exponential terms, one with a τ close to τL and one with a τ close to τD, in an experiment where a solvatochromic

J. Phys. Chem. B, Vol. 114, No. 34, 2010 11327 fluorophore is embedded in a nonpolar dielectric of irregular shape and immersed in a polar liquid. The hydrophobic core of a protein can play the role of the nonpolar dielectric, however, an irregular shape made of polyethylene or diamond can play this role equally well. An exponential term with a τ close to τD is expected only if the nonpolar structure contains at least one narrow channel filled with the polar solvent and if the length of the channel significantly exceeds its diameter or if the channel connects a large-area internal pocket with the external solvent or with another large-area internal pocket. If the structure contains no internal pockets and no channels, then no slow relaxation terms with characteristic times close to τD will be observed, but there can still be more than one exponential term with τ between τL and τD. When the shape of the dielectric structure approaches a perfect sphere, the number of exponential terms in the TDSS reduces to one with τ ) τF given by eq 13. 2.4. Linear-Response Method. Nonequilibrium MD is commonly triggered by switching the partial charges of the fluorophore atoms from Qj0(0) to Qj1(0) at t ) 0. TDSS is then calculated using eqs 1-4 directly from a nonequilibrium MD trajectory; this approach will be called the direct-response method to distinguish it from the linear-response method. TDSS calculated using the direct-response method from just one MD trajectory is overwhelmed by random noise, therefore averaging over many trajectories is usually necessary. An alternative approach, known as the linear-response method, makes it possible to obtain the same information from the autocorrelation of random fluctuations in just one very long equilibrium trajectory. The method is described below. Consider some interaction with energy ∆E(t) that can be instantaneously turned on. For example, Coulombic interaction between the additional charges ∆Qj(0) on fluorophore atoms and the common partial charges of other atoms instantaneously turns on when the fluorophore jumps from the ground state to the excited state. Equation 2 gives the energy of this interaction. According to the fluctuation-dissipation theorem, in the linearresponse approximation, the ensemble average response 〈∆E(t)〉 averaged over an infinite number of nonequilibrium trajectories, in which the interaction is turned on at t ) 0, is connected to the autocorrelation function C(t) defined below by a simple relation:

〈∆E(t)〉 - 〈∆E(∞)〉 ) βC(t)

(14)

Here β is the inverse temperature, β ) 1/(kBTA), where kB is Boltzmann constant and TA is absolute temperature. The autocorrelation function C(t) can be defined either as ensemble average, that is, the average over infinite number of trajectories, or as a time average over one infinitely long trajectory. We will use the second approach, because in practice the autocorrelation function is always obtained from one trajectory.30,31,40,41,43,44 Time-averaged C(t) is defined as the limit at T f ∞ of the autocorrelation function obtained from a finite-length trajectory,

C(t) ) lim CT(t) Tf∞

(15)

Here CT(t) is calculated from a trajectory of length T; it is defined as follows:

CT(t) )

1 T-t

∫0T-t [∆E(t′) - ∆ET][∆E(t′ + t) - ∆ET]dt′ (16)

11328

J. Phys. Chem. B, Vol. 114, No. 34, 2010

Toptygin et al.

Figure 2. Random noise patterns (panels 1A, 2A, and 3A) and corresponding autocorrelation functions CT(t) (panels 1B, 2B, and 3B). The black line in each #B panels represents C10 ns(t) calculated from the 10 ns trajectory in the corresponding #A panel using eqs 16 and 17. The red line represents C10000 ns(t) calculated from the 10000 ns trajectory generated by the same random process as the 10 ns trajectory in the corresponding #A panel. The green and the blue lines in panel 2B depicts C5 ns(t) calculated from for the first 5 ns and from the last 5 ns of the trajectory in panel 2A, respectively.

∆ET )

1 T

∫0T ∆E(t′)dt′

(17)

∆E(t′) must be calculated from an equilibrium trajectory, that is, the interaction should not be turned on or turned off during the MD simulation that generates ∆E(t′) to be used in eqs 16 and 17. The results obtained using the linear-response method contain both random errors and systematic errors. The random errors originate from substituting CT(t) instead of C(t) into eq 14. For a short trajectory the difference between CT(t) and C(t) can be substantial, which can easily lead to false interpretations. The most common artifacts in CT(t) resulting from averaging over a finite-length trajectory are illustrated in Figure 2. The random noise in panel 1A of Figure 2 is a sum of the random fluctuations in 16 noise-driven damped harmonic oscillators with different frequencies and different damping ratios (Langevin dynamics). The frequencies, damping ratios, and amplitudes for each of the 16 oscillators can be found in Supporting Information. This random noise visually resembles ∆E(t) from a real MD simulation, yet, computing this random noise is millions of times faster than a real MD simulation for a real protein, which makes it easy to study the evolution of CT(t) up to T ) 10000 ns. Panel 1A shows only the first 10 ns of ∆E(t), however, the “trajectory” was actually generated all the way up to t ) 10000 ns. In panel 1B of Figure 2, the black line depicts C10 ns(t), calculated from the part of the trajectory in panel 1A. The red line depicts C10000 ns(t), which is practically identical to C(t), because we have established that the increase in T from 5000 to 10000 ns produces no visible changes in CT(t). The difference between C10 ns(t) and C10000 ns(t) is the lowfrequency noise, which is a common artifact associated with the use of the linear-response method. Theoretically, it should

be possible to decrease the amplitude of the low-frequency noise to a desired level; however, in practice, this would require a very long MD trajectory because the amplitude of noise in CT(t) decreases as 1/T. The low-frequency noise, also known as pink noise, is potentially more dangerous than the white noise in direct-response TDSS because the low-frequency noise is not subjectively perceived as noise. The noise pattern in panel 1A of Figure 2 resembles nearequilibrium atomic fluctuations in a protein that does not change its conformation at all. If the entire protein or just one side chain in close proximity to the fluorophore jumps between two discrete isomers, then the noise pattern may look like the ones shown in panels 2A or 3A. ∆E(t) in panels 2A and 3A were generated using a random switch that jumps between two discrete levels; for details, see Supporting Information. The switch in panel 2A spends equal time (on average) in both states; the average time between jumps for this switch equals 5 ns. The switch in panel 3A spends 90% of time in the lower energy state and 10% in the upper energy state, with the residence times of 0.33 ns for the upper energy state and 3 ns for the lower energy state. The random noise from the 16 damped harmonic oscillators was also added to the signals depicted in panels 2A and 3A. ∆E(t) curves were generated for a total of 10000 ns in all cases. Panel 2A of Figure 2 shows a selected 10 ns window during which the protein spent the first 5 ns as isomer 1 (low ∆E) and the following 5 ns as isomer 2 (high ∆E). The autocorrelation C10 ns(t) calculated from the selected window (black line in panel 2B) and the autocorrelation C10000 ns(t) calculated from the full trajectory (red line in panel 2B) are close to each other, and this is for a good reason: the average times the protein spends as isomer 1 or isomer 2 equal 5 ns, and during the selected 10 ns window it spent exactly 5 ns as isomer 1 and 5 ns as isomer 2. Random switching between isomers makes a large contribution to C(t) and it is solely responsible for the difference between the red line in panel 2B and the red line in panel 1B. This switching increases the full amplitude of the autocorrelation function from 1.0 in panel 1B to 7.2 in panel 2B. It also increases the correlation time from less than 0.5 ns in panel 1B to more than 2.5 ns in panel 2B. Thus, switching between two (or more) isomers could be the main source of the slow TDSS observed experimentally on the time scales of nanoseconds or even tens of nanoseconds.11-13,19,24,27,34 The green and blue lines in panel 2B depict C5 ns(t) calculated from for the first 5 ns (during which only isomer 1 was present) and from the last 5 ns (during which only isomer 2 was present) of the trajectory in panel 2A. Because there were no transitions between the isomers during the time windows used to calculate either C5 ns(t), the slow TDSS is not observed in this case. An example of splitting MD trajectory in two parts can be found in the work of Li et al.,31,44 who separately analyzed the trajectories for the isomers 1 and 2 using the linear-response method. No nanosecond-scale relaxation was found in this work.31,44 This example clearly shows how insufficient trajectory length T can radically alter the results obtained using the linear-response method. The example shown in panels 3A and 3B illustrates a possible but unlikely situation. During the 10 ns window selected for panel 3A the protein jumped from the low ∆E state (isomer 1) to the high ∆E state (isomer 2) and back three times in a row. ∆E(t) does not behave like this during every 10 ns window, thus, the piece of trajectory shown in panel 3A is atypical. This explains the big difference between the black and the red line in panel 3B. The black line represents C10 ns(t) calculated from the selected 10 ns window. The red line represents C10000 ns(t) calculated from the full trajectory, and it closely approximates

Origin of Time-Dependent Spectral Shift in GB1

J. Phys. Chem. B, Vol. 114, No. 34, 2010 11329

C(t). While the red line approaches the zero level aperiodically, the black line shows oscillations characteristic of an underdamped oscillator. Golosov and Karplus43 found the underdamped behavior in just one of the eleven linear-response MD trajectories they simulated. This can be the case of a real underdamped oscillator or the case where the low-frequency noise resembles the behavior of an underdamped oscillator, similar to that shown in panels 3A and 3B. To find out which is the case, one would have to generate a much longer linearresponse MD trajectory. The results obtained using the linear-response method may also contain systematic errors. The systematic errors have their origin in deviations from linearity. Equation 14 is valid only in the linear-response approximation; therefore, in the case where β|∆E| > 1, that is, |∆E| > kBTA, the use of eq 14 can result in significant systematic errors. The systematic errors will be further investigated here using a very simple model system, which consists of just one TIP3P water molecule in a uniform electric field that is turned on at t ) 0. The orientation of the water molecule is described here using Euler angles (φ,θ,ψ).60 The Z′ axis of the molecular frame is chosen parallel to the permanent electric dipole moment µ of the water molecule, and the Z axis of the laboratory frame is chosen parallel to the electric field61 E. With this choice of axes, the interaction energy depends on only one Euler angle:

∆E(t) ) -µ · E ) -|µ||E| cos θ

(18)

Using the fact that in the absence of the electric field the orientation of the water molecule is random and therefore 〈cos2θ〉 ) 1/3, we can calculate the value of C(t) for t ) 0:

C(0) ) 〈∆E2(t)〉 ) |µ| 2 |E| 2〈cos2 θ〉 )

1 2 2 |µ| |E| 3

(19) According to the linear-response method,

1 〈∆E(0)〉 - 〈∆E(∞)〉 ) βC(0) ) β|µ| 2 |E| 2 3

(20)

Taking into account that 〈∆E(0)〉 ) 0 and 〈∆E(∞)〉 ) |µ||E|〈cosθ〉, one can obtain the following expression for the equilibrium value of 〈cos θ〉 in the presence of the electric field:

1 〈cos θ〉 ) β|µ||E| 3

(21)

The estimate in eq 21 is based on the linear-response approximation, therefore, it is accurate in weak electric fields only. A general expression for the equilibrium value of the firstrank order parameter 〈cos θ〉 in the electric field of any strength can be obtained using the equilibrium orientational distribution of the molecule, which is essentially a Boltzmann distribution in Euler angles:

d3P(φ, θ, ψ) ) A exp(β|µ||E| cos θ)dφ sin θdθdψ

(22) The constant A in eq 22 must be chosen so that the integral over all orientations equals unity. Multiplying the distribution in eq 22 by cosθ and integrating over all orientations yields

〈cos θ〉 ) coth(β|µ||E|) - (β|µ||E|)-1

(23)

In Figure 3 the estimate obtained using the linear-response method, eq 21, is depicted by the red line, and the value from eq 23, which is accurate at any level of the electric field, is depicted by the blue line. The curves in Figure 3 were calculated using β ) 1/(kB × 300 K), which corresponds to the temperature of 300 K, and |µ| ) 7.829 × 10-30 C · m, which is the magnitude of the permanent electric dipole moment of a TIP3P water molecule. The relative deviation of the linear-response estimate from the actual value reaches 10% at the electric field strength of 7 MV/cm. The additional excited-state charges ∆Qj(0) on the Trp side chain (see Methods) generate an electric field of 7 MV/ cm magnitude at distances up to 8.5 Å from the center of the fluorophore. Thus, if the distance between the center of the indole moiety and the center of at least one water molecule is less than 8.5 Å, the amplitude of the TDSS obtained using the linear-response method is likely to be overestimated. A water molecule in van der Waals contact with the Trp side chain may experience electric fields up to 50 MV/cm; in electric fields of this magnitude, the estimate obtained using the linear-response method is 3.5-fold greater than one obtained using the directresponse method. An overestimated amplitude is always accompanied by a distorted shape of the TDSS; therefore, renormalization does not solve the problem, it just makes it more difficult to be acknowledged. Maroncelli and Fleming62 questioned the validity of the linear-response method when they discovered that in direct-response MD simulations the solvation response to a charge jump in ST2 water does not vary linearly with the magnitude of the charge jump. Comparing the theoretical estimates in eqs 21 and 23 makes it possible to determine in advance whether the results of a future linearresponse MD simulation will be valid or not. 3. Methods 3.1. Trp Atomic Charges. In MD simulations and in calculating TDSS from MD trajectories we used Trp partial atomic charges listed in Table 1. These charges were obtained by Callis and co-workers as described in refs 39 and 46 and represent the average over the last 4 ps of a hybrid QM-MD trajectory of 3-methylindole in 1La excited state in a drop of 1100 explicit TIP3 waters. The charges were originally calculated in accordance with the definition of Mulliken55 and then scaled by a factor of 0.80, as described in ref 39. The charges were kindly made available to us by Callis on May 14, 2007. 3.2. MD Simulations. All-atom MD simulations in explicit TIP3P solvent were performed in a 40.36 × 40.36 × 40.36 Å rectangular box with periodic boundary conditions using the CHARMM program. CHARMM22 all-hydrogen force field was modified slightly to allow unequal charges on the carbon and hydrogen atoms of the Trp side chain. Modified CHARMM22 topology files top_all22_prot_PCallis_G.inp (with ground-state Trp charges Qj0(0) from Table 1) and top_all22_prot_PCallis_E.inp (with 1La excited-state Trp charges Qj1(0) from Table 1) are available, see Supporting Information. Protein Data Bank file 1PGB.pdb63 served as the initial structure of GB1 protein. The protein was solvated using 1834 TIP3P water molecules, energyminimized, gradually heated to the final temperature of 278 K (experimental fluorescence measurements were made at this temperature), and equilibrated with the ground-state charges for an average time of 5000 ps (progressively increasing equilibration time was used to obtain the initial conditions for 100 different trajectories, with an additional 100 ps equilibration between two consecutive trajectories). During the production

11330

J. Phys. Chem. B, Vol. 114, No. 34, 2010

Toptygin et al.

Figure 3. First-rank order parameter 〈cos θ〉 for a TIP3P water molecule in a uniform electric field E. Red line: the estimate obtained using the linear-response method, eq 21. Blue line: the exact value from eq 23. The curves were calculated for 300 K temperature and 7.829 × 10-30 C · m dipole moment.

TABLE 1: Trp Side Chain Atom Charges in Units of Positive Electron Charge e Trp atom

Qj0(0)

Qj1(0)

∆Qj(0)

CB HB1 HB2 CG CD1 HD1 NE1 HE1 CE2 CD2 CE3 HE3 CZ3 HZ3 CZ2 HZ2 CH2 HH2

-0.003 +0.013 +0.016 -0.046 +0.011 +0.030 -0.130 +0.125 +0.052 -0.009 -0.027 +0.027 -0.033 +0.005 -0.039 +0.017 -0.026 +0.017

+0.004 +0.017 +0.015 +0.184 +0.122 +0.047 +0.017 +0.132 +0.082 -0.087 -0.214 +0.028 -0.077 -0.003 -0.243 +0.021 -0.063 +0.018

+0.007 +0.004 -0.001 +0.230 +0.111 +0.017 +0.147 +0.007 +0.030 -0.078 -0.187 +0.001 -0.044 -0.008 -0.204 +0.004 -0.037 +0.001

phase each trajectory was run for 100 ps in the ground state and then for 2000 ps in the excited state (i.e., with the excitedstate Trp charges). MD time step was set to 0.002 ps via the use of “shake bonh”. Nonbonded interactions were limited using cutnb 14.0, ctofnb 12.0, and ctonnb 10.0. Because we used explicit solvent, the dielectric constant ε was set to 1.0. The coordinates of all atoms were saved after each 0.1 ps during the ground-state run, each 0.004 ps during the first 5 ps of the excited-state run, and then each 0.1 ps during the remaining 1995 ps of the excited-state run. One of the 100 trajectories was continued for 30 ns instead of 2 ns in the excited state; the last 26 ns of this long trajectory were used to calculate the autocorrelation CT(t) to be used in the linear-response method. The first 2 ns in the excited state from this trajectory and the other 99 trajectories were used to calculate the ensemble mean in accordance with the direct-response method. 3.3. Calculation of the Total TDSS from MD Trajectories. The total TDSS was calculated using the atomic coordinates in the laboratory reference frame, that is, the coordinates used in the MD simulations. For every fluorophore atom j the electric potential φ(rjlab,t) was calculated using eq 4. The summation in eq 4 was carried out over all protein atoms (except the

fluorophore atoms) and all solvent atoms in the main rectangular box and also in 26 adjacent boxes that were created from the main box by the symmetry transformations (periodic boundary conditions). The electric potentials φ(rjlab,t) were substituted in eq 2 to calculate ∆E(t). For the charges ∆Qj(0) in eq 2, we used the atomic charges from the last column of Table 1, multiplied by e to convert them to SI units. To obtain the direct-response TDSS the energy ∆E(t) was ensemble-averaged over 100 nonequilibrium MD trajectories and then divided by hc to convert the result from energy units to wavenumber units. To obtain the linear-response TDSS we used only one 30 ns long excited-state trajectory. The first 4 ns of this trajectory was discarded. The remaining 26 ns of this trajectory was used to calculate C26 ns(t) as defined in eqs 16 and 17. The latter was multiplied by β and then divided by hc to convert the result from squared energy units to wavenumber units. 3.4. Reference Frame Attached to the Protein. To separate the contributions of fluorophore, solvent, and protein to the total ∆E(t) in accordance with eqs 5-7, it is necessary to know the mean equilibrium coordinates 〈rnmol〉 of each protein atom n in the reference frame attached to the protein molecule. If it was possible to define the reference frame attached to the protein prior to determining the mean coordinates in that frame, then the calculation of the mean coordinates would be straightforward. Unfortunately, the reference frame attached to the protein can be defined only in terms of the mean coordinates 〈rnmol〉 and the mean square deviations from these mean coordinates. In the present study, the reference frame attached to the protein was determined and the coordinate transformations in eqs 5 and 6 were carried out by a Fortran subroutine Match3D, which is described in Supporting Information. For each point in each trajectory, subroutine Match3D finds the rotation matrix R and the translation vector u that minimize the sum

σ)

mol 2 ∑ wn|rlab n - (R〈rn 〉 + u)|

(24)

n

In eq 24 the weight wn equals the inverse mean square deviation for protein atom n, rnlab is the column vector of the current coordinates of that atom in the laboratory reference frame, and 〈rnmol〉 is the column vector of the mean coordinates of that atom in the reference frame attached to the protein molecule. The rotation matrix R and the translation vector u that minimize σ in eq 24 uniquely define the reference frame attached to the protein. Furthermore, for a protein that has a rigid core and some disordered regions, this algorithm naturally attaches the frame to the rigid core and disregards the motion of the atoms in disordered regions via the use of the inversemean-square-deviation weighting. The mean coordinates 〈rnmol〉 and the weights wn can be found in a series of successive iterations, as described below. The iterative procedure operates on a set of protein structures obtained in the course of a series of MD simulations. In this work we used all structures saved between 1000 ps and 2000 ps after excitation along each of the 100 nonequilibrium MD trajectories. Since the coordinates were saved each 0.1 ps, the set consisted of exactly one million structures. The zeroth approximation represents a rough estimate of the mean coordinates 〈rnmol〉 and the weights wn prior to the first iteration. For the zeroth approximation mean coordinates we take the first structure in the set, however, any other structure chosen at random could play this role equally well. For the

Origin of Time-Dependent Spectral Shift in GB1 zeroth approximation weights we use wn ) 1 for all backbone atoms and wn ) 0 and for all side chain atoms. Iteration number i (i ) 1, 2, 3, ...) consists of two steps. During the first step the reference frame attached to the protein is calculated using the (i-1)-th approximation mean coordinates and the (i-1)-th approximation weights. The coordinates of every atom in the attached reference frame are averaged over all structures and the mean coordinates become the i-th approximation mean coordinates. During the second step, the reference frame attached to the protein is calculated using the i-th approximation mean coordinates and the (i-1)-th approximation weights. The mean square deviation (averaged over all structures) of the coordinates of every atom in the attached reference frame from the i-th approximation mean coordinates is calculated and the inverse of this mean square deviation becomes the i-th approximation weight wn. Iterations are repeated until the weights stop changing. When the iterative procedure described here was applied to the set of one million structures of GB1 protein in the relaxed excitedstate configuration, three iterations turned out to be sufficient (the weights after the third iteration were practically the same as those after the second iteration). 3.5. Avoiding Singularities in Component TDSS. The expression for the electrostatic potential φ(r,t) in eq 4 is not applicable if the distance from the point where the potential is calculated to the center of at least one atom is less than the radius of that atom. Indeed, the electric field inside the atom cannot be calculated by placing the net charge of all the electrons and the nucleus at the center of the atom. This does not present a problem as far as the calculation of the total ∆E(t) is concerned, because during a tenable MD simulation the distance between the centers of two nonbonded atoms never becomes smaller than the radius of one of these atoms. For example, in the MD simulations described here the distance between the centers of any two nonbonded atoms was never equal to or less than 1.2 Å at any time in any of the 100 trajectories simulated. This means that if instead of |r - ri| there was max(1.2 Å, |r - ri|) in the denominator of eq 4, then the total time-dependent spectral shift calculated using eqs 1-4 would not change at all. On the other hand, replacing |rimol - 〈rjmol〉| with max (1.2 Å, |rimol - 〈rjmol〉|) in the denominator of eqs 10 and 11 can make a big difference. Once in a while the current coordinates of an atom i may be very close to the mean coordinates of an atom j; when this happens, the contribution of atom i to the spectral shift calculated without limiting the denominators in eqs 10 and 11 experiences a sharp spike. These spikes may be positive and negative with equal probability, and theoretically they should cancel out when the results are averaged over an infinite number of trajectories. In practice, however, we do not have the luxury of an infinite number of trajectories and, therefore, limiting the minimum distance between two points to 1.2 Å in the denominators of eqs 4, 10, and 11 presents a simple practical way to reduce the amplitude of noise in the contributions of fluorophore, solvent, protein, and individual protein atoms to the calculated TDSS. As was stated earlier, this operation does not affect the mean value and/ or the noise in the total TDSS, it only affects the noise in component TDSS. 3.6. Calculation of the Component TDSS from MD Trajectories. The component TDSS were calculated using the atomic coordinates in the molecular reference frame. This frame was individually defined for each time point and each trajectory as described in section 3.4. Atomic coordinates from the trajectory files rilab and rjlab were converted to the coordinates

J. Phys. Chem. B, Vol. 114, No. 34, 2010 11331 in the molecular frame rimol and rjmol using eq 6 and then substituted in eqs 10 and 11. The denominators in eqs 10 and 11 were limited to avoid singularities, as described in section 3.5. The summation in eq 10 was carried out over all protein atoms (except the fluorophore atoms) and all solvent atoms in the main rectangular box and also in 26 adjacent boxes that were created from the main box by the symmetry transformations (periodic boundary conditions). The electric potentials φ(〈rjmol〉,t) obtained from eq 10 were substituted in eq 9 to calculate ∆Ee(t). ∆Ef(t) was obtained in accordance with eq 7, as the difference between ∆E(t) and ∆Ee(t). ∆Es(t) was calculated using the same operations that were used to calculate ∆Ee(t), except that the summation in eq 10 was carried out over solvent atoms only and not over the protein atoms. ∆Ei(t) was obtained from eq 11. ∆Ep(t) was calculated by summing ∆Ei(t) over all protein atoms i except the atoms of the fluorophore and then adding ∆Ef(t) to the sum. 3.7. Fitting Nonequilibrium Relaxation Dynamics by Sums of Exponentials. The total TDSS as well as the contributions to it from the fluorophore, solvent, protein, and individual protein atoms were fit by the following model function, Nexp

∆Ec(t) ) Rc,0 +

∑ Rc,n exp[-max(t, 0)/τn]

(25)

n)1

In eq 25, c represents the name of the component (f for fluorophore contribution, s for solvent, p for protein, i for the contribution of the specific protein atom i, and blank for the total effect), n is the number of the exponential term, Rc,n is the corresponding amplitude, and τn is the corresponding relaxation time. The values of all fitting parameters Rc,n and τn were determined simultaneously in the course of a global weighted nonlinear least-squares minimization using the Fortran program Grelaxfit. The program uses the model function defined in eq 25 and is universally applicable for fitting nonequilibrium relaxation dynamics of any origin. The data for the fitting were obtained as the ensemble averages of the component TDSS over all available 100 trajectories. The standard deviations of the ensemble averages were calculated as described in ref 64. Inverse squares of these standard deviations were used as the weights during the global weighted nonlinear least-squares minimization by the program Grelaxfit. 4. Results and Discussion 4.1. Linear-Response Results. The autocorrelation function C26 ns(t) calculated from the last 26 ns of the 30 ns long excitedstate trajectory is depicted in Figure 4. The system is believed to have approached the excited-state equilibrium during the first 4 ns after excitation. The first 4 ns were not involved in calculating C26 ns(t). In accordance with eqs 1 and 14, C(t) · β/h represents the TDSS on the frequency scale, and C(t) · β/hc represents the TDSS on the wavenumber scale, which is commonly used in experimental work.27 The amplitude of the TDSS in Figure 4 equals 4626 cm-1, which exceeds the upper limit based on experimental data, see below. The t f ∞ limit of the mean65 emission wavenumber for GB1 in water was experimentally determined in ref 27 and equals 27957 cm-1; therefore, in accordance with the linear-response method, the t f 0 value of the mean65 emission wavenumber should be 32583 cm-1, which causes a problem, because this falls within the absorption band of the Trp residue in GB1 (the extinction coefficient at 32583 cm-1 equals 315 M-1 cm-1), and no

11332

J. Phys. Chem. B, Vol. 114, No. 34, 2010

Figure 4. Autocorrelation function C26 ns(t) calculated from the last 26 ns of the 30 ns long trajectory. Inset shows the first 5 ps of C26 ns(t) on the expanded time scale.

fluorophore absorbs light at its mean65 emission frequency.66 This argument suggests that the amplitude obtained using the linear-response method is exaggerated, as it is predicted for those proteins that have water molecules within 8.5 Å from the Trp side chain, see section 2.4. Indeed, the crystal structure of GB1 shows that the Trp residue in this protein is partially accessible to the solvent.63 The linear-response TDSS in Figure 4 also reveals significant levels of low-frequency noise. The noise can be visually recognized on the 13 ns time scale (main panel) but not on the 5 ps time scale (inset). This does not mean that the lowfrequency noise is absent at t < 5 ps. The noise is present at all times, however, since this noise contains no high-frequency components, within the narrow time window the low-frequency noise results in a constant bias. The level of the constant bias is not reproducible from one linear-response trajectory to another, but if the TDSS is calculated from just one trajectory, then it is not obvious that the contribution of the low-frequency noise is present. 4.2. Direct-Response Results. The TDSS calculated from 100 nonequilibrium trajectories using the direct-response method is shown in Figure 5. On the time scale of Figure 5 it is impossible to show subpicosecond relaxation components, therefore, in Figure 6, early parts of the relaxation curves are shown on a 400-fold expanded time scale. In both figures, dots represent ensemble mean TDSS averaged over 100 MD trajectories, and solid lines represent the best global fits to the dots by the multiexponential model function from eq 25, with Nexp ) 5. When fewer than five exponential terms is used, the quality of the fit becomes significantly worse. Increasing Nexp from five to six produces no improvement in the quality of the fit. This does not prove that the relaxation dynamics of GB1 is quintuple-exponential. It shows only that at the present signalto-noise level it is impossible to resolve more than five exponential terms in the TDSS. The number of trajectories would have to be increased significantly to reach a signal-tonoise level sufficient for resolving six or more exponentials terms.

Toptygin et al.

Figure 5. TDSS obtained using the direct-response method, ensembleaveraged over 100 trajectories. Different colors are used to depict the total TDSS (red), the solvent contribution to TDSS (blue), the contribution from all protein atoms, including fluorophore (green), and the contribution from the fluorophore only (orange). Dots represent the ensemble mean values of ∆E(t)/hc and solid lines represent the best global fits by the multiexponential model function from eq 25 with Nexp ) 5.

Figure 6. Direct-response TDSS data from Figure 5, shown on a 400fold expanded time scale to reveal the early stages of the relaxation process. The color scheme and the roles of the dots and solid lines are the same as in Figure 5.

The values of the five relaxation times τn resulting in the best global fit to the data are τ1 ) 36.1 ( 1.6 fs, τ2 ) 384 ( 26 fs, τ3 ) 5.63 ( 0.33 ps, τ4 ) 131 ( 5 ps, and τ5 ) 2.58 ( 1.06 ns. The amplitudes associated with the five exponential terms and their decomposition into the contributions from the solvent, protein, and fluorophore are shown in Figure 7 in the form of a histogram. In section 4.3 we discuss each of the five

Origin of Time-Dependent Spectral Shift in GB1

J. Phys. Chem. B, Vol. 114, No. 34, 2010 11333

Figure 8. Representative relaxed excited-state structure of GB1, depicted using stick representation in PyMOL.68 GB1 atoms are colored in accordance with their contributions to the amplitude of the ultrafast relaxation component, τ1 ) 36.1 fs. Green color means no contribution, red color represents a large positive contribution, and blue color represents a large negative contribution. Intermediate colors correspond to small-magnitude contributions (orange and yellow for positive contributions, cyan or light blue for negative). Figure 7. Histogram of the amplitudes corresponding to the five exponential terms recovered in the global fitting of the direct-response TDSS data by the model function in eq 25. Red color denotes the total amplitude Rn. Blue color denotes the solvent contribution Rs,n. Green color denotes the protein contribution Rp,n. Orange color denotes the fluorophore contribution Rf,n. Error bars represent 95% confidence intervals. Corresponding values of τn are shown at the bottom.

exponential terms individually and explain the mechanics of the corresponding motions. The sum of the five amplitudes equals 3130 cm-1 and represents the change in the emission wavenumber between t ) 0 and t ) ∞. The experimental t f ∞ limit of the mean65 emission wavenumber for GB1 equals 27957 cm-1,27 therefore, the t ) 0 value equals 31087 cm-1. The mean65 wavenumber for 1Lb emission of Trp in a nonpolar environment equals 31416 cm-1. The number was calculated from the experimental emission spectrum of 3-methylindole in hexane, see Supporting Information. Because the energy gap 1Lb - 1A is not sensitive to the electric field,46 we can assume that it is approximately the same for the Trp residue in GB1 protein. The electronic state 1La in GB1 (31087 cm-1) is lower than the 1Lb state (31416 cm-1) by about 329 cm-1 at the instance of excitation and also in the ground state. The uncommon ground-state ordering of the singlet excited states 1La and 1Lb in the in GB1 protein is due to the permanent electric field generated by the charges on the protein atoms. The data shown in Figures 5 and 6 indicate that at t e 0 the net ∆E/hc equals -1310 cm-1, of which protein contributes -1700 cm-1 and the solvent contributes +390 cm-1. In other words, at t e 0 a 1700 cm-1 red shift due to the charges on the protein atoms is opposed by a 390 cm-1 blue shift from the charges on TIP3P water molecules. A similar situation with a large red shift from the protein charges and a smaller blue shift from water was found in a hybrid QM-MD simulation of Trp residue W140 in Staphylococcal Nuclease, W58 in CheY protein, and W126 in T4-lysozyme.39 4.3. Exponential Relaxation Components and Corresponding Motions. 4.3.1. Ultrafast Component. The total amplitude of the ultrafast (τ1 ) 36.1 fs) relaxation component equals 1600 cm-1 or 4.57 kcal/mol. Approximately 59% of this relaxation amplitude is due to the motion of the solvent molecules and 41% is due to the motion of protein atoms. The contribution from the fluorophore motion to this relaxation

component is less than 3% and it is not statistically significant. Ultrafast relaxation of bulk water on a time scale faster than 50 fs was found both in experimental studies4 and in hybrid QMMD simulations of 3-methylindole in TIP3P water.46 In both cases, this was the largest-amplitude relaxation component. The 36.1 fs component also has the largest amplitude among the five components found in this MD simulation of GB1 protein, however, only 59% of this amplitude is due to librations of water molecules and 41% is due to small adjustments in the protein structure, which are illustrated in Figure 8. The H and N atoms of Trp-43 backbone NH group are the most significant contributors on this time scale. The NH group is very close to the positive end of the excited-state Trp-43 fluorophore. In the excited state the positively charged H atom (red color in Figure 8) moves a little further away and pulls the negatively charged N atom (blue color in Figure 8) with it. This results in a 168 cm-1 red shift from the H atom and 115 cm-1 blue shift from the N atom, which travels a shorter distance than the H atom. The net spectral shift from Trp-43 backbone NH group is 53 cm-1. Gly-41 residue backbone CO group has a greater net contribution of 90 cm-1, which consists of a 92 cm-1 red shift due to the O atom (orange color in Figure 8) and only a 2 cm-1 blue shift from the C atom (green color in Figure 8). A 72 cm-1 net contribution to the ultrafast relaxation component is from Lys-31 side chain NH3 group, where each H atom (yellow in Figure 8) contributes a 34 cm-1 red shift, and the N atom (lightblue in Figure 8) contributes a 30 cm-1 blue shift. Asp-40 side chain carboxyl group contributes a net 66 cm-1 red shift, and there is a large number of charged or polar groups, each contributing 10 cm-1 or less. It must be emphasized that protein structural adjustments contributing to the 36.1 fs relaxation component are limited to very small changes of bond lengths, bond angles, and dihedral angles. 4.3.2. Second-Fast Component. The total amplitude of the second-fast (τ2 ) 384 fs) relaxation component equals 536 cm-1 or 1.53 kcal/mol. All of this relaxation amplitude is due to the motion of the solvent molecules. The contribution of protein to this relaxation component is slightly negative and it is not statistically significant, see Figure 7. The longitudinal relaxation time τL of bulk water can be calculated by substituting into eq 12 the values ε∞ ) 5.2, εS ) 78.36, and τD ) 8.27 ps, obtained from the frequency dependence of water dielectric constant.67

11334

J. Phys. Chem. B, Vol. 114, No. 34, 2010

This yields τL ) 550 fs. Because the value of τ2 is so close to the longitudinal relaxation time τL of bulk water, it is fair to say that the second-fast relaxation component actually represents the longitudinal relaxation mode of bulk solvent. The fact that τ2 is 30% less than τL calculated from eq 12 can be attributed to the difference between TIP3P water model used in this MD simulation and the real water, which always results in shorter MD relaxation times as compared to the experiment. For example, in a hybrid QM-MD simulation of 3-methylindole in TIP3P water at 300K Muin˜o and Callis46 found a major exponential relaxation component with τ ) 400 fs, which agrees with the value of τ2 ) 384 ( 26 fs obtained here, but is also about 30% less than τL of real water. Experimental TDSS of coumarin 343 in bulk water was reported to contain two close exponentials with τ ) 126 fs and τ ) 880 fs,4 with the amplitude-weighted mean of 606 fs, which is not far from τL ) 550 fs. 4.3.3. Intermediate Component. The total amplitude of the intermediate (τ3 ) 5.63 ps) relaxation component equals 393 cm-1 or 1.12 kcal/mol. All of this relaxation amplitude is due to the motion of the solvent molecules. The contribution of protein to this relaxation component is near zero and it is not statistically significant, see Figure 7. Considering that the value of τ3 is about 30% less than the transverse relaxation time τD ) 8.27 ps67 and that the 30% difference is attributable to TIP3P water model (see section 4.3.2), it is fair to say that this intermediate relaxation component actually represents the transverse relaxation mode of the bulk solvent. A relaxation time close to τD was not seen in MD simulations of 3-methylindole in bulk TIP3P water46 or in experiments with coumarin 343 in bulk water.4 In section 2.3, we explained why an exponential term with τ ≈ τD is absent in the TDSS of a fluorescent molecule in a homogeneous polar solvent but can be present in the TDSS of the same molecule when it is inserted in a structure containing water channels in a solid dielectric material. 4.3.4. Second-Slowest Component. The total amplitude of the second-slowest (τ4 ) 131 ps) relaxation component equals 137 cm-1 or 0.39 kcal/mol. Of this amplitude, the motion of water molecules contributes 206 cm-1 or +150% and the motion of the protein atoms (including fluorophore) contributes -69 cm-1 or -50%. The motion of the fluorophore contributes 83 cm-1 or +60%. Significant contributions from protein atoms indicate that some conformational change takes place on the time scale of τ4. To understand what type of conformational change is responsible for the TDSS on this time scale, the ensemble mean internal coordinates (bond lengths, bond angles, and dihedral angles) were fitted with the model function from eq 25. In this fitting all τn were fixed at the values listed in section 4.2. Careful inspection of the amplitudes corresponding to τ4 reveals that not a single internal coordinate undergoes a large change on the time scale of this relaxation component. There is practically no change in the ensemble mean dihedral angles χ1 and χ2 of Trp-43, indicating that the side chain of Trp-43 does not change its rotameric state. Furthermore, the side chain of no other residue changes its rotameric state on the time scale of τ4 either. Modest changes (not exceeding 3°) in the ensemble mean values of some backbone dihedral angles φ and ψ are observed on the time scale of τ4, and most of these changes are observed in the turn regions, residues 9-12, 20-23, and 36-42, see Figure 9. This indicates that the secondary-structure elements (one R-helix and four β-strands) do not change their conformation, but an adjustment of the tertiary structure takes place on the time scale of τ4. To better

Toptygin et al.

Figure 9. Representative relaxed excited-state structure of GB1, depicted using cartoon representation in PyMOL.68 The color represents the change in the ensemble mean distance between the R-carbon of each residue and the R-carbon of Trp-43 that takes place on the time scale of the relaxation component τ4. A residue that does not move closer or further away from Trp-43 is depicted by cyan color. The residues that move away from Trp-43 are depicted by green, yellow, orange, and red colors, with the red color corresponding to the largest positive change in the distance (about +0.23 Å) on the time scale of τ4. The residues that move closer to Trp-43 are depicted by the shades of blue color, with the dark blue color corresponding to the largest negative change in the distance (about -0.06 Å) on the time scale of τ4.

understand this structural adjustment, the ensemble mean distance between the R-carbon of each residue and the R-carbon of Trp-43 was calculated and fitted using the model function from eq 25 with fixed τn values. The amplitude associated with τ4 are shown by color in Figure 9. The largest positive change in the distance on the time scale of τ4 is only about +0.23 Å (red color in Figure 9). The largest negative change in the distance on the time scale of τ4 is only about -0.06 Å (blue color in Figure 9). The cartoon in Figure 9 makes it clear that on the time scale of the 131 ps relaxation component the entire R-helix (residues 23-36) and especially its N-terminal end (residues 23-25, red color in Figure 9) move further away from the Trp residue, which allows a slightly greater water access to the Trp side chain, and results in the red shift. The net contribution of all protein atoms to the TDSS is negative on the time scale of τ4, therefore, the energy for this adjustment in the tertiary structure is not derived from the Coulombic interactions within the protein but from the Coulombic interactions between protein atoms and water. Specifically, water molecules are attracted to the increased excited-state dipole moment of the Trp residue and they push the R-helix away from the Trp residue as much as the flexibility of the tertiary structure allows. The rate of this tertiary structure adjustment is defined by the elastic and kinetic properties of the protein, but the energy comes from the Coulombic interaction between the Trp side chain and the solvent. 4.3.5. Slowest Component. The total amplitude of the slowest (τ5 ) 2.58 ns) relaxation component equals 466 cm-1 or 1.33 kcal/mol. Of this amplitude the motion of water molecules contributes -703 cm-1 or -150% and the motion of the protein atoms (including fluorophore) contributes +1170 cm-1 or +250%. The motion of the fluorophore contributes 68 cm-1, which is statistically significant, but represents only a small part of the protein contribution. Fitting of the ensemble mean internal coordinates by the model in eq 25, as described in section 4.3.4, revealed that just two dihedral angles undergo large systematic changes on the time scale of the 2.58 ns relaxation component,

Origin of Time-Dependent Spectral Shift in GB1

J. Phys. Chem. B, Vol. 114, No. 34, 2010 11335

Figure 11. Dynamics of Glu-42 side chain rotamer populations. Curve labels C, D, E, F, and H correspond to the tile labels in Figure 10. Rotamers A, B, G, and I have negligible populations. In Figure 8 the side chain of Glu-42 is depicted in the conformation of rotamer E, which is the most populated rotamer at any time during this simulation, but its population decreases with time. Figure 10. Set of transient Glu-42 side chain conformations, depicted by 2.1 × 106 color dots on a plot of χ2 vs χ1. The internal coordinates (χ1,χ2) were saved at 0.1 ps intervals along 100 trajectories. The color of each dot is blue at early times (ground state and the early part of the excited-state trajectory). During the first 1 ns in the excited state, the dot color gradually shifts from blue to green, and during the second 1 ns, the color shifts from green to red. Overlaps between dots are colored using pigment mixing rules rather than intensity addition rules. Variable brightness enhancement prevents areas where too many dots of different colors overlap from going completely black. Broken lines divide the χ1-χ2 space into nine rotamers, labeled with letters A through I. Each rotamer represents a topologically connected area of high dot density on this plot. To prevent splitting of the connected areas at the artificial boundaries where χ1 ) (180° or χ2 ) (180°, the angles χ1 and χ2 have been redefined so that their domain is 0° e χ < 360° rather than -180° < χ e 180°.

and these are the dihedral angles χ1 and χ2 of Glu-42. This shows that the side chain of Glu-42 changes its rotameric state on the time scale of τ5. A conformation of Glu-42 side chain can be depicted as a point on a plot of χ2 versus χ1. Figure 10 depicts as color dots 2.1 × 106 points (χ1,χ2) saved at 0.1 ps intervals along each of the 100 trajectories. The color of the dot shifts from blue at t e 0, through green at t ) 1 ns, and to red at t ) 2 ns. The dots in Figure 10 tend to group into five islands, with each island representing a different Glu-42 rotamer. Islands D, E, and F contain a significant number of blue dots, which shows that the corresponding rotamers are populated in the ground state. Islands C and H contain predominantly red dots, which shows that the corresponding rotamers become populated in the excited state. The time variation of the populations of the five rotamers is depicted in Figure 11. The shift in the Glu-42 rotamer population from rotamers D, E, and F to rotamers C and H is responsible for the 2.58 ns component of the timedependent red shift, because rotamers C and H have their negatively charged carboxyl group much closer to the positively charged (in the excited state) end of the Trp fluorophore. The negative contribution of the solvent molecules to the 2.58 ns component can be roughly explained in terms of the dielectric model originally suggested by Halle and Nilsson.45 Assuming that the carboxyl group and the fluorophore are material points and all space around them is occupied by water (i.e., no space is occupied by the protein) and that water is a continuous dielectric with a linear response and the dielectric constant εS,

the total TDSS from the carboxyl group in water must be equal to 1/εS of that in the absence of water. The contribution of water equals the TDSS in the presence of water minus the TDSS in the absence of water, thus, the contribution of water equals (1/εS - 1) times the TDSS due to the carboxyl group in the absence of water. The latter represents the contribution of the protein. In our calculations, the contribution of water equals -0.60 times the contribution of the protein, therefore (1/εS - 1) ) -0.60 and εS ) 2.5. This is much less than the dielectric constant of water (εS ) 78.36 at 25 °C),67 suggesting that our initial assumptions are too simplistic. Specifically, the assumption that the protein occupies no space and the further assumption that water is a linear dielectric at the electric field levels that exist near the Glu-42 carboxyl group and the excited-state Trp side chain are too coarse. Nevertheless, on the qualitative level the dielectric model of Halle and Nilsson45 does explain very well why the contribution of water to the 2.58 ns component is negative and why it cancels a major part of the protein contribution to this relaxation component. 4.4. Comparison with the Experimental Data. Experimental TDSS of the Trp residue in GB1 in water at +5 °C was reported in ref.27 The time resolution of the instrument in that study was about 50 ps, therefore, the exponential terms τ1, τ2, and τ3 could not be resolved. The number of exponential terms with τ > 50 ps in this MD simulation equals two, as compared to three exponential terms found in the experiment. This difference is easy to explain. The number of Glu-42 rotamers having nonzero populations equals five, see Figure 10. In the general case, the kinetics of a system with five discrete states is described by four exponential terms, plus there is one more exponential term resulting from a minor adjustment of GB1 tertiary structure, see section 4.3.4. Some of the five exponential terms may have equal or close τ values, which would make it impossible to resolve them from the mixture. Three exponential terms (with the τ values of 77 ps, 637 ps, and 5.56 ns and the amplitudes of 520, 148, 38 cm-1, respectively) were resolved in the experiment,27 which shows that some of the five exponentials are indeed too close. In the present MD simulation the time scale ends at 2.0 ns, therefore all exponential terms with τ > 2.0 ns could not be resolved from one another, and this explains why only two exponential terms with τ > 50 ps

11336

J. Phys. Chem. B, Vol. 114, No. 34, 2010

were resolved. A longer MD simulation is not feasible because it already took more than a year to generate one hundred 2.0 ns long nonequilibrium trajectories. The sum of the amplitudes of the fourth and the fifth term in this MD simulation equals 603 cm-1. The sum of the amplitudes of three exponential terms found in the experimental TDSS equals 706 cm-1. The difference between the two values is only 15%, which can be attributed to the random noise in the simulated TDSS (insufficient number of trajectories) and to experimental errors (a significant uncertainty is expected in the amplitude of the 77 ps exponential because it is too close to the time resolution). While the amplitudes of the simulated and experimental TDSS seem to be in a reasonable agreement, the rates of the spectral shifts are not. For example, the most significant exponential term in the experimental TDSS has the amplitude of 520 cm-1 and the relaxation time τ ) 77 ps. This can be compared only to the fifth exponential term in this MD simulation, because the amplitude of the fourth term is too small, and the first three terms are too fast to be resolved in the experiment. The fifth exponential term has the amplitude of 466 cm-1 and the relaxation time τ5 ) 2.58 ns. The experimental relaxation rate is about 33-fold faster than the one obtained in the MD simulation. This shows that there is a 33-fold difference between the experimental and simulated rate of interconversion between the rotamers of Glu-42 side chain. The rotation about the χ2 axis of Glu-42 appears to be too slow in this MD simulation, although Glu-42 side chain is on the surface and its conformational changes are not obstructed by other protein atoms. Time-resolved anisotropy measurements of fluorescence emission from Trp residues on protein surfaces reveal that the rotation about χ1 and χ2 axes for such Trp residues takes place on the time scale of 55 ps69 or 90/270 ps.70 Considering that a Glu side chain is not much heavier than the Trp side chain, one would expect rotation on a similar time scale. It is not clear then why it takes 2.58 ns in this MD simulation to change the rotameric state of the Glu-42 side chain. 5. Conclusions The work described here provides clear answers to four unrelated questions. The first question is whether the relaxation of the bulk solvent (water) can contain relaxation terms with τ g 5 ps. Debye relaxation model predicts two characteristic relaxation times for every polar solvent: the longitudinal relaxation time τL and the transverse relaxation time τD. For real water at 25 °C, τL ) 550 ps and τD ) 8.27 ps; for TIP3P water both τL and τD are ≈30% shorter. This study shows that both τL and τD can be observed in the TDSS of a Trp residue in a protein. The transverse relaxation time τD is observed only if a fluorophore is embedded in a solid dielectric structure containing water channels or pockets. The solid dielectric may be biologically relevant (the hydrophobic core of a protein) or not (e.g., polyethylene). If the same fluorophore is surrounded by water and there is no solid dielectric, then the TDSS is expected to contain only the τL term. The second question is whether the slow relaxation components with τ > τD in the TDSS of fluorophores in proteins are due to water or due to protein. This study shows that the slow components involve changes in the protein structure, such as a minor adjustment in the tertiary structure or a change in the conformational state of one charged side chain on the surface of the protein. Thus, the rates of the slow components reflect protein dynamics, however, the TDSS contains a part attributable to the motion of protein atoms and a part attributable to the motion of water molecules. Furthermore, negative correlation

Toptygin et al. between water and protein contributions to TDSS is observed for both slow relaxation components in GB1 protein. On the time scale of the 131 ps component, the motion of water molecules shifts Trp emission to the red and the motion of protein atoms shifts it to the blue. On the time scale of the 2.58 ns component, the motion of protein atoms results in a red shift and the motion of water molecules results in a blue shift. The third question is whether the linear-response or the directresponse method should be used to simulate the TDSS in the fluorescence emission of Trp residues in proteins. This study shows that the absolute amplitude of the TDSS obtained using the linear-response method is exaggerated, and this systematic error has a theoretical explanation. In addition, the TDSS obtained using the linear-response method contains high levels of low-frequency noise; random errors of this kind are not found in the TDSS obtained using the direct-response method. The fourth question is whether the rates and the amplitudes of the TDSS obtained in MD simulations are in agreement with the experimental values. This study shows that the amplitudes obtained using the direct-response method are in reasonable agreement with the available experimental data, but some of the rates are not. The relaxation of TIP3P water in MD simulation is 30% faster as compared to experimental data, and the conformational change of Glu-42 side chain is about 33 times faster in the experiment than in the MD simulation. The latter may indicate that the forcefield parameters need some finetuning, however, this isolated finding does not prove that the problem is in the forcefield and it does not provide sufficient information to select the parameter(s) that needs to be adjusted. Additional research into this issue is necessary. Acknowledgment. The authors thank Dr. P. R. Callis for kindly providing the ground-state and excited-state atomic charges for 3-methylindole in water. This research was supported by the National Science Foundation Award MCB-0719248 (to L.B. and D.T.) and the National Institutes of Health Award GM064746 (to T.B.W.). Supporting Information Available: Appendix A, Relationship between spectral shifts and electronic energy levels; Appendix B, Effect of the electric field on the electronic energy levels; Appendix C, Nonlinear Stark effect; Description of the random noise generators that were used to simulate the data for Figure 2; Modified CHARMM22 topology files top_all22_ prot_PCallis_G.inp (with ground-state Trp charges Qj0(0) from Table 1) and top_all22_prot_PCallis_E.inp (with 1La excitedstate Trp charges Qj1(0) from Table 1); Description of a Fortran subroutine Match3D; Corrected emission spectrum of 3-methylindole in hexane. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Time-dependent spectral shift (TDSS) is also known as the timedependent red shift, time-dependent Stokes shift, fluorescence dynamic Stokes shift, and dynamic Stokes shift. As outlined in Appendix A (see Supporting Information), the Stokes shift is either time-invariant or it is a function of two different times. Furthermore, the TDSS can be both the red shift and the blue shift. For these reasons we prefer the term timedependent spectral shift. (2) Ware, W. R.; Chow, P.; Lee, S. K. Chem. Phys. Lett. 1968, 2, 356. (3) Maroncelli, M. J. Mol. Liq. 1993, 57, 1. (4) Jimenez, R.; Fleming, G. R.; Kumar, P. V.; Maroncelli, M. Nature 1994, 369, 471. (5) Horng, M. L.; Gardecki, J. A.; Papazyan, A.; Maroncelli, M. J. Phys. Chem. 1995, 99, 17311. (6) Stratt, R. M.; Maroncelli, M. J. Phys. Chem. 1996, 100, 12981. (7) Gardecki, J. A.; Maroncelli, M. J. Phys. Chem. A 1999, 103, 1187.

Origin of Time-Dependent Spectral Shift in GB1 (8) Kovalenko, S. A.; Schanz, R.; Senyushkina, T. A.; Ernsting, N. P. Phys. Chem. Chem. Phys. 2002, 4, 703. (9) Arzhantsev, S.; Jin, H.; Baker, G. A.; Maroncelli, M. J. Phys. Chem. B 2007, 111, 4978. (10) Sajadi, M.; Obernhuber, T.; Kovalenko, S. A.; Mosquera, M.; Dick, B.; Ernsting, N. P. J. Phys. Chem. A 2009, 113, 44. (11) Brand, L.; Gohlke, J. R. J. Biol. Chem. 1971, 246, 2317. (12) Gafni, A.; DeToma, R. P.; Manrow, R. E.; Brand, L. Biophys. J. 1977, 17, 155. (13) Pierce, D. W.; Boxer, S. G. J. Phys. Chem. 1992, 96, 5560. (14) Riter, R. R.; Edington, M. D.; Beck, W. F. J. Phys. Chem. 1996, 100, 14198. (15) Jordanides, X. J.; Lang, M. J.; Song, X.; Fleming, G. R. J. Phys. Chem. B 1999, 103, 7995. (16) Changenet-Barret, P.; Choma, C. T.; Gooding, E. F.; DeGrado, W. F.; Hochstrasser, R. M. J. Phys. Chem. B 2000, 104, 9322. (17) Vincent, M.; Gilles, A. M.; de la Sierra, I. M. L.; Briozzo, P.; Barzu, O.; Gallay, J. J. Phys. Chem. B 2000, 104, 11286. (18) Pal, S. K.; Mandal, D.; Sukul, D.; Sen, S.; Bhattacharyya, K. J. Phys. Chem. B 2001, 105, 1438. (19) Toptygin, D.; Savtchenko, R. S.; Meadow, N. D.; Brand, L. J. Phys. Chem. B 2001, 105, 2043. (20) Cohen, B. E.; McAnaney, T. B.; Park, E. S.; Jan, Y. N.; Boxer, S. G.; Jan, L. Y. Science 2002, 296, 1700. (21) Pal, S. K.; Peon, J.; Bagchi, B.; Zewail, A. H. J. Phys. Chem. B 2002, 106, 12376. (22) Peon, J.; Pal, S. K.; Zewail, A. H. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 10964. (23) Mataga, N.; Chosrowjan, H.; Taniguchi, S.; Hamada, N.; Tokunaga, F.; Imamoto, Y.; Kataoka, M. Phys. Chem. Chem. Phys. 2003, 5, 2454. (24) Lampa-Pastirk, S.; Beck, W. F. J. Phys. Chem. B 2004, 108, 16288. (25) Qiu, W.; Zhang, L.; Kao, Y.-T.; Lu, W.; Li, T.; Kim, J.; Sollenberger, G. M.; Wang, L.; Zhong, D. J. Phys. Chem. B 2005, 109, 16901. (26) Guha, S.; Sahu, K.; Roy, D.; Mondal, S. K.; Roy, S.; Bhattacharyya, K. Biochemistry 2005, 44, 8940. (27) Toptygin, D.; Gronenborn, A. M.; Brand, L. J. Phys. Chem. B 2006, 110, 26292. (28) Qiu, W. H.; Kao, Y. T.; Zhang, L. Y.; Yang, Y.; Wang, L. J.; Stites, W. E.; Zhong, D. P.; Zewail, A. H. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 13979. (29) Abbyad, P.; Shi, X. H.; Childs, W.; McAnaney, T. B.; Cohen, B. E.; Boxer, S. G. J. Phys. Chem. B 2007, 111, 8269. (30) Halder, M.; Mukherjee, P.; Bose, S.; Hargrove, M. S.; Song, X. Y.; Petrich, J. W. J. Chem. Phys. 2007, 127, 055101. (31) Li, T.; Hassanali, A. A.; Kao, Y.-T.; Zhong, D.; Singer, S. J. J. Am. Chem. Soc. 2007, 129, 3376. (32) Zhang, L.; Wang, L.; Kao, Y.-T.; Qiu, W.; Yang, Y.; Okobiah, O.; Zhong, D. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 18461. (33) Abbyad, P.; Childs, W.; Shi, X.; Boxer, S. G. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 20189. (34) Jesenska, A.; Sykora, J.; Olzynska, A.; Brezovsky, J.; Zdrahal, Z.; Damborsky, J.; Hof, M. J. Am. Chem. Soc. 2009, 131, 494. (35) Zhang, L. Y.; Yang, Y.; Kao, Y. T.; Wang, L. J.; Zhong, D. P. J. Am. Chem. Soc. 2009, 131, 10677. (36) Brauns, E. B.; Madaras, M. L.; Coleman, R. S.; Murphy, C. J.; Berg, M. A. J. Am. Chem. Soc. 1999, 121, 11644. (37) Andreatta, D.; Lustres, J. L. P.; Kovalenko, S. A.; Ernsting, N. P.; Murphy, C. J.; Coleman, R. S.; Berg, M. A. J. Am. Chem. Soc. 2005, 127, 7270. (38) Callis, P. R.; Burgess, B. K. J. Phys. Chem. B 1997, 101, 9429.

J. Phys. Chem. B, Vol. 114, No. 34, 2010 11337 (39) Vivian, J. T.; Callis, P. R. Biophys. J. 2001, 80, 2093. (40) Nilsson, L.; Halle, B. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 13867. (41) Bandyopadhyay, S.; Chakraborty, S.; Balasubramanian, S.; Bagchi, B. J. Am. Chem. Soc. 2005, 127, 4071. (42) Hassanali, A. A.; Li, T. P.; Zhong, D. P.; Singer, S. J. J. Phys. Chem. B 2006, 110, 10497. (43) Golosov, A. A.; Karplus, M. J. Phys. Chem. B 2007, 111, 1482. (44) Li, T.; Hassanali, A. A.; Singer, S. J. J. Phys. Chem. B 2008, 112, 16121. (45) Halle, B.; Nilsson, L. J. Phys. Chem. B 2009, 113, 8210. (46) Muino, P. L.; Callis, P. R. J. Chem. Phys. 1994, 100, 4093. (47) Zheng, Y.; Mamdani, F.; Toptygin, D.; Brand, L.; Stivers, J. T.; Cole, P. A. Biochemistry 2005, 44, 10501. (48) Gallagher, T.; Alexander, P.; Bryan, P.; Gilliland, G. L. Biochemistry 1994, 33, 4721. (49) Gronenborn, A. M.; Filpula, D. R.; Essing, N. Z.; Achari, A.; Whitlow, M.; Wingfield, P. T.; Clore, M. G. Science 1991, 253, 657. (50) Kuszewski, J.; Gronenborn, A. M.; Clore, G. M. J. Am. Chem. Soc. 1999, 121, 2337. (51) Bublitz, G. U.; Boxer, S. G. Annu. ReV. Phys. Chem. 1997, 48, 213. (52) Boxer, S. G. J. Phys. Chem. B 2009, 113, 2972. (53) Pan, C.-P.; Callis, P. R.; Barkley, M. D. J. Phys. Chem. B 2006, 110, 7009. (54) Landau, L. D.; Lifschitz, E. M. Quantum Mechanics (NonRelatiVistic Theory), 3rd ed.; Oxford, England: Pergamon Press, 1977. (55) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1833. (56) Hubbard, J.; Onsager, L. J. Chem. Phys. 1977, 67, 4850. (57) Maroncelli, M.; Fleming, G. R. J. Chem. Phys. 1987, 86, 6221. (58) Castner, E. W.; Fleming, G. R.; Bagchi, B.; Maroncelli, M. J. Chem. Phys. 1988, 89, 3519. (59) Kivelson, D.; Friedman, H. J. Phys. Chem. 1989, 93, 7026. (60) Landau, L. D.; Lifschitz, E. M. Mechanics, 3rd ed.; Oxford, England: Pergamon Press, 1976. (61) We use bold typeface E to denote the electric field vector and light typeface italic E to denote energy. To avoid confusion, the magnitude of vector E is denoted as |E| rather than E. (62) Maroncelli, M.; Fleming, G. R. J. Chem. Phys. 1988, 89, 5044. (63) Gallagher, T.; Alexander, P.; Bryan, P.; Gilliland, G. L. Biochemistry 1994, 33, 4721. (64) Hamilton, W. C. Statistics in Physical Science: Estimation, Hypothesis Testing, and Least Squares; New York: Ronald Press Co., 1964. (65) For space economy, the terms “mean wavenumber” and “mean frequency” are used in this work to denote “Franck-Condon-factor-weighted mean wavenumber” and “Franck-Condon-factor-weighted mean frequency”, respectively. The latter is defined in Appendix A, see eqs A6 and A7 in Supporting Information. (66) A fluorophore cannot have appreciable absorption at its mean65 emission frequency νem,c. Because the probability of emission is proportional to ν3, the mean energy of the emitted photons is always greater than hνem,c, but it cannot be greater then the energy of the absorbed photons, which proves that there are no absorbed photons at νem,c. (67) Kaatze, U. Chem. Phys. Lett. 1993, 203, 1–4. (68) DeLano, W. Pymol Home Page. www.pymol.org (accessed 2009). (69) Pal, S. K.; Peon, J.; Zewail, A. H. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 1763. (70) Toptygin, D.; Wen, X.; Barrick, D.; Brand, L. Biophys. J. 2007, 370A, abstract, 51st Annual Meeting of the Biophysical Society.

JP104425T