Cation Selectivity in Biological Cation Channels Using Experimental Structural Information and Statistical Mechanical Simulation.

Cation selective channels constitute the gate for ion currents through the cell membrane. Here we present an improved statistical mechanical model based on atomistic structural information, cation hydration state and without tuned parameters that reproduces the selectivity of biological Na+ and Ca2+ ion channels. The importance of the inclusion of step-wise cation hydration in these results confirms the essential role partial dehydration plays in the bacterial Na+ channels. The model, proven reliable against experimental data, could be straightforwardly used for designing Na+ and Ca2+ selective nanopores.

Introduction

Biological cation channels (Fig 1) allow the selective passage of Na+, K+ or Ca2+ through cell membranes [1-4] at passage rates near the diffusion limit. The cation discrimination occurs in the “selectivity filter” (SF, pink zone in Fig 1), a short, narrow section [5, 6] of the entire ion conduction pore (yellow volume in Fig 1).

Fig 1

Cation channels from the voltage-gated like family [7] (represented here by NavAb; PDB 3RVY) [13] consist of a four domain transmembrane protein with an ion conduction pore (yellow), a cation selective filter at the extracellular pore mouth (SF, pink) and a gate region at the intracellular pore end (green). Inward pointing carbonyl oxygens from the backbone amide groups provide a lining of the SF [5, 6]. The charged side chains present in the SFs of the bacterial Na+ channels and the Ca2+ selective CavAb [18] structure are embedded in the pore lining. Hydrated (large) and fully dehydrated (small) Na+ (blue) and K+ (green) cations are shown to scale as spheres. Black bars are 5 Å.

In the voltage-gated-like (VGL) chanome [7] theory [8, 9], simulation [10] and experiment [11, 12] have shown the critical role that hydration plays in cation channels. It is now accepted that K+ channels transport K+ ions without a hydration shell [6]. Recent X-ray structures [13-15] and corresponding MD simulations [10, 16, 17] suggest that Na+ ions are transported with a partial hydration shell through bacterial Na+ selective ion channels. Additionally, the recent atomic structures of Na+ channels and a Ca2+ channel [18] show that the charged side-chains are embedded within the SF pore wall with only one of the oxygens of the carboxylic acid group exposed.

Here we investigate the role of hydration using simplified models. Simplified models provide complementary information [19, 20] to atomistic models, giving insights on the value of specific factors for selectivity. A major advantage of such methods is that by calculating properties using rigorous statistical thermodynamics, estimates of the value of a property have well-defined convergence due to the Central Limit Theorem. The top-down approach of simplified models means they may intrinsically identify the simplest and most robust physical principles responsible for a phenomenon, thus identifying the critical features needed to intentionally engineer artificial cation-selective nanopores.

The simplified model we focus on is based on the CSC Hamiltonian (Eq 1), [21-34] which has been successfully used to study the properties of mammalian Ca2+ channels; for example explaining the anomalous mole fraction effect [35, 36]. A major feature of the original CSC Hamiltonian is that the exact formulae for the included physics are used and it therefore has no parameters which are adjusted to tune the simulation results. This contrasts with classical molecular dynamics where a force field is parameterized to reproduce a set of representative physical properties. The key output of the CSC MC simulation is a population profile of ions within the simulation cell, and in particular within the SF. The competition between the cations for occupancy of the SF, which can be deduced from these profiles, determines the ion channel selectivity (α). This is measured experimentally as the ratio of the conductance of two ionic currents. We also compare concentration profiles for different cations from a single simulation, avoiding many of the problems associated with comparing absolute values from different simulations [37].

