Ion selectivity in channels and transporters

1. Historical context

The x-ray structure of the KcsA K⁺ channel offered a first chance to visualize the architecture of a highly selective biological channel (Doyle et al., 1998). The structure (Fig. 1 A) showed that K⁺ ions entering the narrowest region of the pore, the selectivity filter, are coordinated by backbone carbonyl oxygens (Fig. 1 B). The x-ray structure appeared to offer direct support for a simple and appealing structural explanation of ion selectivity in close correspondence with the snug-fit mechanism (Bezanilla and Armstrong, 1972), which thereafter became implicitly taken for granted (Hille et al., 1999). The general idea, borrowed from host–guest chemistry, is that the narrow pore is perfectly suited (at the sub-angstrom level) to provide a cavity of the appropriate size to fit K⁺, but is unable—for structural reasons—to adapt to the slightly smaller Na⁺ (Dietrich, 1985). Very early on, however, MD simulations based on all-atom models suggested that the selectivity filter of the KcsA channel might be too flexible to satisfy the requirement of the traditional snug-fit mechanism (Guidoni et al., 1999; Bernèche and Roux, 2000; Biggin et al., 2001). In MD simulations performed at room temperature, the fluctuating filter appeared to distort easily to cradle Na⁺ with little energetic cost, an observation that did not jibe with the proposed structural explanation of selectivity (Hille et al., 1999).

MD simulations and multi-ion potential of mean force (PMF) computations were used to predict the position of all the K⁺ binding sites along the pore, showing that the free energy barriers opposing ion conduction were not larger than a few k_BT and that the pore was selective for K⁺ over Na⁺ (Bernèche and Roux, 2001). All-atom alchemical free energy perturbation (FEP) MD (FEP/MD) simulations were performed to characterize ion selectivity in quantitative terms. The relative affinity of K⁺ and Na⁺, ΔΔG_Na,K, was calculated for the five cation binding sites located along the narrow pore (S0–S4; Fig. 1 B),

where $\Delta G_{\text{Na,K}}^{\text{site}}=\left[G_{\text{Na}}^{\text{site}}-G_{\text{K}}^{\text{site}}\right]\text{and}\Delta G_{\text{Na,K}}^{\text{bulk}}=\left[G_{\text{Na}}^{\text{bulk}}-G_{\text{K}}^{\text{bulk}}\right].$ (By this definition, a binding site is K⁺selective when ΔΔG_Na,K is positive and Na⁺selective when it is negative.) Such FEP/MD computations assume that ion selectivity can be understood on the basis of thermodynamic binding equilibrium, leaving out any kinetic considerations (see Section 3). The computations showed that the cation binding sites in the dynamical selectivity filter observed in the MD simulations were indeed selective for K⁺over Na⁺, particularly the site S2 located at the center of the narrow pore. At the time, the suggestion of a K⁺-selective flexible filter undergoing angstrom-scale fluctuations seemed to be conceptually incompatible with the available high-resolution x-ray structures.

To resolve this dilemma, a series of computer experiments were then performed on a wide range of atomic models of a KcsA binding site (Noskov et al., 2004a). The goal was, first, to present a convincing demonstration that ion selectivity was possible in a loosely restricted flexible binding site in which the conditions for the classical snug-fit mechanism were not realized, and second, to delineate the microscopic elements giving rise to selectivity under such conditions. The analysis led to a surprising and provocative conclusion. Robust free energy differences controlling selectivity were observed to arise principally from the number and type of ion-coordinating ligands, without the need to enforce the position of those ligands to sub-angstrom precision (Noskov et al., 2004a). Selectivity for K⁺ over Na⁺ was significant (i.e., >4–5 kcal/mol) only with eight carbonyl-like ligands. The relative free energy of K⁺ and Na⁺ in the pore resulted from a local competition between the favorable ion–ligand interactions and the unfavorable ligand–ligand interactions. This pointed to the importance of local energetics in flexible binding sites, a rather novel concept.

