The selective exchange of ions across cellular membranes is a vital biological process. Ca2+-mediated signaling is implicated in a broad array of physiological processes in cells, while elevated intracellular concentrations of Ca2+ are cytotoxic. Due to the significance of this cation, strict Ca2+ concentration gradients are maintained across the plasma and organelle membranes. Therefore, Ca2+ signaling relies on permeation through selective ion channels that control the flux of Ca2+ ions. A key family of Ca2+-permeable membrane channels is the polymodal signal-detecting transient receptor potential (TRP) ion channels. TRP channels are activated by a wide variety of cues including temperature, small molecules, transmembrane voltage, and mechanical stimuli. While most members of this family permeate a broad range of cations non-selectively, TRPV5 and TRPV6 are unique due to their strong Ca2+ selectivity. Here, we address the question of how some members of the TRPV subfamily show a high degree of Ca2+ selectivity while others conduct a wider spectrum of cations. We present results from all-atom molecular dynamics simulations of ion permeation through two Ca2+-selective and two non-selective TRPV channels. Using a new method to quantify permeation cooperativity based on mutual information, we show that Ca2+-selective TRPV channel permeation occurs by a three-binding site knock-on mechanism, whereas a two-binding site knock-on mechanism is observed in non-selective TRPV channels. Each of the ion binding sites involved displayed greater affinity for Ca2+ over Na+. As such, our results suggest that coupling to an extra binding site in the Ca2+-selective TRPV channels underpins their increased selectivity for Ca2+ over Na+ ions. Furthermore, analysis of all available TRPV channel structures shows that the selectivity filter entrance region is wider for the non-selective TRPV channels, slightly destabilizing ion binding at this site, which is likely to underlie mechanistic decoupling.

The significance of Ca2+ 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). Ca2+ is now recognized as a versatile signaling agent, with cellular Ca2+ 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 Ca2+ 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 Ca2+ signaling cascade. The question of how Ca2+ channels selectively permeate Ca2+ 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 Ca2+ 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 Ca2+ 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 Ca2+ cations over Na+ cations (PCa/PNa ∼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)P2) in the cellular membrane and play a key role in Ca2+ homeostasis in the body (Van Goor et al., 2017). Despite their characteristic Ca2+ 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 Ca2+ and Na+ cations, even in the presence of high Ca2+ concentrations, although they are still slightly Ca2+ selective, with a permeability ratio PCa/PNa ∼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 Ca2+-permeating channels, the conventional point-charge models used to describe uncoordinated Ca2+ 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 Ca2+ and proteins, hindering an accurate study of Ca2+ permeation in channels. Previous efforts to resolve this overestimation have included polarizable force field approaches (Li et al., 2015) and rescaled Ca2+ charges (Kohagen et al., 2014a; Kohagen et al., 2014b). More recently, Zhang et al. (2020) published a new multisite Ca2+ model specifically optimized for Ca2+–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 Ca2+–protein binding energies comparable to the quantum mechanical and polarizable models. The multisite Ca2+ model has previously been used to investigate ion permeation of Ca2+ 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 Ca2+ 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 Ca2+-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) P2 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 CaCl2, or a mixture of 75 mM NaCl and 75 mM CaCl2. In the case of simulations containing Ca2+, the standard CHARMM36m Ca2+ ions were then replaced with the multisite Ca2+ 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 Ca2+ model used for simulations of CaCl2, however, is incompatible with a 4 fs time step, and therefore any simulations including Ca2+ 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 E0 of −0.03 V nm−1; this resulted in a transmembrane voltage of ∼410 mV with negative polarity in the intracellular region. The temperature was maintained at 310 K using the Nosé-Hoover thermostat (Evans and Holian, 1985) and the pressure was maintained semi-isotropically at 1 bar using the Parrinello-Rahman barostat (Parrinello and Rahman, 1981). Periodic boundary conditions were used throughout the simulations. Long-range electrostatic interactions were modeled using the particle-mesh Ewald method (Darden et al., 1993) with a cutoff of 12 Å. The linear constraint solver algorithm (Hess et al., 1997) was used to constrain bonds with hydrogen atoms. All individual simulations were 250 ns long and repeated five times for each system, as summarized in Table 2. Additional details for all simulations of Ca2+-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.

To characterize the level of cooperativity in the knock-on permeation mechanisms in TRPV channels, we used PENSA to calculate the SSI shared between discrete state transitions in the occupancy distributions of each binding site (Thomson et al., 2021,Preprint; Vögele et al., 2021; Vögele et al., 2022). A timeseries distribution with a time step of 20 ps for each binding site was obtained, whereby for each frame, if an ion occupied the binding site, then this ion’s atom ID number was recorded, whereas if the binding site was unoccupied, an ID of −1 was recorded. The ID numbers were discrete, and changes between ID numbers in each binding site, therefore, represent discrete state transitions. By quantifying the mutual information shared between changes to the ID numbers in each site, we were able to determine whether ion transitions at one site were coupled to transitions at another during a 20 ps time interval. From this, we concluded whether cations are “knocking” each other, or dissociation occurred independently from one another. The time interval was iteratively optimized to keep noise and finite sampling effects to a minimum (see below). We found that both were smallest when we used an interval of 20 ps.

