Channel noise induced stochastic facilitation in an auditory brainstem neuron model

The form of all models we consider is a generalisation of the Hodgkin-Huxley model (see, for example, Hodgkin and Huxley (1952)). It is described by a differential equation of the following form

where ${C_{\rm m}}$ is the membrane capacitance (pF), $V(t)$ is the membrane potential (mV), ${\bar{g}_{\rm a}}$ is the maximum total channel conductance for ion channels of type a, ${E_{\rm a}}$ (mV) is the reversal potential for ion channels of type a, and $I(t)$ represents a constant injected current. Each term in the sum represents a current, $I_{\rm a}(t)=\bar{g}_{\rm a}G_{\rm a}(t)(V(t)-E_{\rm a})$, since the function $G_{\rm a}(t)$ is dimensionless and constrained to the interval $[0,1]$; this function represents the time course of a conductance. For a passive leak current (a=lk) we have $G_{\rm lk}(t)=1$ in this formulation.

In deterministic models, each conductance variable, $G_{\rm a}(t)$ is formed from the product of at least one ‘activation variable’ (or ‘subunit’) and each such variable, $x$, is described by a differential equation representing ion channel kinetics, of the form

where

and

In these equations, ${\alpha_{x}}(V)$ and ${\beta_{x}}(V)$ are voltage dependent rates for the opening and closing of ion channel subunits respectively, as in the standard Hodgkin-Huxley formalism.

For example, the standard Hodgkin-Huxley model Hodgkin and Huxley (1952) has two ionic currents, Sodium (a=Na) and Potassium (a=K), and three activation variables, $x=m$, $x=n$ and $x=h$, such that $G_{\rm Na}(t)=m^{3}h$ and $G_{\rm K}(t)=n^{4}$. The equations for $\alpha_{\rm Na}$, $\beta_{\rm Na}$, $\alpha_{\rm K}$ and $\beta_{\rm K}$ are well-known Hodgkin and Huxley (1952).

Later we will replace the standard deterministic form of each $G_{\rm a}(t)$ with stochastic processes that model the fluctuations induced by random opening and closing of a finite number of ion channels.

The model of ventral cochlear nucleus (VCN) neurons developed by Rothman and Manis Rothman and Manis. (2003a); Rothman and Manis. (2003b, c) is of the form of Eqn. (1). Unlike the standard Hodgkin-Huxley equation, it is instead comprised of three distinct potassium currents, a sodium current and a hyperpolarization-activated cation current, all of which have been found to be present in VCN neurons. These currents were isolated via the voltage clamp technique and kinetic analysis Rothman and Manis. (2003a). This allowed the kinetic schemes of the ion channels responsible for activation and inactivation of the given currents to be experimentally determined for the potassium currents and verified for the others Rothman and Manis. (2003a). Such a model scheme provides an ideal framework in which to study biophysically realistic channel noise beyond the standard Hodgkin-Huxley model and with potential applications to studying auditory neurons.

The model presented in Rothman and Manis. (2003c) is given in five different configurations corresponding to different cell types. Of these, the Type II configuration shows a classic phasic response upon perturbation, i.e. it produces action potentials only at onset of a stimulus change Izhikevich (2004). This configuration has been used to model bushy cells in the VCN Rothman et al. (1993); Manis and Marx. (1991), and favorably compared to experimental recordings of phasic MSO neurons in gerbils Gai et al. (2009). It was this configuration of the Rothman-Manis neuron model that was shown to exhibit SBSR by Gai et al. Gai et al. (2010). A phasic neuron is one that produces a single spike at the onset of a rapidly changing input, but typically does not fire at all for inputs that vary slowly relative to the action-potential generating time-constants of the model. Gai et al. also noted that the model fires in response to a certain range of slopes of the input rather than a range of amplitudes Gai et al. (2009, 2010), meaning that this neuron has a threshold defined by the magnitude of the derivative of the input current.

The particular Rothman-Manis model we use consists of the following ionic currents: ‘fast sodium’ (${I_{\rm Na}}$), ‘high threshold potassium’ (${I_{\rm KHT}}$), ‘low threshold potassium’ (${I_{\rm KLT}}$), ‘hyperpolarization-activated cation current’ (${I_{\rm h}}$) and leak current (${I_{\rm lk}}$). Note that in Rothman and Manis. (2003c), a third potassium current—‘fast transient’ K+ $(I_{\rm A})$—is discussed, but this is not required in the model configurations we use here.

