INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX REFERENCES

Bacteriorhodopsin (bR), a transmembrane protein residing in the purple membrane of Halobacterium salinarum, functions as a light-driven proton pump to produce a proton gradient across the membrane (Oesterhelt and Stoeckenius, 1971). The protein consists of seven transmembrane -helices (Henderson and Unwin, 1975) surrounding an all-trans retinal chromophore covalently bound to Lys-216 through a protonated Schiff base linkage. Fig. 1 depicts the chemical structure of the protonated all-trans retinal Schiff base and its link to the protein. Retinal contains six conjugated double bonds including the Schiff base group.

View larger version (7K): [in this window] [in a new window]   FIGURE 1   Schematic representation of all-trans retinal bound via a protonated Schiff base to Lys-216 in bR and its conventional atom numbering.

Upon absorption of light in the visible spectral region, the retinal chromophore undergoes in bR an isomerization around the C₁₃==C₁₄ double bond, resulting in a 13-cis configuration. The photoisomerization initiates a photocycle during which bR completes an active proton transport event across the membrane from cytoplasm to the extracellular side. The photocycle comprises a series of intermediates, labeled K, L, M, N, and O, characterized by their absorption maxima, and is completed through return to the initial (BR) state (Lozier et al., 1975). The lifetimes of the intermediates range from subpicoseconds to milliseconds. The photocycle involves protein structural changes, which are coupled to proton transfer between such key residues as the retinal Schiff base, Asp-85, Asp-96, and Glu-204 (Braiman et al., 1988; Gerwert et al., 1989; Brown et al., 1995). A sequence of proton transfers between these residues completes the overall unidirectional proton pump through the channel.

The three-dimensional atomic structure of the resting BR state has been determined by cryoelectronic microscopy (Grigorieff et al., 1996; Mitsuoka et al., 1999) and x-ray crystallography (Pebay-Peyroula et al., 1997; Luecke et al., 1998, 1999a,b; Belrhali et al., 1999). These structures reveal the position of key residues facing the interior of the protein and forming the proton channel. The protonated Schiff base group lies in the middle of the channel dividing it into the cytoplasmic and extracellular half-channels. High resolution x-ray structures (Luecke et al., 1999a; Belrhali et al., 1999) have also shown a hydrogen-bond network (HBN) in the extracellular half of the channel including several internal water molecules. The existence of water molecules was first suggested by nuclear magnetic resonance (NMR) (de Groot et al., 1989) and Fourier transform infrared (FTIR) spectroscopic studies (Maeda et al., 1994; Fischer et al., 1994; Kandori et al., 1995; Chon et al., 1996). In Fig. 2 a, the structure of the retinal binding site of the x-ray structure of BR (Luecke et al., 1999a; PDB code, 1C3W) is shown. Fig. 3 illustrates the topology of the HBN in the retinal binding pocket, including the negatively charged Asp-85 and Asp-212, the Schiff base, and three internal water molecules (Wat-401, Wat-402, and Wat-406).

View larger version (23K): [in this window] [in a new window]   FIGURE 2   Comparison of the x-ray crystallographic structures in the vicinity of the Schiff base in (a) BR (Luecke et al., 1999a; PDB, 1C3W) and (b) K (Edman et al., 1999; PDB, 1QKO) states. The -ionone halves of the chromophore are not shown. In the structure of the K intermediate, dislocation of Wat-402 and rotation of the carboxyl group of Asp-85 are observed. Wat-400 in the K model of Edman et al. (1999) corresponds to Wat-406 in the model of BR of Luecke et al. (1999a).

View larger version (21K): [in this window] [in a new window]   FIGURE 3   Schematic representation of the hydrogen-bond network in the binding pocket in the BR state and atom numbering used in the text. Dashed lines indicates possible hydrogen bonds. Distances of the hydrogen bonds are shown in Å.

After photoisomerization, a bathochromic (red-shifted) intermediate, called K, forms in 3 ps. The K state that is trapped at liquid nitrogen temperature 77 K is called K_LT. In K, the energy of the absorbed photon is stored in a high energy structure, the relaxation of which induces a controlled series of proton translocations. The stored energy has been estimated to be 16 kcal/mol by photocalorimetric experiments (Birge and Cooper, 1983; Birge et al., 1989). The structural changes during the K formation have been extensively investigated by various methods. Low temperature resonance Raman (LT-RR) and low temperature FTIR (LT-FTIR) spectroscopic studies have shown significant chromophore conformational deformations in K_LT (Pande et al., 1981; Braiman and Mathies, 1982; Bagley et al., 1982; Rothschild and Marrero, 1982). LT-FTIR studies have also revealed alterations of the HBN in the vicinity of the Schiff base during the K formation (Kandori et al., 1998a,b, 1999, 2000, 2001a; Kandori and Shichida, 2000).

Recently, an x-ray crystallographic structural model of the K_LT intermediate has been reported (Edman et al., 1999). Fig. 2 b shows the retinal binding site of the x-ray structure (PDB code, 1QKO). The KBR difference electron density map reported in that study (Edman et al., 1999) clearly shows a strong negative density at the position of Wat-402 in BR, indicating the dislocation of this water molecule during the formation of K_LT. The density map also reveals a remarkable conformational change of Asp-85, which leads to the complete dissociation of its hydrogen bond with Thr-89 and its closer distance to the Schiff base. The authors have argued that these structural changes play a key role in the primary proton transfer between the Schiff base and Asp-85 in a later (L-to-M) step. However, due to the weak electron density for K_LT, the position of the dislocated W-402 and the conformation of the chromophore in K_LT have not been determined clearly. Moreover, the complete dissociation of the hydrogen bond between Asp-85 and Thr-89 disagrees with recent LT-FTIR results that indicate a stronger hydrogen bond between these groups upon the formation of K_LT (Kandori et al., 1998b, 1999, 2001).