Similar to McClendon et al. (2009), we observed that finite sampling resulted in independent distributions sharing mutual information (McClendon et al., 2009; Pethel and Hahs, 2014). To overcome this, we calculated a statistical threshold for each simulation via randomly permuted copies of the original data. Random permutations of the original data maintained marginal probabilities for binding site occupation in each simulation while at the same time quantifying the effect of finite sampling on the measurement of SSI. SSI was then calculated between two independently permuted versions of the occupancy distribution for the minimum entropy binding site. Since the upper bound of mutual information between two variables is equal to the lowest entropy of those variables, we used the binding site corresponding to the lowest entropy for obtaining the threshold. This ensured that the portion of SSI which could be attributed to random noise between any two binding sites was always less than or equal to the SSI. This measurement was repeated 1,000 times to resolve a Gaussian distribution from which we obtained the 99% confidence threshold. We subtracted this threshold from the measured values to resolve excess mutual information, or excess SSI (exSSI), shared in discrete state transitions. As it is not possible to transfer negative information, negative exSSI values were corrected to a value of 0.

We also derived a maximum SSI value representing a theoretical upper limit for the information that can be shared between two binding sites, where exSSImax is given by subtracting the random threshold from the minimum entropy of the two binding sites in question.

exSSI(A,B)max=min(H(A),H(B))−threshold(A,B).

(2) To quantify the interdependence of all three ion binding sites within the TRPV pores, the total correlation (TotCorr) was obtained using Eq. 3, where H(A), H(B), and H(C) represent the entropy of binding sites A, B, and C, respectively, and H(A, B, C) the joint entropy of binding sites A, B, and C.

TotCorr=H(A)+H(B)+H(C)−H(A,B,C).

(3)

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 CaCl2 or 150 mM NaCl (herein referred to as monocationic solutions) or a mixture consisting of 75 mM CaCl2 and 75 mM NaCl (herein referred to as dicationic solutions). All simulations performed with Ca2+ ions utilized the multisite Ca2+ 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 Ca2+ and 417 for Na+ in simulations of the Ca2+-selective TRPV channels.

In TRPV5 and TRPV6, Ca2+ 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 Ca2+ 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 Ca2+-selective TRPV channels conduct Na+ well in the absence of Ca2+ (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 Ca2+, but otherwise identical conditions, did not exhibit ion permeation; instead, the Ca2+ 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 Ca2+ 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 PCa/PNa values obtained from our simulations overall show lower Ca2+ 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 Ca2+ 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 Ca2+ 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 Ca2+ 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 Ca2+ channels and the development of channel-permeable models for Ca2+ ions, it had been suggested from experimental observations that Ca2+ channels may obtain their selectivity through competition, i.e., by divalent cations, such as Ca2+, 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 Ca2+ 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 Ca2+ solutions, the three binding sites showed Ca2+ 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 (tr) at the three binding sites were markedly lower than those observed for Ca2+, with ratios of tr (Ca2+):tr (Na+) varying between ∼35:1 and ∼3:1 (Fig. 3). These residence times suggest that Ca2+ 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 Ca2+ and Na+ cations are competing for the binding sites. Binding sites A, B, and C in the Ca2+-selective channels all showed high occupancies with Ca2+ in the mixed cationic solutions (Fig. 3). Across all the replicate simulations, we recorded average Ca2+ 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 Ca2+ and Na+ occupancy varies between ∼85:1 and ∼7:1 (Fig. 3). These values indicate a free energy difference of between 11.5 and 5.0 kJ mol−1 for the preferential binding of Ca2+ over Na+.

The observed increased affinity for Ca2+ cations at the pore binding sites compared to Na+ means that, in a mixed solution, Ca2+ will preferentially occupy these binding sites; however, this also implies that Ca2+ ions face a greater energy barrier when they dissociate from the binding sites. In monocationic solutions, this would result in a greatly reduced Ca2+ conductance with respect to Na+. For instance, based on our observed residence times in monocationic solutions, we would expect an ∼12-fold reduced Ca2+ unbinding rate compared to Na+ for binding site A in TRPV5. However, a much-reduced Ca2+ conductance is neither observed in our simulations nor the experimental literature. Due to the divalent charge of Ca2+, 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 Ca2+ 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 Ca2+ in TRPV5 for example (Fig. 4 left), the plot demonstrates that (i) permeating Ca2+ 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 Ca2+ 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 Ca2+-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 Ca2+ 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 Ca2+ 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 Ca2+ 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 CaCl2 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 Ca2+-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 Ca2+ 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 Ca2+-selective channels. By contrast, the occupancy of binding site A was reduced by ∼50% for both Na+ and Ca2+ 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 Ca2+ in the dicationic solutions. Note that, in the non-selective channels TRPV2 and TRPV3, this is coupled with a particularly low residence time for Ca2+ at binding site A, again suggesting higher exchange rates and weaker binding, despite its occupancy. That is, many Ca2+ 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 Ca2+-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 Ca2+-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 reduc