Ion-pulling simulations provide insights into the mechanisms of channel opening of the skeletal muscle ryanodine receptor

The type 1 ryanodine receptor (RyR1) mediates Ca2+ release from the sarcoplasmic reticulum to initiate skeletal muscle contraction and is associated with muscle diseases, malignant hyperthermia, and central core disease. To better understand RyR1 channel function, we investigated the molecular mechanisms of channel gating and ion permeation. An adequate model of channel gating requires accurate, high-resolution models of both open and closed states of the channel. To this end, we generated an open-channel RyR1 model using molecular simulations to pull Ca2+ through the pore constriction site of a closed-channel RyR1 structure determined at 3.8-Å resolution. Importantly, we find that our open-channel model is consistent with the RyR1 and cardiac RyR (RyR2) open-channel structures reported while this paper was in preparation. Both our model and the published structures show similar rotation of the upper portion of the pore-lining S6 helix away from the 4-fold channel axis and twisting of Ile-4937 at the channel constriction site out of the channel pore. These motions result in a minimum open-channel pore radius of ∼3 Å formed by Gln-4933, rather than Ile-4937 in the closed-channel structure. We also present functional support for our model by mutations around the closed- and open-channel constriction sites (Gln-4933 and Ile-4937). Our results indicate that use of ion-pulling simulations produces a RyR1 open-channel model, which can provide insights into the mechanisms of channel opening complementing those from the structural data.

The type 1 ryanodine receptor (RyR1) mediates Ca 2؉ release from the sarcoplasmic reticulum to initiate skeletal muscle contraction and is associated with muscle diseases, malignant hyperthermia, and central core disease. To better understand RyR1 channel function, we investigated the molecular mechanisms of channel gating and ion permeation. The skeletal muscle ryanodine receptor (RyR1) 4 mediates muscle contraction by releasing Ca 2ϩ from the sarcoplasmic reticulum (1). Provided its essential role in muscle function, RyR1 mutations are associated with muscle diseases malignant hyperthermia and central core disease (2,3). RyR1 shares considerable homology with other Ca 2ϩ release channels including other ryanodine receptor isoforms and isoforms of the inositol 1,4,5-triphosphate receptor. All of these receptors release Ca 2ϩ from a Ca 2ϩ store, are regulated by Ca 2ϩ , and form a cation conducting pore (4 -9). It is likely that these receptors share an evolutionarily conserved channel gating mechanism (10). A structural understanding of RyR1 gating will provide valuable insights into mechanisms of Ca 2ϩ signaling (11).
RyR1 forms a homotetrameric complex of 2,200 kDa. Due to the large size of the receptor structural information is limited. Early cryo-electron microscopy (cryo-EM) studies at ϳ30-Å resolution revealed the global architecture of the receptor (12)(13)(14). The channel consists of a large cytoplasmic N-terminal domain that contains interaction sites for regulatory molecules, and a smaller C-terminal segment containing the pore-forming transmembrane domain. Although the N-terminal segment comprises the largest portion of RyR1, the isolated C-terminal segment (residues 3661 to 5037) can form a functional Ca 2ϩ release channel (15), and the two C-terminal transmembrane segments, including the pore helix and connecting loops, can form a homotetrameric assembly capable of conducting K ϩ , Cl Ϫ , and Ca 2ϩ (16).
More detailed structures of the closed-channel RyR1 were later determined at 10-Å resolution (17)(18)(19). Open-and closedchannel structures determined at 10-Å resolution provided valuable information on the overall motions involved in channel opening, but lacked atomic detail (17). Recently, the closed and open RyR1 and RyR2 channel structures were determined at near atomic resolution of 3.6 to 4.8 Å (20 -23). The pore-forming region consists of a long S6 inner helix of 50 residues, a pore helix of 15 residues, and a GGGIGDE selectivity filter. In the RyR1 closed channel, Ile-4937 lines the pore and forms the hydrophobic constriction site with a pore radius of less than 1 Å, rendering the channel impermeable to Ca 2ϩ (20,21). In the RyR2 open channel, Gln-4933 forms the constriction site with a minimum radius of 2.4 Å (23).
Computational methods offer a means to provide structural insights. Due to the large size of the RyR ion channels, it is not possible to reach the simulation time necessary to observe a channel opening event (microsecond to millisecond time scale). One means of addressing this limitation is to use steered molecular dynamics simulations, a technique that includes an additional external force to propel the system along a particular trajectory (24). Steered molecular dynamics simulations have previously been applied to explore channel opening for the voltage-gated potassium channel (25), the mechanosensitive chan-nel of large conductance (26), and the nicotinic acetylcholine receptor (27).
In this work we use a novel computational method to produce an open-channel model for the RyR1 transmembrane domain. We obtain the model using a form of steered molecular dynamics in which we apply a potential to pull a Ca 2ϩ ion through the pore constriction site. From ion pulling alone we obtain an open-channel model exhibiting a substantially larger pore and counter-clockwise twisting of the S6 helix. It was our expectation that these motions increasing the pore diameter would occur along the path of least resistance. We hypothesize that this channel opening pathway is the physiologically preferred opening pathway and that ion pulling simulations can produce an accurate open-channel model. Indeed, after refinement we find that our open-channel model is consistent with the high resolution cryo-EM structures of the open RyR1 and RyR2 channels (22,23). We also present mutational data sup-porting a mechanism of channel gating in which the S6 helix rotates Ile-4937 out of the pore. In our open-channel model, Gln-4933 forms the minimum radius of the pore. Moreover, by modeling channel gating we propose a role for the conserved glycine residue Gly-4934 in channel gating. The data suggest that increasing the side chain volume of the conserved glycine residue Gly-4934 impedes movement of Ile-4937 away from the pore, which provides a plausible explanation for decreased channel gating and Ca 2ϩ conductance.

