Calcium binding and permeation in TRPV channels: Insights from molecular dynamics simulations

Ion channels are important regulators of intracellular and extracellular ion concentration. The opening of ion channels generates ionic currents, which are fundamental signals in the nervous systems of animals (Hodgkin and Huxley, 1952). The selectivity of ion channels enables them to allow specific types of ions to pass through, which is a key feature of their biological function. For example, selective voltage-gated sodium channels (Na_V channels) and voltage-gated potassium channels (K_V channels) in neuronal membranes are involved in action potential generation and conduction. These channels respond to changes in membrane potential and generate currents specifically conducted by Na⁺ or K⁺. If these channels lose their selectivity, inward sodium and outward potassium currents cannot be generated and action potentials cannot spread among neurons (Abdul Kadir et al., 2018). Divalent ion channels, such as calcium channels that are weakly to highly selective for Ca²⁺, are also involved in essential physiological processes, including muscle contraction, neurotransmitter release, and hormone secretion (Ebashi and Endo, 1968; Ghosh and Greenberg, 1995; Wiederkehr et al., 2011). Some calcium channels are highly selective. For example, Ca_V channels can reach a selectivity of P_Ca:P_Na ∼ 1,000 (ion permeability ratio) despite the high concentration of sodium in the extracellular environment. How ion channels distinguish between ions with subtle differences and selectively permeate specific types of ions are fundamental scientific questions in the field of ion channel research.

For ion channels with available structures, molecular dynamics (MD) simulation is a powerful method for studying the microscopic mechanism of ion permeation and selectivity. There has been much success using this method to elucidate the permeation and selectivity mechanism of monovalent cation channels, such as potassium and sodium channels (Bernèche and Roux, 2001; Jensen et al., 2010; Corry and Thomas, 2012; Xia et al., 2013; Köpfer et al., 2014; Flood et al., 2018; Kopec et al., 2018). Dehydration of the first solvation shell of cations and interactions between multiple ions are thought to be pivotal in the selectivity mechanism of these monovalent ion channels (Mironenko et al., 2021). For calcium channels, multiple hypotheses have been proposed to explain their selectivity (Boda et al., 2000; Corry et al., 2001; Sather and McCleskey, 2003). In the early 1980s, Tsien proposed that calcium channels achieve selectivity by forming multiple binding sites that have a higher affinity for Ca²⁺ than other ions in the pore. Ca²⁺ has been proposed to permeate in a knock-off manner, which requires a repulsive force from a second ion to facilitate the detachment of Ca²⁺ from the binding site (Tsien et al., 1987). In the meanwhile, the knock-on permeation mechanism has also been widely discussed in the field of ion channel research, which describes the phenomenon that multiple ions line up to move in a concerted manner within the pore. The two terms are different, with the former focusing on the dissociation of the bound ion and the latter focusing on the collective motion, both of which resulted from ion knocking. Atomistic and dynamic details are required to discuss the difference or consistency between the two. However, MD studies on calcium channels have been hindered by the lack of accurate calcium models in classical MD.

TRPV channels (transient receptor potential cation channel subfamily V) are excellent candidates for investigating calcium channel selectivity. The TRPV family has six members, and their structures share a similar tetramer scaffold but have distinct calcium selectivity (Venkatachalam and Montell, 2007). Among these, only TRPV5 and TRPV6 are thought to be Ca²⁺-selective (P_Ca:P_Na > 100). Other members of this family, including the capsaicin receptor and thermosensor TRPV1, are considered non-selective cation channels (P_Ca:P_Na < 10; Owsianik et al., 2006). Furthermore, the open-state structures of Ca²⁺-selective TRPV6 and non-selective TRPV1 are available (Cao et al., 2013; McGoldrick et al., 2018), providing great opportunities to investigate their ion binding, permeation, and selectivity using MD simulations. In closed-state crystal rTRPV6 structures, three Ca²⁺ binding sites in the selectivity filter (SF) were identified, formed by D541, T539, and M570, respectively (Saotome et al., 2016), which was confirmed by MD simulations (Sakipov et al., 2018). In the open-state cryoEM hTRPV6 structure, only one Ca²⁺ binding site near D542 in the middle of the SF was identified (McGoldrick et al., 2018). The constriction site in the SF of TRPV1 is formed by a simple glycine and has a rather spacious vestibule above it. Compared with non-selective TRPV1, Ca²⁺-selective TRPV6 has a longer and narrower SF formed by four residues, ⁵³⁹TIID⁵⁴². How the differences in SF structures between TRPV6 and TRPV1 determine their distinct Ca²⁺ selectivity remains to be elucidated.

In this study, we conducted MD simulations for TRPV1 and TRPV6 using a recently developed multisite Ca²⁺ model (Zhang et al., 2020), which is more accurate in describing interactions between Ca²⁺ and proteins and has been successfully used to simulate multiple Ca²⁺ channels (Liu et al., 2021; Zhuang et al., 2022,Preprint; Schackert et al., 2023) to reveal the ion binding and permeation mechanisms through TRPV channels. Our results showed that despite the similar overall protein scaffolds, the differences in the SF structure resulted in distinct ion binding and permeation patterns in TRPV6 and TRPV1. The binding competition and the number of binding sites in SF could determine Ca²⁺ selectivity under physiological conditions.

