Ab Initio Molecular Dynamics of the Green Fluorescent Protein (GFP

We present an ab initio molecular dynamics study of the chromophore of green fluorescent protein (GFP) in the four relevant protonation states. Ground...
8 downloads 0 Views 93KB Size
Download PDF
J. Phys. Chem. B 2001, 105, 5797-5803

5797

Ab Initio Molecular Dynamics of the Green Fluorescent Protein (GFP) Chromophore: An Insight into the Photoinduced Dynamics of Green Fluorescent Proteins Valentina Tozzini* and Riccardo Nifosı` Scuola Normale Superiore and Istituto Nazionale per la Fisica della Materia, Piazza dei CaValieri 7, I-56126 Pisa, Italy ReceiVed: January 4, 2001; In Final Form: April 3, 2001

We present an ab initio molecular dynamics study of the chromophore of green fluorescent protein (GFP) in the four relevant protonation states. Ground-state geometries, electronic structures, and vibrational spectra are calculated. Patterns of vibrations are assigned to recently detected Raman bands. These results provide insight into the correlation between vibrational frequencies and optical absorption wavelengths of the GFP chromophore, observed upon mutations of its environment. An intramolecular mode-coupling mechanism is suggested, which rationalizes the coherent dynamics recently revealed by ultrafast pump-and-probe optical spectroscopy. This mechanism is tested by a simulation mimicking the dynamics after the electronic excitation, which shows a coupling between high and low-frequency modes.

I. Introduction The photoactive green fluorescent protein (GFP) and its mutants have become in the past decades an invaluable tool in molecular and cell biology.1 GFPs are ideal fluorescent markers to monitor protein trafficking in living cells, since their fluorescence occurs without the need of external cofactors and they can be fused to other proteins without affecting their native functions.2 As shown by crystallographic studies,3-6 a β-barrel fold envelops the chromophore, which in wild type (WT-GFP) forms auto-catalytically by posttranslational cyclization of Ser65, Tyr66, and Gly67. The resulting structure, a p-hydroxybenzylindene-imidazolinone, is constituted by the phenolic ring from Tyr66 linked by a bridging carbon to an imidazolinone ring (see Figure 1). Two thermodynamically stable states of the chromophore are photoactive in WT-GFP, one carrying a phenol (A) and the other a phenolate (B), which are responsible for the two bands at 397 and 477 nm of the absorption spectrum, respectively. The equilibrium between those two states is influenced by external factors such as temperature and pH, and by mutations that affect the environment of the chromophore.2,4 Spectroscopic measurements7,8 on single GFP molecules revealed peculiar on-off blinking and photobleaching dynamics, which result in an undesired reduction of the fluorescence. To explain the spectral properties and photodynamics of GFP, the following protonation forms of the chromophore have been proposed in the literature: the neutral form (shown in Figure 1), the anionic form (deprotonated on the phenolic oxygen O14), the cationic form (protonated on the heterocyclic nitrogen N1), and the zwitterionic form (deprotonated on O14 and protonated on N1). The attribution of A and B states to protonation forms of the chromophore has been debated.10,11 However, most authors now agree in assigning A and B to the neutral and anionic states, respectively, since to date no experimental evidence for protonation of the ring nitrogen has been found.9 The zwitterionic form was instead recently invoked to explain the blinking behavior,11 being prone to nonradiative conversion through cis-trans isomerization. Few experiments and, up to the present date, no theoretical calculation have been published on the vibrational dynamics * Corresponding author. Email: [email protected]..

Figure 1. Sketch of the GFP chromophore in the neutral state. The anionic form is obtained by deprotonating O14, the cationic form by protonating N1, and the zwitterionic form by deprotonating O14 and protonating N1.

of GFP chromophore.12 FTIR spectroscopy measurements were performed to detect the change in the protonation state of the chromophore and possibly of its environment,13 by probing changes in vibrational properties of GFP upon photoconversion. Recently, the vibrational spectra in the 1000-1660 cm-1 region have been solved by Raman spectroscopy for WT and S65T GFP mutants and for a model-chromophore molecule.14 The authors observed a linear correlation between a specific vibrational frequency (in the range 1610-1655 cm-1) and the optical absorption energy in different model-chromophore states and GFP mutants. They postulate that a shift in the Raman bands may reflect a change in the stabilization of the ground state, which may in turn dominate the optical absorption shifts. However, the nature of these vibrational modes was not determined, preventing an unambiguous interpretation of these observations. Even more recently, a coupling of the electronic excitation to the low-frequency bands of the vibrational spectrum (500-600 cm-1) was detected by ultrafast pump-and-probe optical spectroscopy on a GFP mutant.15 Motivated by the lack of satisfactory interpretations of the previously mentioned experimental findings, in this paper we present an extensive ab initio study of the dynamical properties of the four protonation forms of the GFP chromophore, based on density functional theory (DFT) and Car-Parrinello molecular dynamics (CPMD). We shall assign the bands revealed by Raman spectroscopy to specific vibrational patterns, and give an interpretation for the observed correlation between vibrational

10.1021/jp010052q CCC: $20.00 © 2001 American Chemical Society Published on Web 05/25/2001

5798 J. Phys. Chem. B, Vol. 105, No. 24, 2001 frequencies and optical absorption energies. A mechanism for the intramolecular mode coupling preliminarily reported in ref 15 is extensively discussed and verified by CPMD simulation. Before discussing the vibrational properties, we also report on the optimized structures and ground state and first excited state electronic structures of the chromophores. II. Model Chromophores and Computational Methods The model neutral chromophore is shown in Figure 1. The connections with the protein backbone (at N2 and C3) are cut and saturated with hydrogen atoms. This is the simplest model chromophore preserving the essential features for the fluorescence. Although the chromophore environment is known to produce sensitive frequency shifts in the electronic excitation energies, calculations with larger model chromophores10 did not show any qualitatively significant additional feature in the electronic structure. It is known from X-ray crystallographic measurements4-6 that O14, O13, and N1 are involved in a hydrogen bond network with the amino acids in the environment. We expect that the vibrational modes localized near these sites and their geometry will be significatively affected by the presence of the environment. By contrast, the electrostatic interaction due to charged or polar residues should have a weaker effect. These points will be the subject matter of further investigations, the present paper being focused on the structure and vibrational properties of the isolated chromophore. Within DFT frame, we evaluated the exchange and correlation functional in the local density approximation (LDA). We used the Perdew and Zunger parametrization of the quantum Monte Carlo exchange and correlation energy,16,17 and included gradient corrections in the functional form proposed by Becke and Perdew.18,19 This approximation scheme has proven accurate in describing the structural and vibrational properties of small conjugate systems like retinals.20 The interaction between valence electrons and inner cores is described by soft-first principles pseudopotentials.21 We used a plane wave basis set for the single-particle Kohn-Sham wave functions, with a 25 Ryd cutoff, which ensures that the bond lengths are within 1% from the full convergence value.20 All the calculations were performed in a cubic cell with a 15 Å lattice parameter and the chromophore was oriented in such a way as to ensure a minimum distance above 10 Å between periodic images. Car-Parrinello molecular dynamics simulations (CPMD)22 were performed using the electronic mass preconditioning scheme23 and a time step of 0.145 fs. The vibrational frequencies were extracted from the velocity autocorrelation function on runs of about 1.5 ps. The normal mode coordinates were obtained by the standard FT technique and the multiple signal classification algorithm.24,25 Additional hints for localized modes were obtained from the spectral analysis in local coordinates. To have a more detailed description of the normal modes, in one case (the anion, see below) we calculated the dynamical matrix within density functional perturbation theory (DFPT)26 using plane waves and the same energy cutoff, cell parameters, and pseudopotentials of the CP simulation. This technique gives the full normal mode decomposition in the harmonic approximation. However, since we shall investigate the effects of coupling and anharmonicity, our study will be mainly based on the CPMD simulations. III. Ground-State Properties Starting from a planar geometry, we optimized the structure of each protonation state of the chromophore in the configuration N1-cis with respect to the C5-C6 bond.

Tozzini and Nifosı` TABLE 1: Bond Lengths in the Optimized Geometries of the Model Chromophorea theory

experiment

bond (Å)

neutral anionic zwitterionic cationic neutral anionic

C4-C5 C5-C6 C6-C7 C7-C8 C7-C9 C8-C10 C9-C11 C10-C12 C11-C12 N1-C3 N1-C5 N2-C3 N2-C4 C4-O13 C12-O14

1.484 1.362 1.438 1.412 1.409 1.380 1.384 1.394 1.395 1.301 1.406 1.383 1.406 1.247 1.382

1.455 1.387 1.404 1.429 1.428 1.364 1.365 1.446 1.444 1.304 1.405 1.385 1.413 1.264 1.284

1.445 1.397 1.393 1.437 1.440 1.356 1.357 1.453 1.449 1.334 1.403 1.344 1.441 1.248 1.269

1.463 1.369 1.420 1.417 1.417 1.371 1.377 1.403 1.398 1.330 1.410 1.332 1.447 1.232 1.359

1.476 1.366 1.457 1.401 1.394 1.375 1.385 1.379 1.374 1.315 1.418 1.369 1.380 1.227 1.385

1.47 1.37 1.40 1.43 1.43 1.38 1.36 1.45 1.47 1.32 1.40 1.43 1.39 1.22 1.28

a Lengths are given in angstroms. Experimental data for the isolated neutral chromophore and for the anionic chromophore embedded in the protein environment are taken respectively from refs 27 and 28 (Protein Data Bank code 1YFP). The C-H, N-H, and O-H bond lengths are, respectively, ∼1.09, ∼1.02, and ∼ 0.99 Å.

The resulting bond lengths are reported in Table 1, together with data from the X-ray structure of the isolated neutral chromophore27 and the X-ray structure of a yellow GFP mutant (S65G/H148G/T203Y, pdb code 1YFP),28 prevalently found in state B (anionic). In the case of the anionic chromophore the comparison is not straightforward due to the influence of the environment on the bond lengths. Moreover, the geometry of the chromophore within different solved GFP structures has a certain degree of variability (deviations of about 0.05 Å for corresponding C-C and C-N bonds, and even more for C-O bonds). The following observations can be made: (i) The geometry of the neutral form compares fairly well (within 1% on average) with the experimental data. In the anionic form the agreement is less accurate, especially concerning the the C-N bonds, probably due to the stretching effect of the backbone chain, which is not included in the present calculations. (ii) The bondorder schemes that can be deduced by the bond lengths for the different protonation states (see also Figure 2) are in agreement with those inferable from experimental data. In particular, the phenolic ring in the neutral and cationic forms has an aromatic character, which is lost in the anion and zwitterion, where there are four single and two double bonds. Moreover, the bridging carbon makes one double bond and one single bond in the neutral and cationic forms, while it has two partial double bonds in the anion and the zwitterion. (iii) the C12-O14 bond length reflects the protonation state of O14 and reproduces the experimental data, while the C4-O13 distance is overestimated in the anion. Since strong H-bond interactions are present between O13 and residues Q94 and R96,3-6 and due to the low accuracy of crystal data, no conclusions can be taken on the nature of this discrepancy. A search for other stable structures with the simulated annealing technique reveals the presence of nonplanar local minima for zwitterionic and cationic chromophores. In these structures the dihedrals N1-C5-C6-C7 and C5-C6-C7C8 depart about 10° from the planar value. The energy deviations from the planar structures are too small (order of hundredth of eV) to be considered significative.29 These local minima can reasonably be attributed to the sterical repulsion between hydrogens of N1 and C8. In fact, in none of the X-ray structures of GFPs there is evidence for the presence of

Molecular Dynamics of the GFP Chromophore

J. Phys. Chem. B, Vol. 105, No. 24, 2001 5799 IV. Electronic Structure

Figure 2. HOMO, LUMO, and bond-order schemes for the four protonation states of the chromophore. The bond-order schemes are deduced by a qualitative analysis of the bonding-antibonding character of the orbitals.

nonplanar forms, which may be ascribed either to the actual absence of the cationic and zwitterionic forms in the GFP at equilibrium, or to the effect of chromophore environment in weakening the steric repulsion. At a temperature of about 1800 K, we observed a rotation around the single bond C6-C7 of the neutral chromophore, leading to the same configuration, a part from the direction of the hydrogen linked to O14 (the two configurations do not differ in energy appreciably). We also optimized the N1-trans configuration of the neutral chromophore with respect to the C5dC6 bond, which turns out to be planar (although the C5-C6-C7 angle is larger than that in the cis case) and higher in energy of about 2.4 kcal/mol. In the following we restrict our study to the planar-cis forms, which are thermodynamically stable within the protein environment.

DFT-LDA correctly describes the variation of the gap between HOMO (highest occupied molecular orbital)-LUMO (lowest unoccupied molecular orbital) as a function of structural changes,20,30 even though it is known to substantially underestimate its absolute value. The energy gaps obtained in the present calculation are the following: neutral ) 2.21 eV, anionic ) 1.78 eV, zwitterionic ) 1.58 eV, cationic ) 1.92 eV. These values are within 60-65% of recent optical absorption data on a model chromophore compound in aqueous solution14 (neutral ) 3.37 eV, anionic ) 2.89 eV, cationic ) 3.15 eV; the zwitterionic form is not reported in the experiment). These findings show that the variation of excitation energy with protonation state is correctly reproduced.31 Figure 2 reports the HOMO and LUMO charge density distributions, showing their π character. In each case, a charge transfer from the phenolic to the imidazolinone ring is observed upon excitation, which can lead to a change in the proton affinity of the protonation sites.32 In agreement with previous theoretical studies,10 we find that, upon excitation, the proton affinity of N1 increases in the anionic and in neutral chromophores, while that of O14 decreases in the anionic and zwitterionic states. In particular, the change in proton affinity of N1 may have a relevant role in the blinking/bleaching processes, which may involve moieties protonated on N1. The internal charge redistribution upon excitation is schematically represented by the bond-order schemes in Figure 2, which are deduced by a qualitative analysis the bonding-antibonding character of LUMO with respect to HOMO. The bond orders of the ground state (left column of Figure 2) are consistent with the bond-lengths of the optimized structures (Table 1). In the excited states (right column) they give a qualitative indication of the change in bond strength upon excitation. In systems with conjugated chains containing double-single bond alternation, the main effect of the excitation is to increase the strength of single bonds and to weaken double bonds,20,33 in some cases reverting the bond-order alternation. In this case, the presence of the rings imposes geometrical constraints, making the effects of the excitation less straightforward: i) the C4-C5 and N1-C5 bonds are strengthened, while the N1-C3 and N2-C4 and C4-O13 are weakened; (ii) when the bridge contains a double and a single bond (i.e., in the neutral and cationic forms) an inversion of the bond order can be observed in the excited state. In the other forms, the effect of excitation is to weaken the bridging bonds; (iii) in the anionic and zwitterionic states, the aromatic character of the phenolic ring is enhanced after excitation, due to weakening of double bonds. The observations made are relevant to the discussion of the vibrational dynamics correlated with electronic excitation (see next section). V. Vibrational Properties Figure 3 shows the vibrational spectra obtained from runs of free dynamics at about 300 K.34 The total number of vibrational modes is 60 in the neutral and zwitterionic chromophores, 57 in the anion, 63 in the cation, and they appear with different spectral weight depending on the temperature and on the intrinsic spectral strength. A list of modes significantly contributing to the spectrum is given in Tables 2 and 3 together with a characterization35 based on normal-mode analysis (see section II). Modes with similar displacement patterns are reported on the same line. For comparison, vibrational modes revealed by Raman spectroscopy14 are shown as arrows and reported in the table. Experimental frequencies are assigned to the closest theoretical ones. Because of the low degree of symmetry of the

