Structural Model of Channelrhodopsin*

Background: The channelrhodopsins are widely used for optogenetic application, whereas a structural model was not available before now. Results: Our modeled structure identifies remarkable structural motifs and elucidates important electrophysiological properties of Channelrhodopsin-2. Conclusion: Channel function relies on the unusual properties of helices 1 and 2. Significance: The atom level structure promotes further understanding of function and may guide the engineering of channelrhodopsins for novel optogenetic applications. Channelrhodopsins (ChRs) are light-gated cation channels that mediate ion transport across membranes in microalgae (vectorial catalysis). ChRs are now widely used for the analysis of neural networks in tissues and living animals with light (optogenetics). For elucidation of functional mechanisms at the atomic level, as well as for further engineering and application, a detailed structure is urgently needed. In the absence of an experimental structure, here we develop a structural ChR model based on several molecular computational approaches, capitalizing on characteristic patterns in amino acid sequences of ChR1, ChR2, Volvox ChRs, Mesostigma ChR, and the recently identified ChR of the halophilic alga Dunaliella salina. In the present model, we identify remarkable structural motifs that may explain fundamental electrophysiological properties of ChR2, ChR1, and their mutants, and in a crucial validation of the model, we successfully reproduce the excitation energy predicted by absorption spectra.

The channelrhodopsins ChR1 3 and ChR2 were originally found in the green alga, Chlamydomonas reinhardtii, where they serve as sensory photoreceptors for phototaxis and phobic responses (1)(2)(3)(4). The demonstration that channelrhodopsins act as light-gated ion channels was followed by major advances in neuroscience, where ChRs are used for neuronal depolarization and firing of action potentials simply via delivery of short light flashes (5)(6)(7)(8). ChRs also suggest hope for certain clinical conditions. For example, ChRs might be employed clinically for recovery of photosensitivity in blind mice with degenerated photoreceptor cells, for deep brain stimulation to treat Parkinson disease, for peripheral nerve stimulation and treatment of heart pacing failure, and for modulating synaptic depression and associated fear (9). Interestingly, although most experiments have been conducted with ChR2 because of superior expression and 2-fold higher cation selectivity, ChR1 is in fact the dominant photoreceptor of the alga. It shows a strongly pH-dependent absorption, reduced inactivation, and a faster photocycle. Because of these specific differences, the analysis of both ChRs is of high interest.
ChRs are thought to belong to the protein family of microbial rhodopsins, which include rhodopsins of eubacteria, archaea, cyanobacteria, algae, and fungi. Although much knowledge has been gained by recent experimental studies, a deeper understanding of ChR function, including light gating, ion conduction, and color-tuning mechanisms, clearly requires the knowledge of the three-dimensional structure. The only experimental structure available so far has been resolved by electron crystallography at 6 Å resolution (10) and revealed that ChR2 consists of seven transmembrane helices similar to other microbial rhodopsins. However, because of the low resolution, precise understanding of the position and orientation of each residue is still not available. In fact, the amino acid residues constituting the rhodopsin helices C-G have been identified in the sequences of ChRs. However, the amino acid sequence that constitutes the first two helices has been much less tractable to structural prediction, because sequence alignment of ChRs with other microbial rhodopsins could not be determined uniquely (2)(3)(4).
The disagreement about the arrangement of helices A and B is not only caused by the lower homology with other microbial rhodopsins but also by the charged residues intruding into the transmembrane region. The sequence around the location of the first two helices contains a very high amount of charged amino acids, which are usually thought to be unlikely to be found in a transmembrane helix region. Although such charged residues seem to characterize the channel, they make it difficult to apply a generalized method for structure prediction. In fact, structure prediction programs as SOSUI (11), which is one of the most accurate algorithms for the secondary structure prediction of membrane proteins with accuracy of 80 -90%, predict ChRs as proteins with only four or five transmembrane helices.
To account for the channel function, however, the situation may be quite different. It is interesting to note that there is a sequence pattern, a 7-residue periodicity, which allows the building of helices with the charged/polar side chains aligned along a quasi one-dimensional chain inside the protein, creating a hydrophilic environment that seems very likely to establish the channel. This sequence pattern could determine the structure of the first two helices, and further help for the alignment comes from the recently discovered sequence of a novel channelrhodopsin (DChR) of the halophilic alga Dunaliella salina, which constitutes an evolutionary link between bacterial rhodopsins as Anabaena sensory rhodopsin (ASR) and bacteriorhodopsin (BR), and the ChRs. With the help of these two findings, we are able to propose a three-dimensional structural model of ChR1 and ChR2 based on new experimental data and computational analysis.
In the first step, we build a homology model using the DChR sequence and the characteristic sequence pattern of helices A and B. Homology modeling has been successfully applied for structure prediction when the template sequence has high homology with that of the target protein, as recently shown for G protein-coupled receptors (12,13), which have seven transmembrane helices similar to the microbial rhodopsins. Based on the ASR structural template, we construct a three-dimensional structural model, which is subjected to extensive molecular dynamics (MD) simulations to examine its stability and dynamical features. This model allows interpretation of recent experimental findings and has stimulated further site-directed mutational studies in our group, to examine the role of important amino acids in the active site in more detail. The characteristic three-dimensional structure of ChRs (especially in the region of helices A and B) also allows study of the distribution of internal water molecules, which traverse the ChR itself and may be responsible for channel function. Finally, we are able to localize ion binding sites and propose a model for the active site based on combined quantum mechanics/molecular mechanics (QM/MM) methods that allow computational reproduction of the excitation energy for ChR1 and ChR2 when applying high level post Hartree-Fock methods (multi-reference configuration interaction) for the QM region.

