This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Hayashi, S. Articles by Schulten, K. Search for Related Content PubMed PubMed Citation Articles by Hayashi, S. Articles by Schulten, K.

Molecular Dynamics Simulation of Bacteriorhodopsin's Photoisomerization Using Ab Initio Forces for the Excited Chromophore

Beckman Institute, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801

Correspondence: Address reprint requests to Klaus Schulten, E-mail: kschulte{at}ks.uiuc.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
Retinal proteins are photoreceptors found in many living organisms. They possess a common chromophore, retinal, that upon absorption of light isomerizes and thereby triggers biological functions ranging from light energy conversion to phototaxis and vision. The photoisomerization of retinal is extremely fast, highly selective inside the protein matrix, and characterized through optimal sensitivity to incoming light. This article describes the first report of an ab initio quantum mechanical description of the in situ isomerization dynamics of retinal in bacteriorhodopsin, a microbial retinal protein that functions as a light-driven proton pump. The description combines ab initio multi-electronic state molecular dynamics of a truncated retinal chromophore model (N-methyl--methylpenta-2,4-dieniminium cation fragment) with molecular mechanics of the protein motion and unveils in complete detail the photoisomerization process. The results illustrate the essential role of the protein for the characteristic kinetics and high selectivity of the photoisomerization: the protein arrests inhomogeneous photoisomerization paths and funnels them into a single path that initiates the functional process. Supported by comparison with dynamic spectral modulations observed in femtosecond spectroscopy, the results identify the principal molecular motion during photoisomerization.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
Retinal proteins, often called rhodopsins, are involved in a variety of responses of living cells to light, such as conversion of light to chemical energy and phototaxis in bacteria, or vision in higher organisms. All rhodopsins contain the molecule retinal as their chromophore. The activation of rhodopsins is initiated by a photoinduced isomerization of retinal around a specific double bond, which represents one of the fastest chemical reactions in nature (Birge, 1990). The extremely fast (200–500 fs) process prevents undesired nonspecific relaxation of the chromophore from the electronically excited state and, hence, furnishes pigments of high sensitivity. Bacteriorhodopsin (bR), shown in Fig. 1, functions as a light-driven proton pump in the purple membrane of Halobacterium salinarum (Oesterhelt and Stoechenius, 1971) and is the most investigated member of the family. Photoisomerization of the chromophore from the all-trans to the 13-cis configuration triggers a photocycle during which one proton is pumped from the cytoplasmic side to the extracellular side of the membrane. Numerous femtosecond spectroscopic experiments (Mathies et al., 1988; Gai et al., 1998; Ruhman et al., 2002) have been performed to identify the kinetics of the photoisomerization, which is schematically illustrated in Fig. 1. Time-resolved vibrational spectroscopy has suggested that bond stretching motions are the primary conformational changes of retinal in the excited state (Song and El-Sayed, 1998). Recent progress in such measurements has made it also possible to observe in real time wavepacket motions of molecular vibrations during the isomerization of retinal in bR (Ye et al., 1999; Kobayashi et al., 2001) and in the visual receptor rhodopsin (Wang et al., 1994) through dynamic spectral modulations of photoabsorption and emission spectra. These studies have revealed characteristic vibrational motions along with their time-dependent couplings in the transient processes and demand a detailed interpretation in terms of the related protein function.

View larger version (34K): [in this window] [in a new window]   FIGURE 1  (a) Structure of bR. Seven transmembrane helices represented by green cylinders form a channel that embeds a retinal chromophore highlighted in purple. (b) Chemical structure of the retinal chromophore and its photoisomerization reaction in bR identified experimentally (Mathies et al., 1988; Gai et al., 1998; Ruhman et al., 2002). The retinal chromophore has a polyene chain with six double bonds and is bound to Lys-216 via a Schiff base. In the present simulation, a short chromophore model consisting of the Schiff base half, i.e., three conjugated double bonds from C to C11 enclosed by a red line, is employed as the quantum mechanically treated molecule. In the resting state of bR (BR570), retinal assumes an all-trans configuration. The isomerization process starts with a photoexcitation at 570 nm to the Frank-Condon region. 100 fs after the photoexcitation, a well-characterized excited state, called I460, with absorption at 460 nm and emission in a longer wavelength region forms. The decay of I460 within the next 500 fs coincides with the isomerization process and results in a vibrationally hot 13-cis product in the ground state, J630. Vibrational cooling in 3 ps leads to K610, which is the first intermediate that can be trapped in cryogenic conditions.

To address theoretically the issues described above, three key methodological advances are required. First, since the photoisomerization process involves molecular motions in many degrees of freedom and proceeds on multi-electronic state energy surfaces, methodology is needed that is capable of describing the chromophore's electronic structure and potential surfaces without imposing assumptions. Recently, high level ab initio quantum mechanical (QM) calculations have become applicable to the chromophore models of increasing size (Garavelli et al., 1997, 1998; González-Luque et al., 2000; De Vico et al., 2002; Vreven and Morokuma, 2000; Molnar et al., 2002; Ben-Nun et al., 2002). The calculations have given reaction energy curves along minimum energy reaction paths (MEPs) or potential profiles around the conical intersection, revealing the complex nature of the potential energy surfaces of the photoisomerization process. Moreover, the description of reaction process requires not only accurate potential surfaces, but also the nonadiabatic electronic coupling between different electronic states that induces electronic transitions during the photoisomerization. Recent computational advance in this regard, amenable through sophisticated program packages, enables one to compute the nonadiabatic electronic coupling along all nuclear coordinates degrees of freedom.

Second, the dynamics of the reaction process should be explicitly described. Although ab initio QM analyses of the MEP energy curves have provided a basic insight into the reaction process, actual dynamic paths can deviate from the MEP due to thermal fluctuation and reaction excess energy. Especially, the characteristic incoherent kinetics and dynamic spectral modulations found in the experimental studies can be directly captured only through molecular dynamics (MD) simulations. In this respect, it is fortunate that ab initio QM MD simulations have finally been feasible to compute on the fly energies and forces of electronically excited molecules in isomerization processes (Vreven et al., 1997; Ben-Nun and Martínez, 1998; Doltsinis and Marx, 2002). These studies shed light on the complex photoisomerization dynamics that involves intramolecular vibrational energy transfers among multi-dimensional reaction coordinates.