The core of the original CSC Hamiltonian is the competition between charge and space. The competition for space is modulated by a cylinder (pink regions in Figs 2 and 3) through which the spheres representing the solute ions must pass. The original CSC model used parameters for the cylinder deduced from experiments [11] combined with solute ion spheres with sizes that matched those used to deduce the cylinder parameters. The charge competition arises through charged spheres representing the side-chain charges, located within the cylinder, that interact with the solute ions as they pass through the cylinder (Fig 2A). To get closer to an atomistic representation of the pore (Fig 2B) means the existing representation can not be used directly. We have to either embed the side-chain charge into the cylinder representing the pore or adopt an alternative representation of the pore. Our recent modification of the CSC Hamiltonian [38] allows us reproduce the position and flexibility of any atom in an X-ray structure. Using this scheme we can build up an alternative representation of the pore from the atoms in the experimental structures that actually line the pore wall (Fig 2C). This would include charged atoms from side-chains in the Na+ and Ca2+ ion channels without needing to treat them as a special case. Hence the competition for space is now modulated by the volume defined by these pore wall spheres. Implicit in using such a volume to modulate the competition for space is that the spheres representing the solute ions also have sizes that correspond to the atomistic structure (Fig 2 balls at top). These solute ions vary in hydration state, being completely dehydrated in the KcsA channel [6] and partially dehydrated in the Na+ ion channels [10, 13–17]. Thus, determining the appropriate size of solute ions involves the consideration of hydration and how full and partial hydration will change the size of these ions as they take part in the competition for space within the SF pore (Fig 2C).

Fig 2

The geometry of the SF in the original CSC model (A, pink) was empirically parameterised to reproduce the cavity of the real SF (B, yellow). In the new model we use the atoms that line the pore (B and C, red balls) to reproduce the pore cavity (C, yellow). Using a realistic SF requires considering cation hydration (Hydrated(large) and fully dehydrated (small) Na+ (blue) and K+ (green) cations as spheres).

Fig 3

Schematic view of the simulation cell used here.

The cell is a cylinder divided by a virtual membrane perpendicular to the cylinder’s axis of rotation. The model cation channel forms the only pore through the membrane and is located at the cylinder’s axis of rotation. The blue outline shows the dielectric boundary surface of the protein. The pink region shows where the spheres modelling the atoms lining the SF are located.

Here we add another important ingredient, cation hydration, to our simplified model. We describe our work to combine partial hydration and an atomistic representation of the SF pore wall in a manner conceptually consistent with the original CSC model. The good agreement between the presented model and experimental selectivity for Na+ and Ca2+ channels suggests that the included electrostatics, volume exclusion and hydration capture sufficient physics to explain selectivity for SFs containing charged residues, which are typical for Na+ and Ca2+ selective channels.

Methods

The changes to the CSC model presented here change how the model is constructed but does not change the CSC Hamiltonian used in the simulation. In this paper we present the following:

An atomistic representation of the SF pore wall. This requires modelling the oxygens of the carbonyl groups lining the SF using a partial charge on one particle as a representation of the carbonyl functional group and its dipole.

A specific model of cation hydration tailored to the CSC cation channel simulations. It introduces no change to the CSC Hamiltonian, preserving the property of a model without adjustable parameters. The hydration model consists of three parts.

The population distribution of the various hydrated cations in solution, calculated directly from the Boltzmann distribution and experimental partial dehydration energies.

A Grand Canonical Monte Carlo move that swaps a cation between hydration states.

A spherical approximation for the various hydrated cation states passing through the ion channel. The sphere diameters are Hille ‘gate-size” [11] and are based on geometric considerations and experimental data for hydrated cations in solution.

We chose the CSC Hamiltonian combined with grand canonical Monte Carlo simulations because they provide cation density profiles in the SF with well defined convergence characteristics from which selectivities can be derived. The CSC combines a Hamiltonian (Eq 1) and simulation set up that includes:

A cylindrical simulation cell divided into two equal halves by a virtual membrane (Fig 3). These halves represent the bulk solution either side of the membrane.

A toroid shape, embedded in the virtual membrane, to represent the ion channel protein (Fig 3, blue outline).

The surface of the toroid represents the dielectric boundary between the aqueous media and protein. The model computes the induced surface charges on this boundary surface (U IC) [39].

The pore of the toroid represents the ion channel pore and provides a volume connecting the bulk volumes either side of the membrane.

