INTRODUCTION

Top Abstract Introduction Methods Results Discussion Conclusions References

The protein bacteriorhodopsin (bR) resides in the membrane of the archaebacterium Halobacterium salinarium and utilizes sunlight to drive a transmembrane proton pump in the most basic form of photosynthesis known (Lozier et al., 1975). The 26-kDa protein incorporates a retinal chromophore bound to a lysine residue via a protonated Schiff base linkage and absorbs light around 568 nm. Photoexcitation triggers an isomerization of retinal about its C₁₃==C₁₄ double bond from an initial all-trans to a 13-cis configuration, which is completed within a few picoseconds. The photoreaction induces a vectorial transfer of a proton, leading to the release of a proton to the extracellular side and uptake from the cytoplasmic side. The structure, photocycle, and previous molecular dynamics studies of bR have been reviewed in Oesterhelt et al. (1992), Ebrey (1993), Schulten et al. (1995), and Lanyi (1997).

The bR photocycle proceeds through several sequential intermediates characterized by distinct retinal absorption maxima. The current view of the photoisomerization process, which has emerged from time-resolved spectroscopy and modeling of the bR photointermediates (Mathies et al., 1988; Dobler et al., 1988; Du and Fleming, 1993; Humphrey et al., 1995, 1997; Xu et al., 1996), suggests that the reaction proceeds through the sequential steps

(1)

where photoexcitation of bR generates bR*, which successively decays to the intermediates J₆₂₅ and K₅₉₀, the latter being a ground-state product with retinal in a 13-cis geometry. The observations suggest that substantial torsional motion takes place during the ~200-fs lifetime of the excited-state bR*.

A schematic structure of bR together with retinal bound as a protonated Schiff base to Lys²¹⁶, water molecules, and key amino acid side groups Asp⁸⁵, Asp⁹⁶, Asp²⁰⁴, and Asp²¹² is shown in Fig. 1. Asp⁸⁵ is unprotonated and negatively charged at the beginning of the proton pump cycle; this group has been shown in mutagenesis (Mogi et al., 1988; Subramaniam et al., 1990; Otto et al., 1990) and Fourier transform infrared spectroscopy experiments (Braiman et al., 1988) to accept a proton from the Schiff base after photoisomerization. The D85N mutant forms the K₅₉₀ intermediate in Eq. 1, but exhibits no proton pumping activity (Mogi et al., 1988; Subramaniam et al., 1990; Otto et al., 1990). Similarly, the D212N mutant involving another unprotonated aspartic acid side group near retinal fails to pump protons (Mogi et al., 1988; Otto et al., 1990). After photoexcitation, the excited-state dynamics of these mutants differs markedly from the behavior for native bR. The fluorescence lifetime for D85N is 20 times longer than the 500-fs lifetime of native bR, and for D212N is 4 times longer (Song et al., 1993). The quantum yields of these mutants, nevertheless, are close to the wild-type value of 0.64 ± 0.04 (Schneider et al., 1989; Govindjee et al., 1990; Logunov et al., 1996a, b).

View larger version (71K): [in this window] [in a new window]   FIGURE 1   The protein bacteriorhodopsin, including retinal, water molecules, and residues Lys216, Asp85, Asp96, Asp212, and Glu204, which participate in proton translocation. Helixes A, B, E, F, and G are drawn as cylinders; helices C and D are drawn as thin tubes to reveal the protein interior. This figure was created with VMD (Humphrey et al., 1996).

Previous simulations of the photoisomerization of bR have employed a single potential surface for the excited state, evaluated by means of QCFF/PI methods (Warshel et al., 1991), or have been described by a simple periodic function (Zhou et al., 1993; Humphrey et al., 1995). However, polyenes are known to have two closely lying singlet excited states: S_single, which, for pure polyenes, is optically allowed and involves mainly single excitations; and S_double, which is optically forbidden and involves mainly two-electron excitations (Schulten and Karplus, 1972; Hudson et al., 1982; Tavan and Schulten, 1986; Serrano-Andreas et al., 1993a). In polyenes and unprotonated Schiff bases of retinal, S_double is the lowest state in energy, whereas for protonated Schiff bases S_single is the lowest state (Schulten et al., 1980). Photoexcitation of bR populates mainly S_single, i.e., the lowest excited singlet state; the higher state S_double is responsible for the two-photon absorption reported by Birge and Zhang (1990).

The question arises, how much does a torsion around retinal's C₁₃==C₁₄ bond during photoisomerization alter the energy level ordering of S_single and S_double such that the state S_double becomes involved in the photodynamics of retinal? This question must be pursued because omission of the second, i.e., higher energy, excited-state S_double in such a case could affect the photoreaction. In fact, nonadiabatic torsion around the C₁₃==C₁₄ double bond promotes the two -electrons residing in the ethylenic ground state of this bond to a double excitation, which argues for an involvement of the state S_double in the all-trans  13-cis reaction.

