Proton Transfer in Bacteriorhodopsin - American Chemical Society

10678

J. Phys. Chem. B 2000, 104, 10678-10691

Proton Transfer in Bacteriorhodopsin: Structure, Excitation, IR Spectra, and Potential Energy Surface Analyses by an ab Initio QM/MM Method Shigehiko Hayashi† and Iwao Ohmine* Chemistry Department, Faculty of Science, Nagoya UniVersity, Furo-cho, Chikusa-ku, Nagoya, Japan 464-8602 ReceiVed: April 20, 2000; In Final Form: August 3, 2000

The structures, the optical properties, and the proton-transfer mechanism of bacteriorhodopsin (bR) are investigated by using an ab initio QM/MM calculation. The nature of the hydrogen bond network (HBN) formed in the retinal pocket is examined. The analysis is made for the “environmental” effect from the molecules involved in this HBN on the geometry and the excitation spectrum of retinal. It is found that IR spectra of the O-D stretching of the bound water molecules and the N-D stretching of the Schiff base retinal are very sensitive to the vibration-electron coupling among these molecules in the pocket. The hypothetical protontransfer reaction path in the bR ground state is investigated, and the “environmental effect” on it is analyzed. The present calculation reveals that electronic interactions such as charge transfer and polarization, not included in most empirical methods, are very important to describe these properties of bR, the structures of HBN, the IR activities, and the proton-transfer mechanism.

1. Introduction Bacteriorhodopsin (bR) is a transmembrane protein in the purple membrane of Halobacterium salinarium and functions as a light-driven proton pump.1-5 The protein consists of seven transmembrane helices, A through G, and contains a chromophore, all transretinal bound covalently to Lys216 via a protonated Schiff base linkage. Photoexcitation of the chromophore initiates a photocycle which involves several intermediate states labeled J, K, L, M, N, and O and returns finally to the initial bR state. During the photocycle, bR serves as a unidirectional proton transport across the membrane from the inside (cytoplasmic) to the outside (extracellular) surface. The photocycle is composed of structural changes and elementary proton-transfer processes, and the proton translocation across the membrane is accomplished by a sequence of proton-transfer processes between key groups: the retinal Schiff base, Asp85, Asp96, and Glu204.6-9 In the primary step of the photocycle, bR f J f K, the light absorption triggers an isomerization of retinal around the C13dC14 double bond. After the following thermal activation/relaxation process to the L-intermediate state, the protonated Schiff base releases its proton to the negatively charged counterion, the carboxylic group of Asp85, in the L f M step. The proton transport is completed by further processes of proton transfers and structural changes in the rest of photocycle. The energy relaxation process in the photocycle must cause the structural changes of the retinal and protein environment in a way to control the acidities (pKa) of the key groups at each step so that the unidirectional proton transport is achieved. Internal water molecules have been considered to play an especially key role in the proton transport.10-20 On the basis of the recent X-ray crystallography studies for the bR ground state,10,11 the K-intermediate,12 and the late M-intermediate,13 * To whom correspondence should be addressed. E-mail: ohmine@aqua. chem.nagoya-u.ac.jp. † Reseach Fellow of the Japan Society for the Promotion of Science. Present address: Beckman Institute, University of Illinois, UrbanaChampaign, Illinois. E-mail: [email protected].

an atomic model of the primary proton transfer from the Schiff base to Asp85 has been proposed. In the bR ground state, three water molecules form a hydrogen bond network (HBN) in the chromophore binding pocket, which consists of polar groups, the Schiff base, Asp85, Asp212, and Arg82. On the other hand, one and two of these water molecules are observed to be dislocated from the binding pocket in the K- and M-intermediates, respectively. This suggests that the HBN rearrangement associated with the water molecular dislocation, induced initially by the photoisomerization of retinal, is a source of the pKa control. This feature of the HBN rearrangement can be studied by Fourier transform infrared (FTIR) experiments14-16 detecting how the frequencies, the intensities, and the transition dipole moment angles of the IR bands from O-H and N-H stretching vibrations change during the photocycle. The dynamics of the initial retinal photoisomerization process itself has also been a subject of great interest.21-39 Detailed knowledge of it will provide insight into the following energy relaxation processes in the photocycle. The isomerization in bR proceeds to a vibrationally hot state, assigned as the J-intermediate, in a very short time (∼500 fs)21-23 with a high quantum efficiency (∼0.6),24-25 and yields a highly selective photoproduct (13-cis conformation). Time-resolved spectroscopy28-30 showed that the fluorescence lifetime is also temperature dependent. These facts indicate that a small potential barrier exists in the isomerization process. The origin of the small barrier has been explained by a three-state model developed on the theoretical5,34 and experimental26,27 grounds. The photoreaction proceeds to the electronic ground state of the 13-cis conformation. Quantum chemical calculations of the isomerization process in the gas phase36-39 revealed that the S1-S0 nonadiabatic transitions occur via a conical intersection accounting for the high efficiency of the isomerization process. The proton transport involves various chemical reaction processes, namely, the bond formation/breaking (proton transfers) and the excited-state dynamics (photoisomerization), which can be greatly affected by the protein environment. The binding pocket indeed consists of the both positively and negatively

10.1021/jp001508r CCC: $19.00 © 2000 American Chemical Society Published on Web 10/21/2000

Proton Transfer in Bacteriorhodopsin charged residues and thus produces the strong electronic field. To understand the proton-transfer mechanism, it is very important to use a reliable description for the electronic structures of the binding pocket, such as an ab initio molecular orbital method, and also to fully take into account the protein environmental effect based on a microscopic molecular model. In the present study, we perform an ab initio QM/MM calculation to investigate the structures, the optical properties, and the proton-transfer mechanism of bR. The QM/MM method has been recently widely applied to deal with chemical reactions in liquid phases and biological systems.40 A pioneer work applying the QM/MM method for bR was carried out by Logunov and Schulten to investigate thermal isomerization of retinal in the dark-adaptation process.41 This method treats a part of the system ( the QM segment), which is involved in electronic change in reaction process, quantum mechanically and the rest of the system (the MM segment), which provides steric and electrostatic environmental effect on the reaction process, by means of the molecular mechanics. We perform the QM/MM calculation for the bR system, treating the molecules in the binding pocket (the chromophore, its counterion residues, and water molecules) as the QM segment and the rest of protein as the MM segment. Although the QM/MM approach provides a very efficient way to explore chemical reaction processes in large biological systems, it hinges on a somewhat arbitrary linking of the QM and MM segments when the division between them is set across covalent bonds. We introduce here an effective charge operator42 based on the restrained electrostatic potential charge (RESP) method43 and construct the QM/ MM Hamiltonian, which provides a proper treatment of the interaction at the QM/MM boundary. The details of the theory and the computational method are given in section 2 and the Appendix. The geometry of the bR ground state is determined, and the structures of the chromophore retinal and HBN in the binding pocket are discussed in section 3. The protein environmental effect on the electronic structures of the chromophore in the ground and the excited states is analyzed in section 4. In section 5, the potential energy curve of the hypothetical proton transfers from the Schiff base to Asp85 in the bR ground state is analyzed, and the role of HBN rearrangement in controlling the proton-transfer process is examined. The vibrational structure of the HBN is examined by normal-mode analysis, and the character of the IR signals is analyzed in section 6. The conclusions are given in section 7. 2. Computational Details The initial protein structure of bacteriorhodopsin is constructed by using the coordinate of the X-ray crystallographic structure by Luecke et al.44 (PDB code: 1BRX). The missing structure of the EF loop in 1BRX is supplemented with the coordinate from the electron diffraction study by Grigorieff et al.45 (PDB code: 2BRD). Nine water molecules are located in the vicinities of the chromophore, Arg82, and Glu204, based on the coordinate revealed in the X-ray diffraction studies reported very recently10,11 (PDB codes: 1C3W and 1QHJ). Asp96, Asp115, and Glu204 are assumed to be protonated.46,47 In the present QM/MM method (see Appendix for the detail of the method), the QM calculation is used to describe a segment of the system (the QM segment), namely, the side chain of retinal-Lys216 (RET), those of Asp85 and Asp212, and three water molecules in the binding pocket (3H2O), whereas the MM calculation is used to describe the rest of the protein (the MM segment). The boundary dividing the QM and MM segments is set to be at the Cβ-Cγ covalent bond of Lys216 and at the

J. Phys. Chem. B, Vol. 104, No. 45, 2000 10679

Figure 1. Atomic numbering of the system used in the present work. The dashed lines indicate possible hydrogen bonds deduced from the X-ray crystallographic structure of Luecke et al.10