Third, the photoisomerization simulation has to take into account the protein environment to elucidate the molecular mechanism of the protein environment catalytic effect essential for the function. The success of the ab initio QM MD simulation has been restricted so far to applications for small systems in isolated condition (in vacuo). The photoisomerization in protein environments (in situ) has been investigated only in the framework of less accurate semiclassical MD simulations using empirical potential energy functions obtained from prior QM calculations (Birge and Hubbard, 1980; Birge et al., 1989; Tallent et al., 1992; Ben-Nun et al., 1998; Humphrey et al., 1998). An alternative treatment applied a hybrid quantum mechanical/molecular mechanical (QM/MM) calculation (Warshel, 1976; Warshel et al., 1991; Warshel and Chu, 2001) that combined an MM description of the protein with a semiempirical (Pariser-Parr-Pople) QM calculation of retinal. Although the latter approach has successfully described multi-electronic processes in situ, the calculation is limited in its account of the excited-state electronic wave function and potential surface, neglecting to describe many degrees of freedom properly.

In this article, we present ab initio QM/MM MD simulations of the complete isomerization process inside bR, including the excited-state dynamics of the chromophore, the electronic transition from the excited state to the ground state, and the subsequent relaxation of the ground state. Combining an ab initio description of the Schiff base half of retinal, the N-methyl--methylpenta-2,4-dieniminium cation fragment depicted in Fig. 1, at the complete active space self-consistent field (CASSCF) level of theory with an MM description of the remaining part of the retinal chromophore and of the protein, enabled us to simulate for the first time the multi-electronic state dynamics of photoisomerization in a complex protein environment without assuming empirical potential functions for the core degrees of freedom.

Unfortunately, the very demanding computation of the CASSCF method for the QM part forced us to employ the short chromophore analog. However, this analog includes the essential retinal segment involved in the isomerization in bR, namely the three terminal double bonds, including the Schiff base. But it should be noted that the approximation of the short analog introduces several artifacts of the potential energy profiles and the electronic properties as investigated by ab initio QM studies (Garavelli et al., 1998; González-Luque et al., 2000). First, the shorter conjugation results in an overestimated excitation energy, leading to a steeper potential profile along the isomerizing bond. The altered photoisomerization potential surface also may severely affect vibrational behavior, in particular vibrational couplings between the stretching and torsional modes. Second, the three conjugated double bond system biases the isomerization around the center bond, i.e., C₁₃=C₁₄, which is actually the bond undergoing the isomerization in bR, whereas a longer chromophore analog in vacuo favors isomerization around either C₁₃=C₁₄ or C₁₁=C₁₂ (De Vico et al., 2002). Finally, the short polyene chain prevents delocalization of the positive charge in the ß-ionone ring half in the excited state, altering the chromophore-protein interaction.

Nevertheless, the present novel in situ simulation and its comparison with in vacuo simulation identify the essential role of the protein environment on dynamics and catalysis of the photoisomerization process in bR. The protein environment significantly affects the photoisomerization kinetics and dynamics. The use of the short analog permitted us to calculate 11 trajectories, and to elucidate the characteristic incoherent dynamics of isomerization in bR. The study also shows that the in situ photoisomerization attains high isomerization bond selectivity, as experimentally observed. Furthermore, since ab initio excited-state MD simulations can naturally provide time evolution of the emission spectrum, we were able to directly examine the dynamic spectral modulation phenomena observed in femtosecond spectroscopy (Ye et al., 1999; Kobayashi et al., 2001). Despite the limitations due to the short chromophore model, ab initio description of the vibrational motions successfully reproduced major features of the dynamic spectra observed experimentally, proposing molecular origins underlying the spectral phenomena. Although the approximation of the short chromophore model used in our investigation has to be reassessed in future studies, the results obtained reveal already now surprising aspects of the multi-dimensional and multi-electronic state dynamics of retinal's in situ photoisomerization.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
Photoisomerization trajectories were computed by ab initio QM/MM MD simulations. The QM and the MM regions were described at the ab initio molecular orbital and molecular mechanics levels of theory, respectively. The Schiff base half of the chromophore with three conjugated double bonds was chosen as the QM region (Fig. 1). Electronic structure, isomerization potential profile, and dynamics of a similar system of the Schiff base with three conjugated double bonds have been examined previously (Garavelli et al., 1997; Vreven et al., 1997). The rest of the chromophore, the protein, and water molecules were treated as the MM region. The electronic wave functions of the QM region in the ground (S₀) and the first -^* excited (S₁) states were determined by the state-average CASSCF method (Werner and Knowles, 1985; Knowles and Werner, 1985). All six valence and ^* orbitals with six electrons were included in the CAS. Dunning's double -basis functions (Dunning, 1971) were used. The MM region was described by the AMBER94 force field (Cornell et al., 1995) for the protein, by the TIP3P force field (Jorgensen et al., 1983) for water molecules, and through previously determined parameters for the MM part of the chromophore (Tajkhorshid et al., 2000; Hayashi et al., 2002). A conventional QM/MM scheme using the link atom approach (Field et al., 1990) was employed. Open valences of covalent bonds at the boundaries dividing the QM and MM regions were filled by dummy hydrogen atoms. Electrostatic interaction of electrons in the QM region with the MM effective point charges were taken into account through one-electron integrals. A QM/MM interface combined with our own MD routine was incorporated into the quantum chemistry package, MOLPRO (a package of ab initio programs written by H.-J. Werner and P. J. Knowles, with contributions from J. Almlöf and co-workers). The trajectories were computed by solving the Newtonian equations of motion with energies and forces obtained from the QM/MM calculations. We employed the velocity-Verlet integrator (Swope et al., 1982) with a 0.5 fs time step. Since energy changes rapidly around the electronic transition points, the time step was switched to 0.25 fs whenever the energy difference between the S₀ and S₁ states dropped below 0.01 a.u. The QM/MM MD simulation of 100 fs required 40 h on a Compaq Alpha21264 500 MHz processor.