Creating the open-channel RyR1
To create the open-channel RyR1, we place a Ca 2ϩ ion at the cytoplasmic end of the channel between residues Glu-4948 and Glu-4952 ( Fig. 1, a and b). From this starting position the Ca 2ϩ ion is pulled into the channel constriction site formed by Ile- Top panels (a, c, e) show the side view of the pore, whereas lower panels (b, d, f) show the pore viewed from the cytoplasmic end of the channel. Side chains for amino acids bordering Ca 2ϩ are shown as sticks and colored according to residue type, red, acidic; green, polar; and white, aliphatic. In the lower panels (b, d, f), Gln-4933 and Ile-4937 are shown in green and white, respectively. The position of Ca 2ϩ is shown as a magenta sphere and waters within 4 Å of Ca 2ϩ are shown in red spheres. The initial structure represents the time point in which Ca 2ϩ was first moved to the position between Glu-4948 and Glu-4952, before water was allowed to equilibrate around the ion. 4937, which is completely dehydrated in the closed-channel structure. Using a short 20-ns simulation, we equilibrate this system with the Ca 2ϩ ion constrained in the pore to generate our intermediate open conformation (Fig. 1, c and d). At this point position restraints have only been applied with regard to the location of the Ca 2ϩ ion and no additional restraints have been applied to the protein.
This simulation results in an open-channel pore with a minimum pore radius of 2.8 Ϯ 0.4 Å. The pore constriction site of this channel is formed by Gln-4933 compared with Ile-4937, which formed the constriction site in the closed-channel RyR1 (21). In this intermediate conformation the Ca 2ϩ ion in the pore is surrounded by 7 to 8 water molecules over the 20-ns simulation, consistent with Ca 2ϩ keeping its first hydration shell as it passes through the channel.
However, the process of pulling the Ca 2ϩ ion into the pore resulted in distortion of the S6 helix. Given the energetic cost of breaking a number of hydrogen bonds in the S6 helix, we presumed the distortion of the S6 helix is an artifact of the ion pulling simulation. Indeed, much of the distortion occurred as the Ca 2ϩ ion moved past the negatively charged residues lining the closed channel pore. In the final open-channel these negatively charged residues move away from the channel pore and cannot simultaneously bind Ca 2ϩ . We therefore refined the model by applying a series of restraints to the protein. We applied hydrogen bond restraints to the S6 helix to reform the distortions from Ca 2ϩ pulling. We also applied C␣ inter-helical restraints for residues Ile-4937, Asp-4945, and Glu-4948 to preserve the open channel conformation and encourage the formation of a symmetric channel. Finally, we note that the initial ion pulling caused Ile-4937 to rotate out of the pore toward the hydrophobic amino acids Ile-4931 and Leu-4935. To further encourage interaction of Ile-4937 with Ile-4931 and Leu-4935, we applied restraints between C␤ atoms of Ile-4937 and Ile-4931 or Leu-4935. After simulation with these imposed restraints we obtain a refined open-channel RyR1 model (Fig. 1, e and f).

