Molecular Na-channel excitability from statistical physics

Experimental understanding of biophysical electrical processes in the cell membrane during the action potential has much progressed during the last 60 years, mainly fostered by the seminal works by Hodgkin and Huxley [1]. They performed extensive experiments on the giant squid axon, and constructed a mathematical model that has constituted since then the basis for the interpretation of the behaviour of nerve and cardiac cells. According to the Hodgkin-Huxley (HH) model, action potential is produced by coordinated ionic fluxes crossing the cell membrane, which acts as a capacitor. Then the excitable characteristics of the membrane action potential is the result of a dynamical coupling among the ionic flux, the membrane conductance and the electrostatic potential [1].

It is now also well known that ions flow along some biochemical molecules (channels) embedded in the cell membrane. These channels present two main structural conformations (open and closed), with transitions between these two states controlled by the membrane potential. Action potential is then known to be the result of the synchronized dynamics of a large number of ionic channels [2]. Moreover, much quantitative physical information is also known on the dynamics of the distinct states of single ionic channels [3–6]. In particular experiments on single channels show very strong fluctuations in the intensity (pA, i.e. a few charges in a microsecond) crossing the channel. As a result the observed stochastic behavior has become an active topic in recent studies.

Several theoretical scenarios have been used to address this stochastic phenomenology. Most of these approaches incorporate fluctuations in some of the elements of the HH theory, for instance by using Langevin [7–10] or master equations [11] for the conductance equations or noise terms in the equation for the membrane charge [12]. A more microscopic approach used Langevin equations for the ions with Poisson equation for the potential membrane [13] in order to obtain the effects of fluctuations on the membrane conductivity.

A microscopic modeling of the excitable dynamics of the action potential, treating channels and ions as physical objects, merits attention. This approach would allow for studying single-channel excitable events, and to obtain additional information on some aspects of the action potential dynamics, such as energetic balances. It would also provide the influence of changes of different physical parameters (concentrations, temperature, etc.) on the whole process, without the need of additional parameter fittings. Then one could address questions such as that concerning the minimum elements necessary to produce the action potential or whether the cooperative coupling of a large number of channels is necessary for excitability.

Our aim in this work is thus to place the action potential spikes within the framework of statistical physics to explain these phenomena at the microscopic level. We will propose a minimum model presenting the desired excitable behavior, and accordingly we will not try to get quantitative agreements with any particular cell type. As a result we expect our simplified model to be close to a very primitive channel, presumably much simpler than in modern organisms and containing the minimum set of elements that permit a bona fide excitability behavior.

In this paper we proceed first with the description of the Na channel and the physical elements that constitute our model. Next we present the numerical results for the excitable behavior of a single Na channel, and relate the excitable properties of the model to the physical mechanisms implicit in the classical HH theory. We will show that the minimum scenario to explain the excitable properties is a single gated Na channel in the presence of a leakage of K ions. We end with some conclusions and perspectives.

Our approach consists of treating ions as Brownian particles, and channels as physical pores with mechanical doors that have two steady states (open and closed). All of them are driven by the membrane potential which also depends on the ionic flux. We will focus on the dynamics of a Na channel with two doors, whose states are defined by the variables Y₁ and Y₂, that will evolve according to their respective dynamics controlled by the membrane potential $\Delta V$. Our approach is then closely related with that of ref. [13] but within a more simple scenario. The leak of K ions through the membrane will be modeled as an additional channel with effective parameters without doors (see fig. 1).

Another relevant point for a microscopic description is that fluctuations should be relevant locally due to the very small number of charges involved and the fact that, although they move deterministically under the electrostatic force, they diffuse also by thermal noise. Following standard formulations of nonequilibrium statistical mechanics, the main variables of the model follow overdamped Langevin equations with their corresponding potential energies and thermal noises. The system is autonomous, and the only source of energy is the Gibbs energy associated with the ionic concentrations and the membrane potential. In this way the model is also able to give information about the energetics of any excitable event. In this formulation, the model parameters and other characteristics can be related to biological experimental information.