CR-Cβ covalent bonds of the asparatic acids. To describe the QM segment properly, the linked atom approach (see Appendix) is used. Hydrogen atoms are used as “dummy atoms” in this approach. We employ two basis sets, the 3-21G basis set48 and the Dunning’s split-valence (DZV) basis set,49 for the QM calculations. In addition, the diffuse functions49,50 are added on the Oδ atoms of the asparatic acids (3-21G(O+) and DZV(O+)). The total number of the basis functions used in the QM calculation for the QM segment consisting of RET, Asp85, Asp212, and 3H2O is, for example, 427 in the DZV(O+) basis set. The coefficient of restraint terms in the RESP method are set to be 10-4 au for the atoms near the boundaries and 10-5 au for the other atoms. Although the RESP charges may depend on the restraint coefficients,43 the classical dipole moment derived from the RESP charges with the present restraint coefficients (e.g., 23.19 D for RET in the Hartree-Fock level with 3-21G) reproduce the quantum dipole moment well, i.e., the expectation value of dipole operator (23.15 D for the same system). The AMBER force field51 is used in the MM calculation to describe the rest of the protein and TIP3P52 is used for water molecules. The AMBER force field parameter is also used to express the van der Waals parameters in the QM segment describing the LJ interaction with the surrounding MM segment. For RET, we employ van der Waals parameters of the atom types of AMBER, CT, CA, HC, HA, N2, and H for the sp3 carbons, the sp2 carbons, hydrogen attached to CT, hydrogen attached to CA, nitrogen in the Schiff base, and hydrogen attached to N2, respectively. The classical Coulombic interaction is cut off at 12 Å. The QM/MM code is implemented in the HONDO53 and GAMESS54 program packages. The numbering of atoms in the retinal pocket used throughout this paper is defined in Figure 1. Water molecules in Figure 1 labeled WA, WB, and WC correspond to the water molecules termed W402, W401, and W406 in the structure of Luecke et al.,10 respectively. 3. Geometry Optimization The geometry of bR is optimized, starting from the initial structure (see section 2), by using the Hartree-Fock (HF) QM/ MM method in the DZV(O+) basis set with the QM segment consisting of RET, Asp85, Asp212, and 3H2O (denoted as RET,-

10680 J. Phys. Chem. B, Vol. 104, No. 45, 2000

Hayashi and Ohmine

Figure 3. Hydrogen bond network around the chromophore binding site. The values indicate atomic distances in the optimized structure of the RET,Asp85,Asp212,3H2O/DZV(O+)QM/MM system (model A) (upper) and those in the X-ray crystallographic structure (1C3W) of Luecke et al.10 (lower in the parentheses).

TABLE 1: Dihedral Angles (deg) of the Polyene Chain of Retinal at the C13‚‚‚N16 Regiona C14-C15-N16-Ce H15-C15dN16-H16 C13dC14-C15dN16 H14-C14-C15-H15 C12-C13dC14-C15 C13-C13dC14-H14

-166.1 -176.9 -179.0 -165.9 -160.4 -166.1

a Results of the model A (RET,Asp85,Asp212,3H O/DZV(O+)) 2 system.

Figure 2. Optimized structure for the RET,Asp85,Asp212,3H2O/DZV(O+) QM/MM system (model A): (a) a view from the Schiff base side to the direction of the β-ionone ring side of the chromophore and (b) a view from the Asp85 side to the direction of the Asp212 side.

Asp85,Asp212,3H2O/DVZ(O+); model A). The optimized bR structure around the retinal binding site is shown in Figure 2. The retinal Schiff base, Asp85, Asp212, Arg82, and three water molecules form a hydrogen bond network (HBN), which is shown in Figure 3. Even in the more restricted 3-21G(O+) basis set (model B), almost the same HBN structure is obtained. We can see from the figure that this HBN structure is very similar to that in the original 1C3W structure. There exist, however, certain differences in the HBN structure between 1C3W and the structure refined by the optimization. In the 1C3W structure, WA is hydrogen bonded to the Oδ1 atom of Asp212. In the geometry-optimized structure, the N-H bond of the Schiff base makes a strong hydrogen bond with the lone-pair electron of the WA water molecule, which fixed the orientation of the WA water molecule so that an O-H bond of WA forms a hydrogen bond with atom Oδ2 of Asp212, which is further hydrogen bonded to Tyr185. The O(WA)-Oδ1,2(D212) distance in the refined structure is, therefore, very different from that in 1C3W, as seen in Figure 3. Hatanaka et al.15,56 have proposed a model of the HBN structure based on Fourier transform infrared spectroscopy

experiments. They predicted that the water molecule is bonded to the Schiff base and has a tetrahedral structure forming hydrogen bonds with the Schiff base, Asp85, Asp212, and N-H of the indole of Trp86. Although WA in our model is located near the position where they suggested, WA is hydrogen bonded only to the Schiff base, Asp85, and Asp212. We cannot find a stable structure forming a hydrogen bond between the water molecule and N-H of Trp86. The conjugated polyene chain of retinal is all trans but twisted and bent around C13‚‚‚N16 bonds, and the N-H bond of the Schiff base is tilted toward Asp212. The torsion around the C15d N16 double bond can be specified by two dihedral angles, C14C15dN16-C and H15-C15dN16-H16. We can see in Table 1 that the dihedral angle of C14-C15dN16-C is -166.1°, considerably deviating from the planar structure (-180.0°), as observed in a molecular dynamics study by Tajkhorshid et al.55 On the other hand, the H15-C15dN16-H16 dihedral angle (-176.9°) is closer to a planar structure. Similarly, the torsion around the C14-C15 single bond is characterized by the nearly planar C13dC14-C15dN16 angle (-179.0°; Table 1)55 and the definitely twisted H14-C14-C15-H15 angle (-165.9°). This shows that the polyene is wrenched, twisted, and bent simultaneously along the C14-C15 and C15dN16 bonds (see Figure 2). The rotation around the C13dC14 double bond is characterized by the rather large deviations of both dihedral angles C12C13dC14-C15 and C13-C13dC14-H14, indicating that C13dC14 involves a strong torsion with a small bend of the polyene plane. The twisted and bent conformation of the polyene at the C13‚ ‚‚N16 bonds arises from the conformational restriction at the

Proton Transfer in Bacteriorhodopsin


Figure 4. Optimized structure for the RET/3-21G QM/MM system (model C). The view is from the Asp85 side in the direction of Asp212 side.

TABLE 2: Angles (deg) Related to the Hydrogen Bonding between the Schiff Base and WA

θ[NH‚OH]g θ[NH‚H1O]h θ[NH‚H2O]i

model Aa

model Bb

model Cc

model Dd

model Ee

ideal Tdf

168.6 118.3 106.5

170.9 116.3 105.0

155.8 149.8 95.9

159.6 144.4 90.4

169.6 123.1 105.7

180.0 109.5 109.5

a RET,Asp85,Asp212,3H O/DZV(O+). b RET,Asp85,Asp212,3H O/ 2 2 3-21G(O+). c RET/3-21G(O+). d RET,3H2O/3-21G(O+). e RET,Asp85, 3H2O/3-21G(O+). f The case of ideal tetrahedral coordination of the N16-H16 bond to the WA water molecule. g Angle between N16-H16 and O(WA)-H16. h Angle between N16-H16 and H1(WA)-O(WA). i Angle between N -H 16 16 and H2(WA)-O(WA).

retinal terminal from the Lys216 covalently bonding to retinal and the strong hydrogen bond between the Schiff base and the WA water molecule, and the steric restriction of the β-ionone ring and the head region of polyene from the surrounding proteins. The N-H bond of the Schiff base is tilted toward Asp212 to make the strong hydrogen bond with WA; the polyene is thus twisted counterclockwise around the axis directing from Schiff base to 3-ionone ring. Under the constraint caused upon the photoexcitation during the initial photoisomerization of retinal to the 13-cis configuration, the polyene may be further twisted counterclockwise; it is noted that the polyene excited state is twisted at its minimum. This may occur if the isomerization is nearly barrierless as suggested by the ultrafast isomerization rate (∼500 fs) observed experimentally. To find the key electronic interactions that lead to the correct molecular geometry and the right hydrogen bond network structure in bR, we carry out the QM/MM geometry optimizations by changing the QM segment sizes. The optimized geometry of the system, when only RET is treated as the QM segment (RET/3-21G: model C), is shown in Figure 4. We can see that the hydrogen bond between the N-H of the Schiff base and the lone pair of the WA water molecule is missing and Oδ1 of Asp212 is the partner for the hydrogen bond of WA as in the 1C3W structure, although the other part of the hydrogen bond network structure is similar to that of model A. Table 2 lists the angles related to the hydrogen bonding between the Schiff base and WA, namely, the angle between the directions of N16-H16 and O(WA)-H16, θ[NH‚OH], the angle between N16-H116 and H1(WA)-O(WA), θ[NH‚H1O],