The open-state structure of hTRPV6 (PDB ID: 6BO8; McGoldrick et al., 2018) used in this study was solved using cryo-EM at a resolution of 3.6 Å. The missing loop between S2 and S3 in this structure was modeled using Modeller (Šali and Blundell, 1993), utilizing the existing open structures (PDB ID: 6BO8, 7K4A; Bhardwaj et al., 2020) as the template. The whole protein structure and two truncated structures were used to build the protein–membrane simulation systems, retaining the channel pore domain (residue 475–588) and transmembrane domain (residue 317–608), respectively. For TRPV1, the DkTx (double-knot toxin)/RTX (resiniferatoxin)-bound open-state structure solved using cryo-EM at a resolution of 2.95 Å was used (Gao et al., 2016). A truncated structure containing only the transmembrane domain (residue 423–713) was used for the simulations.

A POPCbilayer in the simulation system was built using the CHARMM-GUI Membrane Builder (Wu et al., 2014). Ions were added to the system to reach different target concentrations using the CHARMM sodium model and the multisite calcium model developed by us recently (Zhang et al., 2020). The multisite calcium model was optimized for accurate calculation of calcium–protein interactions, and a recent QM/MM study further confirmed its reliability (Schackert et al., 2023). The ion concentration is the same on both sides of the membrane. The simulation systems were equilibrated using the standard CHARMM-GUI equilibration protocol before the production simulations (Jo et al., 2008), as summarized in Table S1. For the TRPV6 simulation system with the pore domain only, the system size was 9.76 × 9.76 × 9.55 nm³ and the thickness of the POPC bilayer (the distance between layers of P atoms) was about 3.7 nm. For the TRPV6 simulation system with the transmembrane domain, the system size was 10.62 × 10.62 × 11.07 nm³ and the thickness of the POPC bilayer was about 3.6 nm. For the TRPV1 simulation system with the transmembrane domain, the system size was 11.20 × 11.20 × 9.80 nm³ and the thickness of the POPC bilayer was about 3.6 nm. A transmembrane potential of 500 mV from the extracellular side to the cytosolic side was applied by adding an external constant electric field E, whose values were obtained by dividing the transmembrane voltage 500 mV with the simulation box length in the z direction. Such a high voltage is not ideal for simulating physiological conditions but is necessary for obtaining enough permeation events within accessible simulation time to get better statistics for now. All simulations were performed using the CHARMM36m force field (Huang et al., 2017) and the TIP3P water model with a time step of 2 fs. The v-rescale thermostat (Hess et al., 2008) with a time constant of 0.5 ps and the Parrinello–Rahman pressure coupling (Parrinello and Rahman, 1981) with a time constant of 5 ps were used to maintain the temperature and pressure at 310 K and 1.0 bar during the simulations. The particle-mesh Ewald method (Darden et al., 1993) was used to calculate electrostatic interactions. The van der Waals interactions were smoothly switched off from 1.0 to 1.2 nm. Gromacs 2018.6 was used to run the MD simulations and trajectory analysis (Abraham et al., 2015). Three independent 500-ns MD trajectories were collected for each simulation condition.

Due to the oversimplification of the simulation system, such as the single-component membrane model and the absence of ligand, the open-state structures of both TRPV1 and TRPV6 were found to be unstable in our MD simulations. Therefore, to maintain the open-state conformation, harmonic position restraints were applied to the majority of the α-carbon atoms of the protein with a force constant of 1,000 kJ mol⁻¹ nm⁻². While this might potentially limit the protein’s dynamics, it offers the advantage of faster convergence, which adequately serves our purpose: analyzing ion permeation through the particular open-state conformation. Still, some regions in the SF were kept free without restraints in the simulations. In simulations with the pore domain of TRPV6, the pore loop–forming SF (residue 539–542) and pore helix regions (residue 514–552) were free. In simulations with the transmembrane domain of TRPV6, the loops between the transmembrane helices and pore helix regions were free (residue 350–378, 405–424, 447–449, 470–475, 514–552). In TRPV1 simulations, the loop formed SF (residue 643–656) and the loops between transmembrane helices and pore helices (residue 456–468, 501–510, 533–537, 557–561, 600–656) were free. All the position restraints that are tested are summarized in Table S2.

