Photoinduced Isomerization of the Photoactive Yellow Protein (PYP

ARTICLE pubs.acs.org/JPCA

Photoinduced Isomerization of the Photoactive Yellow Protein (PYP) Chromophore: Interplay of Two Torsions, a HOOP Mode and Hydrogen Bonding Evgeniy V. Gromov,*,†,‡ Irene Burghardt,§ Horst K€oppel,† and Lorenz S. Cederbaum† †

Theoretische Chemie, Physikalisch-Chemisches Institut, Universit€at Heidelberg, Im Neuenheimer Feld 229, D-69120 Heidelberg, Germany ‡ Laboratory of Quantum Chemistry, Computer Center, Irkutsk State University, K. Marks 1, 664003 Irkutsk, Russian Federation § Institute of Physical and Theoretical Chemistry, Goethe University Frankfurt, Max-von-Laue-Str. 7, D60438 Frankfurt/Main, Germany

bS Supporting Information ABSTRACT: We report on a detailed theoretical analysis, based on extensive ab initio calculations at the CC2 level, of the S1 potential energy surface (PES) of the photoactive yellow protein (PYP) chromophore. The chromophore’s photoisomerization pathway is shown to be fairly complex, involving an intimate coupling between single-bond and double-bond torsions. Furthermore, these torsional modes are shown to couple to a third coordinate of hydrogen out-of-plane (HOOP) type whose role in the isomerization is here identified for the first time. In addition, it is demonstrated that hydrogen bonding at the phenolate moiety of the chromophore can hinder the single-bond torsion and thus facilitates double-bond isomerization. These results suggest that the interplay between intramolecular factors and H-bonding determines the isomerization in native PYP.

I. INTRODUCTION Chromophores are ubiquitous in nature. Some of them are elegantly embedded in proteins, resulting in species able to interact with light. Such species, called photoactive proteins, or photoreceptors, are widely used in nature to make living organisms perceptive to light, a property that is largely utilized for the biological function. Examples range from the striking phenomena of vision and phototropism to less known but important functionalities such as photokinesis and photomotility1 (see ref 2 for an overview). Understanding how the “chromophore protein” pair functions and what the key elements are that make its functioning so efficient has been one of the main goals of contemporary photochemistry and photobiology.3 In the present study we address the primary steps in the photoresponse of Halorhodospira halophila, a bacterium that exhibits a negative phototaxis (photomotility) when exposed to blue light.4 The chomophoreprotein pair mediating this response is the photoactive yellow protein (PYP),5 a small protein comprising 125 amino acids and a unique cofactor, the deprotonated p-hydroxycinnamoyl chromophore (Figure 1). In vivo, PYP functions via a photocycle that is initiated by the isomerization of the chromophore’s double bond,69 very similar to retinal photoisomerization in the visual process.10 The characteristic time scale of the photoisomerization in PYP is in the subpicosecond r 2011 American Chemical Society

time range. The overall photocycle extends over a few seconds and involves the protonation of the chromophore (on a microsecond time scale) as well as the partial unfolding of the protein (on a millisecond time scale).7 PYP has been studied extensively over the past few years resulting in a detailed picture of the various steps of the photocycle (for several recent reviews see, e.g., refs 1116). The recent experimental1724 and theoretical2532 studies that specifically focus on the very initial steps of the cycle have addressed different aspects in the mechanism of the chromophore’s isomerization. Substantial efforts were made, in particular, to unravel the influence of the environment (protein or solution)1719,25,27,29,32,33 and the effects of the substituent in the carbonyl group of the chromophore.31,3436 Despite these efforts, there is no coherent picture yet of the chromophore’s isomerization. This concerns, in particular, the detailed mechanism of the process. Most of the earlier studies (see, e.g., refs 7 and 32) conceived the isomerization as being mainly the torsion around the C7dC8 double bond (Figure 1). In our previous work31 (in the following referred to as paper I) we Received: February 4, 2011 Revised: July 6, 2011 Published: July 11, 2011 9237

dx.doi.org/10.1021/jp2011843 | J. Phys. Chem. A 2011, 115, 9237–9248

The Journal of Physical Chemistry A

Figure 1. Chemical structure of the PYP chromophore: a schematic representation of the chromophore in the protein (left), and the pCTMe 3 2H2O (middle) and pCTMe 3 4H2O (right) species considered in the present study. For the chromophore-in-protein, the actual H-bonds occurring between the chromophore and the protein environment (neighboring amino acids) are schematically shown. The two torsions, around the single bond (R-torsion) and around the double bond (β-torsion) are shown for the pCTMe 3 2H2O species. For both torsional angles, the direction of the arrow to the left defines the counterclockwise (CCW) direction of rotation.

established, however, that the chromophore can also undergo torsion around the adjacent C4C7 single bond which was found to occur barrierless. At the same time isomerization around the C7dC8 bond, also considered in I, involves some barrier. In subsequent work, Groenhof and co-workers proposed that both the single-bond and double-bond isomerization occur in the native PYP chromophore, though the isomerization yield for the single bond is smaller than for the double bond.27,29 This singlebond mechanism has not been yet confirmed experimentally and has been the subject of a lively discussion in the literature (see, e.g., refs 17, 23, 27, and 31). Note that the double-bond mechanism was thought to be promoted by the positive charge of the Arg52 amino acid which is in immediate proximity of the chromophore in the native protein.32 For both mechanisms, the H-bonds between the chromophore and environment (protein or solution) were proposed to facilitate the ultrafast S1S0 decay from the twisted single-bond and double-bond minima.27 It should be however noted that the predicted crucial role of Arg52 has not been corroborated experimentally.19 Moreover, a recent neutron crystallographic study revealed that Arg52 seems to be neutral in native PYP.37 It is also noteworthy that the single-bond torsion was completely neglected in the recent theoretical work of ref 26. In our view, such an approximation is questionable when treating the PYP chromophore. Obviously, there are many open questions in the photochemistry of the PYP chromophore in general and in its photoinduced torsion in particular. In the work reported here we therefore aim at scrutinizing the single vs double bond torsion in order to shed some more light on that issue. To this end we carefully analyze the potential energy surface (PES) of the lowest excited state (S1) of the p-coumaric thiomethyl chromophore, hydrogen bonded with two water molecules (pCTMe 3 2H2O, Figure 1). Special attention will be paid to the topology of S1 along the single-bond and double-bond torsions. In particular, the first half of the transcis isomerization pathway (leading to the 90 double-bond twisted structure) will be analyzed in some detail. The part of the S1 PES related to the cis configuration of the

ARTICLE

chromophore is not considered in the present study. In addition to the pCTMe 3 2H2O species, we will also address the chromophore hydrogen bonded with four water molecules (pCTMe 3 4H2O, Figure 1) aiming at elucidating the effect due to the H-bonds. Note that microsolvated complexes of the PYP chromophore are of interest by themselves in view of several recent experimental studies.38,39 As will be apparent below, the S1 PES of pCTMe 3 2H2O has a very complex topology, exhibiting many minima and saddle points, in fact much more complicated than predicted in I. In this context, the present study is a substantial extension of I. It illuminates many previously obtained findings and brings new important results. These include (i) revealing a new coordinate, hydrogen out-of-plane (HOOP) mode, that is intimately involved in the transcis isomerization of pCTMe 3 2H2O and (ii) identifying the stabilization effect of the H-bonds that is anticipated to hinder (or even entirely suppress) the single-bond mechanism in the native protein. The remainder of the paper is organized as follows. In section II we give details of the calculational procedures. In section III, we present and analyze our results. Section IV summarizes and concludes.

II. COMPUTATIONAL PROCEDURES In the present study we make use throughout of the secondorder approximate coupled cluster (CC2) method,40 for calculating the electronic structure and properties of the ground (S0) and first excited (S1) states of pCTMe 3 2H2O and related species. The method proved to properly describe the S1 excited state of the chromophore which is of ππ* character.31 The deviations of the CC2 excitation energy for this state from the results of the more accurate EOM-CCSD and CR-EOM-CCSD(T) methods were found to be fairly small, ∼0.17 and ∼0.09 eV, respectively.41 A full and partial geometry optimization for the S1 state was carried out using the analytical gradient technique, with a threshold of convergence in the geometry optimization procedure (maximum norm of Cartesian gradient) of 104. No geometrical constraints were imposed in the search (optimization) of the transition states (saddle points). At each stationary point a vibrational analysis (at the same CC2 level) was carried out to establish/confirm the type of the point. The required Hessians were computed by numerical differentiation of the gradients. For each stationary point the atomic charges were computed using the natural population analysis (NBA) approach.42 In addition to various partial geometry optimizations, used to sample the S1 PES along certain coordinates, we also calculated an intrinsic reaction coordinate (IRC) path. A special geometry optimization procedure (with a hypersphere constraint), implemented in the program system Molcas43 (module Slapaf), was used to follow the IRC path (see details of the procedure in ref 43). The IRC path was calculated in mass-weighted coordinates. The value of the hypersphere radius was 0.1 a.u. The CC2 calculations were performed with the use of the resolution-of-the-identity (RI) approximation (RI-CC2),44,45 employing the Turbomole program package (V6.0).46 The available version of the package lacks the procedure to perform IRC calculations; therefore, an external program (Molcas) was used for that. An auxiliary I/O interface was employed in the IRC calculations to link the Turbomole and Molcas modules. Unless stated otherwise, the correlation consistent valence polarized double-ζ basis set of Dunning, augmented by a set of 9238

dx.doi.org/10.1021/jp2011843 |J. Phys. Chem. A 2011, 115, 9237–9248


ARTICLE

Figure 2. Diagram of the stationary points on the S1 PES of pCTMe 3 2H2O chromophore. The S1 and S0 energies and the S1S0 energy gap at each stationary point are shown. The CC2 S0 energy at the ground-state equilibrium geometry S0,eq was taken as origin. The oscillator strength (a.u.) for the vertical S0S1 transition (at the S0,eq point) is given in parentheses. For each stationary point a perspective (3D) drawing of the corresponding molecular structure of pCTMe 3 2H2O is given (the water molecules are omitted to avoid congesting the figure); one generalized structure is given for the S0,eq and S1,min points. The dashed line indicates a spurious link between the S1,min minimum and the Sβ1,TS transition state (see text for details).

diffuse functions (aug-cc-pVDZ)47 was used in the CC2 calculations. The 19 core molecular orbitals were kept frozen throughout in the coupled cluster calculations. All computations were performed in parallel, with 64 CPU-cores being utilized on average. It should be emphasized that the use of diffuse functions is important to describe the S1 PES correctly. Although the lowering of the S1 excitation energy due to inclusion of diffuse functions is relatively moderate (∼0.2 eV), there is a pronounced effect on the topology of the S1 PES. This concerns in particular the species H-bonded via the phenolic oxygen of the chromophore, that is, the pCTMe 3 2H2O and pCTMe 3 4H2O complexes considered here as well as the chromophore in the native PYP environment and in aqueous solution. It is worthwhile to note the suitability of the CC2 method for the present study. As amply demonstrated in I, the method has proven successful in describing the S1 PES of pCTMe 3 2H2O, in the FranckCondon region and at the single- and doublebond twisted configurations. Note that at the twisted structures the S0 and S1 states become close in energy, with the S1S0 energy gap being very small (∼0.4 eV) at the double-bond twisted configuration (see section III for details). In such cases there is usually a strong coupling of the S0 and S1 states, which entails a breakdown of single-reference based coupled-cluster methods. (The problem essentially relates to strong multiconfigurational character of the S0 (reference) state due to its coupling with the S1 state).48,49 There are, however, exceptions that occur when S0 and S1 do not interact (for instance, due to different

spatial symmetry, see, e.g., ref 50) or their interaction (coupling) is small. The latter seems to be the case for the pCTMe 3 2H2O chromophore and, expectedly, related species. The reason for the negligible S0S1 coupling is the different charge distribution in the S0 and S1 states, especially at the twisted configurations (see section III for details). Note, however, that the coupling depends on the chromophore’s geometry and can become larger at certain configurations. This concerns in particular the structures exhibiting both the R- and β-twists to a substantial extent. In the present study we avoided calculations of the S1 state at the chromophore’s configurations where CC2 was known to give unreliable results. In those cases, to obtain the data (S1 energy) interpolation or extrapolation was carried out employing the reliable data obtained at other configurations. In this way the 2D cuts of the S1 PES discussed in section III.E were constructed, with the MATHEMATICA software51 used for the interpolation/extrapolation procedure. We emphasize that this work does not address conical intersections between the S0 and S1 PES, but attention is focused on the topology of the S1 PES. Finally, in most our CC2 calculations, the D1 diagnostic52 was 85%, which should ensure robustness of our results. In a few cases along the IRC path a somewhat elevated value of D1 (∼0.2) was observed. We conclude this section with a note on the CASSCF/ CASPT2 approach which has been widely used to study the PYP chromophore (see, e.g., refs 26, 27, 29, 32, and 53). The 9239



ARTICLE

Table 1. Relative Changes of the Ground State (S0) and Excited State (S1) Energies with Account for the Zero-Point Vibrational Energy (ZPVE) at the Different Stationary Points of pCTMe 3 2H2O (All Values in eV) point

type

S0

S1

S1 (ZPVE)

Table 3. Ground and Excited State Integrated Natural Charge Distributions (Calculated from the Natural Population Analysis) Partitioned between the Chromophore Ring (Including the Two Water Molecules) and the Alkene Part of the pCTMe 3 2H2O Chromophore at the Different Stationary Points

S0,eq

min

0.00

2.92

S1,min

min

0.12

2.82

2.70

SR1,TS

saddle

0.14

2.82

2.69

SR1 1,min

min

1.04

2.66

2.56

SR12 1,TS

S0,eq S1,min SR1,TS SR1 1,min chrom. ring 0.56 0.49 0.58 0.48 0.58 0.47 0.73 0.26

saddle

0.95

2.67

2.55

SR2 1,min

alkene part 0.44 0.51 0.42 0.52 0.42 0.53 0.27 0.74

min

1.04

2.66

2.56

Sβ1,TS Sβ1,min

saddle min

0.82 1.98

3.08 2.50

2.91 2.44

Table 2. Selected Bond Lengths (Å), r, β, and γ Torsional Angles (), and the Out-of-Plane Deviations () of H7 and H8 (outp7 and outp8)a for pCTMe 3 2H2O at the Different Stationary Points bond/angle S0,eq S1,min SR1,TS

SR1 1,min

SR12 1,TS

SR2 1,min

Sβ1,TS Sβ1,min

O1C1

1.306 1.326 1.326

1.304

1.302

1.304

1.307

1.320

C1C2

1.448 1.431 1.432

1.446

1.446

1.445

1.443

1.441

C2C3 C3C4

1.393 1.410 1.409 1.430 1.431 1.430

1.395 1.432

1.393 1.436

1.394 1.436

1.396 1.447

1.397 1.436

C4C5

1.430 1.423 1.423

1.436

1.436

1.432

1.442

1.436

C5C6

1.394 1.404 1.405

1.394

1.393

1.395

1.389

1.397

C6C1

1.443 1.441 1.441

1.445

1.446

1.446

1.455

1.439

C4C7

1.436 1.481 1.483

1.478

1.474

1.478

1.422

1.427

C7C8

1.386 1.405 1.404

1.424

1.421

1.424

1.475

1.465

C8C9

1.447 1.439 1.439

1.414

1.412

1.414

1.418

1.452

C9O2 C9S

1.247 1.265 1.265 1.857 1.870 1.871

1.271 1.905

1.273 1.909

1.271 1.904

1.267 1.877

1.242 1.852

SC10 H1 3 3 H2 3 3 —R

1.821 1.821 1.821

1.820

1.820

1.820

1.822

1.823

3 O1 1.729 1.782 1.783

1.839

1.844

1.840

1.763

1.702

3 O1 1.729 1.774 1.776

1.840

1.844

1.839

1.757

1.704

92.5

89.1

86.0

171.6 178.8

180.0 180.1 169.1

—β

180.0 180.0 179.1 165.7 179.8 165.8

133.1

87.3

—γ

180.0 180.0 178.4

172.3

179.9

172.3 140.6

90.3

outp7 outp8

0.0 0.0

20.6 1.5

0.4 0.1

20.6 1.5

1.2 1.5

0.0 0.0

1.8 0.4

11.4 3.1

a

outp7 and outp8 are used as measures of pyramidalization for the C4, C7, C8, H7 and C7, C8, C9, H8 atoms, respectively.

advantage of this method over the CC2 scheme is that it is suitable to treat degeneracies and quasidegeneracies of the S0 and S1 states thus making feasible to study S0S1 conical intersections and related phenomena. However, CASSCF/CASPT2 calculations become a formidable task for polyatomic molecules when a large active space is required. This is very likely the case for pCTMe 3 2H2O and related species when a diffuse basis set is employed. As we found, the use of diffuse functions results in the lowest virtual molecular orbitals (MO) that are lower in energy than the π* MO involved in the description of the S1 state. As a result the active space needed to describe this state in the case of a diffuse basis set will be much larger than employing a standard valence basis set due to the necessity to include the MOs lying below the π* MO.

S0

S1

S0

S1

S0

S1

S0

S1

SR2 Sβ1,TS Sβ1,min SR12 1,TS 1,min chrom. ring 0.74 0.23 0.73 0.25 0.47 0.48 0.25 0.69 alkene part 0.26 0.77 0.27 0.75 0.53 0.52 0.75 0.31

III. RESULTS AND DISCUSSION Figure 2 summarizes results of our calculations for the stationary points on the S1 PES of the pCTMe 3 2H2O chromophore. In the figure, the points are arranged along the lower and upper abscissa axes, with minima lying along the lower axis and saddle points indicated along the upper axis. For each point we give the relative energies of the S0 and S1 states, the S1S0 energy gap and a (3D-like) sketch of the molecular structure of the chromophore (the water molecules are omitted to avoid encumbering the figure). The S0S1 vertical transition at the equilibrium geometry of the ground state S0,eq [FranckCondon (FC) transition] is also shown; here, both the transition energy and oscillator strength are indicated. A detailed discussion of the results illustrated in Figure 2 will be given in the following sections. Complementary information is available in Tables 13, where we present quantitative data, relative values of the S0 and S1 energies (Table 1), selected geometrical parameters (Table 2) and the integrated natural charge distributions in the S0 and S1 states (Table 3) at the stationary points under consideration. Cartesian coordinates and the S0 and S1 total energies for the stationary points can be found in the Supporting Information. In this work we focus on the twists around the C4C7 and C7dC8 bonds as well as on the coupling between these and other coordinates. For the vertical electronic spectrum of the pCTMe 3 2H2O species as well as characteristics of the S1 excited state, in particular, its charge-transfer character, we refer to the discussion in paper I. We also note that we used different methods for calculating charge distributions in the present work (Table 3) and in paper I. In the latter work, the Mulliken scheme was employed, while here we use the more accurate NBA approach (see section II). Although there are some quantitative differences, both schemes are qualitatively consistent with regard to the charge separation at the stationary points under consideration. Following the analysis of paper I, we note that the stationary points in Figure 2 are associated with different charge distributions in the S1 state. At the planar minimum, S1,min, and especially at the single-bond twisted configurations of the chromophore, R12 R2 SR1 1,min, S1,TS, and S1,min, the negative charge is shifted away from the phenolate part toward the alkene moiety. In contrast, the double-bond twisted minimum, Sβ1,min, exhibits an inverted charge distribution, that is, the negative charge is localized on the phenolic ring (Table 3). This “atypical” charge distribution was found to relate to the different electronic structure at the double-bond twisted configuration (paper I). In particular, at 9240


The Journal of Physical Chemistry A Sβ1,min, the HOMO and LUMO, which mainly contribute to the S1 electronic wave function (with HOMO f LUMO being the principal configuration), are localized on the alkene part and the phenolate part, respectively. This is opposite to the pattern of these molecular orbitals at the single-bond twisted configurations. For further details on changes in the electronic structure of S1, see the discussion in paper I. In the discussion below we make use of the following nomenclatures. We denote the twist (torsion) around the C4C7 single bond as R-twist (R-torsion) and that around the C7dC8 double bond as β-twist (β-torsion). We will also use the notations R and β for the designation of the C5C4C7C8 and C4C7C8C9 dihedral angles, respectively (see Figure 1 for a graphical representation). When giving various values related to the R-/β-twists, we use the notations R-/β-twist (R-/ β-torsion) and just R/β. With the former we imply the “extent” of the twist, that is, the deviation of the corresponding angle (R or β) from 180, while the latter indicates the absolute value of the angles. We are also interested in the H7C7C8H8 dihedral angle that we denote γ. Where relevant, we distinguish the clockwise (CW) and counterclockwise (CCW) torsion of R and β, defined as shown in Figure 1. Numerically, the CW (CCW) torsion corresponds to a decrease (increase) of the corresponding angle. Both positive and corresponding negative values are used for R, β, and other dihedral angles larger than 180. The negative value of an angle implies the CCW torsion. A. In-Plane Relaxation of the Chromophore. According to both the results of the present calculations and our earlier work (paper I), the transition to the S1 state is followed by planar geometry relaxation ending at the S1,min point (Figure 2). The relaxation involves alterations of the bonds and valence angles, in particular an elongation of the C4—C7, C7dC8 and C9—S bonds (Table 2), and a notable decrease of the C4C7C8 angle, by ∼5 as compared with the ground state value. Note also a certain extension of the hydrogen bonds, ∼0.04 Å in comparison with the ground-state equilibrium geometry S0,eq. As discussed in paper I, this is indicative of a weakening of the H-bonding, which relates to a transfer of the electron density (negative charge) from the phenol ring toward the alkene part, occurring in the S1 excited state (Table 3). The adiabatic energy lowering for the S1,min point is ∼0.1 eV (Table 1), in agreement with the result obtained in paper I. The vibrational analysis at S1,min revealed no imaginary frequencies, which suggests that S1,min is a true minimum. Surprisingly, this finding does not agree with our prediction in paper I where this stationary point (denoted S1,TS) was assigned to the transition state with respect to the R-twist of the chromophore. Likewise, this result also disagrees with the earlier study of ref 54, where a similar species (i.e., the chromophore H-bonded via the phenolic oxygen) was considered. The different result obtained here was found to relate to the presence of diffuse functions in the basis set employed in the present study. Indeed, the vibrational analysis at the S1,TS point obtained in paper I, performed with the cc-pVDZ basis set (i.e., without diffuse functions), revealed one mode with imaginary frequency, corresponding in character to the R-twist of the chromophore. Since the present basis (aug-cc-pVDZ) is qualitatively improved (in particular for anionic species), we consider the present result correct. Furthermore, the diffuse functions were found to be crucial in describing the stabilization of the S1,min minimum against the R-torsion due to the H-bonding (see discussion in section III.F). For the pCTMe 3 2H2O

ARTICLE

chromophore, the stabilization turns out to be relatively weak. As a consequence, the S1,min well proved to be very shallow, such that it has no bound vibrational levels (see the discussion in the next section). This, however, is expected to change in the protein environment as discussed in section III.F. B. r-Twist of the Chromophore. According to our present results (Figure 2), there are four stationary points related to the R12 R-twist (i.e., torsion around the C4C7 bond), SR1,TS, SR1 1,min, S1,TS, R2 and S1,min, in contrast to only one point obtained in paper I. As seen from Figure 2, the first point along the R-torsional path is SR1,TS. It corresponds to the transition state that connects the planar minimum S1,min and the fully twisted SR1 1,min structure. The transition state character of SR1,TS was confirmed by vibrational analysis that revealed one mode with imaginary frequency of 33 cm1. As follows from Table 2, the structure of SR1,TS is similar to that of S1,min, with an exception for the R angle, which is about 10 smaller for SR1,TS. Energetically, the SR1,TS point is only ∼1 meV higher than the S1,min point (Table 1). This small barrier was found to “disappear” when the zero-point vibrational energy (ZPVE) is taken into account. In other words, the zeroth vibrational level of S1,min lies above that of SR1,TS, which suggests that the R-torsion of the pCTMe 3 2H2O chromophore should proceed in a “barrierless” fashion. Proceeding from the SR1,TS transition state toward larger torsional angles, one reaches three stationary points, SR1 1,min, R2 S R12 1,TS , and S 1,min (Figure 2). According to our calculations R2 R12 S R1 1,min and S 1,min represent minima, whereas S 1,TS is a transition state separating these minima. As can be seen, all three points are very close in energy, with almost identical energies R2 of the minima. (Note that S R1 1,min and S 1,min are not identical because of the presence of the two water molecules.) The barrier between the minima (S R12 1,TS ) is seen to be just 0.01 eV (Table 1). When the ZPVE is taken into account, the zero R2 vibrational levels of S R1 1,min and S 1,min turn out to be above the barrier (Table 1). In this case, the maximum of the groundstate vibrational function should occur at the transition state, that is, at the S R12 1,TS point. As follows from Table 2 (see also the images of the structures in Figure 2) the transition state SR12 1,TS is almost perfectly (90) twisted, whereas the minima deviate somewhat from 90 in different directions. Interestingly, they also exhibit a notable pyramidalization, that is, deviation of the H7 atom out of the C4C7C8 plane, that amounts to ∼20 in both cases (Table 2). (Note that the pyramidalization was found to be the coordinate R2 R12 that connects the SR1 1,min and S1,min minima via the S1,TS transition state.) These findings suggest a similarity of the R-torsion with a classical twist around the double bond in ethylene. As is wellknown, the pyramidalization in ethylene reflects a sp2 f sp3 rehybridization that in turn is a response to the double-bond breakage. The same effect seems to occur for the R-twisted chromophore (though in a less pronounced fashion than in ethylene), which points to some similarity between the R-twist and the twist around the double bond. On the other hand, the S1S0 energy gap at the R-twisted configurations is relatively large, ∼1.7 eV (Figure 2), which is not the case for the doublebond twist in ethylene. The four R-twisted structures occur along the counterclockwise (CCW) R-torsion, as indicated in Figure 2. The clockwise (CW) R-twist should also be characterized by four stationary R1 R2 points, where SR12 1,TS is the same, while S1,min and S1,min should be R exchanged. The transition state S1,TS of the CW-torsion will be of slightly different structure than along the CCW-twist. This 9241


The Journal of Physical Chemistry A relates to the presence of the two water molecules; in their absence, the CCW and CW twists would be equivalent. Anticipating the discussion in section III.D, the CCW R-twisted structures are coupled to the CW β-twist of the chromophore. With reference to paper I, we recall that the R-twist enhances the charge transfer from the phenol ring to the alkene part in the S1 state, with a direction that is opposite to the S0 ground state (Table 3, more details in paper I). The migration of the charge manifests itself in a further extension of the H-bonds as the R-twist increases (Table 2). These findings will be used later to elucidate the impact of the H-bonds on the R-twist (see section III.F). C. β-Twist of the Chromophore. Our calculations predict two stationary points on the S1 PES related to the β-twist of the chromophore, Sβ1,TS and Sβ1,min (Figure 2), where the former is a saddle point, whereas the latter corresponds to a minimum. As follows from Figure 2, the Sβ1,min point represents the fully β-twisted structure which was previously discussed in some detail in paper I. The present results fully agree with the findings of paper I regarding various properties of Sβ1,min, in particular, its geometry, energetics, and charge distribution. Several interesting points emerge from a comparison of the β- and R-twisted structures. In contrast to the R-twist, there is only one β-twisted configuration (both along the CW and CCW torsion), although the twist is not “perfect” (∼87). This difference from the R-twisted case has to do with the absence of pyramidalization at the Sβ1,min minimum, either for the C4, C7, C8, H7 or the C7,C8, C9, H8 atoms (Table 2). The latter in turn seems to relate to a persistence of the sp2 hybridization at the C7 and C8 atoms despite the breakage of the C7dC8 bond. This can be explained by conjugation of the p-orbitals of C7 and C8 with the corresponding p-orbitals of the adjacent, C4 and C9, atoms. This conjecture is confirmed by the lengths of the C4C7 and C8C9 bonds at Sβ1,min (1.427 and 1.452 Å, respectively), which are close to the lengths of these bonds in the electronic ground state S0,eq (Table 2), where such a conjugation is certainly present. In paper I we indirectly ascertained that there is a torsional barrier separating both the FC point and the R-twisted minimum from the Sβ1,min minimum; this barrier was estimated to occur at a β-twist of ∼55. In the present work we have explicitly identified the barrier which is associated with the Sβ1,TS saddle point. The transition-state character of the point was confirmed by a vibrational analysis that revealed one mode with imaginary frequency of 520.5 cm1. As follows from Table 1, the Sβ1,TS point lies about 0.16 eV higher than the vertical transition energy. The S0S1 energy difference at Sβ1,TS is 2.26 eV, relatively large in comparison to the R-twisted and β-twisted configurations (Figure 2). Interestingly, the charge distribution pattern at Sβ1,TS is practically the same for the S0 and S1 states, with the charge almost equally distributed between the phenolic ring and the alkene part (Table 3). Structurally (Table 2), the Sβ1,TS configuration is characterized by an appreciable twist along β, of about ∼47, very close to our prediction in paper I. There is also an R-twist (∼8) as well as some pyramidalization involving the C4,C7,C8,H7 atoms. Both the geometrical data and energetics for the Sβ1,TS structure agree fairly well with the results of ref 26 where this point was also obtained, although for a derivative of the chromophore considered here. However, in contrast to the predictions of that work we consider the Sβ1,TS transition state separating the Sβ1,min (β-twisted) and SR1 1,min (R-twisted) minima (the full line in Figure 2) rather than the β-twisted and the planar minima (the dashed line in Figure 2) as was proposed in ref 26. This issue

ARTICLE

Figure 3. IRC path on the S1 PES of the pCTMe 3 2H2O chromophore from the Sβ1,TS transition state. At the points marked by crosses, the chromophore was found to exhibit an instability with respect to the R-torsion as established from a vibrational analysis (see text for details).

is discussed in the next section. The disagreement with work of ref 26 might be due to differences in the chromophores and due to the fact that the R-twist was not taken into account in work of ref 26. D. Coupling between β- and r-Torsion. As has been mentioned above, the R and β torsions are not independent, but are coupled to each other. As a consequence, the pathway to the β-twisted structure turns out fairly complex. In this section, we analyze the correlations between the R and β twists from several viewpoints. First, the IRC path is constructed, from the Sβ1,TS transition state (section III.D.1). Second, the potential surface S1 as a function of R and β is analyzed in the vicinity of the FC point (section III.D.2). Third, a β torsional minimum energy path (MEP) in the R-twisted domain is analyzed, to shed light on the transition between the R and β twisted domains (section III.D.3). 1. IRC Analysis. In Figure 3, we present the IRC path from the Sβ1,TS transition state computed as described in section II. According to these calculations, the right (“plus”) branch of the IRC profile leads to the Sβ1,min minimum, while the left (“minus”) branch terminates at the S1,min planar minimum. This suggests that the Sβ1,TS transition state separates the planar minimum and the β-twisted structure. However, a full geometry optimization started from the “0.1” point of the IRC path surprisingly resulted in the R-twisted SR1 1,min minimum instead of the expected planar minimum. (Conversely, a full geometry relaxation from the “+0.1” structure led to the Sβ1,min minimum in agreement with the prediction of IRC.) The inconsistency between the minus branch of IRC and the result of the geometry optimization was found to relate to the instability of the chromophore with respect to the R-torsion that occurs along the minus branch. More precisely, the instability is observed for the first five IRC points, which were recorded from 0.1 to 0.5 a.u., indicated by crosses in Figure 3. This was established from a vibrational analysis performed at the IRC points that revealed one mode with imaginary frequency in the case of those points. The character of the mode was found to correspond to the R-torsion of the chromophore, which is different from the pattern of the mode with imaginary frequency at the transition state itself (see section III.E). These findings suggest that a part of the minus branch of the IRC proceeds along a ridge (or in proximity of a ridge) with respect to the R-twist. The IRC 9242



ARTICLE

Figure 5. Cross sections of the S1 and S0 PESs of the pCTMe 3 2H2O chromophore along the β torsional coordinate, for β changing from 180 to 90 and all the remaining coordinates relaxed (for the S1 state). The value of R (degree) is indicated for each point of the cross sections. The abrupt decrease (increase) of the S1 (S0) energy at β ≈ 173 relates to the R-twist of the chromophore (see text for details).

Figure 4. Upper panel: contour plot of the S1 PES with respect to the R and β torsional coordinates (all the remaining coordinates relaxed) in the vicinity of the S1,min planar minimum. The energy is indicated for certain contour lines; for each next line, the energy increment is 2 meV. Lower panel: several cuts through S1 along R at certain fixed values of β (all the remaining coordinates relaxed). Red dots (upper panel) indicate stationary points. Two crosses stand for valley-ridge inflection (VRI) points (see text for details).

procedure fails to predict the instability in such a case.55 Its failure relates to the fact that it follows a minimum gradient path, which in this case proceeds along the ridge rather than a valley and is not the minimum energy path. In the case of such a topology, one might expect branching (bifurcation) of the reaction path and related phenomena, for example, valley-ridge inflection (VRI) points. These are well-known in the literature, see, for example, ref 56 and references therein. In the present work, we refrain from a deeper analysis of the S1 topology along the IRC, which is beyond the scope of this study. In the next subsection we start exploring the S1 PES as a function of R and β in the vicinity of the FC geometry. The latter is important to shed some light on the initial pathway (toward the Sβ1,TS transition state) after excitation to S1. 2. R and β Torsion in the FC Region. To gain a better understanding of the R versus β torsions in the initial phase following photoexcitation, we focus here on the first portion of the path from the FC geometry toward the Sβ1,TS transition state. Our analysis is based on Figure 4 where we present a contour plot of the 2D cut through S1 with respect to the R and β coordinates (upper panel); several cuts along the R coordinate at fixed values of β are shown in the lower panel. In the figure, R and β varies from 160 to 200, which is about one-half of the pathway toward the Sβ1,TS transition state. At larger excursions of R and β, in particular, in the vicinity of the transition state, the pathway

becomes more complicated due to the involvement of another, third coordinate (see discussion in section III.E); in this case, the 2D picture is not meaningful. From Figure 4 (upper panel) one can clearly see the main topological features of S1 around the planar S1,min minimum, the minimum itself and two (CCW and CW) SR1,TS transition states, as discussed above. The topology of S1 is seen to change from a valley (for β = 180) to a slope (for β j 170) along the R coordinate (Figure 4, lower panel). Here we have identified two VRI points (denoted by crosses in Figure 4), for CW and CCW β twists, respectively. For this we used the condition ∂E/∂R = 0 and ∂2E/∂R2 = 0, which is satisfied at (R ∼ 185.7, β ∼ 173.6) and (R ∼ 174.3, β ∼ 186.7), VRI1 and VRI2, respectively, in Figure 4. Note that these points are not stationary (∂E/∂β 6¼ 0). Moreover, VRI1 and VRI2 belong to the domain of the S1 surface where all the coordinates except for R and β are optimized. One can however expect more VRI points at different finite values of the other coordinates. (This conjecture has been not checked due to the lack of the corresponding procedure). An interesting observation that emerges from Figure 4 (lower panel) is that there is an asymmetry for the cuts along R that gets more pronounced as β gets more twisted. As a result of this asymmetry, the direction of the R-twist depends on the direction of the β-twist. More specifically, the CW β twist started from β = 180 (full lines in lower panel of Figure 4), results in “decreasing” the barrier (that eventually changes to a slope) on the right side of the curves, giving rise to the CCW R twist. Vice versa, the CCW β-twist (dashed lines in Figure 4), in the same fashion, leads to the CW R-twist. In other words, the βR twisting motion is concerted. This is reflected in Figure 2 in terms of connections between the CW β-twisted transition state Sβ1,TS and CCW R-twisted minimum SR1 1,min. This correlation in the βR torsion suggests that there is no (or only little) bifurcation of the reaction path when starting from the planar minimum. A more precise and specific answer can be given from dynamical calculations, which we reserve for future work. 3. β-Torsion in the R-Twisted Domain. In the above discussion we have clearly demonstrated that, for the pCTMe 3 2H2O chromophore, a twist along β inevitably induces a concerted torsion along R, which should eventually lead to the R-twisted R12 R2 domain (the SR1 1,min, S1,TS, and S1,min points). Here we emphasize that once the R-twist has taken place the chromophore gets 9243



Figure 6. Left part: the pattern of the transition mode (the mode with imaginary frequency) obtained from the vibrational analysis at the Sβ1,TS transition state. Right part: a Newman projection for the C7dC8 double bond, showing the β and γ torsional angles and the definition of the HOOP coordinate.

“trapped” in the R-twisted domain. More specifically, the β-torsion, in the R-twisted domain, does not induce a back-twist of R. This implies that the Sβ1,TS transition state cannot be reached from the R-twisted configuration by following the minimum energy path (MEP) along β. This can be inferred from Figure 5 where we present a cut (“quasi MEP”) through the S1 PES along the β torsion, with all other coordinates being optimized. The corresponding S0 cross-section is also shown. In the figure, the abrupt energy decrease of S1 (energy increase of S0), at β j 173, reflects the R-twist of the chromophore. From β = 170 the S1 profile can be considered as the MEP along β in the R-twisted domain. As can be seen, further twist along β does not lead to back-twist of R (see values of R in Figure 4), with the chromophore exhibiting both the R and β twists at the last point (β = 90). At the latter point the S1 and S0 are energetically very close to each other (the energy difference is 0.3 eV). This suggests that a concerted R- and β-torsion constitutes a plausible coordinate leading to a conical intersection of the S1 and S0 PES. We now briefly summarize the above findings that revealed many peculiarities for the β-torsion of the pCTMe 3 2H2O chromophore. (i) The β-twist induces the R-torsion of the chromophore, leading eventually to the R-twisted structure; the chromophore is thus strikingly instable with respect to the R-twist. (ii) The induced R-torsion proceeds in a concerted way which is a result of the asymmetry of the S1 PES along the R-twist; no bifurcation of the reaction path is expected in this case. (iii) The β-torsion alone does not lead to the β-twisted configuration (Sβ1,min), for the system gets trapped in the R-twisted well. From the latter follows that the pathway toward the β-twisted minimum, that is, the first half of the transcis isomerization, should involve both the β- and R-twisting coordinates. Anticipating the discussion in the next section, there is also a third coordinate, the hydrogen out-of-plane (HOOP) distortion, that turns out to be crucial for the isomerization pathway of the pCTMe 3 2H2O chromophore. E. HOOP Coordinate and Its Interplay with r- and βTorsion. The vibrational analysis at the Sβ1,TS transition state revealed an interesting detail related to the transition mode, that is, the mode with imaginary frequency. The pattern of this mode is shown in Figure 6. As can be seen, it exhibits essentially displacements of the H7 and H8 hydrogens (displacements for the other atoms are at least 1 order of magnitude smaller, see Figure S1 in the Supporting Information). In fact, the relevant displacements are very similar to the pattern of the so-called

ARTICLE

hydrogen out-of-plane (HOOP) coordinate, identified in the photoinduced isomerization of retinal and related model systems.5760 To elucidate the implications of the HOOP coordinate, we calculated three 2D cuts through the S1 PES involving the two torsional coordinates, R and β, and the HOOP coordinate defined as the deviation of the H7C7C8H8 dihedral angle (γ) from the β dihedral angle, that is, HOOP = β γ (see Figure 6). The cuts are shown in Figure 7ac. Each cut represents the S1 PES as a function of two coordinates. The third coordinate was fixed at the value corresponding to the Sβ1,TS transition state, and all the remaining coordinates were allowed to relax. The value of HOOP at Sβ1,TS is 7.5. Let us first have a closer look at the cut with respect to the R and β twists (Figure 7a). The cut is quite instructive for the main topological features of the S1 PES discussed above. Indeed, one can see here the planar domain associated with the S1,min minimum (denoted P in Figure 7a), the R- and β-twisted R2 β domains associated with the SR1 1,minS1,min and S1,min points, β respectively, and the S1,TS transition state. It is readily seen that from the planar domain to the transition state, for 140 j β J 220 and R ≈ 180, the S1 surface exhibits a flat area (indicated by a red rectangle in Figure 7a) with a slight slope toward the R-twisted domains. This reflects the instability of the chromophore with respect to the R-torsion as discussed in section III. D.1. Obviously, the torsion along β from the planar domain will result in relaxation of the system to one of the R-twisted domains. This instability is seen to persist practically up to the Sβ1,TS transition state which corroborates our previous considerations in section III.D.1. Note that the slope is relatively small. From the dynamical point of view this suggests that the initial dynamics on the S1 PES should involve a significant spreading of the wavepacket over the region comprising the planar and R-twisted domains. The cut with respect to the R and HOOP coordinates (Figure 7b) illustrates that HOOP does play a role in the transcis isomerization of the pCTMe 3 2H2O chromophore. Depending on the value of HOOP, one can end up in the R-twisted or β-twisted domain, where the former is favored for HOOP < 7.5, whereas the latter occurs for HOOP > 7.5. The direction of relaxation from the Sβ1,TS transition state can be thus controlled by the HOOP coordinate. Note that there is a strong similarity between the cuts in Figures 7b and 7a. The origin of this similarity is an implicit dependence of both cuts on the γ coordinate, i.e. on the H7,C7, C8,H8 dihedral angle. More specifically, the cut in Figure 7b can be looked upon as a function of R and γ, because the variation in HOOP = β γ is now entirely determined by γ, at a fixed value of β. The cut in Figure 7a also reflects a change in γ: here, γ = β HOOP simply follows the variation of β, with a fixed value of HOOP. The above considerations suggest that the saddle-point topology seen in both Figure 7a and Figure 7b is essentially due to the γ coordinate. This is consistent with the pattern of the transition mode which is dominated by the latter coordinate (Figure 6). Because the transition state curvature is somewhat more pronounced in Figure 7a, the β coordinate also participates, as expected. One could thus conclude that the β and γ coordinates operate concertedly at the transition state, entailing the variation in the HOOP mode. The role of the HOOP mode is confirmed by the cut for the β and HOOP coordinates (Figure 7c, see in particular the contour 9244



ARTICLE

Figure 7. Perspective (3D) drawings and contour plots for the S1 PESs of pCTMe 3 2H2O (a) with respect to the R and β coordinates (HOOP is fixed at its value at Sβ1,TS, 7.5); (b) with respect to the R and HOOP coordinates (β is fixed at its value at Sβ1,TS, 133.1); and (c) with respect to the β and HOOP coordinates (R is fixed at its value at Sβ1,TS, 171.6). For all PESs shown, the remaining coordinates are relaxed. On the contour plots, the R, β, and planar (P) domains are indicated; the cross stands for the Sβ1,TS transition state. For all PESs shown, the vertical axis is energy (eV), with the CC2 ground-state energy at the S0,eq point taken as the origin. For all plots, R, β, and HOOP are in degrees.

plot). Here one can again readily see that, at the Sβ1,TS transition state (i.e., for β ≈ 133), the system can be directed by the HOOP mode either to the R-twisted domain, if HOOP < 7.5, or to the β-twisted domain, if HOOP > 7.5. In other words,

starting from a certain value of β, the HOOP indeed can control the course of the chromophore’s isomerization. In this way, “transitions” between the R and β domains, “steered” by the HOOP mode, obviously do not involve large skeletal distortions 9245



Figure 8. Minimum energy path (MEP) on the S1 PES of the isolated (pCTMe, red lines) and H-bonded chromophores, with two water molecules (pCTMe 3 2H2O, black lines) and four water molecules (pCTMe 3 4H2O, green lines), along the R torsional coordinate: full lines correspond to the results obtained with the aug-cc-pVDZ basis set and dashed lines are results obtained with the cc-pVDZ (without diffuse functions) basis set. The curves were aligned vertically with respect to the point R = 180 of the curve for pCTMe 3 2H2O (full line).

of carbon atoms and therefore should occur relatively fast. Note, however, that HOOP alone (i.e., without the β-torsion) does not lead to the β-twisted domain, that is a certain β-twist is required. From the above discussion follows that the photoinduced transcis conversion in pCTMe 3 2H2O operates via an interplay of three coordinates, β, HOOP, and R. It is conceivable that the β-coordinate first “prepares” the chromophore for the isomerization, whereas the HOOP mode accomplishes it when a certain excursion in β is attained. The HOOP mode can thus be understood to exert a local influence in the transition state region, while the β coordinate acts as a global isomerization coordinate. This scenario is in line with experimental observations for rhodopsin where the HOOP coordinate was found to come into play ∼70 fs after the beginning of the excited-state dynamics.57 The R coordinate apparently competes with the transcis isomerization, due to the barrierless relaxation to the R-twisted well. As discussed in the next section this “unfavorable” relaxation pathway is expected to be reduced or even eliminated by the stronger H-bonding protein environment. F. Impact of H-Bonding on the Chromophore’s Torsion. This section addresses a significant change of the R-twisting profile that is induced by the H-bonds between the chromophore and the water molecules. Our discussion is based on Figure 8 where we present results of the calculations for the MEP along the R coordinate, obtained for three different chromophore models. These are the isolated chromophore (pCTMe), the pCTMe 3 2H2O species already considered above, and the complex of the chromophore with four water molecules, pCTMe 3 4H2O (see Figure 1 for the structure). The latter species could in principle occur as a transient configuration in aqueous solution, though we employ it here essentially to demonstrate the impact of the number of H-bonds in a pronounced way. It is clearly seen from Figure 8 that the H-bonds have a profound influence on the S1 PES profile along R. Whereas for the isolated chromophore the planar stationary structure is a saddle point with respect to the R-torsion, it becomes a minimum (the S1,min point) in the case of the pCTMe 3 2H2O complex. Moreover the minimum becomes deeper in the case of pCTMe 3 4H2O that has two more hydrogen bonds. It is also

ARTICLE

Figure 9. Cross sections of the S1 PESs of the pCTMe 3 4H2O and pCTMe 3 2H2O chromophores along the β torsional coordinate, for β changing from 180 to 90 and all the remaining coordinates relaxed. The value of R (degree) is indicated for each point of the cross sections. The abrupt decrease of the S1 energy at β = 160 for pCTMe 3 4H2O and β = 173 for pCTMe 3 2H2O relates to the R-twist of the chromophores (see text for details).

expected that in the latter case there is a bound zero vibrational level, in contrast to the pCTMe 3 2H2O chromophore. The trend in Figure 8 is unambiguously related to the stabilization of the negative charge at the phenol ring by the H-bonds. More precisely, the H-bonding tends to “hold” the charge on the phenol ring, thereby preventing its migration toward the alkene part upon the R-twist. (Redistribution of the negative charge upon the R-twist was discussed in section III.B and can be clearly seen in Table 3.) In other words, the R-twist becomes less favorable in the presence of the H-bonds which manifests itself in the formation of the SR1,TS barrier (in the case of pCTMe 3 2H2O). This barrier further increases in the case of pCTMe 3 4H2O, as the H-bond stabilization effect enhances. This can also be looked upon as the formation and stabilization of the S1,min minimum observed for the pCTMe 3 2H2O and pCTMe 3 4H2O species, respectively. From the above discussion follows that the H-bonds can hinder the R-twist of the chromophore, with the effect increasing with a larger number of the bonds. The impact is expected to depend also on the strength of the bonds. This could lead to differences between H-bonding by water, as compared with the native protein environment where the chromophore is hydrogenbonded to the hydroxyl groups of the Tyr42 and Glu46 amino acids. The H-bonds in this case should be obviously stronger than in the pCTMe 3 4H2O complex, which means that the above effect will be even more pronounced. Our preliminary calculations where the two water molecules are substituted by Tyr42 and Glu46 amino acid fragments corroborate this conjecture. In the protein, the chromophore also exhibits an additional H-bond in the alkene moiety, involving the carbonyl oxygen (see Figure 1). This H-bonding might result in some reduction of the effect as the transfer of the charge toward the alkene part is obviously favorable for the carbonyl H-bond. One can, however, expect that the H-bonding at the phenolic oxygen will dominate as the H-bonds here are considered of a special type (so-called low-barrier H-bonds with a quasi-covalent bond feature)37 and are supposed to be stronger than the carbonyl H-bond, which is of a more usual type. Further investigations are required to confirm this conjecture. The most important consequence of the above effect is an alteration of the pathway and mechanism leading to the β-twisted 9246


The Journal of Physical Chemistry A configuration, that is, for the first half of the transcis isomerization pathway. Specifically, for the chromophore in the protein environment, one might conjecture that the β-twist proceeds “directly”, that is, without the relaxation to the R-twisted well, for the R-torsion is hindered due to the H-bonding effect. These conjecture is partly confirmed by the cut along the βcoordinate for the pCTMe 3 4H2O complex shown in Figure 9 (the cut for pCTMe 3 2H2O is also shown again for comparison). Even though the pCTMe 3 4H2O chromophore still undergoes the R-twist, the latter occurs at a larger excursion of the β-torsion (β ≈ 160), as compared with the pCTMe 3 2H2O species where it takes place at β ≈ 173. The trend is thus clear: at a certain number and strength of the H-bonds, the β-torsion becomes stable against the R-twist. To show this explicitly and to confirm our predictions, additional calculations for the chromophore in the protein environment are required. We reserve this study for a future work. Note that an effect similar to that just discussed was previously reported by Groenhof et. al for the Arg52 amino acid in native PYP.32 Due to its positive charge and proximity to the chromophore ring Arg52 was proposed to stabilize the negative charge at the phenolic ring thus favoring the isomerization around the double bond. As mentioned in the introduction this impact of Arg52 has, however, not been confirmed experimentally,19 and the charged state of Arg52 has been disputed.37 In a recent simulation of the S1S0 decay for a PYP chromophore derivative in water, Boggio-Pasqua et. al reported that hydrogen bonds between the chromophore and water can modulate the S1S0 energy difference at the single- and double-bond twisted configurations, thus facilitating the S1S0 decay.27 This effect also relates to stabilization of the chromophore’s charge at the phenolic ring (upon the double-bond twist) or at the alkene part (upon the single-bond twist). Our present findings on the H-bonding effect are thus in line with results of the previous studies. Lastly, it is worthwhile to emphasize that the above H-bonding effect appears only when the H-bonds are properly described. More specifically, it is crucial for the basis set to contain diffuse functions. This can be clearly seen from comparison of the MEPs in Figure 8 obtained with the aug-cc-pVDZ and cc-pVDZ bases. Apparently, the effect cannot be reproduced with the cc-pVDZ basis. The same is expected to hold true for other bases which do not contain diffuse functions.

IV. CONCLUSIONS AND OUTLOOK We have presented a detailed theoretical account of the structural changes of a model compound for the PYP chromophore (p-coumaric thiomethyl, hydrogen bonded with two water molecules, pCTMe 3 2H2O), following S0S1 photoexcitation. Various twisting/bending degrees of freedom, two torsional modes and a hydrogen out-of-plane bending mode have been considered and their subtle interplay been established. This interplay governs the isomerization process which is considered the very first step of the photocycle which occurs on a sub-ps time scale and triggers the photoreceptor’s signaling cascade. In the literature the primary isomerization process is often assumed to involve a CdC double bond. In our previous work, confirmed by the present more extended investigation, we have shown that in the isolated chromophore a competing singlebond isomerization occurs which is barrierless, or very nearly so. The two torsional modes are intimately coupled to each other and also to a third bending coordinate of HOOP type. We have

ARTICLE

demonstrated, for the first time, that the HOOP mode plays a key role in the PYP transcis isomerization process. Along with the double-bond torsional coordinate, which plays the major role in the beginning of the process, the HOOP is predicted to give a crucial “impulse” to overcome the isomerization barrier. The single-bond torsion, in contrast, was found to hinder the isomerization due to the barrierless relaxation of the chromophore to the single-bond twisted minimum (section III.D). The coupling between the single-bond and double-bond torsions was found to result in many peculiarities for the transcis pathway, which are summarized in the last paragraph of section III.D. Interestingly, the isomerization picture changes with an increasing number of H-bonds in the system: for the chromophore with 2H2O molecules mimicking H-bonding to neighboring amino acids in the native protein, a small barrier for single-bond torsion emerges, provided a proper computational scheme is adopted in the ab initio treatment. With 4H2O molecules and a corresponding increase in the number of H-bonds, this turns out to increase, which is expected to give rise to changes in the mechanism and pathway of the isomerization. In future work we plan to replace the H2O molecules by the tyrosine and glutamate amino acid residues that provide H-bonds to the chromophore in the native protein environment. Here, we expect further strengthening of the H-bonding effect. If confirmed, the dominance of double-bond torsion would be an important environment-induced effect in the chromophore’s biological function, although within a more complex picture than anticipated so far. Along with dynamical calculations, which are also planned in future work to corroborate the static picture of isomerization, this would demonstrate a subtle interplay between different torsional degrees of freedom and structural motifs (H-bonds) toward biomolecular function in PYP. Generalization to other photoactive proteins and photobiological processes in general would be an exciting perspective.

’ ASSOCIATED CONTENT

bS

Supporting Information. Total ground-state and excited-state energies, zero-point vibrational energies at the stationary points of pCTMe 3 2H2O; Cartesian coordinates of the stationary points; Cartesian coordinates of the ground-state equilibrium geometry of pCTMe 3 4H2O; The pattern of the transition mode at the Sβ1,TS stationary point showing the minor components of the transition vector. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT Financial support by the Deutsche Forschungsgemeinschaft (DFG) is gratefully acknowledged. We wish to thank the high performance computer centers bwGrid (member of the German D-Grid initiative, funded by the Ministry for Education and Research and the Ministry for Science, Research and Arts BadenW€urttemberg, http://www.bw-grid.de), Helics2 (Heidelberger Linux Cluster System, http://helics.uni-hd.de), and JUROPA (J€ulich Research on Petaflop Architectures, http://www.fz-juelich.de/ jsc/juropa) for generously providing computational resources. 9247



’ REFERENCES (1) Sgarbossa, A.; Checcucci, G.; Lenci, F. Photochem. Photobiol. Sci. 2002, 1, 459. (2) Photobiological Sciences Online, http://www.photobiology.info (accessed July 5, 2011). (3) Hummer, G. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 2381. (4) Meyer, T. E.; Tollin, G.; Hazzard, J. H.; Cusanovich, M. A. Biophys. J. 1989, 56, 559. (5) Sprenger, W. W.; Hoff, W. D.; Armitage, J. P.; Hellingwerf, K. J. J. Bacteriol. 1993, 175, 3096. (6) Getzoff, E. D.; Gutwin, K. N.; Genick, U. K. Nat. Struct. Biol. 2003, 10, 663. (7) Genick, U. K.; Soltis, S. M.; Kuhn, P.; Canestrelli, I. L.; Getzoff, E. D. Nature 1998, 392, 206. (8) Borgstahl, G. E. O.; Williams, D. R.; Getzoff, E. D. Biochemistry 1995, 34, 6278. (9) Baca, M.; Borgstahl, G. E. O.; Boissinot, M.; Burke, P. M.; Williams, D. R.; Slater, K. A.; Getzoff, E. D. Biochemistry 1994, 33, 14369. (10) Yoshizawa, T.; Kuwata, O. In CRC Handbook of Organic Photochemistry and Photobiology; Horspool, W. M., Song, P., Eds.; CRC Press: New York, 1995. (11) Larsen, D. S.; van Grondelle, R. ChemPhysChem 2005, 6, 828. (12) Changenet-Barret, P.; Espagne, A.; Plaza, P.; Hellingwerf, K. J.; Martin, M. M. New J. Chem. 2005, 29, 527. (13) Kyndt, J. A.; Meyer, T. E.; Cusanovich, M. A. Photochem. Photobiol. Sci. 2004, 3, 519. (14) Mataga, N.; Chosrowjan, H.; Taniguchi, S. J. Photochem. Photobiol., C 2004, 5, 155. (15) Hellingwerf, K. J.; Hendriks, J.; Gensch, T. J. Phys. Chem. A 2003, 107, 1082. (16) Cusanovich, M. A.; Meyer, T. E. Biochemistry 2003, 42, 4759. (17) Changenet-Barret, P.; Loukou, C.; Ley, C.; Lacombat, F.; Plaza, P.; Mallet, J.-M.; Martin, M. M. Phys. Chem. Chem. Phys. 2010, 12, 13715. (18) Changenet-Barret, P.; Plaza, P.; Martin, M. M.; Chosrowjan, H.; Taniguchi, S.; Mataga, N.; Imamoto, Y.; Kataoka, M. J. Phys. Chem. C 2009, 113, 11605. (19) Changenet-Barret, P.; Plaza, P.; Martin, M. M.; Chosrowjan, H.; Taniguchi, S.; Mataga, N.; Imamoto, Y.; Kataoka, M. Chem. Phys. Lett. 2007, 434, 320. (20) Wilderen, L.; Horst, M.; Stokkum, I.; Hellingwerf, K. J.; van Grondelle, R.; Groot, M. L. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 15050. (21) Lee, I.; Lee, W.; Zewail, A. H. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 258. (22) Heyne, K.; Mohammed, O. F.; Usman, A.; Dreyer, J.; Nibbering, E.; Cusanovich, M. A. J. Am. Chem. Soc. 2005, 127, 18100. (23) El-Gezawy, H.; Rettig, W.; Danel, A.; Jonusauskas, G. J. Phys. Chem. B 2005, 109, 18699. (24) Changenet-Barret, P.; Espagne, A.; Charier, S.; Baudin, J.-B.; Jullien, L.; Plaza, P.; Hellingwerf, K. J.; Martin, M. M. Photochem. Photobiol. Sci. 2004, 3, 823. (25) Wang, Y.; Li, H. J. Chem. Phys. 2010, 133, 34108. (26) Coto, P. B.; Roca-Sanjuan, D.; Serrano-Andres, L.; Martín-Pendas, A.; Martí, S.; Andres, J. J. Chem. Theory Comput. 2009, 5, 3032. (27) Boggio-Pasqua, M.; Robb, M. A.; Groenhof, G. J. Am. Chem. Soc. 2009, 131, 13580. (28) Virshup, A.; Punwong, C.; Pogorelov, T. V.; Lindquist, B. A.; Ko, C.; Martínez, T. J. J. Phys. Chem. B 2009, 113, 3280. (29) Groenhof, G.; Sch€afer, L. V.; Boggio-Pasqua, M.; Grubm€uller, H.; Robb, M. A. J. Am. Chem. Soc. 2008, 130, 3250. (30) Ko, C.; Virshup, A. M.; Martínez, T. J. Chem. Phys. Lett. 2008, 460, 272. (31) Gromov, E. V.; Burghardt, I.; Hynes, J. T.; K€oppel, H.; Cederbaum, L. S. J. Photochem. Photobiol., A 2007, 190, 241. (32) Groenhof, G.; Bouxin-Cademartory, M.; Hess, B.; de Visser, S. P.; Berendsen, H. J. C.; Olivucci, M.; Mark, A. E.; Robb, M. A. J. Am. Chem. Soc. 2004, 126, 4228.

ARTICLE

(33) Gromov, E. V.; Burghardt, I.; K€oppel, H.; Cederbaum, L. S. J. Am. Chem. Soc. 2007, 129, 6798. (34) Espagne, A.; Paik, D. H.; Changenet-Barret, P.; Plaza, P.; Martin, M. M.; Zewail, A. H. J. Photochem. Photobiol. Sci. 2007, 6, 780. (35) Espagne, A.; Changenet-Barret, P.; Baudin, J.-B.; Plaza, P.; Martin, M. M. J. Photochem. Photobiol., A 2007, 185, 242. (36) Espagne, A.; Changenet-Barret, P.; Baudin, J.-B.; Plaza, P.; Martin, M. M. In Femtochemistry VII: Fundamental Ultrafast Processes in Chemistry, Physics and Biology; Castleman, A. W., Jr., Kimble, M. L., Eds.; Elsevier: Amsterdam, 2006; p 204. (37) Yamaguchi, S.; Kamikubo, H.; Kurihara, K.; Kuroki, R.; Niimura, N.; amd Y. Yamazaki, N. S.; Kataoka, M. Proc. Natl. Acad. Sci. U.S.A. 2009, 109, 440. (38) Rajput, J.; Rahbek, D. B.; Aravind, G.; Andersen, L. H. Biophys. J. 2010, 98, 488. (39) de Groot, M.; Gromov, E. V.; K€oppel, H.; Buma, W. J. J. Phys. Chem. B 2008, 112, 4427. (40) Christiansen, O.; Koch, H.; Jørgensen, P. Chem. Phys. Lett. 1995, 243, 409. (41) To validate the results of the CC2 method a series of single point calculations for the S1 excitation energy of the pCTMe 3 2H2O chromophore were carried out using the CC2, EOM-CCSD, and CR-EOM-CCSD(T) coupled-cluster scheme and the cc-pVDZ basis set. The calculated S1 excitation energy amounts to 3.11, 3.28, and 3.02 eV at the CC2, EOM-CCSD, and CR-EOMCCSD(T) level, respectively. From this follows that the CC2 method predicts the S1 energy of pCTMe 3 2H2O to a good accuracy, whereas the EOMCCSD method gives a somewhat overestimated value. (42) Reed, A. E.; Weinstock, R. B.; Weinhold, F. J. Chem. Phys. 1985, 83, 735. (43) Aquilante, F.; de Vico, L.; Ferre, N.; Ghigo, G.; Malmqvist, P.-Å.; Neogrady, P.; Pedersen, T. B.; Pito nak, M.; Reiher, M.; Roos, B. O.; et al. J. Comput. Chem. 2010, 31, 224. (44) Eichkorn, K.; Treutler, O.; H., O.; H€aser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 240, 283. (45) Vahtras, O.; Alml€of, J.; Feyereisen, M. W. Chem. Phys. Lett. 1993, 213, 514. (46) Ahlrichs, R.; B€ar, M.; H€aser, M.; Horn, H.; K€olmel, C. Chem. Phys. Lett. 1989, 162, 165. (47) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007. (48) Krylov, A. I. Annu. Rev. Phys. Chem. 2008, 59, 433. (49) Bartlett, R. J. Int. J. Mol. Sci. 2002, 3, 579. (50) Gromov, E. V.; Trofimov, A. B.; Gatti, F.; K€oppel, H. J. Chem. Phys. 2010, 133, 164309. (51) Mathematica, version 6.0; software for technical computation; Wolfram Research: Champaign, IL, 2007. (52) K€ ohn, A.; H€attig, C. J. Chem. Phys. 2003, 119, 5021. (53) Coto, P. B.; Martí, S.; Oliva, M.; Olivucci, M.; Merchan, M.; Andres, J. J. Phys. Chem. B 2008, 112, 7153. (54) Yamada, A.; Yamamoto, S.; Yamato, T.; Kakitani, T. J. Mol. Struct.: THEOCHEM 2001, 536, 195. (55) Quapp, W.; Hirsch, M.; Heidrich, D. Theor. Chem. Acc. 2004, 112, 40. (56) Quapp, W. J. Mol. Struct. 2004, 695696, 95. (57) Kakitani, T.; Akiyama, R.; Hatano, Y.; Imamoto, Y.; Shichida, Y.; Verdegem, P.; Lugtenburg, J. J. Phys. Chem. B 1998, 102, 1334. (58) Kukura, P.; McCamant, D. W.; Yoon, S.; Wandschneider, D. B.; Mathies, R. A. Science 2005, 310, 1006. (59) Weingart, O. Chem. Phys. 2008, 349, 348. (60) Sumita, M.; Ryazantsev, M. N.; Saito, K. Phys. Chem. Chem. Phys. 2009, 11, 6406.

’ NOTE ADDED AFTER ASAP PUBLICATION This article posted ASAP on August 3, 2011. Paragraph six, sentence four of the Results and Discussion section and the caption of Figure 8 have been revised. The correct version posted on August 18, 2011. 9248


Photoinduced Isomerization of the Photoactive Yellow Protein (PYP

Recommend Documents