and the angle between N16-H16 and H2(WA)-O(WA), θ[NH‚ H2O]. The N16-H16 bond forms an ideal hydrogen bond with the lone-pair orbital of the O atom of WA when the N16-H16 bond and the WA molecule are in a tetrahedral arrangement; in that case, θ[NH‚OH] should be 180.0° and both θ[NH‚H1O] and θ[NH‚H2O] should be 109.5°. In the optimized geometry of the model A system, θ[NH‚OH], θ[NH‚H1O], and θ[NH‚ H2O] are found to ensure the values 168.6°, 118.3°, and 106.50, respectively; they are within 15° of the angles of the ideal tetrahedral structure, indicating that a nearly ideal hydrogen bond is formed between the N16-H16 bond and the WA molecule in bR. On the other hand, the angles in model C deviate by more than 20° from the tetrahedral values. The respective angles are also calculated for the QM/MM systems with RET,3H2O/321G (model D) and with RET,Asp85,3H2O/3-21G(O+) (model E). The optimized angles for these models are listed in Table 2. Although the three water molecules and RET are treated as the QM segment in model D, a well-oriented hydrogen bond between the N-H bond and the lone pair orbital of WA is not yet formed and a partner of the hydrogen bond of the WA molecule is still Oδ1 of Asp212 as in model C. In model E, further including Asp85 in the QM segment, the angles approach those of model A. The strong hydrogen bond is now established between the N-H bond and WA, yielding the tetrahedral structure, and then Oδ2 of Asp212, instead of Oδ1, becomes the partner of hydrogen bond of WA as in model A. This suggests that the charge-transfer interactions among the Schiff base, WA, and Asp85 play a key role in establishing the distinct hydrogen bond network in this very polar moiety; the charge transfer interactions are not accounted for if the molecules are not fully treated by the QM method. The importance of this chargetransfer interaction will be seen also in section 6 for the IR spectrum of the O-H and N-H vibrations. 4. Effect of Protein Environment on Chromophore The energetics and electronic properties of the chromophore retinal in bR are analyzed. We examine the solvation of the ground electronic state of RET, which is described by HF wave function with DZV basis function, by first employing the model C type QM/MM division. In this division, treating RET as the quantum segment and its environment as the MM segment, the chromophore-environment interaction is described with the atomic-site-wise Coulombic (electrostatic) and Lennard-Jones terms so that the environment effect can be analyzed simply in terms of individual molecular contribution. For consistency, we use the model C (RET/3-21G) geometry for the analysis; the analysis should be made in principle for a geometry optimized by a same- or lower-level calculation (see below). We then analyze the quantum mechanical effects of the environment on the RET solvation by employing the model A (RET,Asp85,Asp212,3H2O/DZV(O+)) and its optimized geometry. In this level of calculation, the interactions among the molecules included in the quantum segment cannot described by the pairwise interactions. We here use the energy decomposition scheme developed by Kitaura and Morokuma59-61 to examine, particularly, the electronic interaction of the WA water molecule. Excited electronic states in bR are also calculated by the stateaveraged complete active space self-consistent field (CASSCF) and multireference Møller-Plesset (MRMP) perturbation methods for the model C system. Unfortunately, the calculation of the excited electronic states for the model A system is computationally too demanding. To assess the effect of the quantum electronic interaction on the excitation energies, therefore, we carry out the CASSCF calculation for the model


Hayashi and Ohmine TABLE 3: Solvation Energies of RET from Individual Protein Residues (kcal/mol) model C/model Ca Arg82 Asp85 Arg134 Glu194 Asp212 WA

Figure 5. Coulombic interactions of RET with individual protein residues in themodel C type division QM/MM system at the model C optimized geometry. The HF level with DZV basis functions is included.

E system where the hydrogen bonds character in WA is similar to that in the model A system, as shown in the previous section. 4.1. The Solvation Energy of the Ground State. By employing the model C type QM/MM division, the solvation energy of RET in bR (Esolv) is expressed as

Esolv ) VES + VLJ + Vbond + EERO

(1)

where VES and VLJ are electrostatic (Coulombic) Lennard-Jones interaction energies between RET and the protein environment, respectively. VES includes a gain of the electrostatic stabilization by polarization of electronic wave function induced by the electrostatic field of the protein environment. On the other hand, the electronic reorganization energy, EERO, represents the destabilization of the gas-phase energy component of the QM segment (expressed by Eele 0 in eq A-4) caused by the polarization of the electronic wave function.57 Vbond is the intramolecular bonding contribution of the MM force field in the QM/MM boundary region of retinal-Lys216. The solvation energy, Esolv, of the ground-state RET is calculated to be -162.7 kcal/mol at the HF level with a DZV basis in the model C geometry; RET is strongly stabilized in bR. The electrostatic interaction, VES ) -102.5 kcal/mol, is found to be the most dominant contribution to the solvation energy. (Note that the HF values are different from those CASSCF and MRMP values listed in Table 5.) Figure 5 demonstrates how individual protein residues in bR contribute to VES in the ground state. We can see in the figure that five charged residues give large Coulombic interactions. The negatively charged residues in the vicinity of the Schiff base, namely, Asp85 and Asp212, strongly stabilize RET. Arg82, Glu194, and Arg134 also give large Coulombic contribution. The solvation energies from those residues are listed in Table 3. Glu194 locates near Glu204 and participates in the HBN of the extracellular region of the proton channel as well as in the HBN of Glu204, Arg82, and water molecules. Arg134 locates near Glu194, although it does not directly face the proton channel. The interaction between RET and those residues may play a certain role in gating the extracellular channel in the proton translocation. The WA molecule, which is directly hydrogen bonded to the Schiff base, also yields a large Coulombic interaction energy, -15.8 kcal/mol. It is noted that the dipole moment of RET in bR, 23.1 D, is larger than that in the gas phase, 18.3 D (in Hartree-Fock with DZV basis calculation). The electronic ground state is significantly polarized under the strong electrostatic field of the protein environment. Mulliken population analysis indeed reveals that the strong anionic environment around the iminium cation of the chromophore induces a strong polarization of the electrons,

model C/model Ab

totalc

Coulombd

LJe

totalc

Coulombd

LJe

28.3 -71.4 16.0 -21.6 -65.9 -13.6

28.6 -69.7 16.1 -21.4 -63.0 -15.8

-0.3 -1.7 -0.1 -0.2 -2.9 2.2

27.4 -69.4 15.8 -20.9 -61.9 -10.1

27.7 -67.8 15.9 -20.7 -59.5 -18.0

-0.3 -1.6 -0.1 -0.2 -2.4 7.9

a Solvation energies evaluated for the model C type QM/MM division system at the model C optimized geometry. DZV basis functions are used. b Solvation energies evaluated for the model C type QM/MM division system at the model A optimized geometry. DZV basis functions are used. c Total solvation energy is given as a sum of Coulombic and Lennard-Jones contributions. d Contribution of the Coulombic interaction. e Contribution of the Lennard-Jones interaction.

thus reducing the delocalization of the positive charge along the polyene chain and destabilizing the resonance structure of the protonated Schiff base retinal.58 The non-negligible electronic reorganization energy by the electrostatic field, EERO ) 8.9 kcal/mol (in Hartree-Fock with the DZV basis calculation), is thus found in the ground state. The Lennard-Jones interaction energy, VLJ, is -69.4 kcal/mol. The bulky residues which constitute the retinal pocket are found to give a large contribution to VLJ. Vbond is only 0.3 kcal/mol. As seen in section 3, the quantum electronic interaction plays a crucial role in formation of the strong hydrogen bonds around RET; the HBN structure, especially that involved in WA, changes (to the model A geometry) when the quantum effect is accounted for the RET environment. We calculate the RET solvation energy in this model A geometry by using the model C type QM/MM division, which is listed in Table 3. We can see in the table that the RET solvation contributions from individual protein residues in the model A geometry are almost identical to those in the model C geometry. One may expect that the model A geometry calculation provides a more accurate description of a strong hydrogen bond between the Schiff base and WA than the model C geometry. We can see in Table 3, however, that the Schiff base and the WA interaction in the model A geometry (-10.1 kcal/mol) is weaker than that in the model C geometry (-13.6 kcal/mol), although their distance is shortened in the model A geometry. The Coulombic interaction between RET(QM) and WA(TIP3P) at the model A geometry (-18.0 kcal/mol) is stronger than that at the model C geometry (-15.8 kcal/mol), but the L-J repulsion energy at the former geometry (7.9 kcal/mol) is much higher than that at the latter geometry (2.2 kcal/mol). This shows that using the TIP3P model, which lacks the ability to describe the quantum electronic interaction, for WA gives an inconsistent result. 4.2. Effect of Quantum Electronic Interaction in HBN. To analyze the interactions among RET and its neighbor for the model A geometry properly, we must use the model A QM/ MM division, treating RET,Asp85,Asp212,3H2O quantum mechanically, and employ an energy decomposition scheme. Here we use the Kitaura-Morokuma (KM) energy decomposition analysis to deal with the interaction between the WA water molecule and the other RET neighbors. Two monomer molecules are defined in the model A system, i.e., the WA water molecule and a supermolecule including the other part of the QM segment (RET, Asp85, Asp212, WB, WC), and the interaction between them in the protein environment is analyzed. To reduce basis set superposition error (BSSE), we add polarization basis functions49 on WA and the N-H group of the Schiff base.



TABLE 4: The KM Analysis of Interaction Energy between the WA Water Molecule and the Rest of the QM Segment (RET, Asp85, Asp212, WB, WC) in the Model A System (kcal/mol)a ∆E (WA)b

EESc

EEXd

EPLe

ECTf

TABLE 5: Excitation Energies, VESa, EEROb (kcal/mol), and Dipole Moments of the Chromophore (D)c excitation energyd gas

EMIXg

-31.2 (-27.4) -64.7 61.5 (62.8) -14.1 -19.8 (-18.2) 5.9 (7.9) a DZV(O+) and additional polarization basis functions are used (see text). Values in parentheses are those with the BSSE correction. b Total quantum electronic interaction energy. c Electrostatic energy. d Exchange repulsion energy. e Polarization energy. f Charge-transfer energy. g High-order coupling energy.

