Operation of the Proton Wire in Green Fluorescent Protein. A Quantum

Apr 9, 2008 - Departament de Química Universitat Autònoma de Barcelona, 08193 Bellaterra (Barcelona), Spain, Institut de Biotecnologia i de Biomedic...
7 downloads 0 Views 652KB Size
5500

J. Phys. Chem. B 2008, 112, 5500-5511

Operation of the Proton Wire in Green Fluorescent Protein. A Quantum Dynamics Simulation Oriol Vendrell,† Ricard Gelabert,*,‡ Miquel Moreno,‡ and Jose´ M. Lluch‡,§ Departament de Quı´mica UniVersitat Auto` noma de Barcelona, 08193 Bellaterra (Barcelona), Spain, Institut de Biotecnologia i de Biomedicina, UniVersitat Auto` noma de Barcelona, 08193 Bellaterra (Barcelona), Spain, and Theoretische Chemie, Physikalisch-Chemisches Institut, UniVersita¨t Heidelberg, Im Neuenheimer Feld 229, 69120 Heidelberg, Germany ReceiVed: February 8, 2008

A nuclear quantum dynamical simulation of the proton shuttle operating in the green fluorescent protein has been carried out on a high-quality, high-dimensionality potential energy surface describing the photoactive ππ* excited state, and including motion of both the three protons and of the donor and acceptor atoms of the hydrogen bonds in a closed proton wire. The results of the simulations show that proton transfer along the wire is essentially concerted, synchronous, and very fast, with a substantial amount of the green fluorescent species forming within several tens of femtoseconds. In this regard, analysis of the population of the fluorescent species indicates that at least two dynamical regimes are present for its formation. Within the first hundreds of femtoseconds, dynamics is very fast and impulsive. Later on, a slower pace of formation appears. It is discussed that the two largest decay times for the protonated chromophore reported experimentally (Chattoraj, M.; King, B. A.; Bublitz, G. U.; Boxer, S. G. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 8362-8367) might correspond to some irreversible process occurring after formation of the fluorescent species, rather than to cleavage of the chromophore’s phenolic O-H bond.

1. Introduction Few biochemical systems can claim to have experienced the exponential growth of interest as green fluorescent protein (GFP) in the last 15 years. This protein, naturally occurring in North Pacific jellyfish Aequorea Victoria, stands out because of its intense bright-green fluorescence.1 This fact would be a mere curiosity were it not for the discovery that its chromophore, p-hydroxybenzylideneimidazolinone,2 is formed after an autocatalytic intramolecular cyclization process affecting a wellknown hexapeptide sequence.3,4 Any protein whose genetic transcription is enlarged with that of GFP develops GFP-like fluorescence properties upon expression in living cells. This has opened up vast numbers of applications: Nowadays GFP is a well-known biological marker in fields so diverse as molecular biology, cell biology, and medicine because the presence and concentration of proteins can be monitored by measuring the fluorescence due to the GFP part of the target protein.1 Much attention has been placed on the analysis of its structure in order to understand the causes of the peculiar photophysical behavior of GFP. The amino-acid sequence of wild-type GFP (wt-GFP) was first established 15 years ago,5 and shortly afterward the crystallographic structure of wt-GFP was finally solved.6 Nowadays the Protein Data Bank hoards scores of different GFP structures, including mutants and different crystallizations. Even though several of these mutants display markedly different spectral properties, they are remarkably similar structurally. In short, GFP apparently can occur as either a dimeric or monomeric protein, even though the photophysical * Corresponding author. E-mail: [email protected]. † Universita ¨ t Heidelberg. ‡ Departament de Quı´mica, Universitat Auto ` noma de Barcelona. § Institut de Biotecnologia i de Biomedicina, Universitat Auto ` noma de Barcelona.

properties do not apparently depend on this fact. Structurally, GFP can be envisaged as a rigid cylinder or barrel formed by 11 β sheets, inside of which an R-helix traverses diagonally holding the chromophore roughly at the center of the cavity. The detailed crystallographic information allowed the determination of the precise surroundings of the chromophore in wtGFP (Cro). It is accepted nowadays that the chromophore is hydrogen bonded to an internally caged water molecule in the cavity (Wat25), which in turn is bonded in the same manner to a serine residue (Ser205), which finally connects to a glutamate residue (Glu222).7 Other nearby residues contribute to stabilize the chromophore (see Figure 1 for details). How does the structure correlate to the spectroscopic properties displayed by GFP? wt-GFP presents strong absorption at 397 nm (band A) and a smaller absorption band at 477 nm (band B), whose intensity depends on pH (being six times smaller at pH ) 8.0). Fluorescence emission occurs at 510 nm with high quantum yield (0.72-0.85) and a weak unstructured fluorescence peak around 460 nm.8-10 It has been proposed that the absorptions at 397 and 477 nm have their origins in neutral and ionized forms of the phenolic group of the chromophore.4 Interesting data was also harvested from the study of specific mutants of wt-GFP. The Ser65Thr (S65T) mutant does not display the absorption peak at 397 nm while the absorption at 477 nm is increased sixfold and red-shifted.11 Comparison of the protein structures of monomeric wt-GFP7 and the monomeric S65T mutant12 revealed differences due only to the presence of the extra methyl group in the 65th residue or lack thereof. Specially revealing was the fact that the S65T chromophore was deprotonated. Because of the sharp contrast in the spectra, the absorption peak at 397 nm was then ascribed to the neutral form of the chromophore, whereas the absorption at 477 nm to its ionized form.7

10.1021/jp801169z CCC: $40.75 © 2008 American Chemical Society Published on Web 04/09/2008

Proton Wire in Green Fluorescent Protein

Figure 1. Schematic representation of the prevailing interactions involved in the GFP photocycle. W are water molecules, and the dotted lines depict hydrogen bond interactions. Labels in brackets identify the location of amino acids before the autocatalytic cyclization reaction that results in the synthesis of the chromophore.

Figure 2. Representation of the three-state model presumed to operate within GFP: A indicates a neutral chromophore, and I1, I2, and B are different forms including a ionized chromophore. See text for a detailed description of the structures.

