Transient low-barrier hydrogen bond in the photoactive state of green fluorescent protein.

PCCP View Article Online

Published on 24 April 2015. Downloaded by Freie Universitaet Berlin on 09/05/2015 16:50:07.

PAPER

Cite this: DOI: 10.1039/c5cp01067e

View Journal

Transient low-barrier hydrogen bond in the photoactive state of green fluorescent protein ab Marc Nadal-Ferret,a Ricard Gelabert,*a Miquel Morenoa and Jose ´ M. Lluch

In this paper, we have analyzed the feasibility of spontaneous proton transfer in GFP at the Franck–Condon region directly after photoexcitation. Computation of a sizeable portion of the potential energy surface at the Franck–Condon region of the A structure shows the process of proton transfer to be unfavorable by 3 kcal mol1 in S1 if no further structural relaxation is permitted. The ground vibrational state is found to lie above the potential energy barrier of the proton transfer in both S0 and S1. Expectation values of the geometry reveal that the proton shared between the chromophore and W22, and the proton shared between Ser205 and Glu222 are very close to the center of the respective hydrogen bonds, giving support to the claim that the first transient intermediate detected after photoexcitation (I0*) has characteristics similar to those of a low-barrier hydrogen bond [Di Donato et al., Phys. Chem. Chem. Phys., 2012, 13, 16295]. A quantum dynamical calculation of the evolution in the excited state shows an even larger probability of finding those two protons close to the center compared to in the ground state, but no formation of the proton-transferred product is observed. A QM/MM photoactive state geometry optimization, initiated using a configuration obtained by taking the A minimum and moving the protons to the product side, yields a minimum energy structure with the protons transferred and in which the His148 residue is substantially closer to the now anionic chromophore. These results indicate that: (1) proton transfer is not possible if structural relaxation of the surroundings of the chromophore is prevented; (2) protons H1 and H3 especially Received 20th February 2015, Accepted 24th April 2015

are found very close to the point halfway between the donor and acceptor after photoexcitation when the

DOI: 10.1039/c5cp01067e

zero-point energy is considered; (3) a geometrical parameter exists (the His148–Cro distance) under which the structure with the protons transferred is not a minimum, and that, if included, should lead to the fluor-

www.rsc.org/pccp

the triple proton transfer reaction can explain the dual emission reported for the I0* intermediate of wtGFP.

escing I* structure. The existence of an oscillating stationary state between the reactants and products of

1. Introduction Since its discovery, green fluorescent protein (GFP) has been the target of enduring scientific interest, fuelled by its application as a fluorescent marker in biomedicine.1–4 The component behind the photoactivity of GFP is its chromophore 4-( p-hydroxybenzylidene)-imidazolidin-5-one (HBDI), formed autocatalytically in a post-translational reaction involving residues Ser65, Tyr66 and Gly67. The absorption spectrum of GFP presents two maxima, a very prominent one at 395 nm (band A) and a smaller one at 475 nm (band B). The relative intensities of both peaks depend on the pH and a single clear isosbestic point has been reported.5 The pH dependence of the ratio of intensities of the A and B peaks was taken as a first indication that they are linked to different

a

Departament de Quı´mica, Universitat Auto`noma de Barcelona, 08193 Bellaterra, Barcelona, Spain. E-mail: [email protected]; Fax: +34 935812477; Tel: +34 935811669 b Institut de Biotecnologia i de Biomedicina, Universitat Auto`noma de Barcelona, 08193 Bellaterra, Barcelona, Spain

This journal is © the Owner Societies 2015

protonation states of the chromophore. Excitation of the A band results in fast emission from a different species (band I) at 510 nm with high intensity (quantum yield B 0.8). In contrast, the chromophore in solution shows no fluorescence. The dynamics of the photocycle of GFP have been studied via fluorescence upconversion by Chattoraj et al.,6 who found that the increase in the fluorescence of I* at 508 nm takes place very quickly after excitation of the A band at 398 nm. The increase in the fluorescence of I* can be fitted to a biexponential function with time constants of 2.2 ps and 8.1 ps. These values increase to 8.2 ps and 46.0 ps for the protein equilibrated in D2O, even though values as large as 69 ps have been reported for the second time constant in this case.7 The large value of the resulting kinetic isotope effect (KIE) is indicative of the involvement of a proton-transfer process in the rate determining step of the process in the excited state. Crystallographic structures of GFP have started appearing and they will bring molecular resolution to the discussion on the nature of the involved species.8–10 Brejc et al. obtained the crystallographic structures of wild type GFP (wtGFP) and

Phys. Chem. Chem. Phys.

View Article Online


Paper

mutant S65T.10 The latter variant displays a much altered absorption spectrum with no A band and a much larger B band. The crystal structures provided strong support that the proton transfer mechanism is the process responsible for the fluorescence of GFP. The proton is relayed from the chromophore at Tyr66 to a captive water molecule (W22), then to Ser205 and finally to Glu222. A photocycle was then proposed, which, with some variations, still continues to be accepted today. Species A consists of a neutral chromophore in a relaxed protein environment, species B consists of an ionized chromophore in the corresponding relaxed environment, and I contains an ionized chromophore surrounded by an unrelaxed A-type environment. The transformation A* - I* is very fast, while the transformation I* - B* is slow. Experimental evidence supporting this proton relay was obtained in sophisticated time-resolved transient infrared experiments where the protonation of Glu222 was observed.11–13 The biexponential dependence of the fluorescence of the I* species on time has been investigated as a way to access detailed knowledge of the dynamics of the photocycle. Some time ago, Fang et al., using femtosecond stimulated Raman spectroscopy to monitor the decay of certain A* modes, determined the existence of an essentially KIE-less component at the onset of the process (0.6–0.7 ps), suggesting that H(D) motion is not dominant at the beginning of the reaction coordinate in the excited state.14 The interpretation was in terms of a necessary reorientation of the chromophore, as it is not in the optimal arrangement to shuttle its proton to W22 straight after the excitation. After this impulsive phase, the proton transfer would take place. In fact a KIE of 4.3 was measured, which was attributed to a real step with proton transfer characteristics. More recently, Di Donato et al. performed a detailed study on the spectral evolution of the region between 1450 and 1700 cm1 for deuterated GFP, tracking the vibrational modes of neutral and ionized chromophores over time.15 Analysis of the evolutionassociated difference spectra (EADS) revealed that, at least for a fraction of the population, Glu222 receives the proton (deuteron) before the chromophore is ionized. The analysis in terms of the species-associated difference spectra (SADS) provided support for a sequential mechanism of decay of A* into I* via an intermediate (I0*), in contrast with previous work supporting parallel decay.12,16 This intermediate was identified as the excited state counterpart to the ground state intermediate in the process of the formation of A from I.7 Intermediate I0* was shown to have an emission spectrum with features common to both A* and I* (dual emission). This intermediate was suggested to have characteristics of a very short hydrogen bond (SHB) or even a low-barrier hydrogen bond (LBHB)17,18 in view of the above, and they proved that a LBHB signature was present when a broad absorption over 2000 cm1 was observed for this intermediate. This claim was supported by the fact that in the crystallographic structure of the A form, O–O distances along the proton-wire are anomalously short, in the range of 2.5–2.7 Å.10 Other cases connected to the GFP family have been reported where dual emission intermediates have been found.19 GFP itself has been reported to show transient dual emission between 3–10 ps after


PCCP