EXPERIMENTAL PROCEDURES
Generation of Three-dimensional Structural Model-In a first step, we constructed a homology model of ChR2 using the SWISS-MODEL web interface (14) based on the three-dimensional structural template of the ASR x-ray structure (Protein Data Bank code 1xio) (15). The residue correspondence between ChR2 and ASR was based on the sequence alignment shown in supplemental Fig. S1. Although the cRNA of ChR2 encodes in total 737 amino acids, we made use only of the amino acid residues 49 -273 of ChR2, i.e. only those who have a counterpart in the ASR sequence. After the replacement of internal water molecules, the protein and water complex was embedded in 1-palmitoyl-2-oleoylphosphatidylcholine lipid bilayer. The modeling details are described in the supplemental Methods.
pK a Calculations-Because the pK a calculations do not necessarily lead to unambiguous results, we used two different models, the multi-conformer continuum electrostatics (16) model using two different values of the dielectric constant (⑀ ϭ 4.0 and ⑀ ϭ 8.0) and the empirical predictions from the PROPKA server (17). The determination of the protonation state of the titratable residues is detailed in the supplemental Methods and supplemental Table S1.
Molecular Dynamics Simulations-The resulting structural model was subjected to extended MD simulations for structural sampling. To avoid an accidental sampling bias, we performed three 50-ns MD simulations with different initial velocities. As a result, we effectively obtained a 150-ns trajectory for the structural analysis and the subsequent QM/MM simulation. All of the MM calculations were performed using GROMACS 4.5 (18), which adopted CHARMM27 potential energy function (19). The protocol and the setup of the MD simulations are detailed in the supplemental Methods and supplemental Fig.  S10.
QM/MM Simulation-In the QM/MM simulations, all lipid and bulk solvent molecules were removed, and harmonic constraints were imposed on the backbone C␣ atoms to suppress global structural disorder.
We employ the CHARMM27 force field as implemented in the CHARMM 36a2 (20) program for the MM part and apply the approximate density functional theory method SCC-DFTB (21) for the QM region. Self-consistent-charge density-functional tight-binding (SCC-DFTB) has been carefully benchmarked and widely applied in our previous studies to retinal proteins (see the supplemental Methods). To mimic the screening effect of bulk solvent and membrane, the Poisson-Boltzmann charge scaling scheme proposed by Dinner et al. (22) was adopted, where the partial charges of the atoms exposed to bulk phase were scaled. The QM region consists of the retinal, the side chains of Lys-257, Glu-123, Asp-253, and Lys-103 and the closest three water molecules to the Schiff base. As a starting point for the QM/MM geometry optimization and MD simulations, we chose the average structure from the classical MD simulations (see the supplemental Methods).
Excitation Energy Calculation-For the vertical excitation energy calculations, we applied multi-reference configuration interaction method called spectroscopy-oriented configuration interaction (SORCI) method (23) implemented in ORCA 2.8 program. In the calculation, the first order interaction space is divided automatically into weakly and strongly interaction spaces. The former is subjected to a variational treatment, whereas the latter is taken into account using perturbation theory. In the SORCI calculations, Ahlrich's SV(P) basis set was employed, and the threshold parameters T pre ϭ 10 Ϫ3 , T sel ϭ 10 Ϫ6 , and T nat ϭ 10 Ϫ6 were used.
Electrophysiology-DNA mutations were generated by sitedirected mutagenesis (QuikChange II; Agilent Technologies) on an hChR2 (amino acids 1-309)-T159C-mCherry vector and confirmed by sequencing. HEK293 cell culturing and whole cell patch clamp measurements were done as earlier described by Yizhar et al. (24). Transient transfection was performed with FuGENE HD (Roche Applied Science) 20 -28 h prior to the electrical recordings. The external measurement solution contained 140 mM NaCl, 2 mM CaCl 2 , 2 mM MgCl 2 , 2 mM KCl, 10 mM HEPES (pH 7.2). The internal solution was composed of 110 mM NaCl, 10 mM EGTA, 2 mM MgCl 2 , 1 mM CaCl 2 , 5 mM KCl, 10 mM HEPES (pH was adjusted to 7.2 either using CsOH or HCl). HEK cell whole cell patch clamping was performed with an EPC 7 (HEKA, Elektronik GmbH, Lambrecht, Germany) amplifier, whereas a Polychrome V unit (TILL Photonics, Planegg, Germany) was used as activation light source, resulting in a final photon density of ϳ1 ϫ 10 22 photons m Ϫ2 s Ϫ1 at 460 nm on the coverslip. For recording the action spectra, only 50% of the light intensity was used. Excitation wavelength for amplitude determination was 460 nm for all constructs.
Two-electrode voltage clamp measurements on X. laevis oocytes were performed under conditions described before (26) by using a Turbo Tec-03X (NPI Electronic, Tamm, Germany). Data acquisition and light triggering were controlled with pCLAMP software via DigiData 1322A interfaces (Molecular Devices, Sunnyvale, CA). The microelectrodes were fabricated by pulling borosilicate glass capillaries (1.5-mm outer diameter and 1.17-mm inner diameter) using a micropipette puller (model P-97; Sutter Instruments, Novato, CA) and filled with 3 M KCl. The resistances of microelectrodes were 0.5-1.5 M⍀. A polychrome2 (TILL Photonics, Gräfelfing, Germany) was used as tunable light source. The light was applied to the oocytes by using a light guide (diameter of 2 mm). The light intensity was 1.5 ϫ 10 22 photons s Ϫ1 m Ϫ2 at the surface of the oocyte. The data obtained were averages of more than three experiments. The solution contained 100 mM NaCl, 1 mM MgCl 2 , 0.1 mM CaCl 2 , and 5 mM citrate (pH 5.0).

