A cooperative knock-on mechanism underpins Ca2+-selective cation permeation in TRPV channels

TRPV5 and TRPV6 are unique among TRP channels due to their high Ca2+ selectivity, while most other members of this ion channel family do not select for a specific cation type. Ives et al. used biomolecular simulations and in silico electrophysiology to determine the mechanism underlying this unusual Ca2+ selectivity.


Introduction
The significance of Ca 2+ in cellular function was first recognized by Sydney Ringer in 1883, who demonstrated that minute amounts of calcium were required for the contraction of cardiac muscle (Ringer, 1883). Ca 2+ is now recognized as a versatile signaling agent, with cellular Ca 2+ concentrations impacting a broad array of physiological processes ranging from cell proliferation to cell suicide (Carafoli, 2002;Berridge et al., 2003;Clapham, 2007;Patel, 2019). However, the cytoplasmic concentration of Ca 2+ ions is usually kept low due to cytotoxic consequences (Bagur and Hajnóczky, 2017). Therefore, the controlled opening of channels in cellular and organellar membranes is one of the required mechanisms to allow the influx of this ion from the exoplasm and internal storage compartments into the cytoplasm; this subsequently initiates the Ca 2+ signaling cascade. The question of how Ca 2+ channels selectively permeate Ca 2+ in low concentrations over vastly more abundant Na + ions and yet conduct them at high rates has been a longstanding matter of fascination for ion channel researchers (Corry et al., 2001;Hille, 2001).
A key example of ion channels that mediate Ca 2+ permeation across the cytoplasmic membrane is the transient receptor potential (TRP) channel superfamily. In their open state, these polymodal signal-detecting TRP channels allow the transmembrane flux of cations down their electrochemical gradient, thereby increasing the intracellular Ca 2+ and Na + concentration (Ramsey et al., 2006). The malfunction of TRP channels underlies a wide range of diseases, and they are therefore of immense biomedical importance, serving as drug targets for a variety of existing and candidate drugs (Moran, 2018).
TRP channels assemble primarily as homotetramers to form functional ion channels. A conserved structural feature across all TRP channels is the presence of six transmembrane helices (S1-S6) per subunit, forming two distinct transmembrane domains; a four-helix bundle comprising of helices S1-S4 forming the voltage-sensing like domain; and the pore-forming domain consisting of helices S5 and S6 (Hilton et al., 2019).
A four-residue ion selectivity filter (SF) is located at the entrance of the channel pore (Fig. 1). In addition to this conserved transmembrane architecture, members of the TRP superfamily display highly diverse extramembrane loops and Nand C-terminal domains between the different subfamilies (Van Goor et al., 2020).
While they are all cation selective, most TRP channels electrophysiologically characterized to date show only limited discrimination between cation types, as well as between divalent and monovalent cations. However, the TRPV5 and TRPV6 channels are unique due to their high selectivity for Ca 2+ cations over Na + cations (P Ca /P Na ∼100:1 from reversal potential measurements; Vennekens et al., 2000;Yue et al., 2001). Phylogenetic analysis has demonstrated that the TRPV5 and TRPV6 channels of vertebrates originated from an ancestral TRPV5/6 gene, which then diverged to form TRPV5 and TRPV6 from a duplication event after speciation (Flores-Aldama et al., 2000). Both of these channels are constitutively active due to basal levels of phosphatidylinositol 4,5-bisphosphate (PI(4,5)P 2 ) in the cellular membrane and play a key role in Ca 2+ homeostasis in the body (Van Goor et al., 2017). Despite their characteristic Ca 2+ selectivity, both channels have been shown to permeate monovalent cations such as Na + when divalent cations are absent Vennekens et al., 2000;Voets et al., 2003;Bödding and Flockerzi, 2004). By contrast, the remaining members of the TRPV subfamily, TRPV1-4, permeate both Ca 2+ and Na + cations, even in the presence of high Ca 2+ concentrations, although they are still slightly Ca 2+ selective, with a permeability ratio P Ca /P Na ∼10:1. These channels gate in response to a number of stimuli, including raised temperature-in particular the archetypal member TRPV1 (Caterina et al., 1997;Tominaga et al., 1998;Caterina et al., 2000), which has led to TRPV1-4 being referred to as thermoTRPV channels-as well as endogenous and exogenous ligands.
In recent years, MD simulations have been successfully employed to shed light on ion channel function and the mode of action of channel-acting drugs in atomistic detail, for instance on K + channels (Köpfer et al., 2014;Kopec et al., 2018), Na + channels (Ulmschneider et al., 2013;Ke et al., 2014), Cl − channels Figure 1. Structure of the truncated construct of TRPV5 of Oryctolagus cuniculus used in this work, from the extracellular side (top left) and in the plane of the lipid bilayer (bottom left). In this study, the pore is defined as the region between the constrictions of the channel, namely the top residue of the SF (referred to as the α-position of the SF) and the hydrophobic lower gate (right). EC, extracellular side; IC, intracellular side of the membrane. (McKiernan et al., 2020;Chavan et al., 2020), and ligand-gated ion channels (Sauguet et al., 2013). However, in the case of Ca 2+permeating channels, the conventional point-charge models used to describe uncoordinated Ca 2+ ions have historically been inaccurate due to the neglected effects of electronic polarization (Kohagen et al., 2014a). This has resulted in an over-estimation of the binding energies between Ca 2+ and proteins, hindering an accurate study of Ca 2+ permeation in channels. Previous efforts to resolve this overestimation have included polarizable force field approaches (Li et al., 2015) and rescaled Ca 2+ charges (Kohagen et al., 2014a;Kohagen et al., 2014b). More recently, Zhang et al. (2020) published a new multisite Ca 2+ model specifically optimized for Ca 2+ -protein interactions (Zhang et al., 2020). This model correctly replicated experimental values for the hydration-free energy and the number of coordinated water molecules in the first solvation shell and showed Ca 2+ -protein binding energies comparable to the quantum mechanical and polarizable models. The multisite Ca 2+ model has previously been used to investigate ion permeation of Ca 2+ in a range of channels, including the type-1 ryanodine receptor (Zhang et al., 2020;Liu et al., 2021), α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptors (Schackert et al., 2022), and recently the E protein of SARS-CoV-2 (Antonides et al., 2022 Preprint).
In the present work, we set out to elucidate the molecular basis of Ca 2+ selectivity and permeation in the TRPV channel subfamily. We conducted atomistic molecular dynamics (MD) simulations of TRPV channels under transmembrane voltage and compared the cation permeation mechanism observed in the Ca 2+ -selective TRPV5 and TRPV6 channels to the permeation mechanism in two exemplar non-selective TRPV channels, TRPV2 and TRPV3. In total, we observed 2,851 full ion traversals from 17.25 µs of MD simulations, allowing us to decipher the permeation mechanisms and principles of ion selectivity in the TRPV family with statistical power. Our findings suggest that ion conduction in TRPV channels proceeds via a cooperative knock-on mechanism involving multiple ion binding sites. The degree of cooperativity in ion permeation, linking the multiple binding sites, determines the degree of ion selectivity in the channels.