Zoom In Zoom Out Reset image size

To make explicit our microscopic approach we will formulate our approach according to the following assumptions:

• i)  
  The first assumption is the use the capacitor equation for the membrane potential $\Delta V$,

  $$\begin{equation} -C_{\text{M}} \frac{\text{d} \Delta V}{\text{d} t} = I_{\text{Na}} + I_{\text{K}} + I_0, \end{equation} \tag{ 1 }$$

  where $C_{\text{M}}$ is the membrane capacity assumed to be constant and adjusted to a single Na channel. The intensities $I_{\text{Na}}$, and $I_{\text{K}}$ correspond to the flux of Na and the leak of K ions. I₀ is a perturbative current pulse that will trigger the spike. This equation will be integrated as

  $$\begin{equation} \Delta V(t+\Delta t) =\Delta V(t) - \frac{\Delta Q(\text{Na}^+) + \Delta Q(\text{K}^+)}{2 C_{\text{M}}}, \end{equation} \tag{ 2 }$$

  where $\Delta Q(\text{Na}^+)$ and $\Delta Q(\text{K}^+)$ are the balance of charges crossing any of the channel boundaries during the interval of time $\Delta t$ which are obtained from the trajectories of ions. The divisor "2" takes into account that charges cross both boundaries of the channel (then being counted twice) when crossing the membrane from one side to the other one.
• ii)  
  The second assumption has to do with the calculation of the ionic flux. Ions inside the channel are described by point-like particles with electrical charge $+q$ moving in one dimension. Their positions $x_i(t)$ obey overdamped Langevin equations, $\gamma_i {\dot x_i} = - \partial_{x_i} U + \xi_i(t)$, where $\gamma_i$ is the effective friction and $\xi_i(t)$ is a thermal noise of zero mean and intensity $\gamma_i k_{\text{B}}T$. The interaction potential U is the addition of the interactions with the doors (see below) and the electrostatic membrane potential $V_{\text{e}}(x_i,\Delta V)$,

  $$\begin{equation} V_{\text{e}}(x_i,\Delta V)= q \frac{\Delta V}{L} ( x_i- L), \quad 0 < x_i < L. \end{equation} \tag{ 3 }$$

  This Langevin equation has to be complemented with boundary conditions of concentration values $\rho_0 = A c_{\text{in}}$ at x = 0 and $\rho_2 = A c_{\text{out}}$ at x = L, being A the effective section of the channel and $c_{\text{in/out}}$ the bulk ion concentration, interior and exterior to the cell respectively (note that any ion affinity of the pore could be accounted for by changing the value of A in these relations). Boundary conditions are implemented in the following way: ions disappear when hopping out of the channel due to their Brownian motion, and they appear into the channel according with a probability depending on the concentration at this boundary. There is no need of any explicit assumption about the form of the conductances, but nevertheless the fact that it is formulated consistently with statistical mechanics guaranties that this model evolves towards the correct steady-state membrane potential without further parameters or fine tuning. That is, it provides the Nernst potential when a single-ion species can cross the membrane, and the Goldman-Hodgkin-Katz theoretical prediction [14] when different ions compite. Note also that we are explicitly neglecting any ion-ion interaction inside the channel. This is justified by the small number of ions simultaneously present in the system and the screening of the aqueous medium.
