NMR Chemical Shifts of the Rhodopsin Chromophore in the Dark

J. Phys. Chem. B 2008, 112, 1267-1274

1267

NMR Chemical Shifts of the Rhodopsin Chromophore in the Dark State and in Bathorhodopsin: A Hybrid QM/MM Molecular Dynamics Study Ute F. Ro1 hrig† and Daniel Sebastiani*,‡ Ludwig Institute for Cancer Research and Swiss Institute of Bioinformatics, Molecular Modeling Group, Genopode Building CH-1015 Lausanne, Switzerland, and Max-Planck Institute for Polymer Research, Ackermannweg 10, 55128 Mainz, Germany ReceiVed: July 19, 2007; In Final Form: September 26, 2007

We investigate nuclear magnetic resonance (NMR) parameters of the rhodopsin chromophore in the dark state of the protein and in the early photointermediate bathorhodopsin via first-principles molecular dynamics simulations and NMR chemical shift calculations in a hybrid quantum/classical (QM/MM) framework. NMR parameters are particularly sensitive to structural properties and to the chemical environment, which allows us to address different questions about the retinal chromophore in situ. Our calculations show that both the 13C and the 1H NMR chemical shifts are rather insensitive to the protonation state of Glu181, an ionizable amino acid side chain located in the vicinity of the isomerizing 11-cis bond. Thus, other techniques should be better suited to establish its protonation state. The calculated chemical shifts for bathorhodopsin further support our previously published theoretical structure, which is in very good agreement with more recent X-ray data.

1. Introduction The photoreceptor rhodopsin (Rh) is located in the retina of vertebrates and serves for dim-light vision. The crystal structures of rhodopsin1-5 provide the first high-resolution structure of the large and pharmaceutically important class of G-protein coupled receptors. Moreover, they present a major step toward a better understanding of the first steps of the visual cascade. Vision is initiated by light absorption, which triggers 11-cis to all-trans isomerization of the Rh chromophore, the retinal protonated Schiff base (RPSB).6 The primary photoproduct, photorhodopsin, is formed on an ultrafast time scale (200 fs)7 and with high efficiency (quantum yield 0.65).8 Photorhodopsin thermally relaxes within a few picoseconds to bathorhodopsin, which displays a very distorted all-trans conformation.9,10 The photocycle proceeds via formation of a blueshifted intermediate, lumirhodopsin, and meta-rhodopsin I, before attaining metarhodopsin II, the active conformation for G-protein coupling and signal transduction.11 Whereas the basic features of the rhodopsin photocycle have been known for some time, the details of the structural changes between the dark state and the photointermediates as well as the mechanism of intramolecular signal transduction remain very active fields of research. We have previously investigated the photoreaction of rhodopsin by hybrid quantum mechanical/molecular mechanical (QM/ MM)12 excited-state molecular dynamics (MD) simulations of the full protein in a membrane-mimetic environment, treating the chromophore at the density functional theory (DFT) level.13 No restraints, neither on any structural part of the protein, nor on the reaction coordinate, were imposed. We have shown that * Corresponding author. E-mail: [email protected]. Fax: +49-6131-379-100. † Ludwig Institute for Cancer Research and Swiss Institute of Bioinformatics. ‡ Max-Planck Institute for Polymer Research.

the protein binding pocket selects and accelerates the isomerization exclusively around the C(11)-C(12) bond via preformation of a twisted structure. In ≈100 fs, a very twisted all-trans chromophore can be formed within the protein binding pocket by minor rearrangement of the nuclei and a clockwise isomerization sense. Our proposed bathorhodopsin structure was later confirmed by X-ray crystallography (Protein Data Bank accession code 2G87).14 Other models of bathorhodopsin have been developed on the basis of classical MD simulations,15 static excited-state QM/MM calculations,16 static ground-state QM/MM calculations,17 and resonance Raman experiments and calculations.18,19 It is difficult to compare the proposed structures, since their coordinates are usually not published except for selected internal coordinates such as some dihedral angles. However, an important parameter is the preservation of the Glu113-RPSB saltbridge, which has been observed experimentally.20 Experimental data for bathorhodopsin includes vibrational spectra as well as NMR data for selected carbon centers of the polyene chain.21 An interesting feature of the X-ray structures of Rh in the dark state1-5 is the presence of an ionizable side chain, Glu181, from the second extracellular loop close to the isomerizing chromophore bond. It has been suggested that Glu181 replaces Glu113 as a counterion for the RPSB at a late stage of the photocycle.22,23 However, there is an ongoing debate about the protonation state of this residue in the dark state.22-27 A protonated, neutral side chain is suggested by UV-visible spectroscopy with mutated proteins,22,25 pKa calculations,27 and two-photon spectroscopy.28 On the other hand, vibrational spectra fail to show an appropriate band of a protonated carboxylic acid,29,30 and molecular dynamics simulations24,31 as well as FTIR studies23 further support a deprotonated side chain. Nuclear magnetic resonance (NMR) experiments predict a negative charge close to the C(12)32,33 or the C(13) atom,34 which could be provided by Glu113 or by Glu181.

10.1021/jp075662q CCC: $40.75 © 2008 American Chemical Society Published on Web 01/05/2008

1268 J. Phys. Chem. B, Vol. 112, No. 4, 2008

Ro¨hrig and Sebastiani

Figure 1. Left: Numbering scheme for the RPSB. Right: Illustration of the QM subsystem (shown in ball-and-stick representation). The RPSB chromophore as well as a crystallographic water molecule and the sidechains of Lys296, Glu113, and Glu181 are described quantum mechanically.