RESULTS
Alignment-Throughout this study the residue numbering of ChRs follows that of ChR2. Correspondences among ChRs based on homology and alignment are shown in supplemental Fig. S1.
In the last years, several proposals for sequence alignments of ChRs with BR have been put forward (1,3,4,27). Although the position of the helices C-G can be determined quite well because of reasonable sequence homology, the position of the first two helices has remained unresolved because the alignment is highly underdetermined when based on the BR sequence. To include more information for a reasonable modeling, we included the sequences of sensory rhodopsin of the cyanobacterium Anabaena PCC 7120 (ASR) and the only recently identified channelrhodopsin sequences of Mesostigma (28) and the halophilic alga D. salina. 4 Cyanobacteria are considered evolutionary precursors of green algae and as an evolutionary link between eubacteria and eukaryota (29). Moreover, the intracellular part of ASR is more hydrophilic and contains more water than other microbial rhodopsins (15), which might serve as a structural link between light-driven pumps and lightgated channels. DChR in turn might serve as an evolutionary link between algae and halophilic archaea, which lives in symbiosis with Dunaliella supplemental Table S4. Finally, Mesostigma ChR is of interest because Mesostigma is evolutionary distant from Chlamydomonas, Volvox, and Dunaliella and therefore could highlight residues that are essential for channel function.
The most interesting feature of the ChR sequences is the large number of charged residues in the potential region of the 4 P. Hegemann, F. Zhang, and K. Deisseroth, unpublished observation. first two helices, which are Glu-82, Glu-83, Glu-90, Lys-93, Glu-97, His-100, Glu-101, Asp-103, and Glu-104 in ChR1 and ChR2. Glu-82, Glu-90, and Glu-101 are completely conserved in all six known ChRs (Fig. 1). These residues posed a large problem for many of the earlier alignments, because they were either placed in a loop or inside the helices, being partially oriented toward the membrane. However, their spacing shows a very interesting periodicity: charged residues appear at every seventh amino acid site in the series Glu-83, Glu-90, Glu-97, and Glu-104 (and for Lys-93 and His-100 as well).
In an ␣-helix this periodicity locates the residues on top of each other, forming a chain of charged residues; because the ␣-helix contains 3.6 residues per turn, each amino acid residue rotates 100°in the helix, and the side chains of every seventh residue in an ␣-helix roughly point in the same direction. This remarkable feature strongly suggests a functional significance. First, it allows all charged side chains to have an inward orientation that is clearly favored from an energetic perspective. Charges pointing toward the hydrophobic lipid environment can be assumed to be highly unfavorable. This is highlighted in supplemental Fig. S2, where a helical wheel representation is given. Second, the inward oriented chain of charged groups creates a hydrophilic environment, which is likely to be solvated with water molecules. In fact, as will be shown below, such a structure is readily flooded by waters in a MD simulation, forming a water channel that connects the extracellular with the intracellular bulk phase and is only interrupted in the dark state by the retinal moiety. Note, that such a periodicity can also be observed for the aromatic residues, Tyr-43, Tyr-57, Tyr-64, and Phe-71 in BR, and these first three residues belong to helix B. These residues in BR can then be thought of preventing such a water file in BR, which can be seen as a major difference between the pump BR and the channels ChR1 and ChR2.
To build a structural model for the ChRs based on this finding, we identified several amino acids in helix B of the BR and ASR template structures that point toward the protein interior and are separated by six residues. These are Glu-36, Tyr-37, Ala-40, Ile-43, Pro-44, Ser-47, Ala-50, Tyr-51, and Met-54 in the ASR structure and Phe-42, Tyr-43, Thr-46, Val-49, Pro-50, Ala-53, Met-56, and Tyr-57 in the BR structure. Interestingly, only two alignments satisfy the condition that the charged sites in ChRs coincide with the inside-orientated sites in ASR (Fig. 1,  models A and B). However, it is difficult to decide which of the two alignments, which differ by a sequence shift of seven residues, is a better candidate for helix B. Therefore, we have constructed two model structures using both alignments for further investigation.
In a second step, we considered a model for helix A that shows a reduced homology between DChR and ASR compared with helix B. Although it contains no prominent charged residues, one still can identify a seven-residue periodicity for polar amino acid residues between residues 52 and 73 of ChR2, which we propose to constitute helix A. There are two groups exhibiting such a periodicity: the first one consisting of Ser-52, Ala-59 (threonine in ChR1), Leu-66 (glutamine in DChR2), and Gln-73 and the second one consisting of Gln-56, Ser-63, and Tyr-70 (see supplemental Fig. S1). The residues should prefer an inward orientation similar to the charged ones in helix B. The side chains of helix A in BR with inward orientation are numbered as: 9, 13, 16, 20, 23, 24, 27, and 30 (4, 8, 11, 15, 18, 19, 22, and 25 in ASR). With this information, the alignment of helix A can be determined uniquely because the two groups of residues in ChRs can be nicely matched with the corresponding ones from BR and ASR.
Three-dimensional Structure of the First Two Helices-Based on the ASR crystal structure as template, we have built two three-dimensional structural models of ChR2 based on the two different alignments for helix B (Fig. 2 and supplemental Fig.  S3).
Model A at first sight seems favorable compared to model B because the gap region between helix B and helix C is smaller, which is generally preferred by conventional alignment tech-  MARCH 2, 2012 • VOLUME 287 • NUMBER 10 niques and all glutamic acid residues from Glu-82 to Glu-104 reside in the helical region, which is not the case for alignment B (Fig. 2 and supplemental Fig. S3). Further, Glu-82, Ile-89, and Asp-103 would be conserved for all ChRs (except VChR2 with Val-89) and ASR when using model A.