• iii)  
  The dynamical equations for the channel doors are the kernel of our approach. There is strong experimental evidence that the Na channel has two active doors or barriers [2], and that they open and close stochastically according to the value of the electrostatic membrane potential and thermal fluctuations [6]. This hypothesis and the use of Langevin equations are the original parts of our approach. Then we describe a channel door as a physical barrier controlled by the dimensionless variable Y which behaves as a nonlinear spring with two steady states: $Y \sim 0$ (closed) and $Y \sim 1$ (open). These door states are controlled by the elastic potential $V_j(Y_j, \Delta V)\ (j=1,2)$, given by

  $$\begin{eqnarray} V_j(Y_j, \Delta V)&= V_0 \left[ -a \ln (Y_j(1-Y_j)) - b (Y_j-0.5)^2 \right] \nonumber\\ &\quad +\, Q_j(\Delta V - \phi_{\text{ref}}) Y_j, \end{eqnarray} \tag{ 4 }$$

  where Q_j is the charge of each door sensor and $\phi_{\text{ref}}$ is a reference potential. Their values are specific of each door, $Q_1= 12\ \text{e}$, $Q_2=-8\ \text{e}$, whereas we take common values for the other parameters: $V_0= 7 k_{\text{B}}T$, $a=0.2$, b = 9, and $\phi_{\text{ref}}=-40\ \text{mV}$. This potential presents two minima near $Y_i\sim 0$, 1 corresponding to the closed and open states respectively. These minima interchange their relative metastability by changing the value of $\Delta V$. For smaller voltages the door 1 (2) is in the closed $Y_1 \sim 0$ (open $Y_2 \sim 1$) state, and for larger voltages we have the oposite behavior. At $\Delta V = -40\ \text{mV}$ both states are equally probable in both doors.

Since barriers are physical entities, when ions interact with them they interchange momentum and energy. Thus variables Y and x_i have to obey physical laws expressed in terms of dynamical (Langevin) equations constructed from a common potential. With this requirement the potential $V_{\text{I}}(Y,x_i)$ corresponding to the interaction between particles and internal barriers is modeled as

$$\begin{equation} V_{\text{I}}(Y,x_i) = V_{\text{d}}f(Y) \exp \left(-\frac{(x_i - x_{\text{c}})^2}{2 \sigma^2} \right), \end{equation} \tag{ 5 }$$

where V_d is the barrier height, $x_{\text{c}}$ is the position of the barrier center inside the channel and σ is its width (see fig. 2). The function f(Y) modules the aperture of the doors according to the Y_j variables. The parameter values for the two doors have been taken as $V_{\text{d}}(1)=200\ \text{meV}$, $V_{\text{d}}(2)=250\ \text{meV}$, $x_{\text{c}}(1)= 1\ \text{nm}$, $x_{\text{c}}(2)= 3\ \text{nm}$, and $\sigma= 0.283\ \text{nm}$.

Zoom In Zoom Out Reset image size

For the modulated function f(Y) in eq. (5) we have taken the function,

$$\begin{equation} f(Y) = \frac{1}{2}(1 + \cos \pi Y ), \end{equation} \tag{ 6 }$$

which has the values $f(0) = 1$ for the closed state, and $f(1)= 0$ for the open state, and has relative extrema at these points. This property reduces the sensitivity against thermal fluctuations of Y around the steady states.

Regarding the K leakage through the membrane, we consider the motion of K ions as equivalent to moving in an additional (K-selective) channel without any door, and with effective parameters. This provides a charge leakage that restores the membrane potential at the end of the action potential. Then we have not considered for K a gated channel in the spirit of the HH theory (see below), since we are seeking a minimal model and as we will show such a gate is not necessary for excitability.

According to the former assumptions our approach has a set of equations that need to be numerically simulated. Our variables are the position x_i of the ions (Na and K) inside the channel, the Na channel doors Y₁ and Y₂ and the membrane electrostatic potential $\Delta V$.

The whole system can be characterized by the potential energy,

$$\begin{eqnarray} &U(x_i,\Delta V,Y_1,Y_2) =\nonumber\\ &\sum_i V_{\text{e}}(x_i,\Delta V) + \sum_{i,j} V_{\text{I}}(Y_j,x_i) + \sum_j V(Y_j,\Delta V),\quad \end{eqnarray} \tag{ 7 }$$

and accordingly the set of Langevin dynamical for our mechanical variables is