We used umbrella sampling to calculate the potential of mean force (PMF) of Ca²⁺ and Na⁺ in TRPV1 and TRPV6. Systems with only one Ca²⁺/Na⁺ were prepared to calculate one-ion PMF. We applied flat-bottomed cylinder restraints with a radius of 0.5 nm and a force constant of 10,000 kJ mol⁻¹ nm⁻² to keep the selected ion around the pore axis. We first pulled the selected ion from the bottom of the box across the channel to the top of the box. Then a series of frames were extracted from the pulling trajectory with a window spacing of 0.1 nm for the ion along the z direction. For each umbrella window, a 0.1-ns equilibration was performed, followed by 1-ns production simulations. The selected ion was restrained by harmonic potential with a force constant of 1,000 kJ mol⁻¹ nm⁻² along the z direction. The PMF analysis was performed by the weighted histogram analysis method in GROMACS (Kumar et al., 1992) with the cyclic reaction coordinate, and the standard deviations were estimated from 100 bootstraps. For TRPV6, we calculated the PMF of Ca²⁺ in the presence of a Ca²⁺/Na⁺ at the upper binding site. The Ca²⁺/Na⁺ was restrained around the upper binding site with a force constant of 10,000 kJ mol⁻¹ nm⁻². The selected ion was pulled from the lower solution region to the upper binding site and the window spacing was 0.05 nm for umbrella sampling.

The coordination number of an ion is defined as the number of oxygen atoms in the first solvation shell of the ion. The radius of the first solvation shell used for both Na⁺ and Ca²⁺ is 3.0 Å in this study. In bulk solution, the coordination number is 7 for Ca²⁺ and 6 for Na⁺ in our simulations. The error of coordination number is the standard deviation of the mean values calculated from three independent trajectories, respectively.

To calculate the permeation transition probability among the binding sites in the SF, we monitored the trajectories of all the ions that came from the SF and left binding site S2 to enter the cavity. The number of ions that followed the permeation path S1→S2 was designated as N_s1, the number of ions that followed the permeation path P→S2 was N_P, the transition probability of S1→S2 was calculated as $NS1Ns1+NP$⁠, and the transition probability of P→S2 was $NPNs1+NP$⁠.

Fig. S1 shows the structural regions of transmembrane domains of TRPV6 and TRPV1. Fig. S2 shows the analysis of SF conformation at different conductance states of TRPV6. Fig. S3 shows the ion density convergence analysis. Fig. S4 shows the evolution of the z coordinates of permeating Ca²⁺ in the SF of TRPV channels. Fig. S5 shows the Na⁺ binding and permeation pattern in the pore domain of TRPV6. Fig. S6 shows the Na⁺ binding and permeation pattern in the pore domain of TRPV1. Table S1 shows the equilibration protocol and parameters generated by CHARMM-GUI. Table S2 shows the restraint conditions that were tested in our MD simulations.

In our atomistic simulations of the TRPV systems (Fig. 1 A), we applied position restraints to maintain the open conformation of the TRPV channels and found that various position restraints led to distinct permeabilities, particularly those restraints on the SF (Fig. 1 B and Fig. S1). In three 500-ns simulations of the pore domain of TRPV6, when the entire pore domain was restrained to the open-state cryo-EM structure (termed “restrained protein”), the narrow and rigid SF, with the entrance region formed by D542 having a radius of <1 Å (Fig. 1 C), was not permeable to cations and no single Na⁺ permeation was observed (Fig. 1 D). Removal of the restraints on the SF loop (termed “flexible SF”) gave limited flexibility to the SF and the radius was not significantly dilated (Fig. 1 C). No Na⁺ permeation was observed in three 500-ns simulations under such conditions for TRPV6 (Fig. 1 D). To further relax the SF region, the loops and pore helix between S5 and S6 were made fully flexible (termed “flexible pore helix”). Under such a condition, the continuous permeation of Na⁺ was observed (Fig. 1 D). The radius of the pore was dilated to ∼1.5 Å at the SF and the fluctuation increased compared with that in the more restrained situations (Fig. 1 C), allowing the passage of sodium through the SF (Fig. 1 D). The conductance calculated from the simulation with the flexible pore helix (∼62 pS) is comparable with the experimental value of 42–58 pS (Yue et al., 2001). These results highlight the importance of SF flexibility for ion permeation.

In TRPV1, the flexibility of the pore loops and helix (Fig. 1 E) also influenced the radius at SF during the simulations, thus impacting cation permeation. The open-state structure used in this study was solved using the spider toxin peptide DkTx and the vanilloid agonist RTX, whose binding to TRPV1 induces pain in humans. DkTx binds to the extracellular part of the channel and induces a conformational change in the loop between S5 and S6 (Gao et al., 2016). Contrary to TRPV6, the less restrained SF of TRPV1 showed a slightly smaller radius (Fig. 1 F), whereas the number of permeation events dropped by more than half (Fig. 1 G), resulting in a conductance of ∼66 pS, close to the experimental value (62.8 pS; Samways and Egan, 2011).

Simulations of TRPV6 and TRPV1 with different restraints applied to the pore loops and pore helix revealed that the flexibility of this region is essential for cation permeation, which is probably due to the resulting differentiable pore radii when the SF adopts permeating ions (Fig. 1, C and F). Based on the above benchmarking, position restraints were only applied to the transmembrane helices in the following production simulations, leaving the loop region and pore helix between the transmembrane helices fully flexible. It should also be noted that our simulations were conducted under a high voltage of 500 mV, so the resulting conductance may be overestimated when compared with electrophysiological measurements under physiological conditions if the I-V curves were non-linear. This was also observed in previous studies on potassium channels, and it was encouraging to see that the ion binding and permeation mechanisms were preserved (Jensen et al., 2013).

