Exploring the Effects of Intramolecular Vibrational Energy

Sep 30, 2008 - Marc Nadal-Ferret , Ricard Gelabert , Miquel Moreno , and José M. Lluch. Journal of Chemical Theory and Computation 2013 9 (3), 1731-1...
0 downloads 0 Views 1MB Size
Download PDF
J. Phys. Chem. B 2008, 112, 13443–13452

13443

Exploring the Effects of Intramolecular Vibrational Energy Redistribution on the Operation of the Proton Wire in Green Fluorescent Protein 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: June 9, 2008; ReVised Manuscript ReceiVed: July 28, 2008

The operation of the proton wire in Green Fluorescent Protein has been simulated by quantum dynamics and considering the coupling to the protein environment by means of a bath of harmonic oscillators. The simulation consists of 36 explicit and fully quantum degrees of freedom: 6 degrees of freedom represent the configuration of the proton wire, which are coupled to 30 bath coordinates. Regimes of weak and strong coupling have been studied. It is found that presence of the bath induces a fast energy transfer from the proton wire to the bath, with characteristic times under 400 fs. This internal vibrational redistribution happens at the expense of the potential energy content of the proton wire, deformed through the interaction to the bath from its uncoupled state. Strong coupling induces a slowing-down of the operation of the wire because it hinders to some extent the approaching of donor and acceptor atoms to distances in which proton transfer can occur. Internal vibrational energy redistribution affects the dynamics, but from our simulations we conclude that it cannot be the only cause responsible for the experimentally reported fluorescence rise times. 1. Introduction Because of its important use as fluorescent biological marker, Green Fluorescent Protein (GFP) has attracted attention of researchers worldwide in the field of biochemistry and structural biology. GFP is a natural protein encountered in the North Pacific jellyfish Aequorea Victoria. In its natural form, GFP absorbs strongly in the blue region of the electromagnetic spectrum and produces instead bright green fluorescence with high quantum yield (Φ ∼ 0.72-0.85), which is at the core of its technological uses.1 Quite naturally, the photophysical properties of GFP are tightly interwoven with its structure. The crystallographic structure of wild-type GFP (wt-GFP) was determined over a decade ago.2 GFP occurs either as a monomeric or dimeric protein, even though apparently its photophysical properties do not depend on this fact. Structurally, GFP resembles a rigid barrel whose walls are made up of 11 β-sheets. This rigid barrel hosts an R -helix that holds the chromophore (Cro) roughly at the center of the cavity. The consensus nowadays is that the chromophore is stabilized through hydrogen bonding interactions to a captive water molecule (Wat25), which in its turn is bound to a serine (Ser205), which finally connects to a glutamate residue (Glu222).3 Certainly, other interactions exist within the cavity that stabilize the chromophore and keep it fixed in place, like those induced by residues Thr203 and Ser65 (the latter belonging actually to the chromophoric region). Figure 1 shows the most relevant interactions around the chromophore in GFP. A strong correlation exists between the structure of the protein around the chromophore and the photophysical properties of * To whom correspondence should be addressed. E-mail: ricard.gelabert@ uab.cat. † Universita ¨ t Heidelberg. ‡ Departament de Quı´mica, Universitat Auto ` noma de Barcelona. § Institut de Biotecnologia I de Biomedicina, Universitat Auto ` noma de Barcelona.

Figure 1. Main stabilizing interactions playing a role in the GFP photocycle. W are water molecules and dotted lines depict hydrogen bond interactions. Labels in brackets identify the location of amino acids before the autocatalytic cyclization that brings the chromophore into existence.

GFP. Most eloquently, the chromophore has been shown to produce no fluorescence when synthesized and studied in solution.4 wt-GFP absorbs at 397 nm (band A), this being due to absorption of the neutral form of the chromophore, has a smaller absorption at 477 nm (band B) whose intensity depends notably on pH, and fluoresces strongly at 508-510 nm. The identity of the fluorescent species has been inferred from the study of the Ser65Thr (S65T) mutant of GFP, which shows no absorption at 397 nm but shows instead a band B more intense and red-shifted.5 The fact that the S65T mutant was found to contain a deprotonated chromophore provided evidence to identify the fluorescent species with a ionized chromophore.6 Currently, the most widespread model for the photophysics of GFP is supported by time-resolved fluorescence experiments7

10.1021/jp805049c CCC: $40.75  2008 American Chemical Society Published on Web 09/30/2008

13444 J. Phys. Chem. B, Vol. 112, No. 42, 2008

Figure 2. Model used in ref 11 to study the dynamics of the proton wire.