A model of bR was constructed based on high resolution x-ray crystal structures (Luecke et al., 1999; Belrhali et al., 1999), as described elsewhere (Hayashi et al., 2002). To generate an ensemble of initial coordinates and velocities to be used in photoisomerization trajectories, a classical MD simulation of the resting state of bR (BR₅₇₀) at constant temperature (300 K) was carried out (Hayashi et al., 2002). A 12-Å residue-based cutoff was used for nonbonding interactions. The system was equilibrated for 500 ps with 1 fs time step while constraining bond distances that included hydrogen atoms, followed by a 50 ps equilibration with 0.5 fs time step without any bond constraints. The 11 initial configurations used for photoisomerization simulations were then sampled from the following 200 ps equilibrium MD simulation. The on-the-fly QM/MM MD trajectories on the S₁ excited state were started from the initial configurations simulated in an NVE ensemble. The total energy of the system was well conserved and its standard deviation was 0.3 kcal/mol. Each excited-state trajectory was continued for 600 fs or until the trajectory achieved 75% product (see below), resulting in simulation times of 90–600 fs. Along the excited-state trajectories, semiclassical nonadiabatic coupling elements between the S₀ and S₁ states, , were computed using

(1)

where and are adiabatic (Born-Oppenheimer) electronic wave functions, and and are nuclear coordinates and velocities, respectively. At electronic transition points between S₁ and S₀ determined by a semiclassical calculation (see below), the S₀ state trajectories were started to obtain the products. The total number of S₀ state trajectories spawned from the 11 trajectories of S₁ was 65. Excess energies upon the transitions from S₁ to S₀ were distributed to nuclear velocities along the nonadiabatic coupling vectors (Herman, 1984).

To evaluate the transition probabilities from S₁ to S₀ at the transition points, the quantum mechanical probabilities of S₁ and S₀, and , respectively, were computed by solving a semiclassical time-dependent Hamiltonian,

(2)

and are potential energies of S₁ and S₀ along the trajectory on the S₁ state, respectively, and is the nonadiabatic coupling element given in Eq. 1. Energy levels of S₁ and S₀ along the S₀ state trajectories separate rapidly after the electronic transitions (see below), indicating fast decoherence between the S₀ and S₁ nuclear wave functions and prompt decay in the S₀ state. To incorporate the decay effect in Eq. 2, an absorbing term in the S₀ state, , was introduced. was set to be the decoherence time constant, 3.8 fs, estimated based on a Gaussian wavepacket approximation at short time and high temperature limits (Schwartz et al., 1996). The differential equation (Eq. 2) was solved numerically by the fourth order Runge-Kutta method with a time step of fs. For the Hamiltonian matrix elements, spline-fitted functions of time were used to interpolate values from the trajectory calculations. The calculated population of the S₁ state exhibits stepwise decays to either 13-cis or all-trans products at the transition points. The branching ratio between the cis and trans products was then estimated by accumulating the stepwise transition probabilities.

We also performed ab initio MD simulation of photoisomerization of the retinal analog in vacuo to compare the dynamics in vacuo with that in bR (in situ). The retinal analog (Fig. 1) and the method used were the same as in the in situ simulation, except for absence of the QM/MM coupling. An ab initio MD study of the photoisomerization of a similar retinal analog with three conjugated double bonds (Vreven et al., 1997) has shown the time constant of photoisomerization in vacuo to be 70 fs. However, as described below, six trajectories in situ exhibit remarkably longer lifetimes of the excited state (>180 fs). We therefore computed the six trajectories in vacuo starting with the same initial coordinates and velocities as in situ to clarify the molecular origin of the difference. The resultant timescale of the photoisomerization in vacuo was consistent with the previous ab initio QM studies (Garavelli et al., 1997; Vreven et al., 1997).

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
Fig. 2 a shows the time evolution of potential energies of the first excited state (S₁) and the ground state (S₀) along a typical photoisomerization trajectory. Before photoexcitation from the resting BR₅₇₀ state to the Frank-Condon (FC) region at t = 0, the all-trans planar chromophore assumes a slight twist around the C₁₃=C₁₄ bond, as suggested in previous MD (Tajkhorshid et al., 2000) and ab initio QM/MM (Hayashi and Ohmine, 2000) simulations. Immediately after the photoexcitation, the in-plane C₁₃=C₁₄ bond relaxation takes place in 50 fs, as suggested experimentally (Song and El-Sayed, 1998; Kobayashi et al., 2001) and theoretically (Vreven et al., 1997; Garavelli et al., 1997, 1999; Warshel and Chu, 2001), but the chromophore retains its all-trans configuration until t = 130 fs, when the torsion around C₁₃=C₁₄ starts to evolve. Once the dihedral angle reaches a value of about 90°, the twisted conformation remains stable as long as the trajectory proceeds on the S₁ state surface. The overall energy slope in the S₁ state from the FC region toward the 90° twisted conformation is rather gentle, but the energy shows large fluctuations due to the chromophore's intramolecular vibrations and its interaction with the protein environment. On the other hand, energy of the (unpopulated) S₀ state increases steeply as the C₁₃=C₁₄ bond is twisted and approaches the S₁ state energy near the 90° twist, where accordingly frequent electronic transitions occur. The orbitals on carbon atoms C₁₃ and C₁₄ are almost perpendicular at these transition points, even though the carbon center C₁₄ assumes a highly pyramidalized structure (see below), and the nonadiabatic coupling element between S₁ and S₀ exhibits large peaks (Weiss and Warshel, 1979; Molnar et al., 2002), as seen in Fig. 2 a. A semiclassical analysis of the electronic transition process (see Methods) indicates that a significant fraction of the S₁ population decays to S₀ at each crossing point. After the transition to S₀, the highly twisted conformation relaxes quickly to a planar form, either to a 13-cis product (the J₆₃₀ state) or back to the all-trans one (the BR₅₇₀ state).