In this framework, the ion coordination by the carbonyl ligands was described as “dynamic” and “liquid-like,” as a contrast to the snug-fit picture with its geometry enforced to sub-angstrom positional precision. This has apparently led to some unfortunate confusion, as some have taken the semantic too literally. By liquid-like, it was meant that the carbonyls coordinating the ion are fluctuating like the ligands in the first coordination shell of an ion solvated in a liquid (i.e., the first peak in the radial distribution function of the ligands around the ion have similar widths). Furthermore, the relative free energy of K⁺ and Na⁺ in the pore arises from a balance of dynamical local interactions, which bears some similarities with ion solvation in a liquid (see Section 4). Describing the dynamics of proteins as being fluid-like (McCammon et al., 1977; Sheinerman and Brooks, 1998) or liquid-like (Caspar et al., 1988; Kneller and Smith, 1994; Lindorff-Larsen et al., 2005) is actually quite common. Nevertheless, semantic shortcuts can generate confusion because drawing a clear distinction between a material that is stiff and rigid versus one that is dynamic and flexible is partly subjective and context dependent. Notions of rigidity and flexibility in the case of ion channels have quantifiable meaning only once scales of lengths and energies have been specified (Allen et al., 2004). In physical chemistry, there is a long tradition that relies on the Lindemann melting criterion (Lindemann, 1910), according to which a material is declared to be more liquid-like than solid-like when the root-mean-square fluctuations of the atoms exceed ∼10–15% of the nearest neighbor distance. This type of analysis has often been used to contrast the folded and the so-called molten globule states of a protein (Zhou and Karplus, 1997, 1999; Jang et al., 2002). With root-mean-square fluctuations being ∼1 Å and ion-carbonyl distances ∼3 Å, the carbonyl ligands are liquid-like on atomic length scales according to the Lindemann criterion. But the analogies with ion solvation in a liquid should not be overstated. One should never lose sight that the atomic fluctuations from MD are at most on the order of ∼1 Å (Noskov et al., 2004a), and that the protein unequivocally plays an essential structural role by keeping the ion and the carbonyls restricted within some small region.

To further illustrate how size selectivity could be supported in a flexible binding site in the absence of any precise structural geometry, simple models (“toys”) were introduced with carbonyl-like groups confined to a spherical region (Noskov et al., 2004a). By construction, no restraint forces prevent the collapse of the ligands onto the ion. Yet, FEP/MD calculations showed that robust selectivity for K⁺ could occur. That such a simple model would yield selectivity for K⁺ over Na⁺ was unanticipated and quite surprising to many. After this work, several theoretical studies have examined analogous reduced toy models to explore how they could display size selectivity (Asthagiri et al., 2006; Bostick and Brooks, 2007, 2009, 2010; Thomas et al., 2007, 2011; Varma and Rempe, 2007, 2008; Dudev and Lim, 2009). Although differing in emphasis, the emerging consensus from these and our own studies has broadly confirmed that, depending on the number and type of ligands, robust ion selectivity is possible in a loosely restricted flexible binding site for which the classical notion of snug-fit does not apply.

2. The functional concept of ion selectivity

It is useful to note that the notion of ion selectivity can mean different things operationally, depending on the molecular system that is being considered (e.g., an ion channel or a transporter) and on the experimental conditions that are used to probe the system (e.g., equilibrium binding or nonequilibrium flux and ionic current measurements). In the case of a transporter, it is reasonable to expect that ion binding selectivity will be primarily governed thermodynamically by equilibrium binding, because the lifetime of the various conformational states of the protein is extremely long. This situation is unambiguously amenable to FEP/MD computations. In the case of ion channels, selectivity can mean different things, although ultimately it is understood that the wrong kinds of ions should encounter more difficulty than the right kind of ions to permeate.

