The Ca²⁺ permeation mechanism of the ryanodine receptor revealed by a multi-site ion model

As an essential messenger in cells, calcium ions (Ca²⁺) regulate many physiological processes, including neurotransmitter release, muscle contraction, and hormones secretion¹. The concentration of Ca²⁺ in the cytoplasm and organelles is precisely controlled by multiple calcium channels, including the voltage-gated calcium channels in cell membranes and the ryanodine receptors (RyR) in the endoplasmic reticulum (ER) membrane. Ca²⁺ also induces conformational changes of a wide range of Ca²⁺-interacting proteins, such as calmodulin and Ca²⁺-activated ion channels, to trigger downstream signal transduction². Although the pathways of calcium signaling are extensively studied, the molecular interaction details between calcium and proteins have yet to be fully elucidated.

To study the detailed interactions between ions and proteins, we can often use molecular dynamics (MD) simulations to provide microscopic and quantitative insights, thereby obtaining the specific functional mechanisms of the relevant proteins^3,4,5,6. However, the conventional models of Ca²⁺ are far from accurate in calculating the interaction energies between Ca²⁺ and proteins^7,8,9, therefore inadequate to study the precise Ca²⁺-protein interactions. As a consequence, K⁺ and Na⁺ channels have been widely studied by MD simulations, and their detailed permeation mechanisms were revealed^{10,11,12,13,14,15,16,17}, but computational studies of Ca²⁺ channels are rather limited. Notably, several structures of Ca²⁺ channels were resolved recently^18,19,20,21, which provided a solid basis to study their detailed function mechanism further. In particular, the open-state RyR1 channel provides an excellent opportunity for studying Ca²⁺ permeation and selectivity²¹, which makes a reliable Ca²⁺ model even more desirable.

There have been enormous efforts in trying to develop a more accurate Ca²⁺ model. The polarizable force field is theoretically appealing⁹, but its implementation and validation still need further work before being widely accepted and utilized for membrane protein simulations. Kohagen et al. proposed to scale the partial charges on the Ca²⁺ to account for the charge transfer and polarization^22,23. Another strategy is to represent an ion by distributing electrostatic and Lennard-Jones (LJ) interactions on multiple sites, which can introduce a much larger parameter space and therefore make the ion model more tailorable^24,25,26. Unfortunately, none of the existing Ca²⁺ models in the non-polarizable classical force field are able to describe the interactions between Ca²⁺ and protein quantitatively. Therefore, in the present work, we developed a new multi-site Ca²⁺ model particularly optimized for Ca²⁺-protein interaction (inset of Fig. 1), and then utilize this model to study the detailed Ca²⁺ permeation mechanism through the RyR1 channel, which showed distinct features from the widely studied K⁺ and Na⁺ channels.

The labels on the x-axis represent the PDB IDs of the ten Ca²⁺-bound proteins for which we performed the binding energy calculations. The legend indicates the methods of the calculation, QM+10: corrected quantum mechanics calculation; CAM: our multi-site Ca²⁺ model; Drude: the polarizable Drude force field; C36: the CHARMM36 force field. Our multi-site Ca²⁺ model consists of a central atom (C) and six dummy atoms (D) located at the vertices of an octahedron, as shown with the inset. Source data are provided as a Source Data file.

Full size image

The multi-site Ca²⁺ model behaves well in water and protein

We designed a seven-site ion model, as shown in the inset of Fig. 1. There are six adjustable parameters, including b_CD, Q_C, ε_C, σ_C, ε_D, and σ_D, where b_CD is the distance between dummy atoms and the central atom, Q_C is the charge on the central atom, and ε and σ are the LJ parameters of the central (C) and dummy (D) atoms, respectively. We further distinguish the LJ interactions of Ca²⁺ with water and non-water by replacing (ε_C, σ_C) with two sets of (\(\varepsilon _{\mathrm{C}}^{\mathrm{W}}\), \(\sigma _{\mathrm{C}}^{\mathrm{W}}\)) and (\(\varepsilon _{\mathrm{C}}^{{\mathrm{NW}}}\), \(\sigma _{\mathrm{C}}^{{\mathrm{NW}}}\)), respectively. The charges on all the dummy atoms are the same and determined so that the total charge of the model is +2e. By adjusting the aforementioned parameters (Supplementary Fig. 1), we obtained a Ca²⁺ model that can quantitatively reproduce the energetical and dynamic properties of Ca²⁺ in water as well as the Ca²⁺-protein interactions (please refer to the method section and Supplementary Tables 1 and 2 for details). The resulting Ca²⁺ properties in water are shown in Table 1. As can be seen, the hydration free energy (ΔG_h), the first-peak position of the radial distribution function of water around Ca²⁺ (R₁), and the number of coordinated water molecules in the first solvent shell (N_C) have all reached the target values of experiments. In addition, the residence time of water molecules in the first solvation shell (τ_R) can be optimized to below 100 ps, which was consistent with the previous experimental estimation^27,28, and solved the common problem of Ca²⁺ being too sticky. We also performed ab initio MD simulations, and observed multiple exchange of water molecules between the first and second solvation shells within 200 ps, providing additional support that the residence time of our model is reasonable (Supplementary Fig. 2). During the model optimization procedure, it was interesting to see that the water residence time in the first solvation shell can be optimized to a considerable extent with the multi-site model, as models with different parameters can show very close hydration energies while residence time differing by two orders of magnitude.

