Unveiling How an Archetypal Fluorescent Protein Operates

Subscriber access provided by NATIONAL SUN YAT SEN UNIV

Article

Unveiling How an Archetypal Fluorescent Protein Operates: Theoretical Perspective on the Ultrafast Excited-State Dynamics of GFP Variant S65T/H148D Pau Armengol, Ricard Gelabert, Miquel Moreno, and José M. Lluch J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/jp506113g • Publication Date (Web): 21 Aug 2014 Downloaded from http://pubs.acs.org on August 25, 2014

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Unveiling How an Archetypal Fluorescent Protein Operates: Theoretical Perspective on the Ultrafast Excited-State Dynamics of GFP Variant S65T/H148D Pau Armengol,† Ricard Gelabert,∗,† Miquel Moreno,† and José M. Lluch†,‡ Departament de Qu´ımica, Universitat Aut` onoma de Barcelona, 08193 Bellaterra, Barcelona, Spain, and Institut de Biotecnologia i de Biomedicina, Universitat Aut` onoma de Barcelona, 08193 Bellaterra, Barcelona, Spain E-mail: [email protected]

Abstract

and products of the proton-transfer reaction in ground and most importantly, photoexcited state, in terms of potential energy profiles associated to the proton migration. Under these conditions, proton transfer can occur in accord to experimental data available. This set of characteristics are possibly common to a host of other proton-transfer based fluorescent proteins, and helps promoting GFP S65T/H148D to a case of archetypal significance. Thus, our results can be useful to understand the way many fluorescent proteins work and, more generally, the molecular basis for proton-transfer reactions in proteins.

Green Fluorescent Protein variant S65T/H148D has been reported to host a photocycle involving the photoinduced proton-transfer reaction between the chromophore and residue Asp148 under 50 fs and without measurable kinetic isotope effect, and experimental evidence is suggestive of the existence of a highly delocalized proton between these residues. The blinding speed at which this biological system undergoes proton transfer has been ascribed to the extreme increase of acidity of the GFP chromophore in the electronic excited state where proton transfer takes place. This work strives to present a coherent, complete and balanced description of the dynamics of this specific variant of GFP in which it will be shown that this increase of acidity is insufficient to explain the behavior observed. This study tracks the behavior of this photosystem to the delicate interplay between structure and dynamics shown in presence of solvent. In this way, it has been found that the dynamics of this protein intertwines its structure with the intervening solvent to give rise to effectively degenerate situations in what concerns the reactants

Keywords: QM/MM, Molecular Dynamics, Quantum Dynamics, Ergodic Hypothesis, Proton Transfer

INTRODUCTION Green Fluorescent Protein (GFP) can be considered a paradigm of fluorescent proteins. In part, this is so because of the intrinsic interest of the system, as it has served to establish a model of photochromism in fluorescent proteins. 1,2 But it is maybe because of the many applications in biomedicine that depend on the properties of GFP and fluorescent proteins in

∗

To whom correspondence should be addressed Departament de Qu´ımica ‡ Institut de Biotecnologia i de Biomedicina †

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

general 3–6 that these have become famous beyond the science community. 7 Fluorescent proteins nowadays are indispensible resources in imaging techniques. This has resulted in enormous efforts being put into the study and development of novel natural or engineered fluorescent proteins with improved characteristics better suited to imaging applications. The basic photocycle of GFP has been studied extensively and there is consensus about its main aspects. 8 At room temperature, the absorption spectrum of GFP presents two bands, the most intense one is conventionally called A (398 nm) and the less intense B (471 nm). Evidence exists that A and B correspond to the neutral and ionized forms of the chromophore, p-hydroxybenzylideneimidazolinone (p-HBDI). This chromophore is formed in the post-traslational cyclization reaction of the Ser65-Tyr66-Gly67 tripeptide, followed by a 1,2-dehydrogenation of Tyr66. Upon excitation of the A band, strong green fluorescence is registered at 508 nm (species I*) within picoseconds, which makes the process at the core of GFP rank among the fastest processes known in biological systems. It is agreed that I/I* corresponds, like B, to an anionic form of the chromophore, but differs in that the protein environment of the chromophore is not yet adapted to the now anionic species. Transformation of I* into B* is more difficult, and already in the ground state, converting B into A is a process spanning hours. The mechanism by which fluorescence from I* occurs involves excited-state proton transfer (ESPT). This has been established through the determination of a large (∼ 5) kinetic H/D isotope effect in the rise of fluorescence from I*. The proton-transfer has been proposed, based on evidence from the crystal structure of GFP, to involve the net transfer of a proton through a three-link proton wire that involves the chromophore, a water molecule, and residues Ser205 and Glu222. Time-resolved transient infrared experiments have validated this pathway in capturing the protonation of Glu222 side chain. 9–11 Mutagenesis belongs to the toolbox of useful techniques in the subtle art of elucidation of structure-property relationships in proteins.

Page 2 of 32

In the case of wtGFP, a mutant that has received particular attention is generated by substitution of Ser65 by a threonine (S65T). The steric hindrance caused by the extra methyl group manifested itself in a remarkably different absorption spectrum, consisting solely of a peak at 489 nm. 12 When the crystal structure of wtGFP 13,14 and S65T 15 were published, it was seen that both differed mainly in the protonation state of the chromophore, neutral in wtGFP and ionized in S65T. Thus, S65T represented one of the first milestones down a path that has eventually led to the current knowledge on the mechanism of fluorescence in wtGFP. Almost twenty years ago, Remington and coworkers suggested that fluorescent protein BFP (the double mutant Y66H/Y145F) could be used as pH reporter in cells and cellular compartments. 16,17 This suggestion spurred studies to analyze response to pH of GFP-related proteins and the development of novel variants to undertake the measurement of pH in cytoplasm and the inside of organelles with non-invasive techniques. 18,19 Even though GFP’s chromophore has acidbase properties, it was noted long ago that wtGFP is not very responsive to changes of pH of the medium. If novel GFP variants could be devised in which the chromophore were not so unresponsive to bulk pH values, and on top it were possible to modulate the pKa value of titratable groups in the chromophore, it would be possible to monitor pH values in cell compartments with widely different acidity conditions. In this regard, the absorption and fluorescence spectra of the mutant GFP S65T show a remarkable pH dependence which is a desirable property of such potential pH indicators. 17 In the attempt to prepare GFP variants retaining the pH dependence of S65T but with different chromophore pKa values for potential use as indicators, Elsliger et al. targeted the effect of His148, which in both wtGFP and S65T is directly H-bonded to the chromophore and hinders access of the solvent to it, 17 and mutated it to aspartate to explore the effect of a proximal negative charge. The absorption spectrum of GFP S65T/H148D displays two peaks,


2

Page 3 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for the specific observed behavior. Naturally enough, the very short contact between Asp148 and the chromophore has received a lot of consideration. Absorption of the A band is red shifted by about 20 nm with respect to its counterpart in wtGFP, and it was suggested that this red shift was caused by the short hydrogen bond between chromophore and Asp148, 20 which could promote further delocalization of the π-system of the chromophore over the carboxylate group. The extreme velocity at which I* is seen to generate has motivated speculations as to the influence of the very short distance present between Asp148 and the chromophore. In fact, it was already suggested by Remington and coworkers 20–22 that a situation in which the proton were already evenly shared could explain the velocity of the process: If a low-barrier hydrogen bond (LBHB) were to exist between the chromophore and Asp148, the proton would be more evenly shared in the ground state and photoexcitation might cause an almost instantaneous transfer when the hydrogen bond interaction were suddenly weakened in the excited state. To explore this possibility, 1 H NMR analysis of GFP S65T/H148D focusing on the lowfield region of the spectrum was done, but no signal indicative of an anomalously deshielded proton (one of the tell-tale indications of possible involvement of an LBHB) could be found. 20 To further explore the excited state dynamics with enough resolution to unveil the formation and cleavage of individual bonds, Stoner-Ma et al. investigated with time-resolved infrared spectroscopy GFP S65T/H148D. 23 Yet another unexpected result surfaced, as within their experimental time resolution (∼400 fs) no change in the ionization state of a carboxylate residue could be measured as I* was being formed, and only dissipative dynamics were detected during the process, even though, as noted, green fluorescence from I* was obtained instantaneously within instrument response. The authors rationalized this observation arguing that if the transferrable proton is already shared between donor and acceptor before photoexcitation, no complete change in the carboxylate could be observed when the transfer takes place. This sit-

A (415 nm) and B (487 nm), both red-shifted with respect to their wtGFP counterparts, with good pH dependence and a pKa value of 7.8, about two units higher than its parent S65T (pKa =6.0). Nevertheless, this variant had more in stock than could be anticipated: both absorption peaks could be excited to trigger green fluorescence at 510 nm (species I*) with high quantum yield (0.21, to be compared to 0.02 for the parent S65T protein, in which the protonated species has a low yield and does not undergo excited-state proton transfer). Somehow, the proton-transfer mechanism operative in wtGFP, which was disrupted in GFP S65T, had been rewired. Elsliger et al. hinted at the potential uses of S65T/H148D as a fluorescent ratiometric pH indicator for in vivo monitoring in real time and with spatial resolution. 17 A complicated photophysical behavior was uncovered as GFP S65T/H148D was scrutinized in a series of thorough mutagenic and spectroscopic studies by Remington and coworkers. 20–22 The most striking discovery was that the rise of green fluorescence at 510 nm was so fast as to lie below the experimental time resolution (< 170 fs), and on top with rise time insensitive to H/D substitution. This is in sharp contrast to what is found for wtGFP, where green fluorescence appears within picoseconds and is very sensitive to deuteration. 8 The fact that low-temperature steady-state emission experiments detected no measurable rise time for the I* fluorescence even at very low temperatures gave support to the picture of an essentially barrierless ESPT mechanism. 22 The crystal structure of S65T/H148D was resolved 20 and is shown in Figure 1. The structure shows a large cleft separating the β-sheet supporting Asp148 from the neighboring one which might indicate a route of access of the solvent to the chromophore area and explain the larger pH sensitivity of this mutant. The structure was shown to contain a very short contact between one of the carboxylic oxygens of Asp148 and the phenolic oxygen of the chromophore (2.32 ˚ A). This evidenced that Asp148 is the final acceptor of the chromophore’s proton and that the old, wtGFP-like proton-wire, though present, could not be the main reason


3


Tyr145

4.3

0

W1

Cro 2.32

2 .9 5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 32

Asp148 W3

*

2.72

Ser205

Figure 1: Detail of the X-ray diffraction structure of GFP S65T/H148D. Left panel: detail of the region of the protein harboring the chromophore and Asp148, where residues are color-coded: yellow corresponds to the original wtGFP proton-wire member Ser205, brown to the main actors in S65T/H148D (chromophore, Asp148 and on top, Tyr145), and red isolated atoms to water molecules. Right panel: enlargement of the X-ray structure (PDB ID: 2DUF) with most important distances noted in ˚ A. The asterisk denotes the water molecule that belongs to the “old” proton-wire in wtGFP. uation is what could be expected if an LBHB were involved in the ground state, for instance. In 2010, and with much improved time resolution due to changes in the experimental setup (established now at ∼50 fs), time resolved fluorescence up-conversion experiments were carried out by Kondo et al. to try to capture the motion of the proton as it moves from chromophore to Asp148 by monitoring emission under the I* band at different wavelengths. 24 But yet, again no rise in I* fluorescence could be measured even in the red edge of I* and no effects on the decay rates could be observed upon deuteration. The ESPT in this GFP mutant occurs with a portion of the long-lived ionized I* state formed under 50 fs, and slower relaxation processes were also observed with time constants in the range 1-10 ps. The proposal of Kondo et al. contemplates barrierless ultrafast ESPT within 50 fs followed by slower vibrational cooling along the proton-transfer coordinate and relaxation of the protein environment to the new situation of the chromophore. The question of what can unleash such an extremely fast process was rationalized in the possible existence of a short hydrogen bond, but especially in the strong photoacid character of

the chromophore, estimated to be as low as 0.1 in the excited state by Scharnagl and RauppKossmann in a theoretical study based on a statistical-thermodynamical approach 25 (to be compared with a value of 8.2 for the isolated chromophore in solution in the ground state 26 ). The extreme increase of acidity of the chromophore in the excited state: (a) would place the proton equilibrium position right next to Asp148, not the chromophore, after photoexcitation; (b) this would cause the transfer process to be very exoergic; and (c) together with the tight distance between Asp148 and chromophore, would cause the process to be barrierless and very fast. Paraphrasing the authors in Ref. 23, the picture just discussed is, in effect, compatible with experimental evidence available. But direct allusion has been made to concepts like low-barrier hydrogen bonds and certain characteristics of the ground and excited electronic states that “... need to be supported by detailed high-level quantum chemical calculations. Even if such surfaces become available (in itself, a challenging objective), the subsequent calculation of the ultrafast rate coefficients will require nonequilibrium dynamics simulations. 23 ”


