In silico assessment of the conduction mechanism of the Ryanodine Receptor 1 reveals previously unknown exit pathways

The ryanodine receptor 1 is a large calcium ion channel found in mammalian skeletal muscle. The ion channel gained a lot of attention recently, after multiple independent authors published near-atomic cryo electron microscopy data. Taking advantage of the unprecedented quality of structural data, we performed molecular dynamics simulations on the entire ion channel as well as on a reduced model. We calculated potentials of mean force for Ba2+, Ca2+, Mg2+, K+, Na+ and Cl− ions using umbrella sampling to identify the key residues involved in ion permeation. We found two main binding sites for the cations, whereas the channel is strongly repulsive for chloride ions. Furthermore, the data is consistent with the model that the receptor achieves its ion selectivity by over-affinity for divalent cations in a calcium-block-like fashion. We reproduced the experimental conductance for potassium ions in permeation simulations with applied voltage. The analysis of the permeation paths shows that ions exit the pore via multiple pathways, which we suggest to be related to the experimental observation of different subconducting states.


On the modeling and equilibration
Modeling Table 1 provides a list of loop segments that were unresolved in all of the referenced cryo-EM datasets and thus needed to be fitted using the Modeller software [1][2][3][4] . All loop models went through multiple successions of conjugate gradient optimization and simulated annealing. The final model was chosen among 100 independent runs based on the overall score. Table 1. Number of loop segments modelled with the Modeller software and their respective sizes. E.g. 11 loop segments with lengths between 1 and 5 residues were modeled in total. Of those, 2 are part of the reduced model and none is located in the pore region.

Number of residues
Number of loops full model reduced model pore region In addition to loop segments, homology models for two larger unresolved parts were constructed externally. Those segments are given by H1296-N1418 and Q3570-L3644.
The templates used for the full model are shown graphically in Fig. 1.       For the first 30 ns of the simulation run, open-state restraints were kept on every 4th Cα atom. Later the restraints were released, causing an RMSD increase of ≈ 0.1 nm. The umbrella sampling simulations as well as the applied voltage simulations were started from the structure at 60 ns and the equilibration was continued to a total of 200 ns to assess its convergence. The pore RMSD leveled off at ≈ 0.35 nm.
Note that the RMSD does not exactly start at 0, because after the steered simulation run, the model is very similar but not perfectly identical to its template due to thermal fluctuations.

On the Gromacs implementation of the energy-step method
The energy-step, used to restrain ions in their respective compartment, was implemented into Gromacs by expanding upon the set of existing flat-bottom potentials. Fig. 7 shows an example of a energy step. The potential is flat in front and behind a transition region -i.e. there will be no force acting on the particles subject to the potential. A constant force f is exerted in the transition region with width d. The height of the step is given by d · f . Just as any other position restraint, energy-steps can be added to a moleculetype in the topology file. In contrast to the existing implementation of flat-bottomed potentials, the reference points are not determined by the initial particle positions, but the position of the step can be chosen freely as a parameter. Our implementation allows for energy steps in x,y and z directions. Multiple steps can be combined; for instance six steps could be used to restrain atoms in a cuboid-shaped volume.

5/9
Adding the lines below to a moleculetype-section in the topology file restraints the movement of the first atom (i=1) in z-direction (g=11).

[ p o s i t i o n r e s t r a i n t s ]
; i f u n c t g z0 f 1 2 11 0 . 2 5 10000 The step is located at z0=0.25 nm and the force in the transition region is provided by f=10000 kJ/mol/nm. Steps facing in x and y-direction are parameterized by g=9 and g=10 respectively. Smaller values of the g parameter lead to the stock-Gromacs flat-bottom potentials, as described in the section 'Flat-bottomed position restraints' of the Gromacs manual.
Our implementation of the energy-step is available on request.

On the convergence of umbrella sampling
We address the quality of convergence of the umbrella sampling simulations on the example of Ca 2+ . The histogram, provided in figure 8, shows that the entire interval along the reaction coordinate is well sampled and that neighboring windows strongly overlap. To assess the convergence in a more quantitative fashion, we calculated the potential of mean force while every 2nd, 3rd and 4th sampled windows were removed from analysis respectively. The results are summarized in Fig. 9, where the PMF based on the full set of windows is shown as well. The maximal deviation from the full PMF amounts to 6.6 kJ/mol. This is in agreement with the 1σ error range of ±5.4 kJ/mol as obtained by bootstrapping. Furthermore we shorten the analyzed length of each umbrella window by up to 400 ps and compare it to the full data set in Fig. 10.  To further address the question whether 0.6 ns for each of the 209 umbrella windows are sufficient to converge the PMF as well as relax the channel properly, we increased sampling for Ca 2+ ions and observed the corresponding changes to the PMF. 10 umbrella windows, located in the core of the channel, were substituted by longer simulation runs, which lasted 3 ns.
Comparison between the old PMF and the new PMF with increased sampling, shown in Fig. 11, shows no significant change.

Calculating hydration numbers
We calculated hydration numbers for of both first and second hydration shells for all ion types along the pore axis. Hydration numbers were determined by counting the number of water molecules around the ion closer than cut off radii, corresponding the first two hydration shells, in the final frames of each of the umbrella sampling windows. Local averages along the pore axis were taken using a Gaussian filter with a standard deviation of six windows. The cut off radii were determined based the radial distribution functions (RDFs) of a single ion in bulk water, which was simulated for 10 ns. The first and second local RDF

On the classification of the ion permeation trajectories
A heuristic approach was adopted, in order to distinguish and to quantify the different observed modes of exiting the channel, such as visualized in Fig. 5 of the main text.
Classification was done using a scoring function which scores the likelihood that ion i passed through exit j during the course of its trajectory { x Here, N f is the number of captured frames and r j is the exit point of pathway j. The parameter d determines up to which distance to an exit point a trajectory still contributes significantly to the score. Reasonable classification was obtained by choosing d = 1 nm.
Trajectory i was considered unclassifiable if the second largest entry in {S i1 , S i2 , S i3 } was greater than 80% of the largest entry. In the other case, the ion was classified as having passed through the highest scoring exit.