The radius and length of the narrowest part of the pore of the toroid are chosen to model the SF of the ion channel (more information in the supporting material). The dimensions of the SF are now fitted from the X-ray structure rather than using empirical estimates from literature.

Selected atoms from the protein are placed in the narrow pore of the toroid to represent the pore wall and charge sites in the SF of the ion channel protein.

These atoms are localized (U mob) with positions and flexibility taken directly from the coordinates and B-factors of the X-ray crystal structure [38].

Charge-charge interactions (U C) are modeled using the Coulomb equation for dielectric media.

Spatial interactions are handled using hard-body overlap (U overlap). With charged particles not being able to overlap each other or the volume contained within the protein and membrane.

The concentration of the solute ions in the simulation is maintained by a fixed chemical potential (U cp; grand canonical μVT ensemble).

Representation of the pore wall

Including charged and uncharged atoms that line the SF means we can embed charge sites into the pore wall and at the same time impose the same volume constraints that are present in the atomistic structure. These pore wall spheres are surrounded by a cylinder (shaded region in Fig 3), whose role in the competition for space is now to ensure solute ions can not go “outside” the pore wall spheres. For the ion channels studied here, the SF pore wall is composed of polar carbonyl oxygen atoms of the protein backbone [6, 13–15]. These could be represented using two spheres with opposite charges to form a dipole. However, as a first approximation we tested representing these groups using a single partial negative charge sphere. We considered that this approximation is a reasonable starting point as the influence of individual charges and asymmetry in charge reduces very rapidly with distance because our model ions exist in the strong dielectric media of water.

The embedding of the charged side groups in the SF pore wall also caused us to reconsider the earlier approach used in the CSC model of representing carboxylic acid groups as two half charged oxygen atoms. As an approaching cation would cause the charge to localize on the oxygen atom nearest to the pore wall, a single ion appeared to better capture these influential short-range electrostatic interactions better than a two ion representation. We therefore chose to represent the carboxylic group using a single charged atom in the position of the oxygen nearest to the pore wall.

Hydration model

Each cation type passing through the channel may lose differing numbers of waters of hydration in order to fit through the narrowest part of SF. Our main assumption is that we can separate the hydration and dehydration reaction from the passage through the ion channel. For Na+ and K+ cations, this assumption is reasonable since the experimentally measured residence time of waters in the hydration shell of these cations (10−10s) [40] is two orders of magnitude shorter than cation passage through an ion channel (10−8s) [11]. This is not the case for divalent Ca2+. However the much higher energetic cost of removal of waters from the Ca2+ compared to monovalent ions (Table 1) means we can exclude dehydration from consideration when examining Ca2+ selective channels. Indeed, the CSC model without hydration has already provided accurate predictions of Ca2+ channel behavior [29, 30, 41].

Table 1

Model population factor ([A⌈x⌉]bulk[A]bulk) and cation diameters (d ) for partially hydrated cations.

Partial hydration energies used are experimental gas-phase data [42]. The base diameters of the hydrated cations were calculated from the experimental cation–water oxygen (r ) and water oxygen–water oxygen (r ) radial distribution functions in ionic solutions [40, 43]; see Fig 4, .

	K+		Na+		Ca2+
	[K^+⌈x⌉]_bulk[K^+]_bulk	d K+,x	[Na^+⌈x⌉]_bulk[Na^+]_bulk	d Na+,x	[Ca^2+⌈x⌉]_bulk[Ca^2+]_bulk	d Ca2+,x
x = 6	1	8.32	1	7.44	1	7.44
5	2.2×10−2	≈ 8.32	5.8×10−3	≈ 7.44	1.8×10−16	≈ 7.44
4	5.5×10−5	6.21	1.7×10−5	5.44	4.4×10−34	5.44
3	2.5×10−8	5.54	3.7×10−10	5.10	8.5×10−53	5.10
2	7.9×10−13	2.8	6.3×10−17	2.0

Fig 4