Atom Level Structure and Electrophysiological Findings
Considering an evolutionary relationship between ASR and DChR, the latter may be considered as a link between the ChRs and ASR. This relation would support model B: DChR then would have a high homology with the ChRs but also substantial homology with ASR identified in three doublets Leu-91/Val-92, Ile-98/Trp-99, and Ala-106/Met-107 in DChR and Leu-38/Val-39, Ile-45/Trp-46, and Ala-53/Met-54 in ASR (Fig. 1).
Global Structure and Active Site-From the homology model, we generated a three-dimensional structural model embedded into a model membrane surrounded by bulk water, and MD simulations have been performed to probe the dynamics and stability of the model structures (for details see the supplemental Methods). The most obvious residue to study is certainly Glu-123, which is close to the Schiff base in both alignments. Glu-123 is the substitute for D85 in BR, which is the proton acceptor during deprotonation of the Schiff base. Glu-123 has been replaced in earlier studies by the protonated glutamine (E123Q), by the hydrophobic alanine (E123A), or by the threonine (E123T). These mutations accelerate the photocycle and shift the action spectrum to longer wavelengths by 20 nm for Thr and Gln and 9 nm for Ala (30,31), which is in line with both arrangements of Fig. 2. Next to Glu-123, a positively charged residue is located close to the Schiff base but different at different places in both models. It is Lys-93 in the model A and His-100/Lys-103 (in ChR1/ChR2) for model B.
In model A (Fig. 2), the active sites are equal in ChR1 and ChR2, consisting of Glu-90, Glu-123, Lys-93, and, a little further apart, Glu-97, whereas they are very different according to model B (for ChR1, see supplemental Fig. S3). Here, for ChR1 the protonated His-100 in the center is surrounded by the amino acids Glu-97, Glu-101, Lys-103, and Glu-104 on helix B and Glu-123 on helix C, as shown in supplemental Fig. S3. This structure can be seen in the three other channelrhodopsins ChR2, VChR1, and VChR2 as well, although Glu-104 is mutated into serine in VChR1, and the positively charged site shifts from His-100 to Lys-103 in ChR2. This highly polar structural motif would be critical for channel gating and ion conductance but may also contribute to the blue shifted absorption of both ChR1 and ChR2 when compared with BR. The latter is investigated in more detail below.
To understand this residue network in more detail, we mutated several amino acid residues of ChR2 and tested the photocurrents in HEK cells. Mutations were made on the ChR2-T159C background because the currents were very large (31), and 159 is far from the active site (discussed later in Fig. 5).
First, we considered three mutations, namely F100H, K103D, and N104E. According to model B, this would convert the ChR2 active site into the ChR1 counterpart (Fig. 2, model B), whereas in model A these residues would be located at the lower part of helix B, being closer to the potential channel entrance. The protein still showed 55% activity (Fig. 3), and although the mutation of the charged residue in the vicinity of the chromophore generally has a significant impact on the spectrum, the observed max shift was only several nm (supplemental Fig.  S4). This indicates an important role of these residues for conduction, but not for color tuning, which at a first glance seems to be more consistent with model A than for B. Note, however, that the total charge in model B ChR1 and ChR2 in the active site is conserved, which could account or the missing color shift upon the triple mutation. Nevertheless, model A seems to be favorable in this respect.
In separate experiments, we removed the charges from Glu-90, Glu-101, and Glu-123 and replaced them with Gln. The peak current amplitude dropped by 64, 76, and 50% respectively (Fig. 3). In E123Q, the action spectrum was red-shifted by 25 nm (Fig. 3), whereas the other mutants showed only small shifts of a few nm (supplemental Fig. S4). The data suggest that all three residues support function but are individually not of critical importance. However, multiple replacements in the form of E97Q/E101Q or E97Q/E101Q/123Q dropped the activity to 5 or 2% (Fig. 3 and supplemental Table S2). Fig. 4 shows the cytoplasmic side of the three-dimensional structure based on the alignments A and B. Interestingly a number of titratable residues are present at that location, including Glu-90 and His-134, which have been intensively studied by several groups. They are known to be involved in H ϩ selectivity and H ϩ release. Substitution of Glu-90 by Gln or Asp or His-134 by Arg or Asn inhibited H ϩ conductance and promoted Na ϩ selectivity (27,32,33), suggesting that both residues are responsible for H ϩ transfer and release into the medium or part of the internal selectivity filter. Ritter et al. (34) and Radu and co-workers (35,36) monitoring FTIR differences between the dark state and the intermediate P500 reported that ChR2 has two indicative band patterns around 1700 cm Ϫ1 , implying the existence of two protonated carboxyl residues at least. The first band with higher frequency, 1747(ϩ)/1739(Ϫ) cm Ϫ1 , is assigned to Asp-156 (35,36), being identical to protonated Asp-115 in BR. The second one, 1729(ϩ)/1739(Ϫ) cm Ϫ1 , is assigned to Glu-90 because it is not observed in the E90Q mutant (34,35).
The lower frequency of Glu-90 implies stronger hydrogen bonding, which could be to a hydrophilic exposure of the carboxyl side chain. This would support alignment B, where Glu-90 is located on the edge of helix B. However, it could also mean that Glu-90 is strongly interacting with water molecules supposedly located in the channel. As shown for model A in the supplemental Fig. S6, water distributes in the vicinity of Glu-90 during the MD simulations. Therefore, it is likely that Glu-90 forms hydrogen bonds with the inner water or the backbone atom of helix G, lowering the frequency of this vibration.
Sugiyama et al. (32) substituted the glutamates of helix B individually by alanine and found strong inhibition of photore-ceptor currents for Ala-82, Ala-83, Ala-97, and Ala-101 but not for Glu-90. The same finding has been reported by Ruffert et al. (27), proposing a model for helix B that is very similar to our alignment A, because they propose Glu-90 to be located in the middle of helix B. Replacing Glu-90, Glu-97, and Glu-101 with lysines resulted in a large reduction of conductance, highlighting the importance of these residues for channel function. Note that in alignment B, Glu-82/Glu-83 are directed toward the bulk water phase, where they could act as antenna, attracting cations toward the channel. On the other hand, the model building based on the ASR x-ray structure leads to a shorter helix B than when using the BR structure. For this structure, Glu-82/83 could readily be part of helix B. However, to recover the helix during classical MD, longer sampling times would be required, which is beyond the scope of the present manuscript.
Of interest might be K93Q, which does not show Na ϩ currents nor inhibition of the H ϩ conductance by Na ϩ (50). In both alignments, this residue could be part of the potential water channel, with model A being located in the extracellular side, whereas in model B it is located in the cytoplasmic side.
Note that in the present model B, His-265 is located between Glu-90 and His-134, both of which are responsible for cation selectivity (Fig. 4, inset). Therefore we inferred that His-265 might also play an important role. To verify this proposal, we mutated this histidine and measured the photocurrent, which is found to be drastically reduced with regard to the wild type  MARCH 2, 2012 • VOLUME 287 • NUMBER 10 (supplemental Figs. S5 and S6). In model A, His-265 would be still located in the channel but would not interfere directly with Glu-90 and His-134.