The result of the energy decomposition analysis is listed in Table 4. We can see that the total quantum electronic interaction energy (∆E(WA)) is large, -31.2 kcal/mol (-27.4 kcal/mol with the BSSE correction). Although the absolute value of the BSSE correction (3.8 kcal/mol) is not small, it is only 12% of the total interaction energy. The present decomposition analysis is thus expected to provide at least a qualitative picture of the interaction. In the KM scheme, the total interaction energy is decomposed into several components: electrostatic (EES), exchange repulsion (EEX), polarization (EPL), charge transfer (ECT), and high-order coupling (EMIX) energies. In the present system, a considerably large electrostatic stabilization (EES ) -64.7 kcal/mol) is found. In addition, the large stabilization arises from the polarization (EPL ) -14.1 kcal/mol) and charge transfer (ECT ) -19.8 kcal/mol) interactions. These large interactions are due to the strong polar nature of the charged electronic states of the Schiff base, Asp85, and Asp212. On the other hand, EEX gives a large destabilization contribution (61.5 kcal/mol). Asp85 and Asp212 are negatively charged and thus have diffusive anion orbitals. As there exists the positively charged group of the Schiff base, the anionic molecular orbitals of Asp85 and Asp212 expand to the Schiff base. The WA molecule lies between those groups, and the occupied orbitals of WA overlap with the anionic molecular orbitals of Asp85 and Asp212; thus, it gives the large exchange repulsion energy. It is noteworthy that EES and EEX nearly cancel out mutually at the stable geometry and the Heitler-London interaction (EESX )EES + EEX ) -3.2 kcal/mol) is much smaller than the polarization and charge-transfer interactions; the interaction resulting from occupied virtual orbital mixing gives a dominant contribution at the stable molecular geometry. Nevertheless, the range of EES (long-range, 1/r) is different from that of EEX (shortrange, exponential decay), and thus, their magnitudes can be mutually quite different in some unstable structures. Hence EES and EEX, in addition to EPL and ECT, may play an important role in the dynamics of the WA molecule, e.g., in the photoisomerization of RET. 4.3. Excitation Energy. The excitation energies of the RET low-lying states in bR are calculated by using the ab initio QM/ MM method, and then, an analysis is made to elucidate how the protein environment affects the properties of these states. The electronic wave functions of three low-lying π-π* states (S0, S1, and S2) are calculated by the state-averaged CASSCF method with the DZV basis functions. To improve the flexibility of basis functions, polarization functions49 are added on the C and N atoms in the conjugated polyene chain part of the Schiff base retinal. The CASSCF calculation with total valence π active space, which consists of 12 electrons in 12 orbitals (12,12), is too computationally demanding. Instead, the CASSCF with (12,9) active space is employed in the present calculation, since it is found to provide well-balanced excitation energies and orbitals. It is well-known that the energies of the ionic π-π* states of

g

S0 S1 S2

bR

dipole momentf VES

EEROd,e

0.0 0.0 -104.8 8.0 (5.3) 47.4 63.5 -92.3 11.6 (17.6) (66.1) (90.9, 103.5h) 70.8 82.6 -98.9 13.9 (14.8) (97.8) (113.1, 120.8h)

gas

bR

18.4 1.4

24.5 12.5

1.6

18.3

a Electrostatic interaction energy between the QM and the MM segments. b Electronic reorganization energy. c Results with the model C type QM/MM division at the model C optimized geometry. DZV and additional polarization basis functions are used (see text). d Results by the MRMP calculation. Values in parentheses are the CASSCF results. e VLJ and Vbond are -69.4 and 0.3 kcal/mol, respectively, for all states. f Results by the CASSCF calculation. g The absolute energy (the correlation energy) of the ground state in gas phase is -948.6848 au (-2.6091 au) in MRMP. h Results with the model E type QM/MM division at the model E optimized geometry. DZV(O+) and additional polarization functions are used (see text).

conjugated polyenes are significantly improved by taking into account the σ-π dynamic correlation.62 We therefore carried out further refinements on the CASSCF result for the excitation energies by the MRMP method developed by Hirao.63-66 The excitation energies and the properties of the low-lying states calculated by this CASSCF-MRMP method at the model C geometry are listed in Table 5. The excitation energies of the S1 and the S2 states in the gas phase are estimated to be 47.4 and 70.8 kcal/mol, respectively, by the MRMP calculation; these values are significantly reduced from their CASSCF excitation energies, 66.1 and 97.8 kcal/mol, respectively. The electronic structures of these low-lying states, S1 and S2, are ionic (noncovalent).67 This ionic feature is responsible for the large σ-π dynamic correlation in these states. The dipole moments of the chromophore in the S1 and S2 states are considerably smaller (∼1.5 D) compared with that in the ground state (18.4 D). The electronic structure of the S1 state is dominantly described with the HOMO f LUMO singly excited configuration. The S2 state is also characterized by a singly excited configuration, but of the next HOMO f LUMO.67 Due to the positive charge around the Schiff base, the HOMO and the next HOMO are more localized on the β-ionone ring side of the polyene chain and the localization of electron density gives the large dipole moment of the ground state. As the LUMO is localized on the Schiff base side, the HOMO f LUMO and the next HOMO f LUMO excitations effectively lead the delocalization of the positive charge of the Schiff base along the polyene chain. The dipole moments are therefore reduced in the excited states. In bR, the excitation energies of these low-lying excited states are largely blue-shifted from the gas-phase values; the MRMP excitation energies of the S1 and S2 states are calculated to be 63.5 and 82.6 kcal/mol, respectively (see Table 5), which are 16.1 and 11.8 kcal/mol larger than those in gas phase, respectively. In the ground state, the Schiff base of the chromophore is strongly stabilized by the negatively charged protein residues in the vicinity of the Schiff base, namely, Asp85 and Asp212, and by the WA water molecule (see Figure 3 for their locations near the Schiff base). The low-lying excites states, S1 and S2, have the delocalized positive charge distribution, as just discussed, and so are stabilized much less than the ground state by the surrounding protein-water environment; their total solvation energies are still large, for example, Esolv ) -149.8 kcal/mol for S1 and -154.1 kcal/mol for S2, which are 10% and 7% smaller than the ground-state solvation energy, respec-


Figure 6. Coulombic energy differences between the S1 and S0 states, VES(S1) - VES(S0), in terms of individual protein residues in the model C type division QM/MM system at the model C optimized geometry. The CASSCF level with DZV and additional polarization functions is included (see text).

tively. It is noted that the electrostatic field by those protein residues and water induces the large electronic polarization of these excited states, since their electronic structures are noncovalent in nature; for example, the dipole moment of the S1 state becomes 12.5 D in bR, which is significantly larger than its gas-phase value (1.4 D). The difference of VES between the ground and S1 states is estimated to be 12.5 kcal/mol (Table 5). It is the main contribution that increases the S1 excitation energy (by 16.1 kcal/ mol) from the gas phase to bR. Figure 6 shows how individual protein residues and water molecules contribute to the differences, VES(S1) - VES(S0). We can see in the figure that there are large contributions from Asp85, Asp212, and WA. There exist sizable contributions from the charged residues, Argl34 and Glu194, and from the neutral residues in the vicinity of RET, namely, Trp86, Met145, and Trp182. We have seen in Figure 5 that Arg82 has a large contribution to VES (for example, 28.6 kcal/mol in the S0 state), but its difference between the S1 and S0 states is quite small (0.2 kcal/mol). This can be explained as follows. Upon the excitation, the dipole moment of RET decreases significantly along the polyene chain. Its magnitude becomes only half of the ground state value. Arg82 is located, however, nearly perpendicular to the direction of the RET chain, and thus, the charge (Arg82)-dipole (RET) interaction is almost unaffected by the change of the RET dipole moment. It should be noted from Table 5 that the electronic reorganization induced by the strong electronic polarization also contributes to the S1 state destabilization; EERO in the S1 state is found to be 3.6 kcal/ mol larger than that in the ground state. The larger excitation energy of the S2 state in bR than that in gas phase is explained almost in the same manner. The S2 state having a more polarizable wave function is, however, more stabilized by the protein environment state and has a larger dipole moment (18.3 D) than the S1 state. As a result, the energy separation between the S1 and S2 states becomes smaller in bR as compared to the gas phase. Spectroscopic experiments have been carried out for mutants of bR. The absorption maxima of D85N and D212N mutants of bR were observed to be shifted to the longer wavelengths (615 and 585 nm, respectively) than that of native bR (568 nm).68,69 The mutants D85N and D212N yield 20 and 4 times slower fluorescence decay, respectively, than the 500 fs decay of native bR.70 These experimental results for the mutants can be explained well by the results obtained from the present QM/ MM calculation. As discussed above, the electrostatic interactions among the Schiff base and the negatively charged residues such as Asp85 and Asp212 stabilize the ground state more than