Accordingly, we employ in the present study model potential surfaces involving three electronic states, as suggested by Schulten et al. (1995). The three-state model implies that the system is required to pass two crossing regions to return, after photoexcitation, to the ground state. Level crossing is a quantum mechanical process and requires a quantum mechanical description. In this study we will employ for this purpose the so-called density matrix evolution method, also known as mean field approach, in which a 3 × 3 density matrix accounts for the occupancy of the three participating electronic states and the nuclear motion is represented in an approximate fashion through a single classical trajectory (Mavri et al., 1993). More exact description is provided in Ben-nun et al. (1998).

In the following we introduce first the quantum chemical and combined quantum/classical mechanical simulation methods employed in our study. We present then the simulation results accounting for participation of three electronic states in the photodynamics of retinal in wild-type bR, and for the D85N and D212N mutants of bR. From these simulations we compute quantum yields and excited-state lifetimes.

	METHODS

Top Abstract Introduction Methods Results Discussion Conclusions References

The potential surfaces governing photoisomerization processes have been evaluated in vacuo for the retinal analog [CH₂---(CH)₃---(C₂H₃)---(CH)₂---NH---CH₃]⁺. This retinal analog is known to be the shortest analog that exhibits proton pumping activity upon reconstitution with the apoprotein (Steinberg et al., 1991). The choice of this analog was dictated by the fact that retinal itself, which contains six double bonds, proved too large for a sufficiently extensive calculation.

The electronic excitations and potential surfaces of conjugated polyenes have been described by means of extended ab inito calculations by a number of groups. Previous studies have focused on the electronic excitations of octatetraene (Serrano-Andreas et al., 1993b) and retinal (Du and Davidson, 1990), as well as on the ground- and excited-state potential surfaces of butadiene, hexatriene (Ollivuci et al., 1993, 1994), and a polyene Schiff base of two double bonds (Bonacic-Koutecky et al., 1987). In this investigation we extend the earlier calculations of potential surfaces to a polyene Schiff base of four double bonds, i.e., eight -electrons. We employ the program MOLPRO, which has been specifically designed to provide fast and accurate ab initio calculations of excited-state potential surfaces (Werner and Knowles, 1985; Knowles and Werner, 1985).

Calculations were carried out at the CASSCF(8,8)/6-31G level. After an initial self consistent field (SCF) calculation to determine all orbitals, the procedure chosen further optimized the -orbital, allowing all possible  electron configurations (in the eight  orbitals). Such calculations are difficult because of convergence problems and were only feasible because of the advanced optimization algorithms of MOLPRO (Werner, 1987). The MOLPRO program employs an efficient second-order Monte Carlo SCF method for long configuration expansions. The program also provides the feature of so-called state-averaged calculations in which the average energy of the states under consideration is optimized; this feature was used extensively for the present study.

As argued above, one expects the three lowest electronic states (S₀, S₁, and S₂) to be implicated in retinal's photoisomerization. To determine the overall shape of the potential surfaces (cf. Fig. 2), we averaged over the three electronic states S₀, S₁, S₂. In regions of close approach between any two surfaces, additional state-averaged calculations were performed to investigate the nature of the "crossing" between the respective states; in such calculations molecular orbitals were optimized with respect to the average energy of the two particular states of interest (either S₀ and S₁, or S₁ and S₂).

View larger version (21K): [in this window] [in a new window]   FIGURE 2   Computed ground- and excited-state potentials of [H2C==(CH)---(CH)==(CH)---(CH3---C*)==(C*H)---(CH)==NH+---CH3], an analog of the protonated Schiff base of retinal. The dependence of the energies of the ground state (S0) as well as the first (S1) and second (S2) excited states on the torsional angle  of the bond C*==C* is shown. The bond C*==C* corresponds to the C13==C14 double bond of retinal. State-averaged calculations for the S0, S1, and S2 electronic states were performed at the CASSCF(8,8)/6-31G ab initio level (see text). 1 and 2 indicate the two regions of close interaction between the two excited states, and between the first excited state and the ground state, respectively.

The initial geometry of the retinal analog was minimized at the Hartree-Fock level, using a 6-31G basis set. The potential surfaces were computed for a gradual torsion in steps of 15° around the C==C double bond adjacent to the Schiff base of the retinal analog; all other degrees of freedom were kept constant. Basis set improvements beyond the 6-31G level were found not to significantly affect the potential surfaces.

From the adiabatic surfaces in Fig. 2 we change in our simulations to a presentation in terms of a -dependent 3 × 3 Hamiltonian matrix

(2)

where  represents the torsional angle about the C₁₃==C₁₄ double bond (the isomerization reaction coordinate), E_j() (j = a, b, c) are nonadiabatic potential surfaces, and V_ij are the nonadiabatic coupling constants, which we choose to be real and -independent. As discussed further below, the nonadiabatic surfaces E_j and coupling constants V_ij were chosen such that the diagonalized form of H() in Eq. 2 approximates the ab initio surfaces of the states S₀, S₁, and S₂ in Fig. 2. In the new representation the regions of nonadiabatic interaction between the states S₀, S₁, and S₂ are represented by explicit crossing points along the E_j() surfaces.