View larger version (36K): [in this window] [in a new window]   FIGURE 2  (a) Potential energies and populations of the S1 and S0 states along a typical photoisomerization trajectory. The top panel shows time evolution of the potential energies defined as the sum of the electronic energy of the QM region and interaction energy between the QM and the MM regions (see Fig. 1 and Methods). Solid and dashed purple lines indicate the potential energies of S1 and S0, respectively, along the S1 state trajectory. The energies along the trajectories on the S0 state after electronic transitions from the S1 state are drawn by solid lines in blue and orange for the 13-cis (J630) and all-trans (BR570) branches, respectively; the dashed lines in blue and orange indicate the energies of S1 along the S0 trajectories. The chromophore structures shown in the top panel are snapshots at the FC region (t = 0 fs), at a transition point (t = 190 fs), and of a 13-cis product (t = 240 fs). The middle panel shows the nonadiabatic coupling element between the S1 and the S0 states () along the S1 state trajectory. Peaks of indicate the transition points. The bottom panel displays the S1 population along the trajectory. The population decreases stepwise at the transition points, and five distinct transitions during the last 50 fs furnish 75% of the overall transition. S0 trajectories obtained from the first three transitions result in 13-cis products with 52% yield, and the next two transitions lead to all-trans products with 23% yield. (b) Energy differences between S1 and S0 from 11 photoisomerization trajectories and the time evolution of population of the products. In the top panel, the energy differences along the S1 state trajectories (emission energies) are drawn by solid lines in purple. The negative energy differences along the S0 state trajectories (negatives of absorption energies) after the transitions are indicated by solid lines in blue and orange for the 13-cis and all-trans branches, respectively. The dashed lines connect the S1 and S0 states at the transition points. In the bottom panel, populations of the 13-cis and all-trans products resulting from the 11 trajectories are drawn in blue and orange, respectively. A dashed line in gray indicates the residual population of S1, which did not yet decay to S0.

Despite the presence of three isomerizable double bonds in the analog used in our calculations, isomerization proceeds in all trajectories only at the C₁₃=C₁₄ bond, in agreement with experimental evidence (Birge, 1990). Moreover, the twist around this bond occurs in only one rotational sense, namely the one for which the Schiff base hydrogen approaches Asp-212. The observed unidirectional isomerization is consistent with the structures of the early intermediate K₆₁₀ (Fig. 1) recently predicted theoretically (Hayashi et al., 2002) and determined by x-ray crystallography (Schobert et al., 2002; Matsui et al., 2002), where the Schiff base N–H bond points indeed toward Asp-212.

The bond selectivity and unidirectionality of the isomerization in bR are in striking contrast with those of the isolated chromophore (in vacuo). Fig. 3 shows the time evolution of emission energies along in vacuo trajectories and snapshots of the retinal analog with 90° twisted dihedral angles around the isomerizing bonds. Repeating the simulations for six trajectories in the absence of the protein environment, but using the same initial coordinates and velocities as in the in situ simulations, leads to two isomerizations around the C₁₅=N bond in either rotational sense. The results clearly demonstrate the role of the protein environment in the high bond selectivity and directionality of the photoisomerization in bR. At room temperature the retinal chromophore undergoes large thermal fluctuations so that the initial conformation and motion at the moment of photoexcitation can lead to isomerization around the C₁₅=N bond, and to a twist in any direction. The protein environment blocks alternative paths and funnels the trajectory to the isomerization that is coupled to the biological function of the system, i.e., to the C₁₃=C₁₄ isomerization in the case of bR.

View larger version (31K): [in this window] [in a new window]   FIGURE 3  (a) Time evolution of emission energies along six trajectories in vacuo (top), and of the product formation (bottom). Solid lines in red and blue in the top panel indicate emission energies along isomerizations around C13=C14 and C15=N, respectively. A blue dashed line shows the emission energy along an isomerization around C15=N in the direction opposite to that seen in the other trajectories. The product formation drawn by a red line in the bottom panel was computed with the semiclassical method (see Methods) for four trajectories that resulted in the isomerization around C13=C14. A line in gray indicates the residual population of S1 that did not yet decay to S0. (b–d) Snapshots of molecular structures of the retinal analog with 90° twisted dihedral angles around the isomerizing bonds. (b) Isomerization around C13=C14, (c) around C15=N, and (d) around C15=N in the opposite direction.

On the other side, a previous MD simulation (Tajkhorshid et al., 2000) and an ab initio QM/MM calculation using the full retinal chromophore model as the QM region (Hayashi and Ohmine, 2000) have shown that the chromophore in resting bR assumes a significant twist in the C₁₃=C₁₄ bond toward the Asp-212 direction, which is indicative of protein forces favoring isomerization around C₁₃=C₁₄. As suggested in a previous study (Hayashi et al., 2002), the unidirectionality of the isomerization is also due to the protein forces. The chromophore undergoes the unidirectional torsion around C₁₃=C₁₄ to accommodate itself in the binding pocket. Key interactions of the chromophore with Trp-86 and a proximate water molecule that induce this torsion have been identified (Tajkhorshid et al., 2000; Hayashi and Ohmine, 2000).

It is noteworthy that the lifetime of I₄₆₀ ranges from 50 fs to more than 500 fs for the in situ trajectories; the isomerization in some trajectories takes place immediately after formation of I₄₆₀, whereas some trajectories stay much longer in this state. The wide variation of the lifetime of I₄₆₀ in bR is another feature that is absent from the faster and more coherent isomerization dynamics of the retinal Schiff base in vacuo; the isomerization in vacuo completes within 50–100 fs as shown by Vreven et al. (1997) and confirmed by our simulations as shown in Fig. 3. In bR, some trajectories stay and wander in the I₄₆₀ state for a significant amount of time to find an escape path toward the 13-cis isomer. The escape from I₄₆₀ requires matching of the rotational motion around the C₁₃=C₁₄ bond with protein thermal motions as well as other torsions of retinal to minimize the chromophore's conformational change necessary for isomerization in a narrow binding pocket. Key protein motions coupled to the isomerization involve Lys-216, which binds the chromophore. The incoherent behavior of the escape from I₄₆₀ toward the product formation in bR can explain femtosecond spectroscopic observations, such as silent vibrational bands of pump-probe signals (Kobayashi et al., 2001) and the statistical characteristics of stimulated emission decay (Ruhman et al., 2002) after the formation of I₄₆₀.