Classically, ion channel selectivity has been characterized by determining the permeability ratio from the reversal potential (zero net current) under bionic conditions. Computational studies are able to simulate this situation and explain the origin of the reversal potential in an asymmetric ionic concentration with near-quantitative accuracy for OmpF porin (Im and Roux, 2002) and α-hemolysin (Noskov et al., 2004b; Egwolf et al., 2010; Luo et al., 2010). However, measurements of reversal potentials become impractical for highly selective channels such as the K⁺ channels. In this case, alternate methods such as Ba²⁺ blockade relief (Neyton and Miller, 1988a,b) or Na⁺ punchthrough (Nimigean and Miller, 2002) provide more effective strategies to characterize selectivity with quantitative accuracy. Ba²⁺ block is more sensitive to the depth of the free energy minima of the binding sites (i.e., equilibrium binding), whereas Na⁺ punchthrough is more sensitive to the height of free energy barriers (i.e., nonequilibrium rates).

Implicitly, it is assumed that the phenomenon of ion permeation and selectivity through an open conductive pore can be discussed independently from gating and inactivation. In reality, the situation is more complicated. For instance, the filter of some K⁺ channels is unable to remain in a conductive conformation when exposed to an ionic solution with little or no K⁺ ions. Under such conditions, the filter rapidly converts to a nonconductive state because K⁺ ions are needed to maintain the conductive state of the channel, giving rise to some kind of conformationally driven super-selectivity. An example of this is provided by wild-type KcsA channel (Lockless et al., 2007). Of possible relevance to these observations, the selectivity filter of the bacterial KcsA channel recapitulates most of the hallmarks of C-type inactivation displayed by eukaryotic voltage-dependent K⁺ channels (Cordero-Morales et al., 2006, 2007; Cuello et al., 2010). The complex phenomenon of super-selectivity must be distinguished from the simpler and perhaps more traditional notion of conductive-state selectivity, which is the primary subject of this Perspective. Examples of channels that unambiguously display the phenomenon of conductive-state selectivity by remaining in the conductive state even in the complete absence of K⁺ are the MthK channel (Ye et al., 2010) and the synthetic KcsA^D-Ala77 (Valiyaveetil et al., 2006). Computational studies of selectivity in K⁺ channels have mainly focused on the phenomenon of conductive-state selectivity (Bernèche and Roux, 2001; Luzhkov and Åqvist, 2001; Noskov et al., 2004a; Thompson et al., 2009; Egwolf and Roux, 2010). Less attention has been devoted to the characterization of the complex coupling between selective ion binding and the conformational stability of the filter involved in super-selectivity using computational methods. Previous studies of the relative stability of the conductive and so-called low-K⁺ structure of the wild-type KcsA channel and how it depends on ion occupancy provide examples of the challenges and difficulties that might lie ahead (Domene and Furini, 2009; Boiteux and Bernèche, 2011). Addressing these issues systematically will require advanced computational techniques able to describe complex conformational changes in macromolecules (Domene et al., 2008; Pan and Roux, 2008; Pan et al., 2008; Maragliano et al., 2009).

3. Selectivity from the perspective of the multi-ion free energy landscape

The simplest view of conductive-state selectivity starts with the assumption that there exist well-defined binding sites that can be occupied by either K⁺ or Na⁺. Within this context, the problem can be treated with FEP/MD simulations to compute the relative free energy of K⁺ and Na⁺ in the binding site of interest via Eq. 1. But conductive-state selectivity can also be a result of kinetic factors caused by differences in free energy barriers (Thompson et al., 2009; Nimigean and Allen, 2011). Furthermore, it is good to recall that ion channels are not really designed to keep ions in place, but to allow them to diffuse rapidly across the membrane. Therefore, it is likely that binding sites inside the pore, even when they exist, should be rather shallow and separated by small free energy barriers.

A more complete analysis of conductive-state selectivity requires the knowledge of the entire free energy surface governing the multi-ion permeation process, with all its free energy wells and barriers (Bernèche and Roux, 2001; Roux et al., 2004; Furini and Domene, 2009). This can be achieved by computing the PMF for the ions in the pore, W(Z₁,Z₂,Z₃), using umbrella sampling MD simulations (Bernèche and Roux, 2001, 2003). For example, the multi-ion PMF of K⁺ ions translocating through the KcsA selectivity filter displays a series of cation binding sites separated only by small free energy barriers (Bernèche and Roux, 2001), yielding a maximum channel conductance of ∼300–500 pS, in good accord with experiments (Bernèche and Roux, 2003). The predominant knock-on ion conduction mechanism extracted from the PMF calculations is supported by the results from long unbiased MD trajectories (Jensen et al., 2010).