The resulting model requires 6 activation variables, $m,h,n,p,w,z$ and $r$. The particular products that form the conductance time courses, and corresponding constants are shown in Table 1, as are the parameter values of the peak conductances, reversal potentials, leak current and capacitance that we use in the Results section below. The Type I-c, Type I-II and Type II configurations of the model as in Rothman and Manis. (2003c) are used in this study. These have been adjusted to account for temperature as in Gai et al. (2009). The functions describing the kinetics of each activation variable are given in Appendix A.

In the deterministic model given by Equation (1), each conductance time course, $G_{\rm a}(t)$, corresponding to ion channel type a, physically represents the mean fraction of open channels of that type. Mathematically, this is equivalent to the probability that a single ion channel is open (and hence in the conducting state).

In turn, the activation variables are representative of the fraction of ion channel subunits that are in the conducting (open) state. For example, the combined quantities of $m^{3}h$ and ${w^{4}z}$ (see Table 1) are representative of the Na⁺ and K⁺ channels having individually gated subunits for activation and inactivation, with there being three of type $m$ and four of type $w$. Subunit state transition diagrams can be found, for example, in Goldwyn et al. (2007); Sengupta et al. (2010).

In order to introduce a realistic model of channel noise, fluctuations in the fraction of open channels need to be modeled explicitly. As in Goldwyn and Shea-Brown (2011), we write the fraction of open channels as the sum of a mean value, and a stochastic term, i.e.

where $G_{\rm a,det}$ is the deterministic part of the conductance time course, and is equal to the mean fraction of open channels. For the Rothman-Manis model, $G_{\rm a,det}(t)$ is listed in Table 1 for each channel type. This mean is equivalent to the probability that a channel is open at time $t$, given activation variables $\alpha(t)$ and $\beta(t)$. The term $\xi_{\rm a}(t)$ represents channel noise; in an ensemble, the channel noise will be conditionally independent in all neurons, given the membrane potential in each neuron.

We now describe two existing stochastic models of channel noise. The first is the standard Markov model that accurately reflects conceptual models of ion channel subunit openings and closings. The second model is an approximation to the first model that enables greatly reduced simulation runtime.

To begin, we introduce $N_{\rm a}$ to represent the total number of channels of type a. The most conceptually simple way of incorporating channel noise into neuron models consists of tracking the state of each of the $N_{\rm a}$ channels of each type. The probability that a channel will change state is determined by the voltage dependent transition rates $\alpha_{\rm x}(V)$ and $\beta_{\rm x}(V)$ Dayan and Abbott (2005)—see Eqns (2)–(4). The transition between states in a channel are assumed to follow a Markov process and are therefore memoryless. In such a scheme, an open channel occurs when all subunits of a particular ion channel type are in the open state. For example, a sodium channel is in the open state when both the $m$ variable is in state $3$, and the $h$ variable is in state $1$. We introduce the random variable $\mathcal{N}_{\rm a,o}(t)$ to represent the number of open channels of type a at time $t$. We therefore have

Thus, from Eqn. (5), the fluctuations in this model can be expressed as

Direct computational simulation of this Markov chain model is an accurate method, but it is computationally expensive, since it requires numerous random numbers to be generated at each time step, as well as tracking the state of $N_{\rm a}$ channels for each ion channel type.

Recently, there has been considerable interest in developing new methods capable of approximating the full Markov chain description accurately, but which may be implemented in much faster simulations Goldwyn et al. (2011); Goldwyn and Shea-Brown (2011); Linaro et al. (2011); Schmandt and Galán (2012); Dangerfield et al. (2012); Orio and Soudry (2012); Huang et al. (2013). These methods generally represent channel noise using stochastic differential equations (SDEs) that statistically approximate the fluctuations present in the Markov chain description by adding derived noise terms to the deterministic system of equations. In a recent review Goldwyn and Shea-Brown (2011), it was concluded that a simulation method based on much earlier mathematical theory due to Fox and Lu (1994) is accurate and efficient. Therefore, here we use that method in the form described by Goldwyn and Shea-Brown (2011), and compare its performance when applied to the problem of identifying SBSR, versus simulation of the Markov model. However, for our ISR studies, although we found that for small numbers of ion channels that the Markov model is consistent with the SDE approximation, a very large number of channels is needed for the effect to be observable. Consequently, a complete comparison of the Markov and SDE approaches for ISR was computationally infeasible, and our results are confined to the SDE approach.