The excited-state dynamics of GFP has been studied through time-resolved fluorescence.13 From this study a three-state model was proposed, which can be summarized as follows: The A speciessessentially a neutral chromophorescan be photoexcited at 397 nm to species A*, which evolves very rapidly to an intermediate species I* with decay times in the range of a few picoseconds. I* can either decay to the ground state I (and later revert to A) or it can evolve further to B*, eventually relaxing to B. Interconversion of A and B is slow in the ground state. Assuming existence of the intermediate I/I* forms was necessary because A* was known to disappear very quickly but the formation of B* from A* was also known to be slow.13 This model is summarized in Figure 2. Especially significant was the discovery of a very large kinetic isotope effect (KIE) of about 5 on H/D substitution of all exchangeable protons, which naturally suggested proton transfer as the process in the excited state. Accordingly, species A (absorbing at band A) is believed to contain a neutral chromophore, whereas species I and B (the latter absorbing at band B) are supposed to contain the chromophore in anionic form.

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5501 Species I and B must then differ in the environment of the chromophore.13,14 Similar results were obtained through picosecond-resolved fluorescence and femtosecond absorption experiments by other authors.15 The existence of a significant normal KIE with H/D substitution for the nonradiative decay of A* is in agreement with proton-transfer reactions in the excited state (ESPT). On the basis of the structure around the chromophore of wt-GFP and its mutants, a nonradiative deactivation channel was proposed for species A* involving the operation of a three-proton relay, starting with the phenolic proton of the chromophore, and including the captive water molecule Wat25, and the Ser205 and Glu222 residues.7,14 The product would be a non-relaxed form of the I* intermediate (henceforth I/1) having a deprotonated chromophore. I/1 would evolve to a relaxed form I/2, which also contains a deprotonated chromophore, through the loss of hydrogen bonding between Ser205 and Glu222 and the change of conformation of the (just transferred) proton on Glu222 from anti to syn conformation, presumed to be more stable.7,16 Finally, a certain conformational change was expected to occur to turn I/2 into B*. This change implies the rotation of residue Thr203 to stabilize the charged chromophore.7,14 Green fluorescence from the chromophore would then be from the I/2 form. See Figure 3 for details. Evidence of the A* f I* proton pathway and the nature of the I* state has finally come from transient infrared absorption spectroscopy, which has allowed to detect protonated Glu222 after photoconversion.17-19 This mechanism and interpretation is nowadays the most widely accepted. It is not the only proposal there is: Proton wires in GFP have been shown to be more extended than previously known: for instance, there is evidence pointing to a proton wire linking Glu222 near the chromophore with Glu5 on the protein surface.20 On the basis of the asymptotic behavior of fluorescence, Agmon has proposed a mechanism involving a conformational change enabling the rotation of Thr203, which eventually allows the proton to escape to the exterior solution.21,22 So far, the most accepted explanation of the fluorescent properties of GFP, its mechanism and intermediates, has received much attention from theoretical scientists. Much interest has been placed on possible rapid deactivation pathways of the excited state due to transit through non-adiabatic crossings in GFP and its mutants.23-25 Specific studies on the proton relay energetics and dynamics are more scarce.16,26,27 A dynamical simulation of the operation of the proton wire would provide insight into matters such as feasibility of operation of this mechanism, its rate, and nature (that is, whether it is concerted or stepwise, and in such a case, what is the order of motion). Theoretical work by Wang and Smith on the proton chain transfer in the ground electronic state indicates that the reaction coordinate involves sequential motion of the protons.28,29 However, proper dynamical study is extremely complicated in this case for several reasons: (1) the basic system is very large, (2) after photoexcitation the system must evolve in an electronic excited state where transfer to a different electronic state cannot be disregarded,27 and (3) it involves motion of very light atoms, likely to display strong quantum effects (for instance, tunneling). To our knowledge, to date there is only one such dynamical study due to Lill and Helms.16 In their study, a nearly classical molecular dynamics method was used, in which a modified force field to describe the excited state of the chromophore and some corrections to account for non-classical behavior of the transferring protons were employed. The main conclusions were that the proton-wire operation is triggered by

5502 J. Phys. Chem. B, Vol. 112, No. 17, 2008

Vendrell et al.

Figure 3. Schematic representation of the operation of the proton wire in GFP (A* h I/1) and further evolution of the system in the photoactive state. The simulations presented in this work comprise only the A* h I/1 step.

the proton transfer from the chromophore to the captive water molecule Wat25 and that afterward the remaining proton transfers would occur very rapidly, within some tens of femtoseconds through a barrierless process. They also concluded that the conformational change of Glu222 from the anti to syn conformation did not occur within the simulation time (100 ps). However, in their simulation the proton transfer from the chromophore to Wat25 was forced; that is, the starting point of the simulation was an already deprotonated chromophore: This ruled out the possibility to explore a concerted three-proton transfer. Recently, we have published a detailed study of the dynamics of GFP in the electronic ground state based on sophisticated molecular mechanics simulations.26 In this study, it was shown that for long times of simulation the proton wire formed by the chromophore, Wat25, Ser205, and Glu222 was well-formed and remained approximately planar without appreciable deviations. On the basis of these results, a thorough quantum chemical study of a reduced model that included modelizations of all residues in the proton wire was done at the Complete Active Space SelfConsistent Field (CASSCF) level of theory including corrections up to second-order perturbation theory (CASPT2), within Cs symmetry. A large set of energies was obtained for structures where all donor-acceptor distances and the position of protons were varied (See Figure 4). Even though the effect of the protein environment could not be included in this model, the study represented one of the most sophisticated approaches to the photochemistry of GFP, provid-

Figure 4. Reduced model used in this work: The chromophore is substituted by p-hydroxybenzylideneimidazolinone, Ser205 by a methanol molecule, and Glu222 by an acetate.

ing a consistent picture of the potential energy landscape in the photoactive electronic state (we point out nonetheless that certain proposals have been made to include the effect of the protein environment in single-point multiconfigurational calculations30). Our study ruled out the possibility of electronic state crossing between the photoactive ππ* and a dark πσ* state whose energetic proximity was detected in a previous study with a model that included only the chromophore and water molecules.27 From this static picture, prospective mechanisms of

Proton Wire in Green Fluorescent Protein

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5503

Figure 5. Scheme depicting the dynamical model (colored atoms). The potential energy has been computed for structures that include the colored atoms and the groups in parentheses (see Figure 4). The labels at the top of the figure refer to the chemical groups in the real system. The distance between donor and acceptor atoms is represented by Ri (i ) 1...3, depending of the hydrogen bond). The proton coordinates ri are defined as the (signed) distance from proton i to the point halfway between the donor and acceptor atom: in the arrangement shown in the Figure, all ri values are negative.

operation of the proton wire could be investigated on a purely static basis. The most favorable in energetic terms required a sequence of motion of the protons that seems counterintuitive: first transfer of the Ser205 proton, followed by transfer of the Wat25 proton, and finally the proton of the chromophore. However, a hypothetical concerted motion process for the three protons would present a potential energy barrier only 4 kcal mol-1 higher at the same level of theory. Hence, the difference between a stepwise and concerted mechanism could not be established in the absence of real motion of the protons and all of the quantum effects associated with it. The large number of structures studied at the CASPT2 level of theory, spanning a sizable part of the configurational space describing the operation of the proton wire in this model of the GFP opens up the way to comprehensive quantum dynamical simulations of the operation of the proton wire. In this paper, we take on the task of carrying these simulations out and thus provide valuable insight into the short-time dynamics on the ππ* electronic excited state of the proton wire within GFP using state-of-the-art theoretical tools. With this, we will simulate the evolution of the real system from the A* species to the I/1 species, which comprises, strictly speaking, the operation of the proton relay. 2. Computational Details Detailed dynamical simulations of the proton wire are necessary in order to elucidate the questions that are the target of this study. Special attention should be placed on the fact that the proton wire includes heavy atoms whose motion certainly will modulate the dynamics of the wire but especially that it contains three protons. Such light particles are at the borderline between classical and quantum behavior, and explicit consideration of their quantal nature is a must in any balanced simulation. Even though other approaches have been used for GFP in which approximate methods have been employed,16 we think that a detailed and theoretically sound quantum dynamical simulation on a reasonably precise potential energy surface is needed to address the issue of sequence of motion and rate of transfer of the protons in the wire. In what follows, we will provide a description of the dynamical model, the potential energy surface, the kinetic energy operator, and the methodology used to integrate the timedependent Schro¨dinger equation.