photoexcitation.20 Besides, the S65T–H148D mutant of GFP, where the proton-shuttle has been rewired and links the chromophore directly to Asp148, has been found to harbor a very short contact between the donor and acceptor oxygen atoms (about 2.4 Å)21 and some of its properties strongly suggest the existence of a LBHB.22,23 The LBHB concept is tightly connected to terms like the vibrational energy level of the proton(s) being transferred and its(their) expected position along the line linking the donor and acceptor atoms of the hydrogen bond(s). It is desirable to compare the tentative characteristics of this elusive intermediate, I0*, with the corresponding predictions drawn by means of precise and balanced theoretical studies, where these magnitudes are accessible. In this regard, several factors combine to make a meaningful theoretical study of the characteristics of I0* difficult to attain. The structure of this intermediate has been proposed in general terms but no definite structure has been published.15 Moreover, I0* describes an ensemble of structures of the protein in the excited electronic state, and it is not simple to compute reliable excitation energies of biological-sized molecules. Last but not least, to prove that I0* contains a LBHB or SHB, the zero-point energy (ZPE) has to be computed and verified to lie at or above the adiabatic electronic potential energy barrier, separating the reactants and products of the proton transfer in the photoactive state. This computation for an anharmonic potential energy profile is not a trivial issue even for the ground electronic state, taking into account the high dimensionality of the process. In a nutshell, quantum dynamical calculations of the vibrational eigenstates are needed for the excited electronic potential energy surface. Some theoretical work has been published on different aspects of GFP. A few concern themselves with dynamical aspects of the proton-relay process, but this intermediate has not been reported in any. Recently, Grigorenko et al. published a detailed study of the ground electronic state of wtGFP using high-level QM/MM methodology in which they located single structures corresponding to the A, I and B species.24 This study managed to describe very satisfactorily the energy differences between these species, and also provided excitation energies from these structures in good agreement with experimentally registered absorption peaks. These authors tried to locate an excited state structure akin to I0* but their attempts led them invariably to a form with a neutral chromophore. The failure to locate a structure with characteristics matching those ascribed to I0*, however, does not necessarily mean that this intermediate does not exist. If this intermediate has real LBHB character, it is expected that the average position of the protons affected is distant from that dictated by the geometry of the local stationary point found. This is the natural consequence of the anharmonicity of the potential energy surface and the relative values of the ZPE and adiabatic potential energy barrier of the proton transfer(s), which cause the vibrational wave function of the proton to peak closer to the energy barrier summit. In a recent paper,25 we discussed, in these terms, the local geometry of the LBHB reported to exist in the chromophore binding socket of photoactive yellow protein, as determined by


View Article Online


PCCP

neutron diffraction.26 Thus, if an intermediate like I0* exists and is suspected to have LBHB character, it might prove more elusive to find than a regular intermediate as its characteristics depend on the nuclear wave function and it is not possible to take for granted that a stationary point must exist in the potential energy surface in the corresponding area of configurational space. Quantum dynamical studies are then necessary. A decade ago we published a series of articles dealing with the dynamics of the proton-wire.27–30 These studies targeted high dimensional potential energy landscape exploration of ground and photoactive electronic states and quantum dynamical simulation of the translocation of the three protons in the wire. However, they suffered from some limitations. The need to consider a near dark excited state,31 combined with the unavailability of accurate density functionals able to correctly describe excited states with charge-transfer characteristics ¨dinger made the choice of methodology to solve the electronic Schro equation fall on multiconfigurational methods (the nature of the excitation in the chromophore involves partial charge transfer from the phenolic to the imidazole moieties), and it is difficult to consistently scan the entire span of configurational space of a chemical reaction with a single modest set of active orbitals. The most compromising approximation was likely the fact that the protein environment was missing, because no affordable QM/MM methodologies were available to cope with large portions of the potential energy surface in excited states. Within these limitations, a set of conclusions were obtained from the static studies:27,29 in topographical terms, the process of translocation of the 3 protons was exoergic. Full transfer of the three protons via sequences of individual transfers indicated that Ser205 - Glu222 was the proton-transfer to lead the way.27 A 6-dimensional quantum dynamical calculation, involving explicit consideration of protons, donor atoms and acceptor atoms, showed a very fast, essentially concerted, motion of the three protons that finalized within a few hundreds of femtoseconds,28 and vibrational

Paper

coupling to foreign modes could not account for the experimentally observed timescale of the process.30 No other quantum dynamical studies of the excited state have been published to our knowledge. Advances in algorithmics over the last decade have furnished theoretically oriented chemists with density functionals that can deal correctly with charge-transfer excited states in a cost efficient way. Moreover, substantial advancement has taken place in the development of QM/MM codes that enable the use of varied levels of theory in complex systems in a balanced and controlled way, making it possible to finally determine potential energy landscapes of processes involving excited states in proteins. Due to everything mentioned above, we plan to re-visit the process of proton translocation involved in the A* - I* step, computing a very large portion of the potential energy surface of the ground and photoactive states (large enough to encompass I0* as described by Di Donato et al.15), and perform a quantum-dynamical simulation on it to detect signatures of a transient LBHB.

2. Computational methods 2.1

Starting structure of the system

The starting point was PDB entry 1EMA.8 The PROPKA server was used to assign the protonation state of titratable residues at pH 8.0.32–35 Special attention was paid to the chromophore, where the S65T mutation present in the PDB structure was reversed. The system was centered and a sphere of water molecules of a radius equal to 40 Å was used to solvate the system. Fig. 1 shows a representation of the solvated system, which comprised 27 394 explicit atoms. 2.2

QM/MM single point calculations

Calculations of the electronic potential energy were carried out within a QM/MM framework, where the most chemically

Fig. 1 Left panel: model used in the QM/MM calculations including a 40 Å radius water droplet. Right panel: detail of the atoms belonging to the QM part of the system.



View Article Online


Paper

significant portion of the system was treated with quantum mechanics methods and the rest of the residues and solvent were described according to a MM force field. The effect of the MM part on the QM one was introduced using electronic embedding (EE), which lets the QM part be polarized owing to the presence of the cloud of point charges on the MM part. The QM part in this study comprised the chromophore (Cro), the side chains of residues Ser205, Glu222, and W22 (these being the supporters of the proton-wire) as well as the side chains of nearby residues Arg96, Thr203 and His148, plus another water molecule. Fig. 1 shows a representation of the QM part of the system isolated from the environment, which included 111 (including link atoms) QM atoms. This is essentially the same selection as that taken by Grigorenko et al. in their QM/MM study.24 The rest of the atoms of the protein and the water molecules were included in the MM part. Assembly of the QM and MM parts was done using the link atom approach. QM/MM calculations were carried out using the ChemShell software package.36 In all cases, the CHARMM2237,38 force field was used to compute the energy of the MM part. In the ground state calculations, the QM energy was computed using density functional theory (DFT),39 with functional B3LYP40,41 and basis set 6-31G+(d,p), using the Gaussian09 software package.42 For the excited state calculations, time-dependent density functional theory (TDDFT) was used, with the Coulombattenuated functional CAM-B3LYP to compute excitation energies.43 Coulomb-attenuated functionals are convenient as they avoid underestimation of excitation energies in excited states with charge-transfer characteristics. Current density functionals yield excitation energies about 0.5 eV too large for chromophores of the same size and type as the one in GFP,44,45 so this should be taken as a measure of the error bar in the absolute position of the excited state. As ChemShell does not interface with Gaussian09 in excited state calculations, we used a house version of some ChemShell scripts to achieve this. Both in the ground and excited state calculations, final energies included the polarization of the QM part due to the environment, and were computed with ChemShell36 with the charge-shift scheme. To find a candidate for the A structure, a ground-state optimization was done using DL-FIND,46 another component of the ChemShell package. This optimization started off from the structure set up as detailed before. All residues with at least one atom within 15 Å of the chromophore were defined as active system (i.e. atoms allowed to move during optimization). An analogous procedure was followed to optimize the I* structure in the excited state, this time with our own ChemShell script and starting from the A ground state structure, but with the protons within bonding distance of the respective acceptor atoms.

2.3

Quantum dynamics calculations

Dynamical simulations should describe correctly the behavior of protons, which makes it necessary to take into account their


PCCP