Under these conditions, we observed continuous Ca²⁺ permeation events in both channels and measured Ca²⁺ conductance as ∼18 pS for TRPV6 and 12 pS for TRPV1 (Fig. 2). These are close to experimental single-channel conductance of TRPV channels, which are summarized in Table 1. It is also noticeable that both Na⁺ and Ca²⁺ showed burst-like permeation patterns (Figs. 1 and 2), which is more obvious for Ca²⁺ than Na⁺. When breaking down the concatenated trajectory, we observed that the high-permeability burst corresponds to a more dilated SF than the low-permeability trajectory and the side-chain conformational change of D542 is directly involved in the ion permeation through TRPV6 (Fig. S2). This suggests that the ion permeation is highly sensitive to the specific SF conformation.

Our MD simulations showed that the Ca²⁺ binding within the pore achieved convergence within 200 ns (Fig. S3) and that the binding sites of Ca²⁺ in the SF regions of TRPV6 and TRPV1 were distinct. In TRPV6, two types of Ca²⁺ binding sites were identified: site S at the center of the SF and the peripheral binding site P (Fig. 3, A and B). Binding site S can be divided into two subsites: S1 and S2. S1 is located at the entrance of the SF and Ca²⁺ interacts with the D542 side chain, which was occupied 52% of the simulation time. S2 is located in the lower SF, which is a rather extended binding site that can be occupied by only one Ca²⁺ at a time. Ca²⁺ mainly interacts with the carboxyl group of T539 at S2, which was occupied 73% of the simulation time. Ca²⁺ at the binding site P mainly interacts with the negatively charged residue E535. Covering the above binding sites (both S and P), we defined the SF region of TRPV6 as the region formed by residues T539–D542. Notably, the binding sites observed in our simulations were not identical to the ion densities identified in previous structural studies. Compared with the closed-state crystal rTRPV6 structure (Saotome et al., 2016), our simulations showed two binding sites at the SF with higher electron density in the crystal structure, but the third binding site with lower electron density in the cavity was not obvious in our simulations. The open-state cryo-EM hTRPV6 structure solved the binding site located in the middle of the SF, but no Ca²⁺ density was observed at the entrance of the SF (McGoldrick et al., 2018). These differences between our simulations and existing structures can be partly explained by the gating state and flexibility of SF during ion permeation. The crystal structure was in a closed state, whereas it was difficult to capture the less stable ion binding sites in the cryoEM structure; therefore, our simulations provide more details on the ion distribution while permeating. In addition, the transmembrane potential applied in our simulations may cause a shift in the ion binding site compared with the experimental observation under transmembrane potential-free conditions.

In TRPV1, only one stable Ca²⁺ binding site at SF was observed, which was located above the narrowest site formed by G643 (Fig. 3, C and D) and occupied 86% of the simulation time. Ca²⁺ bound to this site mainly interacts with G643, M644, and G645. No peripheral Ca²⁺ binding site corresponding to site P in TRPV6 was observed in TRPV1. To be comparable with TRPV6, we defined the TRPV1 SF region in a similar way to that of TRPV6, including the wide extracellular vestibule and the narrowest site at G643.

Different Ca²⁺ binding sites in the SF region led to different modes of permeation of TRPV6 and TRPV1. In both TRPV channels, a Ca²⁺ transport process consistent with a previously proposed knock-off mechanism was observed (Tsien et al., 1987). Continuous permeation of Ca²⁺ requires the entry of Ca²⁺ into the SF region to knock off Ca²⁺ at the binding site. However, in TRPV6, since two Ca²⁺ ions are bound at the SF, permeation is characterized by a third incoming Ca²⁺, sequentially knocking the two bound Ca²⁺ in the SF. As shown in Fig. 4 A and Fig. S4, when an incoming Ca²⁺ approaches the binding site S1 of TRPV6 from the extracellular site, it knocks the Ca²⁺ bound at S1 by electrostatic repulsion. Sequentially, the Ca²⁺ bound at S1 knocks the nearby Ca²⁺ bound at S2, leading to the dissociation of Ca²⁺ to enter the channel cavity for permeation. Meanwhile, the Ca²⁺ bound at S1 moves to S2 and then S1 is occupied by a newly entering Ca²⁺ from the extracellular site. The three Ca²⁺ ions move in concerted and synergic ways to permeate, while the binding site P is not directly involved in the permeation. In contrast, for TRPV1, the calcium ion permeation process involves only one bound Ca²⁺ within the SF and one incoming Ca²⁺ from the extracellular side, without sequential knocking or concerted movement between multiple ions bound in the SF (Fig. 4 B and Fig. S4).

The collective ion permeation was also observed for Na⁺ in TRPV6 and TRPV1 in the single-cation simulations. At maximum, four and two Na⁺ can occupy the SF of TRPV6 and TRPV1, respectively (Figs. S5 and S6). A more continuous and larger ion density was observed for Na⁺ than for Ca²⁺ within both channels, partially explaining the more efficient Na⁺ permeation than Ca²⁺ in the single-cation simulations. It should be noted that here we only focus on the ions in the SF, while the ions in the cavity would probably be involved in the concerted permeation as well, as observed in the work by Ives et al. (2023). Therefore, multiple ion binding and concerted permeation appear to be a universal phenomenon in the TRPV channels.