Materials and methods
TRPV system construction Truncated TRPV simulation systems consisting of the membranedomain of the channels were constructed as described in Table 1. The systems were built using the CHARMM-GUI server (Jo et al., 2008). The charged N-and C-terminal residues were neutralized by capping with acetyl (ACE) and N-methylamide (CT3) groups, respectively, and all missing non-terminal residues were modeled . In the case of the TRPV5 system, the parameters for PI(4,5) P 2 were generated using the CHARMM General Force Field (Vanommeslaeghe, 2010) through the ligand reader and modeler in CHARMM-GUI (Kim et al., 2017).
The structures were aligned in the membrane using the PPM server (Lomize et al., 2012), inserted into a 1-palmitoyl-2-oleoylsn-glycerol-3-phosphocholine bilayer of 150 × 150Å size using the CHARMM-GUI membrane builder (Jo et al., 2007;Wu et al., 2014), and then solvated. Ions were added using GROMACS 2020.2 (Abraham et al, 2015;Lindahl et al, 2020) to neutralize any system charges and add ions to a concentration of either 150 mM NaCl, 150 mM CaCl 2 , or a mixture of 75 mM NaCl and 75 mM CaCl 2 . In the case of simulations containing Ca 2+ , the standard CHARMM36m Ca 2+ ions were then replaced with the multisite Ca 2+ model of Zhang et al. (2020). Hydrogen mass repartitioning of the system was used to allow the use of 4 fs time steps in simulations of NaCl solutions. The multisite Ca 2+ model used for simulations of CaCl 2 , however, is incompatible with a 4 fs time step, and therefore any simulations including Ca 2+ cations were performed with hydrogen mass repartitioning but at a time step of 2 fs. The protein was restrained in the openstate by applying harmonic restraints on the α-carbon atoms of the lower gate residues (see Table 1).

MD simulations details
All simulations were performed using GROMACS 2020.2 (Abraham et al, 2015;Lindahl et al, 2020) and the CHARMM36m force field for the proteins, lipids, and ions (Huang et al., 2017). The TIP3P water model was used to model solvent molecules (Jorgensen et al., 1983). The system was minimized and equilibrated using the suggested equilibration inputs from CHARMM-GUI (Lee et al., 2016). In brief, the system was equilibrated using the isothermal-isobaric (NPT) ensemble for a total time of 1.85 ns with the force constraints on the system components being gradually released over six equilibration steps. The systems were then further equilibrated by performing a 15 ns simulation with no electric field applied. To prevent closing of the lower gate of the pore, harmonic restraints were applied to maintain the distance between the α-carbon atoms of the lower gate residues of each respective chain (Table 1). To drive ion permeation, an external electric field was applied by using the method of Aksimentiev and Schulten (2005) to production simulations with an E 0 of −0.03 V nm −1 ; this resulted in a transmembrane voltage of ∼410 mV with negative polarity in the intracellular region. The temperature was maintained at 310 K using the Nosé-Hoover thermostat (Evans and Holian, 1985) and the pressure was maintained semi-isotropically at 1 bar using the Parrinello-Rahman barostat (Parrinello and Rahman, 1981). Periodic boundary conditions were used throughout the simulations. Long-range electrostatic interactions were modeled using the particle-mesh Ewald method (Darden et al., 1993) with a cutoff of 12Å. The linear constraint solver algorithm (Hess et al., 1997) was used to constrain bonds with hydrogen atoms. All individual simulations were 250 ns long and repeated five times for each system, as summarized in Table 2. Additional details for all simulations of Ca 2+ -selective and non-selective TRPV channels are described in Tables S1 and S2, respectively, and details for additional control simulations are described in Table S3.

Calculating conductance and selectivity from in silico electrophysiology experiments
The conductance of the channels (G channel ) was calculated according to Eq. 1, where N p is the number of permeation events, Q ion is the charge of the permeating ion in Coulomb, t traj is the length of the trajectory, and V m is the transmembrane voltage. The mean conductance and standard error were calculated from overlapping 50 ns windows of the trajectory.
The selectivity (P Ca /P Na ) was calculated as the ratio between the total sum of Ca 2+ permeation events and the total sum of Na + permeation events from fivefold replicated 250 ns simulations in dicationic solutions with 75 mM CaCl 2 and 75 mM NaCl.
Identification of cation binding sites from MD simulations of TRPV channels Cation binding sites were identified by plotting timeseries of each permeating ion with respect to their position along the pore axis. To further validate this, the Featurizer function of PENSA (Vögele et al., 2021;Vögele et al., 2022) was used to identify the 12 most occupied ion binding sites, as determined by 3-D density maxima within 7Å of the protein. This analysis was performed on a trajectory of concatenated fivefold replicated 250 ns simulations with a 200 ps time step from monocationic simulations.
Calculating ion occupancy probabilities and residence times in the identified cation binding sites The ion occupancy (O ion ) of the identified cation binding sites was calculated by dividing the number of frames (N occupied ), in which an ion's center of geometry is within 3.5Å of the center of geometry of the ion coordinating binding site atoms by the number of frames in the time window (N frames ). These atoms for the respective binding sites were: the carboxylate oxygen atoms of the α-position residue of the SF for binding site A, the carbonyl oxygen atoms of the βand γ-position residues of the SF for binding site B, and the terminal carbon atoms of the hydrocarbon side chain of the isoleucine of the hydrophobic lower gate (I575 in TRPV5) and the amide oxygen atoms of neighboring asparagine residue (N572 in TRPV5). A cut-off distance of 3.5Å was chosen based on the maximum reported distance for calcium-oxygen interactions (Nayal and Di Cera, 1994;Dudev and Lim, 2003;Deng et al., 2006). The mean ion occupancies and standard error were calculated from non-overlapping 50 ns windows of the fivefold replicated 250 ns simulation trajectories with a 20 ps time step.
The ion residence times (t r ) were calculated by averaging the amount of time an individual ion was located within 3.5Å of the center of geometry of the ion coordinating binding site atoms. The mean t r and standard error were calculated from fivefold replicated 250 ns simulation trajectories with a 20 ps time step.
Characterizing permeation cooperativity through mutual information using state-specific information (SSI) from PENSA To characterize the level of cooperativity in the knock-on permeation mechanisms in TRPV channels, we used PENSA to calculate the SSI shared between discrete state transitions in the occupancy distributions of each binding site (Thomson et al., 2021 Preprint;Vögele et al., 2021;Vögele et al., 2022). A timeseries distribution with a time step of 20 ps for each binding site was obtained, whereby for each frame, if an ion occupied the binding site, then this ion's atom ID number was recorded, whereas if the binding site was unoccupied, an ID of −1 was recorded. The ID numbers were discrete, and changes between ID numbers in each binding site, therefore, represent discrete state transitions. By quantifying the mutual information shared between changes to the ID numbers in each site, we were able to determine whether ion transitions at one site were coupled to transitions at another during a 20 ps time interval. From this, we concluded whether cations are "knocking" each other, or dissociation occurred independently from one another. The time interval was iteratively optimized to keep noise and finite sampling effects to a minimum (see below). We found that both were smallest when we used an interval of 20 ps. Similar to McClendon et al. (2009), we observed that finite sampling resulted in independent distributions sharing mutual information (McClendon et al., 2009;Pethel and Hahs, 2014). To overcome this, we calculated a statistical threshold for each simulation via randomly permuted copies of the original data. Random permutations of the original data maintained marginal probabilities for binding site occupation in each simulation while at the same time quantifying the effect of finite sampling on the measurement of SSI. SSI was then calculated between two independently permuted versions of the occupancy distribution for the minimum entropy binding site. Since the upper bound of mutual information between two variables is equal to the lowest entropy of those variables, we used the binding site corresponding to the lowest entropy for obtaining the threshold. This ensured that the portion of SSI which could be attributed to random noise between any two binding sites was always less than or equal to the SSI. This measurement was repeated 1,000 times to resolve a Gaussian distribution from which we obtained the 99% confidence threshold. We subtracted this threshold from the measured values to resolve excess mutual information, or excess SSI (exSSI), shared in discrete state transitions. As it is not possible to transfer negative information, negative exSSI values were corrected to a value of 0.
We also derived a maximum SSI value representing a theoretical upper limit for the information that can be shared between two binding sites, where exSSI max is given by subtracting the random threshold from the minimum entropy of the two binding sites in question.