¨dinger quantum nature. This means that the time-dependent Schro equation must be solved, @ ^ i h Cðq; tÞ ¼ HCðq; tÞ @t

(1)

where q is a vector consisting of a convenient set of dynamical ˆ variables that represent the configuration of its particles, and H is the Hamiltonian (total energy) operator. The solution of this equation renders the time-dependent wave function C, where all information on the state of the system is enclosed. The ˆ depends on the choice of variables on which expression of H the dynamics is to be studied. Our quantum dynamics workhorse is the Heidelberg Multi-Configurational Time-Dependent Hartree software package (version 8.4),47 which conditions the choice of variables and the form of the Hamiltonian. 2.3.1 Dynamical model. To make the calculation feasible, the number of degrees of freedom needs to be kept small, mostly because of the effort required to compute a relevant portion of the potential energy surface. Moreover, in making the choice of this small set of variables, attention needs to be paid to the well-known coupling between the motion of the proton and the distance between the donor and acceptor atoms.48 The dynamical variables used here are the same as those used in ref. 28, and are depicted in Fig. 2. The model consists of 6 variables: Ri (i = 1–3) and ri (i = 1–3). R1, R2 and R3 correspond to the three O O distances in each of the hydrogen bonds (Cro–Wat22, Wat22–Ser205 and Ser205–Glu222). r1, r2 and r3 correspond to the three distances between each proton and the midpoint of the respective hydrogen bonds. Thus, negative values of ri indicate a situation where Hi is bound to the donor, and positive values indicate that it is bound to the acceptor. In this model it is assumed that the protons move exclusively along the straight line linking the respective donor and acceptor atoms. For this dynamical model, the kinetic energy operator is the one published in ref. 28. To prepare the potential energy

Fig. 2 Set of dynamical coordinates used to study proton and heavy atom motion. Coordinates R1, R2 and R3 correspond to the distances between the donor and acceptor atoms of the first, second and third hydrogen bonds, respectively. The position of the hydrogen atom in each bond is given by ri (i = 1–3) as the distance of the proton to the center of the respective hydrogen bond.


View Article Online


PCCP

Paper

operator, we proceeded as follows: starting from the minimum energy structure resulting from the QM/MM optimization in the ground state, a large set of 1426 distinct individual structures were generated keeping all the atoms fixed (this includes the QM part and the environment) but modifying the positions of the oxygen and hydrogen atoms in the proton-wire to correspond to unique arrangements (i.e. values of Ri and ri). Each of these structures was computed for S0 and S1 using the QM/MM approach described above. The resulting computed geometries and energies were fitted into an 8-state Empirical Valence Bond (EVB)49 to obtain a convenient analytical functional form to facilitate interpolation. The fit was done with a least-squares algorithm, minimizing the residual sum vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uN uP u wi V EVB ðqi ; pÞ V ref ðqi Þ 2 t i (2) sðpÞ ¼ N where VEVB(qi;p) is the energy given by the EVB interpolator using parameter set p on the geometrical parameters of reference point qi, Vref(qi) is the reference energy for the same point, and N is the number of reference points. wi is a weighting factor given by wi = min{1,exp[(vi vc)/a]}

(3)

where vi is the DFT/TDDFT potential energy at point qi relative to the lowest in that electronic state, vc is a cutoff value of the potential controlling the domain where the weight is active, and a controls the rate at which the weight approaches 0. Values of vc = 15 kcal mol1 and a = 10 kcal mol1 were used (as in our previous work).29 This procedure was done twice, once per electronic state. In this approach, each diagonal element of the EVB matrix describes one different protonation state of the proton-wire. Each of them is an explicit function of the dynamical variables {ri,Ri}i=1–3 and depends (parametrically) on a set of 55 parameters (a total of 110 parameters). The explicit formulas for the EVB matrix elements are the same as those found in ref. 29. The resulting EVB functions were turned into a MCTDH-usable form using the POTFIT program.47 2.3.2 Quantum dynamics simulation protocols. The Heidelberg Multi-Configurational Time-Dependent Hartree software package (version 8.4) was used to integrate the time¨dinger equation efficiently.47,50–55 Quantum dependent Schro dynamics simulations were used for two purposes: finding stationary (vibrational) states of the proton-wire system and determining the time evolution of the system after photoexcitation. To find the stationary states, propagation in imaginary time (i.e. relaxation) is done until invariance in the (imaginary) timeevolved state. In this way, the ground vibrational state is found. To converge excited vibrational states, the improved relaxation algorithm was used in an analogous way.53 To determine time evolution after photoexcitation, we mimicked the photoexcitation process as follows: first the ground vibrational state of the ground electronic state S0 was obtained through a relaxation calculation. Photoexcitation was then simulated by promoting this ground vibrational state


(the most populated) to the photoactive state S1, where it is not ¨dinger equation was stationary. Next, the time-dependent Schro integrated in real time with the Hamiltonian operator of the ˆ+V ˆ =T ˆ S1. Time evolution of the expectation photoactive state: H values of the observables or operators were computed by first implementing the corresponding operator (xˆ) in MCTDH form and then computing hx(t)i = hC(t)|xˆ|C(t)i

(4)

where |C(t)i denotes the quantum state of the system at time t. Except for the relaxation simulations, energy and unitarity conservation were checked at all times.

3. Results and discussion 3.1

QM/MM optimized structure of the A state

A minimum energy structure corresponding to the A state was found, as detailed in Section 2. A detailed view of the chromophore and its environment in this minimum energy structure is shown in Fig. 3. The structure we found corresponds to a conformer with a neutral chromophore that agrees reasonably well with the structure found by Grigorenko et al.,24 as well as with crystallographic data for the structure of the proton-wire. We found that our A state structure shows an interaction with His148 (the distance between the phenolic oxygen of the chromophore and atom Nt of His148 is 3.01 Å). Residue Thr203 has also been described as being involved in the stabilization of the chromophore. In our structure, Thr203 is placed such that an interaction takes place between Wat22 and the chromophore

Fig. 3 Chromophore and neighboring residues at the potential energy minimum found, corresponding to the A species. Yellow dotted lines denote the proton-wire. White dotted lines identify other close contacts. Numbers indicate the atomic distances in Å (yellow numbers are distances between atoms in the proton-wire, and white numbers are other distances).


View Article Online


Paper

PCCP

on one side, and the alcohol group of Thr203 on the other. Grigorenko et al. also found that the interaction with Thr203 was weaker than with His148, and that the orientation of Thr203 was not critical to stabilize the I intermediate, finding Thr203 normally facing outwards, as we do. The vertical excitation energy of this structure was computed to be 81.8 kcal mol1, which corresponds to a wavelength of approximately 349 nm. This value is off by 46 nm from the experimental A band maximum, which, in this spectral region, corresponds to about 0.4 eV. This value is within the reported accuracy of the TDDFT methodology for chromophores of a similar size.44,45 We conclude that the structure we found is a reasonable representation of the A state of GFP. 3.2

QM/MM potential energy surface

Our goal is to carry out a dynamical analysis of the different species involved in the proton-transfer in GFP, as well as to analyze the A* - I* transformation. The direct involvement of hydrogen atoms makes it necessary to account for their quantum nature. We opted to recover the quantum dynamical approach we used in our previous study of this system without an environment in ref. 28. An in-depth description of the procedure is given elsewhere.28,29 Fig. 2 shows graphically the definitions of the six variables that are used to represent the protonation state of the protonwire. Following the procedure described in Section 2.3.1, a large number of single-point energies for different {Ri}i=1–3 and {ri}i=1–3 6-tuples were computed at DFT (S0) and TDDFT (S1) levels of Table 1 a

Name DA bA d0A DB bB d0B DB+ bB+ d0B+ DC bC d0C DC+ bC+ d0C+ DD bD d0D e1B s1B e2C s2C e3D s3D eA1 sA1 eB2 sB2 eC3 a