$$\begin{eqnarray} \gamma_i {\dot x_i} &= - \partial_{x_i} U(x_i,\Delta V,Y_1,Y_2)+ \xi_i(t), \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \gamma_{Y_1}{\dot Y_1} &= -\partial_{Y_1} U(x_i,\Delta V,Y_1,Y_2) + \xi_{Y_1}(t), \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} \gamma_{Y_2}{\dot Y_2} &= - \partial_{Y_2} U(x_i,\Delta V,Y_1,Y_2) + \xi_{Y_2}(t), \end{eqnarray} \tag{ 10 }$$

where thermal noises fulfill

$$\begin{equation} \langle \xi_a (t) \xi_b(t') \rangle = 2 \gamma_a \,k_{\text{B}} T\, \delta_{a,b}\, \delta (t- t'). \end{equation} \tag{ 11 }$$

Note that the first Langevin equation is for all ions: both Na and K. The simulation of these equations determines the state of the doors and the trajectories of the ions. The numbers of particles entering into and leaving the channels through each boundary are used to evaluate the potential $\Delta V$ through eq. (2). The final output is $\Delta V (t)$ which has to be compared with the known experimental results.

We have considered a single Na channel and the leak of K ions, and we have simulated the whole system of equations (2) and (8)–(11). The parameter values of the Na channel in table 1 have been selected to fulfill the experimental observations [6]. For the K leakage the effective parameter values are: $\gamma_{\text{K}^+} = 200\ \mu \text{s meV/nm}^2$, and $\rho_0^{\text{K}}, \rho_1^{\text{K}} = 20, 0.36\ \text{charges/nm}$, respectively.

Table 1:.  Physical parameter values used in the simulations for a single Na channel.

$\gamma_{\text{Na}^+}$ particle friction	$2\ \mu \text{s meV/nm}^2$
$\gamma_{Y_1}$ door 1 friction	$1000\ \mu \text{s meV/nm}$
$\gamma_{Y_2}$ door 2 friction	$4000\ \mu \text{s meV/nm}$
$k_{\text{B}}T$	$25\ \text{meV}$
L channel length	$4\ \text{nm}$
$\rho_0^{\text{Na}}, \rho_1^{\text{Na}}$	0.01, 1.2 charges/nm
$C_{\text{eff}}$ effective capacity	1.25 charges/mV

More specifically, as in a real experiment, we follow the dynamical evolution of the membrane potential when small and instantaneous discharges $\Delta Q$ of positive ions, corresponding to depolarizing voltage perturbations $\Delta V_0$ of $+80$ or $+70\ \text{mV}$, are applied to the membrane with a period of 5 ms. In fig. 3 we show a typical time interval with 10 of these perturbation events. In order to show the characteristics of the perturbations, we show, on top of this figure, how these pulses are seen when they are applied to the membrane without the presence of the Na channel, i.e. with only the K leakage. We see the expected response of the system as a sudden increase of $\Delta V_0$ followed by a slow relaxation towards the steady value of the membrane potential. We can also appreciate the size of the voltage stochastic fluctuations.

Zoom In Zoom Out Reset image size

In the middle and bottom figures we show the response of the system under these perturbations. The middle graph of the figure corresponds to perturbations equivalent to instantaneous increases of $\Delta V$ of $+70\ \text{mV}$, and the bottom graph to increases of $+80\ \text{mV}$. At each perturbation event the value of the potential membrane $\Delta V(t)$ presents narrow and larger excursions towards positive values. This high increase is due to the fast flux of Na ions into the cell when both channel doors are opened. Then the door 2 of this channel closes suddenly and the outward K flux starts to restore the initial steady state of the membrane potential but in a larger time scale. Although most of the peaks are real excitable events (their height is around twice larger than that of the perturbation), a few of them have some imperfections. In the middle panel we see some failed (f) or missing events when the Na channel door Y₁ does not open, and double peaks (d) when door 2 opens again before the closing of door 1. Also at the bottom graph we see small (s) pulses, in which the door 2 closes very fast and the channel has been active a very short time. One appreciate that for pulses of $+70\ \text{mV}$ (middle graph) the number of errors is larger.