2.1. Proton-Wire Model. The purpose of this study is to carry out a detailed study of the operation of the proton wire as described in the Introduction. The questions for which answers are sought concern the position of the protons and bonded atoms over time. Given the sheer size of the system, it is necessary to make a model: the number of atoms exacts a heavy toll on all electronic quantum calculations of the electronic excited-state energies, and the number of degrees of freedom (geometrical variables) considered explicitly brings about an exponential growth in the number of single-point calculations needed to represent the potential energy, and besides, on the computational costs incurred in determining the time evolution of the system. This work builds up on the results of a previous published paper26 and thus the chemical model of the proton wire is the same, depicted in Figure 4. The quantum chemical model includes all atoms explicitly appearing in the GFP proton wire, while the side chains out of the wire itself are modeled to reduce computational cost. As such, Ser205 has been substituted by a methanol, and Glu222 by an acetate. On the basis of previous results,26 the proton wire in the chemical model has been kept linear and the rest of the side chains to conform to global Cs symmetry. 2.2. Dynamical Model. The specification of the dynamical model consists of elliciting the exact choices of geometrical variables whose time evolution is of interest, along with the specification of the interactions that control the system in the form of its potential energy. The proton wire in GFP consists of three transferrable protons that shift positions from a donor oxygen atom to an acceptor oxygen atom. We draw the attention of the reader to the similarities of this system to that used by one of us to successfully describe proton translocation along a chain of water molecules and capture quantum effects therein.31 Therefore, this dynamical model has also been applied in the current work. Figure 5 displays a scheme of the dynamical model and variables used. As in the work in ref 31, proton and heavy atom motion is restricted to be monodimensional, which is an approximation that simplifies notably the simulation (effectively halving the number of degrees of freedom at least and simplifying notably the Hamiltonian Vide infra). As defined, the dynamical variables R1, R2, and R3 are the distances between each donor-acceptor atom pair, starting at the chromophore and going all the way to glutamate. r1, r2, and r3 represent the signed distance between

5504 J. Phys. Chem. B, Vol. 112, No. 17, 2008

Vendrell et al.

each proton and the point halfway between its donor and acceptor atoms. This specific choice is convenient because it makes it easy to detect when a specific proton has been transferred by simple monitorization of the sign of its coordinate. As it stands, the dynamical model comprises six independent degrees of freedom. The specification of the dynamical model also requires the definition of the interaction governing the system, that is, the potential energy field as a function of all geometrical variables of the system. For practical reasons, the potential energy cannot be given as a set of energy values (this would require over 54 million single-point energy calculations) but rather must be a function ready to be computed for any set of geometrical variables of the system. A total of 338 distinct geometries have been computed at the CASPT2 level of theory with the same active space of six electrons in six molecular orbitals used in our previous work,26 in which all donor-acceptor distances have covered the range between 2.3 and 2.8 Å, and the three protons have explored all possible values. The geometries computed sample all 12 possible single-proton motion processes that can happen in the GFP proton wire.26 As for the rest of the atoms in the model, their relative positions have been held fixed at their values in the minimum in the ground electronic state, S0. When oxygen atoms were displaced along the wire (change of a Ri coordinate), the potential energy was computed for a structure where the distance between the oxygen atoms and their supporting residues was kept fixed (that is, the residues were traslated alongside with the oxygen atom with fixed relative coordinates). Using this procedure, the potential energies of the ground electronic state and the photoactive ππ* state have been obtained. In order to make the computation of the potential energy for a given set of geometrical parameters possible, we fitted an interpolating eight-state Empirical Valence Bond (EVB)32 potential energy function of the six geometrical variables VS0 to the energy values of the electronic ground state, and analogously another function VS1 has been fitted to the energies of the photoactive ππ* excited state. Each EVB state corresponds to a different protonation state of the wire. A similar approach has been used previously to represent proton transfer in excited electronic states.33,34 In the present work, each potential energy function is an explicit function of the 6 geometrical variables R1, R2, R3, r1, r2, and r3, and each depends on a different set of 55 parameters (110 parameters total). The process of fitting of the potential energy surface and its analysis will be described elsewhere.35 2.3. Choice of Dynamics. The fact that the atoms that traslocate in the wire are protons makes it necessary to account for the possibility of quantum behavior. To this end, we have chosen to carry out quantum dynamical simulations. Thus, it is necessary to solve the time-dependent Schro¨dinger equation

ip

∂ Ψ(q,t) ) H ˆ Ψ(q,t) ∂t

(1)

where q is a vector defining the configuration of the system and H ˆ is the Hamiltonian (total energy) operator. The first step consists of defining the Hamiltonian operator. The kinetic energy operator in use in this work is derived from that published by one of us.31 The kinetic energy operator reads

( (

Tˆ ) -

) ( ) (

)

1 1 1 ∂2 1 ∂2 + + 2 2mH1 4M ∂r 2mH2 4M ∂r 2 1 2 2

) )

1 ∂ 1 ∂ ∂2 ∂2 1 + + + 2mH3 4M ∂r 2 M ∂R 2 ∂R 2 ∂R 2 3 1 2 3

( (

2

1 ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + 2M ∂r1 ∂R2 ∂r2 ∂R3 ∂r2 ∂R1 ∂r3 ∂R2

) (

)

∂ ∂ 1 ∂ ∂ ∂ ∂ 1 ∂ ∂ + + + (2) 4M ∂r1 ∂r2 ∂r2 ∂r3 M ∂R1 ∂R2 ∂R2 ∂R3 where mHi is the mass of hydrogen isotope i and M is the mass of the oxygen atom. In doing this it is assumed that during the operation of the proton wire the residues that form it remain essentially static and only the donor and acceptor oxygens move. The potential energy operators have been described in Section 2.2. Quantum dynamical calculations in the six-dimensional model presented above have been performed by means of the multiconfiguration time-dependent Hartree (MCTDH) method.36-39 In particular, the Heidelberg MCTDH package of programs has been used.40 This method has been applied recently and successfully in different aspects of multidimensional intramolecular proton-transfer systems.31,41-43 A brief description of the method is presented here. The MCTDH method is a general algorithm to solve the timedependent Schro¨dinger equation. The MCTDH wave function is expanded in a sum of products of so-called single-particle functions (SPF). The SPFs φ (Q,t) may be one- or multidimensional functions and, in the latter case, the coordinate Q is a collective one, Q ) (qk,...,qj). Because the SPFs are timedependent, they follow the wave packet over time and often a rather small number of SPFs suffices for convergence. The Ansatz for the MCTDH wave function reads

Ψ(q1,...,qf,t) t Ψ(Q1,...,Qp,t) ) n1

np

∑j ∑j ...

1

p

p

Aj1,..., jp

φj (κ)(Qκ,t) ∑ κ)1 κ

(3)

where f denotes the number of degrees of freedom and p is the number of MCTDH particles, also called combined modes. There are nκ SPFs for the κth particle. The equations of motion for the coefficient vector A and for the SPFs are derived from a variational principle. It is important to note that MCTDH uses, in a variational sense, optimal SPFs as this ensures the fastest convergence. The equations of motion are complicated, but because there are comparatively few equations to be solved the MCTDH method can be very efficient. In the concrete simulations that are presented in this work, the MCTDH wave function is expanded as a series of 27 000 time-dependent Hartree products (number of independent Aj1,...,jp coefficients in eq 3). Except in the relaxation dynamics, unitarity and energy conservation have been monitored to guarantee accuracy of the simulation. Naturally enough, because the time propagation of the dynamics is done and the complete information on the wavepacket is stored on disk, it is possible to compute the expectation values of different operators (that is, observables) of use in interpreting the dynamics. Relevant details on these operators will be discussed in the main text as needed. 2.4. Simulation Protocol. If the dynamical study presented in this work is to be meaningful and resemble what occurs in the proton wire in GFP, then the process of photoexcitation has to be simulated reasonably. A traditional approach in quantum