theory and considering the environment. The results were fitted to an EVB functional form (one for each electronic state), as described in Section 2.3.1. The parameters of the fit and the residual fit errors are given in Table 1. The fitted surfaces can be considered as a comprehensive exploration of the A (A*) region of the full potential energy surface, as well as of the A " I (A* " I*) transformation. To discuss the energetics and geometries, we use the same naming scheme as in ref. 29, where the protonation state of the wire is identified by a three-digit binary code. In this code, each position in turn refers to the state of the first (Cro W22), second (W22 Ser205) and third (Ser205 Glu222) hydrogen bond, with ‘‘0’’ meaning that the proton is ‘‘not transferred’’ or bound to the donor atom, and ‘‘1’’ meaning that the proton is ‘‘transferred’’ or bound to the acceptor atom. Thus, the A state (with all the protons bound to their donor atoms) is denoted by ‘‘000’’ and the final state of the transfer by ‘‘111’’. It is instructive to analyze here, in some detail, the topography of the resulting potential energy surfaces. We located the local potential energy minima in the ground and photoactive electronic states, and the resulting values are shown in Table 2. In the electronic ground state, two different minima were found, one corresponding to the absolute reactant (000) and another to the final product of the transfer (111). We found the process to be endoergic by 5.9 kcal mol1. This contrasts with our previous findings in the gas phase and using CASPT2 methodology, where only the 000 structure was located.27,29 In the photoactive state a pair of minima were also located, and

Parameters of the potential energy surface for ground and photoactive states and fit residuals

Value (S0)

Value (S1)

141.24 5.56 0.96 213.66 1.18 0.75 4451.22 0.015 0.54 420.51 1.06 0.88 271.89 0.32 1.39 311.27 1.23 0.96 0.0035 1.84 29.20 0.87 25.24 0.83 0.57 1.70 1.13 1.10 58.26

Parameter names as in ref. 29.


188.64 4.04 0.97 188.24 1.11 0.66 2277.62 0.021 0.28 613.47 0.71 0.86 1305.88 0.046 1.42 412.71 1.07 0.95 4.70 1.05 19.71 0.87 35.08 0.82 0.060 2.05 0.53 1.18 56.02 b

Namea

Value (S0)

Value (S1)

Units

mol1

sC3 eAB sAB eBC sBC eCD sCD H12 H13 H14 H25 H26 H35 H37 H46 H47 H58 H68 H78 D1 D2 D3 D4 D5 D6 D7 D8

0.84 22.40 2.03 0.045 3.47 22.03 2.37 87.13 141.91 106.97 162.03 103.14 122.39 83.18 94.36 72.71 138.59 78.58 116.13 63.45 142.17 84.86 20.01 274.97 84.46 78.86 110.30

0.85 3.14 1.70 0.0017 4.55 16.12 2.43 82.42 173.85 129.44 278.70 111.63 122.74 99.45 93.83 81.68 208.86 60.36 117.09 110.67 201.11 168.34 73.57 486.36 75.94 138.46 158.21

Å kcal Å kcal Å kcal Å kcal kcal kcal kcal kcal kcal kcal kcal kcal kcal kcal kcal kcal kcal kcal kcal kcal kcal kcal kcal

mol1

s (eqn (2))b

0.75

0.86

Units kcal Å1 Å kcal Å1 Å kcal Å1 Å kcal Å1 Å kcal Å1 Å kcal Å1 Å kcal Å kcal Å kcal Å kcal Å kcal Å kcal

mol

1

mol1

mol

1

mol1 mol1 mol1 mol1 mol1 mol1 mol1

mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1 mol1

kcal mol1

Residual fit computed with vc = 15 kcal mol1, a = 10 kcal mol1.


View Article Online

PCCP Table 2 PESes

Paper Geometries and energies of the stationary points in the EVB

Stationary Eb (kcal mol1) r1 (Å) pointa


S0, S0, S0, S1, S1, S1,

000 111 TS 000 111 TS

0.0 5.9 7.2 0.0 3.0 5.6

0.128 +0.162 +0.111 0.0659 +0.182 +0.121

r2 (Å)

r3 (Å)

R1 (Å) R2 (Å) R3(Å)

0.275 +0.168 0.0135 0.302 +0.200 0.0544

0.203 +0.165 +0.0168 0.205 +0.190 0.00949

2.309 2.385 2.367 2.193 2.415 2.377

2.539 2.428 2.403 2.587 2.464 2.430

2.475 2.488 2.462 2.488 2.490 2.461

a Explanation of the naming scheme is given in the text. b Energy origin is taken as the absolute minimum in each electronic state.

again the 000 structure is more stable than the 111 structure by 3.0 kcal mol1. Our previous gas phase results also showed two minima, but the 111 structure was more stable by 6.8 kcal mol1. Grigorenko et al.24 found that the I structure was 1 kcal mol1 above the A structure. A comparison between their and our results needs to be done with care, as our 111 structure was not optimized (relaxed) outside the space of the six dynamical variables selected for the quantum dynamics study, while their result is the result of a more or less extensive minimization. We sought and located the transition states in the 6-dimensional variable space by minimization of the gradient norm. The potential energy barriers in S0 and S1 are 7.2 and 5.6 kcal mol1 respectively. Thus, the portrait of the potential energy landscape is markedly different compared to the one found with our previous methodology. The reduction of the endoergicity and barrier height in S1 is compatible with the increased acidity of the photoactive state. However, we found that these results agree well with our recent findings on GFP mutant S65T/H148D,56 where the pKa decrease of the chromophore is found to be much smaller than the estimated pKa of the isolated chromophore in solution,57 as DE changes only by about 3 kcal mol1 when going from S0 to S1. One way or another, simple translocation of the three protons without further structural changes in the chromophore and/or environment still seems to be unfavorable.

The reaction paths transit through a single transition structure in each state, so the proton translocation processes are concerted in topographical terms. Nevertheless, Fig. 4 shows that the process in S0 and S1 is now quite asynchronous: for negative r3 values (i.e. H3 not yet transferred), the iso-r3 slices have a pronounced L-shape where the first leg of the journey is mostly motion of H1 (the potential energy valley has gentler slopes in the direction of the r1 increase) and only when r1 Z 0 does the valley undergo a turn and aim towards increasing r2. This can already be seen in Table 2, as for the transition states for both surfaces, it is found that r1 4 r2, r3. As r3 accedes to its positive region, the L-shape shifts and a second minimum can be seen, corresponding to the complete transfer of the three protons until eventually even the reactant’s minimum vanishes. All this, from a static point of view, can be summarized as follows: the extent to which H1 and H2 are able to move along their respective coordinates is very limited until H3 reaches its product region. In other words, H3 opens the way to the complete proton-transfer. Of course, dynamical effects are altogether missing and can change this static image noticeably. 3.3

Quantum dynamics

The static picture provided by a topographical analysis can be deceiving as no consideration of the kinetic energy is included. The most accurate, though costly, description of the microscopic world we have is quantum mechanics (QM). In QM, the value of a geometrical parameter that will be determined in a measurement of a quantum system corresponds to the expected value of the operator associated with the geometric parameter, xˆ, via eqn (4). Insofar as the system lies in the ground state, the ground state is localized over a potential energy minimum, and this minimum can be adequately described as harmonic. It is reasonable to expect that the result of eqn (4) will match the value of the parameter at the topographical minimum. This is commonplace and the reason behind the standard procedure of identifying coordinates of a minimum in standard computational QM/MM studies of chemical systems with results of a

Fig. 4 Constant r3 potential energy surface slices of the ground (left) and photoactive (right) electronic states computed at geometries where the heavy atom coordinates are held at the values of the transition state (see Table 2). The five slices correspond (from bottom to top) to values of r3 of 0.2, 0.1, TS, +0.1 and +0.2 Å, where TS means the value of r3 at the transition state (+0.0168 Å in S0 and 0.00949 Å in S1). Color codes of the energy contours are common to both graphics. The black dots denote the positions of the topographical transition structure and correspond to the coordinates shown in Table 2 for S0(TS) and S1(TS).