Hayashi and Ohmine the S1 (and S2) state and thus make the S0-S1 excitation energy larger. In the D85N and D212N mutants, these negatively charged residues are replaced by neutral ones, and thus, this excitation is reduced. Furthermore, the energy separation between the S1 and S2 states is expected to become larger in the mutants due to the weaker electrostatic interactions. The increases of the S1 and S2 separation in the mutants must induce then the larger potential energy barrier on the RET photoisomerization path in the S1 state and thus lead to the slower decay rate of the fluorescence from the S1 state.34 As shown in the previous section, there exist large quantum electronic (e.g., charge transfer and polarization) interactions among the molecules in the binding pocket in the bR ground state. To examine the quantum effect on the RET excitation energies, we carried out the CASSCF calculation for the model E system including not only RET but also Asp85 and 3H2O as the QM segments; the CASSCF treatment of the full model A QM/MM system is too demanding at present. The calculation is carried out with DZV(O+) and the additional polarization basis functions (see above). The resultant CASSCF excitation energies in bR are also listed in Table 5. Taking into account the quantum electronic interaction, both S0-S1 and S0-S2 excitation energies (103.5 and 120.8 kcal/mol, respectively) increase from those in the model C system. As seen above, the positive charge of the chromophore localizes around the Schiff base in the ground state, whereas it delocalizes along the polyene chain in the S1 and S2 excited states. The charge transfer and polarization effects from the anionic residues and the water molecules around the Schiff base, i.e., Asp85, Asp212, and WA, thus stabilize the ground state more than these excited states. Note that the energy difference between the S1 and S2 states in the model E case (17.3 kcal/mol) is smaller than that in the model C case (22.2 kcal/mol; CASSCF), preferable for the early stage dynamics of the photoisomerization process based on the three-state-model.5,26,27,34 The excitation energies of the protonated RET in the gas phase obtained by the present MRMP calculation are significantly smaller than the previously reported ones by the CASSCF71 and MRSD-CI67 methods. There is no experimental spectrum of the protonated RET in the gas phase to be compared with these calculated values. The large reduction of the excitation energies is also seen in bR by improving the level of calculation from the CASSCF to the MRMP calculation (Table 5). The bR excitation energies calculated by the present MRMP calculation (63.5 kcal/mol for the S1 state and 82.6 kcal/mol for the S2 state) are, however, still largely overestimated in comparison with the experimental values, 50.0 and 58.6 kcal/ mol,72 respectively. This discrepancy might arise from the fact that we have used basis functions that are too small and l active spaces that are too small to describe the electronic wave functions of the QM segment, and we have also have neglected both the electronic polarization of the protein environment73 and the charge transfer between the QM and the MM segments in the QM/MM calculation evaluating the solvation energy of the chromophore. Further improvements are thus required for more quantitative comparisons. 5. Potential Energy Curve of Proton Transfer The Schiff base is protonated until the L-intermediate state in the photocycle is reached, and then Asp85 accepts the proton from the Schiff base in the M-intermediate. There exists, however, a path of the concerted proton transfer from the Schiff base to Asp85 through the WA water molecule even in the ground-state geometry shown in Figure 3, although the proton transfer does not take place in the real system until a later stage.


Figure 7. Schematic representations of two states in a proton transfer from the Schiff base to Asp85 in the ground-state bR. The protontransfer process changes the charge states of the Schiff base and Asp85 from the zwitterionic (proton non-transferred state) to the neutral character (proton-transferred state). The values r1 and r2 indicate the distance between H2 of WA and Oδ1 of Asp85 and the distance between H16 of the Schiff base and O of WA, respectively (see also Figure 1). The value r1 decreases along the proton transfer.

We analyze here this hypothetical proton pathway and its environmental effect and see why the proton transfer does not occur in this stage. This may give us a clue for finding the forces which drive the proton transfers in the later stages of the photocycle. Along this hypothetical proton transfer, the electronic structure of the Schiff bases moiety is altered drastically, that is, from the zwitterionic structure, consisting of the protonated Schiff base and the unprotonated Asp85, to a neutral structure (Figure 7). The electrostatic field from the protein environment thus must significantly alter and affect the potential energy along this proton-transfer process. The importance of this protein environmental effect was demonstrated by Hadzˇi et al.74,75 by using a theoretical model using an external point charge and dielectric continuum picture, which is, however, not a strict microscopic model. In the present study, we attempt to construct a microscopic molecular model of the proton transfer in bR by using the QM/MM method. The distance (r1) between the WA molecule and the Oδ1 atom of Asp85 is taken as the coordinate describing the proton transfer (Figure 7); r1 decreases along the proton transfer. Using the QM/MM method with the QM segment consisting of RET, 3H2O, Asp85, and Asp212, the total system geometry is optimized with the 3-21G(O+) basis (equivalent to model B level calculation in section 3) at a given set of r1 values, and then the potential energy of the total system is evaluated with a DZV(O+) basis along this “minimum energy path” (MEP) 76 (RET,3H2O,Asp85,Asp212/DZV(O+)/3-21G(O+) denoted as model I). The resultant potential energy curve along MEP of the proton transfer is shown in Figure 8. There exist two potential minima at r1 ) 1.59 and 1.07 Å, corresponding to the zwitterionic state and the neutral state, respectively. Figure 9shows how the distance between the hydrogen atom of the Schiff base and the O atom of the WA water molecule, r2 defined in Figure 7, changes along MEP. At the zwitterionic state (r1 ) 1.59 Å), r2 is 1.50 Å; that is, the Schiff base has the proton whereas Asp85 is unprotonated, and thus, its carboxyl group is negatively charged. The value r2 decreases slowly with the decrease of r1 and then changes suddenly near r1 ) 1.2 Å, where the Schiff base releases a proton and Asp85 accepts another proton through WA in a concerted fashion and then the moiety becomes charge neutral. The system has a minimum of the neutral state at r1 ) 1.07 Å and r2 ) 1.01 Å. We can see in Figure 8 that the zwitterionic state is 4.9 kcal/mol more stable than the neutral state in the bR ground state, which is consistent with the experimental result.


Figure 8. Potential energy curve of the proton-transfer process along the minimum energy path (MEP). The value r1 is the distance between H2 of WA and Oδ1 of Asp85, which decreases along the proton transfer. The systems is of the zwitterionic state at r1 ) 1.59 Å and of the neutral state at r1 ) 1.07 Å. The thick line with + represents the potential energy of the total bR system along the MEP (model I); the thin line with × represents the potential energy of the QM segment only (in the gas phase) calculated for its geometries in the MEP of model I (model II); the thin line with 0 represents the full potential energy of the bR system without WB for the geometries in the MEP of model I (model III). The symbols +, ×, and 0 indicate the potential energies of model I, II, and III, respectively, calculated for the geometries of the MEP at the r1 values. The lines are drawn by interpolating the potential energies.

Figure 9. The distance between H16 of the Schiff base and O of WA (r2) along the proton-transfer process of model I. The value r1 is the distance between H2 of WA and Oδ1 of Asp85, which decreases along the proton transfer.

It is noted that the abrupt change of r2 occurs near r1 ) 1.2 Å when the electronic character of the QM segment changes from the zwitterionic to the neutral state. This electronic character change seemingly suddenly occurs when MEP is plotted along coordinate r1 only. It is smooth when it is plotted two-dimensionally with r1 and r2 or if the intrinsic reaction coordinate (IRC)77 is employed; IRC yields a large curvature around this point but is smooth. The reaction coordinate of this proton transfer thus must be expressed in terms of IRC or the two coordinates, r1 and r2. The searching IRC or the full plot of two-dimensional potential energy surface is, however, computationally impractical for this large system at present. To clarify the importance of the protein environment for this proton transfer, we evaluated the potential energy curve of the QM segment only (i.e., in the gas phase) for its geometries along MEP of model I (model II). The potential energy curve in the gas phase, plotted in Figure 8, is significantly different from


Figure 10. Coulombic interaction energy between the QM and MM segments (VES) along the proton-transfer process of model I.

Figure 11. Coulombic energy differences between the neutral and the zwitterionic states, VES(neutral) - VES(zwitterionic), in terms of contributions from individual residues. The residues giving large absolute contributions (g1.0 kcal/mol) are Met20 (1.0 kcal/mol), Val49 (-1.3 kcal/mol), Met56 (1.5 kcal/mol), Tyr57 (1.0 kcal/mol), Ala81 (-1.8 kcal/mol), Arg82 (2.1 kcal/mol), Trp86 (1.6 kcal/mol), Thr89 (5.3 kcal/mol), and Asp212 (1.1 kcal/mol).

that in bR. As contrasted with the potential energies in bR, the neutral state is 5.3 kcal/mol more stable than the zwitterionic state in the gas phase (model II). The abrupt change of the potential energy curve occurs at r1 ) 1.2 Å; the difference between this energy curve of model II and that of model I arises from the interactions with the protein environment (i.e., the MM segment) in bR. Figure 10shows how the Coulombic interaction energy VES between the QM and the MM segments (see eq 1) change along MEP in model I (in bR). We can see that the zwitterionic state yields much larger electrostatic stabilization (by 11.7 kcal/mol) than the neutral state. This electrostatic stabilization alters the potential energy of the zwitterionic state to be lower than the neutral-state energy in bR. Figure 11 shows how individual protein residues contribute to the Coulombic energy difference between the neutral and the zwitterionic states, VES(neutral) - VES(zwitterionic), in model I (in bR). We can see that the largest contribution to the zwitterionic state stabilization arises from Thr89. Asp85 has the anionic character in the zwitterionic state and thus is largely stabilized by the hydrogen bond of the OH group of Thr89, seen in Figure 7. In the neutral state, this hydrogen bond is significantly weakened and its distance increases by about 0.02 Å. The second dominant contribution arises from Arg82 (2.1 kcal/mol). Although the absolute interaction energy of Arg82 in VES (-92.5 kcal/mol at the zwitterionic state) is much larger than that of Thr89 (-21.0 kcal/mol), the Arg82 interaction is not much different for the zwitterionic and neutral states because