and is summarized here. Species A (containing a neutral chromophore) is photoexcited at 397 nm to A*, which can fluoresce back to the ground electronic state or deactivate while in the photoactive state through the operation of a 3-proton relay which commences with the phenolic proton of the chromophore, and includes Wat25 and Ser205, to end up at Glu222.3,8 The intermediate species obtained in this way is usually referred to as I*. However, this species forms very quickly: Time resolved fluorescence measured at 508 nm (I* form) could be fitted to a biexponential rise function with τ1 ) 2.2 ( 0.2 ps and τ2 ) 8.1 ( 0.5 ps.7 This fact, combined with the arrangement of residues around the chromophore after arrival of the proton to Glu222, makes it likely that this species is formed first in a nonrelaxed form (I/1) and would evolve to a relaxed form (I/2), which differs from the former by 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.3 Green fluorescence after photoexcitation of A is thought to come from the I/2 form.9,10 Further and slower conformational changes are needed to turn I/2 into the final photoproduct B*, consisting of the rotation of residue Thr203 to stabilize the ionized chromophore. Theoretical studies provide strong interpretative tools to understand the physics that underpine complex processes, without which a reasonable understanding of the latter is very difficult and only of approximative value. Recently we have published a comprehensive, high-dimensional quantum dynamical study of the operation of the proton wire in the photoactive state,11 using state-of-the-art quantum dynamical simulation techniques on an accurate potential energy surface (PES) of the ground (S0) and first excited singlet (S1) electronic states, describing motion of the protons and donor and acceptor atoms for a modeled version of the chromophore and proton wire, depicted in Figure 2 and computed at the Complete Active Subspace Self-Consistent Field level of theory with secondorder perturbative corrections to the energy (CASPT2).12,13 The wave packet propagation scheme used mimicked the process of photoexcitation and ulterior time evolution for a total of 15 ps, restricted to the A* a I/1 process. In this study it was discovered that motion of the three protons is almost simultaneous and very fast. Two different dynamical regimes in the formation if I/1 were determined: within the first femtoseconds dynamics is impulsive and very fast; later on a slower pace of formation sets in, which within the model is identified with the asymptotic decay into equilibrium: being a ”closed” system the ultimate fate of the wire is I/1 (and not B*, nor even I/2). In these simulations, the perprotio proton wire had attained equilibrium before 8 ps had elapsed.

Vendrell et al. Why does a quantum dynamics simulation of the operation of the proton wire on a high quality PES turn out to be so much faster than experimentally measured? Definitely, the model used (Figure 2) does not include important residues and thus is missing hydrogen bonding interactions: Thr203, which stabilizes Wat25, is missing, as is Ser65, which stabilizes Glu222. Their absence can have an effect on the energetics of the wire operation. The cost of computing a sizable part of a highdimensional PES at a high level of theory, and later, do quantum dynamics, is a sufficient justification not to have enlarged the model even more. Besides, there are almost no alternative methods to enable inclusion of protein environment in costly multiconfigurational computations.14 None of these residues interacts with the protons that are being shuttled directly, and, in any case, they would have, as an effect, a likely stabilization of the products of the individual proton transfers. In this way, it is unlikely that their effect would be a substantial slow-down of the rate at which the process occurs. A different possibility can be conceived that might influence the rate at which the proton wire in GFP operates. The proton wire is supported by the protein backbone, to which it is attached in many places and, besides, by side-chain interactions such as those described before and known to stabilize the proton wire. Upon photoexcitation the chromophore is promoted to the photoactive electronic state and because of the topology of the potential energy surface as seen in the model,12,13 the system evolves exoergically by transferring the three protons. However, it is reasonable to assume that the proton wire is also dynamically coupled to the protein backbone, and thus, able to exchange energy with it. Part of the excess energy on the proton wire might be transferred to the different vibrational modes on the backbone and used to vibrationally excite those, depleting the energy content in the wire and thus slowing down the protontransfer process. This process is generally termed intramolecular vibrational energy redistribution (IVR) and is known to plague other photosystems, like for instance bipyridyl systems.15,18 IVR can be a very fast process for chromophores of the size of that in GFP,16,17 and thus it might compete successfully with the generation of I/1 species. In this work we explore the effects of IVR upon the dynamics of the proton wire in GFP. Because of the sheer size of the system, strong modeling is needed. Given that reasonably accurate dynamical simulations have been performed on a model of the proton wire which, in essence and by itself, is a model in vacuo, in this work the effect of dynamical coupling to other degrees of freedom in the system is going to be explored taking the aforementioned simulation as starting point.11 This dynamical coupling requires these additional degrees of freedom to have full dynamical entity, and thus they appear as additional degrees of freedom in the simulation. Besides, these additional degrees of freedom should represent the other nuclear degrees of freedom corresponding to molecular deformations of the backbone, essentially localized in bonds situated close to the proton wire. This would be very costly to parametrize. Instead, we will reproduce the qualitative effect of these degrees of freedom representing them as dynamical shifted harmonic oscillators (bath) coupled to the proton wire and will examine what changes occur to the dynamics of the proton transfer when the coupling regime between proton wire and this reservoir is varied from weak to strong. This should provide an indicative picture of what kind of effects a large number of coupled degrees of freedom is expected to have on the proton wire, and thus, examine the possibility that IVR sits behind the rate of transfer experimentally reported for the proton wire.

Operation of the Proton Wire

J. Phys. Chem. B, Vol. 112, No. 42, 2008 13445

Figure 3. Dynamical model for the uncoupled proton wire. Distances between donor and acceptor atoms are represented by Rj (j ) 1...3) depending on the hydrogen bond. Coordinates rj are defined as the (signed) distance from the jth proton to the point halfway-between its donor and acceptor atom.