exSSI(A, B) max min(H(A), H(B)) − threshold(A, B). (2)
To quantify the interdependence of all three ion binding sites within the TRPV pores, the total correlation (TotCorr) was obtained using Eq. 3, where H(A), H(B), and H(C) represent the entropy of binding sites A, B, and C, respectively, and H(A, B, C) the joint entropy of binding sites A, B, and C.
Characterizing the architecture of the SF of TRPV channels To determine the area formed between residues in the SF, the area of the quadrilateral between the adjacent chains was calculated on the X and Y axes. For this, we used the carboxylate oxygen atoms of the adjacent chains for the α-position residue, and the carbonyl oxygen atoms for the β-, γ-, and δ-position residues. To quantify the SF areas from our MD simulations, the mean and standard error were calculated from nonoverlapping 50 ns windows of the fivefold replicated 250 ns trajectories with a 200 ps time step. Furthermore, the SF areas of all TRPV structures deposited in the Protein Data Bank (PDB), available as of February 4, 2022, were determined. A total of 101 structures were analyzed, with non-tetrameric structures or structures without all the atoms of interest modeled not included. The mean and SEM were calculated for all the available structures for any particular TRPV channel.
Online supplemental material Fig. S1 shows voltage-dependence of occupancy and ion selectivity. Fig. S2 shows additional cation binding sites in TRPV5. Fig. S3 shows pore radii and hydrophobicity in simulations of TRPV2, TRPV3, TRPV5, and TRPV6. Fig. S4 shows occupancy and residence times at postulated W583 ion binding site. Fig. S5 shows permeation traces at low ion concentration. Fig. S6 shows permeation traces in low voltage simulations. Fig. S7 shows impact of binding site affinity differences on exSSI. Fig. S8 shows total correlation of cation permeation between all binding sites in TRPV2, TRPV3, TRPV5, and TRPV6. Fig. S9 shows SF backbone dihedral distributions. Fig. S10 shows multiple sequence alignment of TRPV SF sequences. Fig. S11 shows effect of PI(4,5)P 2 on pore radii. Tables S1, S2, and S3 show MD simulation details. Table  S4 shows permeation times. Table S5 shows SSI details. Table S6 shows selectivity in dicationic solutions. Table S7 shows rootmean-square fluctuation (RMSF) of the SF backbone. Supplemental text at the end of the PDF provides additional information. In TRPV5 and TRPV6, Ca 2+ ions traversed the entire pore within average time spans of 28.4 ± 3.9 ns (TRPV5) and 12.0 ± 1.0 ns (TRPV6; Table S4). The calculated Ca 2+ and Na + conductances from our simulations are shown in Table 3. The considerable Na + conductances we observed agree with the experimental finding that the highly Ca 2+ -selective TRPV channels conduct Na + well in the absence of Ca 2+ Vennekens et al., 2000;Voets et al., 2003;Bödding and Flockerzi, 2004). Notably, these conductances are in quantitative agreement with published values measured for Na + in vitro (Yue et al., 2001;Cha et al., 2007). By contrast, control simulations of TRPV5 using the default CHARMM36m force field parameters for Ca 2+ , but otherwise identical conditions, did not exhibit ion permeation; instead, the Ca 2+ ions remained tightly bound to the protein ion binding sites for the entire course of the simulations. This observation is reflective of the shortcomings of standard parameters for divalent cations in fixed-point charge force fields and highlights the improved accuracy of multisite Ca 2+ models in simulating divalent cation permeation and reproducing in vitro conductances. In addition, no Cl − anions were observed to permeate TRPV channels in any of our simulations.