5800 J. Phys. Chem. B, Vol. 105, No. 24, 2001

Tozzini and Nifosı`

Figure 3. Vibrational spectra of the chromophores. Lines: Fourier transforms of the velocity autocorrelation function from a microcanonical molecular dynamics run at T = 300 K. Arrows: Raman data from 14. The accurate comparison between theoretical and experimental data is reported in Tables 2 and 3.

TABLE 2: Vibrational Frequencies Corresponding to Hydrogen Stretching Modes, Extracted from CPMD Simulationsa neutral

anionic

3585 3430

3445

3035 3060 3030 3020 3010 2955

3030 3015 3005 2985 2945 2965

zwitterionic

cationic

description

3435 3425 3110 3030 3015 2990 2980 2985

3575 3450 3415 3160 3100 3045 3010 2975 2940

O14-H N2-H N1-H C3-H Cph-H Cph-H Cph-H Cph-H C6-H

a Frequencies are reported in cm-1. C are the carbons in the phenolic ph ring. In the “description” Column, we have reported the bonds which are mainly involved in the stretching.

molecule, a detailed description of each mode would be lengthy.35 However, we will discuss in detail the modes that have some relevance in the physical behavior of the protein. Hydrogen Modes. The modes lying between 2800 and 3700 cm-1 are the C-H, N-H and O-H stretching modes (Table 2). Although we find a certain degree of mixing, in particular among the N-H and the C-H in the phenolic ring, these modes can be described in terms of local coordinates, i.e., the bond lengths. The differences between the frequencies of corresponding modes can be attributed to variations in the chemical environment of the atoms bonded to the hydrogen involved in the mode. For example, a comparison of the frequencies of C3-H shows that the chemical environment of C3 is similar in the anionic and neutral chromophore, and in cation and zwitterion separately. We remark that the interaction with the environment will have considerable influence on hydrogen stretching modes, especially that localized on O14-H, which is known to be involved in hydrogen bonds.3-6 Modes in the Region 1000-2000 cm-1. The 1000-2000 cm-1 region (Table 3) is characterized by modes of stretching

vibration of C-C, C-O, and C-N bonds and contains the vibrational bands revealed by Raman spectroscopy (Bell et al.14). The agreement between theoretical and experimental frequencies is within 5%. Even though, to avoid photodegradation problems, the experiments are not strictly performed on resonance, the assumption can be made that the modes associated with the electronic excitation are selectively excited. Thus, since the primary effect of the excitation is to change the strength of the bonds between heavy atoms (see previous section), it is not surprising that the modes detected by Raman spectroscopy are stretching modes (see Figure 5 and below). This behavior is common to other photoactive molecules with a delocalized π system, such as retinals.33 A separate discussion must be devoted to the first of these modes, which corresponds to the C4-O13 stretching vibration. To ensure that in our CPMD simulations we did not miss any mode in the gap around 1700-3000 cm-1, we calculated the harmonic dynamical matrix within DFPT for the anionic chromophore, obtaining the full decomposition in normal modes. As a byproduct, we were able to estimate the effect of anharmonicity in a softening of 1-2% depending on the mode (data not shown35). The C4-O13 stretching mode turns out to be the highest immediately below the 1700-3000 cm-1 gap and can be confidently assigned to the highest frequency of the experimental spectra. By measuring the frequency of this mode in different GFP mutants and in a model chromophore compound in different solvents, Bell et al.14 found that it linearly correlates with the optical absorption energy. Since the vibrational spectrum probes the ground-state properties, they explain the observed correlation in terms of a change in the stabilization of the ground state and deduce that this mode is delocalized over the whole chromophore, being sensitive to variations in the environment. According to our calculations, only a weak degree of delocalization is present, due to the mixing with collective stretching vibrations of the C-C and C-O bonds. However, our results rationalize Bell et al. hypothesis, suggesting the following picture (schematically sketched in Figure 4): the lower frequency of the C4-O13 stretching mode in the anionic form with respect to the neutral form (1630 vs 1675 cm-1) is due to the weaker character of the C4-O13 bond. This is, in turn, determined by the degree of electronic delocalization around C6-C5-C4-O13, which is higher in the anionic with respect to the neutral form (see Figure 4). Correspondingly, the energy gap decreases in the anion (1.78 eV vs 2.21 eV in the neutral form), giving a possible explanation of the correlation found in ref 14. The chromophore environment will, in general, weaken the C4-O13 bond, since possible donors from Q94 and R96 are found in the X-ray structures at hydrogen bonding distance to O13.3-6 This explains why the frequency of the C4O13 mode is lower in the GFP mutants than in the model chromophores14. Moreover, hydrogen bonds with O13 tend to delocalize the electronic system around C6-C5-C4-O13, decreasing the energy gap. The overall effect is that the optical transition energy and the vibrational frequency of the C4-O13 bond are correlated. Because of the low accuracy on the energy gaps, it is not appropriate to compare the linear correlation parameters inferable from our results with those found in ref 14. However, our calculations indicate that the moieties protonated on N1 (zwitterion and cation) do not fit on the same linear relation as the neutral and anionic forms. Indeed, the cation is found out of the correlation by Bell et al. In principle, this fact could be used to detect the presence of absorbing states protonated on N1,

Molecular Dynamics of the GFP Chromophore

J. Phys. Chem. B, Vol. 105, No. 24, 2001 5801

TABLE 3: Frequencies (in cm-1) Extracted from CPMD Simulationsa neutral theor

anion expt

1675 1615 1585

1645 1603 1567

1525 1475 1415 1330

1525 1448 1329

1265 1225

1257 1226

1145

1179

1075 1025

1098

theor 1630 1605 1555 1510 1490 1445 1405 1330 1310 1245 1145(sh) 1125 1060(sh) 1020

960 910 895

960

850 825

840

zwitterion expt 1630 1577 1556 1535 1501 1439 1315 1263 1229 1172 1095

1670 1635 1575 1520 1495 1475 1375 1365 1310 1265 1230 1145 1105 1065 1040

940

510(w)

425 385 290 280(sh) 240 220(sh) 190

theor

expt

1000-2000 cm-1 1725 1651 1595 1585 1584 1540 1539 1510 1455 1452 1355 1290

1350 1270 1225

1155

1182

1090

1077

0-1000 cm-1 975 940

895 880

770 730 710 660(w)

cation

theor

910 850

825 810 785 750 670 610 575 505 470 460 430 425 385

785 730 710 660 580 540 475 450

820 780 710 660 615 560 505 440 (w) 410

430 360 280

270 210

description C4 d O13, C4=O13 CdC, C=C, C-O CdC, C=C, C-C, C-O CdC, C=C, C-C, C-O NdC, C=C C-C C-C C-N C-N, C-C H-rocking, ang def coll str + ang def +H rocking coll str + ang def +H rocking coll str + ang def +H rocking coll str + ang def +H rocking coll str + ang def +H rocking coll str + ang def +H rocking planar def of phen + str of C-C trans def of phen planar def of imidaz trans def of phen planar ang def of the rings trans def of phen ring planar def of the rings planar ang def of rings trans ang def of rings trans ang def of rings trans ang def of rings planar ang def of phen + trans def of imidaz planar ang def planar ang def phen + coll motion trans def trans def planar ang def phen + coll motion and rock of O13 coll planar def + rock O14 coll trans coll trans def coll trans def coll planar def coll trans def coll planar def coll trans def

a A gap line is left where the full decomposition in normal modes from the dynamical matrix indicates the presence of one or more modes which are lacking in the vibrational spectra. In the “description” column, we have reported the bonds which are mainly stretched in the mode (), = , meaning double, partial double and single bond, respectively) or an indicative description (def ) deformation, phen ) phenolic ring, imidaz ) imidazolinone ring, ang ) angular, coll ) collective, trans ) transverse, str ) stretching, rock ) rocking). Cph are the carbons in the phenolic ring. Data labeled with (w) or (sh) are present as weak peaks and shoulders, respectively. Additional information is available on request.

once the vibrational properties were measured together with the optical properties. The three subsequent modes cannot be easily described in terms of local coordinates, since they are delocalized along the principal axis of the molecule. In particular, the second (1590 cm-1 for neutral) is the most delocalized and corresponds to a vibration along the axis of the molecule, involving stretching of the bonds O13-C4-C5-C6-C7, of C8-C10 and C9-C11 in the phenolic ring and of C12-O14. The pattern of vibrations of this mode best matches the displacements induced upon electronic excitation (see sec. IV), explaining why it mostly contributes to the experimental spectra14 (1567 cm-1 in the neutral chromophore). Modes below 1000 cm-1. The modes in this region are angular (planar or transverse) deformations of the rings, and collective wavelike (planar or transverse) vibrations (Table 3). In ultrafast pump-and-probe optical spectroscopy experiments on a GFP mutant15 a (ground state) vibrational band of the anionic chromophore was revealed at 593 cm-1. This technique

is able to detect vibrational modes coupled to the electronic excitation. However, this frequency does not correspond to a stretching mode, and a coupling mechanism between vibrational modes must be invoked to explain its excitation. We propose that the electronic excitation directly couples to stretching modes which in turn are coupled to lower energy angular modes. To verify this hypothesis we performed a CPMD run for the anionic chromophore (which is the stable state for the specific mutant studied in ref 15) simulating the dynamics following electronic excitation. The starting configuration was chosen in order to mimic the pattern of displacement produced upon electronic excitation. This was achieved by scaling each bond length by a factor reflecting the change of bond strength in the excited state (see Figure 2).36 The physical meaning of this procedure is to put the system in the so-called “hot state” in which it is found after excitation.37 The followed procedure is nonstandard and, to a certain extent, arbitrary. However, the results shown in Figure 5 are elucidating, and confirm the coupling hypothesis. Only a limited number of the sixty modes are selected and, in

5802 J. Phys. Chem. B, Vol. 105, No. 24, 2001

Tozzini and Nifosı` of the simulation correctly mimicked the “hot molecule” situation. Moreover, associated with the stretching modes, we find two groups of modes around 800 cm-1 and around 500 cm-1, which correspond to planar angular deformations of the rings. These results confirm the coupling of the stretching vibrations to these low-frequency angular deformations, and explains the 593 cm-1 band detected in the experiment.15 This coupling is presumably due to the geometric constraints introduced by the rings. In a very recent theoretical work38 it has been pointed out that the vibrational mode coupling is an important mechanism for energy transfer and dissipation in proteins. By means of a procedure similar to ours, in that paper the authors give excess energy to a specific vibrational mode coupled to the photodissociation of heme group in myoglobin, and follow the rather complex energy transfer mechanism by a Classical molecular dynamics simulation. VI. Conclusions