2. Computational Details The purpose of this work is to explore the effect on the dynamics of the proton wire in a dynamical model of the GFP chromophoric site, of coupling to a certain large number of degrees of freedom, trying to reproduce at least qualitatively the dynamical effect of the protein environment. Because the physical process involves motion of light atoms, quantum dynamics has to be employed. To solve the time-dependent Schro¨dinger equation (TDSE), the multiconfigurational timedependent Hartree (MCTDH) method is used,19-21 using the Heidelberg MCTDH package.22 As such, this work is a continuation of our study of the dynamics of the uncoupled proton wire,11 and naturally, the dynamical model here is an extension of the model used in that work. We summarize it here. 2.1. Uncoupled Model. The quantum chemical model used to compute the potential energy surface of the S0 and S1 electronic states is represented in Figure 2. The model includes all atoms explicitly appearing in the GFP proton wire, but the side chains out of the wire itself are modeled to reduce the cost of multiconfigurational computations. Thus, Ser205 is turned into a methanol and Glu222 into an acetate. The hydrogen bonds in the model are kept linear and global Cs symmetry is enforced.12 The dynamical model adopts a linear configuration for the complete proton wire, restricting the motion of the atoms in the wire to take place along the line defined by the four oxygen atoms. The model requires the following coordinates to specify the state of the proton wire: R1, R2, and R3 are the distances between each donor-acceptor atom pair, starting at the chromophore and going all the way to Glu222, and r1, r2, and r3 are the signed distances between each proton and the point halfway-between its donor and acceptor atoms. Thus, the uncoupled model of the proton wire requires six independent degrees of freedom to be completely specified: see Figure 3. The potential energy for the uncoupled model was computed for a total of 338 structures covering all relevant areas of configurational space for these six parameters at the CASPT2 level of theory.12,13 To make the computation of the potential energy possible, an interpolating 8-state Empirical Valence Bond (EVB)23 potential energy function of the six geometrical s0 s1 variables was fitted for S0 (VPW ) and S1 (VPW ). Technical details are given elsewhere.13 To solve the TDSE an explicit form of the Hamiltonian operator is needed: σ ˆ PW(q) ) TˆPW(q) + VˆPW H (q) σ VˆPW

(1)

Here, PW refers to “proton wire”. is the potential energy surface for electronic state σ (σ ) S0, S1), TˆPW is the kinetic energy operator for the coordinate system shown in Figure 3

Figure 4. Dynamical model for the coupled proton wire. Besides the coordinates of the proton wire itself, a set of harmonic oscillators is coupled to each hydrogen bond, Q. Only coupling between the heavy atom motion (Rj) in the proton wire and these oscillators is considered. Each hydrogen bond is coupled to a set of harmonic oscillators.

(also given elsewhere),11,24 and q is the following six-coordinate vector:

q ) (r1, r2, r3, R1, R2, R3)

(2)

which completely specifies, within the dynamical model, the geometry of the proton wire. 2.2. Coupled Model. To qualitatively introduce the dynamical effect that the immediate surrounding of the proton wire has on the behavior of the latter, we make use of a bath of harmonic oscillators coupled to the proton wire’s degrees of freedom. This kind of model has been used before to represent dissipation in mixed quantum-classical studies of proton-transfer and in semiclassical approaches.25,26 The proton wire in GFP is supported by interactions between the protein side chains and the donor and acceptor atoms of the proton wire: Thr203 interacts with the oxygen atom of Wat25, as does Asn146, and Ser65 interacts with one of the oxygens of Glu222. In view of this we have decided to consider explicit coupling between the bath modes and only the heavy atom coordinates of the proton wire, Rj. Each hydrogen bond is coupled to a set of N harmonic oscillators (thus, a total of 3N oscillators are considered). This is represented schematically in Figure 4. Within the model, the geometry of the coupled system is known once q and Q are specified, where q is given by Equation 2 and Q is the 3N component vector

Q ) (Q11, Q12, ... , Q1N, Q21, Q22, ... , Q2N, Q31, Q32, ... , Q3N)

(3) where the first subindex indicates the hydrogen bond to which each oscillator is coupled and the second subindex identifies the oscillator. Given that the desired effect of the bath is to act as an energy reservoir/sink, the actual values of the bath coordinates are irrelevant to the study. Take for instance the ith oscillator coupled to the jth hydrogen bond, Qji. Let us assume that it is assigned a frequency νji. Consider the following transformation of the bath coordinate



˜ )Q Q ji ji

mjiji p2

where mji is the mass assigned to the bath mode and

(4)

13446 J. Phys. Chem. B, Vol. 112, No. 42, 2008

ji ) hνji

Vendrell et al.

(5)

˜ ji can be seen to be an adimensional has units of energy. Q quantity. The Hamiltonian of the coupled system can be written as follows:

˜) ˆσ)H ˆ σ(q, Q H σ ˜) + Vˆ ˜ ) TˆPW(q) + VˆPW (q) + TˆBath(Q Coupling(q, Q) (6) σ (here again, σ ) S0, S1). Here TˆPW and VˆPW are the same as in the uncoupled system. With the above definition of the bath coordinates, the kinetic energy of the bath has the form

3

N

∂2 ˆ (Q ˜) ) - 1 T  Bath 2 j)1 i)1 ji ∂Q ˜2

∑∑

(7)

ji

while the potential energy term is 3

N

∑∑

ˆ ˜) ) 1 ˜ - c (R - R0)]2 VCoupling(q, Q  [Q ji j j 2 j)1 i)1 ji ji

(8)

3

˜, t) ) ΨBath(Q

N

∏ ∏ ψji(Q˜ji, t) j)1 i)1

(11)

Such an approximation should remain good as long as the number of bath coordinates is large enough so that each single bath coordinate is weakly coupled to the system, also in the case that the overall effect is large. It has the further advantage that the presence of the bath does not add new configurations to the total wave function. To analyze the results of the dynamical propagation, the expectation values of some observables have been computed along the dynamics. If one is interested in obtaining the expectation value of an observable X, the following quantity must be computed