Schematic for the circumcircles that give the partially hydrated cation diameters, d , used in the CSC simulations (not drawn to scale).

They follow the replacement of four waters of hydration (dashed circles) in a plane perpendicular to the channel axis. Removal of one water does not change the diameter (0, 1 = orange). Loss of more waters reduces the diameter (2 = blue, 3 = pink) to a minimum when four (or more) waters are lost (4 = black). The actual dimensions are calculated from the experimental cation–water oxygen distance (r(MO)) and water oxygen–water oxygen distance (r(OO)).

Populations of partially hydrated cations

The major output of the CSC MC simulation is a population profile of ions within the simulation cell, and in particular within the SF. The separation of timescales between hydration and cation passage for Na+ and K+ ions means we can divide each cation population into hydration state subpopulations, all of which are present in solution: where the total concentration of ions [A] and [B] are the sum of hydration state subpopulations. We can calculate these relative populations (see Table 1) using the Boltzmann equation and the relative dehydration free energies. This can also be applied to the SF to give: where [A⋅X] is the sum of the interactions of hydration state subpopulations of A at some point within the SF (X).

A complete model calculation would involve a simulation containing all hydration states. In practice two factors limit the practicality of such simulations. Firstly, the sub-populations B ⋅ mH2O and A ⋅ nH2O decrease by several orders of magnitude as the cation hydration number, m and n, decreases resulting in 〈[(B ⋅ mH2O) ⋅ X]〉 and 〈[(A ⋅ nH2O) ⋅ X]〉 rapidly converging to zero due to the limits of finite sized simulation. On the other hand, only values of m and n, that allow B ⋅ mH2O or A ⋅ nH2O to fit within the SF and therefore pass through the ion channel will have non-zero contributions to [B ⋅ X] and [A ⋅ X]. We therefore perform preferential sampling of the critical hydration states by ignoring values of m and n for which the values of 〈[(B ⋅ mH2O) ⋅ X]〉 and 〈[(A ⋅ nH2O) ⋅ X]〉 will be zero. This involves determining the largest hydration numbers, ⌈m⌉ and ⌈n⌉, that allows the hydrated cation to fit in the SF. These maximum numbers and consecutively lower values of m and n which have populations within three orders of magnitude of the populations of ⌈m⌉ and ⌈n⌉ are the hydration states included in the final simulations.

In these simulations the sum of the partially hydrated ions in the bulk is still that of the original [A]bulk. This enhanced sampling of the critical hydration states then requires a rescaling factor of the ratio of the sum of the sub-populations used in the simulations to the full population of A to be applied:

Grand Canonical Monte Carlo swap move

The difference in time scales between cation passage and dehydration for Na+ and K+ cations implies that not only can we model dehydration as a set of individual hydration states, but we can also include transitions between these states as a grand canonical Monte Carlo move. We therefore add a new grand canonical move into the CSC MC model that swaps cation hydration states to model variation in the hydration state profile. Such a move is a concatenation of two simple grand canonical moves; a grand canonical deletion of the particle in the original hydration state and a grand canonical addition of the particle in the new hydration state at the location of the original particle.

Model of partially hydrated cations

The CSC Hamiltonian represents a cation’s spacial component as a sphere with the cation’s charge component being a point charge in the centre of the sphere. Therefore, the only size parameter is the cation diameter (d and d ). We choose, as a first approximation, to ignore the dipolar interactions between the cation and its waters of hydration and model hydrated ions using this simple spherical model. This model also eliminates any consideration of the rotational entropy due to the potentially asymmetric shape of the partially hydrated cation cluster. We anticipate, however, that since the cluster has very little freedom of rotation while in the SF due to the size constraints of the narrow pore, a spherical model can be a reasonable first approximation. With a sphere, two choices for modeling the size seemed appropriate (i) ignore the shape of the partially hydrated cations and model the sphere based on the occupied volume, or (ii) maintain some shape information and ignore the occupied volume. The most critical spacial component will be the space the cation occludes during motion along the pore axis because the SF is such a narrow tube ions have limited lateral mobility. For example, if we consider a cation with zero, one or two waters of hydration, these water of hydration could orient along the channel axis and all three hydration states would effectively occlude the same minimal cross-sectional area along the direction of motion. This suggests that the resistance to motion along the pore axis will be in proportion to this cross-sectional area and not to the total volume of the cation complex. We thus take the critical size parameter as the area the partially hydrated cation water cluster occludes along the direction of motion. These diameters are thus Hille’s “gate size” [11] rather than, for example, an “effective” diameter determined from cation diffusion experiments in water.

