Adjacent channelrhodopsin-2 residues within transmembranes 2 and 7 regulate cation selectivity and distribution of the two open states

Channelrhodopsin-2 (ChR2) is a light-activated channel that can conduct cations of multiple valencies down the electrochemical gradient. Under continuous light exposure, ChR2 transitions from a high-conducting open state (O1) to a low-conducting open state (O2) with differing ion selectivity. The molecular basis for the O1 → O2 transition and how ChR2 modulates selectivity between states is currently unresolved. To this end, we used steered molecular dynamics, electrophysiology, and kinetic modeling to identify residues that contribute to gating and selectivity in discrete open states. Analysis of steered molecular dynamics experiments identified three transmembrane residues (Val-86, Lys-93, and Asn-258) that form a putative barrier to ion translocation. Kinetic modeling of photocurrents generated from ChR2 proteins with conservative mutations at these positions demonstrated that these residues contribute to cation selectivity (Val-86 and Asn-258), the transition between the two open states (Val-86), open channel stability, and the hydrogen-bonding network (K93I and K93N). These results suggest that this approach can be used to identify residues that contribute to the open-state transitions and the discrete ion selectivity within these states. With the rise of ChR2 use in optogenetics, it will be critical to identify residues that contribute to O1 or O2 selectivity and gating to minimize undesirable effects.

Channelrhodopsin-2 (ChR2) 2 is a light-activated cation channel endogenous to the green algae Chlamydomonas reinhardtii. In vivo ChR2, together with channelrhodopsin-1 (ChR1), is implicated in phototaxis (1,2). Similar to other microbial-type rhodopsins, ChR2 contains a seven-transmembrane domain (TM) structural motif with the chromophore all-trans-retinal bound to a single lysine on TM 7 through a protonated Schiff base linkage (Lys-257 in ChR2) (1,3). Activation of ChR2 occurs upon photo-isomerization of retinal to 13-cis with blue light ( max ϳ 470 nm) (1). Although structurally similar, ChR2 differs from other microbial-type rhodopsin ion pumps because ChR2 undergoes larger conformational changes in the protein backbone upon photo-isomerization. Like microbial rhodopsins, ChR2 can act as an inefficient proton pump in the absence of an electrochemical gradient (4,5). However, ChR2 functions mainly as a cation channel (6,7). ChR2 conducts a wide array of mono-and divalent metals (Na ϩ , K ϩ , and Ca 2ϩ ) and alkylated amines with a high selectivity for protons (1). The minimum pore diameter of ChR2 has been quantified to be 6.4 Å (8). The unique ion channel properties of ChR2 have made it a popular tool in the emerging field of optogenetics (9 -15).
Activation of ChR2 with blue light gives rise to a transient peak current (I p ), which decays to steady-state value (I ss ) under continuous illumination (Fig. 1A). Once the light is turned off, I ss decays biexponentially back to zero. Subsequent light pulses result in reduced I p but equivalent I SS , suggesting that ChR2 enters a desensitized state in the dark. Initially, the ChR2 photocycle was represented by a three-state model with closed, open, and desensitized states (1). Although the three-state model qualitatively described ChR2 photocurrents, it could not quantitatively describe ChR2 kinetics. More recently, a fourstate model developed under continuous illumination conditions containing two open (conducting) and two closed (non-conducting) states in equilibrium has been able to quantitatively reproduce ChR2 photocurrents (16,17). In this fourstate model, photoactivation transitions ChR2 from the ground (closed) state, C1, to the initial high-conducting (open) state, O1 (Fig. 1B). Under continuous illumination, O1 rapidly decays to the lesser conducting (open) state, O2. When illumination ends, O2 decays to the desensitized ground (closed) state, C2. In prolonged darkness, C2 undergoes a transition to C1, which is due to retinal re-isomerization as well as protonation of one or more residues in or near the retinal binding pocket (3,9,15). Within this four-state model, it is important to note that I p is a combination of mostly O1 with a fraction of ChR2 already in the O2 state, whereas I SS contains a more transiently stable population of O1 and O2. Finally, after light is turned off, ChR2 is primarily in the C2 state. Upon extended time in the dark, ChR2 will be predominantly in the C1 state.
Although the four-state model can quantify the transition states for ChR2 under prolonged activation, the mechanism of light-induced ChR2-mediated cation conductance is unresolved. Analysis of structural and functional experiments has led to the proposal that ion permeation occurs between TMs 1, 2, 3, and 7, which form a large extracellular vestibule filled with water and lined with negatively charged residues (18 -20). More specifically, TM 2 contains five glutamate residues that affect the ion selectivity and function of ChR2 (18,21). Cysteine scanning mutagenesis followed by labeling with methanethiosulfonate derivatives on TM 3 and TM 7 provided experimental evidence that residues comprising these TMs partially mediate ion conductance (8,20). Among these residues is Leu-132, which when mutated to cysteine (CatCh mutant) has large stationary currents and increased Ca 2ϩ permeability compared with WT ChR2 (22). Mutation of Gln-117, which resides at the interface of the extracellular domain and TM 3, has been proposed to disrupt a potential hydrogen bond with Glu-101 at the extracellular interface on TM 2. This hydrogen bond has been suggested to play a role in ion conductance and permeability of ChR2 (20). Furthermore, residues located on the extracellular side of TM 7 and facing TMs 1 and 3 affect ion conductance and kinetics and serve to funnel cations into the pore region (23).
Equally significant for understanding ChR2 function is the gating mechanism for activation. Initially, it was thought that a hydrogen bond was required between Cys-128 and Asp-156 for proper channel function and kinetics (24,25). This hydrogen bond is consistently found throughout microbial rhodopsins, such as the proton pump bacteriorhodopsin, where it contributes to overall protein stability (26 -28). Mutagenesis of Cys-128 exhibited an increased conduction lifetime, which suggested alteration of the potential hydrogen bond between Cys-128 and Asp-156 (24,25,29). However, inspection of the channelrhodopsin chimera C1C2 crystal structure revealed that the distance between Cys-167 and Asp-195 (Cys-128 and Asp-156 in ChR2) is Ͼ 4 Å. This suggests that no direct hydrogen bonding occurs between these residues (19). Furthermore, it has been demonstrated that Asp-156 is the internal proton donor for the Schiff base (30). This indicates that Cys-128 has an undiscovered role in ChR2 function. Currently, gating has been suggested to be controlled, at least partially, through the large conformational changes of the protein and subsequent hydration of helices (31,32).
ChR2 has emerged as a highly successful tool in the field of optogenetics. However, the lack of control over ion selectivity in either open state can produce undesirable indirect effects for in vivo experiments. To rationally design optogenetic mutants to reduce or remove non-selectivity, a better understanding of the molecular basis for the discrete transition between open states and how ChR2 temporally alters ion selectivity is required.
To date, experiments designed to elucidate the unique ion permeation pathway of ChR2 have mainly focused on residues known to be important for homologous retiylidene-containing proteins or through targeting of specific transmembrane domains. Here, our objective was to identify and quantify the role of spatially proximal residues along adjacent TMs that form a barrier to ion permeation. To achieve this objective, we applied steered molecular dynamics (SMD) using a sodium ion to an equilibrated closed-state ChR2 model, based on the fully dark-adapted C1C2 crystal structure, representing the C1 state, to identify residues that comprise a TM barrier. We then generated conservative mutations at identified positions. These mutations were based on inspection of sequence alignments with homologous algal opsin proteins. Photocurrents under extended illumination were measured for WT and mutant constructs. We applied the four-state kinetic model to quantify ChR2 transition parameters from our experimental photocurrent traces (16,33). Finally, photocurrent traces were analyzed to quantify ion conductance, kinetics, and reversal potentials of each mutant ChR2. Using this approach, we have identified proximal residues that contribute to the two conducting states as well as the kinetics and cation selectivity of ChR2.