The dehydration of permeating Ca²⁺ in the SF region has similarities but also shows discrepancies between TRPV6 and TRPV1 (Fig. 5). Ca²⁺ desolvation occurs within the SF region of both TRPV6 and TRPV1, but the desolvation in TRPV6 occurs in a longer region than that observed in TRPV1 (Fig. 5, A and B). In TRPV6, at most 1.5 water molecules in the first solvation shell were removed, and this dehydration was mainly replaced by oxygen atoms on the side chains of charged residues (D542 and E535) in the SF region (Fig. 5 A). In TRPV1, desolvation of permeating Ca²⁺ at the SF in the TRPV1 channel occurred only at the narrowest point formed by G643, where Ca²⁺ was coordinated by the oxygen atoms of the protein backbone (Fig. 5 B). It should be noted that although partial dehydration of Ca²⁺ occurs at negatively charged residues E648 and D646 above the SF of TRPV1, the pore radius at this location is larger than the radius of the first solvation shell formed by coordinated waters around Ca²⁺ (∼4 Å). Therefore, we consider the desolvation of the calcium ion here not due to spatial constraints but due to the stronger electrostatic attraction between the Ca²⁺ and the negatively charged residues. In addition, the calcium ions here do not necessarily participate in the permeation, so the ion desolvation here was not considered as permeation dehydration in this study. Therefore, the dehydration of permeant Ca²⁺ mainly occurs around the upper SF of TRPV6, whereas in TRPV1, this mainly occurs around the constriction site in the lower SF. When dehydration occurs at these constriction sites, there are always oxygen atoms of proteins that can replace water oxygens to coordinate with the permeant Ca²⁺ (Fig. 5).

In our simulations with high-concentration dicationic solutions (150 mM Ca²⁺ and 150 mM Na⁺), the blockage effect of divalent ions on monovalent currents was observed for both TRPV6 and TRPV1, although Na⁺ permeation was not fully precluded. The conductance of the channel decreased significantly when Ca²⁺ was added to the Na⁺ solution (Table 1). The stronger blockage effect observed in simulations with TRPV6 than with TRPV1 (19.1-fold decrease in TRPV6 versus 2.5-fold decrease in TRPV1) is consistent with experimental observations (Vennekens et al., 2000; Samways and Egan, 2011), which may be related to the stronger Ca²⁺ selectivity of TRPV6 over TRPV1 (Gillespie and Boda, 2008; Table 1). Despite the blockage effect of Ca²⁺, permeation events of Na⁺ were still observed in our simulations. The Na⁺ permeation patterns differed between TRPV6 and TRPV1. In TRPV6, Na⁺ is not permeable when Ca²⁺ binds to the SF. Only at the time interval between Ca²⁺ permeation when one Ca²⁺ leaves the binding site S1 can the nearby Na⁺ take the opportunity to bind to this site and cut in line to permeate (Fig. 6, A and B). Since the binding Ca²⁺ moves in concert during permeation, another Ca²⁺ will occupy the binding site soon after one Ca²⁺ leaves, so the time interval available for Na⁺ to cut in line is rather limited. Thus, in TRPV6, calcium ions show a strong blocking effect on monovalent currents. In contrast, the vestibule above the SF region in TRPV1 is rather spacious; when one Ca²⁺ is present in this region, Na⁺ can bypass it and permeate (Fig. 6, C and D). This bypass permeation pattern of Na⁺ in TRPV1 is more likely to occur than the cut-in-line pattern in TRPV6. Therefore, the blockage effect of Ca²⁺ on monovalent currents is weaker in TRPV1 than in TRPV6.

In the high-concentration dicationic systems, TRPV6 and TRPV1 showed slight valence selectivity in our MD simulations. As shown in Table 1, TRPV1 showed a slight preference for Ca²⁺ permeation over Na⁺ (∼7:1), which is consistent with the previous conclusion that TRPV1 is not Ca²⁺ selective. However, for TRPV6, we did not observe any significant Ca²⁺ selectivity either. The ratio of permeation events was ∼2:1 (N_Ca:N_Na) in our simulations with 150 mM Ca²⁺ and 150 mM Na⁺, which is inconsistent with the prediction that the permeability ratio P_Ca:P_Na is > 100:1 under physiological conditions. This is discussed further in the next section.