To describe more explicitly the dynamics of the model during the action potential, we show in fig. 4 a detailed view of a single pulse (the 9th pulse in fig. 3 bottom) with numerical results for other observables. The top frame in this figure is an amplification of the membrane potential, the middle frame is the plot of the ionic intensities during the same pulse, and in the bottom part we find the evolution of the two Na channel doors, Y₁ and Y₂. In these frames we have marked five different times: t₀ is the perturbative trigger time, where the potential is instantaneously increased in an amount $\Delta V_0 = 80\ \text{mV}$. This is followed by the opening of door Y₁ at t₁. Then at t₂ the door Y₂ closes. In the interval $(t_1, t_2)$ both Na-channel doors are open and Na ions enter into the cell producing the rise of the $\Delta V$ pulse. This is manifest in the middle figure where we see the corresponding inward (negative) Na intensity. This interval corresponds to the so-called open state [2]. After t₂ the Na flux is stopped, due to the closing of Y₂, which corresponds to what is known as the inactivated state of the channel. Here an eventual additional perturbation would not induce any channel opening. Now K leak starts to dominate tending to restore the initial or standby state by an outward (positive) K intensity, as seen in the middle frame. Then at t₃ the Y₁ closes and at t₄ the Y₂ opens. The refractory time corresponds to the interval $(t_2,t_4)$ when Y₂ remains closed. After t₄ the channel recovers the steady (excitable) closed state. In this figure we can also see the fluctuations of the door variables and their almost instantaneous transitions following the membrane potential. In the middle frame data of ion intensities have been filtered by using an averaging filter with a window of $62.5\ \mu \text{s}$ to improve the signal from the sea of statistical fluctuations.

Zoom In Zoom Out Reset image size

Thus, our model exhibits some of the fluctuations and imperfections observed in experiments. These figures could be refined by further adjustment of the system parameters to a specific experiment, or by introducing a second kind of K channel with a door, but the excitability properties of the model are clearly manifest.

It is interesting to relate the assumptions of this model to the main elements of the classical Hodgkin-Huxley theory. In this theory the crossing of ions through the membrane is described by time-dependent currents, representing the total of a large number of channels. These currents produce changes in the membrane potential according to the capacitor equation. Our first assumption, eq. (1), is exactly this, but applied to discrete charges (ions) instead of to currents. Moreover, according to HH, the charge intensity crossing many Na channels depends on the membrane potential according to a generic law,

$$\begin{equation} I_{\text{Na}} = g_{\text{Na}}(t) ( \Delta V - V_{\text{Na}}), \end{equation} \tag{ 12 }$$

where $V_{\text{Na}}$ is the Na Nernst potential and $g_{\text{Na}}(t)$ is the ionic conductance. Membrane conductances represent thus the average of the states of a large number of channels, each of them either open or closed. We have substituted this Ohm-type law by the Langevin dynamics of ions along a single channel. The Nernst potential is not a parameter of our model, but instead it is reached automatically (in a single-ion species situation) since it corresponds to the equilibrium state of our model. Analogously the Goldman-Hodgkin-Katz law is verified in the steady state corresponding to more general situations.

Moreover in eq. (12) the membrane conductance depends on other variables subjected to dynamical equations [1]. Namely this conductance depends on the so-called activation and deactivation functions m and h,

$$\begin{equation} g_{\text{Na}}(t)={\bar g_{\text{Na}}} m^3(t)\,h(t), \end{equation} \tag{ 13 }$$

where $\bar g_{\text{Na}}$ is a constant. Activation and deactivation functions are interpreted in the context of our model as the average state of each of the two channel doors, i.e. of our variables Y₁ and Y₂, for a large number of channels. The way these functions m and h are built is the kernel of the HH theory. They obey deterministic linear differential equations, chosen in such a way that each variables m, h have a single steady state that, depending on the value of $\Delta V$, ranges continuously from 1 (all doors open) to 0 (all doors closed). On the contrary our variables Y_i, representing the doors of a single channel, present two steady state (open and closed) in such a way that a stochastic dynamics permits transitions between both states.