Hayashi and Ohmine it is a long-range Coulomb interaction, whereas the Thr89 interaction is short-range (hydrogen bond). We can also see in Figure 11 that there exist non-negligible contributions from several other residues, that is, the interactions of the backbone CdO and the N-H groups of the residues with the QM segment and those of the Ne-H group in the indole of Trp86 with the QM segment. In addition to Thr89, the WB water molecule is hydrogen bonded to Asp85 (see Figure 3). In the above analysis, WB does not appear explicitly because it is treated as part of the QM segment. To assess the role of WB in the proton transfer, we calculate the potential energies of the bR system without including the WB water molecule, at the geometries of MEP in bR (denoted as model III). The potential energy curve is then significantly changed, as seen in Figure 8. The zwitterionic state is destabilized by lack of the hydrogen bond between Asp85 and WB, and the formation of the neutral state becomes exothermic. This shows that the water molecules in the retinal binding pocket play a crucial role in controlling acidity of Asp85. By the photoisomerization of the retinal, the Schiff base is displaced, causing a significant rearrangement of the HBN of the water molecules in the retinal binding pocket, as suggested by Edman el al.12 and Luecke et al.13 based on their X-ray crystallographic structures of the K- and M-intermediates. This rearrangement of the HBN is expected to be one of major factors which control the proton-transfer process in bR. This point must be further investigated by the detailed potential energy analysis and the trajectory study for the photocycle states after the L-intermediate. 6. Vibrational Structure In the preceding section, we discussed the importance of the HBN rearrangement to control the proton-transfer processes in the photocycle of bR. One of the most powerful experimental techniques to detect this rearrangement of the HBN is Fourier transform infrared (FTIR) spectroscopy. Recently, Kandori et al.16 have observed the difference in IR spectra between the bR and the K-intermediate states in a whole mid-infrared frequency region and measured the angles of the IR intensities. By substituting D2O, they found that O-D or N-D stretching vibrations, which are sensitive to the D2O substitution, are widely distributed in the frequency region 2100-2700 cm-1. This shows that various kinds of hydrogen bonds are involved in these O-D and N-D bonds. Although precise structural information of the water molecules78 and threonines in bR was obtained from IR signals,79 no N-D signal of the Schiff base has been assigned. The angular dependence of the IR signals showed that the angles of almost all these bands are larger than the so-called magic angle (52.8°) (the angles here are defined relative to the normal of the membrane surface); the IR intensities decrease when the surface of the sample (the purple membrane) is tilted. The O-D bands of some threonines have an angle near 0°. In the bR structure obtained in section 3, the direction of the N-D bond of the Schiff base is shown to be about normal to the membrane surface, and thus, if the angles of the bond and the IR intensity mutually match, the angle of the N-D IR intensity is expected to be small. In the experimental spectrum, there remains, however, no distinct band yielding a small angle after assigning a small angle band to be a O-D mode of Thr89.79 The vibrational structure of bR is analyzed by using the QM/ MM method. Here we focus on the N-D band of the Schiff base and the O-D bands of three water molecules in the binding pocket because these bands are expected to change most



TABLE 6: IR Spectroscopic Properties of the Water O-D and the Schiff Base N-D Stretching Vibrationsa

mode

assignment

frequency (cm-1)

V1 V2 V3 V4 V5 V6 VND

O-D1 (WB) O-D1 (WA) antisym (WC) sym (WC) O-D2 (WB) O-D2 (WA) N-D

2721 2618 2601 2468 2320 2257 2113

intensityb

angle of intensityc (deg)

angle of bondc (deg)

1.8 12.2 5.7 3.8 18.2 12.1 46.8

18 81 51 53 48 54 52

45 57 50d 74e 26 65 19

a Results of the model A system. b IR intensities are scaled by a value of the D2O symmetric mode in gas phase obtained by the quantum chemistry calculation with the DZV basis set (the intensity ∝ d2/ω ) 3.80 × 10-2 au, where d is the dipole derivatives and ω is the frequency). c Angles are relative to the normal of the membrane surface. For example, the O-D angle is 90° when the O-D bond is parallel to the membrane surface. d Angle of O-D1(WC). e Angle of O-D2(WC).

significantly upon the photoisomerization of RET. The QM/ MM method with the QM segment consisting of RET, Asp85, Asp212, and 3H2O, in the DZV(O+) basis set, is employed for the vibrational analysis. The geometry of bR is optimized with the QM/MM method with DZV(O+) (i.e., model A geometry in section 3), and the normal-mode analysis is performed using the Hessian matrix obtained by numerical differentiation of the energy gradients. To reduce the computational effort, the nuclei taken into account in constructing the Hessian matrix are limited; we calculate only the Hessian elements of C, N16, H16, and C15 atoms of RET; all H and O atoms of the water molecules; and 0δ and Cγ atoms of Asp85 and Asp212, while the coordinates of the other atoms are kept fixed. We have checked the effect of using this limited Hessian matrix by employing the larger matrix sizes in smaller basis sets 3-21G(O+) and found that the effect on the normal modes of the O-D and N-D stretching vibrations (their frequencies, components, and intensities) is small. It is well-known that the vibrational frequencies tend to be ∼10% overestimated by HF calculations as compared with experimental values. We therefore correct all the HF vibrational frequencies by scaling with a factor, 0.9008, which is determined to reproduce the experimental value of the antisymmetric mode of D2O in the gas phase (2788 cm-1). The IR intensities of the modes are calculated by taking the dipole derivatives with respect to the normal coordinates, and their angles are determined (see also table caption). The frequencies, the intensities, and the assignments of the normal modes of the water molecules and N-D of RET in bR are listed in Table 6. The IR intensities are scaled by a value of the D2O symmetric mode in gas phase obtained by the quantum chemistry calculation (DZV basis set). We can see that the six O-D stretching modes of the three water molecules, V1-V6, yield a wide frequency distribution, 2257-2721 cm-1. The N-D stretching vibrational mode of the Schiff base has a low frequency, 2113 cm-1. It is found from analyzing the mode components that there exist small but non-negligible intra- and intermolecular couplings among the O-D and N-D stretching motions. We can see in Table 6 that the angles of IR intensities are significantly different from the angles of the O-D bonds themselves. This is because the IR intensity is very sensitive to intra- and intermolecular vibrational mixing and also to the electron-vibration interaction as discussed below. The V1, V2, V5, and V6 normal modes are mainly described as the decoupled O-D vibrations of the WA and WB water molecules. The two modes, V2 (2618 cm-1) and V6 (2257 cm-1), are the O-D vibrations of the WA water molecule. The

frequency of an O-D stretching is a measure of its hydrogen bond strength. The higher frequency V2 mode is mainly localized on the O-D bond, which is making only a weak hydrogen bond with Asp212. On the other hand, the O-D of the lower frequency V6 mode is strongly hydrogen bonded to Asp85. The IR intensities of those (V2, V6) modes are found to be about 12 times larger than that of the symmetric O-D vibration of a water molecule in the gas phase. The WA molecule is bound to the positively and negatively charged sites, i.e., to Asp85, Asp212, and the Schiff base, and, thus, yields the strong electrostatic and charge-transfer interactions from/to those residues. These intermolecular interactions are responsible for the extensive IR intensities amplification of these WA O-D modes. The V1 (2721 cm-1) and V5 (2320 cm-1) modes are assigned to the O-D vibrations of the WB molecule. The frequency and IR intensity of the V1 mode are comparable to those in gas phase because the O-D bond corresponding to the V1 mode does not make a hydrogen bond with any residues (Figure 3). The other O-D bond of the WB molecule is, on the other hand, strongly hydrogen bonded to Asp85, so the V5 mode has a low frequency and a large IR intensity (18.2). The V3 (2601 cm-1) and V4 (2468 cm-1) modes are assigned to be the antisymmetric and the symmetric vibrations of the WC molecule, respectively. The O-D vibrations of the WC molecule are strongly coupled to each other. These frequencies are reduced by the hydrogen bonds with WB and Asp212. It is possible that the V3 and V4 modes are poorly described in the present calculation because Arg82 attached to the O atom of WC is only treated as a part of the MM segment. The IR intensities of those modes are calculated to be rather small (5.7 and 3.8) in the present calculation but might increase considerably when the electronic (charge transfer and polarization, etc.) interactions of the Arg82 residue is properly accounted for. Kandori et al.78 experimentally assigned the O-D modes of D2O molecules in the K minus bR spectra and also found that all the angles of the IR intensities are larger than the magic angle (52.8°). The present calculation also provides the angle values larger than 48°, in very good agreement with the experimental results, except for the V1 mode (18°). The O-D bond corresponding to the V1 mode is not hydrogen bonded with any residues, and thus, the angle of the O-D bond might thermally fluctuate to a great extent and, therefore, cannot be well described with a static model, like a normal mode calculation with an optimized structure, as done in the present calculation. The N-D stretching vibrational mode of the Schiff base is found to be almost localized on the N-D bond. We can see in Table 6 that the frequency of this mode is very low (2113 cm-1) because this N-D bond is strongly hydrogen bonded to the WA molecule. Moreover, the IR intensity of this band is especially large (46.8) and the angle of IR intensity (520) is totally different from the angle of the bond (190). This anomalous behavior of the IR intensity arises from the electron-vibration interaction.80 It is found that the dipole derivatives with respect to the local N-D stretch mode is directed along the total dipole moment of the QM system, not along the N-D bond itself. As seen in the previous section, the N-D stretching motion is strongly coupled with the motion of the proton transfer from the Schiff base to Asp85 (Figure 9), which involves the distinct change of the electronic structure from the zwitterionic to the neutral state (Figure 7). The N-D stretching motion thus induces the considerable electronic structure reorganization of the Schiff basis moiety, yielding the strong IR intensity and the significant difference in the angle between the N-D bond and the IR