A full treatment of the SDE method can be found in Fox and Lu (1994); Goldwyn et al. (2011); Goldwyn and Shea-Brown (2011) and the references therein. Briefly, the method determines the deterministic component of the conductance time course in the same way as the original deterministic model, but also derives approximations to the noise term, $\xi_{\rm a}(t)$, based on a system size expansion (SSE) of the Markov chain (see Fox (1997); Fox and Lu (1994); Pakdaman et al. (2010)). The noise for each channel type is assumed to be a time-varying Gaussian process that is dependent on the associated activation variables, and therefore also on the membrane potential. Following the approach described in Goldwyn and Shea-Brown (2011) we obtained $\xi_{\rm a}(t)$ by solving a set of $M_{\rm a}$ coupled SDEs for ion channel type a with $M_{\rm a}$ states, of the following form

where ${\bf x}_{\rm a}$ is a length-$M_{\rm a}$ vector, ${\bf A}_{\rm a}$ and ${\bf D}_{\rm a}:={\bf S}_{\rm a}^{2}$ are $M_{\rm a}\times M_{\rm a}$ voltage-dependent drift and diffusion matrices respectively. Noise enters the equations via ${\bf W}_{\rm a}$, which is a vector of $M_{\rm a}$ independent standard Brownian processes. The desired outcome $\xi_{\rm a}(t)$ is given by the $M_{\rm a}$–th element of the vector ${\bf x}_{\rm a}$. Once this is found, the membrane potential equation is given by

Solving (8) requires first finding ${\bf D}_{\rm a}$, and then obtaining its matrix square-root. The drift and diffusion matrices may be constructed for an arbitrary kinetic scheme with $M_{\rm a}$ states, comprised from multiple activation variables (e.g. $m$ and $h$ for sodium channels in the Hodgkin-Huxley model) as follows. The $M_{\rm a}\times M_{\rm a}$ drift matrix ${\bf A}_{\rm a}$ is identical to the transition matrix from the master equation representation of the Markov chain, as noted by Goldwyn et al. (2011), i.e.

where the $M_{\rm a}\times 1$ vector ${\bf p}_{\rm a}$ corresponds to the state occupancy probabilities Dayan and Abbott (2005). These probabilities can be determined at time $t$ from the activation variables (e.g. $m$ and $h$ for a potassium channel), and thus the $M_{\rm a}$–th element of ${\bf p}_{\rm a}$ is equal to the probability that a channel of type a is open (conducting).

In Goldwyn and Shea-Brown (2011), the diffusion matrix ${\bf D}_{\rm a}$ is written for the specific cases of sodium and potassium channels in the Hodgkin-Huxley model. However, here we provide a general method for constructing ${\bf D}_{\rm a}$ from any given drift matrix ${\bf A}_{\rm a}$, and the number of channels of type a, $N_{\rm a}$. The relationship can be expressed as

where ${\bf 1}_{M\times 1}$ is a $M_{\rm a}\times 1$ column vector with all entries equal to $1$, ${\bf I}_{M\times M}$ is the $M_{\rm a}\times M_{\rm a}$ identity matrix, and $\circ$ indicates the Hadamard operator (term by term multiplication). Note that all rows and columns of ${\bf D}_{\rm a}$ sum to zero, and ${\bf D}_{\rm a}^{\top}={\bf D}_{\rm a}$.

Finally, although Eqn (11) shows the theoretical dependence of ${\bf D}_{\rm a}$ on ${\bf x}_{\rm a}$, in practical simulations it is useful to replace the stochastic vector ${\bf x}_{\rm a}$ with ${\bf p}_{\rm a}$, to ensure matrix square roots of ${\bf D}_{\rm a}$ can be found, as discussed in Goldwyn and Shea-Brown (2011).

Very small amplitude white noise has been shown to cause stochastic switching between a periodic limit cycle (indicating action potential generation) and a resting state in Hodgkin-Huxley model neurons Tuckwell and Jost (2012); Gutkin et al. (2009); Tuckwell et al. (2009); Tuckwell and Jost (2010). This occurs for constant suprathreshold input currents that are close to the firing threshold for tonic firing in the deterministic version of the model. Moreover it was shown that there exists a noise amplitude that will minimize the average firing rate over many trials for a given input current. The occurrence of this minima, which is below the firing rate for the corresponding deterministic model, has led to the effect being labeled ISR Tuckwell and Jost (2012).