Open pore conformation of RyR1
After our refinement simulations we obtain an open-channel model of RyR1 with a minimum channel radius of 2.8 Ϯ 0.3 Å. This radius is moderately smaller than the experimentally derived minimum radius of 3.3 to 3.5 Å (28). However, the radius is consistent with those of the open-channel RyR1 The structure previously presented had the side chain of Gln-4933 rotated out of the pore. However, the other two structures had conformations of the Gln-4933 side chain rotated into the pore, causing the average pore profile to reflect the constriction at Gln-4933. The results not only support our model, but also suggest a significant amount of flexibility in the Gln-4933 side chain to allow the passage of Ca 2ϩ .
We assess the accuracy of our RyR1 open-channel model by comparing the alignment of the structural model to the 3.8-Å resolution cryo-EM densities of the open-channel RyR1 (22) (Fig. 3, a and b). We note that the most significant difference between the open-and closed-channel structures occurs at the pore-lining S6 helix, particularly at the cytoplasmic end of the channel. Comparing the alignment of the S6 helix (residues 4910 to 4957) of our open-channel model to the open-channel densities we observe a correlation of 0.93 Ϯ 0.02 (Table 1). In contrast, the alignment of the closed-channel structure (PDB code 3J8H) with the open-channel densities results in a correlation of 0.84 Ϯ 0.02 (Fig. 3, c and d,

Conformational changes of RyR1 channel opening Conformational changes opening the channel
To analyze the mechanism of RyR1 channel gating, we compare global features of the S6 helix to the RyR1 closedchannel structure (PDB code 3J8H) (21). We perform principal component analysis of the closed-channel structure and the 10 open-channel models used to calculate the pore profile to identify the major modes of motion leading to channel opening. The first principal component explains 70% of the motion between closed and open RyR1. The major motion of the first principal component leading to channel opening is an outward rotation of the S6 helix away from the 4-fold channel axis, particularly at the cytosolic end of the channel ( Fig. 4 and supplemental Video S1). We quantify this change by measuring radial and lateral tilting angles for the upper portion of the S6 helix comprised of residues 4930 to 4950. The radial tilting angle is the angle of helical tilting away from the 4-fold channel axis, whereas the lateral angle is the tilting angle between the helix and channel axis in the direction perpendicular to the plane defined by the radial angle (Fig. 4a). The outward radial tilting for the upper portion of the S6 helix is 7 Ϯ 1°, which is in excellent agreement with values of 7 Ϯ 1°and 7°observed for open-channel RyR1 and RyR2 structures, respectively ( Fig. 4b) (22,23). In contrast, we observe less movement for the luminal end of the S6 helix. The change in radial tilting for the lower portion of the S6 helix (residues 4910 to 4930) is 4 Ϯ 3°compared with values of 1 Ϯ 1°observed in the experimental structures. We also observe a relatively small amount of lateral tilting averaging 3 Ϯ 1°tangent to the pore compared with a lateral  We also examine the twisting motion of the S6 helix accompanying channel opening. We observe that backbone atoms for the upper portion of the S6 helix undergo a 13 Ϯ 1°clockwise rotation viewed from the cytosolic side with respect to the closed-channel RyR1 (Fig. 4a) Tilting angles of the open model after the initial ion pulling simulation were 5 Ϯ 10°and 5 Ϯ 1°for radial and lateral tilting angles, respectively. Moreover, we observe a 13 Ϯ 4°clockwise rotation (as viewed from the cytoplasm) after the initial ion pulling simulation, before helical and inter-subunit restraints were introduced. We find that the restraints introduced during refinement effectively reduced the variance observed for rotation and radial tilting of the S6 helix without introducing much new motion. To provide independent evidence for the proposed open model, we generated and determined the functional properties of 7 mutations at positions 4933, 4934, and 4937 ( Table 2). The mutants were expressed in HEK293 cells, and their expression levels and retention of function were determined by immunoblot analysis, caffeine-induced cellular Ca 2ϩ release assay, and in vitro binding assay using the highly RyR-specific ligand [ 3 H]ryanodine. Ryanodine is known to only bind to the openchannel RyR and is therefore commonly used to assess channel function in isolated membrane fractions. Single channel recordings were performed to determine Ca 2ϩ -dependent gating and ion permeation properties of the mutants.

Altering open-channel contacts in the S6 helix affects channel gating
Mutation of residue Ile-4937 to alanine does not change the number of cells that show a caffeine-induced Ca 2ϩ release response compared with wild-type RyR1 (Fig. 6). However, following isolation from HEK293 cells we observe reduced [ 3 H]ryanodine binding, lack of Ca 2ϩ regulation, and lack of Ca 2ϩ conductance ( Table 2, Fig. 7), which suggests that the Ile-4937 to alanine mutation results in a metastable channel in HEK293 cells. As the I4937A mutation only reduces the side chain size, but does not change the polarity, any loss of function resulting from alanine substitution should result from removing nonpolar contacts. Increasing the side chain volume at position 4937 by mutating to valine restores channel function. The I4937V mutation restores [ 3 H]ryanodine binding, Ca 2ϩ regulation, and Ca 2ϩ conductance in single channel experiments. The permeability ratio is also restored to a value that is not appreciably different from wild-type RyR1. Mutation of Ile-4937 to threonine also results in functional channels, however, with Ca 2ϩ currents and permeability ratios that are significantly (p Ͻ 0.05) decreased compared with RyR1-I4937V. As the threonine side chain volume (122 Å 3 ) is similar to the side chain volume of valine (139 Å 3 ) (29), we expect that the reason for the decreased Ca 2ϩ selectivity observed for RyR1-I4937T is the result of disrupting non-polar interactions that alter the openchannel structure. The arrow lengths are scaled according to magnitude of the motion. Likewise, the structure is colored by the magnitude of the motion from smallest (blue) to largest (red). Black arrows in a denote the planes within which radial and lateral helical tilting angles were determined for the helix the arrows are centered on. The black arrows in b highlight the region over which the helical tilting angles were computed. Mutation of Gln-4933 to alanine resulted in caffeine-induced Ca 2ϩ release in HEK293 cells (Fig. 6) but loss of ryanodine binding, increased mean close times, and decreased Ca 2ϩ currents and permeability ratios (Table 2, Fig. 7). Comparable results were obtained with the corresponding mutation in the cardiac ryanodine receptor (RyR2-Q4863A), which showed a caffeine response in HEK293 cells and the presence of functional channels in single channel measurements but a pronounced decrease in [ 3 H]ryanodine binding affinity (30). Mutation of Gln-4933 to valine had an even greater effect resulting in channels with significantly decreased caffeine response, Ca 2ϩ currents, and permeability ratios compared with either wild-type RyR1 or RyR1-Q4933A. Large variability in mean open and close times suggested that channels with a different gating behavior were incorporated in lipid bilayers.
We previously investigated the single channel properties of purified RyR1-G4934A and RyR1-G4934V mutants (31). We here examine the single channel properties of these mutants in membrane fractions avoiding exposure to detergent during purification of RyR1s. Similar to our previous findings the G4934A mutation resulted in functional channels with altered gating behavior, reduced Ca 2ϩ current, and permeability ratio ( Table 2, Fig. 7). However, the G4934V mutation also resulted in functional channels with reduced Ca 2ϩ current and permeability ratios when measured in membrane fractions. The number of channel events decreased and mean close times increased. The result suggests that the G4934V mutation possibly destabilizes the channel so it cannot retain function throughout the process of purification with detergent.

Discussion
We report an open-channel model derived by ion pulling simulations. As we did not include any EM data from an openchannel RyR (22,23) in our modeling, the significant agreement between our model and the open-channel cryo-EM density validates the accuracy of both our model and approach. The rationale for our approach was that moving a Ca 2ϩ ion through the channel constriction site would force structural rearrangements creating a pore capable of allowing Ca 2ϩ to pass through. We assumed that the RyR1 structure would open via the path of least resistance and that this pathway would also be preferred physiologically. The agreement between our model and the experimental structures also supports our hypothesis that channel opening occurs via the most energetically favorable pathway.
Our model accurately reproduced the ϳ12°clockwise rotation of the backbone atoms in the S6 helix viewed from the cytoplasmic end of the channel. The twisting motion rotates C␣ atoms of Ile-4937 and Gln-4933 away from the pore while arranging C␣ atoms of Gly-4934 to face the pore. Rotation of Ile-4937 out of the pore creates a much larger pore radius than could be achieved by widening the inter-helical distances by tilting the S6 helix. Helical rotation has been observed as a pore opening mechanism in other channels including the potassium channel from Streptomyces lividans (KcsA) (32), the NaK channel (33), and the nicotinic acetylcholine receptor (nAChR) (34) suggesting that twisting may be a common pore opening mechanism among ion channels.
Our functional measurements also support rotation of Ile-4937 out of the pore. In the closed-channel structure Ile-4937 lines the pore and forms the hydrophobic constriction site. If Ile-4937 faced the pore in the open channel conformation, we would expect that reducing the side chain volume (I4937A) or making the residue more hydrophilic (I4937T) would make the channel more open. However, we find that the I4937A and I4937T mutations abolish and reduce calcium current, respectively. Moreover, the I4937A mutation markedly reduces ryanodine binding, which suggests that this mutation removes contacts necessary for maintaining a stable channel when removed from its cellular environment. Although the side chain volumes of threonine and valine are similar, we find that the function of RyR1-I4937V is more similar to wild-type RyR1 than RyR1-I4937T. We suggest that it is the polar nature of the I4937T mutation that causes the reduced function of the RyR1-I4937T mutant. Taken together, the results suggest that in the openchannel structure Ile-4937 no longer faces the pore but rather interacts with hydrophobic residues, likely in neighboring S6 helices.
In the open-channel model, Gln-4933 faces the pore. The molecular volumes for alanine (90 Å 3 ) and valine (139 Å 3 ) are, respectively, smaller and similar as compared with the molecular volume of glutamine (151 Å 3 ) (29). On the basis of molecular volume we would expect a more open channel pore, which seems to contradict the experimental finding. However, both the Q4933A and Q4933V mutations also decrease the polarity of the side chain at this site. We compute the pore profiles and apply a Poisson-Boltzmann electrostatic model (35) to determine the Ca 2ϩ solvation energy for wild-type and mutant RyR1 channels (Fig. 8). Our model shows that the Q4933V and  (7) 793 Ϯ 21 (7) 742 Ϯ 14(7) I Ca (pA) (10 mM trans Ca 2ϩ ) Ϫ2 Q4933A mutations increase the pore radius at position 4933 from 2.8 Ϯ 0.3 to 2.9 Ϯ 0.1 Å and 4.0 Ϯ 0.1 Å, respectively. Both mutations also significantly increase the Ca 2ϩ solvation energy, suggesting that the polarity of the residue at the 4933 position is important for Ca 2ϩ conductance. However, we cannot rule out the importance of an interaction between Gln-4933 and Val-4926 to maintain an open-channel conformation (Fig. 5).
Although the radial tilting angle and the twisting angle of the S6 helix are accurately reproduced in our model, we note a difference in the lateral tilting angle for the S6 helix (   In our open-channel model, the direction of rotation brings Ile-4937 past Gly-4934. We previously showed that the G4934A mutation reduced the Ca 2ϩ current, whereas the G4934V mutation resulted in loss of function of purified channels (31). In the current study, we find that RyR1-G4934V can form functional channels with reduced Ca 2ϩ current, when avoiding detergent exposure throughout the process of purification. The result affirms our previous prediction that the G4934V muta-  Table 1. tion reduces channel stability. Previously, we explained the effects of these mutations in terms of increased interaction with residues in the adjacent S6 helix (Ile-4936, Ile-4937, and Phe-4940) within the closed-channel structure. Our model and the structures of the open-channel RyR1 suggest that the G4934A mutation could also affect channel function by impeding the rotation of Ile-4937 to produce a fully open channel. This hypothesis is supported by analysis of the open-and closedchannel times. The G4934A mutation primarily lengthens mean closed times (38.7 Ϯ 9.5 versus 11.0 Ϯ 3.5 ms for WT, n ϭ 8 -11) rather than affect mean open times (0.64 Ϯ 0.08 versus 0.71 Ϯ 0.14 for WT, n ϭ 8 -11), which suggests that the mutation decreases the probability of closed to open transitions. A decrease in closed to open transitions for the G4934A mutation would be consistent with our hypothesis that rotation of Ile-4937 toward Gly-4934 in the adjacent subunit is required for channel opening.

Conclusion
Our modeling of the open-channel RyR1 reveals two primary motions responsible for creating an open-channel pore. The first is a radial tilting of the S6 helix away from the pore along the 4-fold channel axis, which results in a general widening of the pore at the cytoplasmic end of the channel. The second major motion is a twisting of the S6 helix, which rotates the hydrophobic residue Ile-4937 out of the pore. Our results indicate that this rotation brings Ile-4937 past Gly-4934 and into contact with other hydrophobic residues in the S6 helix, which may be responsible for stabilizing the open-channel conformation. Our model also indicates that the constriction site for the open-channel RyR1 is Gln-4933 and our single channel data suggest that the polarity of this residue plays an important role in Ca 2ϩ conductance. Altogether, our modeling reveals both the major motions of the RyR1 transmembrane domain responsible for channel gating and plausible roles of RyR1-Gln-4933, RyR1-Gly-4934, and RyR1-Ile-4937 in channel gating.

Materials
[ 3 H]Ryanodine was obtained from PerkinElmer Life Sciences, protease and phosphatase inhibitors from Sigma, and phospholipids from Avanti Polar Lipids.

Molecular dynamics simulations and pore opening of the RyR1 closed-channel
The transmembrane domain of the 3.8-Å resolution cryo-EM structure (PDB code 3J8H) (21) was inserted into a lipid bilayer with a lipid composition comparable with the one used in single channel measurements (5:3:2 molar ratio of 1palmitoyl-2-oleoylphosphatidylethanolamine, 1-palmitoyl-2oleoylphosphatidylserine, and 1-palmitoyl-2-oleoylphosphatidylcholine) using the GROMACS g_membed tool. The bilayer was generated using the CHARMM-GUI membrane builder (36). The RyR1 transmembrane segment modeled included residues 4545 to 5033 based on the sequence present in the cryo-EM structure. The 33-amino acid loop (residues 4588 to 4621) missing from the reported structure (21) was replaced with 3 glycine residues. Each system was solvated using the TIP3P water model and ionized with 0.02 M CaCl 2 and 0.27 M KCl. All simulations were performed using GROMACS 4.6.1 (37) with the CHARMM36 force field (38,39). The simulation was performed at a constant pressure and temperature of 1 atm and 298 K. All bonds were constrained using the LINCS algorithm (40). The integration time step used for all simulations, unless otherwise noted, was 2 fs. Particle mesh Ewald was used for long-range electrostatic interactions. A 10-Å cutoff was used for non-bonded interactions.
To equilibrate the system, energy minimization was performed for 10,000 steps with harmonic position restraints of 10 kJ/(mol Å 2 ) on the protein backbone atoms followed by 3 ns of equilibration, during which position restraints on protein backbone atoms were gradually reduced from 10 to 0 kJ/(mol Å 2 ). Channel opening was performed using the molecular pulling code implemented in GROMACS. A Ca 2ϩ ion was inserted at the cytosolic side of the channel and pulled through the channel gate at a rate of 0.1 Å/ns over 100 ns of simulation using a force constant of 10 kJ/(mol Å 2 ). The Ca 2ϩ ion was fixed to the channel constriction site for an additional 20 ns of simulation with a harmonic potential of 10,000 kJ/(mol Å 2 ) along the channel axis. The ion was free to move within the membrane plane.
To further refine the open-channel inter-subunit distance, restraints were added between C␤ atoms of Ile-4937 and Ile-4931 in the previous subunit and between C␤ atoms Ile-4937 and Leu-4935 in the previous subunit. The inter-subunit distance restraints used were 7.8 and 7.0 Å for residues 4931 and 4935, respectively. The restraints were designed to further rotate Ile-4937 out of the pore in the direction observed after the initial pulling. As the additional torque of twisting was likely to distort the helices, helical distance restraints were added between backbone nitrogen and oxygen atoms for residues 4930 to 4953. The upper boundaries used for the distance restraints were obtained from distances observed in the cryo-EM structure plus 0.5 Å. Lower boundaries were all 2.7 Å. Inter-subunit distance restraints of 12, 15, and 16 Å were added between C␣ atoms of adjacent subunits for residues Ile-4937, Asp-4945, and Glu-4948, respectively. The listed distances were obtained from the mean distance among four chains observed after the initial ion pulling simulation. The restraints were added gradually over the course of 14 ns of simulation using a time step of 1 fs. Restraint potentials were all 1000 kJ/(mol Å 2 ). A final production run of 20 ns was performed that was used for data analyses.

Alignment with EM density
Initial alignment of structures with cryo-EM maps was performed using the colors function in Situs (41). The last structural snapshot from the pore opening calculations was aligned to the cryo-EM map using a Gaussian density map generated at 4-Å resolution. Rigid body fitting was performed sampling the rotational search space with a 5°step size. The resulting fit was refined using the fit-in-map function implemented in the Chimera software package. The Chimera fit-in-map function was also used to evaluate the correlation between the structural models and the experimental EM data. For all structures we fit the map of the backbone and C␤ atoms in the structure generated at 4-Å resolution for residue ranges 4820 to 4957 and 4910 to 4957. Residues in loops (residues 4861 to 4879 and residues 4891 to 4909) were excluded from the fitting. The atoms included and the resolution for fitting were chosen to maximize the correlation of the experimentally derived structure with the experimental data. The correlations of the open structure (PDB code 5TAL) to the open densities (EMD-8376, EMD-8377, and EMD-8378) and the closed structure (PDB code 3J8H) to the closed density (EMD-2807)(23) were 0.98 and 0.97, respectively.

Structural analyses
For all structural analyses 10 structures sampled every 2 ns from the last 20 ns of simulation were used. Principal component analysis was performed on the 10 open-channel structures from simulation combined with the closed-channel cryo-EM structure (21) using ProDy (42). The eigenvectors of the first principal component were visualized using the NMWiz plugin in VMD (43). Radial and lateral tilting angles were calculated for both upper (4930 to 4950) and lower (4910 to 4930) portions of the S6 helix according to the previous method (44). The mean Ϯ S.E. are reported in the results. The twist angle was measured as the average angle between the respective backbone atoms of the open-channel model and closed-channel cryo-EM structure within the membrane plane. The angle was calculated after centering both over the residue range 4930 to 4950. The reported value is the mean Ϯ S.E. from 10 structures and four chains.

Pore profile and Poisson-Boltzmann calculations
All pore profiles were calculated using the HOLE program (45). To generate pore profiles for RyR1-WT, RyR1-Q4934A, and RyR1-Q4934V, we extracted five structures sampled evenly over the last 20 ns of the pore opening simulation. Amino acid substitution and rotamer optimization were performed using Eris (46,47). For each structure the Eris program was used to perform 200 independent Monte Carlo optimizations allowing side chain repacking within 10 Å of the mutation site for both wild-type and mutant proteins. From each set of calculations we selected the minimum energy structure for both the mutant and wild-type protein (five structures each) and calculated the pore profile for each of these structures.
The same structures were used in the Poisson-Boltzmann calculations performed using APBSmem (48). The charge model used for calculations was the Swanson charge model (49), which was applied using the PDB2PQR software (50). Two focusing layers were used. The finest grid included 97 grid points in the x, y, and z dimension and had a grid length of 90 Å. Ion concentrations were set to 0.1 M. Non-linear Poisson-Boltzmann method was implemented in the Adaptive Poisson-Boltzmann Solver (35). Upper and lower exclusions were set at 16 Å to exclude the membrane from the pore region. Membrane thickness was 42.5 Å with a head group thickness of 7 Å similar to values previously used to calculate the Ca 2ϩ solvation energy profile of the TRPV1 receptor (51). Dielectrics for membrane and protein were set to 2.0 and dielectrics for solvent and head groups were set to 80.0. The ion step size was 1.0 Å along the channel axis and the ionic radius used for Ca 2ϩ was 1.03. The temperature was set at 298 K.

Expression and preparation of wild-type and mutant channels
pCMV5-RyR1 mutants were prepared as described previously (52) or by a gene synthesis method using a proprietary protocol (Genewiz Inc., South Plainfield, NJ). Wild-type and mutant expression vectors were transiently expressed in HEK293 cells with jetPRIME (Polyplus, New York), according to the manufacturer's instructions. Cells were maintained in Dulbecco's minimum essential medium (DMEM/F-12 (1:1)) containing 10% fetal bovine serum and 1ϫ Antibiotic Antimycotic Solution at 37°C and 5% CO 2 and plated the day before transfection. For each 10-cm tissue culture dish, 10 g of DNA was used at a DNA/jetPRIME ratio of 1:2. Following transfection, cells were kept at 35°C and harvested 66 -72 h after transfection. Keeping transfected cells at 35°C instead at 37°C previously used (53) increased the number of functional RyR1s, as determined using [ 3 H]ryanodine ligand binding method. Membrane fractions were prepared in the presence of protease inhibitors and 1 mM oxidized glutathione (GSSG) as described (53).

Cellular Ca 2؉ release
Retention of the function in HEK293 cells was studied by cellular Ca 2ϩ release assay using the Ca 2ϩ releasing drug caffeine described in the previous publication (52). Briefly, HEK293 cells grown on glass coverslips were incubated with 5 M Fluo 4-AM. Cellular Ca 2ϩ release was induced by the addition of ϳ8 mM caffeine and measured in individual cells using EasyRatioPro (Photon Technology International, Lawrenceville, NJ).

SDS-PAGE and immunoblot analysis
Proteins in membrane preparations were separated using 3-12% gradient SDS-PAGE, transferred to nitrocellulose mem-branes, probed with primary rabbit anti-RyR1 polyclonal antibody 6425, and quantified using Image Lab 4.1 Analysis Software (Bio-Rad) (31).

[ 3 H]Ryanodine binding
Expression of functional RyR1 mutants in membrane preparations was assessed using [ 3 H]ryanodine-binding assay. Ryanodine binds with high specificity to RyRs and is widely used to probe for RyR activity and content (54). Expression levels were determined by measuring maximum binding capacity (B max ) of [ 3 H]ryanodine binding of membrane preparations for 4 -5 h at 24°C with a near-saturating concentration of 20 nM [ 3 H]ryanodine in 20 mM imidazole, pH 7.0, 0.6 M KCl, protease inhibitors, 1 mM GSSG, and 0.1 mM Ca 2ϩ . Nonspecific binding was determined with a 1000-fold excess of unlabeled ryanodine. Bound [ 3 H]ryanodine binding was determined using a filter assay (52).

Single channel recordings
Single channel measurements were performed by the fusion of membrane preparations with Mueller-Rudin type planar lipid bilayers containing a 5:3:2 mixture of bovine brain phosphatidylethanolamine, phosphatidylserine, and phosphatidylcholine (25 mg of total phospholipid/ml of n-decane) (52). Single channels were recorded in symmetric 0.25 M KCl, 20 mM K-HEPES, pH 7.4, buffers with additions as indicated in the text. The trans side of the bilayer was defined as ground. Electrical signals were filtered at 2 kHz (0.5 kHz for Ca 2ϩ currents at 0 mV), digitized at 10 kHz, and analyzed at 50% threshold setting (52). Data acquisition and analysis were performed with commercially available software package (pClamp, Axon Instruments, CA), using 2-min recordings for analysis. To determine permeability ratios, single channel activities were recorded in symmetrical 250 mM KCl solution with 10 mM Ca 2ϩ on the trans side and the reversal potential (E rev ) was measured. The permeability ratio of Ca 2ϩ versus K ϩ (P Ca /P K ) was calculated using a modified form of the Goldman-Hodgkin-Katz equation.

Biochemical assays and data analysis
Free Ca 2ϩ concentrations were established by adding varying amounts of Ca 2ϩ to 1 mM EGTA solutions and determined with the use of a Ca 2ϩ selective electrode. Differences between samples were analyzed using Student's t test. p Ͻ 0.05 was considered significant.