Here, we present first-principles calculations of NMR parameters of the RPSB in rhodopsin. NMR is a widespread analysis tool in many areas of chemistry and biology. One of the key quantities in this context are NMR chemical shift spectra, which allow the characterization of the chemical environment of individual atoms. While the dominant factors which determine the location of NMR resonance lines are covalent bond parameters (hybridization state, bond lengths and angles), weaker interactions such as hydrogen bonds and other electrostatic effects are also known to have a sizable influence.35-38 The latter are often more important from a structural perspective, in particular, when the configuration of a molecule is already known. With the recent advances in computational methodology as well as computer hardware, the first-principles prediction of such noncovalent effects on experimentally observable spectra has come into reach for many systems of technological and fundamental scientific interest.39-42 Several methods exist to compute the influence of the chemical environment on NMR parameters within electronic structure approaches. The explicit incorporation of a large number of neighboring molecules is in principle most accurate but computationally very demanding and thus only applicable in simple cases41,43,44 Alternatively, one can resort to embedding schemes, which can treat the environment at various levels of approximation. Recently, the frozen-density embedding scheme within density functional theory has been extended to allow for the calculation of NMR shieldings.45 Polarizable continuum methods have become very popular, although they do not directly reflect the atomistic structure of the chemical environment.46-48 Hybrid QM/MM methods improve on this aspect, but in turn require more effort regarding the sampling of finite-temperature ensemble averages and the description of the QM/MM interface.49-52 Here, we present first-principles calculations of NMR parameters of the RPSB in rhodopsin. The goal of this work is twofold. First, we investigate the influence of the protonation state of Glu181 in the dark state upon the NMR spectra of the RPSB. Second, we compare the NMR spectra computed for our bathorhodopsin structure with the experimental data21 and with theoretical data for a different bathorhodopsin model. The main difference of the Gascon model17,53 with respect to ours is the observation of a stretched salt bridge between the RPSB and the Glu113, which results in a less favorable electrostatic energy. We have performed extensive hybrid QM/MM molecular dynamics simulations, and subsequently sampled the trajectories with QM/MM NMR chemical shift calculations. The RPSB as well as the Lys296, Glu181, and Glu113 sidechains and a crystallographic water molecule hydrogen bonded to Glu181 and Ser1862 are treated quantum mechanically, while the

remaining atoms were included via a molecular mechanics description (Figure 1). 2. Methods and Computational Details 2.1. Hybrid Quantum-Classical (QM/MM) Molecular Dynamics Simulations. We apply a QM/MM molecular dynamics technique54 based on the Car-Parrinello MD (CPMD)55-57 code. The QM subsystem is treated within density functional theory58-60 with the BLYP61,62 functional. We use standard norm-conserving pseudopotentials63-65 and a 70 Ry energy cutoff for the plane-wave expansion of the wave function. The inherent periodicity in the plane-wave calculations is circumvented solving Poisson’s equation for nonperiodic boundary conditions,66 while periodic boundary conditions are retained for the classical simulation box. The wave functions are evolved according to the Car-Parrinello MD scheme57,67 with a time step of 4 au (≈0.1 fs) and a fictitious electron mass of 400 au. The MM subsystem is described by the AMBER6/parm9968-70 force field; electrostatic interactions are taken into account by the P3M method,71 and all bonds involving hydrogen atoms are constrained by the SHAKE72 algorithm. Steric interactions between the QM and the MM part are modeled by a LennardJones potential as prescribed by the classical force field, whereas electrostatic interactions between the QM system and the classical point charges are taken into account by a Hamiltonian electrostatic coupling scheme.54,73 All simulations are performed at a temperature of 300 K. The rhodopsin simulations are based on a model of the protein in a membrane-mimetic environment that has been shown to be stable without restraints for several nanoseconds.24 It is based on the crystal structure of bovine rhodopsin (Protein Data Bank accession code 1HZX, chain A)2 and contains about 24 000 atoms. For details of the classical simulation setup, see refs 24 and 74. For the MD simulations, an orthorhombic quantum cell with edge lengths of 44 × 30 × 40 bohr3 is sufficient to converge the energies and geometries with respect to the cell size. The quantum cell contains the RPSB and the Lys296 side chain up to atom CD, as well as the side chains of Glu113 and Glu181 and a crystallographic water molecule (74 atoms, Figure 1). Where the boundary between the classical and the quantum system cuts through a covalent bond, the quantum system is saturated by a monovalent pseudopotential. 75 For each setup, 500 fs of simulation are carried out for equilibration and 4-5 ps for analysis. Three simulations were carried out: DARK181- is a simulation of Rh in the dark state Glu181 in the unprotonated, negatively charged form. DARK181H is a simulation of the same state, but with Glu181 in the protonated,

Hybrid QM/MM Molecular Dynamics Study

J. Phys. Chem. B, Vol. 112, No. 4, 2008 1269

Figure 2. Left: Computed and experimental34,93 13C NMR chemical shifts of the sp2-hybridized carbons of the rhodopsin chromophore in the dark C state. Right: Chemical shift differences ∆δC ) δcalc - δCexp for the different theoretical approaches.