View Article Online

Paper


Table 3

PCCP Expectation values of the dynamical coordinates in certain vibrational states

Isotopologue

Statea

Eb (kcal mol1)

hr1i (Å)

hr2i (Å)

hr3i (Å)

hR1i (Å)

hR2i(Å)

hR3i (Å)

HHH

cS0 0 cS0 1 cS1 1 cS2 1 cS3 1 cS4 1

13.27

0.0879

0.214

0.134

2.339

2.518

2.480

13.30

0.0450

0.224

0.128

2.304

2.550

2.485

14.09

0.0474

0.232

0.129

2.307

2.565

2.486

14.53

+0.00238

0.149

0.0755

2.338

2.521

2.483

14.87

0.0518

0.243

0.133

2.308

2.581

2.487

17.12

+0.0558

0.0552

0.000635

2.395

2.498

2.484

cS0 0

10.53

0.103

0.238

0.157

2.332

2.531

2.481

cS0 1

10.49

0.0606

0.260

0.159

2.272

2.572

2.488

DDD a

cYX

is the Xth vibrational state in electronic state Y.

b

Energy origin is taken as the absolute minimum in each electronic state.

geometry measurement. Nevertheless, when these conditions do not hold, substantial deviations from the values at the minimum should be expected. One such case in a biologically relevant system is represented by the chromophore in photoactive yellow protein (PYP), where electronic structure methods relying on the minimization of the potential energy alone (without dynamics) fail to account for the observed geometry of the two hydrogen bonds around the chromophore,58 but a quantum dynamical treatment including the effect of ZPE could succeed.25 E We determined the vibrational ground states in S0 cS0 0 E via a relaxation propagation of the S0 and S1 and S1 cS0 1 potential energy surfaces, as described in Section 2, and derived the expected values of the dynamical variables. These results are shown in Table 3. Let us focus, first of all, on the vibrational ground states E E S0 S c0 and c0 1 . The corresponding energy levels (that is, the lowest lying possible stationary states in each electronic adiabat) are almost equal to 13.3 kcal mol1 above the most stable minimum for each surface, while in the case of the perdeutero species, these figures drop to approximately 10.5 kcal mol1 for each electronic state. All these values are large, but they correspond to the zero point energy (ZPE) levels of a system with three bound hydrogen atoms. In any case, the most relevant issue is that they are all above the corresponding adiabatic potential energy barrier, and are therefore candidates from an energetic point of view for being described as LBHBs. Indeed, when computing the expectation values of the positions of the protons, we noticed that they were non-negligibly shifted from their potential energy minima coordinates and almost halfway to the values of the geometry of the TS. When the perdeutero system is considered the effect, even though visible, it is less noticeable as now the ZPE level is lower. Fig. 5 shows the 1-dimensional probability densities for the proton positions for the ground vibrational states in S0 and S1, ð 3 Y Y drj dRj rðri Þ ¼ jhc j cij2 jai

(5)

j¼1

So, although the expectation values of the geometry indicate a position for H1 that is markedly close to the center of the bond (hr1i = 0), the spread of the density makes it evident that


Fig. 5 One-particle probability density functions for the ground vibrational state for each proton of the perprotio isotopologue in S0 (bottom) and S1 (top).

this proton can be found over a large region of configurational space, including the products. For S0, the HHH (DDD) isotopologue in the ground vibrational state has a probability of 33% (29%) of being detected bound to the acceptor atom (hr1i 4 0). These figures drop to 20% (9.7%) for H3. The proton that is most clearly located ‘‘classically’’ bound to its donor atom E is H2, where the figures are 4.1% (0.7%). Thus, in cS0 0 , Cro H1 W22 satisfies the main criteria to be defined as a LBHB. To a lesser extent, the same can be said about the Ser205 H3 Glu222 hydrogen bond. What about S1? In static terms, the situation is even more evident as it corresponds to a potential energy surface closer to degeneracy and with a small barrier. So, ‘‘if’’ the system in S1 was allowed to relax within the S1 adiabatic state (in the sense of vibrational cooling) without spontaneously falling to S0, the ‘‘stationary state’’ reached would definitely show the characteristics shown in Table 3, namely an even stronger LBHB character for H1 and H3. The fractions of density over the product region of protons H1, H2 and H3 are 52% (44%), 4.7% (0.45%) and 22%


View Article Online

PCCP

Paper

(9.5%), respectively, for the perprotio (perdeutero) isotopologue. However, this is not what happens in the first picoseconds after the system is excited to S1. Following the Franck–Condon principle, excitation of the system does not land on a single vibrational state of S1 but rather on a superposition of states S X D S1 S0 E S1 E c 1 ðt ¼ 0Þ ¼ (6) ci c0 ci


i

where the bracket in eqn (6) is the Franck–Condon factor between the S0 ground vibrational state and each of the vibrational states in S1. |cS1i is not a stationary state of the S1 Hamiltonian. Time evolution is provided by S X D S S E iE t=h S E c 1 ðtÞ ¼ ci 1 c0 0 e i ci 1 : (7) i

Table 3 displays some properties of the lowest-lying vibrational states of S1 for the perprotio isotopologue – all of them with similar energies and mostly showing that H1 is expected in each of them to be close to the center of the bond. We expect, therefore, that, although not in a stationary state (whereby the LBHB definition is not applicable), eqn (7) will mix several states that, taken individually as if it was stationary, would describe H1 as satisfying the condition. The same applies to H3. To describe quantitatively the behavior of the system after photoexcitation in our model, we integrated the time¨dinger equation as described in Section 2, dependent Schro E taking cS1 ðt ¼ 0Þ ¼ cS0 as the initial condition. The resulting 0

time evolution of the dynamical coordinates is shown in Fig. 6. In essence, after an initial adjustment at the onset of the simulation (t o 20 fs) where H1 shifts towards larger r1 values, the system enters a phase of periodic oscillations around a stable average for the position of the proton. This is reminiscent of a Rabi oscillation59 entered when the system is no longer in a stationary state and the now time-dependent state ‘‘resonates’’ between a set of the actual stationary states of the Hamiltonian. In our case, hr1i oscillates around 0.03 Å and hr3i around 0.1. Heavy atom modes are activated at t = 0 with a compression of the Cro W22 distance and a stretch of the W22 Ser205 distance. Distance Ser205 Glu222 remains

approximately constant. Summing up: W22 shifts closer to the chromophore as a result of the photoexcitation to pick up H1. In summary, although we think that the term LBHB should not be used strictu senso because the system is not in a stationary state in S1, the general idea that the system has a large probability amplitude of being found in areas of configurational space where H1 is approximately at the center of its hydrogen bond is true, as is the fact that this character is more noticeable in S1 than in S0. We think that this confirms the hypothesis made by Di Donato et al. as to the nature of the I0* intermediate.15 Actually, Di Donato et al. observed protonation of Glu222 in about 3 ps (10 ps for D-GFP) to form I0* and then, in 10 ps, ionization of the chromophore. I0* should have a neutral chromophore with a LBHB-like interaction between Cro and W22. As has been mentioned before, a dual emission is detected experimentally, mostly corresponding to I* with a residual A* component.15 In our simulation, H1 reaches a state at the beginning of the simulation (t Z 20 fs) where it remains at almost hr1i = 0, and H3 follows, to remain at hr3i = 0.1. In the final configuration of the proton-wire, this means that the distance between H3 and Ser205 is approximately 1.14 Å and between H3 and Glu222 is approximately 1.34 Å. It looks like H3 does not form a bond to Glu222 in our simulation because a stable minimum is not found there. This state would explain the existence of a dual emission, but one with a dominating A* feature instead of I*. We think that this is rooted in the limitations of the model, as it has been built by freezing all atoms except the oxygen and hydrogen atoms directly involved in the proton-wire in their positions in the A structure. It might seem unlikely that a lot of heavy atom motion can take place to significantly change the energetics of the proton-transfer surface in the 3–10 ps of the initial stages where changes have been measured. However, Fang et al.14 determined that an impulsive phase lasting about 0.6 ps directly after photoexcitation was KIE-independent, which was interpreted as involving mainly heavy atom motion and identified finally as a wagging motion of the chromophore to adapt to the situation in the excited state. We computed our PESes from the A structure (and not A*) as this kind of motion lies outside the scope of our model. Even in this ‘‘suboptimal’’ arrangement, the endoergicity in S1 is only 3 kcal mol1 (see Table 2), so it is