The use of a truncated retinal analog and the neglect of the stabilizing effect of the protein environment resulted in an overestimate of the S₀  S₁ excitation energy. This deficiency was corrected through a linear scaling of E_a, E_b, and E_c such that the energy difference between the states S₀ and S₁ at the all-trans position measured 50 kcal/mol, a value close to the excitation energy of retinal in bR. The resulting nonadiabatic surfaces E_j() are shown in Fig. 3.

View larger version (19K): [in this window] [in a new window]   FIGURE 3   Nonadiabatic potential energy surfaces used in molecular dynamics simulations of bacteriorhodopsin and bR mutants. The minimum of the surfaces Ec,bR, Ec,D85N, and Ec,D212N are matched to the measured excitation energy of bR and its respective mutants.

The minimum energy separation between S₀ and S₁, as well as between S₁ and S₂, allows one to determine the coupling constants V_ij, i.e., the coupling constants should be half of the minimum energy splitting. After state average calculations (described in the next section) for the ₂ region in Fig. 2, the splitting between S₁ and the ground state is 0.92 kcal/mol, and the splitting between S₁ and S₂ is 4.4 kcal/mol. These values correspond to coupling constants of 0.47 kcal/mol for V_ab and 2.2 kcal/mol for V_ac and V_bc. For simplicity, V_ab was chosen to be 0.5 kcal/mol, and V_ac (as well as V_bc) was chosen to be 1 kcal/mol and 3 kcal/mol in two sets of simulations. It is worth mentioning here that the quantum yield depends sensitively on V_ab: a value of V_ab = 0.5 kcal/mol within the density matrix evolution description adopted here leads to a quantum yield that is close to the experimentally observed value 0.64 (Schneider et al., 1989; Govindjee et al., 1990), whereas a choice of V_ab = 1 kcal/mol would result in a quantum yield of 0.34. Simulations showed that the time for E_c  E_a does not depend sensitively on V_ac and V_bc: a value of V_ac = V_bc = 1 kcal/mol yields a crossing time of 332 fs, and a value of V_ac = V_bc = 3 kcal/mol yields a crossing time of 300 fs. Only simulation results from a coupling constants set (V_ab, V_ac, V_bc) = (0.5 kcal/mol, 1.0 kcal/mol, 1.0 kcal/mol) are shown in the next section.

Hamiltonians similar to Eq. 2 were constructed for the bR mutants D85N and D212N absorbing at 615 nm (Mogi et al., 1988) and at 585 nm (Needleman et al., 1991), respectively. For this purpose the minima of the E_c surface at the all-trans and 13-cis positions were lowered to match the respective excitation energies, whereas E_a, E_b, and the couplings V_ij were assumed to be the same as for wild-type bR. The resulting surfaces E_c,bR, E_c,D85N, and E_c,D212N are compared in Fig. 3.

The photodynamics simulations used as a starting point the refined bR₅₆₈ structure reported by Humphrey et al. (1994), which is based on a structure of bR obtained by electron microscopy (Henderson et al., 1990). The refined structure includes the loop regions between the seven transmembrane helices of bR, as well as 16 water molecules within the protein interior. Five of these water molecules are located within the vicinity of retinal's binding site, in hydrogen-bond contact with the Schiff base proton and nearby amino acids, including Asp⁸⁵ and Asp²¹². As shown by Stuart et al. (1995), the Schiff base binding site is either neutral or positively charged. In the present structure, the counterion of retinal is a hydrogen-bonded complex involving several water molecules, in which Asp²¹², Asp⁸⁵, Tyr⁵⁷, Tyr¹⁸⁵, Arg⁸², and Thr⁸⁹ correspond to a neutral binding site region (Humphrey et al., 1994). The possibility of calcium binding sites near the Schiff base suggested by Stuart et al. (1995) seems to be ruled out by recent x-ray diffraction measurement of crystallized bR (H. Luecke, personal communications). The structure in the current simulations has a retinal binding site similar to bR structures recently reported (Grigorieff et al., 1996; Kimura et al., 1997; Pebay-Peyroula et al., 1997). Structures for bR D85N and D212N mutants were constructed through modification and equilibration of the structure of native bR as previously described by Humphrey et al. (1997). A series of simulations was carried out for native bR as well as for its D85N and D212N mutants. Each simulation started with a different set of (random) velocities and a brief 100-fs equilibration phase, which was followed by simulation of the photoisomerization reaction employing the density matrix evolution (DME) method (Berendsen and Mavri, 1993) based on the potential energy surfaces in Fig. 3. The photoisomerization simulations continued until retinal returned to the ground state, and overall simulation periods varied from 1 to 12 ps. All simulations that completed an all-trans 13-cis isomerization were continued for a further 2 ps to model the J₆₂₅  K₅₉₀ transformation, assumed here to be a relaxation process (Xu et al., 1996); for this purpose the ground-state potential introduced by Humphrey et al. (1995),