Proton Wire in Green Fluorescent Protein

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5505

Figure 6. Expectation values of the dynamical coordinates describing the positions of the transferrable protons Hi (ri) (left) and donor-acceptor atom distances (Ri) (right), for the perprotio wire simulation. The inset in the left panel figure details the first 50 fs of simulation. For definition of dynamical coordinates see Figure 5.

dynamics simulations mimicks the Franck-Condon transition by promoting the ground vibrational state of the system (i.e., the most populated) in S0 to the photoactive state. It is possible to obtain the ground vibrational state by determining the time evolution in imaginary time (i.e., relaxation) with the S0 Hamiltonian H ˆ ) Tˆ + Vˆ S0, until invariance in the (imaginary) time-evolved state, in which case the ground vibrational state is obtained:

t f -i∞ Ψ(q) f Ψ0(q)

(4)

Needless to say, energy is not conserved in relaxation dynamics. Photoexcitation is then simulated by placing the ground vibrational state, Ψ0(q), on S1, that is, solving the timedependent Schro¨dinger equation (eq 1) with initial condition Ψ0(q) and Hamiltonian H ˆ ) Tˆ + Vˆ S1. To analyze the results of the dynamical propagation, we computed the expectation values of some observables along the dynamics. If one is interested in obtaining the expectation value of an observable X, then the following quantity must be computed

〈X〉 ) 〈Ψ(t)|Xˆ |Ψ(t)〉

(5)

where Xˆ is the quantum mechanical operator associated to observable X, and the wave function Ψ(t) is normalized. 3. Results 3.1. Dynamics of the Perprotio and Perdeutero Proton Wires. The ground vibrational state for the perprotio variant of the proton wire has been prepared through the procedure outlined in Section 2.4. This led to a vibrational state where all protons were located near the respective donor atoms, with an energy of 13.8 kcal mol-1 above the potential energy minimum of S0. Then this state was promoted to S1 and propagated in time for 1.5 ps. To analyze the dynamics, the expectation values of the position of each proton in the wire, 〈ri〉(i ) 1...3), and the distance between each donor-acceptor atom pair, 〈Ri〉(i ) 1...3), have been computed along time. The results are shown in Figure 6 for the initial 500 fs of the simulation. Immediately after photoexcitation a sudden, impulsive motion is detected in the donor-acceptor distances, especially in R1, which corresponds to the Cro-Wat25 hydrogen bond, which changes from 2.55 Å (as in the vibrational ground state in S0) to about 2.40 Å in only about 20 fs. A comparable effect has been described previously for other ESPT reactions, for instance from solvatochromic shifts of photoacids.44,45 In the present case, this strong change turns later into damped oscillations that stabilize to values of about 2.50 Å in the long time limit. R2

and R3 do not behave analogously: if at all in the initial stages of the simulation they tend to increase their values slightly, even though some recurrences can be seen in both. Overall, the effect of photoexcitation on the wire skeleton can be described as a marked shortening of the donor-acceptor distance in the CroWat25 hydrogen bond. At longer times, R2 increases to a value of 2.6 Å, describing a certain loosening of the structure of the proton wire, which ends up composed of two units: Cro + Wat25 and Ser205 + Glu222. The effects of photoexcitation on the location of the protons are visible in the left panel of Figure 6: The same impulsive motion experienced by R1 occurs here, and it can be seen that all protons evolve in a more or less concerted way from reactants (that is, proton closer to the donor atom, or negative ri values) to products. This displacement is fast: each proton is closer to its acceptor atom within 20 fs of photoexcitation. Of the three protons, the one shared between Cro and Wat25 lags slightly behind the other two. If the fluorescent species corresponds to the anionic chromophore, then its formation is somewhat delayed. The plots of ri show a general trend toward increasing values, but also display some rapid oscillations of small amplitude. The origin of these “vibrations” already in the product region can find an explanation in one of the features of the model used. As described, the dynamical model is closed: once the wave packet reaches the product region, where the potential energy is lowest, it has no other possibility but to bounce back and revert in the direction leading toward reactants. In the real system there would be other possibilities open, such as activation of other degrees of freedom, or even other more complicated structural changes like the breaking up of the wire through displacement of residue Glu222. Because this is not possible in our model, the wave packet bounces back and interferes with the incoming waves, thus giving rise to this oscillatory pattern in the plots of ri. The long-time behavior of the expectation value of the proton position shows a drift toward more positive values in Figure 6, which describe situations where the corresponding proton is increasingly bonded to the acceptor atom. Figure 7 shows the same data for the full extent of the simulation. It can be seen that after the first 500 fs the average values of each curve stabilize, even though not all at the same value. The process, as far as the final position of the protons is concerned, reaches its end at times around 0.6 ps. Complete isotopic substitution of protons by deuterons, the perdeutero proton wire, has also been done to explore the changes brought about by isotopic substitution on the dynamics. The ground vibrational state appears now 10.7 kcal mol-1 above

5506 J. Phys. Chem. B, Vol. 112, No. 17, 2008

Vendrell et al.

Figure 7. Values of the dynamical coordinate describing the position of the protons in the wire (left) and the distances between donor-acceptor atoms for the perprotio simulation in the long time limit.

Figure 8. Values of the dynamical coordinate describing the position of the protons in the wire for the perdeutero simulation. The inset details the first 50 fs of simulation.

the potential energy minimum. Photoexcitation triggers in this case qualitatively the same sequence of events as in the perprotio case, even though a certain slowdown is detected. Figure 8 displays the time evolution of the expectation values of the positions of the deuterons in the wire. If the kinetic isotope effect was measured with respect to the transfer of the first proton(deuteron), because its transfer brings about the fluorescent anionic chromophore, then the ratio of times at which r1 ) 0 for the perprotio and perdeutero proton wires is approximately 2. 3.2. Discussion on the Concertedness or Asynchronicity of Proton Motion. In our preceding work, the different possibilities for a complete proton transfer along a series of three proton transfers was examined from a static standpoint.26 The most favorable stepwise proton transfer and a hypothetical transition-state structure for the three-proton concerted motion differed only by about 4 kcal mol-1. Even within a static picture, it was clear that a concerted motion of the protons could not be ruled out and that dynamics had to be considered to provide an answer to this question. Insofar the results for the perprotio and perdeutero proton wire have displayed a very marked concerted behavior. Even though not all protons move at exactly the same rate, nor do they move exactly at the same time, their motion happens essentially in a concerted fashion. Interesting information can be obtained from other mixed isotopic substitutions on the wire. In Figure 9 the expectation values for the ri dynamical variables are presented for the proton wires where the first two protons have been turned into deuterons (DDH), and the proton wire where the third proton has been substituted by a tritium atom (HHT).