Fig. 2 b also shows the time evolution of the populations of 13-cis and all-trans products. The products form with a time constant of 200 fs, which is in good agreement with the experimental values of 240–500 fs (Gai et al., 1998); the deviation is most likely due to a narrow plateau of the potential energy curve along the C₁₃=C₁₄ dihedral angle around the FC region and the steeper potential profile toward the 90° twisted region that arise in the short retinal analog used in the present simulations, but not in the complete chromophore (Garavelli et al., 1998; González-Luque et al., 2000). The 13-cis product formation is dominant over return to all-trans retinal, as a result of the first transitions that favor exclusively the cis form (see below). The cis-trans branching ratio calculated from the 11 trajectories is 0.54:0.46, as seen in Fig. 2 b, which agrees well with the experimental value, 0.6:0.4 (Tittor and Oesterhelt, 1990). The slightly underestimated quantum yield of the cis product is likely due to the altered potential profile caused by the use of the short chromophore model and lack of dynamic electronic correlation in the present treatment; there also remains ambiguity in the evaluation of the branching ratio due to poor statistical sampling in the long time region, a limit imposed by the extreme demand on computation time.

At a 90° twist of the C₁₃=C₁₄ bond, each trajectory encounters many transition points (between 2 and 10) until achieving 75% product formation. After the electronic transition, the energy levels of S₀ and S₁, calculated along the S₀ state trajectory, separate rapidly, and recrossing to S₁ does not occur. The first transition in each trajectory furnishes exclusively the cis product, because the transition process involves the large momentum of the rotation around C₁₃=C₁₄ from the trans conformation toward the cis one created by energy funneling from I₄₆₀ to the transition region. In the calculated trajectories, the first transition event contributes on the average 17% cis product, i.e., nearly one-third of all successful isomerization events. Therefore the exclusive transition to the cis conformation at the first transition point enhances the quantum yield of the cis product, as already proposed by Warshel and co-workers (Weiss and Warshel, 1979; Warshel and Chu, 2001). On the other hand, branching to the cis and trans products at transition points after the first one seems to be rather random, and a regular pattern such as alternate transitions to cis and trans as suggested by Warshel and co-workers (Weiss and Warshel, 1979; Warshel and Chu, 2001) was only observed in four trajectories. This random behavior is due to chaotic motions in the transition region originating from highly excited hydrogen-out-of-plane (HOOP) and bending vibrations around the C₁₄ atom, as discussed below.

The electronic transition events were found to be strongly perturbed by the electrostatic field in the protein. The multiple transition in bR is different from transition in vacuo; in the latter case, most of the product is yielded through only one or two transition events (Vreven et al., 1997) (Fig. 3). The relatively small transition probability at each transition point in bR originates from strong interaction of the protonated Schiff base with the counterion groups (Asp-85 and Asp-212) and a proximate water molecule. Fig. 4 displays key molecular orbitals at a transition point, indicating the so-called sudden polarization (Salem, 1979; Bonai-Koutecky et al., 1984) of the electronic structure, i.e., a complete localization of the positive charge in the C₁₄N and C₁₁C₁₃ halves in S₀ and S₁, respectively. Hence, the counterion and water molecules located in the vicinity of the Schiff base stabilize the S₀ electronic state more than S₁, which leads to a shift of the energy surface crossing seam and a larger separation of the energy levels at the transition point, and decreases the transition probability.

View larger version (33K): [in this window] [in a new window]   FIGURE 4  (a) Two orbitals (L and R) characterizing the S0 and S1 electronic structures. In the S0 state, each of the orbitals is occupied by one electron, whereas in the S1 state, L is doubly occupied and R is vacant. Orange dotted lines indicate the isomerizing C13=C14 bond. The localized nature of the orbitals gives rise to localization of the chromophore's positive charge in the C14N (L) and C11C13 (R) halves for S0 and S1, respectively (sudden polarization, see text). (b) Typical conformations of the chromophore for efficient transitions at the electronic transition points in situ. The top and bottom panels show views from side and top, respectively. In the top panel, a high degree of pyramidalization at C14 can be seen, and the wagging angle of C14–H14 (124°) largely deviates from planarity (180°). The doubly occupied L orbital at C14 exhibits an sp3-type hybridization (Salem, 1979) consistent with the pyramidalization frequently found at crossing points. The bottom panel reveals a straightening at C14. The bond angles of C13=C14–C15 (135°) and C13=C14–H14 (98°) largely deviate from the sp2 angle (120°). This is due to a strong hyperconjugation interaction of the vacant R orbital with the C14–H14 orbital providing an sp-type mixing (Ohmine, 1985).

The very demanding computation for this study did not allow us to simulate the complete relaxation of bR to the K₆₁₀ intermediate after retinal's isomerization. However, the observed photoisomerization dynamics corresponds closely to that obtained in a previous study using classical MD simulations (Hayashi et al., 2002), which examined the formation of K₆₁₀. Hence, we expect that the photoisomerization event in the present simulation leads also to the same K₆₁₀ intermediate. In this regard it is noteworthy that Hayashi et al. (2002) included an ab initio QM/MM analysis of K₆₁₀ for the complete retinal chromophore and several groups in the vicinity of the Schiff base, which reproduced the experimentally measured (Birge and Cooper, 1983; Birge et al., 1989) energy storage of 16 kcal/mol in K₆₁₀, and identified the molecular mechanisms of energy storage, namely conformational distortion of retinal and weakened interaction between the Schiff base and the surrounding polar groups. Similar mechanisms were proposed by other researchers (Birge and Cooper, 1983) to be responsible for storing energy after retinal's isomerization in bR.

A stimulated emission pumping measurement by Ruhman et al. (2002) has proved that I₄₆₀ is an excited-state intermediate of the photocycle that is arrested in a specific configuration and decays kinetically to a next photocycle intermediate. The molecular structure and the kinetic decay of I₄₆₀ are discussed above. On the other hand, a key observation of a dynamic spectral modulation in the stimulated emission signal by Ye et al. (1999) detected dynamics within the I₄₆₀ state. Fig. 5 displays the time evolution of emission energy in a trajectory with a long-lived I₄₆₀ state. The spectrum exhibits high frequency oscillations, superimposed on a slow dynamic modulation with a 200 fs period due to molecular vibrations. Indeed, the slow 200 fs component reproduces well a prominent feature in the femtosecond-stimulated emission spectrum observed by Ye et al. (1999). This slow oscillation is assigned to torsional rocking around C₁₃=C₁₄, during which the twist develops up to 40° away from planarity, but a mismatch with neighboring torsion and protein motion prevents the trajectory from proceeding to the 90° twisted conformation.