Two distinct K-like bathochromic intermediates have been resolved by visible (Shichida et al., 1983), FTIR (Rothschild et al., 1985; Sasaki et al., 1993; Weidlich and Siebert, 1993; Hage et al., 1996; Dioumaev and Braiman, 1997), and resonance Raman (Althaus et al., 1995; Doig et al., 1991) spectroscopic measurements. These intermediates appear sequentially in time, and the early and late ones have been conventionally called K and KL, respectively. The formation of KL at room temperature has been observed in ~10 ns by the visible absorption measurements (Shichida et al., 1983), and in 50 to 100 ns by the time resolved FTIR (TR-FTIR) experiments (Hage et al., 1996; Dioumaev and Braiman, 1997). Although the spectroscopic studies have suggested significant structural changes involved in the K-to-KL transition, the KL intermediate has not been well understood at the atomic level.

The structure and dynamics of bR have been investigated by many theoretical studies (Schulten et al., 1995; Logunov and Schulten, 1996; Xu et al., 1996; Humphrey et al., 1998; Ben-nun et al., 1998; Tajkhorshid et al., 2000; Baudry et al., 2001; Warshel et al., 1991; Warshel and Chu, 2001; Nina et al., 1995; Yamato and Kakitani, 1997). Recently, ab initio quantum mechanical/molecular mechanical (QM/MM) calculations have been applied to investigate the structure, the proton transfer reaction, and the electronic and vibrational absorptions of the BR state (Hayashi and Ohmine, 2000).

QM/MM methods have been widely used to study chemical reactions in solution phase and in proteins (Warshel and Levitt, 1976; Warshel et al., 1991; Warshel and Chu, 2001; Field et al., 1990; Gogonea et al., 2001; Maseras and Morokuma, 1995; Merill and Gordon, 1998; Gao and Xia, 1992; Logunov and Schulten, 1996). In QM/MM methods, the part of the system that involves electronic changes is described quantum mechanically (QM), and the rest, which provides steric and electrostatic environmental effects on the QM part, is treated by a molecular mechanical (MM) force field.

Hayashi and Ohmine (2000) computed the proton transfer potential energy profile in bR by QM/MM calculations and showed that the primary proton transfer can be efficiently triggered by a certain rearrangement of the HBN, which can be induced by the chromophore photoisomerization (Hayashi and Ohmine, 2000). The authors have also described quantum mechanically the electronic structures of the chromophore and the HBN in the vicinity of the Schiff base, which is closely involved in the proton transfer process. The electronic nature of the HBN results in anomalous IR signals of the OD and ND stretching modes that have been confirmed by a recent IR measurement (Kandori et al., 2002). The excitation energy of the chromophore in the protein matrix was also calculated and analyzed in terms of the chromophore-protein interactions. The ab initio QM/MM approach has also been applied to identify the molecular mechanism of prominent absorption spectral features of sensory rhodopsin II (Hayashi et al., 2001), for which x-ray crystallographic structures have been determined recently (Royant et al., 2001; Luecke et al., 2001).

In the present paper, we report molecular dynamics (MD) simulations and ab initio QM/MM calculations modeling the early intermediates of bR's photocycle. The photoisomerization and the formation of the photoproducts, i.e., the early intermediates, were simulated by MD. As mentioned above, a proper description of electronic effects of the chromophore and the HBN is crucial for the study of the chromophore-protein interactions and their dynamics. Thus, we introduce a polarizable water model (Rick et al., 1994) for internal water molecules and reparameterized the force field functions by the QM/MM method. The structures of the photoproducts obtained by MD simulations are further refined by QM/MM full geometry optimizations, where the chromophore and the HBN are calculated fully quantum mechanically. At the optimized structures of the intermediates, electronic excitation energies of the chromophore were also calculated to characterize the identified intermediates and to elucidate the molecular basis of their bathochromic shifts. The present study provides detailed atomic structural models of the K and KL intermediates and suggests key structural changes in the proton pump process. The study also addresses discrepancies between the x-ray structural model and spectroscopic observations.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX REFERENCES

In this section, details of preparation of the structural model will be first presented. Then, we will describe protocols and potential energy functions used for MD simulation of the chromophore isomerization and the formation of intermediates. Finally, methods of the QM/MM geometry optimization, energy decomposition analysis, and excitation energy calculations will be described.

Structural model

A structural model of the ground state (BR) was constructed based on the high resolution crystal structures (Luecke et al., 1999a; Belrhali et al., 1999). The model includes nine water molecules: two in the cytoplasmic half of the channel, three in the Schiff base binding site (Wat-401, Wat-402, and Wat-406), one between Arg-82 and Thr-205, and three in the vicinity of Glu-194 and Glu-204. Asp-96, Asp-115, and Glu-204 were protonated, and the other acidic and basic residues were assumed to be charged (Sasaki et al., 1994; Brown et al., 1995).

MD simulation

MD simulations were used to model the chromophore photoisomerization and to obtain starting structural models of the early intermediates, which were refined by the QM/MM optimization described in the next section. The AMBER force field (Cornell et al., 1995) was used for the standard protein residues, whereas force fields of the retinal chromophore and water molecules were reparameterized based on the QM/MM methods, as described below. Time integration of the equation of motion was performed by the leap-frog algorithm with a 1-fs time step. The bond distances including hydrogen atoms were kept fixed by SHAKE. To maintain the overall shape of the protein, harmonic constraints with a force constant of 2.0 kcal/mol/Ų were applied to C atoms further than 12 Å from the retinal chromophore. This approximation is validated by electron microscopic (Bullough and Henderson, 1999) and x-ray crystallographic (Edman et al., 1999) studies, showing that no major conformational changes of helices accompany the formation of K. A 12-Å residue-based cutoff was used for nonbonding interactions. The Berendsen method (Berendsen et al., 1984) was used to control the temperature.

The BR ground state was equilibrated for 200 ps at 300 K starting from the QM/MM optimized structure (Hayashi and Ohmine, 2000), and then an equilibrium simulation of BR for 300 ps was performed to sample initial coordinates and velocities for the simulations of photoisomerization. The chromophore photoisomerization was simulated from the initial configurations using the potential energy function of the excited state, which steers the 13-trans excited state to the 13-cis ground state photoproduct (see below). During the isomerization, the temperature control was turned off. To determine energy minimized structures of the photoproducts from trajectories of the photoisomerization, the system was first cooled gradually to 50 K in 40 ps and then subjected to the QM/MM optimization.