˜, t)|Xˆ|Ψ(q, Q ˜, t)〉 〈X〉 ) 〈Ψ(q, Q

(12)

where Xˆ is the quantum mechanical operator associated to ˜ , t) is taken as observable X, and the wave function Ψ(q, Q normalized. 3. Results and Discussion

where cji are the coupling constants that determine the coupling between bath modes and the main system and R0j are shift parameters, introduced to condition the equilibrium position of the heavy atoms in the main system to retain different values than those from the uncoupled model. In doing this, the mass of the bath modes has vanished from the formulas and the coupling constants can be expressed dimensionally as L-1. We note that even though the proton wire’s potential energy function is definitely different in the ground and photoactive electronic states, it should be a good approximation that the nature of the coupling to the protein environment (in other words, νji, cji, and R0j ) remains approximately constant. Thus, VˆCoupling is taken the same in the simulations in S0 and S1. 2.3. Simulation Protocol. The goal of this study is to explore differences in the dynamical behavior between coupled and uncoupled versions of the proton wire model as it evolves in the photoactive state. A good simulation of the Franck-Condon excitation is needed. A traditional approach in quantum 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 ˆ S0 (eq 6), until invariance in the (imaginary) timeHamiltonian H evolved state, in which case the ground vibrational state is obtained:

t f -i∞

˜) f Ψ (q, Q ˜) Ψ(q, Q 0

(9)

Photoexcitation is then simulated by placing the ground ˜ ), on S1, that is, solving the TDSE vibrational state, Ψ0(q, Q ˜ ) and Hamiltonian H ˆ S1 . with initial condition Ψ0(q, Q ˜ ). In this work We only need to provide a form for Ψ(q, Q we use a Hartree separation of the system and bath variables to define the total wave function, which is given by

˜) ) Ψ (q, t)Ψ (Q ˜, t) Ψ(q, Q PW Bath

(10)

ΨPW is a multiconfigurational expansion including enough configurations for a converged propagation and is the same used previously.11 ΨBath is taken to be a Hartree product on the bath coordinates:

The dynamical study of the reference uncoupled system has been already published.11 Before tackling the quantum dynamical simulations for the coupled system, the bath needs to be specified, which means selection of a set of frequencies, coupling constants and shift parameters. We are aware that the degree of modeling is important and this will prevent us from doing quantitative comparisons to experimental data. We think, however, that the comparison of the results for the uncoupled system and the system coupled to a dissipative bath can provide useful information that can help interpret the behavior of the actual system. We have tried to select these bath parameters in such a way that these comparisons are possible. A reasonable and affordable number of bath modes with frequencies in the range of those present in the surroundings of the proton wire have been selected. A frequency calculation of the proton-wire model has been useful to detect that skeletal vibrations in the proton wire (excluding those explicitly represented in the uncoupled model) cover a spectrum of values between 200 and 700 cm-1. Because of the size of the coupled system the number of bath modes cannot be too large or the simulations become cumbersome. A value N ) 10 has been chosen, yielding for all simulations done in this work a total of 36 quantum degrees of freedom. The frequencies selected are equally spaced with a spacing of 50 cm-1. We have taken 10 distinct frequencies, but these are the same for each hydrogen bond. A molecular dynamics simulation of the complete GFP in the ground-state within a sphere of water molecules produced OsO distances for the three hydrogen bonds between 2.6 and 2.9 Å.12 However, the vibrational ground-state for the quantum chemistry model (which does not include environment) shifted this range to 2.45-2.55 Å. The R0j values in eq 8 have been introduced to modify this equilibrium distance between donor and acceptor atoms of the three hydrogen bonds. In this work all R0j values have been kept fixed at 3.0 Å. As for the coupling constants cji, we have opted for selecting a set of values for these that range from weak to strong coupling. This has been motivated by the impossibility to model the protein environment in detail. Instead, our choice will allow us to see the effects on the dynamics of the proton wire when strong or weak coupling is present, and provide sufficient elements to do an educated guess as to what could be occurring inside GFP. Similar to what has been done for the frequencies and shift parameters, a unique value has been selected for cji for the full proton wire. The values

Operation of the Proton Wire

J. Phys. Chem. B, Vol. 112, No. 42, 2008 13447

S1 Figure 5. Energy analysis of the uncoupled proton wire for the first 100 fs (left) and 1 ps (right). In each graph, values of 〈TˆPW〉 (TPW), 〈VˆPW 〉 (VPW), ˆ S1 〉 (EPW) are shown. Energy origin is the value of each energy type right after photoexcitation. and 〈H PW

corresponding to the simulations presented here are 3 Å-1 for the weak coupling and 6 Å-1 for the strong coupling. For the baths just described a quantum dynamics simulation has been done as explained in Section 2.3. After the groundstate in S0 is found, it is promoted to S1 and the ensuing evolution of the system with time is followed for 1 ps. The wave function is analyzed and the results discussed. 3.1. Energetic Aspects. In order to analyze the changes in the dynamical behavior of the proton wire due to the coupling to a dissipative bath, it is interesting to assess the distribution of energy among different degrees of freedom of the wire. We expand the squared term in Equation 8 and rewrite it as follows: 3

N

∑∑

3

N

∑∑