4

Page 5 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In this contribution, we pick up the gauntlet and have used a wide spectrum of state-ofthe-art methods and techniques in theoretical chemistry to provide a penetrating insight with atomic resolution on the inner workings of GFP S65T/H148D. We take off from realistic representations of the complete protein and environment and using accurate quantum-mechanical molecular-mechanical (QM/MM) methodology to provide a realistic description of the protontransfer coordinate both in the ground state and excited electronic state where protontransfer occurs; then performing long and accurate QM/MM molecular dynamics (MD) simulations to analyze and rationalize the structure and behavior of the environment of the chromophore in constant temperature situations; and finally going all the way to provide estimates of the proton-transfer times using Quantum Dynamics (QD) simulations. We have found that the picture provided so far of the dynamics of GFP S65T/H148D is correct in general lines, 20–24 except, very relevantly, where driving force of the protontransfer is concerned. In this regard, we find that, for the excited state proton-transfer to take place, shadowy characters need to step forth onto the stage: solvent molecules and a companion from inside the protein: residue Tyr145. We will discuss that this GFP variant, maybe because of its intrinsic simplicity in comparison to most other members of its group (as it involves a single proton-transfer), crystallizes some particularities that to a varying extent are also present in the original GFP. In particular, the fact that the system has been “designed” to attain an almost isoergic profile in the ground state and the excited state where proton-transfer occurs, as concerns the protontransfer coordinate, seems to be a common trait of this family of proteins.

vironment) have been performed with the Gaussian09 software suite. 27 Ground and excited state energies have been obtained using Density Functional Theory (DFT) 28 and Time-Dependent DFT (TDDFT), 29,30 respectively, using the 6-31G(d,p) basis set and the Coulomb attenuated functional CAMB3LYP. 31 This functional mixes Hartree-Fock and DFT exchange-correlation contributions according to a distance-dependent ratio, and has been shown to circumvent the dramatic failure of pure functionals when computing excitation energies of excited electronic states of charge-transfer type. It has also demonstrated good performance in describing excitation energies of the GFP chromophore 32 and other systems. 33 Precision of the computed excitation energies are established at ∼0.5 eV with this methodology. 34 When necessary, the effect of solvation has been introduced via the polarizable continuum method (PCM). 35 Optimizations in this work have been done following the gradient in the ground electronic state until a stationary point with all-positive eigenvalues of the Hessian matrix is reached.

Static QM/MM Calculations with Protein Environment Initial coordinates were obtained from the Xray resolved structure deposited in the Protein Data Bank (PDB ID: 2DUF). 20 Missing atoms (mainly hydrogens) were added using the PROPKA server at pH 5.7. 36–39 To realistically represent the system in operational conditions, the protein was placed inside a water droplet of 70 ˚ A diameter. Energies (of ground and excited states) and gradients (for ground-state optimizations) have been computed in a quantum-mechanical molecular-mechanical (QM/MM) framework. Coupling of the QM and MM regions is achieved through the electronic embedding (EE) scheme, in which the force field charges (environment) are allowed to polarize the QM portion. In the current case, the QM part of the system is comprised of the chromophore (formed by the cyclization reaction of Ser65, Tyr66 and Gly67), Asp148 and part of a few of the closest