Atom Level Structure and Electrophysiological Findings
Water Distribution-As a result of the extended MD simulations for both models, a large number of water molecules are found in the interior of ChR2 (see supplemental Fig. S7), which is quite different from the situation in other microbial rhodopsins. The solvent molecules are mainly distributed along helix B caused by the polar and charged residues regularly arranged along this helix. This seems to be reasonable because these two helices differ substantially from those in other rhodopsins and therefore are likely responsible for the channel function. It is also noteworthy that solvent molecules are more populated in the extracellular side than the cytoplasmic side in common with other microbial rhodopsins. Although the distribution is separated by retinal in the dark state, the retinal isomerization triggered by light absorption would cause structural changes at the active site and enable cations to pass through.
DC Gate-In agreement with previous studies (35,36), the present models both locate Cys-128 on helix C and Asp-156 on helix D in the center of the transmembrane region (Fig. 5). Mutation of Cys-128 into Thr/Ser or Ala and mutation of Asp-156 into Ala slows down the photocycle kinetics by 2-4 orders of magnitude (37,38). According to these original electrical studies and subsequent FTIR experiments for these mutants (36), Asp-156 in ChR2 is probably protonated. Thus it has been implied that Asp-156 forms a hydrogen bond with Cys-128, where Cys-128 might be a donor against a carbonyl of the side chain of Asp-156. Furthermore, this arrangement would explain that the DC hydrogen bond would play a significant role in controlling the channel kinetics, justifying the term "DC gate" (36).
However, the present MD simulations exhibit that Cys-128 is closer to Trp-124 than Asp-156 (see supplemental Fig. S8), implying that Cys-128 is probably a proton donor of a hydrogen bond with the oxygen atom of the backbone of Trp-124, conflicting with the FTIR model. To investigate this in more detail, we also performed a 20-ns MD simulation of the C128T mutant, which shows that the threonine mutant forms an even stronger hydrogen bond with Trp-124 than the wild type (supplemental Fig. S8), which again is in conflict with the FTIR model. Interestingly, the BR x-ray structure supports the present computational model, indicating that Thr-90 (Cys-128 in ChR2) forms two hydrogen bonds with Asp-115 (Asp-156 in ChR2) and Trp-86 (Trp-124 in ChR2). Although the hydrogen atoms are not resolved in the x-ray structures, Thr-90 should obviously be the proton acceptor for Asp-115 and a donor to Trp-86. This hydrogen bond between Thr-90 and Trp-86 is equal to a motif discussed by Ballesteros et al. (39), where the side chain of threonine in the g Ϫ conformation ( 1 ϭ gauche Ϫ ) forms the intrahelical hydrogen bond with the i Ϫ 4 carbonyl oxygen. Although there is only a little information about intrahelical hydrogen bonds of cysteines available, we concluded that both the wild type and C128T mutant should form an intrahelical hydrogen bond.
Next, we computed the distance distributions of the sulfur atom of Cys-128 with two oxygen atoms of the side chain of the protonated Asp-156, O␦1 (CϭO), and O␦2 (COH) (supplemental Fig. S8). In the C128T mutant, O␦2 is closer to the side chain of the residue 128 than O␦1, which means that Asp-156 may not be a proton acceptor but a donor in the hydrogen bond. Although it might be also the case for the wild type, the side chain Asp-156 might interact with the backbone of Cys-128 rather than with its side chain, possibly because cysteines are less polar than threonines. In addition, the MD average structure implies that solvent molecules come into the vicinity of Asp-156. In fact, in the BR x-ray structures, a water molecule forms hydrogen bonds with Asp-115 and Leu-87 connecting helix C and D (40). Thus, it is possible that the missing solvent stabilizes the hydrogen bond in the wild type ChR2. However, this issue may be difficult to describe with pure MM methods, because they can fail to reproduce subtle details of complex hydrogen-bonded networks, as reported, for example, for the BR active site (see below). Therefore, more extended QM/MM simulations may be required that are beyond the scope of the present work.
The arrangement with Cys-128 connected to Asp-156 and the Trp-124 backbone would explain why the mutation C128A has a stronger effect than the mutation D156A and why the double mutation C128S/D156A results in a stabilized step function opsin with an open state life time of 30 min (24), which would not be expected for the sole connection of Cys-128 to Asp-156.
Combined QM/MM Simulation-Because classical force field methods can be problematic for complex hydrogenbonded active site structures (see the supplemental Methods), we performed QM/MM simulations using the approximate DFT method SCC-DFTB (22) in the QM region (41).
The active sites of models A and B are quite different (Fig. 6). In the case of model A, there is a larger similarity to the situation in BR, where a prominent pentagonal hydrogen-bonded structure has been found (42), consisting of retinal, Asp-212, Asp-85, and three water molecules. The important structural detail is the interaction of retinal with a water molecule, although there is no salt bridge with one of the counterions. A QM/MM geometry optimization starting from the MD average structure shows that this is similar for model A, although a partial salt bridge with one of the counterions is formed. The active site in model B is not constituted by a strong pentagonal H-bonded network as in BR, because Lys-103 comes into the active site, which seems to change electrostatic environment significantly. The distance between the carboxyl side chains of the counterions in ChR2 (5.7 Å) is smaller than that in BR structure (6.1 Å; Protein Data Bank 1c3w), which may be caused by the replacement of aspartic acid by glutamic acid that has a longer side chain (supplemental Table S3). Note that the Asp/ Glu replacement from BR to ChR also changes the hydrogen bond ASP-85-Thr-89 compared to Glu-123-Thr-127 in ChR2; Thr-127 forms an intrahelical hydrogen bond with the carbonyl oxygen of Glu-123 (Fig. 6) instead of with the Glu-123 side chain.
To probe the stability and the flexibility of the active site, we performed a QM/MM MD simulation (for 1 ns). Although the direct interaction with the Schiff base and the closest water molecule are conserved, Glu-123 seems to be more flexible than D85 in BR, whereas Asp-253 appears as similarly rigid as Asp-212 in BR. For Glu-123, we observe a change in the hydrogen bonding pattern during the MD, i.e. the positions of O⑀1 and O⑀2 interchange during the simulations highlighting the floppiness of the hydrogen-bonded network.
Excitation Energy Calculation-The optical absorption spectrum of ChR2 exhibits several peaks ( max ) of 447, 470, and 410 nm in the dark-adapted state at neutral pH (34). Because it is difficult to define the exact maximum in the absorption spectrum of ChR2 at neutral pH, instead we adopt the max of the action spectrum consisting of a single peak, namely 2.74 eV (452 nm) (Fig. 3 and Table 1) as the reference for the excitation energy calculation (43).
The computed value for the QM/MM optimized structure is 2.64 eV (469 nm) for model A and 2.67 eV (464 nm) for model B. Therefore, both models led to very similar excitation energies, meaning that they cannot be distinguished by such calculations. Because helices C-G are the same in both models, it seems that the different positions of the charged residues in helix B does not have a large impact on the excitation energies, which is a quite reasonable result.
The results have to be compared with the values for BR and SRII as shown in Table 1. QM/MM estimates of excitation energies rely to some degree on an error cancellation, even when using very accurate QM methods like the SORCI method applied in the present study, as we have discussed in detail in our previous work (41, 44 -46), leading to some systematic trends in the computed values. In short, we expect excitation energies to be slightly blue-shifted, mostly because of the neglect of protein polarization.
As can be seen from the examples of BR and SRII from Table  1, the computed value of 2.64/2.67 eV seems to be an excellent estimate of the ChR2 max . This indicates that our computed structure gives excitation energies in semi-quantitative agreement with the experimental ones, although the calculations cannot be used to distinguish between the two models. In general, proximity of the negative charge to the Schiff base stabilizes the electronic ground state and destabilize the excited state leading to blue shift and vice versa. Therefore, the mutational max shifts agree well with the present models where both glutamic acids are supposed to be deprotonated.