For the BR ground state, we used the force field parameters of dihedral angles of the main polyene chain of retinal (set B in Tajkhorshid et al., 2000) determined previously based on density functional calculations (Tajkhorshid et al., 1999). In addition, improper torsional functions with V_n/2 = 0.3 kcal/mol were applied to represent sp² planarities in the conjugated polyene.

For the isomerization, the potential energy function of the C₁₃==C₁₄ dihedral angle, ₁₃₌₌₁₄, was expressed by seven Fourier components listed in Table 1. Fig. 4 shows the energy profile of the potential function. The internal conversion to the 13-cis ground state is described by downhill dynamics on a single potential curve. In this description, the electronic details of the transition process from the excited state to the ground state during the isomerization are neglected. We assumed that the reactive trajectories to the 13-cis configurations are not drastically perturbed by the electronic transition process, because the internal conversion likely takes place through conical intersections, which are expected to be highly efficient doorways to the ground state. The potential energy curve displays a plateau region around the 13-trans configuration, as suggested by quantum chemical calculations (Garavelli et al., 1997, 1998; Humphrey et al., 1998; González-Luque et al., 2000). Unfortunately, the accurate excited state potential profile of the chromophore with six double bonds in the protein matrix has not been calculated yet. We therefore repeated the photoisomerization with several trial potential functions that differ with respect to the width of the plateau region. Those different potential functions produced almost the same photoproducts, although the isomerization rate was sensitive to the width. The other dihedral parameters of the chromophore remained the same as those in the ground state, except for C₁₅==N and C₁₄C₁₅, for which we used single well potentials with the same curvatures at the trans configurations to inhibit their rotations.

View this table: [in this window] [in a new window]   TABLE 1   Force field parameters for the C12C13==C14C15 torsional coordinate of protonated retinal Schiff base in the S1 excited state*

View larger version (19K): [in this window] [in a new window]   FIGURE 4   Energy profile of the potential function of the C13==C14 dihedral coordinate, 13==14, of the retinal chromophore used in the MD simulations of the photoisomerization. In the electronical ground state, the dihedral potential energy function possesses a minimum at the 13-trans configuration, represented by a dashed line. At the time photoabsorption takes place, the potential function is switched to the reactive one, represented by a solid line, which leads to the isomerization to the 13-cis configuration. There are two possible directions of the isomerization, rotation of the NH bond of the Schiff base toward Asp-85 or toward Asp-212, because the dihedral potential function is symmetric. The vertical energy difference of the dihedral potential functions at the planar 13-trans configuration is set to be 39.1 kcal/mol, which together with the contribution from the chromophore-protein electrostatic interaction, V = 11.2 kcal/mol, reproduces the experimental absorption energy (50.3 kcal/mol).

Recently, the reaction potential energy surfaces of retinal analogues from the Franck-Condon region to the product state in isolated form have been precisely investigated by means of high-level quantum chemical calculations (Garavelli et al., 1997, 1998; Humphrey et al., 1998; González-Luque et al., 2000; Molnar et al., 2000). These studies showed the involvement of the stretching motions of the polyene skeleton in the isomerization process, especially in the initial 100 to 200 fs (the H-to-I process) and during the electronic transition at the conical intersections. The quantum dynamics of the electronic transitions in the protein matrix has also been examined by semiclassical techniques (Humphrey et al., 1998; Ben-nun et al., 1998; Warshel and Chu, 2001), suggesting mechanisms accounting for the high quantum yield (0.64) of the isomerization process in bR. The present classical treatment of the isomerization by a one-dimensional potential function is a rather crude approximation and unable to reproduce prominent characteristics of the photoisomerization dynamics observed in ultrafast spectroscopic studies. Nevertheless, it seems to be unlikely that lack of the details of the actual photodynamics leads to a qualitatively irrelevant structural model of the primary photoproduct, because the conformational change of the chromophore during the isomerization is apparently dominated by the rotation around the C₁₃==C₁₄ bond. Moreover, the quality of the structures of the photoproducts, which may be compromised by the applied simple potential functions of the chromophore, is improved through the QM/MM structure refinements that follow the MD simulations in the present modeling scheme.

It is well known that the charge distribution of the retinal chromophore in bR largely changes both during the photoexcitation and isomerization. Therefore, we used effective charge parameters for the chromophore in the ground and excited states, q and q, respectively, derived from QM/MM complete active space self-consistent field (CASSCF) calculations performed at the ground state energy minimum geometry (Hayashi et al., 2001). These charges are listed in Table 2. The change of charge distribution during the isomerization was then expressed by a function that analytically interpolates q and q smoothly at ₁₃₌₌₁₄ = 90° and 270°,

(1)

(2)

a was chosen to be 0.1, which gives an abrupt switch of the charges around ₁₃₌₌₁₄ = 90° and 270° to mimic the so-called sudden polarization at the conical intersection. It is noted that the change of the charge distribution upon the excitation results in a difference of the chromophore-protein electrostatic interactions, V and, therefore, contributes to the S₁S₀ excitation energy (see QM/MM excitation energy calculations section). In the BR state, the V contribution for this effective charge parameter set was estimated to be 11.2 kcal/mol by our preliminary MD simulations. Therefore, the vertical energy difference of the internal potential functions at ₁₃₌₌₁₄ = 180° was set to be 39.1 kcal/mol (=50.3 kcal/mol (568 nm)  11.2 kcal/mol) as shown in Fig. 4.

Rearrangement of the HBN in the proton channel plays an important role for the pump process. Therefore, a proper description of the HBN is also crucial in our MD simulation. As shown by a previous QM/MM study (Hayashi and Ohmine, 2000), there are strong electronic interactions, such as polarization and charge transfer, between the Schiff base, Asp-85, and Wat-402. The Kitaura-Morokuma energy decomposition analysis (Kitaura and Morokuma, 1976) showed that Wat-402 has a polarization energy of 14.1 kcal/mol and a charge transfer energy of 18.2 kcal/mol due to its interaction with the Schiff base, Asp-85, and other surrounding moieties. These electronic interaction energies are considerably larger than those in water-water interactions; for example, the polarization and charge transfer energies in a water dimer are approximately 1 to 2 kcal/mol (Umeyama and Morokuma, 1977). In our preliminary MD simulations, we found that the TIP3P-AMBER force field, in which nonbonding interactions are described by only classical pair-wise potential functions, cannot even reproduce the conformation of Wat-402 observed in the x-ray ground state structure; after the MD equilibration, the Schiff base NH bond directly attaches to O₂ of Asp-85, and Wat-402 no longer bridges the two. We also observed in the QM/MM optimized ground state structure a strong interaction between NH of the Schiff base and a lone pair orbital of Wat-402 (Hayashi and Ohmine, 2000); the TIP3P water model cannot represent this hydrogen-bond correctly because of the poor description of the lone pairs and electronic interactions.