ˆ ˜) ) 1 ˜ 2+ 1 VCoupling(q, Q jiQ  [c 2(R - Rj0)2 ji 2 j)1 i)1 2 j)1 i)1 ji ji j ˜ (R - R0)] 2cjiQ ji j j ˜) + Vˆ (q, Q ˜) )ˆ VBath(Q int

(13)

where VˆBath represents the harmonic oscillator’s energy and Vˆint (interaction potential energy) collects the effects of the shifting terms and the nonseparable terms in the potential. Now, the Hamiltonian for S1 can be rewritten as

ˆ S1(q, Q ˜) ) T ˆ (q) + VˆS1 (q) + T ˆ (Q ˜) + Vˆ (Q ˜) + H PW PW Bath Bath ˜) Vˆint(q, Q S1 ˆ (Q ˜) + Vˆ (q, Q ˜) ˆ PW )H (q) + H Bath int

(14)

S1 ˆ PW ˆ Bath, and Vˆint have been The expectation values of H , H computed over time to quantify the distribution of energy between bath and proton wire. Figure 5 shows the time evolution of energetic values for the uncoupled proton wire, and Figures 6 and 7 show the same information for the weakly and strongly coupled baths, respectively. Note that the origin of energies in the graphs is the value of each kind of energy at t ) 0, that is, their value just after the Franck-Condon transition. The uncoupled proton wire (Figure 5), quite obviously has a constant expectation value for the total energy. It is interesting ˆ s1 〉 instead. After photoexcitation the to analyze 〈TˆPW〉 and 〈V PW proton wire is in a nonstationary state and a very fast decrease of the potential energy is observed, along with a concomitant rise in kinetic energy. In about 20 fs this trend stops and the

values of kinetic and potential energy remain approximately constant with small oscillations around a central value; if at ˆ s1 〉 drops from about -1.8 to all, over the following 500 fs 〈V PW -1 -2.2 kcal mol , and the complementary change is observed for the kinetic energy. The ultrafast decrease in potential energy of the first 20 fs, and the fast decrease afterward up to times of about 500 fs, roughly corresponds to the two dynamical regimes already described,11 perceived now from an energetic point of view. We highlight that after the first 6-7 fs both energy types ˆ s1 〉 passes through a minimum undergo a recurrence (that is, 〈V PW and returns to the original value at the Franck-Condon transition). This initial “heartbeat” must be due to the initial motion of the S0-relaxed wave packet and its first bounce with a potential energy obstacle before leaving for the products area, and is a hint of the existence of a possible, if small, potential energy barrier somewhere between reactants and products of the proton wire. Let us turn now to the weakly coupled model (Figure 6). A glance at the energy content of the proton wire reveals that there is a decrease, basically at the potential energy’s expense. Again, there is an ultrafast energy transfer from the proton wire to the bath in the initial stages of the simulation (∼20 fs) followed by a sustained, but slower, decrease for the following 500 fs. From the point of view of the proton wire this behavior is the same as for the uncoupled proton wire, except for the fact that there is an actual, if small, energy flow out of the proton wire to the bath. Things look different for the strongly coupled model (Figure 7): in the first 100 fs a sudden decrease of total energy in the proton wire of ∼10 kcal mol-1 occurs, that must be transferred to the bath. Thus, the bath is here acting as a very efficient energy sink. It is not that in the strongly coupled simulation the proton wire visits much more stable areas of the potential energy surface than in the other simulations; rather, the cji2(Rj - R0j )2 terms in Vˆint in eq 13 force the Rj values to be large at the beginning of the simulation, and this causes the proton wire to start, at t ) 0, with large values for the potential energy. As soon as the simulation starts in S1, after the impulsive motion of R1 this excess potential energy on the wire is converted into kinetic energy which in its turn is transferred to the bath. This can be observed mainly in the fact that this draining effect happens again at the expense of the potential energy of the proton wire, while the kinetic energy actually ranges between 0 and 2 kcal mol-1 from its value at the Franck-Condon point. The bath conditions here the behavior of the proton wire after these first 100 fs: a series of periodic oscillations in the total and

13448 J. Phys. Chem. B, Vol. 112, No. 42, 2008

Vendrell et al.

S1 ˆ S1 〉 (EPW) for the first 100 fs (left) and 1 Figure 6. Energy analysis of the weakly coupled proton wire. Top: values of 〈TˆPW〉 (TPW), 〈VˆPW 〉 (VPW), and 〈H PW ˆ Bath〉 (left) and 〈Vˆint〉 (right) per bath mode. Energy origin is the value of each energy type right after photoexcitation. ps (right). Bottom: values of 〈H

ˆ s1 〉 (VPW), and 〈H ˆ s1 〉 (EPW) for the first 100 fs (left) and Figure 7. Energy analysis of the strongly coupled proton wire. Top: values of 〈TˆPW〉 (TPW), 〈V PW PW ˆ Bath〉 (left) and 〈Vˆint〉 (right) per bath mode. Energy origin is the value of each energy type right after photoexcitation. 1 ps (right). Bottom: values of 〈H

potential energies of the proton wire match completely the ˆ Bath〉. Thus, strong coupling is oscillations over time of 〈H causing here the proton wire to evolve following very tightly

the bath and not the other way around. The dissipation of energy from the proton wire is very fast: in 100 fs it is gone from the wire and does not return again to it.