K+ and Na+ cations are normally hexa-coordinated in the first hydration shell [40, 42]. We assume that two of the waters of hydration can always be oriented along the SF pore axis and that they will be the last to be removed. This sets the four remaining coordination sites from which water is removed in a plane perpendicular to the SF axis (Fig 4). Using experimental cation–water oxygen (r) and water oxygen–water oxygen (r) distances [40, 43] we calculate d and d as the diameter of the smallest circumcircle around the cation and waters of hydration that occupy from zero to four of these planar sites (see Fig 4 and Table 1).

Note that this coordination constraint only applies to the water hydration shell. The model places no limitations on the coordination of ions with other species. In particular, the cations are free to coordinate with eight atoms in the SF.

Simulation parameters

A major criterion when extending the model has been to keep the model as simple as possible and to use parameters directly derived from external sources, favoring experimental work. A major feature of the unmodified CSC Hamiltonian is that it had no parameters which are adjusted to tune the simulation results. The parameters that are in the CSC Hamiltonian are the fixed dielectric constant of water (80) and protein (10), the charge on an electron, Avogadro’s number, temperature and the Boltzmann factor. Input into the current model includes the atomic position and the B-factor of atoms from X-ray, the experimentally derived atomic radii and the geometry of the dielectric boundary surface between the aqueous and protein media. The dielectric boundary is modeled as a hard surface upon which an induced charge can develop. Previous work [30] has shown that the presence of the dielectric boundary is important but variation of the dielectric constant of the protein media (4–20) [44, 45] and the exact radius of the SF boundary surface [45] are less influential in determining selectivity.

In this work, we can not take the dimensions for the flattened toroid of the dielectric boundary representing the protein (blue line in Fig 3) directly from experiment. The most critical parameter is the diameter of the cylindrical part of the boundary surface in the SF (pink region in Fig 3). Previously, the dimensions of this boundary where critical for both control of pore volume and confinement of the electrostatic field. As we now represent the pore wall explicitly, the boundary surface no longer plays a crucial role in controlling pore volume. Additionally, previous work [45] showed that small changes in the diameter of the surface has a very small influence on the magnitude of the result. Our approach was therefore to maximise the freedom of the spheres representing the pore wall while avoiding solute cations being able to fit between these spheres and the boundary surface. Our initial estimate was a diameter slightly greater than the base atomic position plus particle radius of the spheres representing the pore wall of the SF. Test simulations were then performed to check if the diameter was too small and significantly truncated movement of the spheres of the pore wall or too large and allowed solute cations between these spheres and the boundary surface.

The parameter not taken directly from experiment or calculated through test simulations is the partial charge on the pore wall oxygen atoms. The impact of this charge will be most significant for the K+ channel model, where it is the only endogenous charge in the SF, compared to the Na+ and Ca2+ channel models that have fully charged residues in the SF. We therefore chose the value of −0.1 taken from the partial charge on the pore wall oxygen atoms seen from QM/MM MD simulations of the KcsA K+ channel [46]. This value was used universally for all pore wall oxygen atoms with no attempt to tune the parameter for specific selectivity.

Software Availability

The software used to perform these simulations is available from https://github.com/charge-space-competition/ion-channel.

Results