Continuous permeation of Ca
We note that the P Ca /P Na values obtained from our simulations overall show lower Ca 2+ selectivity than the reported literature values (Table S6). We surmised that this might be, at least partially, due to the higher voltages used in our simulations to enhance the sampling rate. Supplementary simulations performed at a lower voltage demonstrated that, indeed, the selectivity for Ca 2+ increases with lower voltages across the membrane (Fig. S1). Below a voltage threshold of ∼205 mV, however, the sampling of permeation events in the simulations became very poor, such that we were not able to reliably probe the precise voltage range used in the experiments.
We can, additionally, not rule out the contribution of force field inaccuracies. Our simulations with Na + show a remarkable agreement between experimentally recorded and simulated conductance. Even though the Ca 2+ model we used has been carefully parameterized (Zhang et al., 2020), modeling divalent cations is a far from trivial task, and this is the first multisite Ca 2+ model with which simulations of ion channel current have become possible. It can therefore not be excluded that further iterations of model refinement may eventually be required to not only reflect experimental solvation-free energies and protein affinity (Zhang et al., 2020) but also accurately reproduce experimental conductance values. This includes kinetic factors such as the correct characteristics of ion (electro-)diffusion in bulk solvent.
Pore cation binding sites and their preference for Ca 2+ binding Prior to the determination of the atomic structures of Ca 2+ channels and the development of channel-permeable models for Ca 2+ ions, it had been suggested from experimental observations that Ca 2+ channels may obtain their selectivity through competition, i.e., by divalent cations, such as Ca 2+ , binding more tightly to their ion binding sites than monovalent cations, such as Na + (Corry et al., 2001;Hille, 2001).
By analyzing the individual traces of permeating Ca 2+ cations along the pore axis z of TRPV5 and TRPV6 over time, we identified three cation-binding sites inside the channels (Fig. 2). We refer to these cation binding sites as sites A, B, and C, viewed from the extracellular entrance of the channel SF to the hydrophobic lower gate. The three cation binding sites were further confirmed by 3-D density analysis using PENSA (Vögele et al., 2021;Vögele et al., 2022;Fig. S2). The PENSA analysis also identified a number of cation-binding sites outside of the pore within the extracellular loops of both TRPV5 and TRPV6 (Fig. S2), in line with the previous suggestion that TRPV6 contains negatively charged "recruitment sites" that funnel cations toward the entrance of the pore (Saotome et al., 2016;Sakipov et al., 2018).
Of the three binding sites we observed, binding site A is formed by the carboxylate oxygen atoms of the ring of acidic residues at the SF entrance (referred to here as the SF α-position); binding site B is formed by the carbonyl oxygen atoms of the bottom two SF residues (SF γand δ-positions); and binding site C is formed jointly by the hydrophobic gate consisting of a ring of isoleucine residues (I575 in TRPV5) and the amide oxygen atoms of the neighboring asparagine residues (N572 in TRPV5) near the cytoplasmic exit of the pore (Fig. 2). The location of these binding sites coincides with constrictions in the pore profile, as determined using CHAP (Klesse et al, 2019;Fig. S3; see Fig. S11). The distance between binding sites A and B is ∼5Å, and that between binding sites B and C is ∼14Å. We note that Hughes et al. (2018) reported a further constriction below the hydrophobic gate (binding site C) formed by W583 in the TRPV5 structure, and an analogous constriction at W583 can be observed in TRPV6. However, our simulations do not suggest that the side chains of W583 constitute a functionally important ion binding site, as shown in Fig. S4.
In monocationic Ca 2+ solutions, the three binding sites showed Ca 2+ occupancy probabilities of 0.69 ± 0.05, 0.67 ± 0.04, and 0.57 ± 0.03 in TRPV5, and 0.43 ± 0.04, 0.54 ± 0.05, and 0.29 ± 0.02 in TRPV6 (from A to C, respectively; Fig. 3). In monocationic Na + solutions, similar occupancies were observed. However, the Na + residence times (t r ) at the three binding sites were markedly lower than those observed for Ca 2+ , with ratios of t r (Ca 2+ ):t r (Na + ) varying between ∼35:1 and ∼3:1 (Fig. 3). These residence times suggest that Ca 2+ ions have a greater affinity for these binding sites than Na + . This observation was further substantiated when the occupancy of binding sites in dicationic solutions was analyzed, in which Ca 2+ and Na + cations are competing for the binding sites. Binding sites A, B, and C in the Ca 2+ -selective channels all showed high occupancies with Ca 2+ in the mixed cationic solutions (Fig. 3). Across all the replicate simulations, we recorded average Ca 2+ occupancies of 0.51 ± 0.05, 0.48 ± 0.06, and 0.37 ± 0.03 for binding sites A, B, and C in TRPV5, respectively, and of 0.68 ± 0.04, 0.61 ± 0.04, and 0.40 ± 0.04 for binding sites A, B, and C in TRPV6 (Fig. 3). By contrast, the Na + occupancy of each of the three binding sites under these conditions was found to be below 0.07, both in TRPV5 and TRPV6; that is, the ratio between Ca 2+ and Na + occupancy varies between ∼85:1 and ∼7:1 (Fig. 3). These values indicate a free energy difference of between 11.5 and 5.0 kJ mol −1 for the preferential binding of Ca 2+ over Na + . TRPV5 53 ± 7 49 ± 6 92 ± 10 59 ± 6 (Na + ) (Cha et al., 2007) TRPV6 117 ± 12 61 ± 6 29 ± 6 58 ± 4 (Na + ) (Yue et al., 2001) Mean inward conductances and SEM were calculated from overlapping 50 ns windows from fivefold replicated 250 ns simulations. Permeation of monocationic Ca 2+ or Na + cations was simulated in a solution of 150 mM CaCl 2 or 150 mM NaCl, respectively. Permeation of a dicationic mixture of Ca 2+ and Na + cations was investigated in a solution of 75 mM CaCl 2 and 75 mM NaCl.
A highly cooperative knock-on mechanism between three cation binding sites underpins selective Ca 2+ permeation in TRPV channels The observed increased affinity for Ca 2+ cations at the pore binding sites compared to Na + means that, in a mixed solution, Ca 2+ will preferentially occupy these binding sites; however, this also implies that Ca 2+ ions face a greater energy barrier when they dissociate from the binding sites. In monocationic solutions, this would result in a greatly reduced Ca 2+ conductance with respect to Na + . For instance, based on our observed residence times in monocationic solutions, we would expect an ∼12fold reduced Ca 2+ unbinding rate compared to Na + for binding site A in TRPV5. However, a much-reduced Ca 2+ conductance is neither observed in our simulations nor the experimental  literature. Due to the divalent charge of Ca 2+ , increasing the affinity to cation binding sites, this dichotomy had previously been suggested to exist, and it was hypothesized that this paradox could be resolved by assuming cooperativity between successive unbinding events such as in a knock-on mechanism (Corry et al., 2001;Hille, 2001).
In the classic knock-on mechanism, which for example underpins K + channel function, ions transition into and out of multiple ion-binding sites in a highly correlated fashion (Armstrong, 1971;Hille and Schwarz, 1978;Neyton and Miller, 1988;Morais-Cabral et al., 2001;Köpfer et al., 2014). For example, early experiments by Hodgkin and Keynes and later fluxratio measurements established that 3-3.4 K + ions moved in lockstep with each other during permeation through K + channels (Hodgkin and Keynes, 1955;Stampe and Begenisich, 1996).
For each permeating ion in a simulation of TRPV5, Fig. 4 shows the association and dissociation of Ca 2+ and Na + ions at binding sites A, B, and C from top to bottom as color code (bound to A, red; bound to B, orange; bound to C, yellow; transiting within the pore but not bound to a binding site, blue; located in extracellular solvent, dark gray; and located in intracellular solvent, light gray). As can be seen for Ca 2+ in TRPV5 for example (Fig. 4 left), the plot demonstrates that (i) permeating Ca 2+ ions spend the vast majority of their time within the pore at the three binding sites (reflected in the scarcity of blue boxes vs. red, orange, and yellow), (ii) dual and triple occupancy of the three sites, A, B, and C, with Ca 2+ frequently observed (horizontal slices across plot: triple occupancy is observed in 27.2% of the simulation frames, dual occupancy in 49.7%), and (iii) transitions between states show a high degree of correlation, i.e., the ions frequently move in concert into and out of their respective binding sites (horizontal slices; binding state transitions). By contrast, during Na + permeation (Fig. 4 center), the ions are predominantly transiting across the pore without occupying particular binding sites for extended time spans (blue, on average 53% of the traversal time for each ion).
To go beyond visual inspection of the trajectories and to assess the cooperativity of ion permeation in a quantitative way, we developed a new approach based on mutual information, taking into account the "state" of each ion binding site. To achieve this, we assigned a specific binding state (unoccupied or occupied with a specific ion) to binding sites A, B, and C and used our recently developed approach, SSI (Thomson et al., 2021 Preprint), on pairs of adjacent sites to quantify the degree of coupling between ion binding transitions at each of these sites (see Materials and methods). This analysis yields a coefficient quantifying the cooperativity between ion binding and unbinding at neighboring or more distant binding sites, where a greater coefficient signifies a higher degree of coupling. This coupling suggests that when an ion transitions from one site it is more likely that there is a transition at the other. To correct for the non-zero mutual information that samples of two completely independent variables can display due to finite-size effects, we followed the approach of McClendon et al. (2009) and Pethel and Hahs (2014) to yield excess mutual information, or exSSI. We also determined a theoretical upper limit for the maximum mutual information that can be shared between two binding sites by using the minimum state entropy among the two sites. Note that this quantity represents an absolute upper limit; reaching it would require both binding sites to exhibit idealized . Permeation state plots of permeating Ca 2+ (left) and Na + (center) cations through the Ca 2+ -selective TRPV5. Permeation state plots show the state of each permeating cation (columns) at a given time point (rows) by assigning a state to the ions to indicate whether the cation is bound to a binding site, transitioning between binding sites, or in the bulk solution: bound to A, red; bound to B, orange; bound to C, yellow; transiting within the pore but not bound to a binding site, blue; located inextracellular solvent, dark gray; located in intracellular solvent, light gray. Comparison of permeation state plots for Ca 2+ and Na + cations shows that Ca 2+ permeation proceeds in a well-ordered manner with three Ca 2+ cations within the pore and knocking adjacent cations to the next binding site. Na + permeation, on the other hand, is far less ordered, with regularly more than three Na + cations within the pore at a given time. Each plot (left and center) shows exemplars from a single 250 ns simulation of TRPV5 performed in monocationic 150 mM CaCl 2 or 150 mM NaCl, respectively. The structure of TRPV5 shows the colors used in the permeation state plots and the location of the residues constituting the three binding sites (right). Please note, only cations that fully permeate through the pore within the 250 ns simulation are shown in the plots, whereas the binding site occupancies reported above reflect both permeating and non-permeating ions (left and right).
simultaneous states and state transitions throughout the entire simulated time.
The SSI analysis showed that in the Ca 2+ -selective TRPV channels under the simulated conditions, TRPV5 and TRPV6, a high level of information above noise is shared between the transition of ions into and out of binding sites A and B, respectively, both for the permeation of Ca 2+ and Na + (exSSI between 0.8 and 1.6 bits; Table S5 and Fig. 5). This suggests that the ion binding and unbinding processes at each of these binding sites are coupled to one another, constituting a knockon mechanism at relatively short range. We observed three to four water molecules on average between cations bound at binding sites A and B during knock-on, demonstrating a "soft" knock-on mechanism to be in place, as opposed to the "direct" knock-on mechanism between dehydrated K + ions in K + -selective cation channels (Köpfer et al., 2014). As detailed further below, our simulations indicate that only a moderate level of ion desolvation occurs in the SF of the studied TRPV channels, such that the hydration shell of the permeating ions remains largely intact during knock-on.
Similarly, the transitions of ions into and out of binding sites B and C show a large degree of correlation for both Ca 2+ and Na + (Fig. 5). In the case of binding sites B and C, however, this requires a remote knock-on mechanism to be in operation, since these sites are ∼14Å apart. The concept of a remote knock-on event was first proposed by Tindjong et al. (2013) based upon Brownian dynamics simulations and observed by Zhang et al. (2020) in atomistic MD simulations of Ca 2+ permeation in the RyR1 channel. Our SSI analysis suggests that the degree of cooperativity in the remote knock-on mechanism between binding sites B and C (Fig. 5 right) is slightly smaller than the cooperativity in the closer knock-on mechanism between binding sites A and B (Fig. 5 left). Simulations conducted at a lower CaCl 2 concentration (Fig. S5) and reduced transmembrane voltage (Fig.  S6) confirmed that the knock-on mechanism observed occurs also under these milder conditions; however, we did not achieve enough sampling to recalculate the exSSI values at these conditions.