With our model, the calculated binding energies of Ca²⁺ and proteins were also improved to a large extent (Fig. 1). The default Ca²⁺ parameters of the CHARMM force field (C36) led to a significant overestimation of 150–200 kcal/mol⁹, while the average binding-energy discrepancies for ten selected proteins was −0.2 kcal/mol for our multi-site model. Therefore, our model is comparable to the quantum mechanics (QM) and polarizable Drude model in calculating the Ca²⁺-protein binding energies (Fig. 1 and Supplementary Table 1). While in the meantime, our multi-site Ca²⁺ model is compatible with the classical non-polarizable force field, allowing immediate and efficient simulations of complex Ca²⁺ interacting proteins.

We then further validated the thermodynamic property of our Ca²⁺ model around a protein by calculating the potential of mean force (PMF) of pulling a Ca²⁺ away from its binding site in a calpain. As can be seen in Supplementary Fig. 3, our model predicted a ΔG_PMF of ~33 kJ/mol and binding free energy of ~25 kJ/mol, while the default Ca²⁺ parameters of the CHARMM force field resulted in a much larger ΔG_PMF of ~177 kJ/mol and binding free energy of ~168 kJ/mol. The experimental binding affinity is ~35 kJ/mol²⁹, indicating that our Ca²⁺ model is more accurate in reproducing not only the potential energies but also the free energies in the study of Ca²⁺-protein interactions.

The permeability of the open-state RyR1