We simulated the SF from two bacterial Na+ channels (NavAb [13] (Fig 1) and NavMs [14]), a typical K+ channel (KcsA [6, 47]) and the Ca2+ selective mutant of the bacterial NavAb Na+ channel (CavAb) [18], the only Ca2+ selective channel in the VGL chanome [7] for which an X-ray structure has been determined. The simulated concentration profiles of solute cations within the SF have a shape characteristic for each ion channel ([A ⋅ X] and [B ⋅ X], Figs 5 and 6, 7). The location of the minima and maxima for a particular channel occur in similar ranges across cations. Interpreting selectivity from these density profiles requires application of a mechanism of cation passage. A single-barrier mechanism and the assumption the cations have similar average velocities is one of the simplest models and corresponds closely to what Hille describes as a single-ion pore [11]. The calculation of selectivity from this single-barrier mechanism is straightforward and involves comparing cation concentration ratios at the rate determining barrier in the SF (i.e. the minima in density: [A ⋅ X]min and [B ⋅ X]min) to the ratios in the bulk solution (for details see the supporting information): When hydration is included we may have a only subset of partially hydrated cations in the bulk giving: where [m] and [n] are the maximum hydration states that still allow the partially hydrated cations to fit within the SF. Correcting the bulk concentration by the population factors (Eq 2) leads to the following: Our calculated α values (Table 2) show quantitative agreement for the Ca2+ and Na+ channels and qualitative agreement for the K+ channel.

Fig 5

Simulated cation profile (as moving average over 2Å) in the Ca2+ (CavAb) ion channel SF.

[Mn2+·X] green, [Ca2+·X] blue, [Na+·X] pink and [Ba2+·X] orange. Concentration minimum region used to calculate selectivity [−5.5 : −2.5].

Fig 6

Simulated cation profile (as moving average over 2Å) in the Na+ (NavAb) ion channel SF.

[Na+·X] pink and [K+·X] black. Concentration minimum region used to calculate selectivity [−4.0 : −2.0].

Fig 7

Simulated cation profile (as moving average over 2Å) in the K+ (KcsA) ion channel SF.

[Na+·X] pink, [K+·X] black, [·X] red and [Rb+·X] or [Tl+·X] light blue. Concentration minimum region used to calculate selectivity [−6.5 : −2.5].

Table 2

Predicted equilibrium constant values from simulations based on the ion channel X-ray structures. Calculated selectivities are shown relative to the cation selected by the specific channel. (Eq 6: ).

Channel	Ions A,B	Hyd.	α hyd	α^_sim	α sim	Expt
KcsA [6]	K+,Na+	n,m ≤ 2	49,000	0.1	4,900	> 170 [47]
NavMs [14]	Na+,K+	n,m ≤ 4	0.3	30	9	11–18 [10]
NavAb [13]	Na+,K+	n,m ≤ 4	0.3	12	2.4	6–30 [16]
CavAb [18]	Ca2+,Na+	(n,m ≤ 6)	1	600	600	380 [18]

A key finding is that our model reproduces selectivity in the bacterial Na+ channels. The hydration level of cations in the Na+ channels studied is in agreement with suggestions from recent X-ray [13-15] and corresponding MD simulation [10, 16, 17] studies. Our results complement these studies by highlighting the critical role that partial dehydration plays in the selectivity. In our simulations the largest occupancies in the SF are the fourth Na+ hydration state, (Na ⋅ 4(H2O))+, and the third K+ state, (K ⋅ 3(H2O))+. Thus, even though hydration states 4 and 3 are included in our simulations for both Na+ and K+, the (K ⋅ 4(H2O))+ cation is not present in the SF. This allows the value of to overcome the population disadvantage (Na ⋅ 4(H2O))+ has over (K ⋅ 4(H2O))+ in a real aqueous solution (α hyd:Na,K = 0.3), leading to Na+ selectivity overall.

We introduced a simple model of hydration that contributes two key factors into the calculation of selectivity: (i) The population of the partially dehydrated ions as set by their relative partial dehydration energies (α sim,dehyd). (ii) The representation of the partially dehydrated cation as a sphere of a certain radius.