View larger version (32K): [in this window] [in a new window]   FIGURE 5  (a) Emission energy along a trajectory with a long-lived I460 that does not undergo isomerization in 600 fs. The emission energy (red line) exhibits oscillations due to molecular vibrations. The blue line indicates a slow oscillatory component with a 200 fs period obtained by smoothing with a time window function. A Hunning function with a halfwidth at half-maximum of 64 fs was used for the time window. (b) Spectrogram showing the time evolution of characteristic frequencies of the oscillations of the emission energy in (a). The spectrogram was calculated by Fourier transformation of the emission energy gated by the Hunning time window (Kobayashi et al., 2001). The time evolution of the frequencies was computed along the window gate delay time. The frequencies are uniformly scaled by a factor 0.9 (Garavelli et al., 1999). Purple lines indicate traces of the frequency peaks.

Kobayashi et al. (2001) have analyzed the fast oscillatory component of their 5 fs pump-probe absorption spectrum through a spectrogram and revealed remarkable dynamic frequency shifts and vibrational mode couplings. Fig. 5 shows the respective spectrogram calculated here for the long-lived I₄₆₀ trajectory. A strong peak at 1700 cm^-1 is due to the C₁₃=C₁₄ and C₁₅=N stretching vibrations excited during the formation of I₄₆₀. The signal decays quickly, and the peak frequency shifts downward in 100 fs, in good agreement with observation (Kobayashi et al., 2001). The downward shift is attributed to a decrease of the bond order of C₁₃=C₁₄ due to its torsion. After 300 fs, the peak frequency up-shifts, accompanied by an increase of the intensity, which is due to a less twisted conformation on this time domain. A broad peak in the 1000–1250 cm^-1 region includes the HOOP and the in-plane bending of H₁₄ in its low- and high-frequency regions, respectively. Although the peak frequency is almost constant in the first 200 fs, the HOOP component exhibits an upward shift as indicated by a change of the spectral shape. This upward shift originates from a slightly pyramidalized structure at the C₁₄ atom due to the increase of the ionic character of the electronic structure induced by the torsion around C₁₃=C₁₄. After 200 fs, the peak bifurcates to the corresponding HOOP and in-plane modes, as also clearly seen in the experimental spectrogram (Kobayashi et al., 2001). The pyramidalized structure formed in the first 200 fs excites the distinct HOOP motion upon the recovery of the planar conformation.

	CONCLUDING REMARKS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
The novel computational technique based on the ab initio description of in situ isomerization dynamics of retinal unveiled in complete molecular detail the anomalously fast event in the primary process of bR's proton pump function. The simulation reproduced the characteristic kinetics and high selectivity of the photoisomerization in bR observed in experiments, i.e., reproduced that the excited-state process starts with in-plane bond relaxation corresponding to the formation of I₄₆₀, and then undergoes isomerization selectively around the C₁₃=C₁₄ bond in an incoherent manner. Comparison of the isomerization dynamics in situ with that in vacuo proves that the high selectivity in bR is due to the protein environment that blocks inhomogeneous reaction paths and guides retinal to the 13-cis isomerization product. The incoherent dynamics originates from trajectories that stay in the I₄₆₀ state for periods of several hundred femtoseconds until they find an escape pathway to the 13-cis product. The simulations also revealed the essential multi-dimensional molecular motions that enhance the electronic transition from the excited state to the ground one and permitted us to identify the molecular origin of the dynamic spectral modulation phenomena observed in femtosecond spectroscopy. Likewise, simulations based on the available x-ray crystallographic structure of rhodopsin (Palczewski et al., 2000) can describe in complete detail the primary process in vision, which in contrast to bR exhibits apparently highly coherent isomerization dynamics (Wang et al., 1994). Without doubt, simulations similar to the present ones, but incorporating a larger segment of retinal in the ab initio quantum chemical treatment, will be carried out sooner or later, but it is expected that these simulations will alter the detailed picture of the bR photoreaction only in respect to quantitative features such as the photoisomerization timescale.

	ACKNOWLEDGEMENTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
We thank H. Kandori and H. Fujisaki for valuable comments.

The work was supported by the National Science Foundation (MCB-9982629), the National Institutes of Health (PHSS P41RR05969-04), and the Human Frontier Science Program Organization. Computer time was provided by the National Science Foundation National Resource Allocation Committee. Molecular images in the paper were created with VMD (Humphrey et al., 1996).

Submitted on March 15, 2003; accepted for publication May 19, 2003.

	REFERENCES

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
Belrhali, H., P. Nollert, A. Royant, C. Menzel, J. P. Rosenbusch, E. M. Landau, and E. Pebay-Peyroula. 1999. Protein, lipid and water organization in bacteriorhodopsin crystals: a molecular view of the purple membrane at 1.9 Å resolution. Structure. 7:909–917.[Medline]

Ben-Nun, M., and T. J. Martinez. 1998. Ab initio molecular dynamics study of cis-trans photoisomerization in ethylene. Chem. Phys. Lett. 298:57–65.

Ben-Nun, M., F. Molnar, H. Lu, J. C. Phillips, T. J. Martinez, and K. Schulten. 1998. Quantum dynamics of retinal's femtosecond photoisomerization in bacteriorhodopsin. Faraday Discuss. 110:447–462.[Medline]

Ben-Nun, M., F. Molnar, K. Schulten, and T. J. Martinez. 2002. The role of intersection topography in bond selectivity of cis-trans photoisomerization. Proc. Natl. Acad. Sci. USA. 99:1769–1773.[Abstract/Free Full Text]

Birge, R. R. 1990. Nature of the primary photochemical event in rhodopsin and bacteriorhodopsin. Biochim. Biophys. Acta. 1016:293–327.[Medline]

Birge, R. R., and T. M. Cooper. 1983. Energy storage in the primary step of the photocycle of bacteriorhodopsin. Biophys. J. 42:61–69.[Abstract/Free Full Text]

Birge, R. R., and L. M. Hubbard. 1980. Molecular dynamics of cis-trans isomerization in rhodopsin. J. Am. Chem. Soc. 102:2195–2205.