(3)

with k_13-14 = 20 kcal/mol, was used to describe the C₁₃==C₁₄ torsion, replacing the three-level potential in Eq. 2.

The DME description of the torsion around the C₁₃==C₁₄ bond couples the evolution of the quantum mechanical density matrix of the three electronic states involved in the photoisomerization to a single classical trajectory of the entire protein (Berendsen and Mavri, 1993). Using the three-state Hamiltonian H in Eq. 2, we introduce the density matrix (t), the diagonal elements of which represent the occupation of the three states as a function of time. Photoexcitation is accounted for by starting from the pure state,

(4)

During the subsequent simulation, the density matrix evolves according to the Liouville-von Neumann equation,

(5)

The density matrix incorporates the effect of coherence between the quantum states during the dynamics, without the need to calculate explicitly the wave function for the three electronic states.

The DME description is only approximate, as discussed further below. More exact schemes employed in combined quantum/classical molecular dynamics simulations are the surface-hopping algorithm (Tully and Preston, 1971; Tully, 1990) and the multiple spawning algorithm (Martinez et al., 1997; Ben-Nun and Martinez, 1998), which represent the nuclear motion in a less approximate way. These schemes, in the context of bacteriorhodopsin, are extremely expensive computationally, because they require statistical ensembles involving large numbers of trajectories. Simulations of the bR photodynamics using the multiple spawning algorithm are reported in Ben-Nun et al. (1998).

The energy stored in the torsion around the C₁₃==C₁₄ bond, in the framework of the DME description, is given by

(6)

The corresponding forces are obtained after taking the (negative) gradient with respect to atomic coordinates. Because the couplings V_jk are chosen to be independent of , the forces acting on atom k with position _k are

(7)

where _jk = E_j()/_k, i.e., the net force acting on an atom is the average over the forces due to each state, weighted by the occupancy _jj of that state. To the resulting force on atom k are added all other forces acting on the atom from the surrounding protein environment.

Equation 5 was integrated by means of a fourth-order Runge-Kutta scheme to yield (t); rapid oscillations in the occupancy of the states required a 0.1-fs integration time step. For the classical molecular dynamics simulations, the program NAMD (Nelson et al., 1996) was used (after incorporation of the DME method) for integration of the classical equations of motion, and a 1-fs integration time step was employed. The simulations were carried out at a temperature of 300 K. In carrying out the simulations the CHARMm force field (Brooks et al., 1983) was used, except for the torsional motions around the retinal backbone. For the latter we employed torsional potentials as determined by Logunov and Schulten (1996) and provided in Table 1, except for the C₁₃==C₁₄ bond torsion, the potentials of which are characterized in Table 2.

View this table: [in this window] [in a new window]   TABLE 1   Barriers ki in the potential function Eidihe = 1/2 ki[1 + cos(2i + i)] governing torsions around bonds of retinal's backbone, except the 13==14 bond

View this table: [in this window] [in a new window]   TABLE 2   Barriers k(j) and periodicity n in the potential functions Ej = 1/2k(j)[1 + cos(n + j)] + E0 governing torsion around retinal's 13==14 bond

	RESULTS

Top Abstract Introduction Methods Results Discussion Conclusions References

Quantum chemical calculations

The adiabatic, singlet state potential surfaces for torsion around the double bond adjacent to the C==N bond of the retinal analog [CH₂---(CH)₃---(C₂H₃)---(CH)₂---NH---CH₃]⁺ are presented in Fig. 2. For the all-trans [13-cis] geometry the calculations predict, besides the ground state S₀(trans) [S₀(13-cis)], two closely spaced excited states, a state S₁(trans) [S₁(13-cis)], with predominantly single-excited electron configurations, below a state S₂(trans) [S₂(13-cis)], with predominantly double-excited electron configurations. This level ordering is in qualitative agreement with the configuration interaction calculations on a related polyenylic ion reported by Schulten et al. (1980). Both the S₁(trans) and S₂(trans) states are optically allowed from the ground state S₀(trans); the respective transition dipole moments are D_S0S1 = 3.4 a.u., D_S0S2 = 1.3 a.u., a relevant transition dipole moment from state S₁ is D_S1S2 = 0.9 a.u. The doubly excited character of the S₂ state suggests that this state is related to the two-photon allowed ¹A_g state in regular polyenes (Schulten and Karplus, 1972; Tavan and Schulten, 1986). The very close values of S₀  S₁ and S₀ S₂ excitation energies and comparable magnitudes of D_S0S1 and D_S0S2 indicate clearly the importance of both S₁ and S₂ excited states for retinal photodynamics.