Figure 4. Proposed scheme for the correlation between C4-O13 vibrational frequency and optical absorption energy. In the neutral case the electronic charge is localized on the C5-C6 and C4-O13 bonds, while in the neutral form it is distributed around C6-C5-C4-O13. The degree of delocalization around these atoms simultaneously influences the energy gap and determines the strength of the C4-O13 bond.

Figure 5. Vibrational spectrum from the 1.5 ps MD simulation starting from a configuration chosen in order to mimic the effect of electronic transition (see text). As expected, stretching modes (1000-2000 cm-1) are excited, which correspond to the modes detected by Raman spectroscopy (short arrows).14 Due to coupling with stretching modes, lower frequency modes corresponding to angular deformation of the rings are also present. One of these can be assigned to the band observed in ultrafast pump-and-probe spectroscopy (longer arrow).15

particular, some modes in the stretching region are enhanced. As found in the Raman data,14 the theoretical peaks are grouped in two separate sets, at frequencies corresponding to the experimental ones. This ensures that the starting configuration

In this paper we have accurately studied the structural, electronic and vibrational properties of GFP chromophore in the four relevant protonation states by means of DFT based techniques. A charge transfer takes place within the chromophore upon optical excitation, which induces a change in the bond strengths, exciting the vibrational modes detected by Raman spectroscopy. These modes can be identified and attributed to localized or collective stretching deformations, on the basis of our simulations. In particular, the C4-O13 vibration is assigned to the highest frequency of the experimental spectra.14 We suggest that the linear correlation between this frequency and the optical absorption energy found in ref 14 is due to the fact that the C4-O13 bond strength is a “fingerprint” of the electronic delocalization around C6-C5-C4-O13. Hydrogen bonds between O13 and the amino acids in the environment will additionally delocalize the electronic system simultaneously reducing the energy gap and the C4-O13 vibrational frequency. Since this correlation depends also on the protonation state of N1, a joined measurement of the vibrational and optical properties could in principle detect the presence of chromophore states protonated on N1 (zwitterionic or cationic). Our results also show that stretching modes excited by the electronic transition are coupled to lower frequency modes corresponding to angular deformations of the rings. The coupling in the present system is due to geometrical constraints imposed by the presence of the rings. These findings explain the presence of a band at 593 cm-1 observed by ultrafast pump-and-probe optical spectroscopy.15 Moreover, the optical excitation induces a proton affinity change in the protonation site of the imidazolinone ring, which may be important in the blinking and photobleaching behavior of GFPs. MD studies including the chromophore environment will provide further insight in GFP photophysics, in particular regarding the dynamics connecting bright and dark states. Acknowledgment. We thank Francesco Buda, Paolo Giannozzi, and Jorge Kohanoff for making available the codes for Car Parrinello molecular dynamics, DFT calculations, and normal-mode analysis and for useful discussions and technical support. References and Notes (1) Sullivan, K. F.; Kay, S. A. Green Fluorescent Proteins; Academic Press: San Diego, 1999.