In the vestibule above the SF, Ca²⁺ and Na⁺ compete to occupy binding site S1. In our dicationic simulations, a significant decrease in the density of Na⁺ in the SF region of both TRPV6 and TRPV1 was observed compared to the system with only Na⁺ as a cation, indicating that Ca²⁺ repels Na⁺ from this region (Fig. 7, A and C). In TRPV1, Na⁺ density in the entire SF region decreased (Fig. 7 D). In contrast, in TRPV6, Na⁺ density showed a slight increase at the peripheral binding site P (red circle in Fig. 7 B). The differentiation of Ca²⁺ and Na⁺ binding in TRPV6 SF was observed by a detailed inspection of their 3D densities. Ca²⁺ mainly binds to the central binding site S, including S1 and S2, whereas Na⁺ is preferentially bound at the peripheral binding site P (Fig. 7 E). Ion transition probability analysis between these binding sites showed that for both Ca²⁺ and Na⁺, the ions bound at S2, which is the indispensable site during permeation, were all transferred from the S1 binding site in the single cationic systems, indicating that S1→S2 is the major permeation path. In the dicationic systems, all of the Ca²⁺ ions still followed this single permeation path, while for Na⁺, 81% of the permeation events followed the S1→S2 path and 19% followed a new P→S2 path. Most of the permeating ions once bound at site P moved to the S1 binding site for further permeation to S2, but some of them still moved directly from the binding site P to S2 (Fig. 7 F). This inequality in the transition rate to S2 indicates that the ions bound at site S1 have a significantly higher permeation probability, and those bound at site P are less involved in direct permeation. In the TRPV6 dicationic simulation system, Ca²⁺ dominantly occupied the central binding sites, having a higher probability of permeation, while Na⁺ was expelled to the peripheral binding sites P, having a lower probability of permeation. This binding site competition and differentiation between cations of various valences may also contribute to selectivity for Ca²⁺ or blockage of Na⁺ in TRPV6.

To gain quantitative insights on the ion binding competition at the SF, we conducted PMF calculations for both TRPV1 and TRPV6. As shown in the previous section, multiple cations were not only bound at the vestibule or SF but were also directly involved in the collective ion permeation. Calculating the PMF of multiple-ion permeation is very challenging as the calculation is multidimensional and hard to converge. Therefore, we adopted a simplified albeit non-physiological strategy to calculate single-ion PMFs. We removed all the surrounding ions from the protein and calculated single-ion PMFs while one ion translocated through the channels. As shown in Fig. 8, A and B, both TRPV1 and TRPV6 exhibit a strong cation binding affinity at the SF, with a strong preference toward Ca²⁺. For TRPV1, the binding strength to Na⁺ and Ca²⁺ are about −56.3 and −94.6 kJ/mol at the SF. For TRPV6, these values are about −121.9 and −213.5 kJ/mol, respectively. Obviously, the SF of both channels tends to attract cations, and the attraction is much stronger in TRPV6. Such strong binding affinities indicate that it is highly unlikely for a single cation to permeate through the channels, and multiple cations must occupy the vestibule or SF to weaken the binding strength before any permeation can occur. As expected, the binding preference of Ca²⁺ over Na⁺ for TRPV6, as represented by ΔΔE_bind = 91.6 kJ/mol, is significantly stronger than that for TRPV1 (ΔΔE_bind = 38.3 kJ/mol). It should be noted that these are non-physiological conditions without any environmental ions around the protein, so the absolute binding strength from the PMF does not mean much. However, the binding strength difference (ΔΔE_bind) is still qualitatively informative, indicating that the SF of TRPV6 should have a stronger binding preference toward Ca²⁺ than TRPV1 under very low ion concentrations.

As TRPV6 shows a two-ion binding feature within the SF, we examined how the binding of an additional cation at the S1 site affects the PMF of the translocating ion that moves from the S2 site to the cytosolic side (Fig. 8 C). When an additional Ca²⁺ is bound at S1, the depth of the PMF well at the SF changed from −213.5 to −147.9 kJ/mol. This means that the binding of an additional Ca²⁺ at S1 reduced the binding strength of the translocating Ca²⁺ within the SF for about 65.6 kJ/mol. A similar effect was observed if a Na⁺ is bound at S1, but the binding strength of the permeating Ca²⁺ is only reduced for about 40.5 kJ/mol. These results suggest that an incoming cation that binds to S1 can facilitate the permeation of the Ca²⁺ that preoccupies S2 within the SF, and the effect is more significant if the incoming ion is Ca²⁺. Again, these PMF calculations were conducted under highly simplified conditions, so these two-ion PMFs do not correspond to the observed ion permeability in previous simulations or the physiological condition in the cell. The binding strength calculated this way is much stronger than expected from the experimentally observed ion conductance under physiological conditions. However, the difference of reduced binding strength is quantitatively informative if one takes the approximations that the environmental ion distribution was not affected by the S1-bound ion type and the interactions between atoms/ions are pairwise and additive. The difference value of ∼25 kJ/mol indicates that an S1-bound Ca²⁺ indeed exerts a much stronger driving force on the S2-bound Ca²⁺ than an S1-bound Na⁺ does.

Various Ca²⁺-permeable channels showed significantly different Ca²⁺ selectivity, with reported permeability ratios over monovalent cations ranging from ∼1:1 to ∼1,000:1, which are required for these channels to execute their respective biological functions (Sather and McCleskey, 2003). How ion selectivity of various degrees is achieved has long been an intriguing question in the field of ion channel research, which requires detailed structural and dynamic insights at the atomic level. Homologous ion channels that show distinct Ca²⁺ selectivity and whose structures have been resolved, such as the TRPV channels studied in this work, provide an excellent opportunity to address this issue.