In the HHT case, for instance, where the proton between Ser205 and Glu222 has been substituted by a tritium atom, the values for 〈r3〉 differ appreciably in the long time limit and present differences at short time. However, H1 and H2 have not been changed and, despite this, their dynamics present noticeable differences: 〈r1〉 first crosses the zero line at about 15 fs in the perprotio wire (see Figure 6), whereas now does it after 20 fs. A similar picture comes up in the DDH proton wire. In this case, the proton shared between Ser205 and Glu222 has not been changed, but, despite this, 〈r3〉 is noticeably different in the perprotio and DDH cases. Thus, changing a proton’s mass in the wire affects the dynamics of other protons not involved directly in the isotopic substitution. In this regard, the picture of a concerted motion of the protons gets stronger support: at least it now seems clear that the motion of the three protons is strongly correlated. The concertedness of the process that has been found is in contrast to the sequential hopping mechanism established through ultrafast IR measurements of ESPT in water.46 This different behavior can be explained considering two factors: First, a sequential motion of the protons in GFP initiated by the chromophore’s proton would induce important charge separation in the intermediate species that, as pointed out in our previous work,26 implies high-energy reaction paths. When considering aqueous solutions the intermediate and charged species can be stabilized easily. Second, in aqueous solutions two consecutive proton hops are probably separated by a requisite molecular reorientation. This is not the case in GFP, where the proton wire is held approximately fixed in space in an arrangement and relative orientation that eases proton transfer without such rearrangements. However, see ref 47 for a new interpretation of the experiments on ESPT in water, which suggests a mechanism that is more akin to that proposed in this work. 3.3. Time Evolution of Proton-Wire Species. The question of what is the chemical path followed by the system upon photoexcitation remains open: that is, what possible species are present at a given time after photoexcitation, where species means in this context any of the possible protonation states of the wire. As discussed elsewhere,26 there are eight different protonation states for the proton wire, which correspond to the different possibilities open to each proton to have been transferred or not. In what follows, we use the same naming convention as in ref 26, where each protonation state of the wire is referenced by a three-letter code, and in which the first, second, and third letters denote the residue supporting the first,

Proton Wire in Green Fluorescent Protein

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5507

Figure 9. Values for 〈ri〉 for the DDH (left) and HHT (right) isotopic variants of the proton wire.

TABLE 1: Mathematical Description of the Pˆ Operator Used to Analyze the Population of Each of the Eight Protonation States in the Proton Wirea protonation state

operator Pˆ

CWS CWG CSS CSG WWS WWG WSS WSG

[1 - Θ(r1)][1 - Θ(r2)][1 - Θ(r3)] [1 - Θ(r1)][1 - Θ(r2)][Θ(r3)] [1 - Θ(r1)][Θ(r2)][1 - Θ(r3)] [1 - Θ(r1)][Θ(r2)][Θ(r3)] [Θ(r1)][1 - Θ(r2)][1 - Θ(r3)] [Θ(r1)][1 - Θ(r2)][Θ(r3)] [Θ(r1)][Θ(r2)][1 - Θ(r3)] [Θ(r1)][Θ(r2)][Θ(r3)]

a Θ(x) is the Heaviside function, yielding 1 when x > 0 and zero when x < 0. In this definition of Pˆ , we make use of the fact that the ri coordinates are positive only in the product area of the proton transfer, and negative in the reactants’. For the naming convention of the protonation states, see the main text.

second, and third proton of the wire, respectively, according to the following key: C, Cro; W, Wat25; S, Ser205; G, Glu222. To determine the population of a certain protonation state, we need to add up the density in the area of configurational space describing that protonation state. In practical terms, this means computing the expectation value of operator Pˆ (using eq 5) where the form of this operator is given in Table 1. This calculation has been done for the eight protonation states of the perprotio and perdeutero wire. The results are displayed in Figure 10. Generally speaking, the behaviors of the perprotio and perdeutero proton wires are similar. Initially there is a very rapid decrease in the population of the CWS protonation state from 1 to about 0.1 in 30 fs in the perprotio proton wire (∼80 fs in the perdeutero case). After 100 fs the population of the WSG, the final product of the operation of the wire, can be seen increasing steadily. The other six intermediate protonation states are always present, but it is remarkable that most of them are only at very low levels (< 0.1), and among these the most abundant ones are those where the first proton is still bound to the chromophore. The most abundant intermediate species is actually CSG, where the second and third protons have already been transferred. Accepting that what is experimentally measured with the fluorescence of GFP is the appearance of the anionic chromophore, a direct comparison to experiment should be possible if the population of protonated chromophore species is computed. Using the nomenclature used in our previous work26 and in Figure 10, the fraction of A* would correspond to adding up the populations of species CXY, where X can be either W or S, and Y can be either S or G. The result depicts the amount of

protonated chromophore still present at a given time, irrespective of the state of the second and third proton wires. The results are displayed in Figure 11. It becomes quite apparent that there are two different regions or dynamical cycles in Figure 11, and actually in all previous figures where dynamical variables have been depicted. There is an initial phase of short-time dynamics, covering about 200 fs after photoexcitation, in which a fast decay of the population of the neutral chromophore occurs, and a long-time dynamics, whose onset happens after 200-300 fs characterized by a slowly decreasing population of CXY with many small-amplitude recurrences of high frequency. Focusing on the short-time dynamics, it can be seen that the first minimum in the curve depicting the population of CXY species appears around 20 fs for the perprotio and about 30 fs for the perdeutero species of the proton wire. This corresponds in our model to the shortest time in which the A* species (protonated chromophore) evolves to I/1 (deprotonated chromophore). Later on, recurrences of small amplitude and high frequency set on until the end of the simulation. Within this initial phase of the dynamics, a clear isotopic effect of about 1.5 appears for the H/D substitution. This value has led us to postulate that the motion in the short-time dynamics phase is essentially impulsive because of its short duration and the closeness of the H/D isotope effect to that expected from free motion of both isotopes: (mD/mH)1/2 ) 1.41. Thus, for the perprotio proton wire the initial formation of anionic chromophore occurs very quickly within ∼20 fs. Because oxygen atoms do actually move, the isotope effect ought to be smaller than this value. The fact that such a large value is found points in the direction that some quantum tunneling might be present. Alternatively, the long-time dynamics is characterized by a slowly decreasing CXY population with recurrences of small amplitude and high frequency. A visual inspection of Figure 11 reveals that in the perprotio proton wire the period of such recurrences is approximately 40 fs. This matches well with the motion of the protons in the wire from reactants to products, and then back in the direction of reactants, which would be about 2 × 20 ) 40 fs, according to our analysis of the motion in the short-time dynamics phase. Moreover, this period matches that seen in Figure 7 for the 〈ri〉 dynamical variables, which represent the position of the proton in each hydrogen bond. Accordingly, the recurrences in the population of CXY in the long-time dynamics phase of the perprotio proton wire must correspond to bouncing motion of the protons from the product region toward the reactants’ area and back.

5508 J. Phys. Chem. B, Vol. 112, No. 17, 2008

Vendrell et al.

Figure 10. Time evolution of the population of each of the eight different possible states of protonation of the proton wire for the perprotio (left) and perdeutero (right) cases. For the naming convention of the protonation states, see main text.

Figure 11. Time evolution of the population of CXY protonation species, corresponding to all possible protonations of the proton wire where the chromophore is still protonated, for the full dynamical situation (left) and first 100 fs (right). Results are presented for the perprotio (red) and perdeutero (green) versions of the proton wire.

These recurrences are essentially much damped, both in amplitude and in frequency, in the perdeutero proton wire. We think that this is in agreement with the previous explanation of the recurrences in the perprotio proton wire. Deuterons move initially slower than protons in the impulsive phase of the dynamics (30 fs). In both cases, as the first products are formed, the wave packet reaches the end of the products area and, because the model is closed, must bounce back. Along the process, part of the energy the wave packet had in the incoming trip mostly in the directions representing proton motion ends up transferred to other heavy modes (in our model that can be only the donor-acceptor modes, Ri). As a result, the receeding wave packet retains a fraction of the energy it had when arriving in the products’ region. Then the wave packet travels back to the reactants, and a different behavior can be expected for the proton and deuteron: the deuteron, being twice as massive as the proton, has less chances to penetrate into the reactant region. The proton, a much more quantal particle, can penetrate deeper because of its higher zero-point energy and its unmatched capability to enter classically forbidden areas. In the end, the freedom of motion of deuterons once in the product region is constrained to a larger extent than that of protons, showing then less recurrences. We can now discuss the behavior of the simulation in the long-time limit. The process in S1 describing the operation of the wire is exoergic and after the initial stages of the simulation are over is mainly localized in the product region (that is, containing a deprotonated chromophore). However, because of