DISCUSSION
In this work, we proposed a three-dimensional model of ChR2/ChR1 based on the two alignments in Fig. 1, which result from the observation of the seven-residue periodicity of conserved charged amino acids in the ChR sequences. This pattern is therefore likely to have a functional significance, and indeed, it allows for the formation of a hydrophilic channel, which has been confirmed by our MD simulations (supplemental Fig. S7), because the side chains of Glu-90, Lys-93, Glu-97, Glu-101, and Lys-103 in ChR2 are oriented to inside the protein. Note that these residues have been shown by experiments to be significant for the channel function. However, we were not able to determine the alignment unambiguously; therefore we probed the two possible alignments for helix B, which are similar to the previous suggestion from Ruffert et al. (27) regarding helix B.
The experimentally investigated mutations in this work seem to be consistent with both models A and B as well, leaving the determination of the exact structure for future work. The triple mutant F100H/K103D/N104E has only a small impact on the optical spectrum, a finding in favor of model A at first sight, because model B places these residues prominently into the retinal binding pocket, where such drastic changes are thought to lead to sizable effects. However, the strong H-bonded network may allow for charge shifts upon mutation; therefore this issue is not easy to decide.
A further hint results from the QM/MM MD simulations, which show for model A a much tighter bound binding pocket, whereas for model B larger fluctuations and flip-flop motions occur. From a functional perspective, the effect of isomerization, which leads to channel opening, may be easier to understand on the basis of a tighter bound active site. Model B, on the other hand, has been proposed because of the homology of ASR and DChR, considering their evolutionary relationships. Also the location of Glu-90 at the protein surface in model B more easily explains the red-shifted IR absorption. However, this may also be due to interaction with waters in the channel of model A.
Helices C-G are the same in model A and B; therefore, the discussion of the so called DC gate may be independent of a particular alignment. Our MD model leads to different conclusions regarding the interaction of Cys-128 and Asp-156 than the proposed model from FTIR measurements, in particular regarding the C128T mutant. The interactions of threonines at site i with backbone carbonyl oxygens at site i Ϫ 4 are of particular interest, because they may distort the helical structure (39). This is also found with the side chain of Thr-127 on helix C (Fig. 1), which also assumes a g Ϫ conformation forming the intrahelical hydrogen bond with Glu-123 according to our MD simulations, thereby distorting the helical structure leading to lower helicities at residues 126 and 127 in the wild type (supplemental Fig. S9). Furthermore, we also found that the C128T mutation additionally reduces the helicity at residue 126. Because proline generally disrupts a helical structure, Pro-129 also destabilizes helix C (supplemental Fig. S9). As a result, helix C is bent and flexible around the structural motif consisting of Thr-127, Cys-128, and Pro-129.
It should be also emphasized that Thr-127 has an intrahelical hydrogen bond with the backbone of Glu-123 in the active site (Fig. 6). Although Thr-89 in BR (Thr-127 in ChR2) interacts with the Asp-85 side chain (Glu-123 in ChR2) in the dark state (Fig. 6), it has been proposed that Thr-89 forms an intrahelical hydrogen bond with Asp-85 in the O-intermediate state (47). Therefore, it is also possible that the corresponding hydrogen in ChRs may alter the bonding pattern during the photocycle.
Taking into account the flexibility and the connection to the active site, the intrahelical hydrogen bonds may be important for the channel on and off gating during the ChR photocycle. Note that the three residues 127-129 are highly conserved among ChRs (Fig. 1), supporting the present model. In addition, it has been reported that C128X mutations dramatically slow down the photocycle (37,38).
Coming to Thr-159 on helix D, it has been recognized that the mutation T159C causes 10 times larger photocurrents compared with the ChR2 wild type in Xenopus oocytes (31). Moreover, the T159M mutation alters the ion selectivity of the channel in favor of K ϩ with strong consequences on the reversal potential in living cells (31). In the present model, we could not find any polar side chains that could interact with Thr-159. Instead, the side chain of Thr-159 in the g Ϫ conformation forms the intrahelical hydrogen bond with Ser-155, causing a distortion of helix D as shown in Fig. 5. Because Thr-159 is located close to the DC gate, it is likely that Thr-159 and its hydrogen bond act cooperatively with the DC gate. This is in line with experimental findings that modification of the DC gate has a reduced influence in VChR1 or the C1V1 hybrid (24), which naturally contain a cysteine at the 159 position. Furthermore, we can see intrahelical hydrogen bonds between Ser-155-Gly-151 on helix D and Thr-188 -Tyr-184 on helix E in the vicinity of the retinal binding pocket (Fig. 5). Therefore, Ser-155 and Thr-188 may contribute to the function in the same way that the DC gate and Thr-159 do. To prove this, we measured the photocurrent of ChR2 mutants T188A and T188C. Both constructs showed normal targeting and expression levels as evaluated from enhanced cyan fluorescent protein fluorescence. The mutations result in a drastic reduction of the photocurrent showing the importance of this structural motif for the ChR function (supplemental Figs. S5 and S6). T188A showed hardly any photocurrents, whereas photocurrents of T188C were approximately five times smaller than that of wild type. This might be related to the fact that cysteine can still form an intrahelical hydrogen bond, different from alanine, although it can expected to be weaker than that of threonine.
Based on the electron-microscopy structure (10), it has been proposed that the ChR2 channel is located at the dimer interface on the 2-fold axis, lined by transmembrane helices C and D. To investigate this issue in more detail, we performed 30-ns simulations of ChR2 using the dimeric displacement based on the cryoelectron microscopy result.