Fig. 5 compares the interaction energies of Wat-402 with its surroundings in bR, determined in the frame work of the TIP3P-AMBER force field, or by means of the QM/MM calculations. Details of the QM/MM calculations are described in the next section. The energies were estimated for the energy minimum conformation and for 26 other conformations generated by ±15° rotations of the Wat-402 molecule around the Cartesian axes. As clearly seen, the TIP3P-AMBER energies strongly deviate from the QM/MM ones, indicating the poor orientational description of the TIP3P model due to lack of lone pair electrons and electronic interactions.

View larger version (18K): [in this window] [in a new window]   FIGURE 5   Comparison of the interaction energies of Wat-402 with its surroundings, obtained by different empirical force field functions with those obtained by means of the QM/MM method. The energies were evaluated at 27 conformations that were generated by changing orientation of the Wat-402 molecule from the energy minimum geometry. × and  indicate results of the TIP3P and set 1 models (MD simulation section), respectively.

The potential functions describing the internal water molecules were therefore improved as follows. First, the fluctuating charge model (fq-SPC) developed by Rick et al. (1994) was used for all water molecules. This model can take into account the multibody effect due to the electronic polarization induced by the polar environment. Furthermore, we added potential energy terms representing explicitly the angular dependence of hydrogen-bonds and short-range attractions between Wat-402 and the counter-ion aspartate groups. Parameters appearing in the additional potential energy terms were fitted so as to reproduce the QM/MM energies. The reparameterized force field is denoted as set 1. The functional forms and parameters are summarized in an Appendix. Fig. 5 also compares the interaction energies of Wat-402 by using the set 1 force field with the QM/MM values. One can see that the energies of the present model are in good agreement with those of the QM/MM results. Using this model for water, the hydrogen-bond network of the ground state remained stable during a 500-ps MD trajectory. In the simulations of photoisomerization, these additional functions were excluded for the excited state and smoothly switched on during the isomerization with the same switching function used for effective charges (see Eqs. 1 and 2). In the simulations for the modeling of KL, we used the fq-SPC model for Wat-402 but removed the additional functions (see the KL intermediate section). This force field is denoted as set 2.

QM/MM geometry optimization

The early intermediate structures obtained by the MD simulations were refined through QM/MM geometry optimizations. We used a QM/MM method described elsewhere (Hayashi and Ohmine, 2000) that allows one to determine efficiently the optimized geometry of the whole system. One-electron operators based on a restrained electrostatic potential (RESP) method (Bayly et al., 1993) were used to describe the electrostatic interaction between the QM and MM regions, thus, avoiding the time consuming computation of electron integrals of gradient components due to the QMMM interaction during the optimization of the MM region, which usually needs a large number of searching steps. The QM/MM method is implemented in the HONDO (Dupuis et al., 1994) program package. The AMBER force field (Cornell et al., 1995) and the Hartree-Fock (HF) level of theory were used for the description of the MM and QM regions, respectively. For the geometry optimization, the side chains of retinal-Lys-216, Asp-85, and Asp-212, as well as three water molecules (Wat-401, Wat-402, and Wat-406) were included in the QM part. Dunning's split-valence (DZV) basis set (Dunning and Hay, 1977) was used together with diffuse functions on the O atoms of the aspartate side chains. The boundaries dividing the QM and MM regions were set at the CC bond of Lys-216 and the CC bonds of Asp-85 and Asp-212, and the open valences in the QM region were filled by dummy hydrogen atoms following the link atom approach (Field et al., 1990). The interactions between the dummy hydrogen atoms and the MM region were excluded by the RESP treatment (Hayashi and Ohmine, 2000). A 12-Å residue-based cutoff was used for nonbonding interactions.

QM/MM energy decomposition analysis

The QM/MM energies of the BR and intermediate states were analyzed by decomposing them into several contributions. The total energy of the QM/MM system is expressed as

(3)

in which E_QM/MM includes energies of the QM region itself and of the QMMM interaction, and E_MM consists of only interactions within the MM region. Unfortunately, it is rather difficult to consistently determine E_MM of different states in the present scheme because of large fluctuations of polar residues around the surface regions. In this study, therefore, we examined only the E_QM/MM term.

E_QM/MM is decomposed into several terms:

(4)

E⁰ is the QM energy in isolated condition computed at the geometry inside the protein matrix. V and V are electrostatic and Lennard-Jones (LJ) interactions between the QM and MM regions, respectively. V expresses the bonding contributions (bond, angle, and dihedral) at the QM/MM boundary. E is the electronic reorganization (ERO) energy. The external electrostatic field produced by the MM protein matrix polarizes the QM electronic wave function and naturally contributes to V. The induced polarization, however, involves deformation of the electronic wave function, known as ERO, and therefore gives rise to a destabilization effect that to some extent compensates the electrostatic stabilizations (Gao and Xia, 1992; Gao, 1997).

We carried out this analysis for the same QM/MM system used in the geometry optimizations. The QM region includes side chains of retinal-Lys-216, Asp-85, and Asp-212, and three water molecules (see QM/MM geometry optimization section). This QM/MM system is denoted as system A. We also applied the analysis to a QM/MM system in which only the side chain of retinal-Lys-216 is included in the QM region, denoted as system B. In both systems, the QM regions were described by the HF method with the DZV basis set.

QM/MM excitation energy calculations

In the description of the excited state, we employed the QM/MM method using the CASSCF level of theory. The excitation energy (E_ex) was computed for system B (see QM/MM energy decomposition analysis section). Twelve electrons in 11 valence  orbitals of the chromophore were chosen for the CAS space. To improve the flexibility of basis functions, polarization functions were added to the C and N atoms of the conjugated backbone of the retinal Schiff base. The excitation energies were computed by the state-averaged CASSCF calculations for the lowest three -* states (S₀, S₁, and S₂), followed by the state specific CASSCF calculations for the S₀ and S₁ states.