neutral form. The earliest stable intermediate after photoisomerization (BATHO181- , all-trans RPSB), is obtained directly from a QM/MM excited-state isomerization of the chromophore in situ, with Glu181 bearing a negative charge.13 Here, data collection is started after 5 ps of equilibration. 2.2. QM/MM NMR Chemical Shift Calculations. As the basis for the spectroscopic calculations, we have sampled the QM/MM MD trajectories with NMR calculations of a series of snapshots, in the spirit of refs 43 and 44. For consistency, the computational setup in the calculation of the spectroscopic parameters was the same as in the QM/MM MD simulations; that is, the classical environment was represented in the ground state Hamiltonian (and its unperturbed orbitals) via classical point charges. However, a larger quantum cell (50 × 34 × 45 bohr 3) and the Tuckerman solver for the Poisson equation76 were used. The NMR chemical shifts were calculated via density functional perturbation theory77,78 for about 70 snapshots per system (DARK181-, DARK181H, BATHO181-). The set of computed nuclear shieldings was subsequently time-averaged. This averaging yields single numerical values compatible with experiment, where the same averaging process takes place, since the experimental measurement process has a finite duration (in the order of microseconds). This computational averaging is done for periods which are several orders of magnitude shorter than the averaging process in the corresponding experiments. Therefore, we always checked that a satisfying convergence was achieved in our phase space sampling protocol. The standard deviations obtained this way for the spectroscopic parameters characterize mainly the fluctuations in the underlying molecular geometries, which can be surprisingly large on the typical scale of NMR chemical shifts. This has been demonstrated recently on the instantaneous shift distributions in liquid water,44 where the fluctuations cover a range of about 10 ppm for the 1H resonances. Besides the averages of the chemical shifts over these trajectories, their statistical standard deviations ∆δi are shown in the form of error bars in all figures:

∆δi2 )

1

2 (δ(t) ∑ i - δi) N t

(1)

where δ(t) i represents the instantaneous NMR chemical shift value of atom i in a snapshot t, and δi ) 〈δ(t) i 〉. These error bars represent the intrinsic width of the natural instantaneous chemical shift distribution which is due to the thermal motion of the chromophore and its environment. Following the experimental convention, we quote chemical shifts relative to computed nuclear shieldings of standard

reference systems (tetramethylsilane, TMS):

δH(X) ) σH(TMS) - σH(X)

(2)

where σ represents the trace of the nuclear shielding tensor, X is the considered system, and an averaging of these shieldings over the set of sampled configurations is implied. For the bulk susceptibility correction, a sample with spherical shape was assumed.79,80 A special note has to be devoted to the use of pseudopotentials in combination with the nuclear shielding calculations. We use normconserving pseudopotentials to represent the ionic core and the inner electrons for non-hydrogen atoms. Therefore, our computed nuclear shieldings do not account for the effect of core electrons, which is a matter of recent discussions, especially for very heavy atoms.81,82 For first-row atoms, there is a variety of studies in which the performance of the pseudopotential approximation was generally found to be acceptable.39,79,80,83-87 However, it has turned out to be important to reference the sp2 carbon nuclear shieldings to TMS indirectly, that is, via a reference system of the same hybridization state. This is achieved via the identity: 2

δC/sp (X) ) σC(TMS) - σC(X) ) δCexp(C6H6) + 2

[σC(C6H6) - σC/sp (X)]calc (3) while for the sp3 carbons, σC (C2H6) is used. The values for δCexp(C6H6) and δCexp(C2H6) were taken from ref 88. To retain consistency with the finite-temperature NMR calculations for rhodopsin, the reference shieldings were obtained from ensemble averages at room temperature. While the absolute values of NMR chemical shifts for carbons of different hybridization states are typically not as accurate as in all-electron calculations,89 the relative trends for packing phenomena and hydrogen bonding are generally well reproduced and reliable.90-92 3. Results and Discussion 3.1. Carbon NMR Spectra in the Dark State. 3.1.1. 13C NMR of sp2 Carbons. The 13C NMR chemical shift values of the sp2-hybridized carbon atoms of the RPSB in the dark state of Rh are shown in Figure 2. Besides our own results for the protonated and an unprotonated Glu181 side chain, we also show the theoretical pattern obtained by Gascon et al.53 as well as the experimental values.34,93 At variance with the present work, in the earlier calculations53 only a single point calculation for an optimized structure was

1270 J. Phys. Chem. B, Vol. 112, No. 4, 2008

Ro¨hrig and Sebastiani

Figure 3. Left: Computed and experimental34,93 13C NMR chemical shifts of the sp3-hybridized carbons of the rhodopsin chromophore in the dark C state. Right: Chemical shift differences ∆δC ) δcalc - δCexp for the different theoretical approaches.

Figure 4. Left: Computed and experimental93 1H NMR chemical shifts of the RPSB in the dark state of Rh. Right: Chemical shift differences H ∆δH ) δcalc - δHexp for the different theoretical approaches.