Under physiological conditions, the functional consequence of the high selectivity of K⁺ channels is primarily to prevent the entry of Na⁺ ions from the extracellular side. To learn more about this process, all-atom MD simulations with an explicit membrane were used to determine the multi-ion PMFs corresponding to the entry of an extracellular ion, Na⁺ or K⁺, into the selectivity filter of the KcsA channel while it is in its conductive state and occupied by two K⁺ ions (Egwolf and Roux, 2010). The computations, summarized in Fig. 2, were able to capture key differences between the entry of a K⁺ or Na⁺ ion into the pore. Most importantly, the PMFs in Fig. 2 A show that the region corresponding to the binding site S2 near the center of the narrow pore emerges as the most selective for K⁺ over Na⁺. In particular, although there is a stable minimum for K⁺ in the site S2, the Na⁺ ion faces a steeply rising free energy landscape with no local well in this region. Although S2 is a stable location for K⁺, it may be a stable binding site for Na⁺ only in the context of a multi-ion configuration, with other ions bound in the sites S0 and S4. This has important implications for alchemical FEP/MD studies of ion selectivity in the KcsA channel. It is also possible to monitor the average number of coordinating oxygens from backbone carbonyl groups or water molecules around the K⁺ and Na⁺ ion during the entry process to see how it correlates with the free energy landscape. The results are shown in Fig. 2 B. As expected, the number of coordinating water molecules decreases as the ion moves from the bulk phase into the narrow selectivity filter, for Z₁ < 7 Å. Significantly, the increase in selectivity for K⁺ over Na⁺ that is observed near the site S2 appears to be predominantly correlated with changes in the type of coordinating ligands (carbonyls vs. water molecules).

The critical role of the central region near the binding site S2 of KcsA is also reflected in the pore structure of homologous channels. For example, one of the main differences between the nonselective cationic NaK channel and the K⁺-selective KcsA channel is the widening of the pore at the level of the central binding site S2 in NaK (Shi et al., 2006). This led to the suggestion that loss of selectivity at the level of the central site S2 could be the principal reason why the NaK channel is able to conduct Na⁺ unlike KcsA (Zagotta, 2006), which was correlated with MD simulations (Noskov and Roux, 2007). Results from MD about the NaK channel (Noskov and Roux, 2007) were subsequently confirmed by additional x-ray structures (Alam and Jiang, 2009, 2011).

4. Studies of ion selectivity based on reduced models

Although computations accounting for the multi-ion free energy surface are desirable (see Section 3), theoretical studies based on reduced models comprising only the ion and the nearest coordinating ligands can provide meaningful insight when selectivity is concerned with thermodynamic binding equilibrium in a well-defined site. Typically, long-range interactions from distant parts of the protein or the solvent are not expected to directly affect relative free energy differences if, for example, the charge of the ion remains unchanged, as is the case with Na⁺ and K⁺. However, long-range effects can have an indirect influence when, for example, the structural stability of a binding site is affected by a mutation in a distant residue, or it could trigger the onset of some conformation-driven super-selectivity (see Section 2). Nonetheless, describing aspects of ion selectivity from local energetics, even if there are conformational transitions, is formally possible by considering a step-by-step free energy decomposition (Deng and Roux, 2009).

5. Force field, induced polarization, and quantum mechanics (QM)

Conclusions derived from computational studies of ion selectivity depend on the accuracy of underlying computational methods. Simulations of biomolecular systems are essentially limited by the accuracy of the potential energy surface and statistical sampling (Roux, 2010b). To obtain meaningful results, one must seek a balance between more accurate but more expensive representations of the energy surface and the obligation to adequately sample the relevant configurations. Most simulations of biomolecular systems are based on potential functions with fixed effective atomic charges (Åqvist, 1990; Cornell et al., 1995; MacKerell et al., 1998), which account for many-body polarization effects in an average way. Such computationally inexpensive potential functions make it possible to extensively explore conformational space with long MD trajectories (Jensen et al., 2010). But with increasingly powerful computers, simulations based on more sophisticated representations are becoming affordable.