Operation of the Proton Wire Both bath simulations show comparable behavior as far as ˆ Bath〉, is concerned. Although the energy content of the bath, 〈H the intensity of the oscillations (the energy values themselves) ˆ Bath〉 depend on the strength of the coupling, in both cases 〈H increases steeply at the beginning of the simulation, attaining a maximum value in about 100 fs. After photoexcitation, the proton wire has an energy excess and can go to areas of its configurational space with lower potential energy through operation of the wire. However, as seen for the uncoupled system’s simulations,11 for this to occur the system needs to bring closer together the donor and acceptor atoms of each hydrogen bond in the proton wire. This has to happen also here, 1 has not been changed, but the oxygen atoms that donate as VˆSPW or accept each proton intervene also in the interaction potential energy of the bath Vˆint. Reducing the OsO distances (that is, Rj values) occurs with a certain difficulty now: the R0j terms, shift terms, are large (3 Å) and after relaxation in S0 the wave packet has large amplitude in areas of configurational space with Rj larger than in the uncoupled system. The only possibility for the proton wire to operate with ease is to bring the oxygen atoms ˜ ji move to large closer together, which brings down 〈Vint〉 if Q negative values. However, this occurs with energetic penalty ˆ Bath. If there were no other possibility most likely because of H the proton-transfer processes would stop altogether. The model contains an stabilizing term (Vˆint) which allows this process to ˆ Bath〉 occur. It is nonetheles interesting to note that 〈Vˆint〉 + 〈H add up roughly to the same quantity throughout each simulation, that is, the deformation of the bath is paid for by the interaction potential, plus a certain amount of energy flowing out of the proton wire. IVR happens all the time during these simulations, as can be S1 ˆ PW seen by the continuously dropping curves for 〈H 〉. A look at this curve for the weakly coupled case (Figure 6) reveals that the asymptotic value is attained before 1 ps. If a characteristic time could be assigned to the IVR for the weakly coupled bath, a value in the range of τIVR ∼ 300-400 fs would be reasonable. In the strongly coupled model, the behavior is markedly different, and if it could be defined, such a decay time would be much shorter. The bath model used in this work (which incidentally, is the Caldeira-Leggett model)27 has a frictional (dissipative) effect and does not alter the equilibrium position of the proton wire from its uncoupled state substantially. IVR is not negligible, but most of the action in the proton wire occurs in the impulsive phase of the dynamics and is over long before 200 fs have elapsed, so IVR can only very weakly influence the outcome at this stage. After 200 fs have elapsed, the operation of the proton wire is almost at an end, the effects of IVR can be felt, the energy is gone to the bath and its operation is substantially more irreversible than for the uncoupled proton wire. Before moving on toward other aspects of dynamics with dissipation to a bath, there is an unexpected result worth analyzing. The value in all simulations in the asymptotic limit for 〈TˆPW〉 is approximately the same. The actual values at 1 ps are 8.6 (uncoupled proton wire), 8.4 (weakly coupled proton wire), and 10.0 kcal mol-1 (strongly coupled proton wire). In our simulations, the model is closed (considering the proton wire and the bath together). After the proton wire has finished its operation, as discussed elsewhere,11 the system has transferred all three protons and from a chemical point of view no other changes can happen, so it reaches a kind of equilibrium. In this sense, this kinetic energy at long times could be interpreted as coming from the vibrational zero-point energy (ZPE) of the proton wire in its product state in S1. The different range of

J. Phys. Chem. B, Vol. 112, No. 42, 2008 13449 values determined, though close together, is due to the different effects that coupling has on the wire in each case. 3.2. Mechanistic Aspects. We turn now to a different analysis of the simulations presented so far, and focus on the time-evolution of the expectation values of the geometrical parameters of the wire, 〈rj〉 and 〈Rj〉, for each hydrogen bond in it. We note here that, because of the particular definition of the 〈rj〉 coordinates, the jth hydrogen bond has adopted its product form as soon as rj goes from negative to positive values. The expectation values for the geometrical parameters of the proton wire are depicted together in Figure 8. The values of 〈rj〉 for the simulations of the uncoupled and weakly coupled proton wires are very similar: ultrafast creation of the products of the proton wire operation (all 〈rj〉 > 0) within 20 fs of photoexcitation. This is in agreement with the results of the previous section. Proton-transfer is too fast to be affected by IVR. However, IVR is acting throughout for the weakly coupled proton wire: looking at the asymptotic values of 〈rj〉 reveals that the weakly coupled proton wire displays values closer to zero (i.e., protons less extensively transferred) than the uncoupled model. The strongly coupled model presents a much more distinct behavior. The values of 〈Rj〉 at the beginning of the simulation are substantially larger than for the uncoupled and weakly coupled simulations, and resemble more the values of the OsO distances found in the molecular dynamics simulations.12 As we have mentioned before, this may be due to the effect of the cji2(Rj - R0j )2 terms rather than to the bath itself. The motion of the protons toward the products region is delayed: all 〈rj〉 take almost 100 fs to cross to products’ area and interestingly, in these initial stages of the process 〈rj〉 present recurrences (a series of alternate maxima and minima). It is rather straightforward to provide an explanation to this observation. The OsO distances start off very large now and thus proton-transfer is very difficult. The donor-acceptor distances get close together and, in the meantime, the protons move rather fast, but cannot overcome the energetic hindrance of transferring from donor to acceptor when they are too separated. When OsO distances attain values similar to those in the uncoupled or weakly coupled models (∼2.45 Å) proton-transfer occurs with ease. The recurrences observed are mainly (but not only) in 〈r1〉, which makes us think that this is the proton that encounters a larger obstacle when the process starts with the oxygens far from their preferred positions. To better understand these dynamical aspects of the operation of the wire, Figure 9 shows the trajectory that each proton follows on a plane defined by 〈rj〉 and 〈Rj〉. It can be seen that for the uncoupled proton wire the three transferrable protons move straight to the products’ area (〈rj〉 > 0) with slight contraction of the heavy atom’s coordinate R1. Once each proton crosses over this surface, a certain number of back and forth motions from products to reactants and back again are observed (“recrossing”), even though at 200 fs (the final point of the dynamics) the three protons are localized in the products’ region. For the strongly coupled bath, the picture is markedly different. 〈Rj〉 values start very large and thus proton transfer is severely hindered. To occur, 〈Rj〉 must decrease from its initial values (2.65-2.75 Å) to values under 2.55 Å. While this happens, protons bounce many times back and forth trying to reach the products’ region: this is easily seen for H1 in the zigzag line it follows from the upper-left corner of the graph, which undergoes no less than four complete vibrational periods. By the end of this motion, the heavy atom coordinate is already short enough to enable proton-transfer.