Cation permeation in non-selective TRPV channels shows a lower degree of cooperativity
To determine if the remaining, non-selective TRPV channels showed a different permeation mechanism, we next performed simulations of the open-state TRPV2 (Dosey et al., 2019) and TRPV3 (Singh et al., 2019) channels using the same simulation approach as described for the Ca 2+ -selective TRPV channels. These simulations of non-selective TRPV channels also showed continuous ion permeation, with cation conductances, again, in good agreement with published conductance values measured in vitro (Table 4). Overall, we recorded 706 complete inward channel crossings for Ca 2+ and 1,176 for Na + from simulations of the non-selective TRPV channels.
The occupancy of binding sites B and C in the TRPV2 and TRPV3 systems showed no clear difference to the Ca 2+ -selective channels. By contrast, the occupancy of binding site A was reduced by ∼50% for both Na + and Ca 2+ ions in the monocationic solutions (Fig. 3). This suggests that cations are less well coordinated at binding site A in the non-selective TRPV channels, leading to lower affinity binding in the monocationic solutions. All binding sites, however, exhibited a preference for binding Ca 2+ in the dicationic solutions. Note that, in the non-selective channels TRPV2 and TRPV3, this is coupled with a particularly low residence time for Ca 2+ at binding site A, again suggesting higher exchange rates and weaker binding, despite its occupancy. That is, many Ca 2+ ions are observed to diffuse back from A to bulk solution in TRPV2 and TRPV3.
We were therefore curious if the three-site knock-on mechanism described previously for Ca 2+ -selective TRPV channels is also at play in the non-selective TRPV channels. Our SSI analysis confirmed that the cooperativity between binding sites B and C Figure 5. exSSI between ion binding sites quantifies the degree of cooperativity in the knock-on mechanism of cation permeation in TRPV channels. The mean exSSI and SEM between transitions from binding sites A and B (left) and binding sites B and C (right) are shown for Ca 2+ (orange) and Na + (blue) cations in monocationic solutions. For each exSSI, the mean maximum exSSI max and standard error are also shown (gray), and the exSSI norm is reported in Table S5.
in the non-selective TRPV channels was comparable with those calculated for the Ca 2+ -selective TRPV channels (Table S5 and Fig. 5). However, the correlation between ion binding transitions at binding sites A and B was substantially reduced in both of the non-selective TRPV systems (Fig. 5). The nearly complete absence of cooperativity from binding sites A and B demonstrates that a knock-on mechanism is not occurring between these two sites in the non-selective TRPV channels. Instead, our findings suggest that cation permeation in the non-selective TRPV channels occurs via a two-site knock-on mechanism between binding sites B and C.
Since the ion occupancy observed at binding site A is reduced in the case of the non-selective TRPV channels, it is plausible that this lower affinity also impacts the coupling between transitions at binding sites A and B. To test this notion further, the relationship between affinity differences of a pair of binding sites and the knock-on co-operativity was tested systematically by using a toy model with two energy wells (binding sites) possessing a range of different depths (affinities). As shown in Fig. S7, there is a linear relationship between the affinity difference and the observed exSSI. This demonstrates that the diminished affinity of binding site A is likely to be the major reason for the loss of cooperativity in the SF of the non-selective TRPV channels. Our results show that similar binding affinity is a necessary but not sufficient condition for a high degree of cooperativity between two cation binding sites.
Based on our SSI approach to quantify mutual information in binding and unbinding events at different ion binding sites and using the concept of total correlation to evaluate the overall cooperativity in a system across all coupled events, we next calculated the total correlation of ion permeation for all the TRP channels investigated. The reduction in the number of correlated knock-on sites within the non-selective TRP channels can be distinctly seen when comparing the total correlation for each of the non-selective and Ca 2+ -selective TRPV channels (Fig. S8).
The ion binding sites in non-selective TRPV channels preferentially bind Ca 2+ over Na + , as in Ca 2+ -selective TRPV channels, which explains why the non-selective TRPV channels in fact show slight Ca 2+ selectivity (P Ca /P Na ∼10:1; Owsianik et al., 2006). However, our data suggest that the reduced level of coordination at binding site A, and especially the effect on the three-site cooperativity it imparts together with sites B and C, reduces the Ca 2+ selectivity from P Ca /P Na ∼100:1 seen in TRPV5 and TRPV6, and in this way ultimately determines the difference between Ca 2+ -selective and non-selective permeation.

Structural features distinguishing Ca 2+ -selective from nonselective permeation
To understand why cation coordination at binding site A is weakened in non-selective TRPV channels, decoupling its cooperativity, we investigated the area formed between the four subunits of the TRPV channel for each residue in the selectivity filter. The cross-sectional area formed by the carboxylate oxygen atoms at the SF α-position is larger in the non-selective TRPV channels than in the Ca 2+ -selective TRPV channels (Fig. 6). Interestingly, despite the increased average area formed by the carboxylate oxygen atoms at the α-position residue in the selectivity filter, the average area formed by the carbonyl oxygen atoms at both the γand δ-positions are smaller in non-selective TRPV channels than in Ca 2+ -selective TRPV channels (Fig. 6). These differences in selectivity architecture were further confirmed using the pore profile calculated using CHAP (Klesse et al, 2019;Fig. S3). Our simulations showed no differences in the flexibility of selectivity filters between Ca 2+ -selective and nonselective TRPV channels as determined by RMSF calculations of the backbone atoms (Table S7); we did however observe a small difference in the backbone dihedral angle distribution of the β-position residue, which we predict to be due to differences in stabilizing hydrophobic interactions of the residue side-chain (Fig. S9, S10, and S11).
To expand the geometric analysis to all available TRPV channel structures, we also calculated the average area formed by selectivity filter residues from all available TRPV structures deposited in the PDB (Fig. 6). Analysis of these static structures showed the same two trends observed within our MD simulations: (1) Ca 2+ -selective TRPV channels clearly have a smaller average area formed by the carboxylate oxygen atoms of the SF α-position; and (2) non-selective TRPV channels have a slightly narrower constriction at the SF γand δ-position residues. Our MD simulations suggested that the wider opening at the α-position leads to a weaker cation binding interaction at binding site A and cation coordination, which, in turn, decouples binding site A from the co-operative knock-on mechanism that underpins permeation in the Ca 2+ -selective TRPV channels.
The non-selective channels possess a narrower constriction formed by the side chains and carbonyl groups of the γand δ-position residues forming binding site B, whereas the Ca 2+selective TRPV channels exhibit an even narrower constriction at binding site A. Our simulations suggest that the greater occupancy of binding site A is a result of the more confined geometry of this charged binding site in the Ca 2+ -selective TRPV channels. This, in turn, leads to a higher occupancy of the uncharged binding site B in the Ca 2+ -selective TRPV channels (see Fig. 3), as binding site B receives cations from the adjacent binding site A via the knock-on mechanism, rather than having to bind them from bulk solution, as would be the case in the nonselective channels. Thus, despite the constriction not being as narrow, binding site B has a higher occupancy in the Ca 2+ -selective than in the non-selective TRPV channels.  (Zhang et al., 2016) TRPV3 196 ± 10 299 ± 10 222 ± 13 201 (Na + ) (Chung et al., 2004) Mean inward conductances and SEM were calculated from overlapping 50 ns windows from fivefold replicated 250 ns simulations. Monocationic solutions at 150 mM concentration; dicationic mixture of Ca 2+ and Na + at a concentration of 75 mM each.