carried out. Therefore, the authors could afford a larger quantum system (134 atoms) and an all-electron calculation with a hybrid functional (B3LYP). The Glu181 side chain was considered in its neutral, protonated form. In the present simulations, both in DARK181- and in DARK181H, the distance between the Schiff base nitrogen and the closer carboxylate oxygen of the Glu113 counterion amounts to 2.74 ( 0.11 Å. This value is in good agreement with other QM/MM studies,4,53 which generally yield shorter distances than the available X-ray structures (3.1-3.6 Å).1-5 The agreement of our calculated patterns with experiment is generally good and of similar quality as the earlier theoretical results. For the carbon atoms on the ionone ring side of the conjugated chain (C(5)...C(8)), our calculations tend to overestimate the chemical shifts, while for the carbon atoms on the Schiff base side of the conjugated chain (C(12)...C(15)), the opposite is the case. The zigzag pattern of the experimental shift values is somewhat underestimated by our calculations, which results in a similar zigzag pattern of the shift-difference plot. The shift for C(15) is strongly underestimated in both theoretical approaches, which could point to a deficiency of the underlying DFT structures. There is little difference in the spectra of the two models DARK181- and DARK181H. On the ionone-ring side, the values computed for DARK181- are in slightly better agreement with experiment, while for C(10) DARK181H yields a better result. Surprisingly, for C(12) and C(13), both models yield the same chemical shifts. This is startling as the Glu181 side chain is closest to these two sp2-hybridized carbon atoms of the polyene chain, and a negative charge close to these had been postulated based on NMR experiments.32-34 3.1.2. 13C NMR of sp3 Carbons. The chemical shifts of the sp3-hybridized carbon atoms are presented in Figure 3. Our

calculated patterns follow closely the experimental values.34,93 Concerning the deviations from experiment (Figure 3 right), a similar picture as for the sp2 carbons emerges. For C(1)...C(4), the same zigzag shape as found in the earlier calculations is observed, with their values being somewhat closer to experiment. Our results for the methyl groups (C(16)...C(20)) are somewhat less reliable, although the error with respect to experiment is almost constant. Comparing the models DARK181- and DARK181H, we find that the sp3 carbons yield very similar results. The strongest effect is seen on C(16) and C(17), but it is not obvious why these methyl groups, which are relatively far away from the Glu181 group, should feel the charge more than the other carbons. 3.2. Proton NMR Spectra in the Dark State. The NMR signals of all protons of the RPSB in the dark state of rhodopsin are presented in Figure 4. Generally, the hydrogen atoms bound to sp3 carbons are overshielded with respect to experiment, while the hydrogen atoms bound to sp2 carbons are undershielded. The deviations of our proton shifts with respect to experiment are shown in the right part of Figure 4. Consistently with the previous case, the ionone ring (C(1)...C(6)) and its environment exhibit errors of almost 2 ppm, which is surprisingly large. In contrast to this, the protons of the conjugated chain are well described within our QM/MM approach, with a similar accuracy as the setup of Gascon et al. Again, on all sites, the 1H NMR chemical shifts are quite insensitive to the protonation state of Glu181. 3.3. Carbon NMR Spectra in Bathorhodopsin. 3.3.1. 13C NMR of sp2 Carbons. Since our NMR spectra calculations in the dark state of rhodopsin do not yield any significant hint that Glu181 should be protonated, we suppose in line with our earlier work24 that it is more likely to be deprotonated. We therefore continue the comparison of the

Hybrid QM/MM Molecular Dynamics Study

J. Phys. Chem. B, Vol. 112, No. 4, 2008 1271

Figure 5. Left: Computed and experimental21 13C NMR chemical shifts of the sp2-hybridized carbons of the RPSB in bathorhodopsin. Right: C Chemical shift differences ∆δC ) δcalc - δCexp for the different theoretical approaches.

Figure 6. Left: Computed and experimental differences ∆δC ) δCbatho - δCdark of the sp2 carbons. Right: Double differences of the batho-shifts between experiment and calculations: ∆∆δC ) (δCbatho - δCdark)exp - (δCbatho - δCdark)calc.

Figure 7. Computed 13C NMR chemical shifts of the sp3-hybridized carbons of the RPSB in bathorhodopsin. For comparison, the corresponding shifts of the (unprotonated) dark state are also shown; experimental shifts are not available.