The HH theory includes K channels with a different conductance (with a single door) and a ionic leakage. In our model we have only implemented the leak. The K channel with door could be implemented in our model straightforwardly, but it has not been necessary for obtaining excitability.

As a result, both in the HH theory and in our model, the coupling between potential and the channels state triggers a well-synchronized temporal sequence of events, resulting in a sudden discharge of Na ions, the appearance of the spike and the K flux restoring the potential.

We have presented a microscopic physical approach to the excitable properties seen in neuronal cell membranes. The main points of our approach are the following: ions obey classical statistical equations of motion; channels are pores with doors whose dynamics are controlled by elastic nonlinear potentials; and the electrostatic potential of the cell membrane follows the capacitor equation. Moreover, since it is constructed incorporating fluctuations according to fundamental statistical physics (namely according to the fluctuation-dissipation theorem), it provides the correct statistical fluctuations of the diverse variables. This model can then be used to study the dynamics of a small number of channels, and in particular it appears as specially suitable for analyzing single-channel experiments. Note that in global measurements of real neural spikes a large number of channels are involved, and fluctuations will be smoothed out.

We have shown that a single Na channel in the presence of K leakage constitutes an excitable system producing the characteristic spikes in the action potential. Our objective here was not to reproduce the exact form of the action potential for some specific channels or neurons but rather to formulate in terms of fundamental statistical mechanic laws the underlined physical mechanisms in this biomolecular process.

It is worth commenting about the model parameters and their specific values. All of them have a clear physical meaning. Ionic concentrations per length are fixed by the experimental densities and the estimated channel areas. Friction parameters are estimated from experimental time scales, and barrier heights are of the order of a few $k_{\text{B}}T$ as it is expected in the biomolecular scale. Parameters of the doors and the function in eq. (6) have been chosen to fix the door's steady states (open and closed) and their location inside the channel. Other parameters such Q_j and $\phi_{\text{ref}}$ are adjusted to enter in the experimental scale. Since their physical meaning is clear and they are used in physical dynamical equations, the whole model lies within the framework of well-founded physics.

This approach presents several perspectives worth being explored:

• –  
  All model elements are described by standard physical equations based on a single energy functional, and accordingly it is possible to address the energetics of an excitable event. Before and after a pulse the system is in the same thermodynamic state but several (few) charges have changed of reservoir: $\Delta q_{\text{Na}}$ influx of Na and $\Delta q_{\text{K}}$ outflow of K. Thus, it is easy to estimate their loss of Gibbs energy,

  $$\begin{equation} \Delta G = \Delta q_{\text{Na}} \,g(\text{Na}) + \Delta q_{\text{K}} \,g(\text{K}), \end{equation} \tag{ 14 }$$

  where $g(\text{Na}), g(\text{K})$ are the Gibbs energy per particle of Na and K ions.
• –  
  The approach allows for other channels and doors modelizations which could be related to different biochemical structures of the channel proteins. Each door would have specific effective parameters that can be estimated from appropriate experiments.
• –  
  The role of the ionic concentrations on the channel states has not been receiving so far enough experimental attention, but we have observed, in our simulations, important sensitivity due to the ion-door collisions (in this regard see, for instance, fig. 5 in ref. [15]).
• –  
  Our approach allows a new view, from statistical physics, of the well-established Hodgkin-Huxley theory and other models based on it.

Finally, it is worth remarking that we have employed the minimum set of elements that results in the excitable dynamics observed in biological membranes. In this regard, it could also be seen as a modelization of hypothetical primitive channels, which presumably would be much simpler than present biological structures, which are the result of a long evolution and likely much more sophisticated. Thus, our approach opens a complementary scenarioto study ionic channel phenomenology from fundamental physics.

This work was supported by the Spanish DGICYT Project No. FIS2012-37655 and by the Generalitat de Catalunya Projects 2009SGR14 and 2009SGR921. We acknowledge fruitful discussions with Profs. J. García-Ojalvo, F. Giraldez and R. Vicente from Universitat Pompeu Fabra.