Theory and Experiments on Multi-Ion Permeation and Selectivity in the NaChBac Ion Channel

The highly selective permeation of ions through biological ion channels is an unsolved problem of noise and °uctuations. In this paper, we motivate and introduce a non-equilibrium and self-consistent multi-species kinetic model, with the express aims of comparing with experimental recordings of current versus voltage and concentration and extracting important permeation parameters. For self-consistency, the behavior of the model at the two-state, i.e., selective limit in linear response, must agree with recent results derived from an equilibrium statistical theory. The kinetic model provides a good ¯t to data, including the key result of an anomalous mole fraction e®ect . bi ; and only considers the transition from 0 to 1 ions in the ¯lter. We observe a conduction peak, forming in the selective limits for either species.


Introduction
Biological ion channels passively (and stochastically) transport ions through the impermeable cell membrane, but precisely how they do it is a long-standing problem [1][2][3]. Many families of channels exist which can be further subcategorized by their primary conducting ions and gating mechanisms. In this work, we shall focus on a family of narrow voltage-gated channels [4]; and in particular on the prokaryotic Na þ channel NaChBac, from Bacillus halodurans [5]. This represents a model channel providing insight into mammalian voltage-gated channels which play an important role in numerous physiological processes, for example the propagation of the action potential; and have many interesting permeation properties such as high selectivity and blocking phenomena [1][2][3].
The Anomalous Mole Fraction E®ect (AMFE) is a good demonstration of blocking in multi-ion pores. A typical experiment involves considering a pair of conducting species, where one bulk solution is held constant while the ionic concentrations are varied in the other bulk solution. The outcome is a reduction in conductance of the mixed ionic solutions versus pure solutions [6]. It is generally assumed that these channels must conduct ions in single-¯le through multiple binding sites; however, AMFE has also been observed in simulations involving a single site [7]. AMFE can observed in many channels, including but not limited to: K þ and Ca þþ channels [2,[8][9][10] and NaChBac [11].
Prokaryotic channels are an attractive choice because they exist in a simpler protein structure whilst maintaining the permeation properties of their eukaryotic counterparts. The¯rst prokaryotic voltage-gated Na þ channel identi¯ed, was NaChBac in 2001 [5]. A crystal structure is not yet available, but homology models of suggested structure exist (see Fig. 1) [12,13]. These suggest four main binding sites, with two sites playing the dominant role and a knock-on conduction mechanism for Na þ . The selectivity¯lter (pore) has a length L c ¼ 14 # A and average radius R c ¼ 2:8 Å; geometrically constraining the pore to 4 or 6 K þ or Na þ ions in single¯le, respectively. The two dominant sites are located at the level of the carbonyl oxygen Fig. 1. Left, structure of NaChBac [12], made using Chimera [14]. The protein (purple ribbons) is embedded into a lipid membrane (orange lines) and solvated on either side by a NaCl solution (blue and green spheres). Right and top, zoomed in view of the selectivity¯lter representing the pore (highlighted by the yellow box) and right bottom, the lattice approximation of the pore. groups provided by the leucine residues (site 1) and the charged glutamate ring (site 2). The¯rst site has a diameter $ 5:5 # A and therefore is wide enough to accommodate the full primary hydration shell of a transiting Na þ ion. Thus, Na þ is prevented from interacting directly with the dipolar charge from the oxygen atoms. Site 2 is narrower with a diameter of $5 Å; and therefore, the Na þ ion su®ers partial dehydration of its¯rst shell. This energetic cost is o®set by direct interaction with the glutamates. In each case, the depth of the free energy well is almost identical.
A common approach to modeling such stochastic systems involves Eyring rate or kinetic theory. One generally assumes that the pore can be represented by a series of energy wells and barriers, and that occupancies (or states) can be described by a set of master equations from which the net°ux can be calculated. Typically these reproduce saturation of current versus concentration, as observed experimentally [15][16][17], but also more interesting properties such as: the importance of volume exclusion [18], the role of selectivity and interactions between ions [16,19,20] and the consequence of mutagenisis [21]. Generally, this approach faces a persistent criticism over the validity of the transition rates [22][23][24], their relation to structure [25] and self-consistent inclusion of concentration into the energy barrier [26]. To counter this, e®orts have been made to de¯ne rates rigorously using mean¯rst passage time (mfpt) theory [26,27] or equilibrium statistical theory [28].
In this paper, we aim to overcome these issues by developing a multi-species kinetic model within the theoretical framework introduced in earlier work [29,30]. We have the explicit aims of validating the theory by comparison with experimental recordings for NaChBac including the blocking phenomena AMFE, and of explicitly calculating useful parameters including the energy barriers at each binding site and the e®ective di®usion coe±cient in the pore. We start by de¯ning the total system, the state space, the corresponding energy spectra and the statistical ensemble. The statistical°uctuations are then related directly to the°ow of ions at linear response. The model is then extended far from equilibrium by a set of master equations where the transition rates are de¯ned by comparison of the linear responses with equilibrium behavior.
In the work that follows, with SI units: q; k; T , respectively, represent the proton charge, Boltzmann's constant and the system temperature. We use the following bulk di®usion coe±cients: D b K ¼ 1:96 Â 10 À9 m 2 s À1 and D b Na ¼ 1:33 Â 10 À9 m 2 s À1 , for K þ and Na þ , respectively.

Discrete system
To de¯ne the system let us consider a pore, thermally and di®usively coupled on either side to bulk reservoirs [29]. Under typical experimental or physiological conditions the thermal de Broglie wavelength Ã of e.g., the ion Na þ has Ã $ 20 pm. We note that this is much less than the average distance between ions in the bulk and even in narrow pores. Hence to a reasonable approximation the system is classical, and therefore the ions obey Boltzmann statistics. The physical properties of the system are determined by the canonical ensemble assuming constant total volume V and temperature T , and conservation of the total number of particles N i of each species i. However, we note that, in general, the van der Waals forces (that are of quantum origin) still play a fundamental role in ionic permeation and these will be included in future work. We approximate the ion channel as the selectivity¯lter, which we treat as a 1D lattice (see Fig. 1) with a total of M binding sites. Each site in the lattice can hold at most 1 ion and the bulk reservoirs contain mixed solutions of arbitrary concentration. Since the lattice has a¯nite number of sites, its occupancy represents a state of the system, and the total set of states is denoted by fn j g. Each state is then characterized by the set of occupation numbers, for each species i 2 1; . . . ; S at all individual sites m 2 1; . . . ; M labeled from left to right: fn im g ¼ ½n i1 ; n i2 ; . . . ; n iM where P S i¼1 n im 2 0; 1 and P M m¼1 P S i¼1 n im M.

E®ective grand canonical ensemble
In earlier work [29,30] we demonstrated that such a system can be formulated within the e®ective grand canonical ensemble, with the following Gibbs free energy, Here, Eðfn j g; n f Þ denotes the electrostatic contribution which includes the ion-ion and ion-pore long range interactions. The indistinguishability of ions and empty sites is included via W ðfn j gÞ, which manifests via the natural logarithm of the factorial of total number of occupying ions P i n i ! and the factorial of the number of empty sites n es ! (which is equal to ðM À P i n i Þ!). The¯nal term includes contributions from the excess chemical potential di®erence between the bulk and binding site Á" b im , and the natural logarithm of the bulk mole fraction The excess chemical potential di®erence describes the di®erence in non-ideal local interactions between the pore and site, and is largely dominated by the dehydration energy of the ion [31,32].
Within this ensemble (see [29,30] for full details) we can calculate the current I and conductivity of the favored ion in the high selective limit, which is shown to be proportional to the variance in the particle number.
Here, i and c are the electro-chemical and chemical potentials in the pore.

Kinetic theory
To extend the theory far from equilibrium, we introduce a set of master equations describing transitions within our state space. Importantly, this approach requires the speci¯cation of transition rates À, which must be functions of the energy barrier to enter/exit each state. We have seen that each site which can be occupied by at most one ion is described by a transition number n m equal to the sum over all possible species i at the site. The total network can be written in matrix form P : ¼ PÀ, with the row describing each state given by, The right-hand side describes the transiting ion n Ã i , either entering the mth state, or exiting and entering one of the m AE 1 states. For clarity, the conditions H are imposed which require that the state of the pore must remain within the state space and the transitions be physically possible i.e., an ion can only enter a neighboring and empty site, and an ion can only exit from a site that it is already occupying. This set of equations is similar to those commonly used to describe electron°ow in quantum structures [33,34], with the important di®erence that here we do not allow transport into any empty site but instead specify via the conditions H.
Since the time-scale of permeation is short $ 10 ns, it is safe to neglect temporal behavior and consider _ P ¼ 0. The steady-state current can be de¯ned for all possible transitions involving each binding site, where Kircho®'s laws are satis¯ed. Hence the current at each site can be de¯ned from the balance of°uxes, such that its positivity de¯nes ionic°ow from left to right. Hence, the current into the mth site can be written as where we note that the total current can be calculated by summing over all possible states for transitions to the¯rst site. In the equilibrium limit, the transition rates must be constrained by the condition of detailed balance; and the probabilities will reduce to the form given by equilibrium statistical physics [29,30,35].

Application to NaChBac and its Mutant
To describe the permeation of NaChBac and its mutant we employ a toy model with two binding sites; and the potential energy pro¯les shown in Fig. 2. Note that, in general, the depth of the well at either site does not have to be equivalent, and it may change depending on the state. Since anions are electrostatically repelled by the negative pore charge n f , for two conducting species e.g., Na þ and K þ we have nine possible states. We label these A À I for the purpose of notational clarity in Fig. 2 and Eqs. (5) and (6). We note that terms Á" b i depend on the species and hence result in energy level splitting (see Fig. 3).  To agree with the experimental conditions the concentrations of the two species always sum to 0.14 M and hence: x K ¼ ð0:14=c W À x Na Þ. The excess chemical di®erence is taken to be 0 with Á" K ¼ Á" Na ¼ 1 kT , and as such the¯lter is only selective based on the concentration. Along the n f plane levels are parabolic, minimizing when n f ¼ P i n i . At low Na þ concentrations Na þ faces a large concentration energy barrier with its barrier-less transition (the crossing of neighboring levels) occurring at a larger free energy point. As Na þ concentration increases the levels converge in the absence of selectivity, and start to invert beyond this point with Na þ now becoming favored. Right, normalized Na þ and K þ conduction peaks for adding a single ion to an empty pore versus the concentration dependent chemical potential i ¼ kT lnðx i Þ þ " i assuming that D c i ¼ 5 Á 10 À10 m 2 s À1 , E b T ¼ 0:5 kT , ÁE c T ;i ¼ 0:5 kT and ÁE þ " c im ¼ À5 kT . In the selective regions jÁ K À Á Na j ) 0 species speci¯c conduction peaks are observed, maximizing under the conditions ÁG $ 0. Away from these peaks ionic Coulomb blockade occurs and the current is zero [11,29,34,39]. The blocking e®ect AMFE can be observed qualitatively: the orange curve represents the experimental data points.
The corresponding set of master equations can be written as Here, the rates À are de¯ned such that the superscript de¯nes the order of transitions between states i.e., À AB implies transitions between states A and B (see (5)). Diagonal matrix elements are equal to the sum over all transitions involving the site, and the columns sum to zero; and thus we can replace one row with probability conservation. There are two transition types: (i) external i.e., involving a bulk and (ii) internal i.e., between sites. Before we discuss the general form of the rates we note two important properties. In perfectly planar and neutral membranes the transmembrane potential is linear across the pore [32,35], illustrated by the green line in Fig. 2. For simplicity, we shall de¯ne it as follows by taking L ¼ ; c m ¼ m and R ¼ 0, where m describes the relative position of the site m. Consequently, only the exiting rate is in°uenced by the work done in moving against this potential. From the condition of detailed balance we can de¯ne each rate using a functional normalization, which will be computed by comparison of the linear response with Eq. (2), The terms with subscript T represent the contributions from the transition states. We shall consider that the transition state barriers at the bulk-pore interface E b T ;i are identical at either entrance and for both species, although in general this may not be the case. We also introduce the parameter ÁG B im to represent the binding energy at each site, where the terms are as de¯ned in Eq. (1). The internal barrier ÁE c T ;im is given by the di®erence between the internal transition state level and the local interaction at each site: ÁE c T ;im ¼ E T ;i À " c im . To calculate these normalizations A=B, we reduce the system to two-states i.e., a single-site. We thus eliminate internal transitions and apply the condition of selectivity such that we consider each species seperately. By comparison of the linear response, we¯nd that the normalization of A must take the form where ÁG B im describing each site is given by (10). For symmetrical solutions these normalizations are identical at either bulk and reduce to a suppressing (Kramer's type) exponential function in the limit jÁGj ) 0. We allow these normalizations to di®er for asymmetrical solutions (without including the Nernst potential), yet it still reproduces the Kramer's limit if the local energy barrier ÁG B im becomes large. This form is similar to the creation/annhilation probabilities used within the grand canonical Monte Carlo formalism [32,36]. To maintain units and convention we shall also introduce the internal normalization

Comparison with experiment
To compare the theory with experimental recordings and make predictions we consider two experiments. For further details of the experimental methods, including generation of mutant channels and their expression, as well as details of electrophysiological experiments, we refer to [11], and here we only present a concise summary. In the¯rst series of experiments we performed whole-cell current measurements through a NaChBac mutant (LDDWAD) channel, which we assume to have the same conduction mechanism as wild type NaChBac, in di®erent Na þ /K þ concentrations (AMFE experiment). In these recordings the pipette solution contained (in mM) 15 Na-gluconate, 5 NaCl, 90 NMDG, 10 EGTA, and 20 HEPES, pH 7.4 (adjusted with 3mM HCl), meanwhile the bath solution contained (in mM); 137 NaCl, 10 HEPES and 10 glucose, pH 7.4 (adjusted with 3.6 mM NaOH). Permeability to K þ was investigated by replacing the NaCl bath solution with an equivalent KCl solution such that the total ionic concentration was¯xed at 140 mM. Total current across the cell was then normalized and, because one can assume that the total number of channels and their type is conserved in each cell for the duration of the recording, it can e®ectively be modeled as a single channel. Data were collected and the peak current magnitude in the current-voltage relationship was found at a À25 mV voltage drop across the channel. Consequently, there exists strong electrochemical gradients on the ions particularly when Na þ concentration exceeds 0.07 M (the magnitude of the gradient exceeds 2 kT ); and when there is an absence of ions in one side of the membrane e.g., when the bath solution is devoid of Na þ . The second experimental comparison is with single-channel NaChBac channel recordings using an identical bath and pipette solution containing (in mM: 137 NaCl, 10 HEPES and 10 glucose, pH 7.4 adjusted with 3.6 mM NaOH). Single-channel recordings are possible because Na þ is the preferred substrate with su±ciently high conductance as to provide a single-channel current amplitude which signi¯cantly exceeded noise (i.e., a favorable signal to noise ratio).

Linear response
To investigate the permeation process we shall consider the behavior of the kinetic equations under close-to-equilibrium conditions i.e., the linear response. The bulk solutions are considered symmetrical and take the form given by the bath solution in the experiment with a small voltage drop equivalent to 1 kT is applied. In Fig. 3 (left), we consider the energy spectra of the pure states only in the system versus¯lter charge n f and Na þ concentration where the K þ concentration is given by c K ¼ 0:14M À c Na . We also note that the binding interactions of the species are taken to be the same. The electrostatic interaction E is taken from [37,38], and is equal to U c Á ð P i n i þ n f Þ 2 , where U c is the charging energy. We note that for the NaChBac geometry it is relatively small and takes the value $ 6 kT . In the n f plane, the energy spectra form a series of parabolic curves, arranged according to ionic species and number of ions in the pore. The energy spectra can cross at particular values of n f , where the energy di®erence is barrierless ÁG $ 0. This corresponds to degeneracies in the state of the pore and results in resonant conduction of the favored ion (see Fig. 3); the disfavored ion faces an energy barrier which prohibits conduction. This corresponds to Eisenmann selectivity [29,30], and in this example it reduces to the di®erence in concentrations of the species. As the Na þ concentration increases, the selectivity barrier becomes smaller until 0:07 M and then changes sign. Consequently, the free energy spectra of the pure Na þ states become less than their K þ counterparts and the pore now favors Na þ . Figure 3 (right) displays the normalized conduction of the pore in the linear response regime with a small voltage gradient of 25 mV, and the condition of an identical di®usion rate in the pore. It is plotted versus the chemical potentials of both species i ¼ kT logðx i Þ þ " b i ; and only considers the transition from 0 to 1 ions in thē lter. We observe a conduction peak, forming in the selective limits for either species. Peaks maximize under the condition that ÁG $ 0 i.e., a barrierless transition which is governed by ionic Coulomb blockade [11,29,34,39]. The orange curve represents the experimental data from LDDWAD exhibiting the blocking e®ect AMFE. For clarity these are placed on the conduction peak. The qualitative agreement proves that the intersection of the two curves can reproduce this blocking e®ect. To improve agreement we would clearly need to¯t the data under non-equilibrium conditions and this will be undertaken in Sec. 3.3.

Non-equilibrium response
In order to compare with the experimental results we shall now use the kinetic model, implementing the exact experimental conditions. To reduce the¯tting parameters we shall introduce 1 ¼ 0:66, 2 ¼ 0:33 and c T ¼ 0:5, and constant bulk-pore transition state barriers of 1 and 0.5 kT in the NaChBac and AMFE comparison, respectively. The free parameters are the binding energy at each site ÁG B;b im , the internal energy barrier ÁE c T ;im and the parameter . This is introduced via the relation D c i ¼ D b i and necessary because unlike in the bulk, the di®usivity in the pore is unknown. Simulations estimate this parameter to take values of $ 0:01 to 0.5 for pores of the widths considered [40]. In Fig. 4, we display the results of comparison with experimental recordings. In the left panel we illustrate a good¯t to the single-channel I À V data at a concentration of 0.14 M, and make predictions for the shape of the I À V curves at larger concentrations (0.4 M and 0.6 M) and at larger voltages. We¯nd that the pore is symmetrical with respect to voltage and calculate the energy barriers for the¯rst and second ion entering the pore to be: ÁG B Na ¼ 3:75 kT and ÁG B Na ¼ 2:75 kT , respectively, similar to values in [12]. The barriers decrease for the second entry, because the electrostatic attraction to the pore charge is likely to be reduced due to compensation of the total pore charge by the already present ion. Finally, we note that the di®usion coe¯cient takes a reasonable value with ¼ 0:5, and we predict in accordance with single channel currents that saturation occurs after $ 0:2 V [2].
In the right hand panel of Fig. 4, we consider the AMFE whole-cell experiment undertaken for the mutant (LDDWAD). For consistency in¯tting we have treated the values of normalized current as values in [pA], because it is a similar order of magnitude (although di®erent channel) to the single channel data. We are able to model the data quite accurately with the parameters given in the caption. The binding selectivity at each site (S m ¼ ÁG B Na;m À ÁG B K;m ) favors Na þ over K þ ; however, the channel is almost non-selective for the¯rst (m ¼ 1) site when adding Fig. 4. Left, theoretical (T ) current¯tted against single-channel experimental recordings (E) for NaChBac using the following¯tting parameters: ¼ 0:5, for the¯rst and second ions entering the pore respectively (in kT ): ÁG B Na ¼ 3:75 and ÁG B Na ¼ 2:75 and an internal barrier of ÁE c T ;i ¼ 0:5 kT . We note that the error bars were small and hence omitted. This comparison is only made for symmetrical 0.14 M solutions, and so we have predicted the current at the larger 0.4 M and 0.6 M NaCl solutions (orange and green curves respectively). Right, comparison of the¯tted theory (full curve) with the normalised whole-cell current measured for the NaChBac mutant LDDWAD (data points) in mixed K þ /Na þ solutions. The¯tting parameters were: ¼ 1, for the¯rst ion entering (in kT ) ÁG B Na;1 ¼ 4:7, ÁG B Na;2 ¼ 7:3, ÁG B K;1 ¼ 4:5, ÁG B K;2 ¼ 4:8, and for the second ion entering (in kT ) ÁG B Na;1 ¼ 8:2, ÁG B Na;2 ¼ 7, ÁG B K;1 ¼ 5:7 and ÁG B K;2 ¼ 3:7. The equivalent barriers involving mixed ionic states are calculated after subtracting kT lnð2Þ due to the di®erence in permutations term W (see (1)). Finally, the internal transition barriers were (in kT ): ÁE c T ;Na1 ¼ 0:1, ÁE c T ;Na2 ¼ 2:8, ÁE c T ;K1 ¼ 1:5 and ÁE c T ;K2 ¼ 1:8.
the¯rst ion. The total selectivity (given by the sum over m) di®ers between adding the¯rst and second ions taking the values þ2:8 kT and þ5:6 kT , respectively, suggesting that it is harder for K þ to enter if Na þ is already occupying the pore. It is also likely that the conduction mechanism and number of binding sites will di®er in this mutant in addition to protonation e®ects [11]; and so further analysis will be required.

Conclusion
In summary, we have presented a far-from-equilibrium multi-species kinetic theory explicitly tracking the permeating ions through two binding sites. It is directly related to the equilibrium statistical theory, brie°y discussed here but introduced in more detail elsewhere [29,30]; through the calculation of the transition rates and in particular their normalization.
To investigate the permeation of wild type NaChBac and a mutant we have compared the predictions of the kinetic model with two sets of experimental data (see Fig. 4). The¯rst of these contains single-channel I À V curves from which we are able to extract the e®ective Na þ di®usion coe±cient in the pore D c Na $ 6:65Â 10 À10 m 2 s À1 , and the binding energy parameters (ÁG B Na;m ) of 3:75 kT and 2:75 kT for adding the¯rst and second ions, respectively. These are identical in both sites indicating a symmetrical pore. Comparison with the second set of data (describing the mutant) enables us to recover the blocking e®ect AMFE. The barriers varied at each site and always remained selective to Na þ , although this selectivity was small in thē rst (m ¼ 1) site when adding the¯rst ion.
In future, we aim to further compare NaChBac and its mutants by extension of the number of binding sites, analysis of the e®ects of mixed-valence i.e., Na þ /Ca þþ and the e®ects of changing the structure. The enhanced conduction behavior observed within the kinetic theory is also a topic for further research. Finally, we comment that we expect our model to be applicable to the permeation of other voltage-gated ion channels and arti¯cial nanopores.