10688 J. Phys. Chem. B, Vol. 104, No. 45, 2000 intensity. It is noted that the flexibility of the basis function used in present calculation might not be enough to accurately describe the large change of the electronic structure involved in the N-D stretching motion, and thus, its resultant frequency (2113 cm-1) can be underestimated. Nevertheless, it is possible that there exists the hidden N-D mode in a low-frequency region (j r

+

iA

+

ij

∑ A>B

Z AZ B

+

rAB

Zaqp MM κap + V QM-MM + V MM MM (A-2) rap

∑a ∑p

where the indices i and j (A and B) represent the electrons (the nuclei) in the QM segment that includes the dummy atoms and a (p) represents the “interaction sites” in the QM segment (variables for the MM segment in parentheses). The first four terms of the right-hand side constitute the gas-phase electronic Hamiltonian. The Coulombic interaction energies between the electrons (nuclei) in the QM segment and the effective charges in the MM segment are represented by the fifth and sixth terms, respectively. The effective charges of the electrons in the dummy atoms are set to be qdum ) -Zdum so that the electrostatic interactions between a dummy atom and the MM segment mutually cancel out; the terms expressing those interactions do not appear in eq A-2. In the fifth term, qˆ a is the “effective charge operator” of the electrons in the QM segment. The prefactor of the fifth term, γ, is a scaling factor by which the total charge of the QM/MM system is preserved; this factor is determined by the condition

γN ′e + N′Z +

∑

I)

(

qa

∑a ∑a r ωa

)

2

- V′R(P) + 2λe

aR

( ∑ q - N ′) + a

e

a

∑a ga(qa + Za)2

(A-5)

where R is index of the grid points surrounding the QM segment, V ′R(P) is the electrostatic potential of the electronic wave function at the grid point, and ωR is the weight factor. In the present study, the grid points are set to coincide with the interaction sites of the MM segment, and then the square of the corresponding effective charges, q2p, are used as the weight factors, ωp. Here, ωp ) 0 for those p having the 1-2 and 1-3 nonbonding interactions with the QM segment interaction sites. This enables us to estimate the electronic interaction between the QM and the MM segments more reasonably. The electrostatic contribution from the dummy atoms is removed in the electrostatic potential, V R′ (P):

V′R(P) ) VR(P) -

qdum

∑ dum r

(A-6)

dum,R

where VR(P) is the electrostatic potential derived from the electronic density. The total effective charge is constrained to be N e′ with the Lagrange multiplier, λe, by the second term of the right-hand side in eq A-4. N e′ is given by

∑ qdum

N′e ) Ne +

(A-7)

dum

qp ) qtotal

(A-3)

p∈MM

where qtotal is the total charge of the total system and -N ′e (N Z′ ) is the total number of the electrons (the total nucleus charge) in the QM segment without including the dummy atoms. κap in the fifth and sixth terms are the factors for excluding the 1-2 and the 1-3 nonbonded interactions and for scaling the MM 1-4 nonbonding interactions. V QM-MM is the interaction energy between the QM and the MM segments other than the electrostatic interactions. This term is expressed by the MM MM potential functions. The last term, V MM , is the potential interactions within the MM segment. The expectation value of the energy for a HF electronic wave function, E, is expressed as

E ) 〈Ψ|H ˆ |Ψ〉 ) E ele 0 (P) + γ

where Ψ is the HF electronic wave function, E ele 0 (P) is the energy of the MO part whose expression is the same as that in gas phase, and P is the density matrix of electronic state. E MM consists of the terms which do not include the variables of the electrons. The electrostatic interaction between the electrons in the QM segment and the point charges in the MM segment is given by the second term of the right-hand side in eq A-4. In this term, qa(P) is “an electronic effective charge of a interaction site” in the QM segment that depends on the electronic wave function. To define the “effective charge”, we employ the RESP method, where qa is determined by minimizing the following value I:42

qp

∑a ∑p κapr

ap

qa(P) + E MM

(A-4)

where Ne is the number of the electrons in the QM segment. The last term of the right-hand side is introduced to restrain the electronic effective charges. To simplify the SCF procedure, here we used the harmonic penalty function for the purpose. Hence, the RESP effective charges are determined by

q(P) ) a-1 W(P) - a-11

{a}ab )

{W}a )

1

∑R r

ωR

aR

1

∑R r

[

1ta

-1

1 rbR

]

W(P) - N′e(P) -1 1ta 1

(A-8)

+ gaδab

(A-9)

ωRV′R(P) - gaZa

(A-10)

aR

where a and b are the “interaction sites” in the QM segment. The electronic energy is determined by solving the Fock equation derived from eq A-4 with the RESP effective charge


Hayashi and Ohmine

operator. The QM/MM Fock operator consists of two terms:

ˆf ) ˆf 0 + ˆf solv

(A-11)

where the expression ˆf 0 is the same as that in gas phase where and ˆf solv, called the solvated Fock operator,42 gives the electrostatic interaction from the MM segment. The solvated Fock matrix elements in terms of atomic orbitals (AO) are expressed as

∑R ∑a V˜ aΛaRAµν(rR) +

f solv ) - γ

γ

V ˜ t a-11 1ta

-1

1

[∑ ∑ R

-V ˜a )

ΛaR )

]

ΛaRAµν(rR) - Sµν (A-12)

a

qp

∑p r

(A-13) ap

1

∑b (a-1)ab r

ωR

(A-14)

bR

where V ˜ a is the electrostatic potential from the MM segment acting on the site a. The matrix elements have one-electron AO integrals, Sµν and Aµν(rR). The former is the overlap integral, and the latter is the Coulombic interval,

Aµν(rR) )

∫ dr φµ |r -1 rR| φν

(A-15)

Since the electrostatic interaction between the QM and the MM segments is described by the one electron operator, this QM/ MM Hamiltonian can be extended straightforwardly to post HF methods such as MCSCF or MP perturbation theories. The molecular geometry of the QM/MM system is optimized with the analytical energy gradient of the QM/MM system based on the Hamiltonian eq A-2.83 In general, the MM segment is constituted by a large number of atoms compared to the QM segment so it takes many searching steps to optimize the MM segment. The optimization procedure is therefore divided into two steps; the QM optimization with the fixed MM geometry and the MM optimization with the fixed QM geometry and the fixed electronic density, i.e., the fixed effective charges. These steps are carried out iteratively until the total convergence is achieved. Note that the latter step does not include the computational-time-consuming evaluation of derivatives of AO integrals. Hence this iterative procedure enables us to obtain the optimized geometry of QM/MM system very efficiently.76 It is noteworthy that our method ensures the total convergence of the geometry optimization because both the QM and the MM optimization processes obey exactly the same Hamiltonian eq A-2. References and Notes (1) Mathies, R. A.; Lin, S. W.; Ames, J. B.; Pollard, W. T. Annu. ReV. Biophys. Biophys. Chem. 1991, 20, 491. (2) Rothschild, K. J. J. Bioenerg. Biomembr. 1992, 24, 147. (3) Oesterhelt, D.; Tittor, T.; Bamberg, E. J. Bioenerg. Biomembr. 1992, 24, 181. (4) Lanyi, J. K. Biochim. Biophys. Acta 1993, 1183, 241. (5) Schulten, K.; Humphrey, W.; Logunov, M.; Sheves, M.; Xu, D. Isr. J. Chem. 1995, 35, 447. (6) Spudich, J. K.; Lanyi, J. K. Curr. Opin. Cell Biol. 1996, 8, 452. (7) Haupts, U.; Tittor, J.; Bamberg, E.; Oesterhelt, D. Biochemistry 1997, 36, 2. (8) Lanyi, J. K. J. Biol. Chem. 1997, 272, 31209.