As can be seen from Fig. 2, the potential surfaces along the complete [0°, 180°] interval of the torsional angle  exhibit three avoided crossings, two near ₁ = 30°, '₁ = 150°, and one near ₂ = 90°. We examined the surfaces near ₁ and ₂ by means of state-averaged calculations (see Methods) involving only the two states involved in the "crossings," i.e., states S₁ and S₂ near ₁ and states S₀ and S₁ near ₂. These calculations extend results of Schulten et al. (1995), providing an improved description of the potential surfaces close to the crossing regions. In case of the region near ₁, the results revealed that the interaction is of the "avoided crossing" type; the potential curves obtained in state-averaged calculations for two states do not differ considerably from those presented in Fig. 2, implying that the coupling between the S₁ and S₂ states in the ₁ region is strong.

In contrast, examination of the potential surface in the region near ₂, employing averaging over only the two states S₀ and S₁, leads to results that differ considerably from those obtained when all three states are involved in the average. This is demonstrated in the inset in Fig. 2, which shows that the separation between the surfaces becomes less than 1 kcal/mol in this case. Most likely, the crossing between S₀ and S₁ is of the "conical intersection" type (Cederbaum et al., 1981). Clarification of this matter requires consideration of further retinal degrees of freedom in the crossing region, e.g., of the stretch vibration of the C₁₃==C₁₄ bond (Cederbaum et al., 1981). Calculations performed for a smaller retinal analog indicate, indeed, the existence of a conical intersection between S₀ and S₁ (Bonacic-Koutecky et al., 1987).

Following the procedure outlined in Methods, the surfaces shown in Fig. 2 can be interpreted as arising from three nonadiabatic surfaces shown in Fig. 3: one surface (E_b) connecting the all-trans ground state S₀(trans) to the excited state S₂(13-cis); a second (E_a) connecting, conversely, the excited state S₂(trans) to the ground state S₀(13-cis); and a third one (E_c,bR) connecting the excited state S₁(trans) to the excited state S₁(13-cis).

The identification of the two nonadiabatic surfaces E_a and E_b constitutes a key assumption on which the further investigation is based. The surface E_b is consistent with the simple interpretation that the state S₀(trans) contains two -electrons in the bonding orbital of the double bond under consideration, and that a 180° torsion adiabatically lifts these -electrons into the antibonding orbital, i.e., into a doubly excited electron configuration; the occurrence of such configurations is characteristic of the S₂(13-cis) state. The surface E_a can be interpreted in a similar manner.

