Article pubs.acs.org/JPCB
Mechanism by which Untwisting of Retinal Leads to Productive Bacteriorhodopsin Photocycle States Tino Wolter,† Marcus Elstner,*,† Stefan Fischer,‡ Jeremy C. Smith,§ and Ana-Nicoleta Bondar*,∥ †
Institute of Physical Chemistry, Karlsruhe Institute of Technology, Kaiserstrasse 12, 76131 Karlsruhe, Germany IWR, University of Heidelberg, Im Neuenheimer Feld 368, D-69120 Heidelberg, Germany § Center for Molecular Biophysics, University of Tenessee, Oak Ridge National Laboratory, PO BOX 2008 MS6164, Oak Ridge, Tennessee 37831-6164, United States ∥ Theoretical Molecular Biophysics, Department of Physics, Freie Universitaet Berlin, Arnimallee 14, D-14195 Berlin, Germany ‡
ABSTRACT: Relaxation of the twisted-retinal photoproduct state triggers proton-coupled reaction cycle in retinal proteins. Given the crowded protein environments in which the retinal resides, a key open question is whether the retinal relaxation path is governed by the intrinsic torsional properties of the retinal or rather by the interactions of the retinal with protein and water groups. Here we address this question by performing systematic quantum mechanical/ molecular mechanical molecular dynamics computations of retinal dynamics in bacteriorhodopsin at different temperatures, reaction path computations, and assessment of the vibrational fingerprints of the retinal molecule. The results demonstrate a complex dependence of the retinal dynamics and preferred geometry on temperature. As the temperature increases, the retinal dihedral angle samples values largely determined by its internal conformational energy. The protein environment shapes the energetics of retinal relaxation and provides hydrogen-bonding partners that stabilize the retinal geometry.
■
INTRODUCTION Retinal proteins use the light energy absorbed by the retinal chromophore to generate electrochemical gradients, allow ions to pass down electrochemical gradients, or to mediate physiological signals in the cell. Upon absorption of light, the chromophore isomerizes around a specific double bond such that the photoproduct state stores part of the absorbed energy in twisting of the isomerized and other chromophore bonds and in perturbed interactions between the chromophore and its protein environment.1−7 Controlled relaxation of the highenergy photoproduct state8 ensures that the absorbed light energy is used to trigger a productive reaction cycle. How twisted photoproduct states relax toward states that ensure an overall productive reaction cycle is a key issue for all retinal proteins. Here we address this question by performing quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations, reaction path calculations, and computations of excitation energies and vibrational spectra for the bacteriorhodopsin proton pump. The advantage of using bacteriorhodopsin as a model system is that we can compare the results from computations to information from experiments, such as for example, the structures and vibrational fingerprints of the retinal conformers sampled during dynamics and reaction-path computations or the energetics of proton transfers. In the bacteriorhodopsin resting state, bR, the retinal chromophore is all-trans (Scheme 1 and Figure 1A). Upon absorption of light, the retinal isomerizes to 13-cis (Figure 1, panels B−D). The 13-cis, 15-trans K-state photoproduct forms © XXXX American Chemical Society
Scheme 1. Schematic Representation of the Retinal Chain Indicating the Numbers for the Retinal Main Chain Carbon Atomsa
a
The geometry of the retinal was optimized using SCC-DFTB.34
on the picosecond timescale, this being an electronic groundstate intermediate that stores 11−14 kcal/mol of the incident photon energy.6,9 X-ray crystallography at 100 K2,10 and QM/ MM computations6,11 indicated that in K, the retinal Schiff base segment is highly twisted, such that the Schiff base is oriented toward D212 (Figure 1, panels B−C). Such a geometry appears compatible with spectroscopy data, indicating that at low temperatures the K intermediate has a twisted chromophore,12 with the Schiff base NH bond oriented parallel to the membrane plane.13 Special Issue: Photoinduced Proton Transfer in Chemistry and Biology Symposium Received: June 11, 2014 Revised: August 21, 2014
A
dx.doi.org/10.1021/jp505818r | J. Phys. Chem. B XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry B
Article
Figure 1. Retinal geometry in early photocycle intermediates. (A) The all-trans retinal of the bR resting state structure from ref 30. The following atom colors are used: carbon−cyan, nitrogen−blue, and oxygen−red. Water molecules within 6 Å of D85, D212, and K216 are shown as small spheres. For clarity, only part of the retinal and K216 atoms are depicted. (B−D) Illustration of K-state structures from ref 2 (panel B), from ref 10 (panel C), and from ref 29 (panel D). The C12−C13C14−C15 dihedral angle is −56.3° in panel B and −5.5° in panel C. The smaller C13C14 bond twist in the structure from ref 7 as compared to that from ref 2 is associated with an increase by 0.16 Å of the distance between the Schiff base nitrogen and the closest D212 carboxyl oxygen. In ref 29, the retinal polyene chain was modeled in a planar geometry. (E) Protein conformer from ref 22, indicating a highly twisted retinal Schiff base that bridges to D85 via a water molecule. Such a conformer gives an excitation-energy value compatible with the experimental shift of the absorption maximum of L relative to bR25 and proton-transfer energetics compatible with the reaction cycle.22 All molecular graphics were prepared using VMD.85
The molecular picture of a twisted K-state retinal as provided by low-temperature X-ray crystallography2,10 and QM/MM geometry optimizations6,11 may, however, be rather simplistic. Early experiments at 90 K identified at least three different stable conformers of the K-state photoproduct, distinguished by the absorption maximum of the retinal.14 In contrast, at ambient temperatures, time-resolved Raman spectroscopy suggested torsional relaxation and planarization of the retinal polyene chain, 15 and time-resolved FTIR spectroscopy suggested that conversion of an early-K to a late-K-state may involve relaxation of the retinal chain.16 The different geometries of the K-state 13-cis retinal at low versus room temperatures could be due to the protein interfering with the pathway for the rise of K at low temperatures,17 for example, via an irreversible isomerization-induced perturbation,18 the motions of the retinal being more restricted at low temperatures.19
At 150 K, the K intermediate decays into L, the intermediate ready for the transfer of the Schiff base proton to D85. Solidstate NMR data indicate a complex energy landscape for L:20,21 there are four L substates,20 of which only one remains upon increasing the temperature from 150 to 170 K; the other three substates convert back to the ground state.21 The geometry of the retinal Schiff base in the functional L substate derived from the NMR data is characterized by “an exceptionally strong counterion interaction” and torsion around the C15N retinal bond.21 A conformer, the geometry of which is compatible with the NMR data, is that resulting from our previous QM/MM reaction path and molecular dynamics (MD) calculations,22 in which the retinal Schiff base is highly twisted and bridges to the negatively charged D85 via a water molecule (Figure 1E). The water-bridged geometry was derived based on QM/MM computations with second-order Self-Consistent Charge Density Tight-Binding, SCC-DFTB23,24 and confirmed25 via computations with the recently developed third-order SCCB
dx.doi.org/10.1021/jp505818r | J. Phys. Chem. B XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry B
Article
DFTB.26 We demonstrated, using QM/MM approaches, that such a water-bridged geometry gives proton-transfer energetics22 and an excitation energy shift relative to the bR resting state25 that is compatible with experiments. However, the molecular pathway for the relaxation of the twisted K-state intermediate remains unclear. The relaxation path of the twisted K-state is important to determine because knowledge of the structural events along the retinal relaxation path and its energetic determinants may provide clues to general mechanisms for the relaxation of the twisted retinals in microbial rhodopsins. For example, it remains unclear where could a bridging water molecule come from in the preproton-transfer state. As noted before,22 a priori the bridging water may originate from the cytoplasmic bulk, from the cytoplasmic protein region, or from the extracellular side of the retinal Schiff base; movement of a water molecule from its original location in K to the bridging location in L can occur if energetically allowed. QM/MM dynamics at 300 K and reaction path computations22 are compatible with the bridging water molecule being originally located in a cytoplasmic cavity near the retinal Schiff base, as indicated in the crystal structure from ref 27. QM/MM computations (second-order SCCDFTB) using a reactant state with cytoplasmic-oriented retinal Schiff base indicated that translocation of a water molecule from the extracellular to the cytoplasmic side of the retinal Schiff base gives rate-limiting energy barriers of 13−15 kcal/ mol.28 These energy barriers are compatible with the overall millisecond timescale of the reaction cycle, though somewhat large for the microsecond timescale of the K-to-L transition.29 To explore the molecular path for the relaxation of the twisted K-state retinal, here we used QM/MM with third-order SCC-DFTB26 to simulate the dynamics of the K intermediate at nine different temperatures from 50 K to 300 K. For selected structures derived from QM/MM dynamics, we computed excitation energies and vibrational spectra and compared the results to experiments. The QM/MM dynamics simulations indicate that relaxation of the twisted K intermediate proceeds by partial relaxation of bond twists in the retinal Schiff base region such that the Schiff base changes its hydrogen bonding from D212 to T89. Minimum energy computations suggest that reorientation of the retinal Schiff base from D212 to T89 and formation of a water-bridged structure are likely two distinct events, with the Schiff base reorientation occurring faster than water relocation or require relocation of a water molecule from the cytoplasmic side.
groups and waters with atoms further than 14 Å from the retinal Schiff base were constrained with a force, f, inversely proportional to the atomic B factors reported in the crystal structure, f = 0.5(3/2)RT(8π2)/(3B); R is the universal gas constant and T is the temperature in Kelvin. Amino acid residues and waters with atoms between 6 and 14 Å from the retinal Schiff base were restrained harmonically with force constants scaled by the sigmoidal function s(r) = [(r − r1)2(3r2 − r1 − 2r)]/[(r2 − r1)], with r1 = 6 Å and r2 = 14 Å. There were no constraints on groups with atoms within 6 Å from the retinal Schiff base. In the reaction path computations, we considered as the mobile region 827 atoms, including the retinal molecule, one layer of surrounding residues, and water molecules and protein groups whose coordinates likely change during the path. The coordinates of all other atoms were fixed. Potential Energy Function. We used the hybrid QM/MM approach with a hydrogen link atom.36 The QM region consisted of retinal, the side-chains of D85, T89, D212, and K216, and water molecules w401, w402, and w406. The MM region was described with the CHARMM force field31,37,38 and the TIP3P water model,39 using cubic switching with cutoff at 14 Å to smoothly switch off the nonbonded interactions. We set the relative dielectric constant to one. To compute the interactions between the QM and MM parts of the system, we used the QM/MM implementation of the SCC-DFTB method,23,24 with the third-order approximation26 and the 3OB DFTB parameters.40 Test calculations using the previous second-order SCC-DFTB approximation indicated that the method describes well the structural properties of the retinal Schiff base34,41 and the energetics of proton transfers.24,34,42,43 Third-order SCC-DFTB improves the description of the negatively charged ions relative to second-order SCCDFTB.26 The protein models were QM/MM energy-optimized to a RMS energy gradient of 10−3 kcal/mol Å. Reaction Pathway Calculations. To compute the reaction pathway, we used the conjugate peak refinement (CPR) algorithm.44 We refined the CPR path by optimizing its path points with the synchronous chain minimization algorithm45 within the TREK module of CHARMM.31 The computations were performed starting from conformers energy-minimized using QM/MM with third-order SCC-DFTB. Molecular Dynamics Simulations. QM/MM MD simulations for retinal dynamics were performed at 50, 70, 75, 80, 100, 120, 150, 170, and 300 K (Table 1). In each of the simulations, the QM/MM-optimized system was slowly
■
METHODS Protein Models. For the starting protein coordinates in MD simulations of the bR resting state and K-state intermediates, we used the crystal structure from refs 30 and ref 2, respectively. The 9-desmethyl K-state model was prepared from the crystal structure of ref 2 by replacing the 9-methyl group with a hydrogen atom. We used CHARMM31 to construct coordinates for internal atoms that are missing from the crystal structures and capped the protein with the neutral groups −CO−CH3 and −NH− CH3 added to the N (T5) and C (G231) termini, respectively. We used standard protonation states for all titratable residues, except for D96, D115 and E204, which were protonated.32 The protonated E204 models the protonated state of the extracellular proton release group.33−35 In all MD simulations part of the protein was fully flexible and part was subject to harmonic constraints as follows. Protein
Table 1. Summary of the QM/MM MD Simulations simulation
temperature (K)
equilibration (ns)
production (ns)
Simulations on the K-State of Wild-Type Bacteriorhodopsin 1 50 0.3 5 2 70 3 75 4 80 5 100 6 120 7 150 8 170 9 300 Simulations Staring from the K-State Retinal without the 9 Methyl Group 10 50 0.3 5 11 80 C
dx.doi.org/10.1021/jp505818r | J. Phys. Chem. B XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry B
Article
brought to the desired temperature by using a temperature increment of 1K/ps. The heating phase was followed by a first equilibration phase with velocity rescaling and then a second phase with coupling to a Nosé Hoover thermostat.46,47 The simulations were continued with production runs in which the system remained coupled to the Nosé Hoover thermostat, and coordinate sets were saved every 0.1 ps. Computation of Excitation Energies. We computed vertical excitation energies using the ab initio method SORCI48 (spectroscopy-oriented configuration interaction) within the ORCA program.49 The QM atoms were described using the SV(P) basis set50 with additional diffuse s- and p-functions on the carboxylate oxygen atoms of D85 and D212. The usage of SORCI for excitation energies requires several thresholds48 that were set as follows. The weight Tprediag, which gives the configuration state functions built into the reference space, was set to 10−3. The selection threshold Tsel, which differentiates between strongly perturbing versus weakly perturbing subspaces, was set to 10−6 Eh, where Eh is a Hartree (627.5095 kcal/mol). Previous computations on retinal model systems indicated that these Tprediag and Tsel values give a convergence of 0.1 eV for the S1 state.51 Since the polarization of the protein environment can be important,52,53 in computing the excitation energies we replaced the standard CHARMM point charges of the MM atoms with an interactive atomic dipole model which responds to the charge distribution of the S0 and S1 state of the QM region. In this model, the induced dipole moments are first computed from the charge distribution of the S0 and S1 electronic states. The dipole moments are then projected onto a set of monopoles that are added to the standard atomic partial charges of the force field. As a result, the atomic charge distribution is represented by three charges per atom. This computation of dipole moments is done iteratively until the value of the excitation energy no longer changes. A detailed account of the protocol for computing the dipole moments and adjusting the standard CHARMM charges for excitation energies is given in refs 54 and 55. To perform the excitation energy computations, we started from the last coordinate snapshot of each of the nine QM/MM simulations on wild-type bacteriorhodopsin (Table 1) and lowered the temperature to 10 K in 5 K increments; at each temperature value, we performed a brief equilibration of 1.5 ps. The final structure from the dynamics at 10 K was subjected to a geometry optimization to an RMS energy gradient of 10−3 kcal/mol Å2. To assess the differences between geometryminimized structures resulted from quenching of the MD structures, we calculated, for pairs of structures, the backbone root-mean-squared differences (RMSD, in Å) and separately the rmsd for the heavy atoms of the retinal with the covalently bound K216, R82, D85, W86, T89, W138, Y185, W182, and D212. Structures that are within 0.1 Å rmsd are considered the same. This quenching protocol led to three distinct structures that are labeled here as structure-1, structure-2, and structure-3 (Table 2). Computation of Vibrational Spectra. In refs 56 and 57 experimental vibrational spectra for the K intermediate were recorded at ∼80 K and 135 K. To directly compare to these experiments, we computed vibrational spectra as follows. First, starting from structure-1, we performed a set of 100 independent QM/MM MD simulations, 40 ps each, at 80 K. From this set of simulations, we calculated the arithmetic mean of IR intensities. To approximate the nuclear quantum effect,
Table 2. Structural Models Used in Computations of Vertical Excitation Energies structure
original temperature(s)
characteristic structural feature
structure-1
50−120 K
5A
structure-2
150−170 K
Schiff base toward D212 Schiff base toward T89
structure-3
300 K
Schiff base toward T89, w402 displaced
5C
figure
5B
excitation energy 1.96 eV (633 nm) 2.16 eV (574 nm) 2.07 eV (599 nm)
we used a harmonic quantum correction factor as described in ref 58. We followed the same protocol to calculate vibrational spectra for structure-2 and structure-3 at 135 K. To generate reference data and calculate difference vibrational spectra, we used the protocol above to calculate vibrational spectra for the bR resting state at 80 K and 135 K. The difference spectrum between structure-1 and bR at 80 K is labeled here as K-bR-1. Likewise, the two difference spectra at 135 K derived from structure-2 and structure-3 are labeled here as K-bR-2 and KbR-3, respectively (Table 3). Table 3. Summary of the K-bR Difference Spectra Computed spectrum
K-state model
temperature (K)
figure
K-bR-1 K-bR-2 K-bR-3
structure-1 structure-2 structure-3
80 135
5D
Computation of Normal Modes. To assign vibrational modes for retinal in specific conformers, we performed QM/ MM normal mode computations using the VIBRAN module of CHARMM.31 Prior to the normal mode computations, the structures were energy-minimized to an RMSD energy gradient of 10−5 kcal/mol Å. Computation of Histograms and Potential of Mean Force (PMF) Profiles. To characterize the geometry and the dynamics of the retinal at various temperatures tested here, we computed histograms and potential of mean force (PMF) profiles for specific dihedral angles of the retinal polyene chain. All histograms were computed from the last 2 ns of the trajectories (i.e., each histogram was based on 20000 coordinate snapshots); we used a bin width of 0.5°. The PMF profiles for specific dihedral angles were computed from the corresponding histograms according to the equation pmf = −kT ln(y/ymax), where y gives the number of counts for a specific value of the dihedral angle and ymax is the number of counts for the dihedral angle that is sampled most often; T is the temperature at which the simulation was performed.
■
RESULTS AND DISCUSSION In what follows, we first describe the general features of the Kstate retinal at various temperatures. Conformers derived from MD are then evaluated using computations of excitation energies and vibrational spectra. To test whether retinal relaxation may be accompanied by the movement of w402 to a bridging location, we augment the MD simulations by reaction path computations. Temperature Dependence of the K-state Retinal Dynamics. At low temperatures ≤120 K, the retinal chain remains highly twisted, in particular in the segment of the Schiff base (Figure 2, panels A−F); the C13C14 bond is twisted by ∼32°−34° (Figure 3, panels A−F and J−L). But even at these D
dx.doi.org/10.1021/jp505818r | J. Phys. Chem. B XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry B
Article
Figure 2. Temperature-dependent retinal dynamics in wild-type bacteriorhodopsin. (A−I) Coordinate snapshots on the K-state bacteriorhodopsin from simulations performed at (A) 50, (B) 70, (C) 75, (D) 80, (E) 100, (F) 120, (G) 150, (H) 170, and (I) 300 K. From each simulation we show for the retinal and w402 ten equally spaced coordinate snapshots from the last segment of the trajectory, taken each 10 ps. Only the last snapshot (depicted with thicker bonds) is shown for nearby protein and water groups. Hydrogen atoms are shown in white, carbon in cyan, nitrogen in blue, and oxygen in red. At temperatures ≤100 K, the retinal and its surroundings are quite rigid. At 120 K, the retinal chain becomes more dynamic, and at 150 K, it relaxes by reorienting toward T89. At low temperatures, water w402 remains in its location between D85 and D212; at room temperature, it samples a larger region near D85/D212.
cryogenic temperatures, the dynamics of the twisted retinal depends on the temperature: As the temperature increases from 70 to 120 K, the C12−C13C14−C15 dihedral angle undergoes slightly larger fluctuations (Figure 3, panels J−L). At 150 K, the dynamics of the retinal is qualitatively different from that at lower temperatures. Within 320 ps, the retinal relaxes by reorienting toward the cytoplasmic side (Figure 2G), and the twist of the C13C14 bond decreases to ∼14° (Figure 3G). Here, despite the low temperature, the retinal Schiff base is highly dynamic, with relatively large variations of the C12− C13C14−C15 dihedral angle (Figure 3G), such that the Schiff base samples orient toward T89 or toward the cytoplasmic side (Figure 2G). The dynamics at 170 K are similar to those at 150 K (Figures 2H and 3H). At 300 K, the reorientation of the retinal Schiff base toward T89 occurs extremely quickly, during the initial stages of the MD equilibration. The average twist of the retinal C13C14
bond (∼14°) is similar to that at temperatures of 150 and 170 K, though the fluctuations are somewhat larger (Figure 3, panels G−I). W402 becomes dynamic and samples locations on the extracellular side of D85 (Figure 2I). The intriguing temperature dependence of the retinal geometry observed here may be due to the retinal being destabilized due to light-induced twisting in the Schiff base segment,6 its flexibility being facilitated by the small energetic cost for moderate bond twists.59 The reorientation of the twisted retinal toward the cytoplasmic side at temperatures ≥120 K is consistent with the low (