the larger restrictions on deuteron motion, these are more localized in the product region than protons are. This is what is seen in Figure 11, where, as can be seen, the perdeutero species ends up transformed more extensively into products than the perprotio, yielding for our model and simulation time an inVerse KIE for the H/D substitution in the long time range. The higher amount of CXY species at long times for the perprotio proton wire must originate from the sum of systems evolving from reactants to products and systems evolving from products to reactants: otherwise it becomes difficult to explain that in the initial stages the perprotio proton wire reacts faster, thus having smaller amounts of CXY species than the perdeutero. In this regard, we point out that one of the features of the dynamical model in use is that it is closed and cannot describe any ulterior evolution of the wire (such as its breaking). In this sense, WSG represents the most stable minimum in our model, which includes the system up to the formation of the I1 species (see Figure 3). In Figure 10 it can be seen that after the first 200 fs about 50% of the systems correspond to this species. In other words, at long times of simulation, and because there is no possibility in our model for further evolution once the WSG species is formed, our dynamics brings about characteristics more akin to the onset of equilibrium. Several indicia seem to indicate that the proton wire in GFP could break up after proton transfer,7 and some authors have pointed out that the proton-wire network in GFP is actually more extensive than expected and could imply that the shuttled proton does not end its trip at the Glu222 residue, but much farther.20-22

Proton Wire in Green Fluorescent Protein

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5509

Figure 12. Time evolution of the population of CXY protonation species, corresponding to all possible protonation states of the perprotio (left) and perdeutero (right) proton wire, where the chromophore is still protonated. Total time of simulation is 15 ps.

If provisions were made to take into account an irreversible fate for the proton wire in the model (such as breaking of the structure of the wire or escape of the shuttled proton), then we expect that the faster-moving protons would have a more readily accessible exit than the slower-moving deuterons. This exit is missing in the current model. Thus, it is to be expected that the inclusion of a slower and irreversible step would favor the disappearance of WSG species in the perprotio proton wire to a larger extent than in the perdeutero one. In so doing, the amount of remaining WXY species would be substantially reduced and likely show a normal KIE. 3.4. Comparison to Experiment. Maybe the most relevant measurement available related to the dynamics of the proton wire is that of the rise times of the fluorescence of the I* species and the decay of the A* species, determined by Chattoraj et al.13 In the perprotio experiment, the protein was excited at 398 nm (A form) and fluorescence was measured at 460 nm (A* form) and at 508 nm (I* form). The fluorescence responses were fitted to multiexponential decay/rise functional forms. For the decay of A* a triexponential fit was obtained with amplitudes A1 ) 0.50 ( 0.03, A2 ) 0.40 ( 0.05, and A3 ) 0.10 ( 0.01, and decay times τ1 ) 3.6 ( 0.3 ps, τ2 ) 12.0 ( 1.0 ps, and τ3 ) 120.0 ( 15.0 ps. For the rise of I*, only a biexponential fit was reported, with amplitudes A1 ) -0.24 ( 0.02 and A2 ) -0.43 ( 0.03, and rise times τ1 ) 2.2 ( 0.2 ps and τ2 ) 8.1 ( 0.5 ps. A corresponding set of data for the perdeutero proton wire was also determined.13 The process of disappearance of A* was irreversible and has been interpreted as taking place in the range of several picoseconds. As mentioned in the Introduction, there has been a previous theoretical study of the dynamics of the proton wire due to Lill and Helms.16 The study done by Lill and Helms was based on quasiclassical molecular dynamics of the proton wire on a modified potential energy field, which included the effect of the protein environment. In their study, and owing to the difficulty in assessing theoretically the energy barrier for release of the chromophore’s proton in the excited state, the simulations were started after forcing the proton transfer of the chromophore (in our nomenclature, their starting structure would be WWS). Under these conditions, they found that the proton transfer from Wat25 to Ser205 occurred within 10 fs, and the proton transfer from Ser205 to Glu222 within 10-80 fs. Appropriately they defined the transfer of the second and third protons in the wire as ultrafast and concluded that the experimentally measured decay times of the excited-state species A* (at 460 nm) of 3.6

and 12.0 ps 13 ought to correspond to the cleavage of the O-H bond of the chromophore in the excited state. The current study consists of a six-dimensional quantum dynamics simulation of the operation of the proton wire, based on a CASPT2-quality potential energy surface on a model of the proton wire, including motion of heavy atoms. Our simulation clearly supports the fact that the second and third protons in the wire transfer in an ultrafast manner: in Figure 6 the expectation values of the coordinates of all protons are in the product region within 20 fs. However, our simulations do not restrict in any way the motion or location of the phenolic proton in the chromophore, and we have found that it seems to move as fast as the other two. The motion of the three protons appears to be essentially concerted, rather synchronous in nature, and very fast. Because of the quality of the methods used to compute the potential energy surface, it seems likely that the potential energy landscape describing the motion of the phenolic proton has been reasonably accounted for. Consequenly, it seems that the cleavage of the O-H bond in the chromophore cannot be held responsible for the time scale of (many) picoseconds reported experimentally for the process of decay of A*. The question is, then, what process could be responsible for the time scale observed. To suggest an answer to this question, it is worthwhile to analyze further the data from our dynamical simulations. If the fluorescence measured at 460 nm originates in the still protonated chromophore, then the amount of protonated chromophore in our system at a given time can be computed by adding the populations of all CXY species (where X ) W or S, and Y ) S or G). Figure 12 displays the value of the population of CXY species, corresponding to A*, over time for an extended simulation lasting 15 ps. Two intriguing facts become apparent from analysis of Figure 12. First, in the perprotio species significant recurrences in the population of the CXY species disturb the otherwise monotonous decay between 2 and 6 ps. Second, for both isotopomeric variants of the proton wire, the population of the CXY species reaches an asymptotic value after only ∼8 ps, and this value is significantly different from zero (0.13 and 0.07 for the perprotio and perdeutero isotopomers, respectively): after ∼8 ps, in terms of distribution of products and reactants of the process, the operation of the wire is at an end and a kind of equilibrium is reached. Nevertheless, it would be interesting to try and fit our theoretical CXY populations in Figure 12 to some sort of multiexponential decay formula, in order to provide a set of