Our MD simulations explored cation permeation in both Ca²⁺-selective TRPV6 and non-selective TRPV1 channels under non-physiological (high [Ca²⁺] and high voltage) conditions. Unfortunately, despite the higher Ca²⁺/Na⁺ concentration ratio used in simulations than physiological conditions, the experimentally observed ion selectivity could not be reproduced under the simulation conditions. Nevertheless, our simulations provided insights indicating that competition between cations for multiple Ca²⁺ binding sites within the SF is likely the primary step to discriminate between Ca²⁺ and Na⁺ ions during permeation. From the permeation trajectory and PMF analysis, it was observed that Ca²⁺ showed a much stronger binding preference than Na⁺ at the SF entrance of both TRPV6 and TRPV1. In the presence of Ca²⁺, Na⁺ ions were expelled from the central binding site (S1) at the SF entrance (Fig. 7). Although Na⁺ can still occupy the lateral binding site P near the SF entrance of TRPV6, the binding site P does not lead to efficient permeation into the SF. The preferential binding of Ca²⁺ at the SF entrance can generate a blocking effect on Na⁺, which can enhance the permeation probability of Ca²⁺ over Na⁺ for both TRPV6 and TRPV1.

The number of binding sites in SF can generate distinct permeation patterns. With one Ca²⁺ binding site in the SF, Na⁺ can still frequently find their chance to enter the SF, as observed in TRPV1. A possible factor that can enhance Ca²⁺ selectivity is the presence of multiple Ca²⁺ binding sites in the SF, as observed in the SF of TRPV6. When the lower Ca²⁺ binding site (S2) is occupied by Ca²⁺, the upper binding site (S1) can be occupied by either Ca²⁺ or Na⁺, with a stronger binding preference for Ca²⁺. Meanwhile, only when S1 is occupied by Ca²⁺ does the Ca²⁺ at the S2 site in the SF receive a stronger electrostatic driving force to overcome the binding strength, as demonstrated by our PMF calculations (Fig. 8). Therefore, multiple Ca²⁺ ions must line up and march in a concerted manner to permeate more efficiently. Consequently, multiple Ca²⁺ binding sites lead to a concerted permeation manner, which is a characteristic that discriminates the Ca²⁺ permeation in TRPV6 from TRPV1, as also observed by Ives et al. in TRPV5 (Ives et al., 2023). Ives et al. conducted MD simulations on TRPV5 and TRPV2, with distance restraints to keep the gate open and otherwise similar simulation conditions to this study. The consistent results may indicate conserved ion permeation mechanisms in the TRPV ion channel family.

In fact, increasing the number of Ca²⁺ binding sites in channel pores has long been proposed to improve channel selectivity for Ca²⁺ (Tsien et al., 1987). Under physiological conditions, this concerted permeation method only allows Na⁺ permeation when both Ca²⁺ in the SF are away, which is much rarer than the case with only one Ca²⁺ binding site. Therefore, the multiple Ca²⁺ binding sites within SF can reject Na⁺ permeation more efficiently than in TRPV1, resulting in higher Ca²⁺ selectivity. However, in our simulations with high transmembrane potential, although we observed concerted Ca²⁺ permeation, Na⁺ permeation events still occurred in TRPV6. We attribute this to the fact that a non-physiologically high transmembrane voltage (500 mV) was applied in our MD simulations to obtain sufficient statistics, which would generate a much stronger driving force for cations to permeate than under physiological conditions. Thus, monovalent cations can probably still push Ca²⁺ through SF, resulting in a lower valence selectivity than expected. Indeed, Ives et al.’s work observed similar behavior for TRPV5, as well as that the selectivity increases in TRPV5 with decreasing the transmembrane potential (Ives et al., 2023). Hence, caution should be exercised when interpreting the non-equilibrium simulation results on valence selectivity presented in this study. To further evaluate this, we conducted PMF calculations and found that the binding of an incoming ion at the upper SF would facilitate the permeation of the Ca²⁺ at the lower SF of TRPV6, and the effect is more pronounced when the incoming ion is Ca²⁺ than Na⁺. Although the PMF calculation was conducted under non-physiological conditions as well, we think the difference in binding strength reduction of ∼25 kJ/mol is quantitatively informative. This supports the notion that successive occupation of Ca²⁺ in the SF would be more efficient for continuous ion permeation. Nonetheless, we should also point out the possibility that our non-polarizable multisite Ca²⁺ model, although optimized for calculating ion–protein interactions and working well for wide channels like RyRs (Zhang et al., 2020; Liu et al., 2021), may still lack the necessary accuracy for simulating accumulated Ca²⁺ ions confined in a narrow pore. In such cases, the polarizable effect becomes more significant and difficult to account for. It would be interesting to conduct further simulation studies with polarizable force fields in the future.