Birge, R. R., T. M. Cooper, C. M. Einterz, A. F. Lawrence, M. B. Masthay, A. Schick, C. Vasilakis, C. F. Zhang, and R. Zidovetzki. 1989. A spectroscopic, photocalorimetric and theoretical investigation of the quantum efficiency of the primary event in bacteriorhodopsin. J. Am. Chem. Soc. 111:4063–4073.

Bonai-Koutecky, V., J. Kohler, and J. Michl. 1984. Prediction of structural and environmental effect on the S₁–S₀ energy gap and jump probability in double-bond cis-trans photoisomerization. A general rule. Chem. Phys. Lett. 104:440–443.

Cornell, W. D., P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, Jr., D. M. Ferguson, D. C. Spellmeyer, T. Fox., J. W. Caldwell, and P. A. Kollman. 1995. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 117:5179–5197.

De Vico, L., C. S. Page, M. Garavelli, F. Bernardi, R. Basosi, and M. Olivucci. 2002. Reaction path anlysis of the "tunable" photoisomerization selectivity of free and locked retinal chromophores. J. Am. Chem. Soc. 124:4124–4134.[Medline]

Doltsinis, N. L., and D. Marx. 2002. Nonadiabatic Car-Parrinello molecular dynamics. Phys. Rev. Lett. 88:166402.[Medline]

Dunning, T. H., Jr. 1971. Gaussian basis functions for use in molecular calculations. III. Contraction of (10s6p) atomic basis sets for the first-row atoms. J. Chem. Phys. 55:716–723.

Field, M. J., P. A. Bash, and M. Karplus. 1990. A combined quantum mechanical and molecular mechanical potential for molecular dynamics simulations. J. Comput. Chem. 11:700–733.

Gai, F., K. C. Hasson, J. Cooper McDonald, and P. A. Anfinrud. 1998. Chemical dynamics in proteins: the photoisomerization of retinal in bacteriorhodopsin. Science. 279:1886–1891.[Abstract/Free Full Text]

Garavelli, M., P. Celani, F. Bernardi, M. A. Robb, and M. Olivucci. 1997. The C₅H₆NH₂⁺ protonated Schiff base: an ab initio minimal model for retinal photoisomerization. J. Am. Chem. Soc. 119:6891–6901.

Garavelli, M., F. Negri, and M. Olivucci. 1999. Initial excited-state relaxation of the isolated 11-cis protonated Schiff base of retinal: evidence for in-plane motion from ab initio quantum chemical simulation of the resonance Raman spectrum. J. Am. Chem. Soc. 121:1023–1029.

Garavelli, M., T. Vreven, P. Celani, F. Bernardi, M. A. Robb, and M. Olivucci. 1998. Photoisomerization path for a realistic retinal chromophore model: the nonatetraeniminium cation. J. Am. Chem. Soc. 120:1285–1288.

Gonzalez-Luque, R., M. Garavelli, F. Bernardi, M. Merchán, M. A. Robb, and M. Olivucci. 2000 . Computational evidence in favor of a two-state, two-mode model of the retinal chromophore photoisomerization. Proc. Natl. Acad. Sci. USA. 97:9379–9384.[Abstract/Free Full Text]

Hayashi, S., and I. Ohmine. 2000. Proton transfer in bacteriorhodopsin: structure, excitation and IR spectra, and potential energy surface analyses by an ab initio QM/MM method. J. Phys. Chem. B. 104:10678–10691.

Hayashi, S., E. Tajkhorshid, and K. Schulten. 2002. Structural changes during the formation of early intermediates in the bacteriorhodopsin photocycle. Biophys. J. 83:1281–1297.[Abstract/Free Full Text]

Herman, M. F. 1984. Nonadiabatic semiclassical scattering. I. Analysis of generalized surface hopping procedures. J. Chem. Phys. 81:754–763.

Hou, B., N. Friedman, S. Ruhman, M. Sheves, and M. Ottolenghi. 2001. Ultrafast spectroscopy of the protonated Schiff bases of free and C₁₃=C₁₄ locked retinal. J. Phys. Chem. B. 105:7042–7048.

Humphrey, W., A. Dalke, and K. Schulten. 1996. VMD: visual molecular dynamics. J. Mol. Graph. 14:33–38.[Medline]

Humphrey, W., H. Lu, I. Logunov, H. J. Werner, and K. Schulten. 1998. Three electronic state model of the primary photoisomerization of bacteriorhodopsin. Biophys. J. 75:1689–1699.[Abstract/Free Full Text]

Jorgensen, W. L., J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein. 1983. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79:926–935.

Knowles, P. J., and H.-J. Werner. 1985. An efficient second-order MC SCF method for long configuration expansions. Chem. Phys. Lett. 115:259–267.

Kobayashi, T., T. Saito, and H. Ohtani. 2001. Real-time spectroscopy of transition states in bacteriorhodopsin during retinal isomerization. Nature. 414:531–534.[Medline]

Koyama, Y., K. Kubo, M. Komori, H. Yasuda, and Y. Mukai. 1991. Effect of protonation on the isomerization properties of n-butylamine Schiff base of isomeric retinal as revealed by direct HPLC analyses: selection of isomerization pathways by retinal proteins. Photochem. Photobiol. 54:433–443.[Medline]

Logunov, S. L., L. Song, and M. A. El-Sayed. 1996. Excited-state dynamics of a protonated retinal Schiff base in solution. J. Phys. Chem. 100:18586–18591.

Luecke, H., B. Shobert, H. T. Richter, J. P. Cartailler, and J. K. Lanyi. 1999. Structure of bacteriorhodopsin at 1.55 Å resolution. J. Mol. Biol. 291:899–911.[Medline]

Mathies, R. A., C. H. Brito Cruz, W. T. Pollard, and C. V. Shank. 1988. Direct observation of the femtosecond excited-state cis-trans isomerization in bacteriorhodopsin. Science. 240:777–779.[Abstract/Free Full Text]

Matsui, Y., K. Sakai, M. Murakami, Y. Shiro, S. Adachi, H. Okumura, and T. Kouyama. 2002. Specific damage induced by x-ray radiation and structural changes in the primary photoreaction of bacteriorhodopsin. J. Mol. Biol. 324:469–481.[Medline]

Migani, A., M. A. Robb, and M. Olivucci. 2003. Relationship between photoisomerization path and intersection space in a retinal chromophore model. J. Am. Chem. Soc. 125:2804–2808.[Medline]