Results
Channelrhodopsin-2 is part of a distinct subfamily of rhodopsins named channelopsins. These proteins are retinylidenecontaining seven-transmembrane proteins. In contrast to most rhodopsins, which function as ion pumps or receptors, channelrhodopsin-2 can function as both an ion channel and an ion pump. The objective of the experiments described here was to identify and functionally test the role of spatially proximal res-

Channelrhodopsin-2 ion permeation
idues in channelrhodopsin-2 ion conductance by combining steered molecular dynamics and functional studies and analyzing these data using kinetic modeling.
To minimize the system energetic state before SMD experiments, lipid tail groups were melted during a 500-ps equilibration while keeping water and protein fixed. Next, a second 500-ps equilibration was performed with only the protein fixed. Finally, the entire system was released during a 40-ns equilibration. This time frame has been shown previously to be sufficient for performing SMD simulations (34 -38). Fig. 2 shows the overall structure, RMSD of the protein backbone during the 40-ns equilibration, and Ramachandran plot of the equilibrated ChR2 homology model. Following an initial change (Ͻ 0.5 ns) in RMSD, only small fluctuations were observed during the 40-ns equilibration. Based on the RMSD of the protein backbone during the initial 40-ns equilibration, it can be concluded that the homology model of ChR2 is stable. Analysis of residue dihedrals demonstrated that they were within the allowable limits after the 40-ns simulation (Fig. 2C).
Visual inspection of the equilibrated ChR2 model resulted in the observation of bends in both TM 3 and 5, a result of intrahelical hydrogen bonding between a threonine at position i and the carbonyl backbone at i Ϫ 4 (39,40). The hydrogen bond on TM 3 was formed between the hydroxyl side chain of Thr-127 and the backbone carbonyl of Glu-123, whereas the hydrogen bond on TM 5 was formed between the hydroxyl group of Thr-188 and the backbone carbonyl of Tyr-184. Hydrogen bond stability was also observed for Ser-63/Asn-258, whereas no hydrogen bonding between Cys-128 and Asp-156 was observed in our equilibrated structure. The presence or absence of hy-drogen bonds between these sets of residues is consistent with that found in the C1C2 chimera and other ChR2 homology models (19,40,41).