The distinctive features of the surfaces in Figs. 2 and 3 is the existence of regions of nonadiabatic interactions inducing, near ₁ ('₁), the crossing E_a  E_c,bR (E_b  E_c,bR) and, near ₂, the crossing E_a  E_b.

DME simulations

As outlined in Methods, the nonadiabatic potential surfaces in Fig. 3 served as a basis for a combined quantum/classical mechanical simulation, resulting in trajectories that provide coordinates, e.g., the torsional angle of the C₁₃==C₁₄ bond (t), as well as selected observables, e.g., the energy E_13-14 defined in Eq. 6. Fig. 4 presents the [E_13-14(t), (t)] diagram of a typical trajectory resulting from a 1-ps simulation, which ends with 13-cis retinal. Retinal is excited initially in the all-trans geometry by a 568-nm photon, i.e., it is placed on the E_c,bR surface at  = 0 (cf. Eq. 4). After a time ₁ the system reaches the first crossing region and crosses to the surface E_a, where it quickly descends to the second crossing point ₂ at time ₂. Here, after a brief period of meandering (~100 fs) near the crossing point, the system chooses to continue along E_a to complete the photoisomerization and reaches the 13-cis geometry. The energy E_13-14(t), initially follows the potential surfaces closely, but at the end of the photoreaction deviates significantly from the surface E₁[(t)], reflecting the fact that the system in the DME description used here is trapped in a mixed quantum state. For trajectories leading to all-trans retinal, a behavior similar to that shown in Fig. 4 is seen before the second crossing region is reached; however, after the meandering near the crossing point, the system chooses to switch to surface E_b, which leads back toward the all-trans isomer of retinal.

View larger version (21K): [in this window] [in a new window]   FIGURE 4   Energy E13-14 of the 13-14 torsional degree of freedom as a function of the torsional bond angle , i.e., the [E13-14(t), (t)] trajectory, for a sample molecular dynamics simulation. The energy E13-14 is compared to the three nonadiabatic energy surfaces Ea, Eb, and Ec,bR as shown in Fig. 3. Initially retinal is promoted in the all-trans geometry from Eb to Ec,bR; subsequently, retinal crosses to Ea near  = 1. At  = 2 retinal, in the case shown, chooses to continue its motion along Ea to reach the 13-cis isomeric state. Alternatively, retinal may have crossed to Eb to reform the all-trans state (not shown).

The diagonal elements _jj, j = a, b, c, of the density matrix describe the occupancies of the three electronic states corresponding to the surfaces E_a, E_b, E_c,bR. These occupancies are shown in Fig. 5. After excitation, i.e., at t₀ = 0, the occupancies of states corresponding to E_a and E_c,bR oscillate rapidly for ~200 fs because of the energetic proximity of the surfaces. The system is then seen to switch rapidly to the surface E_a, whereupon a fast torsion moves the system in a time (see Table 3) to the second crossing point with the surface E_b. Rapid oscillations are again discernible at this stage, mixing the states corresponding to E_a and E_b; because of the rather weak nonadiabatic coupling V_ab = 0.5 kcal/mol, 70% of the trajectories choose to stay on the E_a surface and proceed toward the 13-cis isomeric state. However, beyond the crossing point the states corresponding to E_b, E_c,bR remain populated to a significant degree. Because of environmental relaxation effects and because of quantum effects on the nuclear motion, not accounted for in our description, the occupancies of these states should vanish within a few femtoseconds, as demonstrated recently for a hydrated electron for which relaxation effects are expected to be stronger, however, than in the present case (Schwartz et al., 1996). The neglect of this decay is a shortcoming of the DME description adopted here, which, as a result, can lead to unphysical "average" final states, when in reality the system would choose between one of the possible "pure" final states with a certain probability distribution. In this respect one should remark that the quantum mechanical nature of the nuclear motion is significant even for molecular systems like retinal and protein, e.g., for nuclear motion with large effective masses (Bornemann et al., 1996). The simulations reported by Ben-Nun et al. (1998) and subsequent fully quantum mechanical simulations address this problem.

View larger version (33K): [in this window] [in a new window]   FIGURE 5   Occupancies jj(t) of the three nonadiabatic potential energy surfaces Ea, Eb, and Ec,bR as a function of time for the [E13-14(t), (t)] trajectory shown in Fig. 4. In the inset are shown the occupancies averaged over all simulations that resulted in 13-cis retinal.

A total of 100 simulations of the photoisomerization of wild-type bR were carried out, resulting in trajectories like the one shown in Figs. 4 and 5. Fifty simulations were also completed for the D85N and D212N mutants of bR, by replacing the potential surface E_c,bR with the surfaces E_c,D85N and E_c,D212N shown in Fig. 3 and given in Table 2. The simulations for the mutants required simulation periods of 1-12 ps because of the longer time required for retinal in D85N and D212N bR to cross to the E_a potential surface. The simulations exhibit asymptotically over 50% occupancy of either of the states corresponding to E_a or E_b. In the example shown in Figs. 4 and 5, this occupancy is almost 90%. The results of the simulations are summarized in Table 3: for each set of simulations, the fraction of trials that completed the all-trans  13-cis isomerization were identified, and the average crossing times ₁ and ₂ were determined.

In the 100 simulations of wild-type bR, 71 13-cis products emerged. Trials that failed to form a 13-cis isomer either did not move beyond the first crossing point ₁ or switched to the surface E_b at the second crossing point ₂; i.e., such trials reformed all-trans retinal. The simulated photoisomerization quantum yield measured 0.71 for wild-type bR, matching closely the observed quantum yield of 0.64 ± 0.04 (Schneider et al., 1989; Govindjee et al., 1990).

Time-resolved spectroscopy of bR has revealed two fast processes, one occurring within 200 ± 70 fs (Dobler et al., 1988) and the other occurring within 500 ± 100 fs (Mathies et al., 1988; Dobler et al., 1988; Du and Fleming, 1993). The observations were interpreted to suggest that photoexcitation is followed by rapid progression along the reaction coordinate from all-trans to 13-cis during the initial 200 fs. More recently, however, femtosecond spectroscopy covering an extremely wide spectral range has indicated a lack of a spectral shift of induced absorption, indicating a flat excited state surface. The observations also revealed a relaxation process with a decay time of 370 fs (Hasson et al., 1996): the E_c,bR  E_a crossing at ₁ = 332 fs can explain the observed fast decay.

Figs. 4 and 5 show the behavior of a single, albeit typical, quantum/classical mechanical trajectory. Of interest is the average behavior of bR. The inset in Fig. 5 shows the time evolution of the density of the states corresponding to the surfaces E_a, E_b, and E_c,bR averaged over all 71 simulations that formed 13-cis photoproducts. The time at which the occupations of E_c,bR and E_a intersect is ₁, whereas ₂ is near the time where the occupation of E₁ exceeds 50%. One can discern that retinal photoisomerizes within ~500 fs to the 13-cis geometry with a final population of the E₁ surface, on average, of ~80%.

For the D85N and D212N mutants of bR, the dynamics of retinal after light excitation is significantly slower in the initial phase. For D85N bR, with 30 of 50 simulations forming 13-cis isomers, the average time to the first crossing point is 5.59 ps: in the case of D212N bR with 37 of 50 completing the 13-cis isomerization, the first crossing occurred at 3.24 ps. The increase of ₁ is due to the feature of our model potentials (see Table 2 and Fig. 3) that the surfaces E_c,D8N and E_c,D212N cross the E_a surfaces at higher (relative to the E_c,x(trans) minima, x = bR, D85N, and D212N) energies than in case of the E_c,bR  E_a crossing. Interestingly, the time between the first and second crossings, , is approximately the same for native bR and for its D85N or D212N mutants; the mutants spent the major part of the overall reaction time "in front of" the first crossing point.

In previous molecular dynamics simulations of the photoisomerization of retinal in bR (Humphrey et al., 1995; Schulten et al., 1995; Humphrey et al., 1997), the isomerization products assumed two distinct structures: case I structures, in which the Schiff base N---H⁺ bond is oriented toward the cytoplasmic side of the protein, and case II structures, in which twists around retinal's single bonds oriented the N---H⁺ bond to point parallel to the plane of the membrane such that N---H⁺ remains connected to Asp⁸⁵ via hydrogen bonds with an intermediate water molecule. The latter product was identified with the K₅₉₀ intermediate, i.e., with the intermediate that actually initiates the proton pump cycle (Humphrey et al., 1995; Schulten et al., 1995; Xu et al., 1995). Fig. 6 compares the retinal geometries and surrounding binding sites for case I (Fig. 6 b), case II (Fig. 6 c), and all-trans retinal (Fig. 6 a). The combined quantum/classical simulations of the present study led to 66% case I products and 34% case II products. This distribution is similar to that in the earlier simulations (58% and 36%, respectively).

View larger version (23K): [in this window] [in a new window]   FIGURE 6   (a) Stereo view of retinal and its binding site in bR, including residues Asp85, Asp212, and internal water molecules. Waters A, C, and D are labeled for better comparison of reactant and product states. (b) Case I 13-cis photoproduct. (c) Case II 13-cis photoproduct. This figure was created with VMD (Humphrey et al., 1996).

	DISCUSSION

Top Abstract Introduction Methods Results Discussion Conclusions References

The three-state model for the photodynamics of retinal in bR has been described in the framework of the DME approach. Our simulations suggest that after light absorption into the first excited state, retinal carries out a rapid torsion around its C₁₃==C₁₄ bond, and in doing so crosses twice between electronic energy surfaces: the first crossing involves a relatively strong nonadiabatic coupling between the first and second excited states at a torsion of ~30°, 332 fs after light absorption; the second crossing occurs between the first excited state and the ground state, in which retinal moves essentially along a steep nonadiabatic surface, completing a 60° motion in another 130 fs. A choice of nonadiabatic couplings on the order of 1 kcal/mol can explain the overall quantum yield. Because of interaction with the protein, two isomerized reaction products arise, one (Case II) with the Schiff base proton in hydrogen bond contact with Asp⁸⁵, the other (Case I) with the Schiff base proton pointing to Asp⁹⁶ on the cytoplasmic side. This approach supports, therefore, the conclusions reached previously regarding a possible mechanism of bR's pump cycle (Schulten et al., 1995; Humphrey et al., 1995, 1997), namely that case II products may actually trigger vectorial proton transfer in bR's pump cycle. In one run, before a case II product finally formed, the system stayed in case I configuration for several hundred femtoseconds. This implies the possibility that case I and case II products interconvert. The fact that this conversion has not been observed in other runs suggests that the interconversion may happen only on a time scale longer than that covered by the present molecular dynamics simulations.

After photoexcitation from the ground-state surface E_b to the excited-state surface E_c,bR, retinal rotates around its C₁₃==C₁₄ bond until it reaches the first crossing point ₁. In this region, interaction with the nearly degenerate E_a surface allows the system to overcome a slight energy barrier before crossing to E_a. For the wild type, this crossing point was reached after 332 fs, which compares well with time-resolved spectroscopy measurements of this process by Anfinrud and co-workers, if one associates the observed 370-fs fast relaxation process (Hasson et al., 1996; Gai et al., 1998) with the crossing event. The near-degeneracy of the surfaces E_a() and E_c,bR() for small angles  implies a continuous interaction that prepares retinal for the crossing and induces a crossing to the second excited-state surface, even for a relatively weak, nonadiabatic coupling; thus the interaction of the system in the region of the first crossing region does not affect the quantum yield. In almost all simulations, retinal eventually surmounted the slight all-trans  ₁ energy barrier of 1-3 kcal/mol. The important characteristic of the initial phase of the excited-state dynamics is that retinal remains very close to its all-trans geometry; rapid torsion sets in only after the first crossing is completed.

Upon reaching the second crossing point ₂, near 90° torsion, retinal makes a nonadiabatic transition to the ground state, characterized by very weak coupling between the respective nonadiabatic surfaces. Essentially, retinal moves along the nonadiabatic surface E_a() (see Fig. 3), halting briefly at the crossing point, as shown in Fig. 4. The brief halt at the crossing point makes it unlikely that this stage of the photodynamics can be identified with the J₆₂₅ intermediate. The crossing determines, however, the quantum yield of formation of 13-cis-products; the observed value of 0.71 (Schneider et al., 1989; Govindjee et al., 1990) is reproduced only for a weak coupling V_ab = 0.5 kcal/mol. The second crossing point is reached within 462 fs for wild-type bR.

The D85N and D212N mutants of bR exhibit strongly increased fluorescence lifetimes without changes in the quantum yield (Song et al., 1993), a behavior that is explained well by the three-state model through an increase in the barrier for the first crossing event. The schematic potential surfaces employed in the simulations cannot be expected to reproduce exact crossing times. Nevertheless, our results provide an explanation of the behavior of bR mutants. Assuming that the crossing times are governed by a Boltzmann distribution at the crossing point ₁ and that the crossing points are at the same absolute height, the ratio of the mutants' excited-state lifetimes to that of bR is

(8)

where E_3,mutant = E_3,mutant(₁)  E_3,mutant(trans). Accordingly, mutations that shift the retinal absorption to the red increase the barrier to photoisomerization and lead to longer excited-state lifetimes. Inserting for E_3,mutant in Eq. 8 the values corresponding to the spectral shifts experienced by the D85N and D212N mutants predicts lifetime increases by factors of 20 and 1.8, respectively. The quantum yield for mutants such as D85N and D212N has been found experimentally to be similar to the wild-type value (Logunov et al., 1996a, b), and the yield obtained for these mutants in our simulations was likewise similar (0.60 and 0.74 for D85N and D212N, respectively, compared to 0.71 for the wild type). The three-state model has been essential for explaining both the quantum yield and the formation time of J/K intermediates.

A similar scheme taking into account state average contributions from different potentials and nonadiabatic couplings has been adopted by Tallent et al. (1992) to simulate the retinal photoisomerization in rhodopsin. The present work differs in the following aspects. Results of the current quantum chemical calculations show that there is a barrier near the 30° and 150° points on S₁, whereas in the study by Tallent et al. (1992) there was no barrier in the first excited state. Present simulations also involved an all-atomic model of the whole structure of bR, accounting thereby for environmental effects.

A deficit of the present study is the limited knowledge of the nonadiabatic couplings V_jk. Unfortunately, one cannot calculate these couplings (in particular, their dihedral angle dependence) in the framework of available ab initio programs. This study estimated the nonadiabatic couplings from the energy splitting of the adiabatic surfaces in the regions of crossings. A second deficit is the neglect of bond length changes in the excited-state nuclear dynamics; such changes are likely to occur on a 10-fs timescale and can affect the S₁-S₂ energy splitting.

The DME description adopted in this paper for the dynamics of classical systems moving on multiple electronic energy surfaces is only approximate. This is apparent from Eq. 5, which relates the force _13-14,k() to a coherent average over all surfaces. The single crossing behavior, as governed by the Landau-Zener formula (Zener, 1932; Landau and Lifshitz, 1977) and its generalization (Crothers, 1981; McDowell and Brandsen, 1992; Kayanuma, 1993), is reproduced well by the DME description (unpublished results), but this approach does not account properly for the "collapse" of the wave function onto wave packets eventually moving incoherently on the individual surfaces. This deficiency is most glaring in Fig. 4, which shows that E_13-14 is significantly above the ground-state energy surface E₁() because of an admixture of state 1 and state 2 character in the asymptotic wave function. As a result, the effective potential surface behind the crossing at ₂ is too shallow. The error introduced by the DME description is tolerable as long as occupancy of a single state is dominant. This is not guaranteed, however, and an extension of the present treatment to an essentially exact description of multiple curve crossing dynamics as provided in Ben-Nun et al. (1998) is desirable.

	CONCLUSIONS

Top Abstract Introduction Methods Results Discussion Conclusions References

Quantum chemical calculations on a retinal analog suggest the involvement of the two lowest excited states and the ground state of the protonated Schiff base of retinal in the primary all-trans  13-cis photoisomerization in bacteriorhodopsin. The calculations suggested a strong nonadiabatic coupling between the two excited states near the all-trans and 13-cis geometries; a weak coupling suggesting a conical intersection was identified between the first excited state and the ground state at a 90° torsion. Combined quantum/classical simulations utilizing three-state potential surfaces in the nonadiabatic representation, which correspond closely to the calculated adiabatic surfaces, reproduce successfully observed crossing times and quantum yield.

The three-state model described here differs from previous models based on ultrafast spectroscopy (Mathies et al., 1988; Dobler et al., 1988; Du and Fleming, 1993). In the new model, photoexcitation of retinal is followed not by rapid movement along the reaction coordinate due to excitation into a repulsive Frank-Condon region, but by motion along a relatively flat excited-state surface and crossing to a region with a second excited state, which exhibits a steep potential gradient inducing rapid isomerization. A proper description of this behavior requires the inclusion of at least the first two excited states. The observed spectral and fluorescent changes occurring within 1 ps after photoexcitation of bR arise because of nonadiabatic transitions between the excited states and an excited state and the ground state. We suggest that a significant fraction of the excited-state lifetime is spent surmounting a small energy barrier before reaching the first nonadiabatic crossing region, and that the fluorescence lifetime of bR is controlled by this barrier and can be altered strongly through mutations that affect the relative stability of retinal's two lowest excited states.