The first test of the hydration model and the pore wall is if it can reproduce the behavior of the original CSC Hamiltonian for Ca2+ ion channels. In our simulations, Ca2+ would not be dehydrated because it has a much larger dehydration energy than the monovalent cations (for example the loss of two water molecules would give Na+ a huge advantage α hyd,Ca,Na = 3×10−29). The Ca2+ α sim,hyd values (Table 3) for fully hydrated ions are in good agreement with experiment and earlier CSC simulations that were based on an empirical model of the SF using cation sizes equivalent to the dehydrated ions. This shows that including atomistic pore wall and hydration in the CSC model continues to reproduce the selectivity behavior of Ca2+ ion channels as seen for simple CSC models without the new features.

Table 3

Predicted selectivity from simulations based on the CavAb [18] ion channel X-ray structures using cations with radii of the fully hydrated () and fully dehydrated () species.

Ions	α sim,hyd	α sim,dehyd	Expt
Na+,Ca2+	600	100	380 [18]
Mn2+,Ca2+	0.4	0.8	Generally Mn2+ > Ca2+
Ba2+,Ca2+	200	3	Generally Ca2+ > Ba2+

To test if using hydrated ions is required we compared these selectivities with ones from simulations with cations at their fully dehydrated radii (Table 3: α sim,dehyd). All selectivities are reduced when using the smaller dehydrated radii cations, with the key selectivity between Ca2+ and Na+ falling below the experimental value. This can be explained from knowing that selectivity decreases as volume increases [45] and we assume the volume defined by the atomistic representation of the pore wall matches the size of fully hydrated cations rather than the size of dehydrated cations. Furthermore, in our model fully hydrated Ca2+ ions can not enter the Na+ selective channels whereas fully dehydrated Ca2+ ions could enter the Na+ selective channels which would lead to (results not shown) these channels being predicted as Ca2+ selective. This shows that including atomistic detail in the SF must be matched with realistically sized hydrated ions.

The second test of hydration and the representation of the pore wall is to see the results for the KcsA channel for ions with hydration energy close to that of K+, the selectivities () [48, 49] follow, in inverse order, the cation hydration energies () [50, 51] rather than the cation size () [52, 53]. From Table 4 we see that α hyd is the dominant term in determining selectivity from simulation, resulting in selectivities that are not in agreement with experiment.

Table 4

Comparison of simulated and experimental [47] selectivity values for the KcsA [6] ion channel with cations of similar full hydration energies [51].

Ions	α hyd	α^_sim	α sim	Expt (PK+PX+)
NH4+K^+	0.02	4.0	0.08	5
Rb+,K+	0.0003	1.4	0.0004	1.2
Tl+,K+	7.5	1.4	10	0.3

The KcsA channel has no charged residues in the SF and so the dominant electrostatic interaction seen in the Na+ and Ca2+ channels is not present. The main attractive interaction is the interaction of the cations with the carbonyl oxygens, the geometry allowing up to eight atoms in the SF to coordinate with a single cation. This is represented in the CSC Hamiltonian only as an electrostatic interaction between the cation and the partial charge on the carbonyl oxygens. The ability of this attractive electrostatic interaction to mimic the attractive interaction inferred from the experimental results was evaluated by comparing simulations with different charges on the oxygens atoms in pore wall of the SF. When going from -0.1 to 0.0 (Table 5, Fig 8) the major impact appears to be to modulate the concentration of ions in the pore, with smaller ions benefitting the most. However, selectivity is derived from the relative concentrations and we find no correlation between the charge and selectivity in our small sample, contrary to earlier suggestions [54]. This points to less significance of this charge for determining selectivity. Therefore the principle of parsimony suggests that these charges should be left out in future simulations. The most telling problem is Rb+ and Tl+, which are modeled identically in our simulations except for their hydration energies. These energies are in inverse order to the experimental selectivities, putting selectivities for these cations in the K+ channel beyond the scope of our Hamiltonian. All this suggests that the non-electrostatic interaction between the cation and carbonyl oxygens that line the SF are important for selectivity in the K+ channel. Handling such interactions requires a more detailed treatment [55-57].