Conduction pathway
Following equilibration, SMD was applied to a sodium ion previously fixed above the center of the electronegative cavity of ChR2. SMD has previously been used to resolve the conduction pathway in the KcsA potassium channel and ion-binding sites of the sodium symporter LeuT and to identify barriers for glycerol permeation in the closed-state aquaglyceroporin, nicotinic acetylcholine receptor, and aquaporin-1 (35,(42)(43)(44)(45). Additionally, SMD has been used to explore binding pathways in monomeric units of higher-order oligomeric proteins, such as bacteriorhodopsin (46). SMD allows for external forces to be applied to molecules on a time scale suitable for molecular dynamics simulations (42). Simulations were performed with monomeric ChR2 in the closed state with retinal remaining in the all-trans conformation. The functional translocation pathway of ChR2 has been shown to be located within the ChR2 monomer (19). SMD was used to explore barriers to ions while closed and was not indicative of the overall single turnover reaction cycle. Sodium entered the pore at the large extracellular vestibule formed between TM domains 1, 2, 3, and 7 (Fig.  3A). This region contains several residues that have been shown to contribute to funneling cations into the pore in addition to negatively charged residues that contribute to ChR2 selectivity (23,47). The Na ϩ -permeation pathway continued in this direction to the side chain of Lys-93, where the applied force increased. From here, the ion entered the TM occlusion site

Channelrhodopsin-2 ion permeation
formed by Ser-63, Glu-90, and Asn-258. Notably, two pairs of stable hydrogen bonds were formed between Ser-63 and Asn-258 and between Lys-93 and Glu-90. When sodium entered this occlusion site, both of these hydrogen bonds were broken and then reformed following sodium permeation. The side chain of Val-86 steered the ion to the inner gate of Tyr-70 ( Fig. 3C) by preventing the ion from permeating out of the pore either through steric occlusion or a hydrophobic interaction. Sodium exited the pore between TMs 2, 3, and 7 and the small intracellular loop between TMs 1 and 2 (IC1), where the applied force reached a second maximum (8 ns in Fig. 3B). IC1 contains highly conserved residues between other microbial-type rhodopsins and channelrhodopsins. No large change in conformation of IC1 was observed.

Force profile of sodium conduction
A normalized force profile was calculated for each individual trial as a rolling average of 500-fs windows and averaged (Fig.   3B). Distinct force maxima were observed at 5.0 and 8.3 ns, which suggests that multiple barriers to ion translocation were encountered during SMD trajectories. Snapshots of this trajectory revealed residues surrounding sodium at time points corresponding to the force maxima. Analysis of our model between 4.5 and 6.0 ns revealed that Val-86, Lys-93, and Asn-258 closely interacted with the pulled sodium ion and prevented it from progressing down the pore (Fig. 3C). The second force maximum was encountered at the interface of the intracellular and TM domains as the ion exited the pore.

Electrophysiology of barrier mutants
We utilized two-electrode voltage clamp electrophysiology to quantify the functional role of the three residues that comprise the TM constriction site. Mutations were guided by a ChR2 sequence alignment with known algal opsin proteins ( Fig.  3D) (48). Val-86 was mutated to alanine and leucine, Lys-93 to asparagine and isoleucine, and Asn-258 to valine and gluta- Figure 3. SMD pulling profile. A, pathway of sodium permeation during a 9.3-ns SMD simulation represented by time points every 1.5 ns with the pulling speed held constant at 5 Å/ns. The trajectory is from a single simulation. B, average force calculated from individual SMD runs using a rolling average window of 500 fs. The force was normalized to the vector used for pulling. The large initial force maximum region corresponded to residues Val-86, Lys-93, and Asn-258. The second peak was observed as sodium was exiting the pore at the interface of the TM domains and intracellular space. C, view down the ChR2 pore from the extracellular side with the three residues of interest shown. The sodium ion is shown in cyan. TM domains of relevant residues are labeled. D, sequence alignment of ChR2 with known algal opsin proteins. Putative transmembrane domains are indicated by black bars, and arrows indicate Val-86, Lys-93, and Asn-258 (ChR2 numbering). Shaded regions indicate sequence conservation. Green, hydrophobic residues; pink, aromatic residues; red, acidic residues; yellow, polar residues; blue, basic residues.

Channelrhodopsin-2 ion permeation
mine. ChR2 mutants were then tested for functional expression in Xenopus laevis oocytes (Fig. 4). The ChR2 proteins containing V86L, K93I, K93N, and N258Q mutations had light-induced current (Fig. 5A). The V86A and N258V mutations resulted in no functional expression. Analysis of Western blotting of ChR2 protein expression revealed that the N258V mutant ChR2 was expressed, whereas the V86A mutant was not (Fig. 5D). The lower molecular weight bands observed may be a post-translational cleaved ChR2 product that is not trafficked to the membrane. All other mutations had expression comparable with WT ChR2. N258V and V86A were not examined further. Interestingly, a complete reduction of peak currents was observed for the V86L mutation upon subsequent light pulses (Fig. 5E).
To quantify changes in ion selectivity for barrier mutants, we measured reversal potentials (E rev ) in 115 mM Na ϩ solutions at pH 7 and 9 and a 115 mM K ϩ solution at pH 9 for both peak and stationary current (Fig. 5 (B and C) and Table 1). N258Q had a significantly more negative E rev at pH 9 when compared with WT. However, E rev was only slightly more negative at pH 7.0, where the photocurrent is carried mostly by H ϩ (47). These large changes in reversal potential demonstrate that N258Q loses permeability for Na ϩ and K ϩ , but not H ϩ . E rev values for both Lys-93 mutants were similar to that for WT and did not appear to alter the selectivity for any of the ions measured. However, when changes in I p and I ss E rev were considered, both mutations at Lys-93 had progressively decreased permeability for Na ϩ compared with WT. This ⌬E rev gives an indication of the permeability of ions as ChR2 progresses through the pho-tocycle (17). Conversely, V86L E rev values were shifted more positive, showing increased permeability for Na ϩ (Ϫ20.3 Ϯ 0.6 mV) and K ϩ (Ϫ17.5 Ϯ 0.8 mV) with respect to WT. Furthermore, we quantified permeability ratios between sodium solutions at pH 7.0/9.0 and sodium and potassium solutions at pH 9.0 ( Table 2). When reducing [H ϩ ] while keeping [Na ϩ ] constant, permeability ratios decreased for all mutants. However, we observed a higher P Na 7.0 /P Na 9.0 for V86L (0.40 Ϯ 0.01) compared with WT (0.28 Ϯ 0.01) and a lower ratio for N258Q (0.11 Ϯ 0.01). This indicated that under steady-state conditions, V86L had a higher H ϩ permeability, whereas N258Q was lower. This correlated with the increased/reduced relative I ss for V86L/N258Q, respectively. A similar trend was observed when comparing P K /P Na at pH 9.0 (WT, 0.82 Ϯ 0.03; V86L, 1.11 Ϯ 0.01; N258Q, 0.72 Ϯ 0.01). K93I and K93N both had higher permeability ratios for Na ϩ at peak currents, but no difference was observed for K ϩ or at stationary current. As a result, these mutants had a significantly lower relative progressive Na ϩ selectivity at pH 9.0 compared with WT.
We also compared steady-state/peak current ratios (I ss /I p ), which are a measure of ChR2 inactivation during prolonged light exposure (49). V86L and K93I had significantly increased ratios compared with WT, indicating less inactivation (Table 1 and Fig. 5A). Conversely, N258Q had small ratios and therefore higher inactivation (18.8 Ϯ 0.8%), whereas K93N had no significant change.
Apparent kinetic rates were extracted from photocurrent traces by exponential fitting ( decay and off ). The decay rate is quantified by fitting photocurrents with a monoexponential open squares, I ss ; green, sodium buffer, pH 7.5; blue, sodium buffer, pH 9.0; purple, potassium buffer, pH 9.0. Blue bars, illumination with 470-nm light. B, peak (blue) and stationary (orange) photocurrents for individual oocytes at V m ϭ 100 in Na ϩ buffer, pH 7.5. Peak currents were not calculated for V86L. Both peak and stationary currents were reduced for ChR2 mutants. No current was observed for the N258V mutant. Black bars, average Ϯ S.E.; n ϭ 7-11.

Channelrhodopsin-2 ion permeation
equation and is largely voltage-independent (49). Similarly, ChR2 undergoes a voltage-dependent biphasic decay when the light is turned off, which is quantified by fitting with a biexponential equation. Both the decay and off rates give insight into the kinetics of ion conductance and global conformational changes occurring during the photocycle. Apparent decay and off rates for all barrier mutants were determined at different holding potentials (Fig. 6 and Table 1). Both K93I and N258Q had faster decay rates compared with WT (Fig. 6A). Likewise, both of these mutants had faster off rates than WT but were less drastic (Fig. 6B). In contrast, V86L decay and off kinetics were severely reduced. Decay rates for V86L were only determined at   Table 2 Permeability ratios of ChR2 constructs Ratios were calculated using differences in reversal potential of the indicated solutions and also between I ss and I p (progressive permeability). Each ratio was calculated using steady-state reversal potentials from individual cells, and results are the average Ϯ S.E. of 7-11 cells. Statistically significant values are denoted with asterisks. *, p Ͻ 0.01. ND, not determined.

Peak current Stationary current
Progressive permeability (I ss ؊ I p ) P Na 9 /P Na 7 P K 9 /P Na 9 P Na 9 /P Na 7 P K 9 /P Na 9 Na pH 7.0 Na pH 9.

Channelrhodopsin-2 ion permeation
Ϫ100 mV because of the loss of I p upon subsequent light pulses (Fig. 5E). These rates were used as guidance for the four-state kinetic model.

Four-state kinetic model
The four-state kinetic algorithm was applied to ChR2 experimental photocurrent traces to quantify the discrete transitions that occur during the photocycle because apparent decay and off rates determined from exponential fitting are dependent on the net result from multiple transitions. Furthermore, this approach can be used to look at the relative occupancy of open and closed states in a time-resolved manner. Here, we used multiparameter optimization to obtain theoretical fits and physical parameters describing ChR2 transition under continuous illumination using Equations 1-5 and 8 (see "Materials and methods") (Fig. 7). The voltage dependence of ChR2 photocurrents is incorporated into the I max parameter of Equation 8. Analysis of electrophysiology experiments was used as guidance for seeding initial parameter values (16,50). Inspection of N258Q parameters identified a faster transition from O1 to O2 and lower ␥ than WT (Table 3). In contrast, the conductance ratio for V86L was much higher than for WT and other barrier mutants, indicating that the large relative stationary currents observed for this mutant were caused by a higher conductance in O2. There were no changes in ␥ for both Lys-93 mutants.
Rates of transitions between states were also calculated. The rates of transition from O1 to O2 (e 12 ) differed for barrier mutants. For N258Q and K93I, this rate was accelerated, whereas for V86L, it was greatly reduced. These theoretical results correlated with the decay rates calculated for each mutant, where N258Q and K93I were faster than the WT and V86L was slower.
Additionally, we investigated the population of conducting states at I p and I ss ( Table 4). The population of the O2 state was similar at maximal currents for WT and all barrier mutants (7-11%). However, large differences in the O1 population were calculated. V86L was made up primarily of O1 at maximal current (77%), whereas WT was significantly less (46%). A similar trend was observed under stationary conditions, where V86L had a much higher contribution from O1 than WT. Thus, the large I ss /I p ratio observed for V86L is a consequence of a higher population of O1 together with higher O2 conductance. In con-trast, N258Q has a small population of O1 (28%) at maximal current. Furthermore, during stationary currents, the O1 population dropped off drastically to 5% with a large rise in O2 population (78%). Therefore, N258Q current is carried largely by the O2 state because of low O1 population at stationary currents and a low ␥ value. Both mutations at Lys-93 did not significantly affect the population of states compared with WT ChR2. Last, we compared I ss /I p ratios with the results of our population analysis. The percentage decrease in O1 population from I p to I ss in our model correlated well with I ss /I p ratios determined from analysis of our electrophysiology results.

Discussion
Channelrhodopsin-2 is a unique dual-functionality protein in that it conducts cations of multiple valencies down the electrochemical gradient as well as pumps protons with low efficiency (0.3 H ϩ /photocycle) in the absence of an electrochemical gradient (5). Activation of ChR2 is the first committed step in phototaxis and has become a facile tool to control excitable cells, such as muscle and neuronal cells, with light. Here, we have used a multimodal approach, including steered molecular dynamics, two-electrode voltage clamp, and kinetic modeling, to identify and quantify the contribution of spatially proximal residues along adjacent TMs to the ChR2 photocycle and cation selectivity.
Trajectories obtained from SMD simulations were consistent with cations entering the putative pore region between TMs 1, 2, 3, and 7. Small resistance was encountered during the beginning stages of pulling, where the ion was located in the large electronegative vestibule. The large force observed from 4.5 to 6.5 ns was correlated to three residues, Val-86, Lys-93, and Asn-258, which directly hinder sodium translocation. A second force maximum was observed at the end of the trajectory and was located at the interface of the intracellular and TM domains.
Protonation changes of residues within the TM domain play a critical role in ChR2 activation, selectivity, and ion conductance (18,30,51,52). Under single-turnover conditions, Glu-90 is deprotonated only during the P480 desensitized state (30). During decay back to the ground state, it becomes reprotonated from an unknown source. Because ChR2 only partially accumulates in P480 under these conditions as part of a branching  K93N; red, N258Q. A, decay rates were determined from monoexponential fitting from I p to I ss at the indicated V m . The decay rate for V86L was determined only at V m ϭ Ϫ100 mV due to the lack of peak current with subsequent light pulses. B, off rates were determined from the biexponential decay after the light was turned off and were dependent on the membrane potential for WT and mutant ChR2s.

Channelrhodopsin-2 ion permeation
photocycle reaction, it remains unclear how the protonation state of Glu-90 affects ion conductance under single-turnover conditions. However, Glu-90 has been shown to be a major determinant for both the selectivity and ion conductance of ChR2 (18).
Inspection of a sequence alignment of homologous algal opsin proteins revealed that Val-86 (in TM 2) encodes hydrophobic residues (Ile, Ala, or Leu) (Fig. 3D). Whereas V86A resulted in a non-expressed protein, the V86L mutation had larger I ss /I p ratios than WT or any of the other mutants presented here, indicating a low level of inactivation and increased O2 ion conductance. This was confirmed through our kineticmodeling data with significantly higher ␥ values than WT. Because there was a complete abolishment of peak currents on subsequent light activation, this is indicative that Val-86 contributes to the recovery back to the ground state.
Analysis of the modeling results also indicates a very slow rate of transition from O2 to C2 (G d2 ) and from C2 to C1 (G r ) when compared with WT ChR2. A similar kinetic effect was observed for the slow-cycling C128T and G224S mutants, where subsequent light pulses had severely reduced I p and had an extended conducting state lifetime, as indicated by the nearly 4-fold increase in off compared with WT (Table 1) (8,25). T off values extracted from photocurrent traces are indicative of the whole-cell ChR2 population closing from two open states. Here, we present the fast values that compose the majority of the fit. It is important to note that the slow and fast values cannot be empirically assigned. Only V86L was determined to have an appreciable difference when compared with WT. We observed that the off rates had a larger dependence on the membrane voltage than WT ChR2. Previously, it has been demonstrated that TM 2 undergoes a large conformational change that is partially responsible for channel gating (31,32). Additionally, Glu-82 has been implicated as a crucial structural element in ChR2 that links movement of TM 2 with the other helices. Combined with our data, it is therefore likely that V86L, which is located one turn above Glu-82, inhibits the movement of TM 2, delaying the opening and closing of the channel. tk;4Interestingly, V86L has nearly identical rates for the forward and back transitions of O1 7 O2, indicating no preference for either state. This position appears to be structurally important for the transition from O1 to O2 under continuous illumination.
Further evidence of the importance of Val-86 to ChR2 structure is the lack of photocurrent from V86A, which is a phenotype similar to the E82A mutant. Not only does TM 2 have strong structural implications in ChR2; it also has a large role in the ion selectivity for ChR2. There are three glutamate residues located on TM 2 (Glu-90, Glu-97, and Glu-101) that are crucial to H ϩ and Na ϩ selectivity (21). V86L had significantly higher permeability ratios for Na ϩ and K ϩ than WT under steadystate conditions. However, under our experimental conditions, we were unable to determine the effects of this residue on the beginning stages of ion flux. Nevertheless, analysis of our

Channelrhodopsin-2 ion permeation
results indicates that charged residues on TM 2 are not the sole contributors to ChR2 selectivity and structural stability. Sequence homology among algal opsin proteins reveals that there is high conservation of asparagine or glutamine at positions homologous to Asn-258 in ChR2. Mutation of Asn-258 to the hydrophobic valine resulted in comparable ChR2 expression to WT but negligible photocurrent under our experimental conditions, suggesting that this mutant is non-functional. However, the N258Q mutation resulted in functional ChR2, but with negatively shifted reversal potentials. The highly reduced permeability for N258Q has also been observed for the homologous mutation in C1C2. At pH 7.0, N258Q is similar to WT permeability, but decreasing the external H ϩ component at pH 9.0 results in hyperpolarized E rev and a reduced permeability ratio. More specifically, when measuring current in K ϩ solution, we would expect E rev to be marginally less than 0 mV if K ϩ is the only permeable ion. However, the internal pH of oocytes is ϳ7.4, and the measured E rev of Ϫ69.5 Ϯ 1.0 mV for N258Q at a pH e of 9.0 indicates outward proton flux with minimal contribution from K ϩ . This is further supported by the similar E rev between Na ϩ and K ϩ at pH 9.0 and the severely reduced permeability ratios between WT and N258Q in identical solutions. We conclude that N258Q severely reduces both Na ϩ and K ϩ conductance and that this position is a major determinant of ChR2 selectivity.
N258Q also displayed accelerated decay kinetics compared with WT, which we attribute to the transition between O1 and O2. Not only does this transition occur much more rapidly than WT, but a larger population of ChR2 occupies O2 in the stationary state. O2 appears to be the kinetically favored state for N258Q. Mutation of the homologous residue in the C1C2 chimera to glutamine also revealed more hyperpolarizing E rev , which was suggested to contribute to chloride permeability (53).
It has been proposed that Asn-258 forms a hydrogen bond with Glu-90 and also Ser-63. Introduction of a negative charge at this position (N258D) results in increased permeability for divalent cations. In our model, we observed a stable hydrogen bond between Asn-258 and Ser-63 and a water-mediated hydrogen-bond network as our simulation progressed. We suggest that loss of Na ϩ and K ϩ flux is a result of the strengthening of the hydrogen-bonding network at the central gate and destabilization of O1. This would prevent larger ions from permeating but allow for H ϩ transport. Furthermore, this residue appears to be critical to ChR2 ion specificity and could be used to predict channel selectivity in channelopsins with unknown function.
Both Lys-93 mutants had reduced photocurrents compared with WT ChR2. Analysis of ChR2 expression revealed that K93N levels were lower than WT, whereas K93I was comparable. The smaller photocurrents of K93I are therefore linked to reduced channel function of single-channel conductance. This is further supported by the lower scaling factor in our kinetic model for K93I, which is representative of single channel conductance. Additionally, Lys-93 mutants had a singular effect on ion permeability for Na ϩ at peak current, leading to a lower progressive selectivity of the channel for Na ϩ at pH 9.0 ( Table 2). This result is surprising because Lys-93 is located in the extracellular vestibule pointing into the pore and is positively charged. Previously, mutation of Lys-132 in C1C2 (Lys-93 in ChR2) to alanine indicated increased K ϩ permeability relative to Na ϩ , but we observed no significant permeability difference for both K93I (P K /P Na ϭ 0.86 Ϯ 0.2) and K93N (P K /P Na ϭ 0.82 Ϯ 0.3) (19). However, under our experimental conditions, there is negligible inward H ϩ flux (pH 9.0), and this may account for the discrepancy. Interestingly, both mutants had decreased progressive selectivity for Na ϩ at pH 9.0. Progressive selectivity can be determined by taking the differences in E rev between stationary and peak currents. Although I ss is not solely populated by the O2 state, these differences in E rev indicate how the ion selectivity changes with respect to time. It appears that removal of a positive charge, not disruption of hydrogen bonding, at this residue is an important factor of sodium selectivity in the early stages of channel opening and ion permeation. Moreover, very little difference for our empirical parameters was observed between Lys-93 mutants and WT. However, the kinetic preference for O2 over O1 was decreased for K93I but maintained for K93N when compared with WT. Lys-93 has been shown to be involved in a complex series of water-mediated hydrogen bonds involving Glu-123, Glu-90, and Asn-258 (40). The K93I mutation effectively removes this hydrogen bond, disrupting the equilibrium between the two conducting states. K93N can retain, partially, the hydrogen-bonding network in this region, which preserves the stability between O1 and O2 observed for WT. Because both K93I and K93N have similar ⌬E rev , the selectivity of ions is probably not related to the hydrogen-bonding network in the region but rather to the charge. Instead, our results indicate that hydrogen bonding in this region is related to the kinetic stability of conformational states.
There is emerging evidence that photoisomerization of retinal upon light activation of channelrhodopsin-2 results in an outward tilt of transmembrane domain 2 (31,54). This conformational change initiates creation of a water-filled pore and subsequent ion conductance. However, the underlying molecular basis of this rearrangement remains unresolved. Analysis of the data described here shows that two residues within TM 2 (Val-86 and Lys-93) are spatially proximal to a residue within TM 7 (Asn-258) in the closed state. Conservative mutagenesis of these results can alter fundamental ChR2 ion conductance properties, including cation selectivity and the relative population of the two open (conductive states). Thus, it is likely that upon photoisomerization, ChR2 conformational changes include separation of Val-86 and Lys-93 from Asn-258, thus enabling ion conductance, and that these spatial proximal residues comprise a pivot point for channel opening/closing.
The use of ChR2 in optogenetics has enabled unparalleled spatio-temporal control of excitable cells both in vitro and in vivo. However, the absence of a true cation specificity leads to many indirect effects in reliable signaling responses. Our experimental approach was able to identify and target residues that contribute to ion selectivity in either O1 or O2. This allows for rational design of ChR2 constructs to minimize indirect effects and establish more effective optogenetic control. Furthermore, this approach can be extended to identifying residues that con-tribute to channel gating and selectivity that would otherwise be overlooked by conventional sequence homology alignments.
Together, by combining steered molecular dynamics with electrophysiological measurements as well as molecular and kinetic modeling, we have identified three spatially proximal residues within the ChR2-conducting pathway that modulate fundamental channel properties, including ion selectivity. Equally, these residues are critical components of the structural integrity of ChR2 and mediate transitions between both open and closed states of the ChR2 photocycle. Analysis of these results can be used to elucidate the molecular underpinnings that define the ion channel functionality of ChR2, potentially predict functionality of homologous proteins that have not been tested directly, and contribute to the development of better constructs for optogenetic approaches.

Building the structural model of channelrhodopsin-2 and molecular dynamics
The ChR2 model was created using the monomeric crystal structure of the channelrhodopsin chimera C1C2 as a template (Protein Data Bank entry 3UG9) and equilibrated as described previously (19,23). All-trans-retinal was added to the structure using the VMD Molefacture plugin using the parameter file described previously (55)(56)(57)(58). Protonation states of titratable residues were determined using the PROPKA plugin (59). All simulations were carried out using NAMD version 2.7 using constant number, pressure, and temperature in addition to periodic boundary conditions (60). Constant temperature and pressure were maintained by utilizing Langevin damping and Langevin piston dynamics. Full-system electrostatics were treated using particle mesh Ewald forces. Short-range van der Waals interactions had a 12 Å cut-off and a switching function at 10 Å. The CHARMM 27 force field was used for all molecular dynamics simulations with 1-fs time steps (61,62). Initially, a 500-ps minimization was performed to melt lipid tail groups while water and protein remained fixed. Minimization was then continued for a second 500-ps simulation in which water was released (34,63). Last, a 40-ns equilibration with the entire system released was performed. The system consisted of 10,056 water molecules and 127 POPC molecules charge-neutralized with 115 mM NaCl for a total atom count of 51,972. The sodium ion used for SMD was fixed during minimization and equilibration. A Ramachandran plot was generated after equilibration to verify that backbone dihedrals were within allowable limits (Fig.  2C) (64).
SMD simulations were initialized from the final time step of the equilibration. Pulling directions were utilized to explore the pathway of sodium conductance through the putative pore of ChR2. Retinal was held in the all-trans conformation during simulations. The sodium ion was tagged for pulling with a spring constant of 5 kcal/mol Å and pulling velocities of 5 and 1 Å/ns. Similar results were observed for pulling velocities of 5 and 1 Å/ns. Forces applied to pull the sodium ion were counteracted by anchoring the protein through the C␣ of select residues that were located at the interface of lipid and water and facing into the lipid to prevent the bilayer and protein from moving along the same axis as the pulled ion. Three trials were performed. All SMD simulations were performed at 293 K.

Microbiology and electrophysiology
Expression and functional measurements in X. laevis oocytes of WT ChR2 and mutants were carried out as described previously (20,23). X. laevis oocytes are a robust heterologous expression system and have been used by a number of groups for ChR2 functional measurements. Oocytes were measured at membrane potentials ranging from Ϫ100 to ϩ40 mV in 20-mV steps in either Na ϩ solution at pH 7.0, Na ϩ solution at pH 9.0, or K ϩ solution at pH 9.0 (115 mM XCl, 2 mM BaCl 2 , 1 mM MgSO 4 , 10 mM HEPES (pH 7.0), or Tris (pH 9.0), where X represents Na ϩ or K ϩ ).

Membrane preparation and Western blotting
Total membrane fractions of Xenopus oocytes were prepared as described previously (8). Briefly, 4 days after mRNA injection, oocytes were homogenized in 200 l of Breaking Buffer (150 mM NaCl, 20 mM HEPES, pH 7.4) supplemented with one tablet of protease inhibitors (Roche Applied Science) and spun at 100 ϫ g for 1 min at 4°C. The supernatant was transferred and spun again at 100 ϫ g for 1 min at 4°C. This process was repeated until no pellet appeared at the bottom. The supernatant was then spun at 16,000 ϫ g for 30 min at 4°C. The supernatant was removed, and the pellet was washed with 200 l of Breaking Buffer and spun a second time at 16,000 ϫ g for 30 min. The supernatant was removed, and the pellet was resuspended in 50 l of Breaking Buffer and 12 l of SDS-loading buffer and incubated for 10 min at room temperature. Samples were spun at room temperature at 16,000 ϫ g for 10 min, and the supernatant was saved.
Samples were separated on an 8% SDS-polyacrylamide gel and transferred to PVDF membrane. The membrane was blocked with 3% BSA in PBST (4 mM KH 2 PO 4 , 16 mM Na 2 HPO 4 , 115 mM NaCl, 0.1% Tween 20, pH 7.4) overnight at 4°C. Membranes were washed with PBST and incubated with anti-HA monoclonal antibody (Cell Signaling Technologies) for 1 h at room temperature (1:1000 in PBST). After washing, the membrane was incubated with 1:1000 HRP-linked anti-rabbit secondary antibody, washed with PBS, developed, and detected by chemiluminescence.

Numerical fitting
The mathematical model was developed based on the fourstate photocycle model (16,33). This model is described by two closed states (C1 and C2) and two open states (O1 and O2) by the rate equations,

Channelrhodopsin-2 ion permeation
dO2 dt ϭ k 2 ϫ C2 ϩ e 12 ϫ O1 Ϫ ͑G d2 ϩ e 21 ͒ ϫ O2 (Eq. 4) C1 ϩ C2 ϩ O1 ϩ O2 ϭ 1 (Eq. 5) where C1, C2, O1, and O2 represent the population of each state. The light-dependent transitions k 1 and k 2 represent the transition C1 3 O1 and C2 3 O2. The back-reactions for these transitions are G d1 and G d2 , respectively. The rates e 12 and e 21 are the forward and back-reactions between the open states, and G r is the thermal conversion from the desensitized C2 to the C1 state. The light-dependent activation rates are determined by the quantum efficiency of photon absorption for retinal, where ⑀1 ⁄ 2 is the quantum efficiency of C1 or C2, F is the photon flux, and p is the sigmoidal activation function described previously (65), with T ChR2 representing the time constant of ChR2 activation. Finally, the current (I) flowing through ChR2 is as follows, where I max ϭ (V m Ϫ E rev ) ϫ g 1 and ␥ is the ratio of channel conductance between O2 and O1 (g 2 /g 1 ). Photocurrent traces were imported, and curve fittings for the four-state model were performed using Mathworks MATLAB R2013. Differential equations were solved using the ODE23 solver in MATLAB. Parameter optimization and minimization was done using the sum of squared errors. For each ChR2 mutant, a global optimization of the parameters was performed using the simulated annealing algorithm (33). Parameters were altered until the fit was a good approximation of the data. Last, a constrained gradient optimization was performed using the fmincon function (50).
Author contributions-R. E. D. conceived and coordinated the study. R. R. and R. E. D. wrote the paper. R. R. performed and analyzed the experiments. Both authors reviewed the results and approved the final version of the manuscript.