Ives et al. Journal of General Physiology
Ca 2+ -selective permeation is not strongly linked to the solvation states of permeating cations A previously reported mechanism of how cation selectivity can be achieved is by desolvation of permeating cations. Differences in desolvation energies between cationic species provide a thermodynamic penalty that can be more favorable for the permeation of one cationic species over another. Such a mechanism has been reported to underpin K + selectivity over Na + in K + channels (Köpfer et al., 2014;Kopec et al., 2018), for example. Alternative mechanisms suggested to yield K + selectivity are based on protein flexibility, especially the plasticity of the K + channel SF (Noskov et al., 2004), the number of stacked ion binding sites in the SF (Derebe et al., 2011), where a reduction from four to three has been shown to abolish K + selectivity, the kinetic model of selectivity (Thompson et al., 2009;Kim and Allen, 2011), and the coordination model (Bostick and Brooks, 2007;Varma et al., 2008). By contrast, Na + selectivity has been suggested to rely chiefly on a "snug fit" coordination of the Na + ion in the SFs of eukaryotic Na V channels and the preservation of its solvation shell while permeating bacterial Na V channels (Dudev and Lim, 2014).
To investigate whether a high degree of desolvation was a dominant factor in Ca 2+ selectivity in TRPV channels, we determined the number of oxygen atoms within a 3Å radius of the cations, representing their first solvation shell (Fig. 7). In the bulk solution, both Ca 2+ and Na + cations showed their expected water coordination number of 7 and 5.6, respectively. As the cations entered the pore, we saw a small degree of partial dehydration of permeating cations at the SF (Fig. 7). In particular, the carboxylate oxygen atoms of the acidic residue at the entrance of the SF coordinated an incoming cation, with these displacing up to approximately two coordinated water molecules from the first solvation shell of the cation. We also observe some desolvation around, or below, binding site C. However, no major differences in the solvation shell of permeating Ca 2+ or Na + cations, or indeed between the Ca 2+ -selective and non-selective TRPV channels were observed. The finding that the degree of desolvation does not differ substantially between Ca 2+ -selective and non-selective channels indicates that their selectivity is not exclusively based on a mechanism of ion dehydration. The largest level of desolvation is seen at ion binding site A in the TRPV5 channel for Ca 2+ . Since the overall desolvation penalty for Ca 2+ is larger than for Na + (Marcus, 1991), one would expect this to lead to the preferential binding of Na + at this site, whereas the opposite is observed (Fig. 3), further arguing against dehydration as a major selectivity mechanism in TRPV channels.

Discussion
Our simulations showed three main cation binding sites in the permeation pathway of TRPV channels, which we term binding sites A, B, and C. Binding sites A and B are formed by the carboxylate oxygen atoms of the α-position residue of the SF and by the carbonyl oxygen atoms of the γand δ-position residues of the SF, respectively. Binding site C is located just above the hydrophobic lower gate of the pore and is formed by the isoleucine residues of the lower gate and amide oxygen atoms of the neighboring asparagine residues.
The identification of these cation-binding sites in our MD simulations is in agreement with previously published TRPV structures and other MD simulations. In their crystal structure Figure 6. Architecture of the four-residue selectivity filter of the investigated TRPV channels. The area formed between the carboxylate oxygen atoms (α-position) and backbone carbonyl atoms (β-, γ-, and δpositions) was calculated for all Ca 2+ -selective and non-selective TRPV systems simulated as part of this study (left), and for all TRPV structures available in the PDB at the time of writing (right). The SF structures of the Ca 2+ -selective TRPV5 and non-selective TRPV3 are shown for reference (center). The mean area between SF residues and SEM was calculated from non-overlapping 50 ns windows from fivefold replicated 250 ns simulations in 150 mM CaCl 2 . The mean area between SF residues and SEM from PDB structures was calculated from all available homotetrameric TRPV structures at the time of writing.
of the Rattus norvegicus TRPV6, Saotome et al. (2016) identified three regions of electron density within the channel pore which they interpreted as cation-binding sites. It should be noted that there are some minor differences in the residues forming binding site C to the binding site reported here, likely due to the structure of TRPV6 of Saotome et al. (2016)  Moreover, numerous structures of TRPV channels deposited within the PDB are resolved with cations bound at one of the cation-binding sites identified in our MD simulations. These include structures from several orthologs of the TRPV channels simulated in this study, as well as of the TRPV1 and TRPV4 channels that were not simulated in this study. For example, structures of the closed-state TRPV1 channel of R. norvegicus (Kwon et al., 2021) and the open-state TRPV4 channel of Homo sapiens in complex with 4α-PDD (Botte et al., 2020 Preprint), both model cations bound at binding site B. This observation further validates the existence of the identified cation binding sites, as well as their conservation among TRPV channels and perhaps the wider TRP superfamily.
In addition to structural data, Sakipov et al. (2018) performed equilibrium MD simulations of Ca 2+ movement in the closedstate structure of R. norvegicus TRPV6 and identified cationbinding sites in the SF (Sakipov et al., 2018). The binding sites from this work are also generally in agreement with our simulation results and the aforementioned structural data. However, Sakipov et al. (2018) reported that their simulation data identified two Ca 2+ cations residing at binding site A, with adjacent chains each occupying an ion. In accordance with previous x-ray crystallographic data, two ions associated with binding site A are not seen in a major population of our simulation ensembles, which we attribute to the use of an optimized multisite model for Ca 2+ ions in the present study (Zhang et al., 2020) and increased sampling for improved statistical analysis.
In monocationic solutions, we observed a high probability of at least two of the three binding sites (A, B, and C) being simultaneously occupied with either Na + or Ca 2+ , respectively, with mostly insignificant differences between the occupancy values for Na + and Ca 2+ at each individual binding site. However, the residence times were markedly reduced at all sites for Na + ions, and overall a tendency toward weaker binding for either ion at binding site A at the extracellular SF entrance was observed. In mixed, dicationic solutions of Na + and Ca 2+ , by contrast, Ca 2+ ions strongly outcompeted Na + ions for association at all pore-binding sites, both in the Ca 2+ -selective and nonselective TRPV channels. Therefore, according to our data, all the binding sites display a much greater affinity for Ca 2+ . This gave rise to the question of how permeation efficiency for conducting Ca 2+ ions is achieved in these channels and how this relates to their varying degrees of selectivity since higher affinity binding is usually expected to lead to slower permeation rates.
By ensuring cooperativity between binding/unbinding events at multiple binding sites, permeation rates can be enhanced (Corry et al., 2001;Hille, 2001). We hypothesized that the level of cooperativity between sites A, B, and C in the pore could underpin the difference between highly and less Ca 2+ -selective TRPV channels. We, therefore, developed a novel method to quantify co-operativity during ion permeation in pores with multiple ion binding sites based on mutual information and total correlation measures using the SSI approach (Thomson et al., 2021 Preprint). We anticipate that this method will be similarly useful for the study of permeation mechanisms and the basis of selectivity in other channels. Our analysis showed that there is a substantial degree of co-operativity for Ca 2+ permeation between binding sites B and C across all TRPV channels, whereas a clear distinction exists between the co-operativity between binding sites A and B in Ca 2+ -selective and non-selective TRPV channels. In the non-selective TRPV channels, binding site A is decoupled from binding site B. By contrast, binding sites A and B are even more strongly coupled than binding sites B and C in the case of the Ca 2+ -selective channels. We suggest that this marked difference in co-operativity mechanistically explains the different levels of Ca 2+ selectivity in TRPV channels.