We have demonstrated, using both sinusoidal and ramp inputs, that slope-based stochastic resonance can be observed due to channel noise. In the case of sinusoidal inputs, we found that channel noise induced spiking in response to low frequency inputs occurs only during the rising phase. However Gai et al. (2009) also reported firing during the falling phase of the sinusoidal input, when using injected current noise. This confirms the comment of Gai et al. (2009), that unrealistically high noise levels were used to generate action potentials on the falling phase of the input.

In the case of a ramp input current, SBSR was demonstrated by production of a sigmoidal increase of firing rate with increasing slope. It has been shown previously that such a sigmoidal shape in the probability of firing may be exploited for increased information transmission in a parallel population of stochastic Hodgkin-Huxley models Ashida and Kubo (2010). This previous demonstration by Ashida and Kubo (2010) of a form of stochastic facilitation, known as suprathreshold stochastic resonance Stocks (2000); McDonnell et al. (2007), due to channel noise is consistent with the central conclusions of the current paper.

We have demonstrated that ISR can occur due to channel noise in the Hodgkin-Huxley and Rothman-Manis Type I-II models. As for the case of injected current noise Tuckwell and Jost (2012); Gutkin et al. (2009); Tuckwell et al. (2009); Tuckwell and Jost (2010), the effect is observed due to stochastic switching between a limit cycle and a fixed point. The minima in the ISR curve occurs when the noise level is such that the probability of initial switch from limit-cycle to fixed point, followed by a long dwell in the fixed point, is highest.

Both the Hodgkin-Huxley and Rothman-Manis models exhibit Class-II excitability Izhikevich (2007), which means that the deterministic forms of the models exhibit high tonic firing rates above some threshold value of constant injected current, $I$, and no spiking below the threshold, such that spike rates discontinuously jump from zero to a large value as $I$ increases. This behavior can be attributed to the existence of a subcritical Andronov-Hopf bifurcation Izhikevich (2007). As noted by Tuckwell et al. (2009), the combination of this dynamics with noise means the model can stochastically switch between the resulting stable limit-cycle and fixed point, thus leading to bursting behavior that causes ISR. Such noise-induced bursting is consistent with the results of a previous detailed study of channel-noise induced bursting in the Hodgkin-Huxley model Goldwyn et al. (2007).

On the other hand, we did not observe ISR in the Rothman-Manis Type I-C model. The Type I-C model has a much lower threshold to tonic spiking in response to constant current, I, in comparison with the Type I-II model, and as mentioned lacks a low threshold potassium current. This is sufficient to change the dynamics of the system such that the deterministic form of the model exhibits Class-I excitability, which means that as the input current $I$ is brought above threshold, spike-rates increase continuously from zero. Such behavior is attributable to the existence of a saddle-node-on-invariant-circle bifurcation, which means a stable fixed point disappears for above-threshold current levels Izhikevich (2007). Thus, the system acts like an integrator (as noted by Rothman and Manis. (2003a)) and there is no scope for bistable switching from a limit-cycle to a rest-state, and hence no opportunity for ISR to be observed.

The extent of the stochastic behavior observed in all models due to channel noise may be qualitatively understood by considering the case where the model neuron’s membrane potential is close to causing action potentials to occur. In this regime, fluctuations in all membrane currents can contribute to whether or not an action potential occurs. The size of the fluctuations in each current are determined by variance in the number of open channels, $\mathcal{N}_{\rm a,o}(t)$, the associated peak conductance, and the difference between the membrane potential and the reversal potential.

Since the variance of an activation variable is high when $\mathcal{N}_{\rm a,o}(t)$ is low Schneidman et al. (1998); Goldwyn and Shea-Brown (2011), if such a current has a large peak conductance, the path variance of the membrane potential may also be significantly affected in this case.

For the models with Class-II excitability in which ISR occurs, the sodium channel has the largest peak conductance and a value of $G_{\rm Na}(t)$ that is close to zero when the membrane potential is close to threshold. Thus, the occurrence of large $P_{R\rightarrow L}$ and $P_{L\rightarrow R}$ that results in rapid state switching for low channel numbers can be attributed primarily to noise in the sodium channels. Indeed, we found that removing the stochastic noise from only the sodium channels meant that much lower channel numbers in remaining channel types was required to see state switching.