To analyze the molecular mechanism of the characteristic bathochromic shifts of the intermediates, E_ex was also decomposed in the same way as described in QM/MM energy decomposition analysis section. Because V and V do not change upon the excitations, E_ex is expressed by only three contributions (Hayashi et al., 2001),

(5)

Here E is the excitation energy in isolated condition; V is the difference between the chromophore(ch)-protein(pr) electrostatic interaction energies in the S₀ and S₁ states; E is the difference between the ERO energies in those states. Note that upon excitation the positive charge of the chromophore, which is mainly localized in the Schiff base region in the ground state, migrates toward the -ionone ring side of the polyene chain. Hence, the negatively charged counter-ion groups in the vicinity of the Schiff base stabilize the ground state more strongly than the excited state, and, consequently, V is a large positive contribution (Hayashi and Ohmine, 2000; Hayashi et al., 2001).

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSION APPENDIX REFERENCES

MD simulations were performed to examine the photoisomerization and structural changes in the early intermediates, K and KL. Conformational changes of the retinal chromophore and rearrangements of the HBN in the retinal binding pocket revealed by the simulations and QM/MM optimizations are described below. The energetics of the K and KL formations are analyzed. The molecular mechanism of the absorption spectral shifts of the intermediates are also discussed.

Photoisomerization

Using different initial configurations (coordinates and velocities) selected from the ground state equilibrium trajectory at time intervals of 5 ps, 40 isomerization MD trajectories, 20-ps-long each, were calculated. After the photoexcitation, the trajectories stayed around the Franck-Condon region transiently before proceeding to the 13-cis configuration. The average lifetime of the excited state was 116 fs. Because in the present model, the vibrational relaxation of the in-plane skeleton modes of the polyene part contributions to the transition from the Franck-Condon region (H) to a fluorescent state (I) is not described, the lifetime of the excited state is expected to be shorter than the experimental one (the formation time of the vibrationally hot photoproduct, J, is ~500 fs). The lifetime is rather sensitive to the potential profile describing torsion around the C₁₃==C₁₄ bond, and precise molecular models are therefore needed to reproduce the experimentally observed lifetime.

The simulations provided an insight into the directionality of the isomerization, i.e., clockwise or counter clockwise around the C₁₃==C₁₄ bond. The two possible directions of the isomerization result in the rotation of the NH bond of the Schiff base toward Asp-85 or toward Asp-212 (Fig. 4), respectively. In our simulations, all isomerizations happened toward Asp-212. The unidirectional isomerizations toward Asp-212 have been also reported in a recent computational study by Warshel and Chu (2001).

Because a symmetric function has been used for the description of the C₁₃==C₁₄ dihedral potential (Fig. 4), the unidirectionality of the isomerization is purely due to forces of the environment. Our previous MD and QM/MM studies showed that the conformation of the C₁₃ ~ N part of the chromophore in the BR ground state is twisted and the NH bond of the Schiff base is tilted toward Asp-212 (Tajkhorshid et al., 2000; Hayashi and Ohmine, 2000). Particularly, the C₁₂C₁₃==C₁₄C₁₅ dihedral angle deviates considerably (19° in the QM/MM optimized geometry) from planarity. A twisted conformation of retinal toward Asp-212 in BR, which was also observed in the present study, originates from the forces of the environment that favor the isomerization in the direction of Asp-212.

The unidirectional torsion arises from two kinds of interactions of retinal with the surroundings. One is the steric interactions with residues constituting the retinal binding pocket. Because the binding pocket exhibits a slightly bent shape, the conjugated chain of retinal adjusts itself to the binding pocket by twisting around the C₁₃==C₁₄ and C₁₅==N double bonds (Tajkhorshid et al., 2000; Hayashi and Ohmine, 2000). Previous MD simulations (Tajkhorshid et al., 2000) have clearly shown that reduction of the bend of the binding pocket by replacing Trp-86 by Ala results in the relaxation of torsion at those double bonds. The other interaction is due to the strong hydrogen bond of the Schiff base NH bond with Wat-402. The resulting twisted structure of the chromophore could not be observed when Wat-402 was described by the TIP3P model, instead of the present water model that includes three body effects and explicit interaction of lone pairs. Therefore, strong hydrogen-bonds between the Schiff base, Wat-402, and Asp-85 are essential for the twisted retinal conformation that ensures a unidirectional photoisomerization.

The K intermediate

To proceed to the first intermediate, eight isomerization trajectories were continued for 200 ps. All runs resulted in almost the same 13-cis photoproduct, in which the chromophore had a strongly deformed structure. Fig. 6 b displays the QM/MM energy minimized structure of this photoproduct, assigned to the K (or K_LT) intermediate. One can clearly see that the chromophore is largely twisted in the C₁₃ ~ N region, and the NH bond of the Schiff base is almost perpendicular to the membrane normal and pointing toward Asp-212. In a fully planar 13-cis conformation of the chromophore, the NH bond is expected to point "up" (parallel to the membrane normal). However, in none of the eight simulations performed at room temperature was a completely relaxed planar 13-cis conformation observed. The NH bond maintains its interaction with Wat-402 and Asp-212 in the K intermediate, thus preventing the chromophore from relaxing into a planar conformation.

View larger version (25K): [in this window] [in a new window]   FIGURE 6   Structural changes in the vicinity of the Schiff base in the early phase of the photocycle, modeled in the present study. Stereo views are shown for the QM/MM optimized structures of (a) BR, (b) K, and (c) KL. The hydrogen-bonds are depicted as dashed lines. Wat-402 in the calculated KL structure corresponds closely to Wat-401 in the x-ray structure model of K (Edman et al., 1999), and Wat-401 in the present model is suggested to be the missing water in the x-ray model (see the KL intermediate section).