5510 J. Phys. Chem. B, Vol. 112, No. 17, 2008 decay times to compare with experiment. However, the behavior of the CXY population for the perprotio isotopomer does not correspond to a monotonic decay and thus does not lend itself to such a fitting process, barring some extensive conditioning of the data set, which would compromise the quality of the results. For this reason, we do not provide a set of decay times quantitatively derived from the time evolution of the CXY population. This does not prevent us from guessing the time scale of different dynamical processes because they can be seen from the dynamical results presented so far. As described in the preceding sections, the dynamical behavior of the proton wire seems to undergo at least two completely different regimes: a very impulsive process occurring during the first tens of femtoseconds in which about 50% of the CXY species is consumed, followed by a slower decay of the CXY species to a certain asymptotic value denoting equilibrium. The process of transferring from reactants to products has been estimated to reach 50% of the extent in around 20 fs in our simulations for the perprotio isotopomer (see Figure 6), and about 30 fs for the perdeutero isotopomer (see Figure 8), thus belonging to the very short-time dynamics regime, taking place in the first moments after photoexcitation. Proton motion during this stage is essentially concerted. In the impulsive phase of the dynamics the process is notably faster for the perprotio than it is for the perdeutero. A similar ultrafast impulsive component has been observed for other systems.48 In the actual experimental setup of Chattoraj et al., the excitation pulse had a duration of 70 fs.13 This could mean that the actual proton transfer could be occurring during these initial phases of the process. Besides, experimental resolution is nowadays established around 50100 fs, which means that all measurements available might have a blur in the initial time where the impulsive motion phase occurs. Later on, a slower decrease of the CXY population sets on, which eventually reaches the plateau at its asymptotic value. Right after photoexcitation, the system has a very large energy excess and it is likely that it evolves impulsively, mainly in the reaction channels that include the cleavage of the chromophore’s O-H bond (and thus being captured in Figures 6 and 8), but also along other secondary channels in which this bond is not cleaved (channels where only the second or third protons are transferred), as can be seen in Figure 10. Once the impulsive phase of the process is over, the intermediate products of these secondary reaction channels will evolve exoergically to the final product of the process, eventually breaking the chromophore’s O-H bond in a slower “leaking” process. This can be seen for the first 500 fs in Figure 10, where for instance species CSGs formally still species A*sdecreases very slowly its population while turning into WSG, which increases almost concomitantly. This overall finalizing process involving the evolution of all CXY protonation states (and including the bouncing of all protons from final products back to reactants) must start where the impulsive process ends, and from Figure 12, it must take about 6-8 ps to complete. Thus, this might justify a decay time in the range of 3-4 ps. It is tempting to identify the experimental decay time τ1 ) 3.6 ( 0.3 ps with this dynamical process. Finally, equilibrium is attained after ∼8 ps in our simulation, with 87% (perprotio) or 93% (perdeutero) WSG species. Overall, our model predicts a dynamical process, initiated by an ultrafast concerted motion of all protons, taking 8 ps to depopulate CXY species, and ending in an equilibrium. Although one can try to identify the first decay constant experimentally measured to the transformation of CWS into WSG, it is clear that two other, larger decay times were experimentally

Vendrell et al. measured that cannot be accounted for in this way. Besides, an equilibrium cannot be deduced from the dynamical data published by Chattoraj et al.13 However, as mentioned before, the model on which the simulations have been computed is closed; that is, it does not provide an irreversible exit to the proton shuttled by the wire after it has been accepted by Glu222. Including such an escape route would have made the computational model significantly more cumbersome at this stage. When the simulation time is increased, chances are increased that the wavepacket initially propagated from the reactants region (CWS) has already reached the products region (WSG) and is effectively being reflected back to reactants. Such behavior can be detected as “recurrences” in the dynamical plots: for instance, a sudden increase of the population of reactant species or concomitantly a decrease in the amount of product species. Take for instance the perprotio case in Figure 12. As can be seen, after the initial smooth decay lasting ∼2 ps, there is a sudden rise in the population of CXY, which can only be interpreted as part of the products reVerting to reactants. This behavior has not been observed experimentally. It would not happen to the same extent if the current model included an irreVersible exit to the product of the wire, such as the splitting of the wire, or the possibility to extend the shuttle further. Both possibilities have been discussed in the literature on GFP. The striking difference with the perdeutero proton wire, in which no large-scale recurrences like that occur, is illuminating. Again, the perdeutero species is less likely to bounce back to reactant’s region because its larger mass goes along with a larger difficulty in entering classically forbidden areas. Should such an exit channel exist and have an energy barrier associated to it, the perdeutero proton wire would be less efficient than its perprotio counterpart in reaching the end of the process. Thus, if the model were enlarged from A* h I/1 to A* h I/1 f I/2 it is likely that a normal and large KIE would be theoretically determined for the process of formation of I/2. In this context, the longer decay times measured experimentally (12.0 ( 1.0 ps and 120.0 ( 15.0 ps)13 would correspond to this irreversible step. We think that the results derived from our theoretical model indicate that, in order to explain the rise times and KIEs experimentally detected, an exit channel to form species I/2, possibly with an energy barrier, is needed. 4. Conclusions In this paper, we have presented the first quantum mechanical simulation of the operation of the proton wire responsible for fluorescence in the green fluorescent protein on a reliable highdimensional potential energy surface that also includes the degrees of freedom corresponding to the distance between donor and acceptor atoms. Both the perprotio and perdeutero versions of the proton wire have been studied. Careful analysis of the dynamical results shows that after photoexcitation the hydrogen bond between Cro and Wat25 undergoes a sudden shortening. This is the main effect of photoexcitation on the structure of the proton wire. Ensuing proton motion is very fast and within 20 fs the expected position of all protons is closer to the respective acceptor atoms. After this sudden burst of motion it slows down, but the proton coordinates continue to approach the acceptor atoms until the end of the simulation. Our study concludes that proton motion in the wire is essentially concerted and quite synchronous. Simulations carried out on specific isotopic substitutions of the wire show also that motion of all protons is tightly coupled.

Proton Wire in Green Fluorescent Protein The population of protonated chromophore has been computed over time. It is possible to discern at least two dynamical regimes operating on the wire. At very short times after photoexcitation, when the system has excess energy, the wire is under the effect of impulsive dynamics lasting a few tens of femtoseconds. Within this short period a large fraction of the systems develops an ionized chromophore, even though not necessarily the final product of the proton-shuttle. Later on, a slower decrease in the amount of protonated chromophore can be seen. Within the simulation time of 1.5 ps, a normal (larger than unity) kinetic isotope effect (KIE) is observed during the initial stages of the dynamics, which turns into an inVerse KIE in the assymptotic limit. This anomaly is discussed in terms of a specific feature of the model, this being that the model in use does not contemplate an irreversible exit of the product of the proton transfer of the wire. Likely processes would be the escape of the proton to the surface of the protein via an extended protonwire network or the effective splitting of the wire. Were such an irreversible step incorporated in the model, it is likely that the KIE measured would probably be normal instead. We think that this provides theoretical support to the existence of such irreversible processes in the operation of the proton wire. The current model is used to explore the long-time limit result of the dynamics. It is found that an equilibrium is attained within ∼8 ps of photoexcitation, with about 10% of protonated chromophore remaining. Overall, after the impulsive phase is over, the decrease of the population of the protonated chromophore occurs within a time frame that could correspond to the first decay constant measured experimentally (3.6 ( 0.3 ps).13 This equilibrium is the logical outcome in a closed model as the one proposed in this study. If the model included an irreversible step after the deprotonation of the chromophore, then the aforementioned equilibrium would not occur. In such a case, the larger decay times reported experimentally (12.0 ( 1.0 ps and 120.0 ( 15.0 ps)13 would correspond to this irreversible step. From the results presented here, it seems that the proton relay must break up after its operation, or otherwise involve some irreversible step. This would correspond to the process I/1 f I/2 (see Figure 3), which has been described in the literature as involving a conformational change of Glu222.7 Work along this line is currently underway in our laboratory. Acknowledgment. We are grateful for financial support from the “Ministerio de Educacio´n y Ciencia” and the “Fondo Europeo de Desarrollo Regional” (project CTQ2005-07115/ BQU) and from the “Generalitat de Catalunya” (project 2005SGR00400). O.V. acknowledges the European Commission for a Marie Curie individual fellowship. References and Notes (1) Zimmer, M. Chem. ReV. 2002, 102, 759-781. (2) Shimomura, O. FEBS Lett. 1979, 104, 220-222. (3) Cody, C. W.; Prasher, D. C.; Westler, W. M.; Prendergast, F. G.; Ward, W. W. Biochemistry 1993, 32, 1212-1218. (4) Heim, R.; Prasher, D. C.; Tsien, R. Y. Proc. Natl. Acad. Sci. U.S.A. 1994, 91, 12501-12504. (5) Prasher, D. C.; Eckenrode, V. K.; Ward, W. W.; Prendergast, F. G.; Cormier, M. J. Gene 1992, 111, 229-233. (6) Yang, F.; Moss, L. G.; Phillips, G. N. Nat. Biotechnol. 1996, 14, 1246-1251. (7) Brejc, K.; Sixma, T. K.; Kitts, P. A.; Kain, S. R.; Tsien, R. Y.; Ormo¨, M.; Remington, S. J. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 23062311. (8) Ward, W. W.; Cody, C. W.; Hart, R. C.; Cormier, M. J. Photochem. Photobiol. 1980, 31, 611-615.

J. Phys. Chem. B, Vol. 112, No. 17, 2008 5511 (9) Morise, H.; Shimomura, O.; Johnson, F. H.; Winant, J. Biochemistry 1974, 13, 2656-2662. (10) Kummer, A. D.; Kompa, C.; Lossau, H.; Pollinger-Dammer, F.; Michel-Beyerle, M. E.; Silva, C. M.; Bylina, E. J.; Coleman, W. J.; Yang, M. M.; Youvan, D. C. Chem. Phys. 1998, 237, 183-193. (11) Heim, R.; Cubitt, A. B.; Tsien, R. Y. Nature 1995, 373, 663-664. (12) Ormo¨, M.; Cubitt, A. B.; Kallio, K.; Gross, L. A.; Tsien, R. Y.; Remington, S. J. Science 1996, 273, 1392-1395. (13) Chattoraj, M.; King, B. A.; Bublitz, G. U.; Boxer, S. G. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 8362-8367. (14) Palm, G. J.; Zdanov, A.; Gaitanaris, G. A.; Stauber, R.; Pavlakis, G. N.; Wlodawer, A. Nat. Struct. Biol. 1997, 4, 361-365. (15) Lossau, H.; Kummer, A.; Heinecke, R.; Po¨llinger-Dammer, F.; Kompa, C.; Bieser, G.; Jonsson, T.; Silva, C. M.; Yang, M. M.; Youvan, D. C.; Michel-Beyerle, M. E. Chem. Phys. 1996, 213, 1-16. (16) Lill, M. A.; Helms, V. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 2778-2781. (17) van Thor, J. J.; Zanetti, G.; Ronayne, K. L.; Towrie, M. J. Phys. Chem. B 2005, 109, 16099-16108. (18) Stoner-Ma, D.; Jaye, A. A.; Matousek, P.; Towrie, M.; Meech, S. R.; Tonge, P. J. J. Am. Chem. Soc. 2005, 127, 2864-2865. (19) Stoner-Ma, D.; Melief, E. H.; Nappa, J.; Ronayne, K. L.; Tonge, P. J.; Meech, S. R. J. Phys. Chem. B 2006, 110, 22009-22018. (20) Agmon, N. Biophys. J. 2005, 88, 2452-2461. (21) Leiderman, P.; Huppert, D.; Agmon, N. Biophys. J. 2006, 90, 10091018. (22) Agmon, N. J. Phys. Chem. B 2007, 111, 7870-7878. (23) Weber, W.; Helms, V.; McCammon, J. A.; Langhoff, P. W. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 6177-6182. (24) Martin, M. E.; Negri, F.; Olivucci, M. J. Am. Chem. Soc. 2004, 126, 5452-5464. (25) Sinicropi, A.; Andruniow, T.; Devico, L.; Ferre, N.; Olivucci, M. Pure Appl. Chem. 2005, 77, 977-993. (26) Vendrell, O.; Gelabert, R.; Moreno, M.; Lluch, J. M. J. Am. Chem. Soc. 2006, 128, 3564-3574. (27) Vendrell, O.; Gelabert, R.; Moreno, M.; Lluch, J. M. Chem. Phys. Lett. 2004, 396, 202-207. (28) Wang, S. F.; Smith, S. C. Phys. Chem. Chem. Phys. 2007, 9, 452458. (29) Wang, S. F.; Smith, S. C. J. Phys. Chem. B 2006, 110, 50845093. (30) Sinicropi, A.; Andruniow, T.; Ferre, N.; Basosi, R.; Olivucci, M. J. Am. Chem. Soc. 2005, 127, 11534-11535. (31) Vendrell, O.; Meyer, H. D. J. Chem. Phys. 2005, 122, 104505. (32) Warshel, A. J. Phys. Chem. 1982, 86, 2218-2224. (33) Vendrell, O.; Moreno, M.; Lluch, J. M.; Hammes-Schiffer, S. J. Phys. Chem. B 2004, 108, 6616-6623. (34) Cembran, A.; Gao, J. Theor. Chem. Acc. 2007, 118, 211-218. (35) Vendrell, O.; Gelabert, R.; Moreno, M.; Lluch, J. M., submitted. (36) Meyer, H. D.; Manthe, U.; Cederbaum, L. S. Chem. Phys. Lett. 1990, 165, 73-78. (37) Manthe, U.; Meyer, H. D.; Cederbaum, L. S. J. Chem. Phys. 1992, 97, 3199-3213. (38) Beck, M. H.; Ja¨ckle, A.; Worth, G. A.; Meyer, H. D. Phys. Rep. 2000, 324, 1-105. (39) Meyer, H. D.; Worth, G. A. Theor. Chem. Acc. 2003, 109, 251267. (40) Worth, G. A.; Beck, M. H.; Ja¨ckle, A.; Meyer, H.-D. The MCTDH Package, Version 8.2, 2000. Meyer, H.-D.The MCTDH Package, Version 8.3, 2002. See http://www.pci.uni-heidelberg.de/tc/usr/mctdh/. (41) Petkovic´, M.; Ku¨hn, O. J. Phys. Chem. A 2003, 107, 8458-8466. (42) Coutinho-Neto, M. D.; Viel, A.; Manthe, U. J. Chem. Phys. 2004, 121, 9207-9210. (43) Ortiz-Sa´nchez, J. M.; Gelabert, R.; Moreno, M.; Lluch, J. M. J. Phys. Chem. A 2006, 110, 4649-4656. (44) Solntsev, K. M.; Huppert, D.; Agmon, N. J. Phys. Chem. A 1998, 102, 9599-9606. (45) Solntsev, K. M.; Huppert, D.; Agmon, N. J. Phys. Chem. A 1999, 103, 6984-6997. (46) Mohammed, O. F.; Pines, D.; Dreyer, J.; Pines, E.; Nibbering, E. T. J. Science 2005, 310, 83-86. (47) Siwick, B. J.; Bakker, H. J. J. Am. Chem. Soc. 2007, 129, 1341213420. (48) Jimenez, R.; Fleming, G. R.; Kumar, P. V.; Maroncelli, M. Nature 1994, 369, 471-473.