Retinal isomerization and water-pore formation in channelrhodopsin-2

Retinal photoactivation is a QM process. However, modeling multimicrosecond channel opening requires a computationally efficient classical description. To describe the protein-embedded and covalently bound retinal, we first developed force field parameters for the chromophore following the QM/MM force-matching approach of Doemer et al. (28). The force field parameters and charges of the chromophore were modified so that the forces in an identical classical simulation would match those of a 96-ps QM/MM MD reference simulation. The force constants of the new set of parameters (available at https://gitlab.mpcdf.mpg.de/MPIBP-Hummer/ChR-2) are somewhat softer compared with the original general amber force field (GAFF). Importantly, the relaxed bond lengths in the QM/MM simulations of the protein-bound retinal show a significant bond length alternation (BLA) between single and double bonds, $\text{BLA}={{\displaystyle \sum }}^{\text{}}{d}_{C-C}-{{\displaystyle \sum }}^{\text{}}{d}_{C=C}$. The QM/MM force-matching force field reproduces the QM/MM simulation (BLA_QM/MM-ff = 0.3), which is smaller than the one in the crystal structure (BLA_X-Ray = 0.67, with some uncertainty at 2.3-Å crystallographic resolution) but in agreement with previous density functional theory (DFT) simulations (28, 31) and larger than the GAFF parameters (BLA_GAFF = 0). As atomic partial charges, we used dynamically generated restrained electrostatic potential (D-RESP) charges averaged over the QM/MM trajectory (32) (SI Appendix and SI Appendix, Fig. S1). Compared with the standard GAFF HF/6-31G* charges, we observed a slightly higher polarization in the C13−C15 bond and around C5−C6, similar to what was reported for retinal in rhodopsin (28). Overall, the changes in bond lengths and charges are indicative of a less efficient π conjugation compared with gas phase, consistent with a blue shift in the absorption spectra of the retinal in the protein environment.

To validate the model, we calculated optical absorption spectra of the RSB, as a highly sensitive reporter on the conformation of the chromophore and its interactions with the protein. Previous studies have identified the pitfalls of not considering the thermal sampling and the protein environment when calculating the absorption spectra in rhodopsin (33). Zerner’s intermediate neglect of differential overlap for spectroscopy (ZINDO/S) and time-dependent DFT computations showed the best agreement with experimental data (34). Therefore, we applied the semiempirical ZINDO/S method to a large ensemble of structures extracted from the classical and the QM/MM simulations, including all residues within 4 Å of the chromophore (SI Appendix). The absorption spectra calculated using snapshots from the QM/MM match perfectly with the new force field (SI Appendix, Fig. S2), indicating that the new parameters reproduce the QM/MM ensemble. By contrast, the GAFF simulations show significantly lower excitation energies, due to a poor description of the BLA and the chromophore−protein interactions. We attribute the discrepancy between the absolute values of the experimental and the QM/MM absorption maximum (470 and 512 nm, respectively) to the limitations of the DFT method in reproducing quantitatively the experimental BLA, which results in a red shift in the absorption spectrum (33, 35–37), and to uncertainties in the ZINDO/S method as denoted previously (28). We thus concentrate on the spectral shifts relative to the dark-adapted state, which are in excellent agreement with experiment in the M state (38), and within 10 nm in the K state (Fig. 3).

Our model of the closed state is based on the C1C2 crystal structure (11) as explained in Materials and Methods. During the 4-µs production run, the structure remained stable, with a root-mean-square deviation of the transmembrane helices reaching a plateau value of ∼1.1 Å (SI Appendix, Fig. S3). Water molecules readily entered both sides of the channel between helices H1, H2, H3, and H7 (Fig. 4A), and, after about 20 ns, the number of water molecules reached the equilibrium value of 60 ± 6 molecules (SI Appendix, Fig. S4). Similar results were reported for previous models of ChR2 (6, 31, 39), and in simulations of related opsins (16).

We also observed structural water molecules far from the channel (SI Appendix, Fig. S5). In the so-called DC gate (40), a water molecule entered spontaneously from solution to bridge an H-bond interaction between residues D156 and C128. This water molecule has been observed in previous MD simulations (31, 39). Similarly, a structural water molecule mediates an H bond between the side chain of W223 and the backbone oxygen of S256. This water molecule is also present in the crystal structure of the C1C2 chimera (11) and emulates the water molecule bridging the side-chain nitrogen atom of W182 and the backbone oxygen of A215 in bacteriorhodopsin (bR) (41). Interestingly, four structural water molecules entered into the ChR2 dimer interface to bridge between helices H5, H2′, H4, and H3′, exchanging repeatedly with bulk.

A wide gap in the water density indicates the gate of the putative ion channel. Residues K257 (and the RSB), N258, E90, E123, and D253 form the bottleneck that blocks the leakage of water, protons, and ions (Figs. 2A and 4A). A dense network of interhelical interactions of helix H7 with helices H2 and H3 preclude further penetration of water molecules and ions. E123 (H3) forms a stable salt bridge with the protonated nitrogen, in agreement with the crystal structure of C1C2 (11) and with previous MD simulations of ChR2 (19) and C1C2 (16). Previous QM/MM MD simulations of ChR2 showed alternative RSB−water and RSB−D253 interactions (31). In our simulation, the closed-state structure gives no indication that later, in the P₁⁵⁰⁰/K state, D253 becomes the proton acceptor. This conundrum will be addressed in the following section.

While the paper was under revision, the crystal structures of the dark (closed) state of the wild type and the C128T mutant of ChR2 were published (42). As described in SI Appendix and shown in SI Appendix, Fig. S13, our model accurately captures protein, cofactor, and hydration structures.

The photocycle of ChR2 begins with light absorption and isomerization of the C13−C14 double bond of the RSB moiety. This local conformational change triggers a cascade of events that leads to the preopen channel, the conducting state, desensitized structures, and, ultimately, back to the dark-adapted state (Fig. 1). Previous models of the 13-cis RSB photointermediate have been constructed by structural superposition of the D⁴⁷⁰ (dark) state with the crystal structure of the K state of bR (43) followed by geometry optimization of the 13-cis RSB coordinates in ChR2 (20); by a stepwise tilting of the C13−C14 torsion and MD equilibration steps—all of the RSB torsions were fixed during the equilibration MDs in a sort of big-steps umbrella sampling simulation (19); or by modeling the RSB in the all-trans configuration (44). Here we adopted a more flexible strategy based on enhanced sampling simulations. Starting from various initial structures from our equilibrium trajectory of the D⁴⁷⁰ state, we ran a series of metadynamics simulations (30) as described in SI Appendix. The metadynamics simulations were performed using the ground-state force field parameters for the retinal. This simplification does not include the effect of photoexcitation on the electronic structure of the retinal and on the isomerization energy barrier. We used the C13−C14 torsion as reaction coordinate, and restrained the neighboring torsions with a smoothed-wall potential to avoid the double-torsion flip that leads to the 13-cis−15-syn conformation. This isomer has also been observed in the D⁴⁷⁰ state (45) and corresponds to a closed (nonconducting) channel state.

Computational isomerization of the chromophore led to a P₁⁵⁰⁰ state with a strained (hot) 13-cis−15-anti conformation (Fig. 2B). During the isomerization, the nitrogen of the Schiff base changed its orientation by ∼90°, away from its dark-state partner E123 to form an H bond with D253. A twist of the retinal ethylene chain (C13 to C15) and the C20 methyl group pushes the side chain of F226 toward H6. Interestingly, the change in the orientation of the nitrogen atom of the RSB toward D253 resembles the conformational change of the RSB in the K and L states of bR, as observed by X-ray crystallography (43, 46–48). Two of these structures (43, 47) also show that the C20 methyl tilts toward Y185 (equivalent to F226 in ChR2), and Matsui et al. (43) reported the tilting of the side chain of Y185 as a response to the >1-Å displacement of the C20 methyl group. By contrast, the earliest structure of the K state in bR (49) does not show any of these changes. Overall, the pathway of retinal isomerization in ChR2 resembles that in bR (7).

As another similarity of bR and ChR2, tilting of the Schiff base remodels the H-bond networks in the gate region. In bR, the Schiff base is H-bonding with a water molecule in the dark state, but this water is disordered in the K state (47, 49), or moved (43, 46). The water displacement facilitates the first proton transfer event (48, 49) and provokes a movement of the Schiff base toward the cytoplasmic side of the channel. In the dark state of ChR2, the RSB is H-bonding with E123. Isomerization of the C13−C14 bond and H-bond partner change of the RSB leave E123 free. This transition is accompanied by an influx of water molecules, which facilitates the cleavage of E90−N258 and formation of E90−E123 H-bond interactions. This implies that, as a result of the photoisomerization, H2−H7 and H3−H7 interhelical interactions are replaced by H2−H3 and H7−H7 (RSB−D253) interactions. This is in agreement with FTIR data showing changes in the H-bonding pattern and strength already in the K-like intermediate (50). The differential absorption spectra of the P₁⁵⁰⁰ photointermediate shows good agreement with experiment (Fig. 3), supporting our model.

As a consequence of H-bond remodeling, a preopen pore forms, as indicated by a continuous band of water density between the two vestibules (Fig. 4B). We defined a cylindrical volume centered between helices H1, H2, H3, and H7 that covers the pore between opposite sides of the lipid membrane (SI Appendix, Fig. S6). Then we monitored the evolution of the water density in the z (transmembrane) direction in time. In the dark state, the water density is zero at the center of the channel (i.e., in the central gate), leaving a gap or bottleneck of unsolvated protein residues. The length of the gap was defined as the maximum extent in the z direction of regions below a minimum threshold of water density. SI Appendix, Fig. S7 shows the evolution of the gap length in time for the states modeled. The closed or dark-adapted state (D⁴⁷⁰) displays a large gap that is persistent during the whole trajectory. By contrast, the first photoactivated state shows that, ∼0.6 μs after isomerization, the gap vanishes and stays closed for the rest of the simulation. The active site after the isomerization is more hydrated compared with the dark state. The influx of water molecules modifies some of the H1−H2, H2−H7, and H3−H7 interhelical interactions, inducing conformational changes in the side-chain orientation of several residues of the pore. Cleavage of these interactions, which form a tight mesh of contacts in the dark state that keep the channel closed, let water molecules flow across the pore after retinal isomerization.

We characterized the water permeation by calculating the potential of mean force (PMF) and the position-dependent diffusion coefficient along the transmembrane direction in the pore. We accomplished this by using a maximum likelihood approach based on the transition matrix (51). SI Appendix, Fig. S8 shows the PMF and the diffusion coefficient evolution in time (using time windows of 250 ns of MD simulation). The PMF shows a high barrier at the central gate for the first time window (0 μs to 0.25 µs) that prevents water flux. However, this barrier disappears ∼0.5 µs after isomerization. The diffusion coefficient shows a nearly flat profile, indicating that the calculated PMF is reliable and that there is no other major variable needed to describe the water flux.

Even though water can pass through the gate, we did not observe a spontaneous ion permeation event in any of our simulations. To rule out that this is an artifact of the simulation conditions, we performed a series of MD simulations with increased ion concentration and with the effect of an external electric field. The absence of ion permeation in any of these simulations agrees with the experimental observation of water influx into ChR2 early in the photocycle (<10 μs) into a nonconducting state (15). According to the experiments, the conducting channel is only observed at later stages of the P₂³⁹⁰ intermediate and upon formation of P₃⁵²⁰. Based on the area of the amide-I band, approximately a third of the total of 10 water molecules present in the conducting state enter the channel in the preopen state (15). To compare directly to the experiments, we calculated the number of water molecules that are in H-bonding interaction with the backbone carbonyls of the helices. We observed an increase by 3.2 and 4.4 water molecules in the K- and the M-state MD simulations, respectively (SI Appendix, Fig. S9 and Table S1). As anticipated, only helices H1, H2, H3, and H7 showed a change in hydration. We conclude that the simulation captures the first helix solvation event observed experimentally (15).

The transition from the P₁⁵⁰⁰ to the P₂³⁹⁰ state reflects the proton transfer from the RSB to D253. Indeed, QM/MM MD simulations (including the RSB and the proton acceptor in the quantum region; see SI Appendix) showed that the proton can jump spontaneously to D253. This proton transfer process is remarkably different in bR. In bR, the proton acceptor is D85/E123_ChR2. Proton transfer in bR is mediated by a bridging water molecule and/or T89_bR (47, 52, 53), which establish a hydrogen bond network from the RSB to D85_bR. After proton transfer, the ionic form of D212_bR is stabilized by H-bonding interactions to Y185_bR, while D85_bR is protonated. In ChR2, the proton transfer takes place in the K state, where the RSB is pointing toward D253_ChR2 and does not require any bridging residue or water molecule. After proton transfer, flipping of the RSB toward the cytosolic side makes the process quasi-irreversible. The protein environment is also different, favoring the neutral form of D253 (which cannot H-bond to F226/Y185_bR) and stabilizing the ionic form of E123/D85_bR by H-bonding interactions with E90 and K93. Therefore, the mechanism and the thermodynamics in ChR2 are substantially different from that of bR.

We simulated the P₂³⁹⁰ M intermediate in ChR2 by manually transferring the proton from the RSB to D253 in seven structures taken from our K-like ensemble and extending classical MD simulations for 1 μs each. In six out of seven trajectories, after 200 ns to 400 ns of MD simulation, the unprotonated retinal underwent a conformational change that oriented the nitrogen atom toward the cytosolic side (Fig. 2C). As expected after proton transfer, the absorption spectrum of this conformation shows a large blue shift compared with the dark-adapted state. The difference spectra between the dark and the M states are in perfect agreement with experiment, as shown in Fig. 3. QM/MM MD simulations of the M state show a nearly identical geometry with a more accurate absorption spectrum due to a better description of the BLA (SI Appendix, Fig. S10). The change in the orientation of the Schiff base makes the proton transfer reaction quasi-irreversible, positioning it for eventual reprotonation. In this conformation, the retinal moiety adopts a nearly planar structure, and steric interactions of the C20 methyl group induce a displacement of W223 toward the cytosolic side. This effect was also observed in the X-ray structures of the M state in bR, where the C20 methyl group pushes on W182_bR/W223_ChR2 (47, 48, 54), and in a theoretical study of the C1C2 chimera (20). In ChR2, this displaces the structural water bound to W223. These motions have been suggested as a mechanism by which the photoisomerization of the retinal induces large conformational changes at the cytoplasmic side of helix H6 (55).

The orientation of the RSB toward the cytosol shortens the distance between the Schiff base and the eventual proton donor D156 from ∼10.5 Å in the D⁴⁷⁰ state to ∼9 Å in the M-like state. This distance is still too large for a direct proton transfer, but it could be bridged by a water molecule and C128, as the distance between the Schiff base and the cysteine is also reduced (SI Appendix, Fig. S11). D253 shares H-bond interactions with E90, E123, K93, N258, and the newly entered water molecules in the active site (SI Appendix, Fig. S12).

In the M-like state the water pore remains open. Water molecules diffuse through the channel in a series of jumps between local water sites. Fig. 5 shows the water connectivity in the channel and the path of a water molecule crossing the central gate between E90 and N258 from the cytosolic to the extracellular sides of the membrane. In P₂³⁹⁰, helix 2 preserves the tilting angle and intercysteine distances of the dark state (SI Appendix, Table S2). The distance changes seen in EPR experiments (9, 21) and the tilting angle change in MD simulations (19) should thus take place later in the photocycle, i.e., in the P₃⁵⁰⁰ state and/or after deprotonation of the E90 residue. The structure of the extracellular side of the central gate shows a complex H-bond network between residues D253 (protonated), E90 (protonated), E123, K93, E97, R120, and H249, providing a suitable route for proton transfer to the extracellular side.

The structure of the closed state of C. reinhardtii ChR2 was constructed by homology modeling based on the crystal structure of the C1C2 chimera—PDB ID 3UG9 (11)—and the Swiss-Model software (57). The protonation states of all residues were chosen based on the physiological pH and their local environment. Based on FTIR experimental data (8), E90 and D156 were modeled in their nonionic form. All crystallographic waters were kept, and the system was embedded in a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) lipid bilayer. A ChR2 dimer was solvated with TIP3P water molecules (58) and 100 mM NaCl and enclosed in an 81 × 80 × 166 ų box. The protein was modeled using the AMBER99-SB (59, 60) all-atom force field with side-chain torsional angle corrections (61), and the lipids were modeled with the GAFFlipid force field (62). The chromophore was modeled with an ad hoc set of QM/MM force-matching parameters (28) (SI Appendix). All classical MD simulations were run with the NAMD 2.9 software (63). The geometry of the initial complex was optimized by a series of energy minimization cycles. Then, the system was thermalized and equilibrated to 300 K, followed by an MD simulation of the D⁴⁷⁰state in the NPγT ensemble for 4 μs.

The PMF and the position-dependent diffusion coefficient D(z) of water permeation through the channel were calculated using a maximum likelihood method (64). First, a transition matrix was calculated by dividing the pore into 11 bins and counting all transitions for all water molecules in windows of 260 ns at different time lags. Then, the equilibrium probabilities P_i and the rate coefficients R_j,i were optimized to maximize the likelihood of the observed transition matrix. Finally, the PMF and D(z) were obtained from probabilities and rates. SI Appendix, Fig. S8 shows the transition matrices, the PMF, and the diffusion coefficient before and after channel opening in the simulation of the photointermediate P₁⁵⁰⁰.