Molecular Dynamics of the GFP Chromophore (2) Cubitt, A.; Heim, R.; Adams, S.; Boyd, A.; Gross, L.; Tsien, R. Y. Trends Biochem. Sci. 1995, 20, 448. (3) Yang, F.; Moss, L. G.; Phillips, G. N., Jr. Nat. Biotechnol. 1996, 14, 1246. (4) Brejc, K.; Sixma, T. K.; Kitts, P. A.; Kain, S. R.; Tsien, R. Y.; Ormo¨, M. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 2306. (5) Ormo¨, M.; Kubitt, A. B.; Kallio, K.; Gross, L. A.; Tsien, R. Y.; Remington, S. J. Science 1996, 273, 1392. (6) Palm, G. J.; Zdanov, A.; Gaitanaris, G. A.; Stauber, R.; Pavlakis, G. N.; Wlodawer, A. Nat. Struct. Biol. 1997, 4, 361. (7) Dickson, R. M.; Cubitt, A. B.; Tsien, R. Y.; Moerner, W. E. Nature 1997, 388, 355. (8) Cinelli, R. A. G.; Ferrari, A.; Pellegrini, V.; Tyagi, M.; Giacca, M.; Beltram, F. Photochem. Photobiol. 2000, 71, 771. (9) Elsliger, M.; Wachter, R. M.; Hanson, G. T.; Kallio, K.; Remington, S. J. Biochemistry 1999, 38, 5296. (10) Voityuk, A. A.; Michel-Beyerle, M. E.; Ro¨sh, N. Chem. Phys. 1998, 231, 13. Voityuk, A. A.; Michel-Beyerle, M. E.; Ro¨sh, N. Chem. Phys. Lett. 1997,272, 162. Voityk, A. A.; Michel-Beyerle, M. E.; Ro¨sh, N. Chem. Phys. Lett. 1997, 296, 269. (11) Weber, W.; Helms, V.; McCammon, J. A.; Langhoff, P. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 6177. (12) During the submission of this paper, a referee brought to our attention the following related paper: Langhoff et al. Chromophore Charge States and the Proton Shuttle Mechanism in Green Fluorescent Protein: Inferences Drawn from Ab initio Theoretical Studies of Structures and Infrared Absorption Spectra. J. Phys. Chem. B 2001. In press. (13) van Thor, J. J.; Pierik, A. J.; Nugteren-Roodzant, I.; Xie, A.; Hellingwerf, K. J. Biochemistry 1998, 37, 16915. (14) Bell, A. F.; He, X.; Wachter, R. M.; Tonge, P. J. Biochemistry 2000, 39, 4423. (15) Cinelli, R. A. G.; Tozzini, V.; Pellegrini, V.; Beltram, F.; Cerullo, G.; Zavelani-Rossi, M.; Tyagi, M.; Giacca, M. Phys. ReV. Lett. 2001, 86, 3439. (16) Perdew, J. P.; Zunger, A. Phys. ReV. B 1981, 23, 5048. (17) Ceperley, D. M.; Alder, B. J. Phys. ReV. Lett. 1980, 45, 566. (18) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (19) Perdew, J. P. Phys. ReV. B 1986, 33, 8822. (20) Bifone, A.; de Groot, H. J. M.; Buda, F. J. Chem. Phys. B 1997, 101 2954. Bifone, A.; de Groot, H. J. M.; Buda, F. Chem Phys. Lett. 1996, 248 165. (21) Vanderbilt, D. Phys. ReV. B 1990, 43, 7892.

