Coulomb blockade model of permeation and selectivity in biological ion channels

Figure 1(a) shows the generic, self-consistent, electrostatic model of the selectivity filter of a calcium/sodium channel whose properties we will analyse. We consider it as a negatively charged, axisymmetric, water-filled, cylindrical pore through the protein hub in the cellular membrane; and, to match the dimensions of the selectivity filters of the Na⁺/Ca²⁺ channels on which we focus, we suppose it to be of radius R = 0.3 nm and length L = 1.6 nm [14, 15, 63]. The x-axis is coincident with the channel axis, with x = 0 in the center of channel. There is a centrally placed, uniformly charged, rigid ring of negative charge $0\leqslant | {Q}_{f}/e| \leqslant 7$ embedded in the wall at ${R}_{Q}=R$. The left-hand bath, modeling the extracellular space, contains non-zero concentrations of Ca²⁺ and/or Na⁺ ions. For the Brownian dynamics simulations, the computational domain length L_d = 10 nm, its radius R_d = 10 nm, the grid size h = 0.05 nm, and a potential difference in the range $0-25$ mV (corresponding to the depolarized membrane state) was applied between the left and right domain boundaries.

Zoom In Zoom Out Reset image size

The mobile sodium and calcium ions are described as charged spheres of radius ${R}_{i}\approx 0.1$ nm for both ions (allowing use of the implicit solvent model [64, 65] with neglible ion radii), with diffusion coefficients of ${D}_{\mathrm{Na}}=1.33\times {10}^{-9}$ m² s⁻¹ [15, 66] and ${D}_{\mathrm{Ca}}=0.79\times {10}^{-9}$ m² s⁻¹ [15], respectively.

We take both the water and the protein to be homogeneous continua describable by relative permittivities ${\varepsilon }_{w}=80$ and ${\varepsilon }_{p}=2$, respectively, together with an implicit model of ion hydration whose validity is discussed elsewhere [58]. We approximate _w, D_Na, and D_Ca as being equal to their bulk values throughout the whole computational domain (see below, section 2.4).

The electrostatic potential U for an ion, and the potential gradients, were derived by numerical solution of Poisson's electrostatic equation within the computational domain shown in figure 1:

$$\begin{eqnarray}&&-\nabla ({\varepsilon }_{0}\varepsilon \nabla \phi )={\rho }_{0}+\sum _{i}{{ez}}_{i}{n}_{i}\quad \quad ({\rm{Poisson's}}\;{\rm{equation}}),\end{eqnarray} \tag{ 1 }$$

where is the relative permittivity of the medium (water or protein), ρ₀ is the density of fixed charge, z_i is the charge number (valence), and n_i is the number density of moving ions.

Self-consistent electrostatics within the narrow, water-filled, channel in the protein differs significantly from bulk electrostatics [49, 67, 68]. The huge gradient between ${\varepsilon }_{w}=80$ and ${\varepsilon }_{p}=2$, the discreteness of ionic charge and the specific channel geometry, lead to permeation, quasi-1D axial behaviour of ions inside the channel [49, 50, 69], and hence to single-filing. Consequently, we use a 1D dynamical model to simulate the axial single-file movement of cations (only) inside the selectivity filter and in its close vicinity.

Figure 1(b) shows axial single-ion potential energy profiles for ${Q}_{f}=-1e$, including the repulsive self-energy barrier U_s and the attraction energy U_a attributable to the fixed negative charge. The dielectric self-energy U_s plays a crucial role in controlling ion permeation through the narrow channel. It will be shown below that U_s defines the strength of Coulomb blockade for an ion of particular valence z, and that it determines the condition for strong valence selectivity [30, 58]. The site attraction is proportional to $z\times {Q}_{f}$, so that a variation of Q_f can significantly change the resultant profile. In particular for ${Q}_{f}\simeq -1e$, the self-potential barrier of the dielectric boundary force can be balanced by electrostatic attraction to the fixed charge Q_f, resulting in almost barrier-less conduction [58].

The ions' motion through a channel can be described as an electro-diffusion process and investigated through Brownian dynamics simulations [15, 30, 44, 45, 57, 70–73] of the ionic trajectories. Taking account of the single-filing required by electrostatics, we can solve the over-damped, time-discretized, Langevin equation numerically. An axial step $\Delta {x}_{i}$ of the ith ion is defined as [74],

$$\begin{eqnarray}&&\Delta {x}_{i}=-{{zeD}}_{i}{\nabla }_{x}\phi ({x}_{i})\Delta t+\sqrt{2{D}_{i}\Delta t}\ {\xi }_{i}(t)\quad \quad ({\rm{Langevin}}\;{\rm{equation}}),\end{eqnarray} \tag{ 2 }$$

where D_i is the ionic diffusion coefficient, $\Delta t$ is the time step, ${\xi }_{i}(t)$ is normalized white noise, z_i is the valence of the ith ion, and the potential $\phi ({x}_{i})$ is calculated self-consistently from (1) at each simulation step.

Brownian dynamics simulations of the ion current J and occupancy P were performed separately for CaCl₂ and NaCl solutions, and also for a mixed-salt configuration, with concentrations $[\mathrm{Na}]=30$ mM and $20\;\mu {\rm{M}}$ $\leqslant \;[\mathrm{Ca}]\leqslant 80$ mM.

Our reduced model obviously represents a considerable simplification of the actual electrostatics and dynamics of moving ions and water molecules within the narrow selectivity filter. The validity and range of applicability of this kind of model have been discussed in detail elsewhere [35, 58, 75, 76]. The most significant simplifications are probably: the use of continuum electrostatics; the use of the implicit solvent model; the use of Brownian dynamics with uncorrelated noise sources for the charged particles; and the assumption of 1D (i.e. single-file) movement of ions inside the selectivity filter. We can partially accommodate these simplifications by use of effective values [58].

Figures 2(a), (b) and 3(a), (b) present the Brownian dynamics results [57–59] for permeation of the model channel figure 1(a) by calcium and sodium ions, respectively, in pure baths of different concentration, plotted in each case as a function of Q_f.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 2(a) exhibits a sequence of narrow conduction bands M₀, M₁, M₂, separated by stop-bands of almost zero-conductance centred on the blockade points Z₁, Z₂, Z₃. Figure 2(b) shows that the M_n peaks in J_Ca correspond to transition regimes where the channel occupancy P_Ca jumps from one integer value to the next, and that the stop-bands correspond to saturated regions with integer ${P}_{\mathrm{Ca}}=1,2,3...$.

The band M₀ corresponds to single-ion low-barrier conduction (see red curve in figure 1(b)). Band M₁ corresponds to double-ion knock-on conduction, which is well-established for L-type Ca²⁺ channels [14, 77, 78]; a similar peak was obtained by Corry et al [15] in Brownian dynamics simulations of the Ca²⁺ channel. Band M₂ corresponds to triple-ion conduction, a process that can be identified with the permeation of ryanodine receptor calcium channels [7] (see table 1 below for further detail). These bands can be considered as examples of self-organization in ion channels [16, 73].

Table 1.  Putative identification of ionic Coulomb blockade model bands with some known wild type (WT) channels and mutants sequences (extended from [58]).

N	Channel/mutant	Site locus	Nom. ${Q}_{f}/e$	Band	Band's ${Q}_{f}/e$
(a)	WT Na+-selective Nav channel [22]	DEKA	$-1$	L${}_{0}$	$-0.5$
	WT Non-selective OmpF [27]	KRRRDE	$-1$	M0	$-1$
	WT Ca2+-selective L-type channel [23]	EEEE	$-4$	M1	$-3$
	WT NaChBac[18], NavAB [28]	EEEE	$-4$	Z2	$-4$
	WT Ca2+-selective RyR channel [7]	DDDD(ED)	$-6$	M2	$-5$
(b)	WT Na+-selective Nav channel [22]	DEKA	$-1$	L${}_{0}$	$-0.5$
	Ca2+-blocked Nav mutant	$\Rightarrow$ DEKE	$-2$	M0	$-1$
	Ca2+-permeable Nav mutant	$\Rightarrow$ DEEA	$-3$	Z1	$-2$
	Ca2+-selective Nav mutant	$\Rightarrow$ EEEE	$-4$	M1	$-3$
(c)	WT Ca2+-selective L-type channel [23]	EEEE	$-4$	M1	$-3$
	Ca2+-blocked Cav mutant	$\Rightarrow$ EEQE	$-3$	Z1	$-2$
	Na+-conductive Cav mutant	$\Rightarrow$ EEKE	$-2$	M0	$-1$
	Na+-selective Cav mutant	$\Rightarrow$ EEKA	$-1$	L0	$-0.5$
(d)	WT Na+-selective NaChBac [18]	EEEE	$-4$	Z2	$-4$
	Ca2+-selective CaChBac mutant	$\Rightarrow$ EEEE+DDDD	$-7$	M3	$-8$
(e)	WT Na+-selective NavAB [28]	EEEE	$-4$	Z2	$-4$
	Ca2+-selective CavAB mutant	$\Rightarrow$ EEEE+DDDD	$-7$	M3	$-8$

Comparison of the J_Ca and P_Ca plots in figures 2(a) and (b) shows that for the M_n points near which conduction occurs, ${Q}_{f}=-{ze}(n+1/2)$; whereas the non-conducting regions of constant P correspond to Z_n points with ${Q}_{f}=-{zen}$, i.e. to the neutralized state [59]. These bands will be interpreted below as strong Coulomb blockade.

Figures 3 (a) and (b) plot the equivalent results for sodium ions in pure NaCl baths of different concentration showing (a) the sodium current and (b) the occupancy as functions of Q_f. The current J_Na exhibits a single-ion peak L₀ that would appear to be analogous to the calcium conduction band M₀ of figure 2(a). For larger Q_f there are weak, strongly overlapping, conduction bands between which the current does not fall to zero, making the sodium conductance relatively independent of Q_f. The separations of the J_Na band maxima are approximately half the size of those for the calcium bands, reflecting the charge difference between Na⁺ and Ca²⁺ ions. We will refer to this scenario as an example of weak Coulomb blockade.

Figures 2(c) and 3(c) plot the ground state energies characteristic of Coulomb blockade graphs [56]. They will be discussed in the next section. Note that we use ground state in the physics sense, implying the state of minimum energy.

We are now in position to introduce the Coulomb blockade model of conduction and selectivity in Ca²⁺/Na⁺ ion channels. We will find that it is able to account for the pattern of calcium and sodium bands seen in figures 2 and 3 in terms of strong and weak Coulomb blockade oscillations [55], respectively .

The discreteness and entity of the ionic charge allow to us to introduce exclusive 'eigenstates' $\{n\}$ of the channel for fixed integer numbers of ions inside its selectivity filter, with total electrostatic energy U_n. The transition $\{n\}\to \{n+1\}$ corresponds to the entry of a new ion, whereas $\{n\}\to \{n-1\}$ corresponds to the escape of a trapped ion. The n ions' eigenstates form a discrete exclusive set of $\{n\}$-states [79] :

$$\begin{eqnarray}&&n=\{0,1,2,...\}\quad \sum _{n}{\theta }_{n}=1;\quad {P}_{c}=\sum _{n}n{\theta }_{n},\end{eqnarray} \tag{ 3 }$$

The electrostatic exclusion principle (3) leads to Fermi–Dirac statistical distributions [80] for θ_n and P_c, as will be derived below.

The total energy U_n for a channel in state $\{n\}$ can be expressed as:

$$\begin{eqnarray}&&{U}_{n}={U}_{n,s}+{U}_{n,a}+{U}_{n,\mathrm{int}}\end{eqnarray} \tag{ 4 }$$

where ${U}_{n,s}$ is the self-energy, ${U}_{n,a}$ is the energy of attraction, and ${U}_{n,\mathrm{int}}$ is the ions' mutual interaction energy.

We approximate U_n as the dielectric self-energy ${U}_{n,s}$ of the excess charge Q_n, based on the assumption that both the ions and Q_f are located within the central part of the selectivity filter, so that application of Gauss's theorem to the n similar ions captured within its volume gives a Coulomb blockade-like quadratic dependence of U_n on Q_f [51]:

$$\begin{eqnarray}{U}_{n} & = & \displaystyle \frac{{Q}_{n}^{2}L}{8\pi {\varepsilon }_{0}{\varepsilon }_{w}{R}^{2}}=\displaystyle \frac{{Q}_{n}^{2}}{2{C}_{s}}\quad ({\rm{}}\;{\rm{Electrostatic}}\;{\rm{energy}}\;{\rm{}})\\ {Q}_{n} & = & {zen}+{Q}_{f}\quad \quad \quad \quad ({\rm{}}\;{\rm{Excess}}\;{\rm{charge}}).\end{eqnarray} \tag{ 5 }$$

Here, C_s stands for the geometry-dependent self-capacitance of the channel, and Q_n represents the excess charge at the selectivity filter for the n ions as a function of Q_f. A binomial expansion Q_n² in (5) gives first approximations for ${U}_{n,s}$, ${U}_{n,a}$ and ${U}_{n,\mathrm{int}}$ that are consistent with the energetics analysis in [58] and with the 1D Coulomb gas model of ion–ion and ion-fixed charge interactions [49, 50]. A more realistic account of the interactions would result in corrections to (5) and to the formulæ derived from it.

With (5) we arrive at an equation identical to that for electronic Coulomb blockade (except for the presence of z), and our further consideration follows standard Coulomb blockade theory [53, 55]. Remarkably, however, the ionic version of the phenomenon exhibits valence selectivity precisely because it contains the valence z.

To interpret the conduction bands in terms of Coulomb blockade, we calculate U_n as a function of Q_f for $n=0,1,2,3$ and examine the minimum (i.e. ground state) energy

$$\begin{eqnarray}&&{U}_{G}({Q}_{f})=\underset{n}{\mathrm{min}}({U}_{n}({Q}_{f},n))\quad ({\rm{Ground}}\;{\rm{state}}\;{\rm{energy}})\end{eqnarray} \tag{ 6 }$$

for the ground state occupancy n_G, and the excess charge Q_G, all as functions of Q_f (see figures 2(c) and 3(c)). The ground state (6) is by definition a stable state in thermal equilibrium.

Coulomb blockade appears in low-capacitance systems on account of quantization of the quadratic energy in (5) on a grid of discrete states (3), provided that the ground state {n_G} is separated from neighbouring $\{{n}_{G}\pm 1\}$ states by large Coulomb energy gaps as function of geometry ($R,L$) and ion valence z:

$$\begin{eqnarray}&&{U}_{s}=\displaystyle \frac{{e}^{2}}{2{C}_{s}}=\displaystyle \frac{{\lambda }_{{\rm{B}}}L}{2{R}^{2}},\quad {U}_{s,z}={z}^{2}{U}_{s}\gg 1\quad ({\rm{Coulomb}}\;{\rm{gap}}),\end{eqnarray} \tag{ 7 }$$

where ${\lambda }_{{\rm{B}}}$ stands for Bjerrum length [49] (0.7 nm for water at T = 298K). This is the applicability condition for the strong electrostatic exclusion principle. Ion channels are extremely small and have tiny capacitance: the dimensionless self-energies of monovalent and divalent ions are ${U}_{s,1}\approx 5$ and ${U}_{s,2}\approx 20$, respectively, providing strong Coulomb blockade effects for divalent ions.

It follows from (5) that U_n versus Q_f for given z is described by an equidistant set of identical parabolæ of period equal to the ionic charge ${ze}$. These patterns are plotted for Ca²⁺ in figure 2(c) and for Na⁺ in figure 3(c). We note that ${U}_{G}({Q}_{f})$ exhibits two different kinds of ground state singular points, marked as M_n and Z_n. The minima of U_n (and correspondingly the blockade regions) appear around the neutralization points Z${}_{n}=-{zen}$ where the net charge at the selectivity filter Q_n = 0 and the occupancy P_c is saturated at an integer value [51, 58].

Figure 2(c) illustrates the fact that for Ca²⁺ the ground state $\{{n}_{G}\}$ Z_n points are separated from neighbouring $\{{n}_{G}\pm 1\}$ states by an impermeable Coulomb gap of $\sim 20$ (see (7)) hence providing strong Coulomb blockade. At the neutralized (Q_n = 0) blockade points Z_n, Ca²⁺ ions are prohibited by the self-energy barrier from entering the uncharged channel. The crossover points M_n (${U}_{n}={U}_{n+1}$) allow low-barrier (almost barrier-less) $\{n\}\leftrightarrows \{n+1\}$ transitions (see figure 1 (b)); they correspond to the P_c transition regions and to the conduction peaks in J [58]. The M_n points are separated from higher energy states by impermeable barriers of $\sim 40$; these points are equivalent to the ion-exchange transitions reported by Zhang et al [51].

The pattern of sodium ground states in figure 3(c) arises from energy gaps U_c that are 4× smaller than for calcium, They are too small to prevent thermally activated transitions between neighbouring states and they allow the coexistence of more then two $\{n\}$-states; correspondingly there is only a weak exclusion principle, in turn giving rise to weak Coulomb blockade.

The positions of the singular Q_f points in figure 2(a) can be written as:

$$\begin{eqnarray}{{\rm{Z}}}_{n} & = & -{zen}\pm \delta {Z}_{n},\quad \quad \quad \quad \;(\mathrm{Coulomb}\;\mathrm{blockade}),\\ {{\rm{M}}}_{n} & = & -{ze}(n+1/2)\pm \delta {M}_{n}\quad (\mathrm{Resonant}\;\mathrm{conduction}),\end{eqnarray} \tag{ 8 }$$

where $\delta {Z}_{n}$, $\delta {M}_{n}$ are possible corrections for the singular parts of the affinity and ion–ion interaction (see above) and possible electrostatic field leakage. Equations (8) describe two interleaved periodic sets of points having periods equal to the ionic charge ${ze}$, very similar to their counterparts in electronic Coulomb blockade [54].

The points Z${}_{n}=-{zen}$ are neutralization points, where the excess charge Q_n = 0 and U_G takes minimal values. Such points are stable and current is prohibited. Charge neutrality is important, but it is not the only factor that influences channel conductance: the term $\pm \delta {Z}_{n}$ accounts formally for the other factors (unaccounted parts of ion–ion/ion-ligand interactions and hydration energy, fields leaks etc).

The points ${M}_{n}=-{ze}(n+1/2)$ are where barrier-less conduction occurs, for which ${U}_{n}={U}_{n+1}$. Similarly, the terms $\pm \delta {M}_{n}$ account formally for any unconsidered perturbations.

We may therefore interpret the Brownian dynamics-simulated calcium conduction bands of figure 2(a) as Coulomb blockade conductance oscillations [55] which appear in this case as $| {Q}_{f}|$ is being increased, and the corresponding occupancy steps in figure 2(b) as a Coulomb staircase [56]. The deviations in the precise positions of M_n and Z_n can be attributed to field leaks and the model simplifications.

Table 1, showing the putative identifications of the bands/singularities in the Coulomb blockade model with real channels, wild type and mutants in the Ca²⁺/Na⁺ family, has been extended from [57, 58].

Table 1(a) describes shows identifications of conduction bands with wild type channels.

The single-ion Na⁺ barrier-less point L₀ can be speculatively identified with the non-selective DEKA sodium channel (${Q}_{f}=-1e$) [1, 17]. The single-ion Ca²⁺ barrier-less point M₀ can be identified with the non-selective OmpF porin (${Q}_{f}=-1e$) [27], or with nonselective Na⁺–K channel [81].

Mammalian calcium channels exist in several modifications. Some of them (L-type, T-type) share the same highly conserved four-glutamate (EEEE) locus at the selectivity filter with nominal ${Q}_{f}=-4e$ [14]. The permeation properties (sharp selectivity, AMFE, double-ion nature of calcium permeation) of these channels are consistent with the double-ion M₁ Ca²⁺ conduction band [57, 58]. The double-ion knock-on mechanism of conduction and selectivity in the L-type calcium channel has been derived from the experimentally observed double-affinity of AMFE [14, 78]. The same M₁ peak can be also identified with the Ca²⁺-selective mutant of the Nav sodium channel (${Q}_{f}=-3e$) [22] and a calcium-selective mutant of OmpF porin (${Q}_{f}=-4e$) [27].

The Ryanodine receptor RyR calcium channel has ${Q}_{f}\approx -6e$ and relatively weak selectivity [7]. We connect it with M${}_{2}\approx -5e$ three-ion conduction point [57, 58].

Bacterial sodium channels, NaChBac [18] and NavAB [28] possess EEEE loci (with nominal ${Q}_{f}=-4e$), typical of mammalian calcium channels but nonetheless exhibit Na⁺-selective features and divalent blockade. We connect them (speculatively) with the Z₂ Ca²⁺ blockade point for ${Q}_{f}=-4e$.

It is evident that the nominal charges of the channel loci exceed the charge of the model bands by nearly $\Delta {Q}_{f}=1e$. For example, at M₁ we have ${Q}_{f}=-3e$ from the Coulomb blockade model, whereas the nominal ${Q}_{f}=-4e$ for the EEEE loci of L-type calcium channels (and there are corresponding discrepancies at the M₂ and M₃ points). We speculate that these systematic discrepancies may be connected to field leaks, to the distant influence of the positive charges of the other ends of the amino acids buried in the rest of the protein, or to possible protonation of the side chains [82]. Molecular dynamics simulations [28, 35, 38] could be particularly helpful in determining the effective values of Q_f for wild type and mutant channels.

The above identification scheme can thus account for many mutation transformations observed in the Ca²⁺/Na⁺ channels family. The main test and validation of the Coulomb blockade model is the correct prediction that single $\pm 1e$ mutations should lead to sharp changes of calcium conductance from resonance to blockade and vice versa.

Similarly, table 1(b) describes the radical mutation-induced transformation of the Nav sodium channel to a calcium channel, when the DEKA locus has been sequentially changed by point mutations [22].

Table 1(c) shows that the downward mutations of EEEE calcium channels described in [23] lead to a moderate decrease of divalent/monovalent (Ba²⁺/Na⁺) selectivity for EEEE $(-4e)\to \mathrm{EQEE}(-3e)$ mutation and a sharp, two-order decrease for the EEEE $(-4e)\;\to$ EKEE $(-2e)$ mutation (though see discussion below).

Tables 1(d), (e) show that the upward mutations of the EEEE loci of bacterial sodium channels (EEEE $(-4e)\;\Rightarrow$ EEEE+DDDD $(-8e)$ calcium channels described in [18, 28] lead to channel transformation to a calcium-selective mutant CaChBac and CavAB, respectively. Ca²⁺ selective mutants demonstrate selectivities of up to $600\times$ in favour of Ca²⁺ ions over Na⁺, and AMFE. We putatively connect these transformations with the jump ${Z}_{2}\Rightarrow {M}_{3}$. Note, that the model predicts rather high numbers of ions inside the channel (between 3 and 4 for M₃) but there is no experimental evidence to support these specific occupancy numbers.

One of the most striking consequences of the Coulomb blockade oscillations is that, for channels having ${Q}_{f}\approx {M}_{n}$, adding one negative charge should dramatically decrease the Ca²⁺ conductance. So far, all mutation chains showed increase of calcium selectivity with growth of negative Q_f which can also, however, be explained by alternative models [32]. On the other hand these experiments are limited to $| {Q}_{f}| \leqslant 4e$ where the Coulomb blockade model also predicts an experimental increase of calcium selectivity with growth of $| {Q}_{f}|$ in the range $| {Q}_{f}| =0-4e$ (see figure 11 of [58)] but predicts a sharp drop in calcium selectivity near the neutralization blockade point Z₂.

We can suppose that the NavAB/NaChBac EEEE locus has an effective ${Q}_{f}={Z}_{2}=-4e$, and hence a 1e mutation of the EEEE ring (or a $\pm 1e$ mutation of the neighbouring ring of residues) should lead to increase of S for both directions of mutation. The currently available data (see tables (d), (e)) were obtained for mutation steps of ${Q}_{s}=4e$ and thus cannot resolve the effect of the smaller 1e jumps.

These identifications can be seen as an initial verification of the Coulomb blockade model predictions, albeit with some discrepancies. Further investigations are needed to confirm the channel identifications, to understand why the nominal Q_f is systematically slightly smaller than the Q_f values of band maxima in the model, and to establish the reason why mammalian and bacterial channels with EEEE-loci have opposite selectivity properties.

We now develop a description of the (interconnected) shapes of the Coulomb blockade oscillations in J, and the Coulomb staircase in P. We consider the case of divalent Ca²⁺ bands arising from strong Coulomb blockade (figure 2).

The equilibrium distribution of P_c in the vicinity of the M_n points (and hence calculation of the shapes of P_c(U) or ${P}_{c}({Q}_{f})$) follows from standard Coulomb blockade theory. As mentioned above, the energy level separation for divalent Ca²⁺ is so large $(\approx 40)$ that the general set of eigenstates (3) reduces to a simple two-state exclusive set:

$$\begin{eqnarray}&&m=\{n,n+1\};\quad {\theta }_{n}+{\theta }_{n+1}=1;\quad {P}_{c}=n+{\theta }_{n+1}.\end{eqnarray} \tag{ 9 }$$

The electrostatic exclusion principle (9) plays the same role here as that of the Pauli exclusion principle in quantum mechanics [83–85]; the standard derivation via the partition function, taking account of exclusion principle (9) leads [86–88] to Fermi–Dirac statistics for θ_n+1 and an excess (fractional) occupancy ${P}_{c}^{*}={P}_{c}\mathrm{mod}1$ :

$$\begin{eqnarray}&&{P}_{c}^{*}={\left(1+\mathrm{exp}\left({U}_{n+1}-\mu \right)\right)}^{-1}={\left(1+{P}_{b}^{-1}\mathrm{exp}\left({U}_{c}\right)\right)}^{-1},\end{eqnarray} \tag{ 10 }$$

where $\mu ={\mu }_{0}+\mathrm{ln}{P}_{b}$ is the chemical potential, P_b is a reference occupancy related to the bulk concentration, and μ₀ is a constant potential assumed here to be ${\mu }_{0}=\langle U\rangle =({U}_{n}+{U}_{n+1})/2$. Hence:

$$\begin{eqnarray}&&{U}_{c}={U}_{n+1}-{\mu }_{0}=({U}_{n+1}-{U}_{n})/2.\end{eqnarray} \tag{ 11 }$$

The Fermi–Dirac equation (10) is equivalent to the Langmuir isotherm [86] and to Michaelis–Menten saturation. A similar Fermi function was obtained earlier [51] for the variation of P_c with concentration.

Note that the Fermi–Dirac distribution needed for cases where the exclusion principle is based on volume exclusion, and the ions have unequal diameters, has been investigated by Liu and Eisenberg [85, 89, 90].

It follows from (5) that U_c is a linear function of Q_f:

$$\begin{eqnarray}&&\quad {U}_{c}={k}_{z}\displaystyle \frac{\Delta {Q}_{f}}{e},\end{eqnarray} \tag{ 12 }$$

where $\Delta {Q}_{f}={Q}_{f}-{M}_{n}$ and ${k}_{z}={{zU}}_{s}$. Hence the Fermi–Dirac function (10) also describes the dependence of P^* on Q_f:

$$\begin{eqnarray}&&{P}_{c}^{*}={\left(1+{P}_{b}^{-1}\mathrm{exp}\left({k}_{z}\displaystyle \frac{\Delta {Q}_{f}}{e}\right)\right)}^{-1},\end{eqnarray} \tag{ 13 }$$

where for our geometry with z = 2 (Ca²⁺ ions) the dimensionless scaling coefficient ${k}_{z}\approx 10$. This is the final result for the shape of the Coulomb staircase of occupancy as a function of Q_f. Figure 2(b) shows qualitative agreement of $P({Q}_{f})$ shapes with (13), including the concentration-dependent shift between curves with different P_b. We will make a more detailed comparison below.

Approximate forms of the current as a function of energy (or fixed charge) can be found via the variance σ of the Fermi–Dirac occupancy P due to thermal fluctuations, as follows from the fluctuation–dissipation theorem and linear response theory [91]. The ability of an energy level to contribute to the current/conductance is then proportional to ${\sigma }^{2}={\rm{d}}P/{\rm{d}}{U}_{c}$ via the Landauer approximation [54, 55]:

$$\begin{eqnarray}&&\displaystyle \frac{{J}_{c}}{{J}_{\mathrm{max}}}\propto \displaystyle \frac{{\rm{d}}P}{{\rm{d}}{U}_{c}}={\mathrm{cosh}}^{-2}\left(\displaystyle \frac{{U}_{c}}{2}\right),\end{eqnarray} \tag{ 14 }$$

where J_max is the barrier-less diffusive current. Taking account of scattering, one reaches the standard Coulomb blockade theory approximation [55, 92]:

$$\begin{eqnarray}&&\displaystyle \frac{{J}_{c}}{{J}_{\mathrm{max}}}={U}_{c}{\mathrm{sinh}}^{-1}\left({U}_{c}\right)\approx {\mathrm{cosh}}^{-2}\left(\displaystyle \frac{{U}_{c}}{2.5}\right).\end{eqnarray} \tag{ 15 }$$

Alternatively (15) can be derived by the quasi-equilibrium (or nonequilibrium reaction rate [93]) approach with explicit solution of the Nernst–Planck equation (i.e. the Goldman–Hodgkin–Katz (GHK) solution ) taking account of the Fermi–Dirac occupancy (10), so (15) can thus be called the Fermi-GHK approximation; a similar result was obtained earlier by Mott [94, 95].

Figure 4 reveals a resonant conductivity as U_c is varied, for both the Landauer (14) and Fermi-GHK (15) approximations: in each case there is a peak coinciding with the maximum in the derivative of P_c, ${\rm{d}}P/{\rm{d}}{U}_{c}$. (Note that, from (5), variation of U_c is equivalent to variation of Q_f.) In practical terms, the difference between the two approximations is small: they both represent double-exponential peaks of half-width ${U}_{1/2}\approx 2.3$ . The form of this current is similar to that of the tunneling current in a quantum dot [54]: an even, double-exponential function of U_c, reflecting the symmetry of the escape and relaxation trajectories [96]. Note that even small asymmetries, easy to overlook, can have disquieting effects [97].

Zoom In Zoom Out Reset image size

For a quantitative comparison of the theory with P_c as obtained from the Brownian dynamics simulations, we calculate the effective (excess) well depth ${U}_{c}^{*}$ as:

$$\begin{eqnarray}&&{U}_{c}^{*}=\mathrm{ln}\displaystyle \frac{1-{P}_{c}^{*}}{{P}_{c}^{*}}={k}_{z}\displaystyle \frac{\Delta {Q}_{f}}{e}.\end{eqnarray} \tag{ 16 }$$

The Fermi–Dirac function (10) predicts that ${U}_{c}^{*}$ should be linear in Q_f, i.e. (16) represents linearising coordinates for the Fermi–Dirac equation (10); the Coulomb blockade model also predicts the geometry-dependent coefficient k_z.

Figure 5(a) shows a sawtooth-like dependence of ${U}_{c}^{*}$ on Q_f, confirming that the ${P}_{c}^{*}$ transitions obey the Fermi–Dirac function (10) of U_c. For M₁, the slope corresponds well with analytic k_z; the discrepancies in slope at M₂ and M₃ are attributable to our neglect of ion–ion interactions.

Zoom In Zoom Out Reset image size

Figure 5(b) compares the predictions of the Coulomb blockade model with the conduction bands M₀, M₁, M₂ obtained from the simulations [57]. The Landauer and Fermi-GHK peaks are calculated using (14) and (15) respectively with the values of ${U}_{c}^{*}$ taken from plot (a), and there are no adjustable parameters. The BD peak shapes and positions are described reasonably well by the model, extending the simpler fitting in [58, 59]. Again, the agreement is very good for M₁, and less good for M₂ and M₃ on account of our neglect of field leaks and the singular parts of the interactions.

In general, the results of Brownian dynamics simulations [57, 58] correspond well with the predictions of the Coulomb blockade model.

Comparison of the Ca²⁺ (z = 2) bands pattern (figure 2) with the Na⁺ (z = 1) picture (figure 3 ) clearly shows that, unlike its electronic counterpart, ionic Coulomb blockade is valence-dependent: the positions of the bands ${{\rm{M}}}_{0}={ze}/2$, and their period ${ze}$, shift in proportion to z as described by equation (8) and they broaden/narrow in proportion to z² (5) [58].

The ionic Coulomb blockade model attributes selectivity to modification of the U_G dependences on Q_f. Valence selectivity is provided by the electrostatic z-dependent shift of ${U}_{G}(n,{Q}_{f})$ curves and the corresponding shifts of the M_n and Z_n singular points. Alike-charge selectivity can in principle be accounted for by hydration-dependent shifts of M_n and Z_n.

The model also predicts that a monovalent (e.g. sodium) ionic current can be effectively blocked by a few divalent (e.g. calcium) ions, an effect well-known in calcium and sodium channels. Calcium blockade is an aspect of AMFE, which is a signature of EEEE calcium channels [14]. The threshold of the divalent block IC₅₀ is defined as the Ca²⁺ concentration which results in a 50% decrease in the monovalent current [14]. We can estimate IC₅₀ on the assumption that the sodium current is proportional to ($1-{P}_{\mathrm{Ca}}$), i.e. to the fraction of time when the channel is not blocked by calcium ionc. Hence equations (12, 13) result in an exponential dependence of IC₅₀ on Q_f:

$$\begin{eqnarray}&&{{IC}}_{50}\propto \mathrm{exp}\left({k}_{2}\displaystyle \frac{\Delta {Q}_{f}}{e}\right),\end{eqnarray} \tag{ 17 }$$

which is a prediction that can be tested.

Mutant studies of the DEKA sodium channel [22, 24, 26] have measured the dependence of AMFE on the selectivity filter locus and its Q_f, and it was shown that log(IC₅₀) can successfully be fitted by a linear function of Q_f (figure 6), in agreement with the prediction (17) of the Coulomb blockade model, or

$$\begin{eqnarray}&&\mathrm{log}({{IC}}_{50})=a+b\displaystyle \frac{{Q}_{f}}{e},\end{eqnarray} \tag{ 18 }$$

where a and b are constants. The best-fit line gives b = 2.29 [24] which does not agree well with the $b={k}_{2}\;\approx$ 10 calculated for our model. This discrepancy can be interpreted as evidence that the geometrical parameters (R and L) of the selectivity filter of a DEKA sodium channel differ from values taken in our model (R = 0.3 nm, L = 1.6 nm). To be compatible with experiments the selectivity filter would need to be shorter and/or wider. We can get good agreement for e.g. R = 3.5 nm and L = 0.4 nm (such a short selectivity filter would in fact correspond well to the model of [70]) for which k₂ = 2.28. Because the exact values are not yet known, our fitting can be used to estimate these parameters. Note, however, that the fitting (18) is not unique to the Coulomb blockade model but follows from the dominance of electrostatics for divalent ions [24].

Zoom In Zoom Out Reset image size

Although an ion moving inside a channel is a room-temperature classical system described by Newtonian dynamics, we see that it exhibits some quantum-like mesoscopic features [52]. They include the appearance of strong Coulomb blockade oscillations in conductance, a Coulomb staircase, and a Fermi–Dirac distribution of occupancy for Ca²⁺ ions. We attribute such behavior to the charge discreteness, to the strong electrostatic interaction in confinement [52] and to the electrostatic exclusion principle. It has been shown rigorously that, in the presence of an exclusion principle, Brownian motion leads to a Fermi–Dirac distribution of the Brownian particles [84].

A parametric study [57] based on Brownian dynamics modelling has shown that, in accordance with (7), strong Coulomb blockade may be expected in channels of radii $R=0.25-0.35$ nm for ions having z = 2, e.g. for Ca⁺ ions in calcium/sodium channels. It is an interesting and significant question whether or not strong blockade and the corresponding oscillations may also arise in the conduction of monovalent ions in K⁺ channels. The latter have a selectivity filter of appropriate radius, $R\approx 0.25$ nm [34]. Furthermore, K⁺ ions are fully hydrated, which should lead to a stronger electrostatic interaction and hence to a marked decrease in the effective _w, bringing it close to ${\varepsilon }_{w}=1$. In such a case equation (7) suggests that we could expect an observable Coulomb blockade effect, even for monovalent ions in potassium channels [40, 41], a speculation that needs to be tested.

The relationship of the stop bands in the model to the subconductance states seen in experiments, and to the noise of closed channels, remains to be determined. Subconductance states and unusual baseline noise are found in a wide variety of natural biological channels [98].