13450 J. Phys. Chem. B, Vol. 112, No. 42, 2008

Vendrell et al.

Figure 8. Time evolution of 〈rj〉 (left) and 〈Rj〉 (right) coordinates for the uncoupled proton wire (top), weakly coupled proton wire (center) and strongly coupled proton wire (bottom).

Figure 9. Trajectory of each proton over the initial 200 fs of the dynamics for the uncoupled (left) and strongly coupled (right) models of the proton wire, represented as a parametric function of time, on a coordinate system representing 〈rj〉 and 〈Rj〉. Scales on both axes have been kept equal to display the different range of motion of all protons in both cases. At t ) 0 the three protons start their motion at the left part of each graph and tend to move toward the right when time increases.

3.3. Effects of IVR on the Dynamics. Inclusion of a bath that interacts with the proton wire does alter the dynamics of

the latter. Even for a weakly coupled bath there is a certain energy draining that brings energy from the proton wire to the

Operation of the Proton Wire bath. Kinetic energy remains roughly the same for the proton wire in all simulations studied, so the energy transfer is done at the expense of the potential energy of the wire. This potential energy is mostly deformation energy, coming from the lengthened distances between donor and acceptor atoms (Rj) induced by the interaction potential, Vˆint. With this proviso, the energy flow is fast, with characteristic times of the order of 200-400 fs for the weakly coupled bath and much faster than that for the strongly coupled model. Coupling has also effects on the expectation values of the structural parameters of the wire. On the other hand, in the strongly coupled proton wire, coupling alters the asymptotic values of the geometrical parameters also, but a different effect is to be noticed at short times: while for the uncoupled proton wire all protons reach the product region during the impulsive phase within 20 fs,11 for the strongest coupling studied this is effectively delayed until 100 fs have elapsed (see Figure 8). This slowing-down occurs as a consequence of the fact that the donor and acceptor atoms start off much farther apart and proton-transfer can only occur quickly when these are within shorter distances. The time increase to obtain the product of operation of the proton wire is related to the time it takes for each donor-acceptor pair to go from 2.7 Å to about 2.5 Å. However this 5-fold increase in the transfer time happens only for a highly coupled bath, able to drain within 100 fs no less than 10 kcal mol-1 of energy from the proton wire, which is roughly equivalent to the total kinetic energy of the wire itself. Such strong coupling causes the wire to follow the time evolution of the bath and not the other way around. Thus, we think that even though IVR has an effect on the operation of the wire, and that this effect is to delay the operation of the wire, it is not intense enough to delay the formation of the product species of the wire to times in the range of picoseconds, as measured experimentally, for any conceivable coupling existing around the proton wire in GFP. If IVR does not stand behind the experimentally measured fluorescence rise and decay times, the question remains as to what is the likely cause of their actual values. In our previous work,11 we found that long-time simulations of the uncoupled proton wire were too fast when compared with the reported experimental decay times. It was also found that for times of about 2 ps, perprotio variants of the wire experienced an increase in the population of protonated (i.e., A* form) of the wire, which was interpreted in terms of part of the product formed reverting to reactants. This was possible because the model of the proton wire is closed, and thus, the wave packet, upon reaching the products’ area had no other option but to bounce toward the reactants again. This behavior could only be avoided if an irreversible exit would exist, such as the splitting of the wire between Ser205 and Glu222 occurred after operation because of the anti-syn inversion of the newly formed OsH bond and attachment of the Glu222 to Ser65. Besides, if such process involved a certain energy barrier, the kinetic isotope effect of about 5 detected experimentally for the fluorescence rise times could be accounted for. We are currently improving the current model to include such possibility in our laboratory. 4. Conclusions In this work the effect of dissipation of energy from the proton wire in GFP to the protein environment has been studied, by means of a simplified model of the proton wire in its photoactive state11-13 and simulating the coupling to the protein environment through coupling to a bath of harmonic oscillators. This system has been studied quantum mechanically integrating the timedependent Schro¨dinger equation for a global system consisting