Fig. 6 Time evolution of the dynamical variables: hrii (i = 1–3) (left) and hRii (i = 1–3) (right). Solid lines correspond to the perprotio isotopologue and dashed lines correspond to the perdeutero isotopologue.



View Article Online


Paper

Fig. 7 Schematic representation of the controlled deformation of the base EVB potential in S1. The solid line is a section across the PES without correction. The dashed line is a corrected profile. DE0 is the natural endoergicity (in our case, 3 kcal mol1). DE is the parameter that regulates the final exoergicity in eqn (8).

conceivable that, if a certain degree of relaxation was allowed in our 111* structure, the endoergicity of the process would be further reduced. Because of the complexity of the calculations, carrying out this relaxation with the aim of running a quantitative quantum dynamical simulation is difficult, as it is possible that the relaxation process cannot be described with a single extra coordinate. This would imply modifying the model thoroughly, recalculating the potential energy surface with this extra coordinate, and deriving the kinetic energy term. We opted instead to modify our model in a way that qualitatively reproduces the effect of this relaxation and observe the changes this induces. A relaxation that sets in after H3 moves to the products could be introduced as a correction that increases the well depth of the XX1 (X = 0, 1) potential energy region, as depicted in Fig. 7. Such an effect can be introduced functionally by including an additive term to the potential of the form: V S1 ðr1 ; r2 ; r3 ; R1 ; R2 ; R3 Þ ¼ V EVB;S1 ðr1 ; r2 ; r3 ; R1 ; R2 ; R3 Þ

(8) 1 1 DE þ tanh a r3 r03 2 2

PCCP

where the term in square brackets is just a switch function that causes the regions where r3 4 r03 to come down in energy by an amount, DE. By choosing r03 close to 0 on the positive side, we guarantee minimal distortion at the reactants’ side (including a transition structure). The values selected in our simulations are a = 7.5 a01 and r03 = +0.25 a0. Because DE0 = 3 kcal mol1 (q.v. Table 2), values of DE similar to this figure will render the proton transfer approximately isoergic. We ran simulations where the potential energy surface of S1 is given by eqn (8) with a range of values for DE starting at 2 kcal mol1 (where the process is still endoergic by just 1 kcal mol1). Fig. 8 shows the time evolution of the expectation values of r1 and r3. It can be seen that the basics of the process (stagnation of hr3i at values around 0.1 Å) are still reproduced for DE = 2 kcal mol1, which still renders an endoergic process. However, as the process becomes isoergic or exoergic (DE Z 3 kcal mol1), hr3i attains values of 0 or positive. Except for the simulation with DE = 10 kcal mol1, hr3i stays positive over 100 fs and then returns to negative. The long time dynamics reveal an emerging pattern of oscillations where hr3i goes from the reactants to products periodically with periods around 250 fs. This is indicative of the wave packet bounding back and forth without obstacles, being an example of a Rabi oscillation between two states, one centered over the 000 structure and the other over the 111 structure. Proton H1 follows a similar pattern. As DE increases over the value of DE0 (3 kcal mol1), hr1i shifts to positive values. This implies that for values of DE 4 3 kcal mol1 it can be said that the chromophore is ionized most of the time, and hence emission from I* will prevail over emission from A*. Also note that hr1i does not enter the positive region until hr3i does so too. Grigorenko et al. performed accurate optimizations on the ground and photoactive states of wtGFP.24 Concerning S1, a vertical excitation energy of 57.8 kcal mol1 (492 nm) was found for the I species and 73.7 kcal mol1 (385 nm) for the A species. The A - I process was found to imply an energy gain of 1 kcal mol1. Overall, an approximate value of 15 kcal mol1 is deduced for the A* - I* process from these data. The process

Fig. 8 Time evolution of hr1i (left) and hr3i (right) for simulations with different values of parameter DE (eqn (8)). The dashed line (for DE = 0 kcal mol1) is the result for the intact EVB potential.



View Article Online


PCCP

in our case implies an energy gain, but we highlight that the process we have studied is strictly A* - 111, by which we mean that the protons undergo translocation but the heavy atoms remain static. Thus, it seems that some amount of heavy atom motion is necessary in order to explain the exoergicity of the process in the excited state. Quantum dynamics simulations predict that, if only proton motion is allowed, protons remain essentially in a configuration similar to the one expected in A but with a large degree of delocalization, especially for H1. The behavior of this proton in our simulation is very similar to what has been reported concerning the existence of a transient LBHB in the photoactive state at the beginning of the process.15 Our QD simulations show that the expected position of the proton is very sensitive to slight modifications of the potential energy profile. Given that after photoexcitation it is expected that some heavy atom motion takes place to improve their arrangement to the conditions of the excited state, this will cause the depth of the ‘‘111’’ potential energy minimum to increase, and this fact would force then a larger localization of the protons to the neighborhood of the acceptor atom, even though this deeper minimum would be reached through a heavy-atom coordinate different from those included in the dynamical model. To explore this point, we carried out a QM/MM geometry optimization starting at the ‘‘111’’ structure on the full configurational space of the S1 electronic state, in which the active region consisted of all residues with atoms within 15 Å of the chromophore. Details of the chromophore and its environment in the S1 minimum corresponding to the product of the proton-transfer are shown in Fig. 9. First and foremost, the structure in Fig. 9 is connected directly to our ‘‘111’’ minimum without an intervening barrier (that is, the optimization drove the geometry from one to the

Paper

other along a monotonous downhill path). This is an indication that our ‘‘111’’ structure is not really a local minimum if the rest of the heavy coordinates are allowed to vary. Comparison of Fig. 3 (structure ‘‘000’’ or A) and Fig. 9 (structure I*) shows that, apart from the natural shift of the three protons toward the acceptor atoms, most distances between interacting residues do not vary significantly. Distance Ser205 Glu222 has increased from 2.85 Å to 3.03 Å, which is compatible with the fact that I* hosts a glutamic acid instead of a glutamate. The situation at the other end is more intricate: W22 is now slightly more distant from the chromophore (2.58 Å in I*, 2.66 Å in A) and closer to Ser205 (2.59 Å in I*, 2.63 Å in A), which seems odd now that the chromophore is anionic. However, this can be understood in terms of the change of position of His148: in A it was at 3.01 Å and in I* it was at 2.85 Å. This substantial decrease of the distance to the residue helps stabilize the nascent charge in the chromophore. No significant changes can be observed in the distances at which Arg96 and Thr203 interact with the chromophore. Grigorenko et al. published optimized structures for the A, I and B states of wtGFP.24 Their I structure shows distances that are longer than those we found (in particular, His148 Cro, for which a distance of 3.62 Å is reported), with the exception being Glu222 Ser65, where the intervening distance amounts to 2.72 Å. Besides the differences in methodology, these structures have to be compared with care as our optimization yields us a candidate for I*, not I. At any rate, Grigorenko et al. found a change of 0.16 Å in the His148 Cro distance in the transformation A - I, while we found the same distance for exactly the same change for the A* - I* process. It seems clear then that the His148 Cro distance is strongly coupled to the generation of the I* species from the 111 structure. Hence, a complete QD study of the A* - I* process must include, besides the set of coordinates describing the position of the protons and the distance between donors and acceptors, the His148 Cro distance at the very least. Finally, Grigorenko et al. also found that the orientation of Thr203 seemed to play no relevant role in the stabilization of the structures participating in the photocycle of GFP.24 In an analogous way to what was reported, we also found that Thr203 retains an orientation in which the side chain interacts with W22 and the chromophore. Besides, Thr203 does neither change this orientation nor vary the distance to these groups significantly along the A - I* reaction pathway.