Conclusion
We have characterized the cation permeation mechanisms in four members of the TRPV channel family. We identified three cation binding sites within the pore, each of which displayed greater affinity for Ca 2+ binding over Na + . A novel application of mutual information between ion binding and unbinding at consecutive binding sites (SSI) enabled us to quantify the degree of knock-on taking place in ion permeation. This showed that the level of Ca 2+ selectivity in TRPV channels is determined by the co-operativity, or coupling, between the transitions at three pore cation-binding sites. The Ca 2+ -selective TRPV channels display a highly correlated three-site knock-on, whereas the cation binding site at the extracellular entrance is decoupled from the mechanism in the non-selective TRPV channels, which reduces the overall preference for Ca 2+ permeation.
Data availability MD simulation inputs and analysis scripts used for this study are deposited in a public GitHub repository, available at: https:// github.com/cmives/Ca_selectivity_mechanism_of_TRPV_ channels.

Supplemental material
Summary of MD simulations within this study Tables S1, S2, and S3 show details of all MD simulations performed in this study.
Average permeation times for Ca 2+ and Na + cations in MD simulations of TRPV channels The average permeation time of cations from monocationic simulations is described in Table S4. The average permeation time was defined as the time taken between the cation binding at binding site A, and dissociating from binding site C into the bulk solution. This data showed that for every TRPV system, the average permeation time was greater for Ca 2+ permeation than for Na + , with the exception of TRPV6 where there was no discernible difference between the average permeation times. This greater permeation time of Ca 2+ cations is likely a result of their greater affinity for the cation binding sites. However, the difference in permeation time is less than one would expect based on the residency times (Fig. 3). For example, the average time for a Ca 2+ cation to permeate through TRPV5 is ∼1.5-fold slower than Na + permeation, whereas the residency time of cations at binding site A of TRPV5 would suggest a 12fold difference. This supports the concept of cooperativity between successive binding sites increasing the unbinding rate of ions from binding sites, as proposed by Hille (2001).

The effect of cation concentration and transmembrane on the knock-on permeation mechanism
In vitro electrophysiology experiments of TRPV5 and TRPV6 in the literature have been conducted with a maximum Ca 2+ concentration of 30 mM Hoenderup et al., 2001) and concentrations of up to 30 mM Ca 2+ were used in studies of other TRPV channels (Watanabe et al., 2003). Moreover, many electrophysiological recordings on TRPV channels were performed at Na + concentrations of 150 mM . We therefore believe that the concentrations we probed in our simulations do not deviate too strongly from realistic concentrations used under experimental conditions. However, to confirm that the cooperativity and knock-on mechanism observed was not a consequence of the higher cation concentrations, we also performed simulations of TRPV5 in a lower Ca 2+ concentration of 25 mM CaCl 2 , rather than 150 mM CaCl 2 , which is well within the range of concentrations used in in vitro experiments. As can be seen in Fig. S5, the knock-on mechanism between Ca 2+ ions is preserved at a lower Ca 2+ concentration. Specifically, two three-site knock-on permeation events can be seen at ∼50 and ∼78 ns.
Pore architecture of Ca 2+ -selective and non-selective TRPV channels from MD simulations To investigate the time-averaged architecture of TRPV channels in our simulations, we used CHAP (Klesse et al., 2019) to create profiles of the pore radius and the hydrophobicity of pore-facing residues (Fig. S3). The CHAP analysis of the simulated TRPV structures showed similar profiles to the respective starting, static structures. These include an upper constriction formed by the SF and a lower constriction formed by a hydrophobic gate.
Identification of cation binding sites using PENSA In addition to permeation traces (Fig. 2), cation binding sites were further identified using the Featurizer function of PENSA (Vögele et al., 2021). This analysis identified cation binding sites from 3-D density maxima in the simulation. As can be seen in Fig. S2, this analysis identified binding sites A, B, and C, as well as cation binding sites within the extracellular loops of the protein. This finding is in agreement with previous suggestions of recruitment sites in TRPV6: negatively charged or polar residues that attract cations and funnel them toward the entrance of the pore (Saotome et al., 2016;Sakipov et al., 2018). In particular, our MD simulations suggested that the following residues act as recruitment sites: E522, D525, T528, F531, S532, E535, Y547, and Y549.

W583 does not constitute a cation-binding site in our simulations
In the structure of TRPV5, Hughes et al. reported a constriction below the hydrophobic lower gate formed by residue W583 (Hughes et al., 2018). An analogous constriction formed by W583 can be observed in the structure of TRPV6 determined by McGoldrick et al. (2018). Therefore, we investigated whether this constriction also constituted a cation binding site in our simulations (putatively referred to as binding site D), in addition to binding sites A, B, and C. Analysis of simulations of TRPV5 in a monocationic solution of 150 mM CaCl 2 showed that the occupancy probability and residence time of Ca 2+ at binding site D is substantially lower than at the other binding sites (Fig. S4). Therefore, we did not consider W583 to constitute a functionally important cation-binding site. To evaluate and quantify the overall cooperativity in a system across all coupled binding site transition events, we used the concept of total correlation. This total correlation analysis showed that the Ca 2+ -selective TRPV channels had a greater total correlation than the non-selective TRPV channels (Fig. S8). The greater total correlation in Ca 2+ -selective TRPV channels is due to these channels consisting of a three-binding site knock-on permeation mechanism. However, non-selective TRPV channels consist of a two-binding site knock-on permeation mechanism due to the loss of cooperativity between binding sites A and B, reducing the total correlation of these channels.