Table 3 compares geometrical parameters of the chromophore in the QM/MM optimized structures of BR and K. In K, the strong torsions are mainly localized in the C₁₃ ~ N region. Especially, the C₁₄C₁₅==NC and C₁₂C₁₃==C₁₄C₁₅ double bonds are largely twisted by 30° and 25°, respectively. Although dihedral angles of the other half of the polyene chain, C₁₃ to C₅, do not significantly deviate from 180°, small rotations around C₉==C₁₀ and C₁₁==C₁₂ result in an overall twist in this half. Consequently, the C₂₀ methyl group attached to C₁₃ is tilted toward Asp-212, and the angle of the C₁₃C₂₀ bond to the membrane normal in K differs by 30° from that in the BR ground state (for the numbering of retinal carbons, see Fig. 1).

Fig. 7 a schematically depicts the HBN in the K intermediate. The HBN does not undergo large conformational changes upon the formation of K, as can be discerned by comparison of Figs. 3 (BR) and 7 (K). Only minor changes in the hydrogen-bond distances were found, except for the hydrogen bonds between H of the Schiff base and Wat-402. The isomerization of the chromophore moves the H atom upward, resulting in a weakened hydrogen bond between the Schiff base and Wat-402. The distance between Wat-402:H₂ and Asp-85:O₁ also slightly increases due to the reduction of charge transfer and polarization effects from the Schiff base.

View larger version (18K): [in this window] [in a new window]   FIGURE 7   Schematic representations of the HBN in the binding pocket in (a) K and (b) KL. Dashed thin lines represent hydrogen-bonds. Distances of the hydrogen-bonds are shown in Å. Dashed thick arrows in gray illustrate conformational changes of the hydrogen-bond network during the K-to-KL transition.

The strongly twisted chromophore in our K intermediate presents a remarkable difference from the x-ray crystallographic model of Edman et al. (1999), in which the chromophore conformation is assumed planar. Furthermore, the K intermediate in the present study involves only minor changes of the HBN, as seen in Fig. 6, whereas in the x-ray model (Fig. 2 b) significant changes of the HBN are observed as discussed in the Introduction section. In the present K model, Wat-402 occupies almost the same position as in BR, and the hydrogen-bond between Thr-89 and Asp-85 is not broken.

The present K model is, however, in keeping with the results of LT-RR (Pande et al., 1981; Braiman and Mathies, 1982) and LT-FTIR (Bagley et al., 1982; Rothschild and Marrero, 1982) spectroscopic measurements, which indicate a strong enhancement of the hydrogen-out-of-plane (HOOP) vibrational peaks of the chromophore upon the formation of K. One of the prominent peaks at ~960 cm¹ has been assigned to the 15-HOOP mode (Braiman and Mathies, 1982). The intensities of the HOOP modes are small in the planar conformation and become large when the structure is distorted. Hence the torsional deformation around the Schiff base in our K model is consistent with the large HOOP peaks in the LT-RR and LT-FTIR spectra.

The IR spectral shape of the HOOP peak at 960 cm¹ is altered by the mutation of Thr-89 to Cys (Kandori et al., 1999). The present model of K can explain this mutation effect. The H₁₅ atom is close to Thr-89:O (2.58 Å), and replacement of the O atom by the S atom, which is larger, in the T89C mutant can affect the spectral shape of the 15-HOOP. Contrary to the x-ray model (Edman et al., 1999), FTIR measurements (Kandori et al., 1998b, 1999, 2001) indicate that the hydrogen-bond between Thr-89 and Asp-85 becomes stronger upon the formation of K. The distance between Thr-89:H and Asp-85:O in the simulated K intermediate is 0.01 Å shorter than that in BR (Figs. 3 and 7 a), implying a stronger hydrogen-bond, and thus supporting the FTIR result. The almost perpendicular orientation of the NH bond of the Schiff base to the membrane normal is also in good agreement with a recent polarized FTIR measurement (Kandori et al., 2002).

The KL intermediate

As discussed above, the chromophore has a strongly twisted conformation in K because the Schiff base maintains its connection to Wat-402. Nevertheless, the hydrogen-bond is rather weak in K and the distance between H of the Schiff base and Wat-402:O (2.7 Å) is remarkably larger compared to other hydrogen bonds (1.5 ~ 1.8 Å), as seen Fig. 7 a. Furthermore, during the simulations, the distance between the Schiff base and Wat-402 frequently becomes longer due to thermal fluctuation and tension by the strongly deformed chromophore in the Schiff base region. Breaking of the hydrogen-bond between the Schiff base and Wat-402 seems to be a prerequisite for the relaxation of the chromophore to a planar conformation in which the NH bond points "up," i.e., toward the cytoplasmic side. Although none of the isomerization trajectories applying the set 1 force field used in the modeling of the K intermediate resulted in the complete cleavage of this hydrogen-bond, it is likely that the dissociation of the hydrogen bond leads to a transition from K to the next intermediate, KL.

Because the lifetime of K can be as long as 100 ns according to a recent TR-FTIR study (Doiumaev and Braiman, 1997), it is not feasible to observe the decay of K in MD simulations. Furthermore, the lifetime of K can increase even further due to the strengthened hydrogen bond between the Schiff base and Wat-402 in the set 1 force field. Unfortunately, it is rather difficult to determine empirical potential functions that can quantitatively describe the dissociation of hydrogen bonds that involve strong electronic interactions. To accelerate the decay of the K intermediate, therefore, we used the set 2 force field, which does not include the additional function for the hydrogen bond (see MD simulations section), and repeated the simulations to examine the dynamics of the photoproducts when this hydrogen bond can be easily broken.

Starting from the initial configurations taken from the BR equilibrium trajectory at time intervals of 10 ps, 32 isomerization trajectories, 80-ps-long each, were computed. The additional hydrogen-bond function in the set 1 force field was turned off at the beginning of isomerizations. The absence of the strong hydrogen bond did not affect the dynamics of the isomerization; all isomerizations took place toward Asp-212, as seen in the simulations of the K intermediate. In addition, the first photoproduct is also similar to the K intermediate obtained in the simulations with the set 1 force field, i.e., the chromophore is strongly deformed in the Schiff base region, and the NH bond is perpendicular to the membrane normal.

After the formation of the K intermediate characterized as the first photoproduct, however, approximately two-thirds of the trajectories (21 in the 32 trajectories) proceeded to a second stable photoproduct in 80 ps. The remaining 11 trajectories stay in K on this time scale. Because the transition is a stochastic process, the 11 trajectories should reach the same second stable state after a longer time.