Table 5

Comparison of change in maximum concentration and selectivity from simulations based on the KcsA X-ray structure [6] with and without “polarization” charges on the SF oxygen atoms.

Ion (A)	diameter /Å	[A ⋅ X](max,0.0)	[A ⋅ X](max,−0.1)	α_(sim,K:A,-0.1)α_(sim,K:A,0.0)
Na+	1.94	0.001	0.05	0.8
K+	2.82	0.001	0.006
Rb+	3.00	0.001	0.003	1.1
NH4+	3.20	0.001	0.001	0.8

Fig 8

Variation of [A ⋅ X]min (at maximum [A]bulk) in KcsA K+ channel simulation with “polarization” charge (red = 0:0, blue = −0.1) as a function of cation diameter.

The simulation results for the KcsA model shows the limit of our Hamiltonian where volume, hydration and electrostatic interactions are assumed to be the most dominant terms and therefore neglecting attractive non-electrostatic interactions. Our results show this assumption is reasonable for ion channels that contain charged residues in the SF, such as Na+ and Ca2+ channels, but fails for the K+ channel which is without these charged residues. The results for the Na+ and Ca2+ channels demonstrate that the dehydration model presented here, while simple, appears to capture the main features of dehydration important for selectivity within the model.

Limitations

We have presented here a model that is specific to the needs of simulating the confined geometry of ion channel SFs within the CSC method. As in any model, ours has advantages and disadvantages. Some of these have been described in previous sections. In addition, a limiting factor we have found are the limited number of published experimental values for partial hydration free energies over the full range of cations used in biological studies (for example the partial dehydration energies are not available for Cs+). Terms omitted from the current model include non-electrostatic attractive forces, a large number of details of the protein structure outside the SF, the bulk dynamics of the protein and the membrane voltage [58]. Explicit water molecules are not included, thus the complex networks of water molecules seen in atomistic simulations [59-66] are not reproduced. The model of partially hydrated cations as static spheres of discrete sizes and without dipoles can obviously be augmented with further terms, for example dynamic rearrangement of the waters of hydration affecting both the cluster shape and electrostatic field could be added.

Concluding remarks

The augmented Hamiltonian leads to satisfactory results not only for Na+ and Ca2+ channels, but also for K+ because the ion-water interactions are critical to ion channel simulations. Close inspection of α (Table 1) reveals that selectivity based purely on hydration could vary over at least eight orders of magnitude. This implies that validation of the description of cation-water interactions must be a primary requirement for all models.

The CSC Hamiltonian used in this paper complements other studies by highlighting and isolating the most important contributions to selectivity. We confirm quantitatively and through robust simulation the critical role hydration, volume exclusion and charge sites play in selectivity, as first proposed by Eisenman [8, 9] in the 1960s. We have shown how partial dehydration plays a critical role in Na+ selectivity for bacterial Na+ ion channels, as only suggested by earlier studies. However, the physics included in the CSC model is insufficient to describe selectivity in K+ channels, possibly due to the absence of attractive non-electrostatic interactions in the current Hamiltonian.

We show that for selectivity, cation hydration states must to be taken into account for any Hamiltonian when using atomistic detail in the SF. We have also shown that the CSC Hamiltonian is limited to pores where the electrostatic interactions dominate. Within this limitation, the elements of the Hamiltonian isolated here gives us a robust basis to consider designing cation selective nanopores. Thus, not only does this model predict selectivity for a subset of biological cation channels, the model can also apply to the significant simplification required in forward-engineering similar functionality.

Supporting Information.

Section A, Derivation of Eq 3 from single-barrier model and GHK equation. Section B, The mathematics of aggregating the results from sets of simulations into the single equilibrium constant. Section C, The mathematics of converting the X-ray B-factors into the harmonic potential terms used in the simulations. Section D, The pore geometry and atom positions used in the simulations. Section E, Tables of complete raw population data for Figs 5–7.

(PDF)

Click here for additional data file.