METHODOLOGY Cluster Calculations Pure quantum-mechanical (QM) calculations of cluster type (that is, without protein en-


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

residues to both (Val68, Ser147 and Asn149), for a total of 58 QM atoms. The rest of the protein, and the crystallization and solvent water molecules went into the MM region, for a total of 18,508 MM atoms. Calculations were done within a QM/MM framework using the ChemShell suite. 40,41 Gaussian 09 was used 27 to describe the QM part of the system with the same choice of functional and basis set as in the cluster calculations, while the CHARMM22 force field was employed 42,43 for the MM part via the DLPOLY ChemShell device. 44 For excited state calculations, TDDFT/CAM-B3LYP was the method of choice in the QM part. Link atoms were used to connect the QM and MM parts of the system. The QM/MM interaction energy was determined with the charge-shift approach, in which the QM/MM electrostatic interactions are worked out by the QM code including polarization originating in the charge distribution of the MM part of the system. MM charges close to link atoms are shifted away and point dipoles are added to compensate. Single-point calculations were done after minimal relaxation to eliminate bad contacts between water molecules and the protein introduced when preparing the solvated model. Optimization in the ground electronic state was performed with the HDLCOpt module of ChemShell. 45 The active region (the atoms allowed to move during optimizations) comprised the full protein and all water molecules within 4 ˚ A of it. A thin layer of water molecules in the outside of the droplet were, then, held fixed to prevent evaporation but the protein was not restricted in any other way.

Page 6 of 32

in the field we must turn to semiempirical Hamiltonians to achieve reasonable accuracy and avoid the superlative costs of more rigorous approaches. In this, Thiel’s orthogonalization model (OMx) methods have been developed in the last decade and proved to provide accurate results for ground-state calculations of biochemical systems. 48,49 The QM method chosen for our study is OM3, 50 as implemented in the MNDO99 package, 51 and the CHARMM22 force field was used for the MM region. 42,43 Our approach requires the computation of a precursor trajectory from which to sample starting configurations (“seeds”) for several independent short trajectories on which averaging is to be performed. The precursor trajectory is obtained from the following protocol: the system starts off from the structure used for point QM/MM calculations. Stochastic boundary conditions are used to control the temperature of the simulation, and SHAKE is used to keep distances between hydrogens and their bound atoms constant. The heating starts by setting the temperature to 10 K and then the equations of motion are solved. Temperature is increased in steps of 10 K for every 100 fs of simulation, until 300 K is reached. After full heating is completed, an equilibration phase begins, in which temperature is kept constant and the trajectory is integrated for 10 ps. After this, equilibration is assumed and the precursor trajectory is finally started and obtained by continuing integration for 20 ps of production phase at constant T . The trajectories for analysis start off independent snapshots from the 20 ps production phase of the precursor trajectory, taken at regular intervals of 800 fs to ensure that these samples are not correlated. Each of these “seeds” is assigned random velocities from a Gaussian distribution corresponding to a temperature of 300 K and sent each to a single processor machine to be integrated. As the velocities assigned do not correspond to an equilibrated trajectory, each undergoes a short equilibration phase to ensure that velocities are properly randomized, lasting 3 ps. After that, the trajectory is assumed to be equilibrated and the production begins, lasting for 40 ps. As 25 such trajectories have been

Dynamical QM/MM Calculations We have devised a simulation strategy that allows us to exploit the ergodic hypothesis 46,47 by averaging over many shorter trajectories, each using a single machine. Refer to the Results and Discussion section for details. The question remains as to what QM methodology is to be used. Ab initio- or DFT-quality QM methods are out of the question as long dynamics are sought, so as is commonplace


6

Page 7 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


started and each has run from 40 ps, the total simulation time is 1.0 ns. On this superset of 25 trajectories dynamical averages can be computed as needed.

the Hamiltonian is time-independent, the timeevolved wave function can be written in equivalent form using the time-evolution operator in spectral form as 53

Quantum-Dynamical Wave Packet Simulations

Ψ(q, t) = hq|e−iHt/~ |ψ0GS i X (2) = hq|ψiPS ihψiPS |ψ0GS ie−iEi t/~

ˆ

i

Two different types of QD simulations have been done in this work: determination of stationary states of the vibrational motion (vibrational energies and wave functions), and timeevolution of the photoexcited wave packet. We deal with them in turn. Vibrational energies and wave functions have been obtained for selected potential energy profiles, each corresponding to the transfer of a proton between the acceptor and donor atom. To do this, the time-independent Schrödinger equation for the motion of the proton has been solved: h i ˆ i = Tˆ + Vˆ ψi = Ei ψi Hψ (1)

where the summation extends over all vibrational states of the excited electronic state (|ψiPS i), and {Ei } are the eigenvalues of the excited state eigenstates. In the summation in Equation 2 the first bracket is the ith vibrational wave function of the excited state, while the second is the overlap between the ith vibrational wave function of the excited state and the vibrational ground state of the ground electronic state (i.e. the Franck-Condon factor). Convergence of the time evolution of the wave packet was verified to be unitary to within one part over 105 or better. Expectation values can be computed on Ψ(q, t) to yield time-dependent observables.

where Vˆ is the potential energy operator (the aforementioned potential energy profile, computed at QM/MM level) and Tˆ is the kinetic energy operator for the motion of H(D) in one dimension. The sync-Discrete Value Representation (sync-DVR) of Colbert and Miller has been used to build the matrix representation of the Hamiltonian in Equation 1, using 50 points as basis set. 52 Eigenvalues ({Ei }) and eigenvectors ({ψi }) were obtained through diagonalization of the Hamiltonian matrix. Convergence of the results with respect to the size of the basis set was established. Time evolution of a wave packet representing the state of the proton (deuteron) in the excited state was achieved in two steps. First, the vibrational eigenstates of the ground electronic state were obtained as detailed above. Then, assuming that only the ground vibrational state in the ground electronic state (|ψ0GS i) is populated at room temperature, it (i.e. the system) is promoted to the excited state. Because it is no longer an eigenstate of the Hamiltonian, the time evolution is determined by solving the time-dependent Schrödinger equation. Because

RESULTS The wealth of data on the ultrafast dynamics in GFP variant S65T/H148D has found a plausible explanation relying on the following assumptions: (1) That the process of protontransfer in the ground state is unfavorable, (2) that a large increase of acidity is induced by photoexcitation of the chromophore, which (3) causes the proton transfer to be favorable in this excited state, to the extent of (4) being barrierless. These assumptions would explain the sudden rise (τ < 50 fs) of fluorescence of the I* band without measurable H/D kinetic isotope effects. Besides, the lack of evolution of the carboxylic bond signature in transient infrared experiments could be explained if the reactant in the ground state had characteristics akin to a very delocalized hydrogen bond, as is the case, for instance, in LBHBs. Much of the proposed explanation on photochemistry of S65T/H148D relies on the suspected increase of acidity of the GFP chro-


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

mophore after photoexcitation. A sudden increase in acidity of the chromophore should stabilize the potential energy minimum of the photoproduct (such that in comparison to the ground state, where the chromophore is found to be neutral, the process of ionization of the chromophore is found to be spontaneous). As pointed out by Kondo et al., 24 this assumption is reasonable in view of the expected electronic redistribution, but has been quantified using theoretical methods on HBDI in solution by Scharnagl and Raupp-Kossmann to drop by more than 8 pKa units, from 8.3 in the ground state to 0.1 in the photoactive state, 25 where the term “photoactive state” is introduced to refer to the excited electronic state where proton-transfer can occur. The proton acceptor in this process is Asp148, and aspartic acid has a solution pKa of about 3.7. Thus, a plunge of the chromophore pKa of this magnitude would cause the process in the excited state to be spontaneous and, very likely, barrierless. Such a dramatic change in value leads us to expect a well-formed potential energy minimum in the ground state corresponding to a neutral chromophore, and also a wellformed potential energy minimum in the photoactive state defining a ionized chromophore, likely missing a secondary minimum describing a neutral chromophore. Potential energy profiles can be computed routinely for most systems of small and medium size within a cluster approach. However, in the present case the environment of the chromophore and Asp148 is highly structured, with groups supporting localized charges that can act polarizing the chromophore. Thus, polarization by the structured environment of the chromophore cavity is likely to be nonnegligible. Because of this, we have chosen to perform QM/MM potential energy scans corresponding to the proton-transfer coordinate in which chromophore, Asp148 and some close residues are treated quantum-mechanically, and the rest of the system is described by means of a force field. With “the rest of the system” we mean the rest of the protein as well as solvent, which is introduced explicitly using a water droplet of 70 ˚ A diameter and consisting

Page 8 of 32

of 4995 water molecules (counting in those of crystallographic origin). These potential energy profiles have been computed twice, first using the geometry of the protein as provided by PDB entry 2DUF without relaxation, except minimal MM optimization of the solvent to eliminate bad contacts, and second, we minimized the energy of the protein and water droplet, and then computed unrelaxed potential energy profiles on the resulting structure. Figure 2 shows the structure of the chromophore cavity after optimization, while the profiles are shown in Figure 3. From now on we will use the abbreviation Cro to refer to the chromophore. Ground state profiles in Figure 3 show a single potential energy minimum with a O(Cro)–H distance of ∼0.98 ˚ A (unoptimized), or ∼0.99 ˚ A if the full system is allowed to relax. If optimization is allowed, the structure relaxes in such a way that the CroAsp148 O–O distance rises to 2.64 ˚ A, the chromophore attains almost perfect planarity and Asp148 adopts an arrangement closer to perpendicularity with respect to the phenolic ring of the chromophore. We will discuss the significance of this minimum energy structure when we present the results of the OM3 QM/MM dynamics simulations. Neither in the crystal structure nor in the optimized one do we observe a secondary minimum corresponding to the proton-transfer product. S0 profiles show a certain degree of anharmonicity but not very large deviations of the position of the proton from the potential energy minimum are expected because of this. Thus, it is reasonable to define the situation as follows: the proton is well localized around the chromophore’s phenolic oxygen. This casts some doubts on the supposition that a low-barrier hydrogen bond be involved in the ground electronic state of this GFP mutant. Photoactive state profiles (in green in Figure 3) also do not show the expected behavior. Both unrelaxed and relaxed (on S0 ) show a single minimum almost exactly on top of the S0 minimum and corresponding also to a chromophore-bound proton. Excitation wavelength is about 408 nm for the unrelaxed profile (crystallographic structure) which is in very


8

Page 9 of 32

4.68

Tyr145

W1

Cro

W2

2.

2.64

85

* Asp148

2.77

W3

Ser205

Figure 2: Details of the DFT QM/MM optimized structure of GFP S65T/H148D. Left panel: pictorial representation of the system object of the QM/MM study, including a solvent droplet of 70 ˚ A diameter. Right panel: zoom in to the chromophore cavity at the QM/MM optimized structure, with most relevant distances indicated in ˚ A. The asterisk identifies the water molecule that belongs to the original proton-wire in wtGFP. 90

90

85

85

80

80

75

75

70

70 E / kcal mol

65

−1

S0 S1 1 ( ππ*) S2 S3

−1

E / kcal mol

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


20

20

15

15

10

10

5

5

0

S0 1 ( ππ*) S1

65

0 0.8

0.9

1.0 1.1 1.2 dCro−H / Å

1.3

1.4

0.8

0.9

1.0 1.1 1.2 dCro−H / Å

1.3

1.4

Figure 3: Potential energy profiles of the proton-transfer in wtGFP variant S65T/H148D in the electronic ground (DFT QM/MM) and excited (TDDFT QM/MM) states. Left panel: unrelaxed potential energy profiles in which the system is kept at the same geometry reported in PDB entry 2DUF, but an explicit droplet of solvent of 70 ˚ A diameter surrounds the protein. The three lowestlying excited electronic states are also shown, but only the 1 ππ ∗ (S2 ) state is optically accessible, the rest having negligible oscillator strengths. A potential diabatic state crossing is detected around A). Right panel: (unrelaxed) potential energy profiles of the A (dAsp148−H ∼ 1.2 ˚ dCro−H ∼ 1.2 ˚ ground and lowest-lying excited states carried on the S0 -relaxed minimum of the protein+water droplet system. For comparison purposes, dotted lines represent the potential energy profiles on the unoptimized structure (plotted solid in left panel) and solid lines the profiles on the S0 -optimized structure. ACS Paragon Plus Environment

9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 32

experimentally measured with non-zero Stokes shift, corresponding to the emission of a ionized chromophore, so the ESPT reaction does definitely take place. Because the QM/MM methods used have been demonstrated to perform accurately in many similar systems, 54–56 we think that a satisfactory explanation must exist able to explain the observed facts, even if it should be more intricate than anticipated a priori. To find this we will critically analyze the hypotheses that have been invoked so far by us and our experimental colleagues.

good agreement with experiment (415 nm), especially taking into account the accuracy of the TDDFT/CAM-B3LYP QM methodology used. Relaxation of the structure in the ground state de-stabilizes the photoactive state by ∼10 kcal mol−1 , which corresponds to a wavelength of 353 nm, a value a little bit too energetic but still within the accuracy reported in the literature for this methodology. 34 In the crystal, we have found that the photoactive state is not the first singlet excited state. A lower energy singlet state with vanishingly small oscillator strength exists, and we have found indications of a state crossing between the latter and the photoactive state in areas of configurational space corresponding to a ionized chromophore and a neutral Asp148. Nevertheless, optimization clearly reduces the energy of the ground electronic state, in this way increasing the energy gap to the photoactive state. In this process, this dark state disappears above the photoactive one. It is not clear whether this state might be playing any role in possible radiationless relaxation channels by which a diminished quantum yield of 0.21 (compared to 0.80 for wtGFP) might arise. Summarizing, both in S0 and in the photoactive electronic state our findings point to the proton being clearly bound to the chromophore, even though a certain anharmonicity is present (especially in the case of relaxed potential energy profile) in S0 , which could cause a certain degree of proton delocalization towards longer O–H distances, but not to the extent of supporting the claim of a delocalized proton in the ground state. This partial picture, especially in what concerns the lack of excited state protontransfer reactivity, flatly contradicts the expectations exposed by Kondo et al. 24 to explain the photochemistry of S65T/H148D. In fact, our results seem to indicate that, in what concerns proton-transfer the system is inert in both ground and photoactive state and fluorescence should occur basically from the A* species with negligible Stokes shift. Our results also fail to appreciate the significant decrease of pKa of the phenolic proton thus far ascribed to the electronic transition to the photoactive state of HBDI. 25 Undeniably, though, fluorescence is

Acidity of the wtGFP Chromophore in the Photoactive State Phenolic protons are known to increase their acidity upon photoexcitation. For instance: In theoretical terms, in our previous studies on Schiff bases we found that photoactive excited states show a pronounced minimum corresponding to the ionized phenolic oxygen. 57,58 For the isolated GFP chromophore (HBDI) Scharnagl and Raupp-Kossmann determined the pKa * value by means of theoretical calculations of thermodynamic type. 25 Concretely, they derived the following equation for the pKa difference between ground and excited state pKa∗ − pKa = ∆pKa∗ + ∆∆pKaS

(3)

where ∆pKa∗ =

1 ∗ [∆EA∗ − ∆EAH ] RT ln 10

(4)

and ∗

∆∆pKaS

1 he−β∆∆EA i0 =− ln −β∆∆E ∗ AH i RT ln 10 he 0

(5)

Equation 4 represents the differences in excitation energies for the ionized and neutral chromophore measured at their corresponding ground state minima, and Equation 5 quantifies the contributions due to the excitation energy differences of the ionized and neutral chromophores in their configurational space areas, and contributes mainly entropic corrections to pKa∗ . The methodology used to compute the terms in Equation 3 was a hybrid quantum-


10

Page 11 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with AM1 for ∆∆pKaS (∆∆pKaS =+2.9) can be used directly in our calculations, we can estimate the pKa∗ value. From Table 1 we get that ∗ ∆EA∗ − ∆EAH =–12.4 kcal mol−1 (AM1 Value: 25 –15.2 kcal mol−1 ), which produces pKa∗ =2.1. This value represents a plunge of more than six pKa units with respect to the value in the ground state. However, the two pKa unit difference with respect to Scharnagl and RauppKossmann’s value is connected essentially to the 9 nm excitation difference for the neutral chromophore, which translates into an energy difference of 2.5 kcal mol−1 in this spectral region. Incidentally, the TDDFT excitation wavelength of neutral HBDI in water we obtain is closer to experiment which means that our pKa∗ estimate might be more accurate. The value of pKa∗ =2.1 for the chromophore in the photoactive state needs to be compared to the pKa of aspartic acid, of 3.7, to predict whether proton-transfer is possible or not. As the chromophore is now about 1.6 pKa units more acidic, reaction should be spontaneous. Then, why do we encounter only endoergic potential energy profiles for this reaction in the photoactive state? As obvious as it might seem, the chromophore (HBDI) differs in two aspects from solution conditions to its state inside the protein. First, there is the structured background due to the different residues in the cavity where the chromophore lies, but also, the chromophore is distorted inside this cavity: the two rings are not coplanar in the crystallographic structure, but they are in the optimized calculations with solvent. It is easy to estimate what is the effect of these two factors (geometrical distortion and protein polarization effects): If the HBDI fragment is taken out of the PDB structure and, without optimizing, the excitation energies of neutral and anionic species in presence of solvent are computed, one gets the effect of the geometry distortion undergone in going from the optimal structure in the solvent to the structure that the protein accomodates without the protein polarization. This is shown in the second-to-last column of Table 1. With these excitation energies and assuming, as before, that ∆∆pKaS of Scharnagl and Raupp-

classical approach, implying the use of the AM1 semiempirical Hamiltonian, introducing solvation effects through a continuum model, because of the large number of calculations involved to compute the averages in Equation 5. 25 The vertical excitation energies in Equations 4 and 5 were computed as the difference between the second and first roots of the CI matrix, using the AM1 Hamiltonian. The use of a semiempirical method to derive the excitation energies in Equations 3-5 is apparently of some concern, as this method was not parameterized to reproduce excited state properties. This is of special relevance here, as excitation energies for the neutral and ionized chromophore are used in Equations 4 and 5, and while in Equation 5 these excitation energies enter as a ratio of averaged quantities (where errors might approximately cancel), in Equation 4 these quantities appear in full, which can have serious effects if systematic errors are expected in the excitation energies. While it is true that CISD excitation energies with the AM1 Hamiltonian have been reported to perform acceptably, 59 we have decided to quantitatively verify how well our TDDFT results perform when looking for the cause of the absence of reactivity in our excited state potential energy profiles. To this end we have computed the excitation wavelengths for the neutral and ionized chromophore in aqueous solution using our TDDFT/CAM-B3LYP methodology and describing the solvent with the PCM model. Results appear in Table 1. The most striking result in Table 1 is the very good agreement between the excitation energies obtained by Scharnagl and Raupp-Kossmann and those obtained with our methodology: excitation of neutral and anionic HBDI are blueshifted by about 40 nm with respect to experimental data with both methodologies This is meaningful not only because it provides a further assessment that AM1 excitation energies give good result, but especially because this tells us flatly that our TDDFT results ought to be compatible with the enormous increase of acidity of the GFP chromophore reported by Scharnagl and Raupp-Kossmann. 25 In fact, assuming that their values computed


11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 32

Table 1: Excitation wavelengths (in nm) of neutral and anionic HBDI in water and of the chromophore in S65T/H148D

d

Protonation Experiment CISD/AM1 TDDFT/CAM-B3LYP Form (water) (water)a (water)b (water)c (protein)d neutral 368e,f ,370g 323 332 388 408 f g e anion 425 ,426 ,428 389 388 436 ∼447 a b Ref 25; S0 -optimized structures in solution, solvent introduced with PCM method; c Unoptimized structure cut from PDB (ID: 2DUF), solvent introduced with PCM method; QM/MM calculation on unoptimized PDB (ID: 2DUF) structure with explicit solvent; e Ref 26; f Ref 60; g Ref 61

Kossmann can be used, we obtain a value of pKa∗ =5.3. The pKa∗ would now be above that of Asp148 and reaction would already not occur. Moreover, if now the chromophore is surrounded by the full protein and solvent as in the crystal, we recover the excitation energies of the neutral form, and a tentative anionic form in which the proton is about the same distance from Asp148 as it was in the case of the neutral chromophore from its oxygen. This is in the last column of Table 1. Comparing these excitation energies with the previous ones reveals what the effect of the protein polarization is. Computing under the same assumptions the value of the excited state pKa∗ yields 6.7. The reason for this change is that the excitation energies of both forms are smaller, but the change is more pronounced for the neutral form, which is smaller by 20 nm while the anionic form by 11 nm (which in the respective regions of the spectrum correspond to 3.6 and 1.6 kcal mol−1 , respectively). At any rate, consideration of geometrical distortion and the effect of the polarization due to the protein affect by reducing the acidity as predicted in water. We think that the trend shown in Table 1 and in the estimated pKa∗ values derived therefrom provides support to the non-reactivity in the photoactive state found in the potential energy profiles computed before. A final assessment on the feasibility of the methodology used so far in this paper can be derived from the excitation energies inside the protein in Table 1. The experimental absorption maximum of S65T/H148D is reported to be 415 nm, which matches extremely well the

408 nm predicted by our method in the protein. Our 447 nm value for the anionic does not compare well with the final photoproduct at 510 nm. 21 However, this would correspond to the final photoproduct emitting when some structural rearrangements might have taken place. In this regard, it is noteworthy to highlight that in their time-resolved fluorescence measurements, Kondo et al. detected 24 transient species with very fast time evolution (within ∼0.1 ps) when excitation happened at 415 nm and fluorescence was monitored at 450 nm (see Figure 1 in Ref 24). In fact, the short time interval between excitation and fluorescence of this transient species precludes any complex structural rearrangement from the original structure (with a neutral chromophore). In other words, this species must differ very little from form A. That a quick change occurs with little geometric rearrangement could be indicative of the motion of just the proton. When we just shift the proton to Asp148 without relaxing, we obtain an excitation wavelength of 447 nm (Table 1) which is again in very good agreement to their measurement. We tentatively conclude that the transfer of a proton without further structural relaxation is well predicted as far as excitation energy is concerned by our calculations, which are compatible with a much-less-acidic than expected chromophore. Unavoidably, the next question in need of an answer is: if the chromophore is not more acidic than Asp148 in the photoactive state, how come the proton gets transferred anyway, and on top, so quickly?


12

Page 13 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


QM/MM Dynamics Simulations

of O(103 ) energies and gradients of the QM part, reaching the nanosecond timescale is very costly. This rules out nowadays effectively QM methods based on ab initio or DFT approaches, as each of these energy+gradient evaluations will amount to at least some minutes with powerful parallel computational infrastructure. But the picture in semiempirical Hamiltonians is not much better, as many of the most powerful programs (like MNDO99 51 ) do not parallelize well. In practice integration of the dynamical equations must proceed on single processors and this limits the maximum simulation time to a few picoseconds per day depending on the size of the system, which makes long simulations in the nanosecond range unfeasible. A possibility to circumvent this problem consists of making good use of the ergodic hypothesis, 46,47 which lies at the basis of statistical mechanics simulations: it is assumed that the time average of a parameter derived from the simulation matches the average of the same parameter over the statistical ensemble. Now: if computing a long trajectory on which to calculate averaged quantities is unattainable, an approach consisting of averaging over many shorter trajectories might do the trick, if ergodicity holds. In a nutshell, the idea is to substitute a very long and costly trajectory with many shorter independent trajectories that belong to the same boundary conditions as the first. As these “independent trajectories” are really independent they can be executed in parallel on different processors, thus achieving effective linear scaling as each processor integrates its own trajectory and adding more processors just means that as many new trajectories are integrated simultaneously. As each of these independent trajectories is much shorter than 1 ns, they can be done in reasonable time even at QM/MM level. The final magnitudes are obtained averaging the value of the observable magnitudes over the full set of individual trajectories. This manual parallelization, or divide-and-conquer approach enables the simulation of nanosecond scale QM/MM dynamics in reasonable amounts of real time. We note in passing that a comparable approach has been recently applied by Koshloff et al. to accumulate over 2 ms worth of

If sheer acidity increase is not behind the reactivity observed in the photoactive state of S65T/H148D, it might be necessary to analyze the dynamical evolution of the protein in physiological conditions. As a matter of fact this is the traditional approach in theoretical studies of enzymes and proteins in general. However, undertaking this approach in the present case is fraught with serious difficulties. S65T/H148D boasts a very short contact between the phenolic oxygen of the chromophore and an acceptor oxygen of Asp148 of only 2.32 ˚ A. It is this closeness which suggested in the first place that maybe the proton is not bound to the chromophore but shared among both groups. Affordable MD methods involve long simulations (usually in excess of 20 ns) based on wellestablished force-fields. However, very short HBs (or alleged LBHBs) do not fit this scheme and have not been parameterized in any forcefields known to us, as they would require a specific new type of interaction to be included, describing a proton effectively shared among two atoms. In cases like these, about the only possibility consists of using QM/MM methods and locate the short hydrogen bond inside the QM region of the system. Needless to say, QM/MM methods are very expensive and cannot commonly be used to perform long simulations. An example of such calculations has been published by some of us on a protein suspected to harbor a LBHB entity (photoactive yellow protein, PYP). 56 Using such QM/MM approach to do molecular dynamics we managed to painstakingly simulate 100 ps worth of trajectory, which sufficed for the purpose it was done. Sometimes the length of a simulation can be critical if a certain arrangement of the protein needs to take place before a certain process is triggered. In view of this we have chosen to take a different approach for our study of S65T/H148D that should provide a simulation of good length with current computational infrastructure. The main problem in considering QM/MM dynamics simulations is that the computational effort rises with the length of the simulation: as every picosecond requires the evaluation


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

classical MD simulations on G-protein-coupled receptors, and used to sample the occurrence of rare events in these complicated systems by means of Markov state models. 62 In practice, ergodicity is ensured by sampling configurational space points connected by a precursor trajectory, and then, spawning one trajectory from each once velocities are randomly reinitialized (and a short equilibration done). See Figure S1 for a scheme detailing how this approach works. Following the procedure described in the Computational Methods section, a “production set” of 1 ns OM3 QM/MM dynamics simulation was accumulated from which statistical averages can be computed, which by virtue of the ergodic hypothesis should match the corresponding averages computed on a single 1 ns trajectory. The only caveat is that no chronological information can be obtained in this way, except for each of the individual member trajectories. Before drawing conclusions of physical nature from this long simulation, it is necessary to verify how well does it describe the key interaction between the chromophore and Asp148. Figure 4 shows the distribution of O–O distances of this tight contact. This graph corresponds to a Gaussian-like distribution centered at 2.35 ˚ A, which corresponds to a well formed, strong hydrogen bond, which is stable over time (as no structures have been found with dO−O > 2.6 ˚ A). Even though this value is in excellent agreement with experimental measurements from Xray diffraction data, a note of caution is necessary here: the QM method of choice for our OM3 QM/MM molecular dynamics simulation, despite ranking among the best methods currently available in the field of modern semiempirical Hamiltonians, is known to somewhat overestimate the strength of hydrogen bond interactions. 48 This makes us think that the actual situation might differ from the result of our simulation in this aspect, such that the Asp148-chromophore interaction is likely a little bit weaker than OM3 predicts, which should translate itself in a longer O(Asp148)–O(Cro) distance than found by us. At any rate, OM3 provides a better description than all semiem-

Page 14 of 32

pirical methods known to us for this kind of interaction. Throughout the simulation, Asp148 stays very close to the chromophore, but also its relative orientation is quite stable over time. Figure 4 also shows the distribution of values of the dihedral angle formed by the carboxylic group of Asp148 with respect to the phenolic ring of the chromophore. This angle shows also a unimodal Gaussian-like distribution of values with an average of 51.6◦ . The distribution is indicative also of a fairly stable relative orientation of both fragments, or in other words, that the “primary” environment around the proton-transfer characters is quite stable over a respectable simulation time, and because this structure is not unlike the one reported in the X-ray diffraction structure (which shows dO−O = 2.32 ˚ A, ∠Asp148-Cro=58.6◦ ) we conclude that apparently the X-ray structure, as concerns the proton-transfer characters, holds well under OM3 QM/MM molecular dynamics at room temperature. We have obtained clear, bell-shaped unimodal distributions for the two parameters analyzed, a fact that can be taken as an indication of both the structural stability of the closest environment of the chromophore in constant T simulations on one hand, and on the other, that the overall design based on the ergodic hypothesis is valid as the structures of all simulations fit in a unique uniform distribution. At this point we would like to comment on the significance of the optimized structure found at the beginning of the Results section. Proteins in solution at room temperature are dynamical entities and over time explore a wide area of configurational space, as should be evident in this case from Figure 4. The optimized structure, although obtained with a different methodology, represents at best a specific configuration of particularly low energy (not even the lowest energy configuration). As such, even if this particular structure should exist, it will be visited (found) an extremely tiny fraction of time due to the vanishing statistical weight of this configuration compared to the rest. It is debatable whether reliable conclusions can be obtained from an optimized structure of a protein in view of this, being much more correct to


14


1200

1200

900

900

600

600

N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

N

Page 15 of 32

300

0 2.1

300

2.2

2.3 2.4 dCro−Asp148 / Å

2.5

2.6

0 0

15

30 45 60 75 ∠ Cro−Asp148 / degrees

90

105

Figure 4: Distribution of structures of the Cro-Asp148 unit as derived from the 1.0 ns OM3 QM/MM dynamics simulation. Left panel: Histogram of the O–O distances between donor (Cro) and acceptor (Asp148). Right panel: orientation of Asp148 with respect to the phenolic ring of the chromophore. The figures are based on the 1.0 ns OM3 QM/MM dynamics simulation. and the same number of excited state TDDFT QM/MM calculations. On these potential energy profiles we have carried out topographical analyses aimed at determining the individual potential energy landscape encountered by the protein along its dynamical traverse in configurational space at room temperature. Thus, we have computed the ergicity of the protontransfer process in the ground state (∆E) as the energy of the ionized chromophore structure of the profile minus the energy of the neutral chromophore in the same profile. The ergicity in the photoactive state (∆E ∗ ) is derived in an analogous way. It is worth recalling (Figure 3) that the situation we encountered in our study of the crystal form (as well as after optimization of the ground state) would correspond, in these terms, to large values of ∆E (ca. 10 kcal mol−1 ) and ∆E ∗ (somewhat smaller), corresponding to a non-reactive process in ground and excited state, with very similar λabs and λem (no Stokes’ shift).

contemplate the plethora of them that arise in a room temperature simulation. After the full ensemble of trajectories is shown to provide a reasonable description of the closest region of the proton-transfer between Asp148 and the chromophore in constant T conditions, it is time to determine whether the dynamical traverse of the protein gives rise to instantaneous configurations that differ significantly from the situation we encountered in the crystal and optimized structures. Because the energetics of the ground, but especially the photoactive, states are critical to determine whether the proton-transfer can take place, we have opted to use high level QM/MM calculations at DFT (for the ground state) and TDDFT (for excited states) levels, and including polarization due to the environment via the QM/MM approach (henceforth DFT QM/MM and TDDFT QM/MM) taking as valid the OM3 QM/MM dynamics simulation. Thus, we have selected 500 structures from the QM/MM dynamics simulation evenly spaced (2 ps apart) to avoid correlation, and for each we have computed a 1-dimensional nonrelaxed potential energy profile corresponding to the proton-transfer, built by changing the position of the proton from donor to acceptor. Each potential energy profile consisted of 10 individual points. This represented a total of 5,000 single-point ground state DFT QM/MM,

Proton-Transfer Energetics Figure 5 displays the results derived from the DFT QM/MM and TDDFT QM/MM potential energy landscapes in terms of the energetics of the ground and photoactive states. Analysis of these data reveals a drastically different situation than in the crystal structure. The


15

10

0

5

∆E* / kcal mol

−5

0

5

−5

0 −10

50

100

∆E* / kcal mol−1

20 15

15

−1

20

−10 −5

0

5

10

15

20

15

20

∆E / kcal mol−1

N 100 N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

10


50 0 −5

0

5

10

∆E / kcal mol

−1

Figure 5: Analysis of the energetics of the proton-transfer in ground and photoactive electronic states, as derived from the individual potential energy profiles of proton-transfer from chromophore to Asp148, computed at DFT QM/MM and TDDFT QM/MM levels. The abscissa represents the exo- or endoergicity of the process in the ground state for a given profile, and the ordinate does the same for the photoactive state. The values are computed as the energy difference of the Asp148-H· · · Cro structure minus the energy of Asp148· · · H-Cro. Red dots correspond to profiles defined as “not reactive (∆E ∗ > 0) ” while blue dots to profiles defined as “reactive (∆E ∗ < 0)”. Green dots correspond to profiles where the proton is found bound to Asp148 directly, in the ground state. The horizontal inset shows the histogram of ∆E and the vertical inset that of ∆E ∗ . Diagonal line represents the points where ergicity is the same in ground and photoactive states.

lower horizontal histogram in Figure 5 represents the distribution of the snapshots according to the energetics of the proton-transfer in S0 , ∆E. A very large number of snapshots lead to quite isoenergetic potential energy processes in the ground state, with about 50% of the structures in the range ∆E < 3 kcal mol−1 . In this way, consideration of the dynamics drastically changes our concept of the energetics of the ground state from a proton clearly bound to the chromophore in the crystal structure to a proton that might be substantially delocalized between donor and acceptor atoms. This is made all the more eloquent by the fact that a small but significant fraction of the potential energy profiles analyzed have negative ∆E values (concretely, 5.7%). In terms of the position of the absolute minima in S0 , this implies that the proton is already transferred even in the ground state, or in other words, that the system is found in 5.7% of the samples harboring in the ground state a neutral Asp148 and a ionized chromophore. These results are indicative of a substantial change in the state of affairs from the crystal structure, all the more remarkable because the raw DFT QM/MM methodology is exactly the same in both cases, so that the changes, if any, can only be due to the fact that the dynamics brings the system across areas of configurational space significantly different from those represented by a potential energy minimum or the crystal structure. Overall, the picture that the dynamical study renders of the protein in solution corresponds to a much flatter (isoergic) profile for the ground state, with about 40% of the potential energy profiles displaying about ±1.5 kcal mol−1 energy difference between reactant and product of the proton transfer. With energy differences of this small entity it is conceivable that in a substantial part of them the ZPE lies at or above the possible potential energy barrier of the proton-transfer. In practice this would describe a situation where the proton is effectively shared by the donor and the acceptor atoms with like strengths. It is noteworthy that such a possibility was already suggested by StonerMa et al. as a possible explanation to the lack of change in the ionization state of a carboxylate


16

Page 16 of 32

Page 17 of 32

residue concomitant to the formation of the I* species within time resolution of 400 fs. 23 This suggested that if the proton were already shared between donor and acceptor in S0 , no possibility of observing the protonation of Asp148 should exist. Our results concerning the flatness of the potential energy profile in S0 give theoretical support to this suggestion. In the photoactive state the situation is also markedly different than in the crystal. The results we showed before when discussing the pKa∗ make us expect a certain increase of the acidity of the chromophore, and thus, it is legitimate to expect that, on average, ∆E ∗ of a snapshot should be less than ∆E for the same profile. Figure 5 shows that almost all profiles comply with this expectation (as most dots are under the isoergicity line). Hence we confirm a general increase of the acidity of the chromophore in the protein in the photoactive state. The values of ∆E ∗ —which in the end will determine whether proton-transfer occurs favorably or not— are clustered around zero, with ∼70% of the snapshots having ∆E ∗ < 1.5 kcal mol−1 . This means that upon photoexcitation the system will find itself most of the time in a situation where proton-transfer can occur. The natural question coming to mind is: Why is it that the dynamics changes the situation so much as to make the proton-transfer potential energy profile so much flatter in ground and photoactive state? The answer to this question is nevertheless connected also to the increase of viability of the proton-transfer in the photoactive state. We have looked in the different configurations that arise along the dynamics of the protein, and found two main factors that contribute to this change in the potential energy landscapes. (a) Tyr145 is a inconspicuous residue situated close to Asp148. In the crystal structure it can be seen with the residue tail projected to the inside of the β-barrel structure and not particularly close to the main characters of the protontransfer. Its phenolic oxygen is 4.3 ˚ A apart from the chromophore’s phenolic oxygen. However, during the OM3 QM/MM dynamics simulation this ligand is seen to approach the chromophore and remain for a while in its vicin-

1000

800

600 N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


400

200

0 2.5

3.0

3.5

4.0 4.5 5.0 dO(Cro)−O(Tyr147) / Å

5.5

6.0

Figure 6: Histogram of the O–O distances between chromophore and residue Tyr145 resulting from the 1.0 ns OM3 QM/MM dynamics simulation. ity. Figure 6 shows the distribution of O–O distances between O(Cro) and O(Tyr145) over the full OM3 QM/MM dynamics. A bimodal distribution is detected with two peaks, one at dO(Cro)−O(Tyr145) ∼ 4.4 ˚ A (containing about 65% of the structures, and which coincides with the distance measured in the crystal structure) in which the residue tail is oriented towards the inside of the β-barrel structure, and another, containing about ∼ 35% of the structures, which corresponds to Tyr145 being close to the chromophore’s phenolic oxygen (peaking at dO(Cro)−O(Tyr145) ∼ 2.8 ˚ A). This second conformation is akin to a hydrogen bond interaction between Tyr145 and the chromophore (See Figure S2 in the Supporting Information for a picture of this kind of alternance in structure). What effect can it have that this residue establishes over time a hydrogen bond interaction to the chromophore? This interaction is established with a formally neutral chromophore, which is stabilized by this interaction. However, the presence of this hydrogen bond to the chromophore will have a larger stabilizing effect on the ionized chromophore. Thus, this interaction might increase the tendency of the chromophore to release its proton (increase its acidity) as it will preferentially stabilize its ionization products. (b) We have also observed that solvent water


17


the distance from the center, and Ni (r) is the number of atoms of type i that are found at distance r from the center. A value of gi (r) = 1 denotes a local density equal to that of the bulk for atoms of type i, and values of gi (r) larger (smaller) than 1 indicate a local density that is higher (lower) than in the bulk. If we take for i the oxygen atom of water molecules, it is possible to spot and quantify the existence of solvation layers around specific centers. Figure 7 shows gO(Wat) (r), that is, the radial distribution function of the oxygen atom of water molecules for the OM3 QM/MM dynamics simulation, taking as centers the oxygen atom of the chromophore, both oxygen atoms of Asp148, and the oxygen of Tyr145. As has been shown before, Asp148 usually maintains a locked orientation with respect to the chromophore (see Figure 4). This means that one of the carboxylic oxygens has been found to point normally towards the outside of the protein (which we denote by O1 ) and the other towards the inside of the β-barrel cavity (which we call O2 ). Figure 7 shows that both oxygens of Asp148 have a well-peaked radial distribution function of water molecules with a pronounced maximum at r ∼ 2.65 ˚ A. Afterwards gO(Wat) (r) drops to almost zero and subsequently a second, less pronounced maximum is appreciated at 4.5 ˚ A for O1 , which indicates the presence of a stable secondary solvation shell. O2 has also this maximum in gO(Wat) (r) a little later (4.8 ˚ A). These features describe a well solvated Asp148, representing the arrangement of water molecules situated on the outside of the protein. The chromophore also has water molecules close, as gO(Wat) (r) peaks at 2.85 ˚ A, but it is less pronounced, indicating that less water molecules are around in interaction with it. Last in this list of interactions is Tyr145, which shows a very low value over the entire r domain for gO(Wat) (r) and that indicates the absence of strong interacting water molecules in its vicinity (peaking at about 2.9 ˚ A). All of this makes good sense, if one remembers that Asp148 is close to the surface of the protein (and because of its locked arrangement, one of its oxygens is closer to the surface than the other systematically), the chro-

1.5 X=O(Cro) X=O1(Asp148) X=O2(Asp148) X=O(Tyr145) gO(X)−O(Water)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1.0

0.5

0.0 2.5

3.0

3.5 4.0 dO(X)−O(Water) / Å

4.5

5.0

Figure 7: Radial distribution functions of the oxygen atom of waters in the 1.0 ns OM3 QM/MM simulation, where the bulk density used is taken to be 1 g mL−1 . The different curves correspond to the radial distribution function of water oxygens from the specific centers indicated. molecules are likely involved in the energetics of the proton-transfer process. In fact, Asp148 is situated, already in the crystal structure, close to the surface of the β-barrel, facing the chromophore but also in contact to the surface of the protein and consequently in a potential situation to interact with water molecules. The relevance of the interaction with water molecules dawned on us in a typical case of serendipity (see Supporting Information) in the initial stages of this research. Analysis of the contribution of water molecules is troublesome because they are not covalently bonded to the protein and can move around, which makes it complicated to trace what water molecule might be involved in interactions with the chromophore or Asp148 along the dynamical situation. A standard method to study this makes use of the radial distribution function g(r), which is defined as the ratio of the local density of a given atom type to the bulk density of the same atom type at a given distance from a specific center: 47 gi (r) =

1 dNi (r) 4πr2 ρi dr

Page 18 of 32

(6)

where ρi is the bulk density of atom type i, r is


18

Page 19 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ton transfer in the ground state.

mophore is buried in the protein but is very close to Asp148, while Tyr145 is deeply buried in the β-barrel structure and thus in hydrophobic environment. It is possible to compute the average number of water molecules that are interacting directly with each of these atoms by isolating dNi (r) in Equation 6 and integrating gO(Wat) (r) up to the first minimum: Z rmin NO(Wat) = gO(Wat) (r)ρO(Wat) 4πr2 dr (7)

Photoactive State Energy Profiles We have analyzed again the distribution of TDDFT QM/MM potential energy profiles in Figure 5 to see how they fare as far as the photoactive state is concerned. As mentioned before, 5.7% of the profiles described a system already with a ionized chromophore in S0 . The rest are described as follows: 64.8% of the profiles have a positive (albeit small) value of ∆E ∗ and thus describe a situation where the excited chromophore is more stable when it is neutral. Nevertheless, the vertical histogram in Figure 5 indicates a very large proportion of profiles of ∆E ∗ < 2 kcal mol−1 , so the endothermicity is very small for the proton-transfer process. Already 29.5% of the TDDFT QM/MM potential energy profiles have ∆E ∗ < 0 and correspond to a ionized chromophore as the most stable species in the photoactive state. This percentile is suggestive, especially when compared to the measured quantum yield of 0.21. To show explicit examples, Figure 8 displays three selected potential energy profiles at DFT QM/MM level for the ground state and TDDFT QM/MM level for the photoactive state, corresponding to a typical ∆E ∗ > 0 (“not reactive”), ∆E ∗ < 0 (“reactive”) and finally a case where ∆E ∗ < 0 and fluorescence is observed at large Stokes shift. It is striking that in the first two cases in Figure 8 (left and center panels) the potential energy profiles are very unimodal, with a single minimum that puts the proton closer to the chromophore than to Asp148 in the ground electronic state, and without a secondary less stable minimum describing a ionized chromophore. This describes that even for non-reactive profiles there is a single minimum that describes a proton effectively shared between Asp148 and the chromophore. It is possible, actually, to see this quantitatively if the ZPE is computed for these profiles. Solving the time-independent Schrödinger equation as described in the Computational Methods, we encounter ZPE values of 2.0 kcal mol−1 for the typical “non-reactive” profile, 2.0 kcal mol−1

0

This renders the following figures: O1 (Asp148) 1.26 waters, O2 (Asp148) 1.10 waters, O(Cro) 0.92 waters and O(Tyr145) 0.37 waters. An average of more than 2 water molecules are found stabilizing Asp148, which is reasonable as it is normally anionic and interacts strongly with them. One of these molecules is inside the protein, while the rest are outside it and affect the external side of Asp148. The chromophore is interacting stably with one water molecule. This is in fact the water molecule that partakes in the old well-known 3-link proton shuttle connecting the chromophore, this water molecule, Ser205 and Glu222. Finally, Tyr145 only interacts with 0.37 water molecules, or rather, with one water molecule about 37% of the time, which matches very well with the 35% residence time of Tyr145 in the neighborhood of the chromophore (the first peak in Figure 6). This means that the only water molecule that Tyr145 interacts with (and then, feebly) is the same that is captive between chromophore and Ser205. Summarizing, water molecules contribute to stabilize the different modes of interaction of Asp148 and the chromophore, and Tyr145 will stabilize better a ionized chromophore and is present about 35% of the time. Besides, while in the static studies Asp148 is found to interact with three water molecules, the snapshots of the OM3 QM/MM dynamics simulations analyzed show that only two water molecules are found at any time interacting with Asp148. We conclude that the finely balanced number of water molecules on one side, and the presence of residue Tyr145 on the other, determine the flatness of the potential energy profiles for the pro-


19

Page 20 of 32 100

95

95

95

90

90

90

85

85

85

80

80

80

E / kcal mol

E / kcal mol

75 20

−1

100

−1

100

E / kcal mol

75 20

75 20

15

15

15

10

10

10

5

5

5

0

0 0.8

1.0

1.2 1.4 dCro−H / Å

1.6

0 0.8

1.0

1.2 1.4 dCro−H / Å

1.6

0.8

1.0

1.2 1.4 dCro−H / Å

1.6

Figure 8: Potential energy profiles obtained by DFT QM/MM (ground state, in red) and TDDFT QM/MM (photoactive state, in green) methods on three selected snapshots of the 1.0 ns OM3 QM/MM dynamics simulation. In the left panel, a situation typical of the “not reactive” profiles in topographical terms (∆E ∗ > 0). In the central panel, a common situation for profiles describing a “reactive” profile (∆E ∗ < 0). In the right panel, a profile that is “reactive” and gives a large Stokes shift. The donor-acceptor O–O distances are 2.44 ˚ A, 2.38 ˚ A and 2.53 ˚ A, respectively. In black, the ground vibrational levels of ground and photoactive electronic states (solid lines) and probability density functions |ψ0 |2 (dotted lines, not to scale) of the protium isotopologue. for the typical “reactive” profile and 2.4 kcal mol−1 for the “reactive” profile leading to large Stokes shift. Figure 8 shows these levels, as well as the vibrational probability density function |ψv=0 |2 . As these probability densities describe where is it most likely to find the proton in a position measurement, it should be obvious that the proton is widely spread between donor and acceptor even when the most stable minimum corresponds to a neutral chromophore in topographical terms. All of this is in complete agreement with the hypothesis formulated by Kondo et al.: 24 the proton will be found about anywhere between the phenolic oxygen of the chromophore and the oxygen of Asp148, and can be consequently described without ambiguity as “shared” between both. Thus, the lack of time evolution of the carboxylic signature in the time-resolved infrared experiments can be explained —now with detailed microscopic description— by the fact that the ground vibrational level is high enough as to lead to a delocalized proton, similarly to what would be expected in the case of an LBHB. To try to discover whether there is any struc-

50 All Reactive Not Reactive Product in S0

40

30 N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

−1


20

10

0 2.2

2.3 2.4 dO(Cro)−O(Asp148) / Å

2.5

Figure 9: Distribution of O–O distances between chromophore and Asp148 on the snapshots studied at TDDFT QM/MM and obtained from the 1.0 ns OM3 QM/MM dynamics simulation. The distribution functions shown correspond to those snapshots that were found to be reactive and not reactive in the photoactive state and those found already in the product region in S0 (ionized chromophore and Asp148 in neutral form).


20

Page 21 of 32

tural factor that favors negative ∆E ∗ values, we have collected a large set of interatomic distances between atoms close to the Cro-Asp148 unit, including solvent waters, for all profiles and statistically compared the profiles with negative ∆E ∗ and those with positive ∆E ∗ . A Student-t test with 95% confidence revealed the two sets to be the same. Then, we conclude that no individual interactions differ between both situations and must thus fall back to the conclusion that the system seems to be callibrated to yield essentially isoergic profiles. The deviation to positive or negative ∆E ∗ values are small and explainable by invoking the fluctuations encountered in the dynamics, which are blatantly absent from a pure static picture, like the one provided by the X-ray structure. Figure 9 shows how “reactive” and “not reactive” profiles correlate with the donor-acceptor distance. Although both types of situations occur throughout all the spectrum of O–O distances, it is possible to notice that the distribution of “reactive” profiles peaks at O–O distance values slightly shorter (2.33 ˚ A) than those of “not reactive” ones (2.37 ˚ A). Again, as OM3 (our semiempirical method of choice for the QM/MM dynamics) seems to overestimate the strength of hydrogen bond interactions, the actual situation might differ from our simulation spreading the O–O values to longer values, which would imply that the proportion of endoergic potential energy profiles would increase. However, in interpreting this extrapolation we must bear in mind that the histograms in Figure 9 are based on TDDFT QM/MM results, so it is very reasonable to expect that a large fraction of systems will be cluttered around ∆E ∗ ∼ 0 and thus, no significant deviation of the previous conclusions should be expected.

function of the O(Cro)–O(Asp148) distance. It is not surprising that for most of the profiles the Stokes shift turns out to be small (ca. 5 nm) in view of the kind of potential energy profiles we have encountered. Nevertheless, a small number of profiles (seen as the offshoot that projects to the top right in Figure 10) reach relatively large values in the range of 10-20 nm. These come from a small set of reactive profiles which resemble the ones in the right panel of Figure 8, and have a large Stokes shift because they display a shoulder of potential energy in S0 corresponding to a frustrated anionic chromophore minimum, and another in the photoactive electronic state corresponding then to a frustrated neutral chromophore. Profiles like this one are encountered especially when the O(Cro)–O(Asp148) distance is long (the case depicted in Figure 8 corresponds to 2.53 ˚ A). 30

20 ∆λ / nm

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


10

0 2.1

2.2

2.3 2.4 dO(Cro)−O(Asp148) / Å

2.5

2.6

Figure 10: Representation of the predicted Stokes shift for the snapshots of the 1.0 ns OM3 QM/MM dynamics simulation that lead to a reactive profile, as defined in the main text. This comprises about 29% of the snapshots analyzed. It is very suggestive that Kondo et al., in their time-resolved fluorescence experiments, detected instantaneous rise of fluorescence (with a time resolution of ∼ 50 fs) on the blue edge of the I* band (450 nm) after excitation at 415 nm. 24 In other words, considering an Stokes shift of 450 - 415 = 35 nm, the chemical process that causes it takes place within 50 fs or less. Our calculations tell that if the proton just moves and everything else remains frozen in

Predicted Stokes Shift The above results allow us to determine, for those profiles with ∆E ∗ < 0 in which we expect spontaneous proton transfer in the photoactive state, what the wavelength shift would be if the system emits fluorescence right when the product (ionized excited chromophore) is formed. This is represented in Figure 10 as a


21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

time, profiles with a large O(Cro)–O(Asp148) distance in our OM3 QM/MM dynamics simulations lead to Stokes shifts of about 20 nm (the largest we found being 24 nm). This again fits really well with the tentative picture that Kondo et al. made of a very fast motion of the proton in the femtosecond time range, that ought to correspond to the ballistic motion of the proton from chromophore to Asp148. They detected different time evolutions when monitoring fluorescence on a wide range of wavelengths (450-520 nm) that should correspond to the relaxation of the chromophore and cavity after the sudden traslocation of the proton, while the ionized chromophore is still in the excited photoactive state. This specific aspect escapes the present simulation and work is in progress along this line in our laboratory. Two new questions arise: The first is that we encounter very few configurations on our OM3 QM/MM dynamics describing this large Stokes shift. We must remember that OM3 overestimates the strength of hydrogen bond interactions. The real interaction is surely weaker and we should expect longer Cro-Asp148 distances on average. This would increase the representation of potential energy profiles like the one in the right panel of Figure 8, which has a large Stokes shift. As is easy to guess, this would flesh out the offshoot in Figure 10 at long O–O distances and large Stokes shift to the detriment of the region with short O–O distances and small ∆λ.

Page 22 of 32

toabsorption obeys the Franck-Condon principle, according to which the electron cloud rearranges much faster than the nuclei. This is simulated by changing the electronic state whilst keeping the nuclear configuration intact, which in terms of the simulation to be done means that the vibrational Hamiltonian is suddenly changed to the one of the photoactive state. As the vibrational ground state is not a stationary state of the photoactive state Hamiltonian, this state will evolve with time under the control of the time-dependent Schrödinger equation. We have solved this wave packet propagation using the procedure described in the Computational Methods section, for a total time of 150 fs. To monitor the state of the system at a given time t we have computed the survival probability, which represents the fraction of the proton probability density that is found in the region of reactants (R): Z Θ(t) = |ψ(q, t)|2 dq (8) q∈R

where R covers the values of the coordinate(s) of the proton that describe it as being closer to the donor atom (O(Cro)) than to the acceptor (O(Asp148)). The survival probability is a function with values between 0 and 1, and it can be interpreted as the probability that, upon measurement of the position of the proton, it be found in the reactant (neutral chromophore) region. The results are shown in Figure 11. The results we obtain indicate recurrent motion of the wave packet. The behavior is appropriately described as a periodic bouncing motion from reactant to product, which is accurately described as a Rabi oscillation, 63 between two limiting cases (proton bound to chromophore, on one hand, and on the other, proton bound to Asp148). The period of this Rabi oscillation ranges from ∼20 fs to about ∼32 fs. The time necessary to achieve the minimal value of the survival probability is half this value (from ∼10 fs to ∼16 fs, and this value should be taken as a theoretical estimate of the proton (H) transfer time in a situation where the environment (that is, everything except the proton that is transferring) remains fixed.

Timescale of Proton Motion The second question concerns the timescale of the proton motion. We have performed quantum dynamics calculations to describe the proton motion. This we have done for the three representative cases depicted in Figure 8, in a way that tries to mimick the photoabsorption process, as set forth in the Computational Methods. This is done by first finding the ground vibrational eigensate for S0 for the potential energy profile (this is the state of the system the instant before the photon is absorbed). Figure 8 displays the corresponding probability densities for each profile. Second: Pho-


22

Page 23 of 32

almost complete) of the proton within times under 20 fs. The estimate of the transfer time fits well with the available experimental estimates found by Kondo et al. 24 Hereafter, the proton “resonates” between donor and acceptor atom. This simple periodic oscillation that is predicted by our QD calculations is not real, though: Relaxation of the structure of the supporting groups should ensue, to accomodate the charge shift from Asp148 to chromophore that the protontransfer carries with it. It is to be expected that this relaxation of the scaffold that supports the proton-transfer groups will stabilize the minimum of the photoactive state with the ionized chromophore preferentially, leading to a substantial red shift of the fluorescence wavelength. The study of this reaction coordinate in the excited state is an absolutely non-trivial task and is currently in progress in our laboratory. We have also solved the QD for the deuterium (D) transfer. The graphical representation of the survival probability is depicted in Figure 11 alongside the protium transfer one. Deuterium transfers slower in all cases when compared to protium transfer on the same profile, yielding values between 16 fs and 22 fs. When compared to transfer of proton on the same profile, the transfer times obey τD /τH & 1.4, in good agreement with ballistic (in the sense of barrierless) motion. Nevertheless, deuterium has been found to transfer under 35 fs, and so well under the experimental resolution time. Thus, our calculations predict a slower motion of the deuteron that will be undetectable with current time resolution available.

Θ(t)

1.0 0.5

Θ(t)

0.0 1.0 0.5 0.0 1.0 Θ(t)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.5 0.0 0

20

40

60

80

100

120

140

t / fs

Figure 11: Time evolution of the survival probability (Equation 8) for three typical potential energy profiles corresponding to a non-reactive situation (red), a typical reactive profile (blue) and a quite favorable profile for excited-state proton transfer (green). Solid lines correspond to the dynamics of protium (H) transfer and dotted lines to the dynamics of deuterium (D) transfer. It is noteworthy that proton-transfer occurs even for profiles with ∆E ∗ > 0 (“not reactive”, where the motion of the proton is unfavorable energetically) even if to a smaller extent: only 60% transfer in the “not reactive” case, to be compared to 80% transfer for the “reactive” case or even 95% for the case with large Stokes shift. The fact of why transfer does occur even in situations where it is energetically uphill is just another eloquent proof of the unique importance of dynamical simulations over discussions based exclusively on topographical assessments (in other words, whether the process is energetically favorable or not). The profiles, in general and as seen in Figure 5, are better defined as very shallow in almost all cases, and the simple energy content of the wave packet after landing in the photoactive state is large enough to ellicit effective transfer. At any rate, even endoergic profiles show appreciable transfer. This is so because, apparently, the system is “designed” to attain configurations that do not deviate substantially from an isoergic (∆E ∗ ∼ 0) situation. It will in all cases show appreciable transfer (in some cases,

CONCLUSIONS This work strives to provide an atomic-scale description of the mechanism operating within the GFP variant S65T/H148D using state-of-theart methodologies and techniques in theoretical chemistry, including non-equilibrium quantum dynamical wave packet propagation to obtain theoretical estimates of the time of transfer. GFP variant S65T/H148D was synthesized


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

in an attempt to modulate the pKa value of the chromophore of its original system, S65T, by substitution of His148 by a proximal negative charge provided by an aspartate residue. This protein was surprising in that it recovered the bright green fluorescence that is lost in S65T, even though it was established that proton-transfer occured now between the chromophore and the new Asp148 group. X-ray diffraction studies found a very short contact between the chromophore and this aspartate residue and not along the conventional proton shuttle route operative in wtGFP. Time resolved spectroscopic studies established that (1) proton (and deuteron) transfer occured under 50 fs with no measurable kinetic isotope effect and (2) within experimental resolution no time evolution could be detected of the IR signature of the carboxylate residue. This set of experimental data found a tentative explanation in the form of a hypothesis that relied heavily on the drastic plunge of pKa of the chromophore predicted for the chromophore in solution, 25 which suggested a barrierless transfer in the photoactive state. Besides, the tight contact between donor and acceptor found in the crystal structure, combined with the impressive velocity of the transfer, was used to propose a very delocalized proton in the ground electronic state, which would explain the lack of measurement of time evolution of the IR signature of a carboxylate group (as the proton would be already half-transferred in the ground electronic state). This work, using QM/MM computations for ground and photoactive state on the crystal structure, has revealed that the proton is clearly bound to the chromophore, and cannot be described as delocalized in the conventional sense of the term. Optimization of this structure does not change this conclusion. In any case, reaction in the photoactive state is not predicted, in clear contradiction to the experimental observations. It is remarkable that appropriate use of theoretical methods manages to confirm the hypotheses that were used by experimentalists to explain this ultrafast process. However, we find beyond reasonable doubt that the so-far

Page 24 of 32

accepted driving force, namely the sudden increase of acidity of the GFP chromophore in the photoactive state, is not enough to trigger the proton-transfer. The chromophore is very sensitive to distortions of its geometry (especially to out-of-plane torsions 25 ) and the presence of the highly structured protein environment. Even though there is a clear increase of acidity, it does not suffice to overcome the acidity of aspartic acid. It is not until the full-blown complexity of the protein, including solvent water molecules, is brought onto the stage that the otherwise unfavorable energetic conditions change drastically into others in which reaction is possible and can occur. It has been noted elsewhere 64 that biological macromolecules might show very different properties in the crystal state and in solution, under biological conditions, and we think that this is another suggestive example of this norm. Sophisticated and costly QM/MM dynamics simulations of the protein in solution reveal the protein to be a dynamic entity involved in a delicate interplay with the solvent. The main characters in the chromophore cavity interact among themselves but also other groups that have the potential to establish hydrogen bond interactions. Tyr145 is a residue close to Asp148 that is oriented towards the inside of the protein and that has been found to stabilize the chromophore (and more importantly, to stabilize better the ionized chromophore) even though in the crystal structure it is found oriented in a different way. Water molecules manage to interact in right the proper proportion with Asp148 and the chromophore. In all, these interactions manage to bring the otherwise unfavorable potential energy profiles for the proton-transfer into a situation of essentially isoergicity both in ground and photoactive states when studied under precise DFT QM/MM and TDDFT QM/MM computational techniques. Under these conditions the proton-transfer is much more favorable than in the crystal structure (or its optimized counterpart). Moreover, this isoergicity situation occurs without occurrence of a potential energy barrier in neither the ground nor the photoactive state. Thus, the vibrational ground state


24

Page 25 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


work and shown to lie within the expectations suggested by experimental data. 65 Moreover, Di Donato et al. studied with timeresolved IR spectroscopy the behavior of wtGFP in the photactive state, tracking the time evolution of the IR signature of carboxylic groups (that is, tracking in the end the change of protonation state of Glu222). 66 In their detailed study they reached the conclusion that the process is initiated in the photoactive state by the transfer of the last proton (Ser205 → Glu222), which gives rise to some intermediate with characteristics of a highly delocalized proton system (concretely, the protons originaly anchored to the chromophore and Wat25). In their view, the proton transfer is initiated by the transfer of the last proton, but reaches a situation of isoergicity in what corresponds to the transfer of the remaining protons, such that the reaction comes to an end when the environment drags it thereto. In all, it looks like S65T/H148D (and because of the above, likely also wtGFP) is involved in a delicate dynamical behavior that is consubstantial to its structure and the nature of the solvent in which it performs its activity, and that carries with it the arising of situations of almost degeneracy for all relevant steps of the photocycle that involve proton migration. This might then occur at high speed and with no intervening barriers (in the case of wtGFP, these barriers are there but must be small). GFP S65T/H148D is special in that (1) it has a single proton-transfer event and (2) it occurs so fast that no KIE for the proton transfer could be observed. It is in this sense of absolute simplicity that we think that GFP S65T/H148D can be considered an archetype for this kind of processes. We hope that our findings in this work can be useful to understand how fluorescent proteins operate at the molecular level. During the peer-review process of this work one of the reviewers brought to our attention the very recent theoretical study by Zhang et al. on the same GFP mutant. 67 In this work the authors perform a Molecular Dynamics simulation to select one candidate structure with adequate hydrogen bonding layout, and the authors use this structure to carry out

for the proton motion describes an essentially delocalized and effectively shared proton in the ground state, thus confirming the hypothesis postulated by Stoner-Ma et al. 23 If representative potential energy profiles derived from the simulation of S65T/H148D GFP are used to compute proton-transfer times upon photoexcitation, we encounter very short transfer times that are well under the reported experimental precision (50 fs) and without measurable kinetic isotope effect. While GFP S65T/H148D ranks among the fastest fluorescent proteins derived from wtGFP, it is also true that other members of the family (especially wtGFP) also operate really fast proton-shuttle mechanisms in the photoactive state. It is remarkable that this is achieved by modulating the energetics of the system in such a way that unfavorable conditions are rendered almost isoergic. This is not a exclusive feature of GFP S65T/H148D: the case of wtGFP is especially revealing. It is well known that the A and B bands in wtGFP correspond to structurally relaxed proteins adapted to a neutral and a ionized chromophore, respectively, whereas the I band belongs to a non-equilibrium structure with a ionized chromophore inserted in an unrelaxed chromophore (i.e. yet belonging to a neutral chromophore). It has been known for some time that absorption maxima of A and B bands are related in a 6:1 proportion, which (provided that molar absorptivities are the same) can be expressed as a free energy difference of ∆ion G0 ∼ 1 kcal mol−1 . In other words, there are almost degenerate situations connecting the migration of a proton accross a long distance. Another piece of experimental evidence points in the same direction: mutation S65T causes the most stable species in S0 to contain the ionized chromophore. 15 The small steric hindrance added by the extra methyl group in S65T is responsible for this difference, so it is natural to conclude that the two forms of the chromophore in wtGFP must be very similar in terms of energy. In a recent computational work by Grigorenko et al., the relative energies of single structures belonging to the A and B species in the ground state were determined within a QM/MM frame-


25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

a QM(CASPT2//CASSCF)/AMBER study (that is, a QM/MM-level simulation) of the topography of the excited state energetics using energy optimizations. Using this procedure the authors determine that proton transfer can occur in the excited state via a barrierless mechanism and that downhill relaxation channels are open to the system. The authors also identify a possible deactivation channel involving a concerted asynchronous hula-twist photoisomerization that explains the lower quantum yield of this mutant with respect to the wtGFP. Even though the theoretical methods used to derive the conclusions in their work and ours are different, the conclusions are not contradictory, but rather complementary. In fact, our study does not explore the excited-state energetics beyond the event of photoexcitation and first stages of proton transfer. Our analysis of the ground state is quite thorough and serves the purpose of establishing the energetics of the system prior to photoexcitation. We find that the protein environment and dynamical behavior of the protein are determinant to explain that most photoexcitation events will leave the system in the excited state in what can be called effective degeneracy, or in other words, a potential energy surface that, in what concerns the sole motion of the proton, is flat or favorable. This can explain why time-resolved measurements detect the presence of deprotonated chromophore already at the experimental limit (50 fs). Of course, once in the excited state the system must evolve somehow and find ways to relax. In effect, experimental data shows time evolution of the excited state populations when monitored further to the red at longer times. Our study does not occupy itself with this. The results obtained by Zhang et al. using S1 -based optimization suggest that the relaxation will involve a substantial increase of the O–O distance between chromophore and Asp148. 67

Page 26 of 32

and 42. This material is available free of charge via the Internet at http://pubs.acs.org/.

Acknowledgments This work was supported by the “Ministerio de Econom´ıa y Competitividad” through project CTQ2011-24292 and by the “Generalitat de Catalunya” through project 2009SGR409. Use of computational facilities at the “Centre de Serveis Cient´ıfics i Acadèmics de Catalunya (CESCA)” is gratefully acknowledged. Protein graphics on left panels of Figures 1 and 2 were prepared with the UCSF Chimera package: 68,69 Chimera is developed by the Resource for Biocomputing, Visualization, and Informatics at the University of California, San Francisco (supported by NIGMS P41-GM103311). Molecular graphics on right panels of Figures 1 and 2 were prepared with the MOLEKEL package. 70

References (1) Zimmer, M. Green fluorescent protein (GFP): Applications, Structure, and Related Photophysical Behavior. Chem. Rev. 2002, 102, 759 – 781. (2) Meech, S. R. Excited State Reactions in Fluorescent Proteins. Chem. Soc. Rev. 2009, 38, 2922–2934. (3) Lin, M. Z.; McKeown, M. R.; Ng, H.-L.; Aguilera, T. A.; Shaner, N. C.; Campbell, R. E.; Adams, S. R.; Gross, L. A.; Ma, W.; Alber, T. et al. Autofluorescent Proteins with Excitation in the Optical Window for Intravital Imaging in Mammals. Chem. Biol. 2009, 16, 1169–1179. (4) Bravaya, K. B.; Grigorenko, B. L.; Nemukhin, A. V.; Krylov, A. I. Quantum Chemistry Behind Bioimaging: Insights from Ab Initio Studies of Fluorescent Proteins and Their Chromophores. Acc. Chem. Res. 2012, 45, 265–275.

Supporting Information Available: A scheme detailing the method used to compute lengthy QM/MM dynamics simulations. Description of a MM simulation on S65T/H148D that illustrates the effect of water molecules. Complete list of authors of References 3, 27, 41

(5) Subach, F. V.; Verkhusha, V. V. Chromophore Transformations in Red Fluo-


26

Page 27 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(14) Yang, F.; Moss, L. G.; Phillips, G. N. The Molecular Structure of Green Fluorescent Protein. Nat. Biotechnol. 1996, 14, 1246 – 1251.

rescent Proteins. Chem. Rev. 2012, 112, 4308–4327. (6) Shcherbakova, D. M.; Subach, O. M.; Verkhusha, V. V. Red Fluorescent Proteins: Advanced Imaging Applications and Future Design. Angew. Chem. Int. Ed. 2012, 51, 10724–10738.

(15) Brejc, K.; Sixma, T. K.; Kitts, P. A.; Kain, S. R.; Tsien, R. Y.; Ormö, M.; Remington, S. J. Structural Basis for Dual Excitation and Photoisomerization of the Aequorea Victoria Green Fluorescent Protein. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 2306 – 2311.

(7) Voss Andreae, J. Protein Sculptures: Life’s Building Blocks Inspire Art. Leonardo 2005, 38, 41–45. (8) Chattoraj, M.; King, B. A.; Bublitz, G. U.; Boxer, S. G. UltraFast Excited State Dynamics in Green Fluorescent Protein: Multiple States and Proton Transfer. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 8362 – 8367.

(16) Wachter, R. M.; King, B. A.; Heim, R.; Kallio, K.; Tsien, R. Y.; Boxer, S. G.; Remington, S. J. Crystal Structure and Photodynamic Behavior of the Blue Emission Variant Y66H/Y145F of Green Fluorescent Protein. Biochemistry 1997, 36, 9759–9765.

(9) Stoner Ma, D.; Jaye, A. A.; Matousek, P.; Towrie, M.; Meech, S. R.; Tonge, P. J. Observation of Excited-State Proton Transfer in Green Fluorescent Protein Using Ultrafast Vibrational Spectroscopy. J. Am. Chem. Soc. 2005, 127, 2864 – 2865.

(17) Elsliger, M. A.; Wachter, R. M.; Hanson, G. T.; Kallio, K.; Remington, S. J. Structural and Spectral Response of Green Fluorescent Protein Variants to Dhanges in pH. Biochemistry 1999, 38, 5296–5301.

(10) Stoner Ma, D.; Melief, E. H.; Nappa, J.; Ronayne, K. L.; Tonge, P. J.; Meech, S. R. Proton Relay Reaction in Green Fluorescent Protein (GFP): PolarizationResolved Ultrafast Vibrational Spectroscopy of Isotopically Edited GFP. J. Phys. Chem. B 2006, 110, 22009 – 22018.

(18) Kneen, M.; Farinas, J.; Li, Y. X.; Verkman, A. S. Green Fluorescent Protein as a Noninvasive Intracellular pH Indicator. Biophys. J. 1998, 74, 1591–1599. (19) Llopis, J.; McCaffery, J. M.; Miyawaki, A.; Farquhar, M. G.; Tsien, R. Y. Measurement of Cytosolic, Mitochondrial, and Golgi pH in Single Living Cells with Green Fluorescent Proteins. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 6803–6808.

(11) van Thor, J. J.; Zanetti, G.; Ronayne, K. L.; Towrie, M. Structural Events in the Photocycle of Green Fluorescent Protein. J. Phys. Chem. B 2005, 109, 16099 – 16108.

(20) Shu, X.; Kallio, K.; Shi, X.; Abbyad, P.; Kanchanawong, P.; Childs, W.; Boxer, S. G.; Remington, S. J. Ultrafast Excited-state Dynamics in the Green Fluorescent Protein Variant S65T/H148D. 1. Mutagenesis and Structural Studies. Biochemistry 2007, 46, 12005–12013.

(12) Heim, R.; Cubitt, A. B.; Tsien, R. Y. Improved Green Fluorescence. Nature 1995, 373, 663 – 664. (13) Ormö, M.; Cubitt, A. B.; Kallio, K.; Gross, L. A.; Tsien, R. Y.; Remington, S. J. Crystal Structure of the Aequorea Victoria Green Fluorescent Protein. Science 1996, 273, 1392 – 1395.

(21) Shi, X.; Abbyad, P.; Shu, X.; Kallio, K.; Kanchanawong, P.; Childs, W.; Remington, S. J.; Boxer, S. G. Ultrafast Excited-


27


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

State Dynamics in the Green Fluorescent Protein Variant S65T/H148D. 2. Unusual Photophysical Properties. Biochemistry 2007, 46, 12014–12025.

Page 28 of 32

(28) Parr, R. G.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford University Press: Oxford, U. K., 1989. (29) Bauernschmitt, R.; Ahlrichs, R. Treatment of Electronic Excitations within the Adiabatic Approximation of Time Dependent Density Functional Theory. Chem. Phys. Lett. 1996, 256, 454 – 464.

(22) Leiderman, P.; Genosar, L.; Huppert, D.; Shu, X.; Remington, S. J.; Solntsev, K. M.; Tolbert, L. M. Ultrafast Excited-State Dynamics in the Green Fluorescent Protein Variant S65T/H148D. 3. Short- and Long-Time Dynamics of the Excited-State Proton Transfer. Biochemistry 2007, 46, 12026–12036.

(30) Dreuw, A.; Head Gordon, M. SingleReference ab Initio Methods for the Calculation of Excited States of Large Molecules. Chem. Rev. 2005, 105, 4009– 4037.

(23) Stoner Ma, D.; Jaye, A. A.; Ronayne, K. L.; Nappa, J.; Meech, S. R.; Tonge, P. J. An Alternate Proton Acceptor for Excited-State Proton Transfer in Green Fluorescent Protein: Rewiring GFP. J. Am. Chem. Soc. 2008, 130, 1227–1235.

(31) Yanai, T.; Dew, D. P.; Handy, N. C. A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP). Chem. Phys. Lett. 2004, 393, 51–57.

(24) Kondo, M.; Heisler, I. A.; Stoner Ma, D.; Tonge, P. J.; Meech, S. R. Ultrafast Dynamics of Protein Proton Transfer on Short Hydrogen Bond Potential Energy Surfaces: S65T/H148D GFP. J. Am. Chem. Soc. 2010, 132, 1452–1453.

(32) List, N. H.; Olsen, J. M.; Rocha-Rinza, T.; Christiansen, O.; Kongsted, J. Performance of Popular xc-Functionals for the Description of Excitation Energies in GFP-like Chromophore Models. Int. J. Quantum Chem. 2012, 112, 789–800.

(25) Scharnagl, C.; Raupp Kossmann, R. A. Solution pKa Values of the Green Fluorescent Protein Chromophore from Hybrid Quantum-Classical Calculations. J. Phys. Chem. B 2004, 108, 477–489.

(33) Leang, S. S.; Zahariev, F.; Gordon, M. S. Benchmarking the Performance of TimeDependent Density Functional Methods. J. Chem. Phys. 2012, 136, 104101. (34) Sobolewski, A. L.; Domcke, W. Ab Initio Potential Energy Functions for Excited State Intramolecular Proton Transfer: A Comparative Study of oHydroxybenzaldehyde, Salicylic Acid and 7-Hydroxy-1-indanone. Phys. Chem. Chem. Phys. 1999, 1, 3065–3072.

(26) Bell, A. F.; He, X.; Wachter, R. M.; Tonge, P. J. Probing the Ground State Structure of the Green Fluorescent Protein Chromophore Using Raman Spectroscopy. Biochemistry 2000, 39, 4423– 4431.

(35) Miertus, S.; Scrocco, E.; Tomasi, J. Electrostatic Interaction of a Solute with a Continuum. A Direct Utilization of Ab Initio Molecular Potentials for the Prevision of Solvent Effects. Chem. Phys. 1981, 55, 117 – 129.

(27) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian 09 Revision C.01. Gaussian Inc. Wallingford CT, 2009.

(36) Li, H.; Robertson, A. D.; Jensen, J. H. Very Fast Empirical Prediction and Inter-


28

Page 29 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


pretation of Protein pKa Values. Proteins 2008, 61, 704–721.

lar Dynamics Simulation Package. J. Mol. Graphics 1996, 14, 136–141.

(37) Bas, D. C.; Rogers, D. M.; Jensen, J. H. Very Fast Prediction and Rationalization of pKa Values for Protein-Ligand Complexes. Proteins 2008, 73, 765–783.

(45) Billeter, S. R.; Turner, A. J.; Thiel, W. Linear Scaling Geometry Optimisation and Transition State Search in Hybrid Delocalised Internal Coordinates. Phys. Chem. Chem. Phys. 2000, 2, 2177–2186.

(38) Olsson, M. H. M.; Sondergard, C. R.; Rostkowski, M.; Jensen, J. H. PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa predictions. J. Chem. Theory Comput. 2011, 7, 525–537.

(46) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (47) McQuarrie, D. A. Statistical Mechanics; University Science Books: Sausalito (California), 2000.

(39) Sondergard, C. R.; Olsson, M. H. M.; Rostkowski, M.; Jensen, J. H. Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of pKa Values. J. Chem. Theory Comput. 2011, 7, 2284–2295.

(48) Tuttle, T.; Thiel, W. OMx-d: Semiempirical Methods with Orthogonalization and Dispersion Corrections. Implementation and Biochemical Application. Phys. Chem. Chem. Phys. 2008, 10, 2159–2166.

(40) ChemShell, a Computational Chemistry Shell. See www.chemshell.org, Last accessed on: August 17, 2014.

(49) Korth, M.; Thiel, W. Benchmarking Semiempirical Methods for Thermochemistry, Kinetics, and Noncovalent Interactions: OMx Methods are Almost as Accurate and Robust as DFT-GGA Methods for Organic Molecules. J. Chem. Theory Comput. 2011, 7, 2929–2936.

(41) Sherwood, P.; de Vries, A. H.; Guest, M. F.; Schreckenbach, G.; Catlow, C. R. A.; French, S. A.; Sokol, A. A.; Bromley, S. T.; Thiel, W.; Turner, A. J. et al. QUASI: A General Purpose Implementation of the QM/MM Approach and its Application to Problems in Catalysis. J. Mol. Struct. (Theochem) 2003, 632, 1 – 28.

(50) Scholten, M. Ph.D. Thesis, Universität D¨ usseldorf, 2003. (51) Thiel, W. MNDO99 Program, Version 7.0. Max-Planck-Institut f¨ ur Kohlenforschung, M¨ ulheim, Germany, 2005.

(42) MacKerell, A. D.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S. et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586–3616.

(52) Colbert, D. T.; Miller, W. H. A Novel Discrete Variable Representation for Quantum Mechanical Reactive Scattering via the S-Matrix Kohn Method. J. Chem. Phys. 1992, 96, 1982 – 1991. (53) Sakurai, S. S.; Napolitano, J. J. Modern Quantum Mechanics, 2nd International Edition; Pearson, 2013.

(43) MacKerell, A. D.; Feig, M.; Brooks, C. L. Improved Treatment of the Protein Backbone in Empirical Force Fields. J. Am. Chem. Soc. 2004, 126, 698–699.

(54) Nadal-Ferret, M.; Gelabert, R.; Moreno, M.; Lluch, J. M. How Does the Environment Affect the Absorption Spectrum of the Fluorescent Protein

(44) Smith, W.; Forester, T. R. DL-POLY 2.0: A General-Purpose Parallel Molecu-


29


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 32

Andersen, L. H. Absorption Spectrum of the Green Fluorescent Protein Chromophore Anion In Vacuo. Phys. Rev. Lett. 2001, 87, 228102.

mKeima? J. Chem. Theory Comput. 2013, 9, 1731–1742. (55) Randino, C.; Nadal-Ferret, M.; Gelabert, R.; Moreno, M.; Lluch, J. M. A Time-Dependent DFT/Molecular Dynamics Study of the Proton-Wire Responsible for the Red Fluorescence in the LSSmKate2 Protein. Theor. Chem. Acc. 2013, 132, 1327.

(62) Kohlhoff, K. J.; Shukla, D.; Lawrenz, M.; Bowman, G. R.; Konerding, D. E.; Belov, D.; Altman, R. B.; Pande, V. S. Cloud-Based Simulations on Google Exacycle Reveal Ligand Modulation of GPCR Activation Pathways. Nature Chem. 2014, 6, 15–21.

(56) Nadal-Ferret, M.; Gelabert, R.; Moreno, M.; Lluch, J. M. Are There Really Low-Barrier Hydrogen Bonds in Proteins? The Case of Photoactive Yellow Protein. J. Am. Chem. Soc. 2014, 136, 3542–3552.

(63) Fayngold, M.; Fayngold, V. Quantum Mechanics and Quantum Information; Wiley-WCH, 2013. (64) Klinman, J. P. Importance of Protein Dynamics During Enzymatic C–H Bond Cleavage Catalysis. Biochemistry 2013, 52, 2068 – 2077.

(57) Ortiz-Sánchez, J. M.; Gelabert, R.; Moreno, M.; Lluch, J. M. Theoretical Study on the Excited-State Intramolecular Proton Transfer in the Aromatic Schiff Base Salicylidene Methylamine: An Electronic Structure and Quantum Dynamical Approach. J. Phys. Chem. A 2006, 110, 4649–4656.

(65) Grigorenko, B. L.; Nemukhin, A. V.; Polyakov, I. V.; Morozov, D. I.; Krylov, A. I. First-Principles Characterization of the Energy Landscape and Optical Spectra of Green Fluorescent Protein Along the A → I → B Proton Transfer Route. J. Am. Chem. Soc. 2013, 135, 11541 – 11549.

(58) Ortiz-Sánchez, J. M.; Gelabert, R.; Moreno, M.; Lluch, J. M. ElectronicStructure and Quantum Dynamical Study of the Photochromism of the Aromatic Schiff Base Salicylideneaniline. J. Chem. Phys. 2008, 129, 214308.

(66) Di Donato, M.; van Wilderen, L. J.; Van Stokkum, I. H. M.; Stuart, T. C.; Kennis, J. T. M.; Hellingwerf, K. J.; van Grondelle, R.; Groot, M. L. Proton Transfer Events in GFP. Phys. Chem. Chem. Phys. 2011, 13, 16295 – 16305.

(59) Voityuk, A. A.; Michel Beyerle, M. E.; Rosch, N. Structure and Rotation Barriers for Ground and Excited States of the Isolated Chromophore of the Green Fluorescent Protein. Chem. Phys. Lett. 1998, 296, 269 – 276.

(67) Zhang, Q.; Chen, X.; Cui, G.; Fang, W.; Thiel., W. Concerted Asynchronous Hula-Twist Photoisomerization in the S65T/H148D Mutant of Green Fluorescent Protein. Angew. Chem. Int. Ed. 2014, 53, 8649–8653.

(60) Niwa, H.; Inouye, S.; Hirano, T.; Matsuno, T.; Kojima, S.; Kubota, M.; Ohashi, M.; Tsuji, F. I. Chemical Nature of the Light Emitter of the Aequorea Green Fluorescent Protein. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 13617 – 13622.

(68) UCSF Chimera. See www.cgl.ucsf.edu/chimera, Last accessed on: August 17, 2014. (69) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E.

(61) Nielsen, S. B.; Lapierre, A.; Andersen, J. U.; Pedersen, U. V.; Tomita, S.;


30

Page 31 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


UCSF Chimera–A Visualization System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25, 1605 – 1612. (70) Varetto, U. MOLEKEL, Version 5.4.0.8. Swiss National Supercomputing Centre, Lugano, Switzerland, 2009.


31


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Graphical TOC Entry


32

Page 32 of 32

Unveiling How an Archetypal Fluorescent Protein Operates

Recommend Documents