The excitation energy of the second photoproduct is red shifted from BR and blue shifted from K (see below). Therefore, we assign this state to the KL intermediate. Fig. 8 shows the QM/MM optimized structure of the KL intermediate. The K-to-KL transition is characterized by significant structural changes in the retinal binding pocket, as illustrated schematically in Fig. 7. The first event is a complete dissociation of the hydrogen bond between the Schiff base and Wat-402, which persists in the K intermediate. The dissociation leads to conformational relaxation of the chromophore to the 13-cis planar form, where the NH bond of the Schiff base points to the cytoplasmic side and starts to establish a weak hydrogen-bond with Thr-89:O. Along with formation of the hydrogen bond, the Schiff base part of the chromophore becomes twisted again but in the opposite direction to the K intermediate. One can discern this whiplash motion in Fig. 6 c, which displays the retinal binding site in KL, the hydrogen-bond between the Schiff base and Thr-89:O and the twist in the Schiff base region. As will be discussed in the next section, the distorted conformation of the Schiff base is stabilized by electrostatic interaction with Asp-85, Asp-212, and Thr-89.

View larger version (34K): [in this window] [in a new window]   FIGURE 8   Structures of the K and KL intermediates. The seven transmembrane helices are shown as cylinders; loops are shown in tube representation; -strands are drawn as ribbons. The retinal chromophore and key amino acids are highlighted (K is shown in gray, blue, and red; KL is shown in green). Water molecules (W-401, W-402, and W-406) are represented as red (K) and green (KL) spheres.

The geometrical features the chromophore in the KL intermediate are listed in Table 3. The C₁₂C₁₃==C₁₄C₁₅ dihedral angle deviates by 26° from the planar conformation. Although the C₁₄C₁₅==NC and C₁₃==C₁₄C₁₅==N dihedral angles are also twisted considerably, the torsions are somewhat relaxed compared with that in the K intermediate. The relaxed torsions are consistent with the LT- and TR-FTIR and TR-RR measurements, where the 15-HOOP mode in KL (~980 cm¹) exhibits a high frequency shift from that in the K intermediate (960 cm¹) (Rothschild et al., 1985; Sasaki et al., 1993; Weidlich and Siebert, 1993; Hage et al., 1996; Dioumaev and Braiman, 1997; Althaus et al., 1995). As previously shown (Tajkhorshid et al., 1999; Tajkhorshid and Suhai, 1999), the H₁₅ atom carries a large positive partial charge. The strong hydrogen bond between H₁₅ and Asp-212:O with 2.52 Å distance further stabilizes the twisted conformation in the Schiff base region. This hydrogen bond can, therefore, also contribute to the aforementioned high frequency shift of the 15-HOOP mode, because the interaction leads to a steeper potential well along a librational motion of the C₁₅H₁₅ bond that corresponds to the HOOP motion.

One also notices in Fig. 6 c and 8, a large movement of the Lys-216 chain during the K-to-KL transition. In particular, the NC single bond rotates by 90.2° (Table 3), and the C atom moves down (Fig. 6 c). The conformational change of the chromophore-Lys-216 moiety causes a significant rearrangement of the HBN in the binding site. As clearly seen in Fig. 6 c, Wat-402 is dislocated from its position in the BR and K states (between the Schiff base, Asp-85, and Asp-212). Due to the complete dissociation of its hydrogen bond with the Schiff base, and steric interaction with the Lys-216:C, Wat-402 is forced to move toward the extracellular side, initiating an overall rearrangement of the HBN. Wat-401 moves and bridges between Asp-85 and Asp-212, and, as a result, the carboxylate group of Asp-85 rotates slightly. The dislocation of Wat-402 also triggers a downward movement of Arg-82 (1.5-Å displacement at Arg-82:C) during the K-to-KL transition, as clearly seen in Fig. 8. A remarkable downward shift of the G helix, induced by the conformational change of Lys-216, is also observed in the KL structure.

Several structural features discussed above for the KL intermediate are actually consistent with the x-ray crystallography model of K (Edman et al., 1999). The dislocation of Wat-402 is in good agreement with the strong negative density in the KBR difference electron density map. In addition, a significant movement of the Lys-216 side-chain around the C atom and the downward shift of the G helix are common to both models. Furthermore, the carboxylate group of Asp-85 in our KL model is rotated as observed in the x-ray K model, although in our KL model this rotation is smaller and the hydrogen bond between Thr-89 and Asp-85 is not broken. Therefore, the x-ray model closely corresponds to the KL intermediate obtained in this study.

The x-ray measurements have been carried out for an intermediate state trapped at 110 K, a temperature 33 K higher than the liquid nitrogen temperature (77 K) used for trapping the K intermediate in spectroscopic studies. This slightly higher temperature might give rise to a mixture of the K and KL states. According to our model, the K formation involves only minor structural changes of the HBN. Therefore, it is likely that the observed difference electron density map is mainly due to structural differences between the BR and KL states. The present models show that the conformations of the chromophore around the Schiff base in K and KL are remarkably different (see Fig. 6, b and c, and Table 3). The inhomogeneity due to the presence of the two K-like species can also account for the absence of a positive density around the Schiff base region of the chromophore in the x-ray difference electron density map (Edman et al., 1999).

The present KL model also indicates the location of Wat-402, which has not been identified in the x-ray K model. According to our model, Wat-402 moves to a position that closely corresponds to the position of Wat-401 in the x-ray model of K. Therefore, the missing water in the x-ray model is the water bridging between Asp-85 and Asp-212 in our KL model (Wat-401 in Fig. 6 c). In fact, there seems to be a positive density at this position in the x-ray difference electron density map (see Fig. 5 in Pebay-Peyroula et al., 2000), although no water molecule is modeled at this position in the x-ray structure (PDB code, 1QKO).