NMR spectra of the RPSB in rhodopsin and in bathorhodopsin on the basis of the unprotonated models DARK181- and BATHO181-. In BATHO181-, the saltbridge between the Schiff base nitrogen and the carboxylate oxygen of Glu113 is very stable, and its length (2.72 ( 0.12 Å) is very similar to the dark state (2.74 ( 0.11 Å), which is in good agreement with experimental data.20 Also the N-H-O angle (167 ( 6° in the dark state) remains practically unchanged (166 ( 8° in bathorhodopsin) and close to linear, unlike in the Gascon model. The only available NMR data for bathorhodopsin stems from low-temperature experiments with a selectively 13C-labeled retinal.21 Comparison of these experimental values with our results and those by Gascon et al.53 is shown in Figure 5. Our overall pattern is very similar to the experimental one. Especially for C(12) and C(14) our calculations are in better agreement with experiment than those of Gascon et al., while for C(8) and C(11) they are worse. These findings give further support to our theoretical bathorhodopsin model,13 despite the remaining statistical and systematic errors inherent in our computational

Figure 8. Computed 1H NMR chemical shifts of the RPSB in bathorhodopsin. For comparison, the corresponding values for the (unprotonated) dark state are added as well; experimental proton chemical shifts are not available for bathorhodopsin.

approach. Again, the chemical shift of C(15) is rather badly reproduced by both DFT-based theoretical approaches. 3.3.2. Isomerization Shifts of sp2 Carbons. Computed NMR frequencies from ab initio methods are subject to numerous errors which are due to limitations of the theoretical approaches and the approximations used therein. However, many of these errors cancel each other when considering changes in spectroscopic properties which occur when the chemical environment is altered. In this spirit, we have looked at NMR chemical shift differences between the DARK181- and the BATHO181states, which represent the isomerization shift of the chromophore. Complementary to the “absolute” shift pattern along the conjugated chain, the influence of conformational changes on that pattern might be more insensitive to computational parameters. The results are shown in Figure 6. Experimentally, these chemical shift differences are only known for selected carbon atoms.21 When comparing the results from the Gascon structure and our calculations for the unprotonated DARK181-, the agreement

1272 J. Phys. Chem. B, Vol. 112, No. 4, 2008

Ro¨hrig and Sebastiani

Figure 9. Computed 13C(sp2) and 1H NMR chemical shifts using different computational setups. The black line includes the effect of atomic motion (QM/MM NMR sampling of the QM/MM molecular dynamics trajectory), while the green and violet data represent QM/MM NMR and QM-only NMR calculations in the annealed conformation.

of our model with experiment is better for carbons C(12) and C(14), while the chemical shift trends of C(10) and C(11) are better reproduced in the calculations by Gascon et al. The deviations of 2...3 ppm for C(10) and C(11) might indicate that the sampling of structural parameters such as the bond length alternation or the dihedral angle distribution are not yet perfectly equilibrated within our QM/MM molecular dynamics trajectory. We find again the same zigzag shape for C(11)...C(15) as already observed in Figure 2. 3.3.3. 13C NMR of sp3 Carbons. For the aliphatic carbon atoms, there is no experimental NMR data for the batho state.21 A comparison of our values is only possible with respect to the work of Gascon et al.,53 which is shown in Figure 7. In addition, we replot the pattern from the dark configuration to illustrate the changes which occur upon isomerization. The agreement between the two computational approaches is again reasonable. However, the changes between the dark state and bathorhodopsin are not very large compared with the natural width of the chemical shift distribution as indicated by error bars in Figure 7. 3.4. Proton NMR Spectra in Bathorhodopsin. Also the proton NMR chemical shifts change somewhat when the chromophore undergoes isomerization, as illustrated in Figure 8. Only H(8) and H(11) (which are both located close to the isomerizing bond) show a sizable difference between the dark state and the batho state. While this is not an unexpected behavior, it confirms that our model systems are consistent in the sense that two different molecular dynamics runs lead to similar results except where the configurations actually differ. In this respect, Figure 8 also illustrates the intrinsic sensitivity of the backbone hydrogens: Most average proton NMR chemical shifts in bathorhodopsin are only a small fraction of a part per million different from those in the dark state. 3.5. Effect of the Protein Environment and Atomic Motion. We have investigated the impact of phase space sampling and the polarization effect of the protein framework on the calculated NMR chemical shift spectra. In Figure 9, we compare the shifts of the hydrogen and sp2 carbons from the most complete protocol used in this work (QM/MM NMR calculations of an ensemble of snapshots from a trajectory generated using QM/MM molecular dynamics simulations) to the corresponding spectra obtained from a simplified computational setup. We have annealed our rhodopsin chromophore from the end of the QM/MM MD run to T ) 0 K and computed the equilibrium NMR chemical shifts with and without the point charges representing the surrounding protein.

It turns out that the incorporation of the polarizing effect of the charges in the protein is very important for an adequate representation of the electronic structure of the chromophore. The violet line in Figure 9 visibly has the largest deviations from experiment. The effect is weaker but still present when comparing the spectrum in the equilibrium geometry (green lines) with the ensemble average obtained from sampling the molecular dynamics trajectory. It is actually suprising that the difference between DARK181- and DARK181H is somewhat more pronounced in the optimized geometry than in the ensemble average (see Figure 2). This might indicate that, despite our relatively long simulation times, the phase space could not yet be completely sampled by our MD simulations. Very similar pictures are obtained for annealed configurations of bathorhodopsin (not shown). 4. Conclusion We have presented a first-principles investigation of the NMR chemical shift spectra of the retinyl chromophore of rhodopsin in the dark state and in bathorhodopsin. We have performed extensive hybrid quantum-mechanical/molecular modeling (QM/ MM) molecular dynamics simulations, which were then sampled with QM/MM NMR chemical shift calculations. Our calculations are based on molecular dynamics simulations, in which the RPSB as well as the Lys296, Glu181, and Glu113 sidechains and a crystallographic water molecule are treated quantum-mechanically, while the remaining atoms are included via a molecular mechanics description. The 1H and 13C NMR results for both possible protonation states of Glu181 are reported for the dark state of rhodopsin. We observe some disagreement between calculated and experimental shifts in certain cases. One reason could be that, for atoms which are in direct contact with the classically described protein environment, present-day QM/MM descriptions are not yet sophisticated enough to reproduce subtle effects on the hydrogen’s orbitals. For instance, the orbital of the proton’s bonding electron must be orthogonal to neighboring orbitals, which can substantially modify the proton’s nuclear shielding. A similar trend had already been observed in previous work50-52 and is subject to ongoing attempts to improve the interactions in the QM/MM contact region. Another more fundamental origin of this problem could be the imperfect description of conjugated carbon-carbon bonds by present-day exchange-correlation functionals. It is known that the degree of bond length alternation can be significantly underestimated, for example, in polyacetylene.94,95 One may suspect that the zigzag shape of the error curves of the sp 2 -hybridized carbon atoms in the dark state as well as in

Hybrid QM/MM Molecular Dynamics Study bathorhodopsin (Figures 2 and 5) is due to an imperfect description of the bond length alternation (BLA) in the conjugated carbon chain by the applied DFT/GGA (BLYP) methodology. However, we have shown previously74 that the maximum bond length deviation in this region for the RPSB in vacuum optimized with BLYP and with B3LYP is only 0.009 Å (C(9)-C(10) bond). For the RPSB model system C5H6NH2+, it has been shown in turn that calculations at the B3LYP, CASPT2, and MP2 level yield a very similar BLA.96,97 Additionally, it should be noted that the apparent correlation between the BLA pattern and the 13C chemical shift pattern is not valid for the whole conjugated chain, since the experimental NMR pattern for the dark state has a point of inversion at C(8) (Figure 2), while the bond lengths are truly alternating over the whole chain (C(5)-C(15), data not shown). We therefore conclude that the zigzag shape of the error curves cannot easily be attributed to an incorrect description of the BLA by the BLYP functional. However, since several DFT studies of the RPSB in different environments using different functionals53,98,99 all show a systematic underestimation of the 13C chemical shift values of C(15), this probably points to a deficiency of commonly used exchange-correlation functionals94,95 The influence of the protonation state of Glu181 upon the NMR spectra of the RPSB seems to be very small. A similarly small influence of the presence or absence of this charge has been found before by investigating the visual spectra of Rh.100 The agreement between the experimental NMR data for bathorhodopsin and the calculated NMR spectra for our model13 is slightly better than the values calculated for the model of Gascon et al.17 However, the differences are too small and not systematic enough to separate the influence of the different structural data from the influence of the different computational approaches. In summary, despite some deviations, we could show that the QM/MM modeling of the NMR resonance lines of the RPSB yields a good overall agreement with experiment. Acknowledgment. D.S. acknowledges support from the Deutsche Forschungsgemeinschaft (DFG) under Grant SE 1008/ 2. References and Notes (1) Palczewski, K.; Kumasaka, T.; Hori, T.; Behnke, C. A.; Motoshima, H.; Fox, B. A.; Trong, I. L.; Teller, D. C.; Okada, T.; Stenkamp, R. E.; Yamamoto, M.; Miyano, M. Science 2000, 289, 739-745. (2) Teller, D. C.; Okada, T.; Behnke, C. A.; Palczewski, K.; Stenkamp, R. E. Biochemistry 2001, 40, 7761-7772. (3) Okada, T.; Fujiyoshi, Y.; Silow, M.; Navarro, J.; Landau, E. M.; Shichida, Y. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 5982-5987. (4) Okada, T.; Sugihara, M.; Bondar, A.-N.; Elstner, M.; Entel, P.; Buss, V. J. Mol. Biol. 2004, 342, 571-583. (5) Edwards, P. C.; Li, J.; Burghammer, M.; McDowell, J.; Villa, C.; Hargrave, P. A.; Schertler, G. F. X. J. Mol. Biol. 2004, 343, 1439-1450. (6) Yoshizawa, T.; Wald, G. Nature 1963, 197, 1279-1285. (7) Schoenlein, R. W.; Peteanu, L. A.; Mathies, R. A.; Shank, C. V. Science 1991, 254, 412-415. (8) Kim, J. E.; Taubner, M. J.; Mathies, R. A. Biochemistry 2001, 40, 13774-13778. (9) Palings, I.; van den Berg, E. M.; Lugtenburg, J.; Mathies, R. A. Biochemistry 1989, 28, 1498-1507. (10) Sekharan, S.; Sugihara, M.; Weingart, O.; Okada, T.; Buss, V. J. Am. Chem. Soc. 2007, 129, 1052-1054. (11) Kliger, D. S.; Lewis, J. W. Isr. J. Chem. 1995, 35, 289-307. (12) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227-249. (13) Ro¨hrig, U. F.; Guidoni, L.; Laio, A.; Frank, I.; Rothlisberger, U. J. Am. Chem. Soc. 2004, 126, 15328-15329. (14) Nakamichi, H.; Okada, T. Angew. Chem., Intl. Ed. 2006, 45, 42704273. (15) Ishiguro, M.; Hirano, T.; Oyama, Y. ChemBioChem 2003, 4, 228231.

J. Phys. Chem. B, Vol. 112, No. 4, 2008 1273 (16) Andruniow, T.; Ferre, N.; Olivucci, M. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 17908-17913. (17) Gascon, J. A.; Batista, V. S. Biophys. J. 2004, 87, 2931-2941. (18) Yan, E. C. Y.; Ganim, Z.; Kazmi, M. A.; Chang, B. S. W.; Sakmar, T. P.; Mathies, R. A. Biochemistry 2004, 43, 10867-10876. (19) Kukura, P.; McCamant, D. W.; Yoon, S.; Wandschneider, D. B.; Mathies, R. A. Science 2005, 310, 1006-1009. (20) Pan, D.; Ganim, Z.; Kim, J. E.; Verhoeven, M. A.; Lugtenburg, J.; Mathies, R. A. J. Am. Chem. Soc. 2002, 124, 4857-4864. (21) Smith, S. O.; Courtin, J.; de Groot, H.; Gebhard, R.; Lugtenburg, J. Biochemistry 1991, 30, 7409-7415. (22) Yan, E. C. Y.; Kazmi, M. A.; Ganim, Z.; Hou, J.-M.; Pan, D.; Chang, B. S. W.; Sakmar, T. P.; Mathies, R. A. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 9262-9267. (23) Lu¨deke, S.; Beck, M.; Yan, E. C. Y.; Sakmar, T. P.; Siebert, F.; Vogel, R. J. Mol. Biol. 2005, 353, 345-356. (24) Ro¨hrig, U. F.; Guidoni, L.; Rothlisberger, U. Biochemistry 2002, 41, 10799-10809. (25) Yan, E. C.; Kazmi, M. A.; De, S.; Chang, B. S.; Seibert, C.; Marin, E. P.; Mathies, R. A.; Sakmar, T. P. Biochemistry 2002, 41, 3620-3627. (26) Lewis, J. W.; Szundi, I.; Kazmi, M. A.; Sakmar, T. P.; Kliger, D. S. Biochemistry 2004, 43, 12614-12621. (27) Periole, X.; Ceruso, M. A.; Mehler, E. L. Biochemistry 2004, 43, 6858-6864. (28) Birge, R. R.; Murray, L. P.; Pierce, B. M.; Akita, H.; Balogh-Nair, V.; Findsen, L. A.; Nakanishi, K. Proc. Natl. Acad. Sci. U.S.A. 1985, 82, 4117-4121. (29) Nagata, T.; Terakita, A.; Kandori, H.; Shichida, Y.; Maeda, A. Biochemistry 1998, 37, 17216-17222. (30) Kandori, H.; Kinoshita, N.; Yamazak, Y.; Maeda, A.; Shichida, Y.; Needleman, R.; Lanyi, J. K.; Bizounok, M.; Herzfeld, J.; Raap, J.; Lugtenburg, J. Biochemistry 1999, 38, 9676-9683. (31) Martinez-Mayorga, K.; Pitman, M. C.; Grossfield, A.; Feller, S. E.; Brown, M. F. J. Am. Chem. Soc. 2006, 128, 16502-16503. (32) Mollevanger, L. C.; Kentgens, A. P.; Pardoen, J. A.; Courtin, J. M.; Veeman, W. S.; Lugtenburg, J.; de Grip, W. J. Eur. J. Biochem. 1987, 163, 9-14. (33) Han, M.; Smith, S. O. Biochemistry 1995, 34, 1425-1432. (34) Smith, S. O.; Palings, I.; Miley, M.; Courtin, J.; de Groot, H.; Lugtenburg, J.; Mathies, R.; Griffin, R. Biochemistry 1990, 29, 8158-8164. (35) Yates, J. R.; Pham, T. N.; Pickard, C. J.; Mauri, F.; Amado, A. M.; Gil, A. M.; Brown, S. P. J. Am. Chem. Soc. 2005, 127, 10216-10220. (36) Schulz-Dobrick, M.; Metzroth, T.; Spiess, H. W.; Gauss, J.; Schnell, I. ChemPhysChem, 2005, 6, 315-327. (37) Ochsenfeld, C.; Brown, S. P.; Schnell, I.; Gauss, J.; Spiess, H. W. J. Am. Chem. Soc. 2001, 123, 2597-2606. (38) Bu¨hl, M.; Kabrede, H.; Diss, R.; Wipff, G. J. Am. Chem. Soc. 2006, 128, 6357-6368. (39) Gervais, C.; Dupree, R.; Pike, K. J.; Bonhomme, C.; Profeta, M.; Pickard, C. J.; Mauri, F. J. Phys. Chem. A 2005, 109, 6960-6969. (40) Murakhtina, T.; Delle Site, L.; Sebastiani, D. ChemPhysChem 2006, 7, 1215- 1219. (41) Bu¨hl, M.; Grigoleit, S.; Kabrede, H.; Mauschick, F. T. Chem. Eur. J. 2006, 12, 477-488. (42) Sekharan, S.; Weingart, O.; Buss, V. Biophys. J. 2006, 91, L07L09. (43) Sebastiani, D.; Parrinello, M. ChemPhysChem 2002, 3, 675. (44) Murakhtina, T.; Heuft, J.; Meijer, J.-E.; Sebastiani, D. ChemPhysChem 2006, 7, 2578-2584. (45) Jacob, C. R.; Visscher, L. J. Chem. Phys. 2006, 125, 194104. (46) Amovilli, C.; Mennucci, B. J. Phys. Chem. B 1997, 101, 10511057. (47) Tomasi, J.; Mennucci, B.; Cance`s, E. J. Mol. Struct. (THEOCHEM) 1999, 464, 211-226. (48) Klein, R. A.; Mennucci, B.; Tomasi, J. J. Phys. Chem. A 2004, 108, 5851-5863. (49) Cui, Q. J. Chem. Phys. 2002, 117, 4720-4728. (50) Sebastiani, D.; Rothlisberger, U. J. Phys. Chem. B 2004, 108, 2807. (51) Cui, Q.; Karplus, M. J. Phys. Chem. B 2000, 104, 3721-3743. (52) Komin, S.; Gossens, C.; Tavernelli, I.; Ro¨thlisberger, U.; Sebastiani, D. J. Phys. Chem. B 2007, 111, 5225-5232. (53) Gascon, J. A.; Sproviero, E. M.; Batista, V. S. J. Chem. Theory Comput. 2005, 1, 674-685. (54) Laio, A.; VandeVondele, J.; Rothlisberger, U. J. Chem. Phys. 2002, 116, 6941-6947. (55) Hutter, J. CPMD, version 3.12; Copyright IBM Corp. and MPIFKF Stuttgart: 1990-2007; http://www.cpmd.org. (56) Hutter, J.; Curioni, A. ChemPhysChem 2005, 6, 1788-1793. (57) Marx, D.; Hutter, J. Ab initio molecular dynamics: Theory and implementation. In Modern Methods and Algorithms of Quantum Chemistry; Grotendorst, J., Ed.; Forschungszentrum Ju¨lich: 2000; pp301-449, Vol 1 of NIC Series. (58) Hohenberg, P.; Kohn, W. Phys. ReV. B 1964, 136, 864-871.

1274 J. Phys. Chem. B, Vol. 112, No. 4, 2008 (59) Kohn, W.; Sham, L. J. Phys. ReV. A 1965, 140, 1133-1138. (60) Jones, R. O.; Gunnarsson, O. ReV. Mod. Phys. 1989, 61, 689746. (61) Becke, A. D. Phys. ReV. A 1988, 38, 3098-3100. (62) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785789. (63) Troullier, N.; Martins, J. L. Phys. ReV. B 1991, 43, 1993. (64) Goedecker, S.; Teter, M.; Hutter, J. Phys. ReV. B 1996, 54, 1703. (65) Hartwigsen, C.; Goedecker, S.; Hutter, J. Phys. ReV. B 1998, 58, 3641. (66) Hockney, R. W. Methods Comput. Phys. 1970, 9, 135-211. (67) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471-2474. (68) Case, D. A.; Pearlman, D. A.; Caldwell, J. W.; Ill, T. E. C.; Wang, J.; Ross, W. S.; Simmerling, C.; Darden, T.; Merz, K. M.; Stanton, R. V.; Cheng, A.; Vincent, J. J.; Crowley, M.; Tsui, V.; Gohlke, H.; Radmer, R.; Duan, Y.; Pitera, J.; Massova, I.; Seibel, G. L.; Singh, U. C.; Weiner, P.; Kollman, P. A. Amber; University of California. (69) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179-5197. (70) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049-1074. (71) Huenenberger, P. H. J. Chem. Phys. 2000, 113, 10464-10476. (72) van Gunsteren, W.; Berendsen, H. Mol. Phys. 1977, 34, 13111327. (73) Laio, A.; VandeVondele, J.; Rothlisberger, U. J. Phys. Chem. B 2002, 106, 7300-7307. (74) Ro¨hrig, U. F.; Guidoni, L.; Rothlisberger, U. ChemPhysChem 2005, 6, 1836-1847. (75) von Lilienfeld, O. A.; Tavernelli, I.; Rothlisberger, U.; Sebastiani, D. J. Chem. Phys. 2005, 122, 14113. (76) Martyna, G.; Tuckerman, M. J. Chem. Phys. 1999, 110, 2810. (77) Putrino, A.; Sebastiani, D.; Parrinello, M. J. Chem. Phys. 2000, 113, 7102-7109.

Ro¨hrig and Sebastiani (78) Baroni, S.; de Gironcoli, S.; del Corso, A.; Giannozzi, P. ReV. Mod. Phys. 2001, 73, 515. (79) Mauri, F.; Pfrommer, B.; Louie, S. Phys. ReV. Lett. 1996, 77, 53005303. (80) Sebastiani, D.; Parrinello, M. J. Phys. Chem. A 2001, 105, 1951. (81) Schreckenbach, G.; Wolff, S. K.; Ziegler, T. J. Phys. Chem. A 2000, 104, 8244-8255. (82) Straka, M.; Kaupp, M. Chem. Phys. Lett. 2005, 311, 45-56. (83) Mauri, F.; Louie, S. Phys. ReV. Lett. 1996, 76, 4246-4249. (84) Mauri, F.; Pfrommer, B.; Louie, S. Phys. ReV. Lett. 1997, 79, 2340. (85) Yoon, Y.; Pfrommer, B.; Mauri, F.; Louie, S. Phys. ReV. Lett. 1998, 80, 3388. (86) Mauri, F.; Vast, N.; Pickard, C. J. Phys. ReV. Lett. 2001, 88, 085506. (87) Sebastiani, D.; Goward, G.; Schnell, I.; Spiess, H. W. J. Mol. Struct. (THEOCHEM) 2003, 625, 283-288. (88) Jameson, A.; Jameson, C. Chem. Phys. Lett. 1987, 134, 461. (89) Pickard, C. J.; Mauri, F. Phys. ReV. B 2001, 63, 245101. (90) Hoffmann, A.; Sebastiani, D.; Sugiono, E.; Kim, K. S.; Spiess, H. W.; Schnell, I. Chem. Phys. Lett. 2004, 388, 164-169. (91) Sebastiani, D. Int. J. Quantum Chem. 2005, 101, 849-853. (92) Sebastiani, D. Mod. Phys. Lett. B 2003, 17, 1301-1319. (93) Creemers, A. F.; Kiihne, S.; Bove-Geurts, P.; DeGrip, W. J.; Lutenburg, J.; de Groot, H. J. M. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 9101-9106. (94) Jacquemin, D.; Perpe`te, E. A.; Ciofini, I.; Adamo, C. Chem. Phys. Lett. 2005, 405, 376-381. (95) Gascon, J. A.; Sproviero, E. M.; Batista, V. S. Acc. Chem. Res. 2006, 39, 184-193. (96) Page, C. S.; Olivucci, M. J. Comput. Chem. 2003, 24, 298-309. (97) Fantacci, S.; Migani, A.; Olivucci, M. J. Phys. Chem. A 2004, 108, 1208-1213. (98) Buda, F.; Giannozzi, P.; Mauri, F. J. Phys. Chem. B 2000, 104, 9048-9053. (99) Touw, S. I. E.; de Groot, H. J. M.; Buda, F. J. Phys. Chem. B 2004, 108, 13560-13572. (100) Schreiber, M.; Buss, V.; Sugihara, M. J. Chem. Phys. 2003, 119, 12045-12048.