A promising strategy to improve the accuracy of MD studies is the development of potential functions that explicitly account for induced polarization (Grossfield et al., 2003; Patel et al., 2004; Lamoureux and Roux, 2006; Whitfield et al., 2007; Yu et al., 2010a,c). Although simulations of full proteins with polarizable force fields are still not routinely performed, there are already some interesting results concerning ion solvation in water (Lamoureux et al., 2006; Yu et al., 2010c) and in liquid amides (Grossfield et al., 2003; Yu et al., 2010a). For example, according to FEP/MD computations, there is an extremely sharp variation in the relative free energy of Na⁺ and K⁺ as a function of the number of N-methylacetamide molecules in reduced models, whereas the corresponding free energy variations with water molecules are much smaller (Yu et al., 2010a). This observation highlights the quantitative differences presented by water and carbonyls as ion-coordinating ligands.

Another area of considerable progress concerns simulations based on QM electronic structure theories (Bucher and Rothlisberger, 2010). Those can be performed in the context of Car–Parrinello MD (CPMD; http://www.cpmd.org/), which is an algorithm for generating a MD trajectory of a periodic system that approximates the Born–Oppenheimer surface of a density functional theory (DFT) method. Typically, CPMD simulations are limited to systems of only a few hundred atoms and a few tens of picoseconds and are several thousand times slower than simulations using polarizable force fields. CPMD has been used extensively to calculate the solution structure of hydrated ions (Whitfield et al., 2007) and has recently been used to model the coordination modes of ions in the KcsA selectivity filter (Bucher et al., 2007, 2010). As a word of caution, conventional DFT methods are not exact and often underestimate intermolecular attractions resulting from van der Waals dispersive forces. This can lead to limitations in representing the solvation of ions such as K⁺ in liquid water; the simulations lead to a K⁺–O radial distribution function that is broader and shifted to larger distances than the experimental and molecular mechanical force field distribution functions (Whitfield et al., 2007). As illustrated in Fig. 4, there is some improvement with CPMD simulations with an empirical correction to account for van der Waals dispersion, but the agreement with the experimental distribution functions estimated from neutron scattering and anomalous diffraction x-ray absorption-fine structure (XAFS) is poorer than that with simulations based on a polarizable force field (Whitfield et al., 2007). Notably, the coordination number of K⁺ (at r = 3.5 Å) is considerably affected by the dispersion correction. The coordination number increases from 6.15 in the uncorrected CPMD simulations to 7.45 in the dispersion-corrected CPMD simulations—a value larger than the one obtained from the non-polarizable force field (6.77) and the polarizable force field (6.80). This illustrates that the approximations currently used in CPMD simulations can substantially affect ion solvation and must be considered with caution. It cannot be assumed that a strategy based on QM methods is inherently superior to one based on a well-calibrated force field for modeling ion channel selectivity, but they remain valuable tools to parameterize and validate these simulations.

Ultimately, what matters is determining the degree of confidence regarding the conclusions that are drawn from any given computation. A good strategy is to assess the validity of a dominant trend using high-level methods. For example, the importance of the carbonyl–carbonyl repulsion in the selectivity filter of KcsA, first noted in FEP/MD simulations based on a non-polarizable force field, was subsequently assessed using high-level ab initio computations to show that the repulsion between water molecules is smaller than that between carbonyl groups (see supplemental material in Noskov et al., 2004a). More recently, the importance of the carbonyl–carbonyl repulsion was demonstrated again in the CPMD simulations of the KcsA channel performed by Bucher et al. (2010). It is also of interest to revisit and validate some previous conclusions about the high K⁺ selectivity of the binding site S2 in KcsA using FEP/MD computations with a polarizable model or a fully QM calculation based on DFT (Fig. 1 C). The results given in Table I indicate that in all cases, the relative free energy of Na⁺ and K⁺ is fairly similar, indicating that meaningful conclusions were drawn from the FEP/MD computations based on the non-polarizable force field.

6. Conclusion

Early MD simulations of the KcsA channel brought to the forefront the provocative idea that ion selectivity could robustly be achieved despite considerable structural flexibility, prompting a reexamination of the microscopic factors giving rise to ion selectivity using a wide range of computational experiments. Some of those computational experiments represented fairly artificial situations (e.g., switching off the carbonyl–carbonyl repulsion in KcsA), whereas others involved reduced models (e.g., a few confined carbonyl-like ligands). Ideas about the need for precise rigid binding sites were reassessed, and issues previously believed to be settled were reopened.

Classical notions about conductive-state selectivity such as the snug-fit mechanism, field strength, and free energy profiles are still important, but they are now sharing the stage with additional and less familiar concepts relevant to flexible and loosely restricted systems. The provocative notion that ion selectivity in a flexible binding site can emerge from local energetics has now been firmly grounded. This analysis has helped delineate the existence of two different idealized mechanisms able to support ion selectivity (Yu et al., 2010b).

Much of the progress has relied on a combination of all-atom MD simulations on known structures with explicit solvent and membrane and the analysis of reduced models. It is our contention that the reduced models, although useful for identifying and illustrating fundamental physical principles, will always be subservient to studies on real protein structures. Although features of reduced models can be a matter of debate, it is important to recall that the results from all-atom MD computations on KcsA, NaK, and LeuT have, thus far, been in good accord with experiments. This relative success of computer simulations suggests that they provide a realistic representation of these molecular systems, and that they can contribute to insight about ion selectivity.

Until now, most efforts have implicitly been focused on the factors governing conductive-state selectivity, with a fair amount of attention devoted to notions of thermodynamic equilibrium. More efforts will be required to better understand the complex phenomenon of conformation-driven super-selectivity, which remains poorly understood. Research into ion channel selectivity will continue to advance, with additional channel structures, along with improvements in computational models.

Appendix

Eq. 3, taken from Yu et al. (2010b), is closely related to Eq. 2 of Varma et al. (2011). The difference between the average ion–ligand interactions for Na⁺ and K⁺, ${\langle U_{\text{ion}-\text{ligand}}\rangle }_{\text{(Na)}}-{\langle U_{\text{ion}-\text{ligand}}\rangle }_{\text{(K)}}$, in Eq. 3 is equivalent to $\left[\Delta U_{\text{K}\to \text{Na}}^{\text{IL}}\right]$, and the difference between the average ligand–ligand interactions for Na⁺ and K⁺, ${\langle U_{\text{ligand}-\text{ligand}}\rangle }_{\text{(Na)}}-{\langle U_{\text{ligand}-\text{ligand}}\rangle }_{\text{(K)}}$, in Eq. 3 corresponds to $\left[\Delta U_{\text{K}\to \text{Na}}^{\text{LL}}+\Delta U_{\text{K}\to \text{Na}}^{\text{intra}}\right]$. Differences in entropy between ions, $-T\Delta S_{\text{K}\to \text{Na}}$, are often small in practice and can be neglected for the sake of simplicity (Roux, 2010a). The remaining term is the average internal energy contribution from the external field imposed by the rest of the protein $\Delta U_{\text{K}\to \text{Na}}^{\text{field}}$. Assuming that the external field from the protein corresponds to some harsh confining potential acting on the ligands, $U^{\text{field}}$, its average is proportional to ${\displaystyle \int d\mathbf{x}}{U}^{\text{field}}\left(\mathbf{x}\right)\exp [-U^{\text{total}}(\mathbf{x})/k_{\text{B}}T]$, which is expected to be small because the Boltzmann factor is negligible at all the points x where $U^{\text{field}}$ is nonzero (x represents all the relevant degrees of freedom).

This Perspectives series includes articles by Andersen, Alam and Jiang, Nimigean and Allen, Dixit and Asthagiri, and Varma et al. (scheduled for the June 2011 issue).