Ca²⁺ permeates through the SF of both TRPV6 and TRPV1 in a slightly dehydrated manner, and dehydration occurs in a more extended region of TRPV6, suggesting that the permeation pathway is more constricted for hydrated Ca²⁺ in the SF of TRPV6 than in TRPV1. Fortunately, the dehydration of Ca²⁺ in TRPV6 occurs around the charged residues in the SF and the oxygen atoms of the charged residues can take the role of coordination with the permeant Ca²⁺, thus lowering the dehydration-free energy barrier. In contrast, for TRPV1, the steric constriction site is below the electrostatic attraction site in SF. Although TRPV6 and TRPV1 utilize electrostatic and steric interactions to attract and bind Ca²⁺, they have different strategies to regulate ion permeation in the SF.

The distinct permeation mechanisms of TRPV6 and TRPV1 are based on their distinct SF structures. In the Ca²⁺-selective channel TRPV6, two adjacent binding sites line up in the narrow and long SF regions, consisting of negatively charged D542 and polar T539. This SF structure, which is ∼2 nm in length, enables Ca²⁺ to bind partially dehydrated and allows two Ca²⁺ ions to bind simultaneously in the SF region. In contrast, the SF of TRPV1 is much shorter and has a wide vestibule of negative electrostatic potential above the SF, with the narrowest region consisting of only one non-charged G643. This structure only allows for one Ca²⁺ binding site, and Ca²⁺ in the upper vestibule is dynamic and less efficient in blocking the permeation of Na⁺.

Although conventionally called knock-off (Hess and Tsien, 1984; Tang et al., 2014), the concerted Ca²⁺ permeation through multiple binding sites in TRPV6 is similar to the knock-on permeation behavior in the selective Na⁺ and K⁺ channels, which possess two Na⁺ and four K⁺ binding sites in the SF region, respectively (Corry and Thomas, 2012; Köpfer et al., 2014). However, there are evident differences in the selectivity mechanisms of these channels. In the selective Ca²⁺ channels, such as TRPV6, slightly dehydrated Ca²⁺ ions occupy two binding sites in the SF that are farther away from each other (with a separation of 11 Å), and the determinative factor for ion selectivity is likely the difference in electrostatic binding affinity and driving force generated by different valences. In this context, we can probably say that the Ca²⁺ knock-off was achieved by multiple Ca²⁺ ions’ remote knock-on. In the highly selective K⁺ channels, fully dehydrated K⁺ ions sit next to each other (with a separation of about 3.4 Å), and the determinant factor for ion selectivity is the dehydration energy difference of various ions when entering the SF. The highly selective Na⁺ channels also showed two loosely packed Na⁺ binding sites in the SF (with a separation of 10 Å), where the determinative factor was the partial dehydration energy when Na⁺ entered the SF. In fact, the difference in the (partial) dehydration energy of K⁺ and Na⁺ entering the SF can also be viewed as the difference in the binding affinities of these monovalent ions to the SF.

Therefore, by studying two representative TRPV channels with distinct Ca²⁺ binding and permeation patterns, as well as relating the results to the existing knowledge of highly selective Na⁺ and K⁺ channels, we conclude that ion binding competition at multiple binding sites in the SF of ion channels is required to generate high-blocking effect or ion selectivity. Ion competition at one binding site can generate a moderate blocking effect and ion selectivity. To achieve high blocking or selectivity while maintaining ion permeation, multiple ion binding sites with the same ion preference in SF are required. In such a scenario, the permeating ions would have to follow a concerted knock-on permeation behavior, which can simultaneously facilitate optimal ion permeability and selectivity. In addition, our results showed that not only the specific state of the channel structure but also the flexibility of the SF is essential for regulating ion permeability, highlighting the importance of protein dynamics.

It should be kept in mind, though, that our simulations were performed under non-physiological conditions (high [Ca²⁺] and high voltage), and therefore some of the results may not be quantitatively accurate. For instance, the selectivity may be transmembrane potential dependent, and a non-physiologically high voltage would abolish valence selectivity as observed in our simulations. Therefore, our conclusions on ion binding and permeation are probably more reliable than those on ion selectivity. To more quantitatively study ion valence selectivity under physiological conditions, complete ion channel structures in a realistic simulation setup, more accurate ion models that better consider polarization, physiological Ca²⁺ concentration and transmembrane voltage, and longer simulation time are desirable in future simulation studies.

The MD simulation input files and analysis scripts used for this study are deposited in a public GitHub repository, available at: https://github.com/songcgroup/Ca_in_TRPV. The other data are available from the corresponding author upon reasonable request.

Crina Nimigean served as editor.

We thank Prof. Ulrich Zachariae at the University of Dundee for helpful discussions.

This work was supported by the National Key R&D Program of China (grant No. 2021YFE0108100) and the National Natural Science Foundation of China (grant No. 21873006).

Author contributions: C. Liu and C. Song conceived the idea and designed the study. C. Liu conducted molecular dynamics simulations and analysis. L. Xue conducted additional simulations and analysis, including umbrella sampling. C. Liu and C. Song wrote the first manuscript. C. Liu, L. Xue, and C. Song revised the manuscript. C. Song acquired funding and supervised the work.