The conventional ion models generally work well in studying ion-protein interactions for K⁺ or Na⁺, but fail consistently for Ca²⁺ ions⁹, and therefore no calcium permeation was observed in previous MD studies of Ca²⁺ channels³⁰. We performed MD simulations on the type-1 ryanodine receptor (RyR1), an intracellular calcium release channel required for skeletal muscle contraction, with our Ca²⁺ model. The open-state structure of RyR1 (PDB ID: 5TAL [https://doi.org/10.2210/pdb5TAL/pdb]) was obtained from des Georges et al.’s²¹ work. The simulation systems consist of the pore domain of RyR1 embedded in a lipid bilayer of 1-Palmitoyl-2-oleoyl-sn-glycerol-3-phosphocholine (POPC) and an aqueous solution of either 150 mM Ca²⁺ or 250 mM K⁺ (Fig. 2a). The protein was restrained to the open-state crystal structure, and a transmembrane potential of 100 mV was applied along the direction from the sarcoplasmic reticulum (SR) lumen to the cytosol.

a The simulation system of Ca²⁺ (gray spheres) permeation through the RyR1 channel. The pore domain of RyR1 (cartoon) is embedded in a POPC membrane (gray surface). Cl⁻ ions are shown with blue spheres. Only part of water molecules (small red and white spheres) are shown for clarity. b The cumulative number of K⁺ ions permeating through the RyR1 channel as a function of simulation time. The dashed gray lines correspond to three independent simulation trajectories, and the solid red line corresponds to the average of the total permeation count. The conductance calculated from these trajectories is 910 ± 39 pS. c Same as b but for Ca²⁺ ions. The conductance from six trajectories is 141 ± 30 pS with our Ca²⁺ model. Source data of b and c are provided as a Source Data file.

Full size image

We first studied the permeation of K⁺ through RyR1 as a test as well as a validation. Three independent 300-ns MD simulations were conducted and generated enough permeation events (~500) for statistical analysis (Fig. 2b). The conductance was calculated to be 910 ± 39 pS, which was close to the experimentally measured saturated K⁺ conductance of ~1 nS and indicates that the RyR1 structure under study is indeed in its open state and that the K⁺ model in the CHARMM force field is reasonably accurate in describing the interactions between ions and proteins. However, our MD simulations of Ca²⁺ permeation through the same open-state RyR1 showed that the channel is not permeable to Ca²⁺ at all with the default Ca²⁺ parameter of CHARMM. Not a single permeation event was observed in three 500-ns MD trajectories (Fig. 2c). Since the channel is in its open state, this indicated that the default Ca²⁺ model gave us qualitatively wrong simulation results here, as also observed by another recent study³⁰. A close inspection of the MD trajectories showed that the Ca²⁺ ions were tightly bound to the protein, again confirming that the binding affinity between Ca²⁺ and protein is too strong with the conventional Ca²⁺ model. In contrast, when our Ca²⁺ model was used for the same simulations, we observed continuous Ca²⁺ permeation as expected (Fig. 2c). The conductance calculated from six 500-ns trajectories is 141 ± 30 pS, which agrees reasonably well with the experimental value of ~172 pS with the same ion concentration³¹. Therefore, we believe that our multi-site Ca²⁺ model is indeed more accurate in studying the permeation behavior of the Ca²⁺ channel.

The Ca²⁺ binding sites in the pore region

We performed an 800-ns simulation without the transmembrane potential to identify the Ca²⁺ binding sites in RyR1. The contour plot of the Ca²⁺ density, ρ(R, z), is presented in Fig. 3a, from which two major binding sites in the luminal vestibule (L) and above the selectivity filter constriction (S), and one minor binding site at the gate constriction (G) can be identified within the transmembrane pore. The selectivity filter and gate constrictions are located near the residues G4894 and Q4933, respectively (Fig. 3b). The corresponding positions in the contour plot are indicated by solid lines labeled with SF and GT in Fig. 3a. The isosurfaces of probability density corresponding to these binding sites are shown in Fig. 3b. By calculating the residence time of carboxylate oxygen of negatively charged residues within a sphere with a radius of 5 Å around the binding site L, we identified that the binding site L is formed by the interaction of Ca²⁺ ions with D4899, E4900, and D4903 residues (Fig. 3b), which agrees well with previous experimental and computational studies^30,32. The large probability of finding Ca²⁺ in the selectivity filter also indicates that the Ca²⁺ can easily accumulate around the luminal vestibule and move into the upper filter, and the rate-limiting step of permeation is the process of passing through the lower selectivity filter or gate constrictions. In addition, a continuous cytosolic binding region (C) was found near residues D4938, E4942, and D4945, which interact with the permeating Ca²⁺ and may influence the ion permeability as well. This is consistent with a previous experimental study showing that D4938 and D4945 determine the ion flux and selectivity³³.

a The contour plot of Ca²⁺ density on the R-z plane around RyR1. Four binding sites within the pore are designated as L, S, G, and C. The positions of the selectivity filter (SF) and gate (GT) constrictions are indicated with solid lines. The ion channel is divided into chambers by the dashed lines at the density saddles. The pore residues in close proximity are labeled on these lines. b, c Side views of the RyR1 channel. Only two chains are shown for clarity. The bottleneck residues (GLY-4894 at SF and GLN-4933 at GT), and the negatively charged residues at the binding sites L and C are shown as ball-and-sticks. Source data of (a) are provided as a Source Data file.

Full size image

The Ca²⁺ ions are fully hydrated during permeation

We calculated the number of oxygen atoms coordinated with the permeating Ca²⁺ and monitored from which residues these oxygen atoms were (water or protein). As shown in Fig. 4, the open-state RyR pore is relatively wide compared to K⁺ and Na⁺ channels. As the first solvation shell is around 2.4 Å from the Ca²⁺ and the radius of a water molecule is usually considered to be 1.4 Å, we consider the pore region with a radius < 4.0 Å to be the narrow pore region that contains the rate-limiting constriction sites. Interestingly, the average number of oxygen atoms coordinated with the Ca²⁺ was almost constant during the permeation process of Ca²⁺, as shown in the right panel of Fig. 4, and nearly all of these oxygen atoms were from water molecules within the narrow pore region. This clearly indicates that the permeating Ca²⁺ ions do not need to dehydrate when permeating through the open-state pore, and therefore the first solvation shell was intact in the narrow pore region. On the other hand, in the wide upper selectivity filter, one of the water molecules coordinated with Ca²⁺ can be replaced by the oxygen from E4900 (Fig. 4), simply because the strong electrostatic attraction between Ca²⁺ and E4900 pushed one of the coordinated water molecules away. A similar phenomenon was observed near D4938. It should be noted that, the water replacement at these sites are not caused by the steric dehydration when ions passing through a narrow pore as observed in K⁺ and Na⁺ channels, but rather due to strong electrostatic attraction from negatively charged residues in a wide vestibule (r > 5 Å), and therefore should not be considered as dehydration due to the permeation.

The left panel shows the pore radii along the pore axis of the open-state RyR1 (black solid), K⁺ channel (red dotted), and Na⁺ channel (blue dashed). The right panel shows the total number (black) of coordinated oxygen atoms around the calcium ions within the pore, and the contributions from protein (red) and water (blue) respectively. The significant contributions from protein oxygen are marked with corresponding protein residue IDs containing the oxygen. The transparent background shows the pore profile within the open-state RyR1 structure, and the narrow pore region is highlighted between the black dashed lines. Source data are provided as a Source Data file.

Full size image

The Ca²⁺ permeation pattern

The narrow region of the channel consists of the lower selectivity filter, the cavity, and the gate, which can be divided into two chambers by the saddle points of the Ca²⁺ density, as indicated by the dashed lines in Fig. 3a. The upper and lower chambers contain the binding sites S and G, respectively. The typical permeation pattern of Ca²⁺ ions is shown in Fig. 5a, from which it seems that permeation through the narrow pore region of the channel occur mainly in a one-at-a-time manner, meaning that there is only one Ca²⁺ residing in this narrow pore region, either in the upper or the lower chamber. The probability of both chambers being occupied by Ca²⁺ was only 2.4% in the trajectories, while the probability of only one chamber being occupied was 68.6%. Therefore, most of the time, only one Ca²⁺ can occupy the narrow pore region when permeating through the open-state channel. However, this does not mean that other Ca²⁺ are not involved in the permeation process. In fact, we can see that the exit of the lower Ca²⁺ from the narrow pore is often accompanied by the (transient) entering of a second Ca²⁺ (blue bars in Fig. 5a). This indicates that, Ca²⁺ ions do not have to be in direct contract or adjacent positions to influence each other in the RyR channel. This permeation pattern is distinct from K⁺ and Na⁺ channels, in which usually both the narrow selectivity filter and cavity can be occupied by multiple permeating ions at the same time^11,13,16. This is probably due to the fact that the electrostatic repulsion is much stronger between divalent ions than monovalent ions (Supplementary Fig. 4), and that RyR1 has a much shorter narrow selectivity filter region than typical K⁺ and Na⁺ channels (Fig. 4). Interestingly, when K⁺ permeating through the RyR1 channel, there is a significant fraction of the time in which two K⁺ ions were found occupying the narrow pore region (Fig. 5b), which is consistent with Miranda et al.ʼs work³⁴.

a For Ca²⁺, the occupied chamber accommodating either the selectivity filter or gate is shaded in gray unless both of them are occupied, which is indicated by blue bars depicted with small arrows. Representative configurations with one of the chambers occupied by a Ca²⁺ ion (in silver) are shown above. b For K⁺, similar to above, the gray shaded areas indicate one-ion occupancy while the blue shaded areas indicate two-ion occupancy. Source data are provided as a Source Data file.

Full size image

Free energy of Ca²⁺ within the RyR1

We calculated the PMF of Ca²⁺ within the pore using the equilibrium MD simulation trajectories. The result is shown in Fig. 6. As can be seen, the binding sites we observed in our simulations correspond well with the free energy wells in Fig. 6, and the free energy barrier between S and G corresponds to the rate-limiting step we found in the simulation trajectories, which is only ~9.8 kJ/mol. This is consistent with our finding that Ca²⁺ ions do not dehydrate its first solvation shell when passing through this region, as the dehydration energy of Ca²⁺ is very high. We believe the energy barrier here mainly represents an electrostatic effect, due to the existence of the multiple negatively charged residues in the upper SF vestibule, as well as the minor Ca²⁺ binding site G, which can accommodate one Ca²⁺ that will hinder other Ca²⁺ from approaching the cavity by electrostatic repulsion. Such a low free energy barrier also explains why we observed a large conductance of Ca²⁺ through RyR1. The energy barrier between the binding sites G and C is even lower (~2.5 kJ/mol), indicating the gate constriction is not the rate-limiting site and Ca²⁺ will not occupy the G site for a long time before exiting the narrow pore.

The binding sites and some relevant residues are labeled. The rate-limiting free energy barrier between S and G is estimated to be ~9.8 kJ/mol. Source data are provided as a Source Data file.

Full size image

The interaction between Ca²⁺ and protein is of great importance in studying Ca²⁺-mediated biological processes. Although the multi-site ion model cannot explicitly represent the charge transfer and polarization effect, the simulation results showed that our new model behaves much better than conventional single-point ion models. Moreover, our model is entirely consistent with the currently widely used non-polarizable force fields, such as CHARMM and AMBER, and therefore can be easily used in MD simulations. As shown in the result section, our seven-site Ca²⁺ model can reproduce the solvation properties of Ca²⁺, including the hydration energy, the first solvation shell size (first-peak position of RDF), the coordination number, and the residence time of water in the first solvation shell, as well as the Ca²⁺-protein binding energies in a more quantitative way that is comparable to quantum chemistry calculations. The optimized Ca²⁺ model was further validated by thermodynamic calculations and then utilized in the simulations of the RyR1 ion channel. In contrast to the conventional Ca²⁺ model, where ions get stuck in the channel in MD simulations, our new model resulted in the continuous permeation of Ca²⁺ through the channel, and the calculated conductance is in good agreement with the experimental electrophysiology result. Therefore, we believe this multi-site Ca²⁺ model can be widely used in simulating many Ca²⁺-involved bio-systems in addition to ion channels, and our simulations of RyR can provide detailed information about the Ca²⁺ permeation mechanism.

The ion binding and permeation mechanisms have been widely studied for K⁺ and Na⁺ channels with molecular dynamics simulations^{10,11,12,14,15,16,17}, but only rarely studied for Ca²⁺ channels with conventional ion models^30,34. Very recently, Heinz et al. performed extensive PMF calculations of various cations in the RyR1, and located the key binding sites within RyR1 as well as the monovalent cation permeation pathways³⁰. Miranda et al. used monovalent cations and Cl⁻ to explore the permeability and origin of the high selectivity of cations over anions of RyR2³⁴. Both work provided highly valuable insights on the ion permeability and selectivity of RyR channels, but the Ca²⁺ permeation details remain unclear due to the absence of a reliable Ca²⁺ model. From our MD simulations with the new multi-site Ca²⁺ model, the major Ca²⁺ binding sites in the open-state RyR1 were determined to be near the D4900 and G4894 residues on the luminal side, which is consistent with previous experiments and simulations^30,32,33,34, and the rate-limiting step of permeation is found to be the step of passing through the selectivity filter constriction, where Ca²⁺ ions dwell for a relatively long time before passing through due to the existence of a free energy barrier here (Figs. 5a and 6). Our simulation indicated that Ca²⁺ can easily accumulate in the wide upper selectivity filter, as shown by the high Ca²⁺ density from our equilibrium simulations (Fig. 3a). On average, there were about five Ca²⁺ or eight K⁺ in the selectivity filter in the presence of a 100-mV transmembrane potential (Supplementary Fig. 5). This indicates that K⁺ cannot fill up the electrostatic energy well in the selectivity filter as effectively as Ca²⁺ do, implying that the selectivity among cations is mainly achieved by the charge-space competition, as hypothesized by a previous experimental study³⁵. Since Miranda et al. already showed that the highly negatively charged SF vestibule is also responsible for the selectivity of cations against anions, we conclude that the highly charged rings in the SF is responsible for not only the high selectivity between cations and anions, but also the low selectivity among divalent and monovalent cations.

Na_V and TRP channels have a shorter SF, thus sharing a higher structural similarity with RyR channels than K_V channels, but their SFs are still longer than RyR1. Previous studies showed that there are multiple Na⁺ occupying the SF of Na_V channels, and they permeate by the so-called “loosely coupled knock-on” mechanism¹³, in which multiple partially hydrated Na⁺ interact each other within the SF and cavity. Darré et al. observed at least three monovalent or divalent cations can co-exist in the narrow pore region of TRPV1³⁶. Recent structural and MD studies of the TRPV6^37,38, a calcium-permeable ion channel, also showed that Ca²⁺ can tightly bind to the SF, and one or two additional Ca²⁺ need to approach it, line up and knock it off from the binding site. There are multiple Ca²⁺ occupying the narrow SF and cavity during this process. This is named the “knock-off” mechanism, which probably determines the ion selectivity of the Ca_V and TRPV6 channels^37,39,40. We did not observe such a well-organized “knock-on” or “knock-off” phenomena in our simulations of RyR1. As shown in Fig. 3a, there is only one Ca²⁺ binding site in the narrow SF, and the other Ca²⁺ ions in the upper SF do not line up, which may still push and knock off the lower Ca²⁺ in the SF, but in a loosely coupled manner.

From our simulations, Ca²⁺ ions permeate through the narrow pore region of the RyR1 channel following a “remote knock-on” mechanism and one-at-a-time manner (Fig. 5a), which is distinct from the well-studied direct knock-on mechanism in K⁺ channel and loosely coupled knock-on in Na⁺ channels^11,13,16. Previous studies have shown that multiple monovalent ions can enter the narrow pore region of the K⁺ and Na⁺ channels, and line up to facilitate the so-called ‘knock-on’ permeation, either tightly or loosely coupled^{10,11,12,13,14,15,16}. This was also observed for K⁺ in RyR³⁴. However, it is not the case for the Ca²⁺ in the RyR1 channel, as only one Ca²⁺ was observed in the narrow pore region (r < 4 Å) during the Ca²⁺ permeation in our MD simulations, and this region contains both the lower selectivity filter and cavity. There may be two reasons for this different permeation pattern, one is that the narrow selectivity filter of RyR is shorter compared to that of K⁺ and Na⁺ channels (Fig. 4), which can hardly accommodate multiple ions, and the second reason is that the electrostatic repulsion between divalent ions is much stronger than monovalent ions, which makes it energetically unfavorable for multiple divalent cations to sit side by side within a certain distance. In fact, we observed that when there was a Ca²⁺ in the cavity, other Ca²⁺ cannot stay in the lower selectivity filter (Fig. 5a). However, this does not mean other Ca²⁺ ions are not interacting with the permeating Ca²⁺. By calculating the electrostatic potential profile (Supplementary Fig. 4), we can see that the Ca²⁺ in the lower binding site can exert a large force on those Ca²⁺ in the SF, and vice versa. Therefore, even though there is only one Ca²⁺ ion in the narrow pore region, those Ca²⁺ in the SF can still electrostatically push it. We term this “remote knock-on”, which was first proposed from a simple model ion channel for the enhancement of ion permeability⁴¹. This unique remote knock-on and one-at-a-time permeation pattern observed in our MD simulations occurs due to both structural features of RyR and strong repulsive interactions of divalent cations, which was not observed in previously studied ion channels.

Another distinct permeation feature of RyR is that Ca²⁺ ions are nearly fully hydrated during the translocation along the narrow pore region. Previous studies have shown that K⁺ is nearly fully dehydrated^11,16 and Na⁺ is partially dehydrated^13,14 during permeation, meaning all of the water molecules or at least several water molecules within the first solvation shell of the ions must be removed or replaced by other residues at the narrowest constriction sites. In fact, this dehydration process was believed to be the key factor responsible for the ion selectivity of the channels^13,17,42, as different ions have different sized solvation shell and selectivity filters with different steric and chemical features can discriminate them by the free energy difference during the dehydration process. Interestingly, we did not observe such a dehydration behavior when Ca²⁺ was permeating through the RyR1 channel. As shown in Fig. 4, the number of water molecules are nearly constant in the narrow pore region and the oxygen atoms coordinating with the Ca²⁺ from the protein is nearly zero throughout, indicating that the first solvation shell of the Ca²⁺ is intact during permeation, and no dehydration occurred when the ion passing through the narrow constriction sites in the selectivity filter and lower gate. This finding provides a clear picture on the exact hydration states of Ca²⁺ ions as they pass through the pore, which is consistent with earlier speculation that ion dehydration may not be a significant component of selectivity or permeation⁴³. From the structural point of view, the open-state RyR1 is wider than K⁺ and Na⁺ channels in the selectivity filter region (Fig. 4), and therefore sterically allows the Ca²⁺ to permeate with its first solvation shell intact. On the one hand, this makes the free energy barrier of Ca²⁺ translocation lower, allowing highly efficient Ca²⁺ permeation to regulate ion concentration in the cytosol for muscle contraction and heartbeat, which is otherwise hard to imagine as Ca²⁺ has a much higher dehydration energy compared to K⁺ and Na⁺. On the other hand, this probably also leads to the weaker ion selectivity of RyR1, compared to the highly selective K_V, Na_V, and Ca_V channels, since the powerful dehydration-based selectivity mechanism cannot be utilized here.

It should be noted that, we adopted a truncated RyR1 in our simulations. Although this covers the whole transmembrane region, it does not represent the entire protein. The large cytoplasmic domain, which is important for the gating and regulation of the RyR permeation, was not considered in the present study. Interestingly, Heinz et al.’s recent study of a less truncated RyR1 showed very similar binding sites to ours, although the specific PMF profile was not identical due to the different ion models and truncation sites adopted³⁰. Therefore, we think the truncation may affect the ion motion after entering the cytoplasmic side, but the permeation behavior through the transmembrane pore and the binding sites along this path should remain valid. We also calculated the conductance of the pore using the fluctuation-dissipation theorem by monitoring the spontaneous Ca²⁺ permeation events and obtained a value of ~139 pS (Supplementary Fig. 6), which is in very good agreement with the value calculated under a transmembrane potential mentioned above (~141 pS), indicating nice convergence of the simulations.

In summary, we developed a multi-site Ca²⁺ model which behaves well when studying its interactions with proteins. It is entirely consistent with the widely used CHARMM force field, and we are expanding to other force fields and divalent ions as well. With the new Ca²⁺ model, we revealed the detailed Ca²⁺ permeation process through the open-state RyR1, and discovered that multiple Ca²⁺ ions can accumulate in the wide upper selectivity filter, and then permeate through the narrow pore region following a remote knock-on and one-at-a-time pattern, in which the permeating Ca²⁺ is fully hydrated with an intact first solvation shell (Fig. 7), distinct from the widely studied K⁺ and Na⁺ channels. The transmembrane pore region is wide enough to allow hydrated cations to permeate and therefore facilitate a large ion conductance. The selectivity filter of RyR1 is shorter and wider than the well-studied K_V and Na_V channels and cannot discriminate cations by their solvation shells, leading to a weak selectivity among cations by the so-called charge-space competition mechanism. Our simulation results shed new light on the high efficiency and weak cation selectivity of the Ca²⁺ permeation through the RyR channels.

There is only one Ca²⁺ in the narrow pore region, with a longer residence time at the selectivity filter constriction site (darker color), and a shorter residence time at the gate constriction (lighter color).

Full size image

Simulations for parameterization

All the MD simulations for the ion model parameterization were performed with OpenMM (version 7.0.1)⁴⁴, as the package is highly flexible and can be easily customized with a Python interface, and the CHARMM force field⁴⁵ (version C36_Aug15) was utilized in this work.

We built an ion-in-water system to optimize the ion properties in the water, which consists of one multi-site Ca²⁺ ion in a cubic water box of 3 × 3 × 3 nm³. The TIPS3P water model was used in consistency with the CHARMM force field. In our simulations, NPT ensembles were generated by integrating the Langevin dynamics with a time-step of 2 fs and a collision frequency of 5 ps⁻¹. The temperature was maintained at 298 K, and the pressure was regulated at 1 bar using a Monte Carlo barostat. Water molecules were kept rigid during simulations, and the cutoff of non-bonded interactions was 1 nm. For calculations of the radial distribution function, coordination number, and residence time, a 20-ns trajectory was generated. For hydration free energy calculations, 1-ns trajectories were generated for each of the 14 alchemical states (please also see below). The hydration free energy was estimated using the python implementation of the multistate Bennett acceptance ratio downloaded from https://github.com/choderalab/pymbar. The properties related to the radial distribution function were calculated using MDTRAJ⁴⁶.

To optimize the ion-protein interactions, we used the Ca²⁺-Protein systems previously investigated by Li et al.⁹, and used OpenMM and CHARMM force field to perform the calculations.

Fitness functions for optimization

We define two fitness functions to optimize the parameters of our multi-site Ca²⁺ model. One is designated as the protein-fitness function \((\lambda _{\mathrm{p}}^2)\) and the other as the water-fitness function \((\lambda _{\mathrm{w}}^2)\). The protein-fitness function is defined against a dataset of quantum-mechanically calculated Ca²⁺-protein binding energies, \(\Delta E_{p,f}^{{\mathrm{QM}}}\), where p identifies the index among N_p (=10) proteins, and f indexes N_f (=21) trajectory snapshots for each protein. The formula for \(\lambda _{\mathrm{p}}^2\) is

where \(\overline {\Delta E} _p^{\mathrm{X}} = \mathop {\sum}\nolimits_f {\Delta E_{p,f}^{\mathrm{X}}/N_{\mathrm{f}}\, ({\mathrm{X}} = {\mathrm{QM}},{\mathrm{MM}})}\), and MM indicates the molecule-mechanical results. Due to limitations of the methodology level and the basis-set size used in quantum-mechanical calculations, a systematic correction of binding energies (ΔE_C = 10 kcal/mol) is added to \(\overline {\Delta E} _p^{{\mathrm{QM}}}\) following Li et al.’s strategy⁹.

In addition to Ca²⁺-protein interactions, we also optimize our model to reproduce its energetical, structural, and dynamic properties in water. Specifically, these properties include the hydration free energy (ΔG_h), the first-peak position of radial distribution function (R₁), the coordination number (N_C), and the residence time of water in the first coordination shell (τ_R). The water-fitness function is defined as

where quantities with superscripts of t, e, and w stand for theoretical, experimental, and weighting values, respectively.

Target properties of Ca²⁺ in water

We followed the approach of Mamatkulov et al.⁴⁷ to determine the target hydration energy of Ca²⁺ (−1504 kJ/mol) as

where ΔG(CaCl₂) is the measured hydration energy of CaCl₂⁴⁸, ΔG_h(Cl⁻) the theoretical hydration energy determined from Smith-Dang parameters⁴⁹ for Cl⁻, ΔG_press the energy needed to compress one mole of ion gas at 1 atm into a liter, and ΔG_surf the energy change of crossing the air-water interface for 1 mol of ions. Marcus assessed the Ca–O internuclear distances in calcium salt solutions from different studies and concluded that the generally consistent result is 0.242 nm⁵⁰, which was used as our target value for \(R_1^{\mathrm{e}}\). A recent neutron diffraction study⁵¹ revealed that the average number of water molecules in the first hydration shell of Ca²⁺ is close to 7, which we took as the target value for N_C^e. It is difficult to experimentally determine the residence time of water molecules in the first hydration shell of the calcium ion⁵², and a nuclear magnetic resonance (NMR) study estimated its value to be <100 ps²⁷. Since the existing Ca²⁺ models generally overestimate τ_R, we set \(\tau _{\mathrm{R}}^{\mathrm{e}}\) as zero in the course of optimization and checked the final result so as to be the same order of magnitude as 100 ps. The numerical values for the targeted experimental properties and their corresponding weights are listed in Supplementary Table 2.

Target binding energies of Ca²⁺ with proteins

We followed Li et al.’s⁹ protocol to calculate the binding energies between Ca²⁺ and proteins. As noted by Li et al., the default single-site model of Ca²⁺ in the CHARMM C36 force field generally overestimates the Ca²⁺-protein binding energies by 150–200 kcal/mol with respect to quantum-mechanical binding energies. They selected 10 high-resolution crystal structures of Ca²⁺-binding enzymatic proteins (PDB IDs shown in Fig. 1) and performed MD simulations of the solvated proteins. From the equilibrated trajectories, 21 conformation snapshots were extracted for each protein and truncated to include atoms within a sphere of ~0.55 nm around the ion. Quantum mechanical calculations were then carried out with the truncated models to obtain a dataset of binding energies, which is used in this work to optimize the multi-site Ca²⁺ model. For the systematic correction of binding energies (ΔE_C) that is added to \(\overline {\Delta E} _p^{{\mathrm{QM}}}\), we use an estimation of 10 kcal/mol as Li et al. did in their work. All the ab initio and Drude binding energies were from the original work by Li et al.⁹.

Calculation of properties from simulations

The theoretical hydration energy, \(\Delta G_{\mathrm{h}}^{\mathrm{t}}\), consists of two terms, i.e.

The first term refers to the free energy change corresponding to alchemically switching off the ion-water electrostatic and LJ interactions, and the second term is a correction due to the finite-size simulation box. The alchemical transformation for turning off the electrostatic interaction was performed with 10 intermediate states, during which the charge of ion was sequentially reduced by a tenth of the total charge. After that, it took four uniform intermediate states with a soft-core potential (α = 0.5)⁵³ to annihilate the LJ interaction. MD trajectories (1 ns) of one ion in a cubic water box (3 nm) were used to estimate the alchemical free energy by the method of multistate Bennett acceptance ratio⁵⁴. For the finite-size correction, we took the formula derived by Hummer et al.⁵⁵,

where Ze is the ion charge, ξ_EW is the Wigner potential, ε_r (=82) is the relative dielectric constant, and L (=3 nm) is the box size.

The first-peak position of the radial distribution function \(R_1^{\mathrm{t}}\) was calculated using the MDTRAJ software⁴⁶ from a 20-ns trajectory, and the coordination number \(N_{\mathrm{C}}^{\mathrm{t}}\) was computed by the integration of the first peak. We followed the definition described by Impey et al.⁵⁶ to calculate the residence time \(\tau _{\mathrm{R}}^{\mathrm{t}}\). First, the residence time distribution \(n_{{\mathrm{ion}}}(t)\) was computed as

where N_t is the number of time frames and N_w is the number of water molecules. P_j(t_n, t; t^*) takes a value 1 if the water molecule j stays in the first hydration shell from t_n to t_n + t without leaving for any period larger than t^*, and takes the value 0 otherwise. Then, n_ion(t) was fitted with \(n_0\exp ( - t/\tau _{\mathrm{R}}^{\mathrm{t}})\) to obtain \(\tau _{\mathrm{R}}^{\mathrm{t}}\).

The molecule-mechanical Ca²⁺-protein binding energies \((\Delta E_{p,f}^{{\mathrm{MM}}})\) were calculated as the potential energy difference between Ca²⁺-bound and Ca²⁺-free configurations. Since our model of Ca²⁺ has a geometry of octahedron, the total energy of a Ca²⁺-bound system depends on the orientation of the octahedron. Thus, we first minimized the energy with respect to the orientation of the octahedron and obtained E_min(Ca²⁺ + Protein), and then calculated the Ca²⁺-protein binding energy using the formula ΔE^MM = E_min(Ca²⁺ + Protein) − E(Ca²⁺) − E(Protein).

Optimization strategy

First, we did a thorough scan and prescreening, using the conventional optimization routines (such as conjugate gradient and basin hopping), random sampling, and sifting the parameter space by the properties described above in sequence. Initially, we focused on six parameters (b_CD, Q_C, ε_C, σ_C, ε_D, σ_D), i.e., treating interactions with water and protein on an equal footing. The searched parameter space as shown in Supplementary Table 3 is determined from previous works^25,26, and we substantially extended the range of Q_C, ε_C, and ε_D based on our optimization results. However, after systematic optimization, the calculated water properties and protein binding energies were still slightly off from the target values. Therefore, we further split (ε_C, σ_C) into (\(\varepsilon _{\mathrm{C}}^{{\mathrm{NW}}}\), \(\sigma _{\mathrm{C}}^{{\mathrm{NW}}}\)) and (\(\varepsilon _{\mathrm{C}}^{\mathrm{W}}\), \(\sigma _{\mathrm{C}}^{\mathrm{W}}\)), thus calculating the interactions of Ca²⁺ with protein and water separately. In fact, we took a two-step strategy to optimize our multi-site model. First, the water-fitness function \(\lambda _{\mathrm{w}}^2\) was minimized to obtain optimal values for (b_CD, Q_C, \(\varepsilon _{\mathrm{C}}^{\mathrm{W}}\), \(\sigma _{\mathrm{C}}^{\mathrm{W}}\), ε_D, σ_D). Second, with the parameters of (b_CD, Q_C, ε_D, σ_D) fixed, the protein-fitness function \(\lambda _{\mathrm{p}}^2\) was scanned with a two-dimension grid to find the best values of (\(\varepsilon _{\mathrm{C}}^{{\mathrm{NW}}}\), \(\sigma _{\mathrm{C}}^{{\mathrm{NW}}}\)) that reproduce quantum-mechanical Ca²⁺-protein binding energies. Due to the roughness of \(\lambda _{\mathrm{w}}^2\) owing to the stochastic nature of the calculation, we developed a simple yet effective steep-descent-like algorithm that alternatingly searches along each dimension of the parameter space instead of the gradient. The parameter searching along each dimension was performed at linearly spaced steps, which were then finely tuned near the local minimum of a finite interval centering at the current minimal point. We used this strategy to further optimize the top parameters obtained in the initial stage of searching, which indeed yielded more satisfactory results. We found that the minimization mainly took place in the dimensions other than \(\varepsilon _{\mathrm{C}}^{\mathrm{W}}\) and ε_D, which indicated that the landscape of the water-fitness function is much smoother along these two dimensions. Therefore, we selected a couple of combinations of (\(\varepsilon _{\mathrm{C}}^{\mathrm{W}}\), ε_D) based on the scanning results and carried out further optimization until the fitness value was sufficiently low. The final optimized parameters are presented in Supplementary Table 4.

Details of permeation simulations

To study the permeation of Ca²⁺ through RyR1, we performed MD simulations with GROMACS⁵⁷ version 5.1.3 with the CHARMM36 force filed⁴⁵ and the TIPS3P water model. The channel pore domain (residue 4820–4956) of the cryo-EM structure of the open-state RyR1²¹ (PDB ID: 5TAL [https://doi.org/10.2210/pdb5TAL/pdb], resolution: 4.3 Å) was used as the starting structure. The resolution is not very high, but should be good enough for the present study as the open-state pore of RyR1 is relatively wide and the subtle configurations within the pore may affect the ion permeation to a minor extent, and the side chains are free to move in our MD simulations to sample low-energy local conformations. The OPM⁵⁸ web server and Membrane Builder in CHARMM-GUI⁵⁹ were used to build four POPC-RyR1 simulation systems.

A POPC-RyR1 simulation system with ~150 mM of calcium ions and a POPC-RyR1 simulation system with ~250 mM of potassium ions were built to calculate the conductance of the calcium ions and potassium ions with a transmembrane potential of 100 mV. The ion concentrations were chosen to match the experimental conditions under which the conductance was measured³¹. As K⁺ exhibits a higher saturation concentration than Ca²⁺, we used a [Ca²⁺] of 150 mM and [K⁺] of 250 mM in our simulations. Six 500-ns trajectories were simulated for the former system, and three 300-ns trajectories for the latter. By counting the number of permeation events (N_p) for a trajectory of certain time length (t_traj) under a transmembrane potential of V_tm, the ion conductance C_ion was calculated with:

where Q_ion is the charge of the permeating ion.

Another POPC-RyR1 simulation system with ~150 mM calcium ions was built to calculate the conductance of the calcium ions with our new calcium model with the same transmembrane potential of 100 mV. As the ion model cannot be simulated as a rigid body with GROMACS, we used bond and angle restraints to make the multi-site Ca²⁺ model as rigid as possible, which affect the thermodynamics of the Ca²⁺ to a negligible extent (please refer to Supplementary Table 5). Six independent 500-ns trajectories were conducted. And, one 800-ns trajectory was run for the same system but without an electric field to obtain the stable binding sites of the calcium ions in the open-state RyR1, which would be more comparable to ion positions in the cryo-EM density maps than with an electric field.

Based on the 800-ns simulation trajectory, we estimated the free energy profile of Ca²⁺ along the pore. We first calculated the Ca²⁺ density in the simulation system with a grid-spacing of 1 Å, and then we determined the maximum density of Ca²⁺ at a certain position along the axis of the pore: \(\rho ^ \ast \left( z \right) = {{\max }_{x,y} }\rho (x,y,z)\). Then the free energy profile is calculated as \(G\left( z \right) = - RT\ {\mathrm{ln}}\ \rho ^ \ast \left( z \right)\).

All the simulation systems were first equilibrated with the standard CHARMM-GUI equilibration protocol followed by the production simulations with the harmonic position restraints applied on the α carbon of the protein (with a force constant of 1000 kJ mol⁻¹ nm⁻²). The position restraints were applied to make sure the open-state conformation was stable in the absence of the large intracellular domain and the caffeine and ATP bound to it, but the atoms other than α carbon, including all the side chains of the transmembrane domain, were fully flexible in our simulations, allowing a certain degree of flexibility of the protein (Supplementary Fig. 7). For all the production simulations, the periodic boundary conditions were used and the time step was 2 fs. The v-rescale algorithm with a time constant of 0.5 ps was used to maintain the temperature at 310 K, and the Parrinello-Rahman algorithm⁶⁰ with a time constant of 1 ps was used to maintain the pressure at 1.0 bar. The Particle-mesh Ewald method⁶¹ was used to calculate electrostatics, and the cut-off length of the van der Waals interaction was 1.2 nm.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.