The small noise amplitude required to increase firing rates above the minimum as noise decreases in ISR corresponds to a large number of channels, showing that stochastic effects due to ion channel noise are present near threshold for large $N_{\rm a}$. This is contradictory to the frequently used assumption that for large $N_{\rm a}$ the fluctuations due to channel noise may be neglected, but is consistent with other previous studies of the influence of channel noise Schneidman et al. (1998); Shuai and Jung (2005); Faisal and Laughlin (2007); Goldwyn and Shea-Brown (2011).

Our simulations took advantage of the fact that the run time of the SSE noise model does not scale with $N_{\rm a}$, as the Markov model does. Therefore large numbers of ion channels were able to be efficiently simulated, allowing this effect to be elucidated.

As mentioned above, our results assume that all channel types have the same number of channels. This is a simplification made for the purposes of assessing whether channel noise can impact on spiking in a similar manner to current noise. Having established that it can, it is of interest to consider the influence of numbers of channels that vary realistically for each channel type.

Various methods can be used to estimate total numbers of channels, e.g. see Buchholtz et al. (2002). Here, however, we make use of existing data on the conductance of individual channels from  Gutman and et. al (2005); Hofmann et al. (2005); Catterall et al. (2005) and then divide the total conductance for each channel type in our models by this value to obtain estimates. Since total conductance is temperature dependent, and the temperature at which the individual channel conductances were measured is not available, we use the total conductances used in Rothman and Manis. (2003c). This procedure is facilitated by discussion in Rothman and Manis. (2003a) of the molecular identities of each channel type in the Rothman-Manis model. Our results are tabulated in Table 2, including references to papers used to obtain the single channel conductances.

The results in Table 2 are consistent with the following analysis. Using the neuronal diameter for the VCN neuron 21 $\mu m^{2}$ Rothman and Manis. (2003c) and assuming channel densities are about 18 per $\mu{\rm m}^{2}$, for potassium channel and about 60 per $\mu{\rm m}^{2}$ for sodium channels as in Hodgkin-Huxley type models Schneidman et al. (1998), the number of channels in a spherical soma is estimated at being about 25000 potassium channels and about 83000 sodium channels.

Our first order estimations place the estimated number of channels in vivo in a region where stochastic effects are strongly observed in these simulations.

When channel noise is accounted for in model neurons, the response for inputs close to a firing threshold (whether an amplitude or slope threshold) may nontrivially differ from the deterministic model. Here we have shown that with the inclusion of intrinsic channel noise in a model of a VCN neuron that stochastic effects significantly alter the close to threshold behavior for both a phasic firing and a tonic firing configuration of the model, showing the effects of SBSR and ISR respectively. For the effect of SBSR both the Markov chain model and the SSE model were used to show the effect. Thus our study suggests that it would be interesting in future work to carry out detailed analysis of whether coding and processing of acoustic stimuli in the auditory brainstem may benefit from channel noise.

The assumption that for a large number of ion channels the fluctuations due to individual ion channels may be disregarded was found to be incorrect for some cases. We emphasise that the path of the membrane potential in phase space may have a large variance due to intrinsic channel noise for neurons in which the current making the strongest contributions to the membrane potential at threshold has low occupancy in the conducting (open) state. In agreement with this, stochastic effects were found to strongly influence behavior of the tonically firing models considered here at channel numbers in the range of those found in biological neurons. The different forms of stochastic behavior displayed by these tonic models were able to be qualitatively explained using considerations of the near-threshold value of $\mathcal{N}_{\rm a,o}(t)$, the associated membrane conductances and the near-threshold behavior of the model. These considerations when applied to the phasically firing Type II model predict that for a stochastic effect to be seen, the number of channels must be quite low to cause variance in dominating low threshold current, as indeed is the case.

In order to perform simulations for large numbers of individual ion channels the SSE model Fox and Lu (1994); Goldwyn et al. (2011) was used. Future work may include comparing the results of these simulations with other stochastic differential equation approximations to the Markov chain description, such as those found in Schmandt and Galán (2012); Dangerfield et al. (2012); Orio and Soudry (2012); Huang et al. (2013). It will also be useful to consider stochastic facilitation effects due to channel noise in more detailed multi compartmental models.