(9) Subramaniam, S.; Lindahl, M.; Bullough, P.; Faruqi, A. R.; Tittor, J.; Oesterhelt, D.; Brown, L.; Lanyi, J.; Henderson, R. J. Mol. Biol. 1999, 287, 145. (10) Luecke, H.; Schobert, B.; Richter, H. T.; Cartailler, J. P.; Lanyi, J. K. J. Mol. Biol. 1999, 291, 899. (11) Belrhali, H.; Nollert, P.; Royant, A.; Menzel, C.; Rosenbusch, J. P.; Landau, E. M.; Pebay-Peyroula, E. Structure 1999, 7, 909. (12) Edman, K.; Nollert, P.; Royant, A.; Belrhall, H.; Pebay-Peyroula, E.; Hajdu, J.; Neutze, R.; Landau, E. M. Nature 1999, 401, 822. (13) Luecke, H.; Schobert, B.; Richter, H. T.; Cartailler, J. P.; Lanyi, J. K. Science 1999, 286, 255. (14) Maeda, A.; Sasaki, J.; Shichida, Y.; Yoshizawa, T. Biochemistry 1992, 31, 462. (15) Hatanaka, M.; Kandori, H.; Maeda, A. Biophys. J. 1997, 73, 1001. (16) Kandori, H.; Kinoshida, N.; Shichida, Y.; Maeda, A. J. Phys. Chem. B 1998, 102, 7899. (17) Xu, D.; Sheves, M.; Schulten, K. Biophys. J. 1995, 69, 2745. (18) Humphrey, W.; Xu, D.; Sheves, M.; Schulten, K. J. Phys. Chem. 1995, 99, 14549. (19) Nina, M.; Roux, B.; Smith, J. C. Biophys. J. 1995, 68, 25. (20) Roux, B.; Nina, M.; Pomes, R.; Smith, J. C. Biophys. J. 1996, 71, 670. (21) Petrich, J. W.; Breton, J.; Martin, J. L.; Antonetti, A. Chem. Phys. Lett. 1987, 137, 369. (22) Mathies, R. A.; Brito Cruz, C. H.; Pollard, W. T.; Shank, C. V. Science 1988, 240, 777. (23) Dobler, J.; Zinth, W.; Kaiser, W.; Oesterhelt, D. Chem. Phys. Lett. 1988, 144, 215. (24) Govindjee, R.; Balashov, S. P.; Ebrey, T. G. Biophys. J. 1990, 58, 597. (25) Tittor, J.; Oesterhelt, D. FEBS Lett. 1990, 263, 269. (26) Hasson, K. C.; Gai, F.; Anfinrud, P. A. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 15124. (27) Gai, F.; Hasson, K. C.; McDonald, J. C.; Anfinrud, P. A. Science 1998, 279, 1886. (28) Alfano, R. R.; Govindjee, R.; Becher, B.; Ebrey, T. G. Biophys. J. 1976, 16, 541. (29) Shapiro, S. L.; Campillo, A. J.; Lewis, A.; Perreault, G. J.; Spoonhower, J. P.; Clayton, R. K.; Stoeckenuis, W. Biophys. J. 1978, 23, 383. (30) Logunov, S. L.; Song, L.; El-Sayed, M. A. J. Phys. Chem. 1996, 100, 18586. (31) Birge, R. R.; Cooper, T. M.; Lawrence, A. F.; Masthay, M. B.; Vasilakis, C.; Zhang, C. F.; Zidovetzki, R. J. Am. Chem. Soc. 1989, 111, 4063. (32) Yamato, T.; Kakitani, T. Photochem. Photobiol. 1997, 66, 735. (33) Xu, D.; Martin, C.; Schulten, K. Biophys. J. 1996, 70, 453. (34) Humphrey, W.; Lu, H.; Logunov, I.; Werner, H. J.; Schulten, K. Biophys. J. 1998, 75, 1689. (35) Ben-Num, M.; Molnar, F.; Lu, H.; Phillips, J. C.; Martinez, T. D.; Schulten, K. Faraday Discuss. 1998, 110, 447. (36) Molnar, F.; Ben-Num, M.; Martinez, T. J.; Schulten, K. J. Mol. Struct.: THEOCHEM, a special WATOC issue in press. (37) Garavelli, M.; Celani, P.; Bernardi, F.; Robb, M. A.; Olivucci, M. J. Am. Chem. Soc. 1997, 119, 6891. (38) Vreven, T.; Bernardi, F.; Garavelli, M.; Olivucci, M.; Robb, M. A.; Schlegel, H. B. J. Am. Chem. Soc. 1997, 119, 12687. (39) Garavelli, M.; Vreven, T.; Celani, P.; Bernardi, F.; Robb, M. A.; Olivucci, M. J. Am. Chem. Soc. 1998, 120, 1285. (40) (a) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 6, 700. (b) Gao, J. Acc. Chem. Res. 1996, 29, 298. (c) Monard, G.; Merz, K. M., Jr. Acc. Chem. Res. 1999, 32, 904 and references therein. (d) Maseras, F.; Morokuma K. J. Comput. Chem. 1995, 16, 1170. (e) Merill, G. M.; Gordon, M. S. J. Phys. Chem. A 1998, 102, 2650. (41) Logunov, I.; Schulten, K. J. Am. Chem. Soc. 1996, 118, 9727. (42) Ten-no, S.; Hirata, F.; Kato, S. J. Chem. Phys. 1994, 100, 7443. (43) Bayly, C.; Cieplak, P.; Cornell, W.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269. (44) Luecke, H.; Richter, H. T.; Lanyi, J. K. Science 1998, 280, 1934. (45) Grigorieff, N.; Ceska, T. A.; Downing, K. H.; Baldwin, J. M.; Henderson, R. J. Mol. Biol. 1996, 259, 393. (46) Sasaki, J.; Lanyi, J. K.; Needleman, R.; Yoshizawa, T.; Maeda, A. Biochemistry 1994, 33, 3178. (47) Brown, L. S.; Sasaki, J.; Kandori, H.; Maeda, A.; Needleman, R.; Lanyi, J. K. J. Biol. Chem. 1995, 270, 27122. (48) Binkley, J. S.; Pople, J. A.; Hehre, W. J. J. Am. Chem. Soc. 1980, 102, 939. (49) Dunning, T. H., Jr.; Hay, P. J. In Method of Electronic Structure Theory; Shaefer, H. F., III, Ed.; Plenum Press: New York, 1977; Chapter 1. (50) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; von Schleyer, P. R. J. Comput. Chem. 1983, 4, 294.

Proton Transfer in Bacteriorhodopsin (51) 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. (52) Jorgensen, W. L.; Chandreskhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1982, 79, 926. (53) Dupuis, M.; Chin, S.; Marquez, A. In RelatiVistic and Electron Correlation Effects in Molecules and Clusters; Malli, G. L., Ed.; NATO ASI Series; Plenum Press: New York, 1992. (54) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. J.; Windus, T. L.; Dupuis, M.; Montgomery, J. S. J. Comput. Chem. 1993, 14, 1347. (55) Tajkhorshid, E.; Baudry, J.; Schulten, K.; Suhai, S. Biophys. J. 2000, 78, 683. (56) Hatanaka, M.; Kashima, R.; Kandori, H.; Friedman, N.; Sheves, M.; Needleman, R.; Lanyi, J. K.; Maeda, A. Biochemistry 1997, 36, 5493. (57) Sato, H.; Hirata, F. J. Phys. Chem. 1998, 102, 2603. (58) Tajkhorshid, E.; Suhai, S. Theor. Chem. Acc. 1999, 101, 180. (59) Morokuma, K. J. Chem. Phys. 1971, 55, 1236. (60) Kitaura, K.; Morokuma, K. Int. J. Quantum Chem. 1976, 10, 325. (61) Chen, W.; Gordon, M. S. J. Phys. Chem. 1996, 100, 14316. (62) Nakayama, K.; Nakano, H.; Hirao, K. Int. J. Quantum Chem. 1998, 66, 157. (63) Hirao, K. Chem. Phys. Lett. 1992, 190, 374. (64) Hirao, K. Chem. Phys. Lett. 1992, 196, 397. (65) Hirao, K. Chem. Phys. Lett. 1993, 201, 59.

J. Phys. Chem. B, Vol. 104, No. 45, 2000 10691 (66) Nakano, H. J. Chem. Phys. 1993, 99, 7983. Nakano, H. MR2D, version 2; University of Tokyo: Tokoyo, 1995. (67) Du, P.; Davidson, E. R. J. Phys. Chem. 1990, 94, 7013. (68) Mogi, T.; Stern, L.; Marti, T.; Chao, B.; Khorana, H. Proc. Natl. Acad. Sci. U.S.A. 1988, 85, 4148. (69) Needleman, R.; Chang, M.; Ni, B.; Varo, G.; Fornes, J.; White, S. H.; Lanyi, J. K. J. Biol. Chem. 1991, 266, 11478. (70) Song, L.; El-Sayed, M. A.; Lanyi, J. K. Science 1993, 261, 891. (71) Yamamoto, S.; Wasada, H.; Kakitani, T. J. Mol. Struct.: THEOCHEM 1998, 451, 151. (72) Birge, R. R.; Zhang, C. F. J. Chem. Phys. 1990, 92, 7178. (73) Houjou, H.; Inoue, Y.; Sakurai, M. J. Am. Chem. Soc. 1998, 120, 4459. (74) Hodosˇcˇek, M.; Hadzˇi, D. Can. J. Chem. 1985, 63, 1528. (75) Hadzˇi, D. J. Mol. Struct. 1988, 177, 1. (76) Zhang, Y.; Liu, H.; Yang, W. J. Chem. Phys. 2000, 112, 3483. (77) Fukui, K.; Kato, S.; Fujimoto, H. J. Am. Chem. Soc. 1997, 97, 1. (78) Kandori, H. Personal communication. (79) Kandori, H.; Kinoshita, N.; Yamazaki, Y.; Maeda, A.; Shichida, Y.; Needleman, R.; Lanyi, J. K.; Bizounok, M.; Herzfeld, J.; Raap, J.; Lugtenburg, J. Biochemistry 1999, 38, 9676. (80) Torii, H.; Tasumi, M. J. Phys. Chem. B 1997, 101, 466. (81) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1986, 7, 718. (82) Field, M. J.; Bash, P. A.; Karplus, M. A. J. Comput. Chem. 1990, 11, 700. (83) Sato, H.; Hirata, F.; Kato, S. J. Chem. Phys. 1996, 105, 1546.

Proton Transfer in Bacteriorhodopsin - American Chemical Society

Recommend Documents