4. Conclusions

Fig. 9 Chromophore and neighboring residues at the local potential energy minimum corresponding to the product of the triple proton transfer in S1. Distances are in Å (yellow numbers are distances between atoms in the proton-wire, and white numbers are other distances).


In this paper we have analyzed the feasibility of spontaneous proton transfer at the Franck–Condon region directly after photoexcitation. Using DFT methodology within a QM/MM framework, a structure corresponding to the A species of wtGFP has been found that agrees reasonably with the one reported recently by Grigorenko et al.24 Computation of a sizeable portion of the potential energy surface in the ground and photoactive states in the immediate surroundings of the Franck–Condon


View Article Online


Paper

region of species A (i.e. without nuclear relaxation) reveals that the process of triple proton transfer is endoergic by 3 kcal mol1 in the photoactive state and must undergo a barrier of B5.6 kcal mol1 corresponding to the transition state of an asynchronous concerted transfer. Determination of the ZPE level of the ground and excited states using quantum dynamics simulations confirms that for both electronic states, but especially the photoactive one, this energy level lies above the barrier. The expectation values of the variables describing the position of the protons indicate clearly that protons H1 (shared between the chromophore and W22) and H3 (shared between Ser205 and Glu222) will be found very close to the central region of the respective hydrogen bonds, giving support to the claim of Di Donato et al. that the first transient intermediate (I0*) detected has characteristics akin to those of a low-barrier hydrogen bond.15 However, the situation we find is not stationary and we think that this terminology is not strictly applicable in this case. Quantum dynamical simulations of the evolution of the system in the excited state show, in general terms, very scant changes with respect to the initial state. While it is true that the degree of localization of the protons is closer to the center of the respective hydrogen bonds, no definitive move to attain a closer vicinity to the acceptor atoms is observed. This is seen as an indication that the proton-transfer process is not operative if no heavy atom relaxation is permitted. The fact that the two minima of the reactants and products of the proton transfer are close to degeneracy suggests that slight rearrangements of the heavy atoms outside of the dynamic model used could change the energy balance of the process. Quantum dynamical simulations, where this energy balance is altered to bring it into degeneracy, or exoergicity, in a controlled way, result in clear transit of the protons to the acceptor side of their respective hydrogen bonds. A photoactive state geometry optimization using QM/MM methodology, including all residues with atoms within 15 Å of the chromophore and initiated at a configuration obtained taking the A minimum and moving the protons to the product side, yields a minimum energy structure with the protons transferred and in which the residue that has undergone the largest change is His148, which is now 0.16 Å closer than before the proton transfer. We interpret these results in the following way: assuming a sequential mechanism of decay of A* into I* via intermediate I0*, we find that: (1) proton transfer is not possible if structural relaxation of the surroundings of the chromophore is prevented; (2) protons H1 and H3 especially are found very close to the point halfway between the donor and acceptor after photoexcitation when the zero-point energy is considered; and (3) a geometrical parameter exists (the His148–Cro distance) under which the structure with the protons transferred is not a minimum, and that, if allowed to vary, should lead to the I* structure found by Grigorenko et al. The existence of an oscillating stationary state between reactants and products of the triple proton transfer reaction can explain the dual emission reported for the I0* intermediate of wtGFP.


PCCP