TABLE 1 Excitation energy (eV) and wavelength in parenthesis (nm)
In SORCI, the QM region consists of the entire chromophore and lysine side chain. With the 2-fold axis between helices C and D, as suggested by cryoelectron microscopy, the simulations do not show any steric hindrance and do exhibit the densely packed hydrophobic residues between the proteins. Note that the alignment about the two helices is uniquely determined. Therefore it is less likely that the channel pore is located at the dimer interface.

SORCI/CHARMM
We have presented two models of the helix A and B arrangements in ChR and have discussed the expected properties on the bases of experimental data and QM/MM calculations. In both models, the key characteristics are the charged residues in helix B that define in conjunction with residues of helix C and G the ion conducting pore that goes all the way through from the extracellular phase to the intracellular bulk phase. This is the fundamental distinction from the pumping rhodopsins where the retinal is never simultaneously connected to the extra-and intracellular phases by a water column and hydrogen bond networks. A combination of modeling and electrical studies of mutants revealed residues that are important for conductance and ion selectivity. The light switch is conserved in ChR from other microbial rhodopsin, but details about early retinal isomerization and gating of the channel remain unknown and will be addressed in a more extended project already in progress. Strictly speaking, we are unable to prefer any of those two models, and we could even imagine that both arrangements may be realized in different ChRs or even in one ChR under variable ionic or voltage conditions.