Molnar, F., M. Ben-Nun, T. J. Martinez, and K. Schulten. 2002. Characterization of a conical intersection between the ground and first excited state for a retinal analog. J. Mol. Struct.: THEOCHEM. 506:169–178.

Oesterhelt, D., and W. Stoechenius. 1971. Rhodopsin-like protein from the purple membrane of Halobacterium halobium. Nat. New Biol. 233:149–152.[Medline]

Ohmine, I. 1985. Mechanism of nonadiabatic transitions in photoisomerization processes of conjugated molecules: role of hydrogen migrations. J. Chem. Phys. 85:2348–2362.

Palczewski, K., T. Kumasaka, T. Hori, C. A. Behnke, H. Motoshima, B. A. Fox, I. L. Trong, D. C. Teller, T. Okada, R. E. Stenkamp, M. Yamamoto, and M. Miyano. 2000. Crystal structure of rhodopsin: A G-protein-coupled receptor. Science. 289:739–745.[Abstract/Free Full Text]

Ruhman, S., B. Hou, N. Friedman, M. Ottolenghi, and M. Sheves. 2002. Following evolution of bacteriorhodopsin in its reactive excited state via stimulated emission pumping. J. Am. Chem. Soc. 124:8854–8858.[Medline]

Salem, L. 1979. The sudden polarization effect and its possible role in vision. Acc. Chem. Res. 119:12687–12688.

Schobert, B., J. Cupp-Vickery, V. Hornak, S. O. Smith, and J. K. Lanyi. 2002. Crystallographic structure of the K intermediate of bacteriorhodopsin: conservation of free energy after photoisomerization of the retinal. J. Mol. Biol. 321:715–726.[Medline]

Schulten, K. 1978. An isomerization model for the photocycle of bacteriorhodopsin. In Energetics and Structure of Halophilic Organisms. S. R. Caplan and M. Ginzburg, editors. Elsevier, Amsterdam, Holland. 331–334.

Schwartz, B., E. Bittner, O. Prezhdo, and P. Rossky. 1996. Quantum decoherence and the isotope effect in condensed phase nonadiabatic molecular dynamics simulations. J. Chem. Phys. 104:5942–5955.

Song, L., and M. A. El-Sayed. 1998. Primary step in bacteriorhodopsin photosynthesis: bond stretch rather than angle twist of its retinal excited-state structure. J. Am. Chem. Soc. 120:8889–8890.

Song, L., M. A. El-Sayed, and J. K. Lanyi. 1993. Protein catalysis of the retinal subpicosecond photoisomerization in the primary process of bacteriorhodopsin photosynthesis. Science. 261:891–894.[Abstract/Free Full Text]

Swope, W. C., H. C. Andersen, P. H. Berens, and K. R. Wilson. 1982. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: application to small water clusters. J. Chem. Phys. 76:637–649.

Tajkhorshid, E., J. Baudry, K. Schulten, and S. Suhai. 2000. Molecular dynamics study of the nature and origin of retinal's twisted structure in bacteriorhodopsin. Biophys. J. 78:683–693.[Abstract/Free Full Text]

Tallent, J. R., E. W. Hyde, L. A. Findsen, G. C. Fox, and R. R. Birge. 1992. Molecular dynamics of the primary photochemical event in rhodopsin. J. Am. Chem. Soc. 114:1581–1592.

Tittor, J., and D. Oesterhelt. 1990. The quantum yield of bacteriorhodopsin. FEBS Lett. 263:269–273.

Vreven, T., F. Bernardi, M. Garavelli, M. Olivucci, M. A. Robb, and H. B. Schlegel. 1997. Ab initio photoisomerization of a simple retinal chromophore model. J. Am. Chem. Soc. 119:12687–12688.

Vreven, T., and K. Morokuma. 2000. The ONIOM method for the S₁ excited state photoisomerization path of a retinal protonated Schiff base. J. Chem. Phys. 113:2969–2975.

Wang, Q., R. W. Schoenlein, L. A. Peteau, R. A. Mathies, and C. V. Shank. 1994. Vibrationally coherent photochemistry in the femtosecond primary event of vision. Science. 266:422–424.[Abstract/Free Full Text]

Warshel, A. 1976. Bicycle-pedal model for the first step in the vision process. Nature. 260:679–683.[Medline]

Warshel, A., and Z. T. Chu. 2001. Nature of the surface crossing process in bacteriorhodopsin: computer simulations of the quantum dynamics of the primary photochemical event. J. Phys. Chem. B. 105:9857–9871.

Warshel, A., Z. T. Chu, and J.-K. Hwang. 1991. The dynamics of the primary event in rhodopsins revisited. Chem. Phys. 158:303–314.

Weiss, R. M., and A. Warshel. 1979. A new view of the dynamics of singlet cis-trans photoisomerization. J. Am. Chem. Soc. 101:6131–6133.

Werner, H.-J., and P. J. Knowles. 1985. A second order multiconfiguration SCF procedure with optimum convergence. J. Chem. Phys. 82:5053–5063.

Ye, T., E. Gershgoren, N. Friedman, M. Ottolenghi, M. Sheves, and S. Ruhman. 1999. Resolving the primary dynamics of bacteriorhodopsin, and of "C₁₃=C₁₄ locked" analogue, in the reactive excited state. Chem. Phys. Lett. 314:429–434.

This article has been cited by other articles:

				V. I. Prokhorenko, A. M. Nagy, S. A. Waschuk, L. S. Brown, R. R. Birge, and R. J. D. Miller Coherent control of retinal isomerization in bacteriorhodopsin. Science, September 1, 2006; 313(5791): 1257 - 1261. [Abstract] [Full Text] [PDF]

				H. Fujisaki and J. E. Straub Chemical Theory and Computation Special Feature: Vibrational energy relaxation in proteins PNAS, May 10, 2005; 102(19): 6726 - 6731. [Abstract] [Full Text] [PDF]

				T. Andruniow, N. Ferre, and M. Olivucci Structure, initial excited-state relaxation, and energy storage of rhodopsin resolved at the multiconfigurational perturbation theory level PNAS, December 28, 2004; 101(52): 17908 - 17913. [Abstract] [Full Text] [PDF]