Acknowledgements This work was supported by the ‘‘Ministerio de Economıá y Competitividad’’ through projects CTQ2011-24292 and CTQ2014-53144-P. Use of computational facilities at the Con`mics de Catalunya (CSUC) is sorci de Serveis Cientı´fics i Acade gratefully acknowledged. Protein graphics on the left panel of Fig. 1, 3 and 9 were prepared with the VMD program.60 Molecular graphics on right panel of Fig. 1 were prepared with the Molekel program.61 Mathematical plots of the data were done with the GNUPLOT open source software.

References 1 2 3 4 5 6 7

8 9 10

11 12 13

14 15

16 17 18 19

20

M. Zimmer, Chem. Rev., 2002, 102, 759–781. S. R. Meech, Chem. Soc. Rev., 2009, 38, 2922–2934. J. J. van Thor, Chem. Soc. Rev., 2009, 38, 2935–2950. S. J. Remington, Protein Sci., 2011, 20, 1509–1519. W. W. Ward, C. W. Cody, R. C. Hart and M. J. Cormier, Photochem. Photobiol., 1980, 31, 611–615. M. Chattoraj, B. A. King, G. U. Bublitz and S. G. Boxer, Proc. Natl. Acad. Sci. U. S. A., 1996, 93, 8362–8367. J. T. M. Kennis, D. S. Larsen, I. H. M. van Stokkum, M. Vengris, J. J. van Thor and R. van Grondelle, Proc. Natl. Acad. Sci. U. S. A., 2004, 101, 17988–17993. ¨, A. B. Cubitt, K. Kallio, L. A. Gross, R. Y. Tsien and M. Ormo S. J. Remington, Science, 1996, 273, 1392–1395. F. Yang, L. G. Moss and G. N. Phillips Jr., Nat. Biotechnol., 1996, 14, 1246–1251. K. Brejc, T. K. Sixma, P. A. Kitts, S. R. Kain, R. Y. Tsien, ¨ and S. J. Remington, Proc. Natl. Acad. Sci. U. S. A., M. Ormo 1997, 94, 2306–2311. D. Stoner Ma, A. A. Jaye, P. Matousek, M. Towrie, S. R. Meech and P. J. Tonge, J. Am. Chem. Soc., 2005, 127, 2864–2865. J. J. van Thor, G. Zanetti, K. L. Ronayne and M. Towrie, J. Phys. Chem. B, 2005, 109, 16099–16108. D. Stoner Ma, E. H. Melief, J. Nappa, K. L. Ronayne, P. J. Tonge and S. R. Meech, J. Phys. Chem. B, 2006, 110, 22009–22018. C. Fang, R. R. Frontiera, R. Tran and R. A. Mathies, Nature, 2009, 462, 200–204. M. Di Donato, L. J. G. W. van Wilderen, I. H. M. van Stokkum, T. Cohen Stuart, J. T. M. Kennis, K. J. Hellingwerf, R. van Grondelle and M. L. Groot, Phys. Chem. Chem. Phys., 2011, 13, 16295–16305. J. J. van Thor, K. L. Ronayne, M. Towrie and J. T. Sage, Biophys. J., 2008, 95, 1902–1912. W. W. Cleland and M. M. Kreevoy, Science, 1994, 264, 1887–1890. C. L. Perrin, Acc. Chem. Res., 2010, 43, 1550–1557. G. T. Hanson, T. B. McAnaney, E. S. Park, M. E. P. Rendell, D. K. Yarbrough, S. Chu, L. Xi, S. G. Boxer, M. H. Montrose and S. J. Remington, Biochemistry, 2002, 41, 15477–15488. I. H. M. van Stokkum, B. Gobets, T. Gensch, F. van Mourik, K. J. Hellingwerf, R. van Grondelle and J. T. M. Kennis, Photochem. Photobiol., 2006, 82, 380–388.


View Article Online


PCCP

21 X. Shu, K. Kallio, X. Shi, P. Abbyad, P. Kanchanawong, W. Childs, S. G. Boxer and S. J. Remington, Biochemistry, 2007, 46, 12005–12013. 22 X. Shi, P. Abbyad, X. Shu, K. Kallio, P. Kanchanawong, W. Childs, S. J. Remington and S. G. Boxer, Biochemistry, 2007, 46, 12014–12025. 23 D. Stoner Ma, A. A. Jaye, K. L. Ronayne, J. Nappa, S. R. Meech and P. J. Tonge, J. Am. Chem. Soc., 2008, 130, 1227–1235. 24 B. L. Grigorenko, A. V. Nemukhin, I. V. Polyakov, D. I. Morozov and A. I. Krylov, J. Am. Chem. Soc., 2013, 135, 11541–11549. 25 M. Nadal-Ferret, R. Gelabert, M. Moreno and J. M. Lluch, J. Am. Chem. Soc., 2014, 136, 3542–3552. 26 S. Yamaguchi, H. Kamikubo, K. Kurihara, R. Kuroki, N. Niimura, N. Shimizu, Y. Yamazaki and M. Kataoka, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 440–444. 27 O. Vendrell, R. Gelabert, M. Moreno and J. M. Lluch, J. Am. Chem. Soc., 2006, 128, 3564–3574. 28 O. Vendrell, R. Gelabert, M. Moreno and J. M. Lluch, J. Phys. Chem. B, 2008, 112, 5500–5511. 29 O. Vendrell, R. Gelabert, M. Moreno and J. M. Lluch, J. Chem. Theory Comput., 2008, 4, 1138–1150. 30 O. Vendrell, R. Gelabert, M. Moreno and J. M. Lluch, J. Phys. Chem. B, 2008, 112, 13443–13452. 31 O. Vendrell, R. Gelabert, M. Moreno and J. M. Lluch, Chem. Phys. Lett., 2004, 396, 202–207. 32 H. Li, A. D. Robertson and J. H. Jensen, Proteins: Struct., Funct., Bioinf., 2005, 61, 704–721. 33 D. C. Bas, D. M. Rogers and J. H. Jensen, Proteins: Struct., Funct., Bioinf., 2008, 73, 765–783. 34 M. H. M. Olsson, C. R. Sondergard, M. Rostkowski and J. H. Jensen, J. Chem. Theory Comput., 2011, 7, 525–537. 35 C. R. Sondergard, M. H. M. Olsson, M. Rostkowski and J. H. Jensen, J. Chem. Theory Comput., 2011, 7, 2284–2295. 36 ChemShell, a Computational Chemistry Shell, see www.chem shell.org, last accessed on April 10, 2015. 37 A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiorkiewicz Kuczera, D. Yin and M. Karplus, J. Phys. Chem. B, 1998, 102, 3586–3616. 38 A. D. MacKerell, M. Feig and C. L. Brooks, J. Am. Chem. Soc., 2004, 126, 698–699. 39 R. G. Parr and W. Yang, Density Functional Theory of Atoms and Molecules, Oxford University Press, Oxford, UK, 1989. 40 A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652. 41 C. T. Lee, W. T. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789. 42 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng,


Paper

43 44 45 46 47

48 49 50 51 52 53 54

55 56 57 58 59 60 61

J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. ¨ . Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski Daniels, O and D. J. Fox, Gaussian 09, Revision C.01, Gaussian Inc., Wallingford, CT, 2009. T. Yanai, D. P. Dew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57. A. L. Sobolewski and W. Domcke, Phys. Chem. Chem. Phys., 1999, 1, 3065–3072. N. H. List, J. M. Olsen, T. Rocha-Rinza, O. Christiansen and J. Kongsted, Int. J. Quantum Chem., 2012, 112, 789–800. ¨stner, J. M. Carr, T. W. Keal, W. Thiel, A. Wander and J. Ka P. Sherwood, J. Phys. Chem. A, 2009, 113, 11856–11865. ¨ckle, H.-D. Meyer, M.-C. Heitz, G. A. Worth, M. H. Beck, A. Ja S. Wefing, U. Manthe, S. Sukiasyan, A. Raab, M. Ehara, C. Cattarius, F. Gatti, F. Otto, M. Nest, A. Markmann, M. R. Brill and O. Vendrell, The MCTDH Package, Versions 8.3, 8.4, 2000–2007, see http://www.pci.uni-heidelberg.de/tc/ usr/mctdh. R. P. Bell, The Tunnel Effect in Chemistry, Springer, 1980. A. Warshel, J. Phys. Chem., 1982, 86, 2218–2224. H.-D. Meyer, U. Manthe and L. S. Cederbaum, Chem. Phys. Lett., 1990, 165, 73–78. U. Manthe, H.-D. Meyer and L. S. Cederbaum, J. Chem. Phys., 1992, 97, 3199–3213. ¨ckle, G. A. Worth and H.-D. Meyer, Phys. M. H. Beck, A. Ja Rep., 2000, 324, 1–105. H.-D. Meyer and G. A. Worth, Theor. Chem. Acc., 2003, 109, 251–267. Multidimensional Quantum Dynamics: MCTDH Theory and Applications, ed. H.-D. Meyer, F. Gatti and G. A. Worth, Wiley-VCH, Weinheim, 2009. H.-D. Meyer, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 351. P. Armengol, R. Gelabert, M. Moreno and J. M. Lluch, J. Phys. Chem. B, 2015, 119, 2274–2291. C. Scharnagl and R. A. Raupp Kossmann, J. Phys. Chem. B, 2004, 108, 477–489. K. Saito and H. Ishikita, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 167–172. M. Fayngold and V. Fayngold, Quantum Mechanics and Quantum Information, Wiley-WCH, 2013. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38. U. Varetto, MOLEKEL, Version 5.4.0.8, Swiss National Supercomputing Centre, Lugano, Switzerland, 2009.


Ground-state proton transfer kinetics in green fluorescent protein.

Vibrational analysis on the revised potential energy curve of the low-barrier hydrogen bond in photoactive yellow protein.

Short Hydrogen Bonds and Proton Delocalization in Green Fluorescent Protein (GFP).

Excited state dynamics of the isolated green fluorescent protein chromophore anion following UV excitation.

A Potential Energy Function for Heterogeneous Proton-Wires. Ground and Photoactive States of the Proton-Wire in the Green Fluorescent Protein.

Are there really low-barrier hydrogen bonds in proteins? The case of photoactive yellow protein.

Incomplete proteasomal degradation of green fluorescent proteins in the context of tandem fluorescent protein timers.

Local fitness landscape of the green fluorescent protein.

ZnO nanoparticles assist the refolding of denatured green fluorescent protein.

Primary structure of the Aequorea victoria green-fluorescent protein.

Neutron crystallography of photoactive yellow protein reveals unusual protonation state of Arg52 in the crystal.

Correction: Multiscale Approach to the Determination of the Photoactive Yellow Protein Signaling State Ensemble.

A photoactivatable green-fluorescent protein from the phylum Ctenophora.

Proton transfer reactions and hydrogen-bond networks in protein environments.

Hydrogen bond stereochemistry in protein structure and function.

Multiscale approach to the determination of the photoactive yellow protein signaling state ensemble.

Regulation of protein-ligand binding affinity by hydrogen bond pairing.

Spatial dependence of protein-water collective hydrogen-bond dynamics.

Enhanced green fluorescent protein-mediated synthesis of biocompatible graphene.

Anaerobic green fluorescent protein as a marker of Bifidobacterium strains.

Transgenic pigs with pancreas-specific expression of green fluorescent protein.

Highly Fluorescent Green Fluorescent Protein Chromophore Analogues Made by Decorating the Imidazolone Ring.

Strong ionic hydrogen bonding causes a spectral isotope effect in photoactive yellow protein.

TDDFT study on the excited-state proton transfer of 8-hydroxyquinoline: key role of the excited-state hydrogen-bond strengthening.