The downward movement of Arg-82 toward the extracellular side in our KL model is not reported in the x-ray crystallographic study. This disagreement may be due to the different temperatures used for generating the intermediate state. In the x-ray study, the low temperature (110 K) was maintained throughout the experiments, whereas in the present study the KL state was obtained by simulations at room temperature (300 K). Therefore, in the simulations we may observe structural fluctuations that cannot be detected in low temperature experiments. It is noted that the downward movement of Arg-82 has been observed in x-ray crystallographic structures of the later states, L (Royant et al., 2000) and M (Luecke et al., 1999b, 2000; Sass et al., 2000; Facciotti et al., 2001). The movement in the present KL model is, however, much smaller than the motions observed in the x-ray structures of L and M. This small movement is induced by the HBN rearrangement, especially the downward migration of Wat-402, upon the formation of KL.

Energetics of the early intermediates

As seen in the previous sections, the chromophore experiences significant conformational changes during the formation of the early intermediates, resulting in the HBN rearrangement in the retinal binding pocket. Understanding the energetics related to these conformational changes provides an insight into the molecular mechanism of the energy storage in the early phase of the pump cycle. We therefore decomposed and analyzed the QM/MM energies, E_QM/MM (Eqs. 3 and 4), calculated for the early intermediate states. E_QM/MM is defined as the sum of the energies of the QM region itself and its interaction with the MM region. Details of the method are described in the QM/MM energy decomposition analysis.

Table 4 lists the differences in E_QM/MM between K, KL, and BR. The difference in E_QM/MM between K and BR for system A is 16.2 kcal/mol. Because the MM region in system A stays almost unchanged upon the formation of K (see the K intermediate section), energy of the MM region, E_MM, in K is expected to be almost the same as in BR. Thus, an enthalpy change, i.e., energy change of the whole system, upon the formation of K can be well described by the energy difference in E_QM/MM alone. The enthalpy change estimated in the present study is in good agreement with the experimental value of 16 kcal/mol (Birge and Cooper, 1983; Birge et al., 1989).

View this table: [in this window] [in a new window]   TABLE 4   Changes of QM/MM energies (EQM/MM*) upon the formation of K and KL, and their additive contributions, E0, V,§ V,¶ V, and E** in kcal/mol

To calculate the contribution of the chromophore's conformational change alone, we also evaluated E_QM/MM for system B, in which only the chromophore is included in the QM region, and Asp-85, Asp-212, and water molecules of the binding pocket reside in the MM region. The difference in E_QM/MM between K and BR for system B is 23.1 kcal/mol, which is 6.9 kcal/mol larger than for system A. The difference is due to the slight changes in the structure and electronic interaction of nonretinal groups participating in HBN (Figs. 3 and 7).

Table 4 also summarizes the decomposed contributions to E_QM/MM for system B (see QM/MM energy decomposition analysis section). The difference in E⁰ is due to the conformational change of the chromophore. The strong torsion around the Schiff base region in K accounts for the large E⁰ difference between K and BR (11.9 kcal/mol). Change of the chromophore-protein interactions also contributes to E_QM/MM. The electrostatic stabilization V in K is significantly smaller (23.3 kcal/mol) than in BR. Upon the formation of K, the H atom of the Schiff base moves away from such groups as Asp-85, Asp-212, and Wat-402, resulting in the weakening of the electrostatic stabilization. The changes of the electrostatic stabilization effects of Asp-85, Asp-212, and Wat-402 are estimated to be 4.4, 5.6, and 16.5 kcal/mol, respectively. On the other hand, the LJ interaction energy, V, significantly decreases (8.5 kcal/mol) upon the formation of K. This drastic decrease is mainly due to the reduced strong short-range repulsion between the Schiff base and Wat-402 (8.8 kcal/mol). The ERO energy, E, also gives a stabilizing contribution (3.4 kcal/mol), consistent with the large destabilizing contribution of V (Sec. 2.4).

Altogether, upon the formation of K, the energy (16.2 kcal/mol) is stored in the form of the conformational distortion of the chromophore (E⁰ = 11.9 kcal/mol) and the weakening of the interaction of the Schiff base with the surrounding polar residues (V + V + E = 11.4 kcal/mol), part of which is compensated by the structural and electronic changes of the HBN in the binding pocket (6.9 kcal/mol). This scheme is comparable with the results of Birge and Cooper (1983) suggesting that roughly one-half of the energy is stored in the form of the changes of the electrostatic interaction.

The analysis of energetics for KL is also summarized in Table 4. The difference in E_QM/MM between KL and BR for system A is 29.0 kcal/mol, showing an energy increase by 12.8 kcal/mol from K to KL. This apparent endothermic change is due to the fact that the energy change of the MM region, which experiences conformational changes during the K-to-KL transition, is not included in our analysis; E_QM/MM does not include some of the energy changes caused by the extensive rearrangement of the HBN, e.g., the movement of Arg-82 toward the extracellular side in the K-to-KL process (see the KL intermediate section). In fact, during the K-to-KL transition, the side chain of Arg-82 was found to gain large electrostatic stabilizations from Glu-194 (10.9 kcal/mol) and Thr-205 (3.4 kcal/mol) located at the extracellular end of the channel in the MM region, and those energies are not included in E_QM/MM. After taking into account such effects, it becomes evident that, despite the apparent increase of E_QM/MM (12.8 kcal/mol), the energy of the KL intermediate is not higher than K.

The increase of E_QM/MM from BR to KL in system B was estimated to be rather small (4.5 kcal/mol) compared with increase from BR to K (Table 4). The increasing stabilization of the chromophore during the K-to-KL transition implies that the energy, which was mainly stored in the chromophore itself and its interaction with the surroundings in K, is now partially used to rearrange the HBN in the extracellular channel, as seen in the KL intermediate section. The contribution of the conformational distortion of the chromophore, E⁰, is still large (7.1 kcal/mol) in KL, but 4.8 kcal/mol smaller than K. This is consistent with the relaxed torsion of the Schiff base group in KL, as seen in Table 3. Compared with the K intermediate, the chromophore-protein electrostatic interaction (6.7 kcal/mol) is significantly smaller in KL. Although the hydrogen bond between the Schiff base and Wat-402 completely dissociates during the transition of BR to KL and the chromophore loses large electrostatic stabilization energy (20.8 kcal/mol), the loss is compensated by stabilization effects of Asp-85 (10.0 kcal/mol), Thr-89 (3.4 kcal/mol), and Asp-212 (