A cooperative knock-on mechanism underpins Ca²⁺-selective cation permeation in TRPV channels

The significance of Ca²⁺ 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²⁺ is now recognized as a versatile signaling agent, with cellular Ca²⁺ 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²⁺ 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²⁺ signaling cascade. The question of how Ca²⁺ channels selectively permeate Ca²⁺ 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²⁺ 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²⁺ 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 N- and C-terminal domains between the different subfamilies (Van Goor et al., 2020).

TRP channels are gated open by a particularly wide range of stimuli, which include temperature, small molecules, transmembrane voltage changes, and mechanical cues (Caterina et al., 1997; Tominaga et al., 1998; Caterina et al., 2000). This superfamily of genes can be divided into seven main subfamilies: TRPA (ankyrin), TRPC (canonical), TRPM (melastatin), TRPML (mucolipin), TRPN (no mechanoreceptor potential C), TRPP (polycystin), and TRPV (vanilloid). It should be noted, however, that several less well-characterized TRP subfamilies have also been reported, including the TRPY (Palmer et al., 2001; Zhou et al., 2003; Chang et al., 2010), TRPVL (Peng et al., 2015), and TRPS (Himmel et al., 2020) subfamilies. The TRPV channels are the most intensely studied channel subfamily.

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²⁺ 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₂) in the cellular membrane and play a key role in Ca²⁺ homeostasis in the body (Van Goor et al., 2017). Despite their characteristic Ca²⁺ selectivity, both channels have been shown to permeate monovalent cations such as Na⁺ when divalent cations are absent (Nilius et al., 2000; 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²⁺ and Na⁺ cations, even in the presence of high Ca²⁺ concentrations, although they are still slightly Ca²⁺ 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 (McKiernan et al., 2020; Chavan et al., 2020), and ligand-gated ion channels (Sauguet et al., 2013). However, in the case of Ca²⁺-permeating channels, the conventional point-charge models used to describe uncoordinated Ca²⁺ 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²⁺ and proteins, hindering an accurate study of Ca²⁺ permeation in channels. Previous efforts to resolve this overestimation have included polarizable force field approaches (Li et al., 2015) and rescaled Ca²⁺ charges (Kohagen et al., 2014a; Kohagen et al., 2014b). More recently, Zhang et al. (2020) published a new multisite Ca²⁺ model specifically optimized for Ca²⁺–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²⁺–protein binding energies comparable to the quantum mechanical and polarizable models. The multisite Ca²⁺ model has previously been used to investigate ion permeation of Ca²⁺ 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²⁺ 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²⁺-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.

Truncated TRPV simulation systems consisting of the membrane-domain 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 (Jo et al., 2014). In the case of the TRPV5 system, the parameters for PI(4,5) P₂ 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-oleoyl-sn-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₂, or a mixture of 75 mM NaCl and 75 mM CaCl₂. In the case of simulations containing Ca²⁺, the standard CHARMM36m Ca²⁺ ions were then replaced with the multisite Ca²⁺ 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²⁺ model used for simulations of CaCl₂, however, is incompatible with a 4 fs time step, and therefore any simulations including Ca²⁺ cations were performed with hydrogen mass repartitioning but at a time step of 2 fs. The protein was restrained in the open-state by applying harmonic restraints on the α-carbon atoms of the lower gate residues (see Table 1).

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₀ of −0.03 V nm⁻¹; 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²⁺-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.

Analysis of MD trajectory data was performed using in-house written Python scripts, utilizing GROMACS modules (Abraham et al, 2015; Lindahl et al, 2020), the SciPy library of tools (Oliphant, 2007; Pérez and Granger, 2007; Millman and Aivazis, 2011; Van der Walt et al., 2011), and MDAnalysis (Michaud-Agrawal et al., 2011; Gowers et al., 2016). Analysis of the pore architecture was performed using the Channel Annotation Package (CHAP; Klesse et al, 2019). All plots were generated in Python using Matplotlib (Hunter, 2007) and Seaborn (Waskom et al., 2018).

The selectivity (P_Ca/P_Na) was calculated as the ratio between the total sum of Ca²⁺ permeation events and the total sum of Na⁺ permeation events from fivefold replicated 250 ns simulations in dicationic solutions with 75 mM CaCl₂ and 75 mM NaCl.

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.

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.

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.

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 non-overlapping 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.

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₂ 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 root-mean-square fluctuation (RMSF) of the SF backbone. Supplemental text at the end of the PDF provides additional information.

We performed MD simulations of the pore domain of open-state TRPV5 (Hughes et al., 2018) and TRPV6 (McGoldrick et al., 2018) channels embedded in 1-palmitoyl-2-oleoyl-sn-glycerol-3-phosphocholine lipid bilayers under transmembrane voltage (approximately −410 mV). The aqueous solutions contained either 150 mM CaCl₂ or 150 mM NaCl (herein referred to as monocationic solutions) or a mixture consisting of 75 mM CaCl₂ and 75 mM NaCl (herein referred to as dicationic solutions). All simulations performed with Ca²⁺ ions utilized the multisite Ca²⁺ model developed by Zhang et al. (2020) unless otherwise stated. In both the monocationic and the dicationic solutions, the applied voltage drove a continuous flow of permeating ions through all investigated open-state TRPV channels. Overall, we recorded 433 complete inward channel crossings for Ca²⁺ and 417 for Na⁺ in simulations of the Ca²⁺-selective TRPV channels.

In TRPV5 and TRPV6, Ca²⁺ 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²⁺ 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²⁺-selective TRPV channels conduct Na⁺ well in the absence of Ca²⁺ (Nilius et al., 2000; 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²⁺, but otherwise identical conditions, did not exhibit ion permeation; instead, the Ca²⁺ 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²⁺ 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.

We note that the P_Ca/P_Na values obtained from our simulations overall show lower Ca²⁺ 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²⁺ 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²⁺ 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²⁺ 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.

Prior to the determination of the atomic structures of Ca²⁺ channels and the development of channel-permeable models for Ca²⁺ ions, it had been suggested from experimental observations that Ca²⁺ channels may obtain their selectivity through competition, i.e., by divalent cations, such as Ca²⁺, 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²⁺ 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²⁺ solutions, the three binding sites showed Ca²⁺ 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²⁺, with ratios of t_r (Ca²⁺):t_r (Na⁺) varying between ∼35:1 and ∼3:1 (Fig. 3). These residence times suggest that Ca²⁺ 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²⁺ and Na⁺ cations are competing for the binding sites. Binding sites A, B, and C in the Ca²⁺-selective channels all showed high occupancies with Ca²⁺ in the mixed cationic solutions (Fig. 3). Across all the replicate simulations, we recorded average Ca²⁺ 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²⁺ 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⁻¹ for the preferential binding of Ca²⁺ over Na⁺.

The observed increased affinity for Ca²⁺ cations at the pore binding sites compared to Na⁺ means that, in a mixed solution, Ca²⁺ will preferentially occupy these binding sites; however, this also implies that Ca²⁺ ions face a greater energy barrier when they dissociate from the binding sites. In monocationic solutions, this would result in a greatly reduced Ca²⁺ conductance with respect to Na⁺. For instance, based on our observed residence times in monocationic solutions, we would expect an ∼12-fold reduced Ca²⁺ unbinding rate compared to Na⁺ for binding site A in TRPV5. However, a much-reduced Ca²⁺ conductance is neither observed in our simulations nor the experimental literature. Due to the divalent charge of Ca²⁺, 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 flux-ratio 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²⁺ 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²⁺ in TRPV5 for example (Fig. 4 left), the plot demonstrates that (i) permeating Ca²⁺ 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²⁺ 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 simultaneous states and state transitions throughout the entire simulated time.

The SSI analysis showed that in the Ca²⁺-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²⁺ 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 knock-on 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²⁺ 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²⁺ 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₂ 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.

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²⁺-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²⁺ 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²⁺-selective channels. By contrast, the occupancy of binding site A was reduced by ∼50% for both Na⁺ and Ca²⁺ 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²⁺ 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²⁺ at binding site A, again suggesting higher exchange rates and weaker binding, despite its occupancy. That is, many Ca²⁺ 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²⁺-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 in the non-selective TRPV channels was comparable with those calculated for the Ca²⁺-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²⁺-selective TRPV channels (Fig. S8).

The ion binding sites in non-selective TRPV channels preferentially bind Ca²⁺ over Na⁺, as in Ca²⁺-selective TRPV channels, which explains why the non-selective TRPV channels in fact show slight Ca²⁺ 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²⁺ selectivity from P_Ca/P_Na ∼100:1 seen in TRPV5 and TRPV6, and in this way ultimately determines the difference between Ca²⁺-selective and non-selective 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²⁺-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²⁺-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²⁺-selective and non-selective 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²⁺-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²⁺-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²⁺-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²⁺-selective TRPV channels. This, in turn, leads to a higher occupancy of the uncharged binding site B in the Ca²⁺-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 non-selective channels. Thus, despite the constriction not being as narrow, binding site B has a higher occupancy in the Ca²⁺-selective than in the non-selective TRPV channels.

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²⁺ 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²⁺ 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²⁺ or Na⁺ cations, or indeed between the Ca²⁺-selective and non-selective TRPV channels were observed. The finding that the degree of desolvation does not differ substantially between Ca²⁺-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²⁺. Since the overall desolvation penalty for Ca²⁺ 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.

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 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) being in the closed-state, rather than the open-state structure of McGoldrick et al. (2018) used in the present study. During the transition from a closed state to the open state, TRP channels undergo a rotation about the pore-forming S6 helix, which changes the pore-facing residues (McGoldrick et al., 2018).

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²⁺ movement in the closed-state structure of R. norvegicus TRPV6 and identified cation-binding 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²⁺ 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²⁺ 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²⁺, respectively, with mostly insignificant differences between the occupancy values for Na⁺ and Ca²⁺ 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²⁺, by contrast, Ca²⁺ ions strongly outcompeted Na⁺ ions for association at all pore-binding sites, both in the Ca²⁺-selective and non-selective TRPV channels. Therefore, according to our data, all the binding sites display a much greater affinity for Ca²⁺. This gave rise to the question of how permeation efficiency for conducting Ca²⁺ 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²⁺-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²⁺ 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²⁺-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²⁺-selective channels. We suggest that this marked difference in co-operativity mechanistically explains the different levels of Ca²⁺ selectivity in TRPV channels.

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²⁺ 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²⁺ selectivity in TRPV channels is determined by the co-operativity, or coupling, between the transitions at three pore cation-binding sites. The Ca²⁺-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²⁺ permeation.

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.

Crina M. Nimigean served as editor.

We thank Professor Chen Song of Peking University for providing parameters and discussion on the multisite Ca²⁺ model. We also thank the University of Dundee I.T. services for maintenance of the School of Life Sciences high-performance computing cluster that is utilized in this research.

C.M. Ives was supported by the Medical Research Council (grant number MR/N013735/1), and N.J. Thomson was supported by the UK Research and Innovation Biotechnology and Biological Sciences Research Council (grant number BB/M010996/1). Open Access funding provided by the University of Dundee.

The authors declare no competing financial interests.

Authors’ contributions: C.M. Ives and U. Zachariae conceived the idea and designed the computational study, C.M. Ives conducted the simulations, C.M. Ives analyzed the simulation data, C.M. Ives and N.J. Thomson developed the use of SSI for knock-on permeation, N.J. Thomson analyzed the SSI data, U. Zachariae supervised the work, C.M. Ives and U. Zachariae wrote the manuscript with contributions from N.J. Thomson, and all authors edited and reviewed the manuscript.