J. Phys. Chem. B, Vol. 105, No. 24, 2001 5803 (22) Car, R.; Parrinello, M. Phys. ReV. Lett., 1985, 55, 2471. (23) Tassone, F.; Mauri, F.; Car, R. Phys. ReV. B 1994, 50, 10561. (24) Kohanoff, J. Comput. Mater. Sci. 1994, 2, 221. (25) Marple, S. L. Jr. In Digital Spectral Analysis with Applications; Prentice-Hall: New Jersey, 1985. (26) Giannozzi, P.; de Gironcoli, S.; Pavone, P.; Baroni, S. Phys. ReV. B 1991, 43, 7231. (27) Kurimoto, M.; Subramony, P.; Gurney, R. W.; Lovell, S.; Chmielewsky, J.; Kahr, B. J. Am. Chem. Soc. 1999, 121, 6952. (28) Wachter, R. M.; Elsliger, M.-A.; Kallio, K.; Hanson, G. T.; Remington, S. J. Structure 1998, 6, 1269. (29) Energy differences of 0.05 eV are usually considered at the limiting confidence level in the present approximation. Moreover, the stabilization/ destabilization introduced by the environment and by the vibrational entropy should be considered. (30) Saitta, A. M.; Buda, F.; Fiumara, G.; Giaquinta, P. V. Phys. ReV. B 1996, 53, 1446. (31) We remark that the HOMO-LUMO analysis is a rather correct description of the electronic excitation in the case of GFP chromophores, since the HOMO-LUMO contribution to the transition from ground state to first excited is 70-90%.10 (32) The change in proton affinity (PA) can be roughly estimated via a Fo¨rster cycle, according to which the variation of PA of S upon excitation is PA(S*) - PA(S) ) gap(S) - gap(HS), were S and HS are the protonated and deprotonated state, respectively. (33) Curry, B.; Palings, I.; Broek, A. D.; Pardoen, J. A.; Lugtenburg, J.; Mathies, R. In AdVances in Infrared and Raman Spectroscopy; Clark, R. J. H., Hester, R. E., Eds.; Heyden, London, 1985; Vol. 12. (34) The thermalization is achieved by adopting the usual procedure of randomizing initial conditions (for example, see: 20). (35) Additional information is available on request. (36) Changing the scaling factors by 10-20% did not produce any qualitative change in the spectrum. Any accurate choice of these factors would be inappropriate, since the description of the change in bond strengths due to the electronic excitation is only qualitative. (37) Blanchet, V.; Zgierski, M. Z.; Seideman, T.; Stolow, A. Nature 1999, 401, 52. Concerning the time length of our simulation, we observe that a time scale of 1.5 ps does not allow a complete thermalization, implying that the “hot-state” situation is described. (38) Moritsugu, K.; Miyashita, O.; Kidera, A. Phys. ReV. Lett. 2000, 85, 3970.