SSI and total correlation of cation binding site transitions in cation permeation
The effect of binding site affinity on the knock-on permeation mechanism To investigate if it is a necessary requirement for the knock-on mechanism for binding sites to have a similar binding affinity, we analyzed the exSSI between two binding sites in a toy model. This toy model consisted of binding site occupancy data generated for two binding sites, binding sites ε and ζ, for the equivalent of a 100-ns simulation with the ion occupation state of each binding site reported every 20 ps, as in the analysis of real simulation trajectories. In this simple model, the occupancy probability of binding site ε was 0.9 and the occupancy probability of binding site ζ was adjusted to investigate the effect of differing binding site affinities. The ratio of occupancy probabilities between ε:ζ was: 1:1, 1:0.75, 1:0.5, 1:0.25, 1:0.1, 1:0.01. exSSI analysis of this system showed a clear, linear, and positive correlation between the ratio of binding site affinities and the exSSI shared between them (Fig. S7). These results show that similar binding affinity is a necessary, but not sufficient condition for a high degree of cooperativity between two cation binding sites, since despite similar affinity, coupling could be absent.
Lower voltage simulations of TRPV5 in a dicationic solution As described previously, TRPV5 and TRPV6 are highly Ca 2+ selective with a reported in vitro experimental permeability ratio (PCa/PNa) of ∼100:1 Yue et al., 2001). The remaining members of the TRPV subfamily, TRPV1-4, are still slightly Ca 2+ -selective however, with a permeability ratio P Ca /P Na of ∼10:1 (Owsianik et al., 2006). Table S6 shows the calculated P Ca /P Na from our MD simulations of TRPV channels in dicationic solutions at a voltage of −410 mV.
These calculated in silico selectivity values are lower than the in vitro selectivity values. We surmised that this may be due to the higher voltages (∼410 mV) used in our simulations to enhance the sampling rate of ion permeation events. To test this, we ran simulations of TRPV5 in a dicationic solution using a reduced voltage of ∼205 mV. These simulations showed an increase in Ca 2+ selectivity, as well as an increase in occupancy probability of Ca 2+ cations in the cation binding sites (Fig. S1). However, these simulations with a lower voltage resulted in a lower number of permeation events, meaning that they were not used to elucidate the mechanism of cation permeation due to reduced sampling.
Furthermore, permeation traces from these simulations showed that a clear and distinct knock-on mechanism between permeating Ca 2+ cations is present in lower voltage simulations. This confirms that the observed knock-on mechanism and cation cooperativity are not a result of the increased voltage used in these simulations as compared with the experiment.
Conformational flexibility in the selectivity filter of Ca 2+ -selective and non-selective TRPV channels In addition to our investigations of the differences in SF architecture between Ca 2+ -selective and non-selective TRPV channels, we also investigated whether there were differences in the flexibility of the SFs. Analysis of the RMSF of the backbone atoms for each residue in the SF showed similar RMSF values for all residues of all TRPV channels (Table S7), with RMSF values ranging between 0.6 and 1.3Å.
In contrast, analysis of the backbone dihedral angles showed some differences in the angle distribution between Ca 2+ and nonselective TRPV channels (Fig. S9). The most obvious difference was in the β-position of the SF, where non-selective TRPV channels showed wider distribution of the φ angle, and principally on the ψ angle.
Comparison of the amino acid composition of the SF of TRPV channels shows that the β-position residue in Ca 2+ -selective TRPV channels is usually an isoleucine residue; however, in non-selective TRPV channels this β-residue is a much smaller glycine residue (Fig. 6). As the side-chain of the β-position residue faces into the hydrophobic pocket between the SF and pore-helix, the hydrophobic side-chain of isoleucine is involved in hydrophobic interactions with this region, increasing the backbone stability of the β-position in Ca 2+ -selective TRPV channels. On the other hand, the glycine residue at the β-position in the non-selective TRPV channels will have greater backbone flexibility. Furthermore, the increased flexibility at this β-position will affect the adjacent αand γ-positions, explaining why these two positions show differences in their backbone dihedral angle distributions between Ca 2+selective and non-selective TRPV channels, but no difference in the δ-position.
The effect of PI(4,5)P2 protonation state on the TRPV5 simulation system The cryo-EM structure utilized of TRPV5 of Oryctolagus cuniculus was solved to a resolution of 4Å and resolved a PI(4,5)P2 molecule bound to each chain of the homotetrameric structure (Hughes et al., 2018). At such a resolution, it is not possible to accurately model the protonation state of the PI(4,5)P2 molecule. Hughes et al. report that PI(4,5)P2 binding induces conformational changes related to channel activity. As our simulation protocol included harmonic restraints on the lower gate of the protein to maintain the open state of the structure, the protonation state of the PI(4,5)P2 would be inconsequential to the stability of the pore during our MD simulations. To confirm this, we performed threefold 150-ns simulations of TRPV5 with deprotonated PI(4,5)P2 molecules (molecule charge = −4) and compared the pore architecture to the production simulations of fivefold 250-ns simulations of TRPV5 with protonated PI(4,5)P2 molecules.
As can be seen in Fig. S11, our simulations show no significant difference in the pore architecture dependent on the protonation state of the PI(4,5)P2 in TRPV5 simulations. This demonstrates that our simulation protocol and use of harmonic restraints on the lower gate reliably constrain the protein structures in their open-state conformations. Figure S1. The effect of high voltage on Ca 2+ -selectivity in simulations of TRPV5 in dicationic solutions. Simulations performed at a lower voltage of ∼205 mV (left) resulted in an increased occupancy probability of Ca 2+ cations (orange) and a reduced occupancy probability of Na + cations (blue) compared to simulations at a higher voltage of ∼410 mV (center). This resulted in increased Ca 2+ selectivity in lower voltage simulations, as summarized in the P Ca /P Na value (right).

Ives et al.
Journal of General Physiology S3 TRPV channel permeation mechanism https://doi.org/10.1085/jgp.202213226 Figure S2. Cation binding sites in TRPV5 identified by PENSA. 3-D density maxima of Ca 2+ cations within 7Å of the protein was analyzed to identify 12 cation binding sites shown as pseudo-atoms (left). This analysis identified binding sites A (red), B (orange), and C (yellow) and several other "recruitment sites" (purple). The location of these recruitment sites (purple) help to attract cations and funnel them toward the pore entrance (top right and top left). Figure S3. Pore architecture of TRPV channels from MD simulations, showing the average radius and hydrophobicity of the channel withrespect to the relative z coordinate, obtained using CHAP (Klesse et al, 2019). The mean radius or hydrophobicity (black) and SD (gray) were calculated from concatenated trajectories of fivefold replicated 250 ns simulations in 150 mM CaCl 2 with a 200 ps time step. The shaded gray region represents the SD. The average position of binding sites A, B, and C are shown as shaded red, orange, and yellow regions, respectively. The dashed line in the pore radius plots indicates the radius of a dehydrated Ca 2+ ion.
Ives et al. Journal of General Physiology S4 TRPV channel permeation mechanism https://doi.org/10.1085/jgp.202213226 Figure S4. W583 does not form a functionally important cation binding site in simulations of TRPV5 in 150 mM CaCl 2 . W583 is located on the S6 helix, below the hydrophobic lower gate (center). Analysis of the occupancy probability (left) and the residence time (right) showed that the constriction formed by W583 does not coordinate Ca 2+ cations as efficiently as binding sites A, B, and C in our simulations.  Figure S7. The effect of similar binding site affinities on exSSI generated from a model of two consecutive binding sites. When the binding site affinity of binding site ε was reduced relative to binding site ζ, the exSSI also decreased in a linear fashion.
Ives et al. Journal of General Physiology S6 TRPV channel permeation mechanism https://doi.org/10.1085/jgp.202213226 Figure S8. Total correlation of cation permeation between cation binding sites from simulations of TRPV channels. Comparison of the total correlation showed that the Ca 2+ -selective TRPV channels have a greater total correlation than non-selective TRPV channels (center; Ca 2+ , orange bars; Na + , blue bars). This greater total correlation in Ca 2+ -selective TRPV channels is a consequence of a knock-on mechanism between three binding sites (left). However, in nonselective TRPV channels, cation coordination at binding site A is reduced, resulting in reduced cooperativity, and a two-binding site knock-on mechanism between binding sites B and C only (right). Provided online are seven tables. Table S1 shows summary of simulation details of Ca 2+ -selective TRPV channels. Table S2 shows the summary of simulation details of non-selective TRPV channels. Table S3 shows the summary of simulation details of additional control simulations of Ca 2+ -selective TRPV channels. Table S4 shows the average time to permeate through the TRPV pore, as defined by the z position between binding sites A and C. Table S5 shows calculated exSSI and exSSI norm values of cation transition from binding sites in MD simulations of TRPV channels. Table S6 shows selectivity ratios of Ca 2+ and Na + permeation events from simulations of TRPV channels in a dicationic solution. Table S7 shows RMSF of the backbone of SF residues of TRPV channels from MD simulations. Figure S11. Pore architecture of TRPV5 channels from MD simulations with protonated (top) and deprotonated (bottom) PI(4,5)P 2 molecules. The plots show the mean radius and hydrophobicity of the channel with respect to the relative z coordinate. The shaded gray region represents the SD. The average position of binding sites A, B, and C are shown as shaded red, orange, and yellow regions, respectively. The dashed line represents the radius of a dehydrated Ca 2+ ion.