J. Phys. Chem. B, Vol. 112, No. 42, 2008 13451 of 36 explicit degrees of freedom. Two different regimes of coupling have been studied: weak and strong coupling, and the results have been compared to the uncoupled simulations. The energetic and mechanistic aspects of the simulations have been explored. It has been found that even for the weakly coupled model there is an effective energy transfer from the proton wire to the bath, which is taken mainly from the potential energy of the proton wire. This excess potential energy arises from the deformation of the proton wire due to its coupling to the bath. This energy transfer is seen to be very fast, occurring even for the weakest coupling with characteristic times under 400 fs. The strongly coupled model presents the proton wire following closely the evolution of the bath, and energy dissipation is even faster in this case. Scrutiny of the expectation values of the structural parameters of the wire reveals that for the strongest coupling studied, there is an actual delay for the transfer times of the protons to the products’ area of configurational space: while in the uncoupled model this happened within 20 fs of photoexcitation, for the strongly coupled model this occurs within 100 fs. The reasons behind this substantial slowing down lie in the fact that for the strongly coupled model the starting point of the simulation presents donor-acceptor distances notably larger than those found in the uncoupled model. Thus, for proton-transfer to occur, the donor and acceptor atom pairs must get closer, this being the time it takes until effective proton-transfer can occur. The coupling to a bath and the corresponding IVR processes can appreciably change the behavior of the wire with respect to the uncoupled case, mostly increasing the formation time of the fluorescent species. However, the reported effect cannot be the only cause behind the larger experimental fluorescence rise times, considered to represent the rates at which the formation of the fluorescent species occur. Work is in progress in our laboratory to shed light in these matters. Acknowledgment. The authors are grateful for financial support from the “Ministerio de Educacio´n y Ciencia” and the “Fondo Europeo de Desarrollo Regional” (Project CTQ200507115/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) Yang, F.; Moss, L. G.; Phillips, G. N. Nat. Biotechnol. 1996, 14, 1246–1251. (3) 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, 2306– 2311. (4) Niwa, H.; Inouye, S.; Hirano, T.; Matsuno, T.; Kojima, S.; Kubota, M.; Ohashi, M.; Tsuji, F. I. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 13617– 13622. (5) Heim, R.; Cubitt, A. B.; Tsien, R. Y. Nature 1995, 373, 663–664. (6) Ormo¨, M.; Cubitt, A. B.; Kallio, K.; Gross, L. A.; Tsien, R. Y.; Remington, S. J. Science 1996, 273, 1392–1395. (7) Chattoraj, M.; King, B. A.; Bublitz, G. U.; Boxer, S. G. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 8362–8367. (8) Palm, G. J.; Zdanov, A.; Gaitanaris, G. A.; Stauber, R.; Pavlakis, G. N.; Wlodawer, A. Nat. Struct. Biol. 1997, 4, 361–365. (9) Lill, M. A.; Helms, V. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 2778– 2781. (10) Weber, W.; Helms, V.; McCammon, J. A.; Langhoff, P. W. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 6177–6182. (11) Vendrell, O.; Gelabert, R.; Moreno, M.; Lluch, J. M. J. Phys. Chem. B 2008, 112, 5500–5511. (12) Vendrell, O.; Gelabert, R.; Moreno, M.; Lluch, J. M. J. Am. Chem. Soc. 2006, 128, 3564–3574. (13) Vendrell, O.; Gelabert, R.; Moreno, M.; Lluch, J. M. J. Chem. Theory Comput. 2008, 4, 1138–1150.

13452 J. Phys. Chem. B, Vol. 112, No. 42, 2008 (14) Sinicropi, A.; Andruniow, T.; Ferre, N.; Basosi, R.; Olivucci, M. J. Am. Chem. Soc. 2005, 127, 11534–11535. (15) Marks, D.; Prosposito, P.; Zhang, H.; Glasbeek, M. Chem. Phys. Lett. 1998, 289, 535–540. (16) De Belder, G.; Schweitzer, G.; Jordens, S.; Lor, M.; Mitra, S.; Hofkens, J.; De Feyter, S.; Van der Auweraer, M.; Herrmann, A.; Weil, T.; Mu¨llen, K.; De Schryver, F. C. Chem. Phys. Chem. 2001, 2, 49–55. (17) Chou, P. T.; Chen, Y. C.; Yu, W. S.; Chou, Y. H.; Wei, C. Y.; Cheng, Y. M. J. Phys. Chem. A 2001, 105, 1731–1740. (18) Gelabert, R.; Moreno, M.; Lluch, J. M. Chem. Phys. Chem. 2004, 5, 1372–1378. (19) Meyer, H.-D.; Manthe, U.; Cederbaum, L. S. Chem. Phys. Lett. 1990, 165, 73–78. (20) Beck, M. H.; Ja¨ckle, A.; Worth, G. A.; Meyer, H.-D. Phys. Rep. 2000, 324, 1–105.

Vendrell et al. (21) Meyer, H.-D.; Worth, G. A. Theor. Chem. Acc. 2003, 109, 251– 267. (22) Worth, G. A.; Beck, M. H.; Ja¨ckle, A.; Meyer, H.-D.; Heitz, M.C.; Wefing, S.; Manthe, U.; Sukiasyan, S.; Raab, A.; Ehara, M.; Cattarius, C.; Gatti, F.; Otto, F.; Nest, M.; Markmann, A.; Brill, M. R.; Vendrell, O. The MCTDH Package, Versions 8.3, 8.4; See http://www.pci.uni-heidelberg.de/tc/usr/mctdh, 2000-2007. (23) Warshel, A. J. Phys. Chem. 1982, 86, 2218–2224. (24) Vendrell, O.; Meyer, H.-D. J. Chem. Phys. 2005, 122, 104505. (25) Kim, S. Y.; Hammes-Schiffer, S. J. Chem. Phys. 2006, 124, 244102. (26) Wang, H. B.; Thoss, M.; Sorge, K. L.; Gelabert, R.; Gime´nez, X.; Miller, W. H. J. Chem. Phys. 2001, 114, 2562–2571. (27) Caldeira, A.; Leggett, A. J. Phys. ReV. Lett. 1981, 46, 211–214.

JP805049C