The What and Where of Adding Channel Noise to the Hodgkin-Huxley Equations

Conductance-based equations for electrically active cells form one of the most widely studied mathematical frameworks in computational biology. This framework, as expressed through a set of differential equations by Hodgkin and Huxley, synthesizes the impact of ionic currents on a cell's voltage—and the highly nonlinear impact of that voltage back on the currents themselves—into the rapid push and pull of the action potential. Later studies confirmed that these cellular dynamics are orchestrated by individual ion channels, whose conformational changes regulate the conductance of each ionic current. Thus, kinetic equations familiar from physical chemistry are the natural setting for describing conductances; for small-to-moderate numbers of channels, these will predict fluctuations in conductances and stochasticity in the resulting action potentials. At first glance, the kinetic equations provide a far more complex (and higher-dimensional) description than the original Hodgkin-Huxley equations or their counterparts. This has prompted more than a decade of efforts to capture channel fluctuations with noise terms added to the equations of Hodgkin-Huxley type. Many of these approaches, while intuitively appealing, produce quantitative errors when compared to kinetic equations; others, as only very recently demonstrated, are both accurate and relatively simple. We review what works, what doesn't, and why, seeking to build a bridge to well-established results for the deterministic equations of Hodgkin-Huxley type as well as to more modern models of ion channel dynamics. As such, we hope that this review will speed emerging studies of how channel noise modulates electrophysiological dynamics and function. We supply user-friendly MATLAB simulation code of these stochastic versions of the Hodgkin-Huxley equations on the ModelDB website (accession number 138950) and http://www.amath.washington.edu/~etsb/tutorials.html.


Introduction
Understanding the role of noise in cellular dynamics and function is a central challenge across computational biology.This is as true in neuroscience as in any field [1][2][3], and a universal source of noise in electrically active cells that has garnered increasing attention is the stochastic activity in ion channels [4][5][6].This channel noise has been studied in a variety of neural systems including electrical stimulation of the auditory nerve by cochlear implants [7,8, e.g.], as well as in entorhinal cortex [9], cerebellar granule cells [10], and hippocampal CA1 pyramidal neurons [11].Modeling studies have suggested that channel noise can influence information processing [12], spike time reliability [13], stochastic resonance [14], firing irregularity [10,15], subthreshold dynamics [9,10], and action potential initiation and propagation in morphologically detailed models [11,16].Channel noise is at work in many other systems such as the activity of cold receptor cells [17], nicotinic acetylcholine receptors [18], and calcium release by Inositol 1,4,5-Triphosphate receptors [19].
Despite widespread interest in channel noise, it has remained unclear what the options are for including this noise source in the canonical model of neurophysiology -the Hodgkin-Huxley (HH) equations for the 2 action potential [20].The direct approach provides a gold standard: each of N channels of a particular type transitions independently and randomly among discrete configurational states.This yields a continuous-time Markov chain with voltage-dependent transition probabilities; see [21] for a recent review.In the limit that N → ∞ for each channel type, the classical Hodgkin-Huxley equations are recovered [22][23][24][25][26][27].For finite N , one simulates the Markov process via a Gillespie-type algorithm [16,[28][29][30].
Is there a simpler approach, where one modifies the HH equations by adding a few well-placed noise terms?Beyond conceptual and computational simplicity, this would provide a direct link to powerful results on the dynamics and geometry of these equations [31,32].This line of research was initiated by Fox and Lu [22,33] who derived candidate sets of stochastic differential equations (SDEs) using a system-size expansion applied to the Markov chain.The past few years have seen increasing interest in this problem, spurred on by the promise yet apparent shortcomings of this SDE approach [10,27,[34][35][36][37][38][39].
As recent work attests [27,38,39], accurate methods for incorporating channel noise into the HH equations are finally emerging in the form of methods both new and old.These works demonstrate that adding noise terms to the HH equations can indeed give a compressed and accurate reproduction of the channel fluctuations.However, the placement of these terms is critical, and -as a decade of research attests -less than obvious.A key focus of our review is a unified presentation of the methods that provide the most accurate approximations to Markov chain models of channel noise.A common feature of these methods is that they introduce noise processes as conductances in the HH equations.

Stochastic versions of the Hodgkin-Huxley Equations
As they are an almost universal point of reference for neuron modeling, we consider the equations introduced by Hodgkin and Huxley to model action potentials in the squid giant axon [20].
Here, V is the membrane voltage, the gating variables m, h, and n represent the fraction of open channel subunits of different types, aggregated across the entire cell membrane.These fractions are combined in the terms m 3 h and n 4 to regulate total conductances for Na + and K + currents.The constant C represents the capacitance of the cell membrane; E Na , E K , and E L are reversal potentials; ḡNa and ḡK are maximal conductances; and g L is the leak conductance.
Comprehensive introductions to this model can be found in many standard texts [23,31,32].We emphasize that our discussion applies to any conductance-based model of excitable cells, including point, compartmental, or spatially extended neurons, and related models of calcium release [40].To model channel noise within this differential equation framework, we seek ways of introducing fluctuations into this (classically) deterministic system.In this review, we consider three approaches, which we classify as follows: • Current noise: Replace Eq. 1 with where ξ V (t) is a Gaussian white noise process.
• Subunit noise: Replace Eq. 2 with where the ξ x (t) are Gaussian processes that may depend on x and V .• Conductance noise Replace Eq. 1 with where the noise processes ξ Na (t) and ξ K (t) are Gaussian processes that may depend on x and V .
Table 1 summaries the differences among these models, which we now discuss in detail.
A summary of the three classes of channel noise models that we discuss in this review, and how they differ from the deterministic HH equations which have no noise.

Current noise
The simplest method for incorporating noise into the HH equations is to add a fluctuating current term ξ V (t) to the dV /dt equation, as shown in Eq. 1*.Here, we assume ξ V (t) is only a function of time.Stochastic currents of this form are frequently used to drive the HH model, often in the context of the diffusion approximation for synaptic inputs [41][42][43].In the present context, however, we emphasize that ξ V (t) is meant to represent the combined effect of the stochastic activity of ion channels on the voltage dynamics of the cell.This approach is appealing due to its simplicity, but since channel noise is generated by the stochastic activity of ion channels in the cell membrane, it seems likely that the fluctuation term ξ V (t) should also depend on V or the subunit variables.Another drawback is that, to date, there is no principled method for determining the intensity of the noise.Nonetheless, there may be cases in which current noise can be justified on empirical grounds.For instance, for a single membrane area and a constant applied current, Rowat compared the interspike interval distribution generated by a Markov chain model to the distribution generated by HH equations with current noise and found remarkably close agreement [15].

Subunit noise
In the HH framework, an ion channel's configuration is determined by the states of its constituent subunits, where each subunit can be either in an open or closed state [6,23,44].For instance, each K + channel is composed of four n-type subunits, all of which must be open in order for the channel to be permeable to K + ions.Each subunit randomly transitions between its open and closed state.This suggests that the most appropriate place to add noise may be to the equations that describe the fractions of open subunits, as in Eq. 2*.Moreover, since all subunits are independent and all subunits of the same type are statistically identical, it is tempting to combine the resulting noisy fractions of open subunits to regulate conductances in the same way as one would in the deterministic HH equations; namely by computing m 3 h and n 4 .The variables m, h, and n represent the aggregated fraction of open subunits, whereas the quantity that influences the membrane potential is the fraction of individual open channels.In the limit of infinitely many channels (and therefore vanishing fluctuation terms), m 3 h and n 4 do correctly model the fraction of open channels.For a finite number of channels, however, there is no guarantee that fluctuations in the these quantities will correctly model fluctuations in the membrane-wide fractions of open channels.
To see this, note that if all channels were gated by a single subunit, then the subunit model would be appropriate -in this case, the (noisy) fraction of open subunits is identical to the (noisy) fraction of open channels.In the HH model, however, four subunits gate each channel.Combining the quantities m, h, and n together to form the quantities m 3 h and n 4 neglects the important fact that each ion channel is composed of a specific package of subunits.The states of the particular subunits within a channel, not the average state of all subunits in the cell membrane, determine whether that channel is open or closed.Thus, random transitions of individual channels among their different configurational states occur with different statistics than predicted by random transitions of the aggregated subunit variables alone [27].This fact leads to quantitative errors produced by the subunit noise approach, as we will review below.
Subunit noise was first proposed in [22] and has been used many times; see [10,14,17,19,[45][46][47][48][49][50][51][52], among others.By applying a system size expansion to the states of populations of subunits, Fox and Lu arrived at a Langevin equation description of the subunit dynamics, precisely of the form of Eq. 2*, where the noise terms ξ x (V, t) (x = m, h, or n) are Gaussian processes with covariance function Here, δ(•) is the Dirac delta function and N represents either the number of Na + channels for the m and h subunits or the number of K + channels for the n subunit.Although the authors acknowledged that the subunit noise approach has no rigorous justification and must be validated empirically, it has been widely used as an approximation to Markov chain ion channel models.However, numerical studies have revealed inaccuracies in this approximation that persist even as the number of channels increases [35,37].Relative to the Markov chain model, the subunit noise models produce weaker conductance and voltage fluctuations [37,53], lower firing rates [12] (and, equivalently, longer mean interspike intervals [35]), less variability in the occurrences and timing of spikes in response to a brief pulse of current [34,36], and they transmit information at a higher rate [12].Furthermore, mathematical analyses of the voltage clamp statistics of the subunit noise model have proven that it does not generate stationary distributions of open channels that accurately approximate those of the Markov chain model [27,38].The analysis in [27] revealed similar inaccuracies in a related model proposed by [19], in which the terms m 3 h and n 4 terms in Eq. 1 are replaced by m 1 m 2 m 3 h and n 1 n 2 n 3 n 4 , respectively, where the subscript denotes independent solutions to SDEs of the form of Eq 2*.Others have proposed simplifying Eq. 3 so that the noise terms do not depend on V , and are simply Gaussian white noise [10].While such approaches may be justified on empirical grounds, in general they should not be considered as systematic approximations to Markov chain ion channel models.

Conductance noise
The remaining possibility is to incorporate fluctuations directly into the fractions of open channels.This seems natural, as the fraction of open channels controls ionic currents.Our intuitive understanding of the HH equations, which can be made rigorous as in [23,25,27], tells us that the mean fractions of open Na + and K + channels are given by m 3 h and n 4 .The most direct approach to adding channel noise to the HH equations, therefore, is to add zero mean stochastic processes to the deterministic values of m 3 h and n 4 .Following this idea leads to Eq. 1**, which is a compact mathematical description of channel noise that preserves the original structure of the HH equations and has the biophysically desirable interpretation that channel noise induces fluctuations in the ionic conductances.We now review three channel noise models [22,27,38] and, with a brief set of calculations, place them in the unified framework of conductance noise.

Conductance noise models based on voltage clamp
Two recent studies have developed conductance noise models based on stationary statistics of channel activity in voltage clamp -called the "quasistationary" channel model in [27] and the "effective" model in [38].Using the standard assumption that all ion channels are independent, the stationary distribution of open channels in voltage clamp is a binomial distribution parameterized by the total number of channels and the probability that any given channel is open.The probability that a channel is open depends on V , and thus a voltage clamp analysis generates a family of binomial distributions indexed by V , which is treated as a fixed parameter.The means of the distributions of open channels are given by familiar terms from the deterministic HH equations: m 3 h for Na + channels and n 4 for K + channels.If these binomial distributions are well approximated by Gaussian distributions, then the stationary distribution of open channels in voltage clamp can be accurately approximated by a family of zero mean, voltage dependent Gaussian processes that are added to the voltage dependent equilibrium values of m 3 h and n 4 .
The effective model of [38], for instance, represents the fraction of open K + channels in voltage clamp as n 4 + ξ K (V, t) where the stochastic processes are sums of independent Ornstein-Uhlenbeck (OU) processes (i.e.Gaussian colored noise).In other words, ξ K (V, t) = i ζ i (V, t), where the ζ i (t) are defined by SDEs of the form: with timescales τ i (V ) and noise amplitudes σ i (V ) [27,38].The quasistationary channel model in [27] produces equivalent Gaussian processes in voltage clamp.The difference between the two methods is that, in [27], there is a single noise process shared by all OU processes: dW i (t) = dW (t) for all i in Eq. 4. While this leads to different values of σ i (V ), our own simulations of these models (not shown) did not reveal any systematic differences in the outputs of the two models.
To simulate these models with a freely evolving membrane potential, one must assume Eq. 4 is valid outside of voltage clamp.In practice, one numerically integrates Eq. 4, where V is updated in each time step according to Eq. 1**.There is no assurance that this approach is valid in the context of a dynamic membrane potential.If V changes on longer time scales than the correlation times in the conductance fluctuations, then such an approximation may be appropriate, but a defining feature of neural dynamics is the rapid change in V during the course of an action potential.Voltage clamp-based methods may be less reliable, therefore, for modeling the spiking activity of neurons.
These channel noise models were developed in [27] and [38] in order to approximate the original Markov chain description of channel kinetics.Their structural details -i.e. the number of ζ i (V, t) processes used to define ξ K (V, t) and ξ N a (V, t) and the values of τ i (V ) and σ i (V ) in Eq. 4 -were defined based on the stationary statistics of Markov chain model.The voltage clamp approach itself, however, can be made general and model-independent.The only necessary ingredients are the autocovariance functions, as a function of the voltage clamp value, for fluctuations in the conductances.Moreover, if these stationary autocovariance functions can be expressed as sums of exponential functions, then the Gaussian representation theory for multiple Markov processes ensures that they can be approximated as a linear combination of OU processes [54].

Conductance noise models based on Fox and Lu's system size expansion
Lacking in all of the previously discussed methods is a direct approach for modeling the dynamics of fluctuations in the fractions of open channels as the voltage V dynamically evolves.Surprisingly, the early work of Fox and Lu addressed this problem, but has apparently gone overlooked ever since.Fox and Lu derived a system of SDEs in which each dynamical variable represents the fraction of ion channels in a specified configuration.This differs from their more widely-used model, the subunit model discussed previously, in which the dynamical variables represent the fractions of open subunits.The resulting system of SDEs does not visibly resemble the HH equations, but with a few calculations we next show that this approach produces a conductance noise model in the form of Eq. 1**.
The starting point of Fox and Lu's analysis are vectors that describe the fractions of Na + and K + channels in each configuration as a function of time.We denote these by y(t) and x(t).For instance, the elements of x represent the fraction of K + channels that have all subunits closed, three subunits closed and one open, etc.The state that will be of most interest is the conducting state, in which all subunits are open.We denote the corresponding elements of y and x as y o and x o , and write the current balance equation as: The dynamics of y o and x o are determined by drift and diffusion matrices, which Fox and Lu obtained from the original Markov chain description through a system size expansion [22,33,55].We omit the details of the system size expansion, which can be found in [22,33].A rigorous discussion of a related method for passing from the Markov chain kinetics to a system of SDEs has been recently presented [26], and we do not reproduce those details here.The result of Fox and Lu's expansion is a coupled system of linear SDEs of the form: The matrices A Na (V ) and A K (V ) are the drift term or deterministic part of the dynamics, and are identical to the transition matrices from the master equation representation of the Markov chains for the Na + and K + channels [22,25,27].The matrices S Na (V, y) and S K (V, x) are matrix square roots of diffusion matrices; they depend on the state variable and the voltage-dependent transition rates.Stochasticity arises via the independent, standard Brownian processes W Na (t) and W K (t).
We now demystify the connection between these equations, in which fractions of open channels are obtained from a high-dimensional system of coupled SDEs, and the standard HH equations, in which the fractions of open channels depend on the subunit variables.The key is to split the equations for x and y into two parts: a deterministic equation that exactly matches the gating variable equation ( 2), and a fluctuation equation for the noise terms.To accomplish this we define new variables x and x, which evolve via: with initial conditions x(0) = x(0) and x(0) = 0.The sum x+x solves Eq. 7, so this is an exact decomposition of x into a deterministic part x and a fluctuation part x.We can also apply a similar decomposition to y.As discussed by a number of authors [23,25,27], solutions to the deterministic equation (Eq.8) can be generated by appropriate combinations of m, n, h, the gating variables from the deterministic HH equations: ȳo = m 3 h and xo = n 4 .This leaves the fundamental structure of the HH equations intact.Eq. 5 can be replaced by the modified HH voltage equation (Eq.1**), where the conductance noise terms ξ K (V, t) and ξ Na (V, t) are defined to be xo (t) and ŷo (t), respectively.In sum, the high-dimensional SDEs derived by Fox and Lu [22] do not modify the deterministic structure of the HH equations.Instead, as shown in Eq. 9, their sole purpose is to shape the fluctuations in the fractions of open channels.An important strength of this method is that it yields a description of channel fluctuations that is equally valid outside of voltage clamp.Furthermore, as shown in [27], the stationary statistics of open channels for this method match exactly those of the Markov chain model, and it accurately replicates spiking statistics for channel numbers as small as 600 Na + and 180 K + channels (membrane area of 10 µm 2 ).
One complication in solving these systems of SDEs is the need to determine S Na (V, y) and S K (V, x), by computing matrix square roots in each time step.In order to guarantee the existence of these matrix square roots, we replace the values y and x in the diffusion matrices with deterministic values obtained from the gating variables, or equivalently the solutions Eq. 8 for x and the corresponding equation for y.

Comparing stochastic versions of the Hodgkin-Huxley equations: simulations
How well do the simplified noise models match the "gold standard" Markov chain model of ion channel kinetics?Extensive comparisons between Markov chain and subunit noise models have been reported in prior studies [12,27,34,35,37].Studies have also compared Markov chain models to a current noise model [15], voltage clamp conductance noise models [27,38], and Fox and Lu's system size derived conductance model [27,39].An exhaustive numerical investigation of these approaches is beyond the scope of this review, but in Figs.1-3 we show simulation results that illustrate key differences among these approaches.All simulations use standard parameter values for the HH equations [20].The voltage clamp conductance noise model is defined as in [38].In all simulations, we used the Euler-Maruyama method with 0.01ms time step for solving the relevant differential equations [56] and a Gillespie-type algorithm to implement the ion channel kinetics in the Markov chain [28,30].To generate Gaussian pseudorandom numbers, we produced uniform pseudorandom numbers with the Mersenne Twister algorithm [57] and then transformed these using the Box-Muller method [58].Simulation code is available upon request, and is based on the work of [38] and [27].Both of these groups have made their code available on the ModelDB website [59], accession numbers 127992 and 128502, respectively.To complement this review, we supply user-friendly Matlab simulation code of these stochastic versions of the Hodgkin-Huxley equations on the ModelDB website (accession number 138950) and at our website http://www.amath.washington.edu/∼etsb/tutorials.html.
We will first compare time-varying distributions of the fractions of open channels.Intuitively, one would expect that the number of open channels (all of which assumed to be independent), should be binomially distributed.For a predefined voltage trajectory, this is indeed the case, as has been proven by [25].The time-varying distributions of the fractions of open Na + and K + channels in a Markov chain model of ion channel kinetics approach an asymptotically stable, voltage dependent binomial distribution with means and variances given by solutions to the deterministic subunit equations Eq. 2: Var[Fraction Open Na Var[Fraction Open Na We can use this result to compare channel noise models outside of voltage clamp.Fig. 1A shows a single voltage trace obtained from a Markov chain model with 6000 Na + channels and 1800 K + channels (membrane area 100µm 2 ) with no applied current (I = 0).Using this sample path as an input to the channel noise models, we compare the statistics of the fractions of open channels for the different models.Panel B shows the mean fractions of open Na + and K + channels, as computed from Eqs. 10 and 11.All channel noise models produced mean values that were in close agreement with these values, so we did not plot those results.
The results for the variance of the fractions of open channels, as shown in Fig. 1C and 1D, tell a different story.The variance in the fractions of open Na + channels are computed from Eqs. 12 and shown in black in Panel C. The variance is accurately captured by Fox and Lu's conductance noise model (red), but misestimated by the subunit noise model (blue) and voltage clamp conductance noise model (green).Of particular note is the fact that the voltage clamp conductance noise model utterly fails to track the Markov chain variance during the spike (right inset of Panel C).This illustrates the point, made earlier, that voltage clamp methods may not be appropriate in regimes when V changes rapidly.The subunit noise model underestimates the variance during the subthreshold period (left inset), and overestimates the variance during the spike at ∼ 70ms (right inset).
Fig. 1D shows variances in the fraction of K + channels.Again, Fox and Lu's conductance noise model is most consistent with the equilibrium binomial distribution result.The voltage clamp model provides a reasonably close approximation, but the subunit noise model alternately undervalues the variance prior to the spike (see inset), and overvalues the variance near the time of the spike.
We also simulated spike trains from the channel noise models, using a current input of the form where ξ I (t) is a Gaussian white noise process with zero mean E[ξ I (t)ξ I (t )] = δ(t − t ).In Fig. 2, we show the mean and coefficients of variance (CV) of interspike intervals (ISIs) obtained from simulations of the Markov chain and SDE models in response to a constant current input (I noise = 0µA/cm 2 ).Similar simulation results have been reported in [12,15,27,35].We present results for different amounts of constant current For Panels C and D, exact variances (black) were computed from Eq. 12 and Eq. 13 and all other variances were estimated from 5000 repeated simulations of the channel noise models.input (x-axis) and two different membrane areas.The magnitude of fluctuations in the current noise model was chosen so that the mean insterspike interval of this model would match that of the Markov chain model: ξ V (t) = 1.94η(t) for a membrane area of 10µm 2 and ξ V (t) = 4.3η(t) for a membrane area of 100µm 2 , where η(t) is a Gaussian white noise process with mean zero and E[η(t)η(t )] = δ(t − t ).
In the left column, we see that all models, with the known exception of the subunit noise model (blue), accurately reproduce the mean ISIs of the Markov chain (black), although there are slight discrepancies apparent for the current noise (cyan) and voltage clamp (green) methods.These discrepancies are even more visible when comparing the coefficient of variation of the ISIs (right column).Fig. 3 presents the results of similar simulations, the only difference being the addition of Gaussian white noise fluctuations to the applied input (I noise = 2µA/cm 2 ).The ISI statistics for the subunit noise model remain the most different from the Markov chain.For all conditions tested, it is clear that the Fox and Lu's conductance noise model (red) generates ISI statistics that are most similar to the Markov chain model.

Discussion
We stand at a promising moment for the study of channel noise in conductance-based models.In recent years, due to a spate of simulation studies drawing attention to discrepancies between subunit noise models and Markov chain ion channel models [12,27,34,35,37,38], there has been a growing sense of pessimism regarding whether SDEs could prove an effective framework for modeling the stochastic activity of populations of ion channels [53, e.g.].However, thanks to the development of novel approximation methods [27,38] and the rediscovery, analysis, and testing of past efforts [22,33], new life has been breathed into the SDE approach.The validity of SDE versions of the HH equations is now more clearly justified, and the door is open for these models to generate insight into how channel effects spike timing, reliability, propagation, and other aspects of neural dynamics.
A central theme of this review is that the addition of fluctuations in conductance terms, or equivalently in the fractions of open channels, should be the preferred way for including channel noise in the HH equations.This approach, which we have termed conductance noise, generates models that preserve the mathematical structure of the HH equations and that accurately approximate Markov chain models.In the case of the high-dimensional SDE model derived by Fox and Lu in [22], this was not obvious at first glance, and may be one reason why this aspect of their work has been overlooked.Through a brief calculation, however, we elucidated the connection between this model and the HH equations by showing how the high-dimensional SDEs can be decomposed into a deterministic part identical to the classical HH equations and a fluctuation part representing channel noise.
Although SDE models for channel noise are generally validated by making comparisons to Markov chain versions of the HH model, there is no guarantee that the Markov chain framework will remain the "gold standard."Indeed, critiques of the Markov chain approach have been articulated [60, cf.] and alternative mathematical models have been proposed [61, e.g.].With this in mind, it is useful to draw a distinction between "derived models" and "empirical models."The subunit and conductance noise models introduced by Fox and Lu [22,33] are in the former category.They are constructed with explicit reference to the conformational states of ion channels and their subunits, as defined by a Markov chain model of ion channel kinetics.In contrast, the current noise model and the voltage clamp conductance noise models can be thought of as "empirical" since they can be constructed from observable quantities.In our simulations, for instance, we used spontaneous firing rate to set the current noise level and the stationary statistics of open channels in the Markov channel model to define the noise processes in the voltage clamp conductance models.In principle, empirical measurements of conductance fluctuations in voltage clamp, without reference to a Markov chain model, could be used to construct conductance fluctuations.Empirical models that can be fit to, or validated against, quantities that are readily available from electrophysiological data are an attractive direction for future research, as they may inspire new methods for incorporating channel noise in conductance-based models.
The effects of channel noise have been a subject of intense interest in computational neuroscience and related fields in computational biology.The stochastic approaches reviewed in this paper represent an important extension of the conductance-based model framework introduced by Hodgkin and Huxley [20].Due to decades of analysis of the HH equations and an abundance of theoretical tools [62] and numerical methods [63, e.g.] for studying SDE models, we believe that appropriate methods for adding noise processes to the HH equations and their cousins throughout electrophysiology will play an important role in the future of computational biology.

Figure 1 .
Figure 1.Analysis of responses of channel noise models for a fixed voltage trajectory.A: Voltage trace obtained from the Markov chain model with no current input, 6000 Na + channels and 1800 K + channels.Dynamics are characterized by a prolonged subthreshold period followed by a spontaneous, channel-noise induced spike at 70ms.B: Means of fraction of open Na + and K + channels for the given voltage trace in Panel A, as computed from Eqs. 10 and 11.C: Variance in the fraction of open Na + channels.D: Variance in the fraction of open K + channels.Left insets in Panels C and D show magnified views of the period preceding the spike.Right inset in Panel C shows magnified view during the spike.For Panels C and D, exact variances (black) were computed from Eq. 12 and Eq. 13 and all other variances were estimated from 5000 repeated simulations of the channel noise models.

Figure 2 .
Figure 2. ISI statistics for DC input.A: Mean of ISIs for a membrane area of 10µm 2 (600 Na + and 180 K + channels).B: CV of ISIs for same membrane area as Panel A. 500 spikes were used to estimate the mean and variance, and error bars indicate standard error in the mean for 10 repeated measurements.C: Mean of ISIs for a membrane area of 100µm 2 (6000 Na + and 1800 K + channels).D: CV of ISIs for same membrane area as Panel C. 500 spikes were used to estimate the mean and variance, and error bars indicate standard error in the mean for 10 repeated measurements for all models except the Markov chain model, for which only 4 repeated measurements were used.

Figure 3 .
Figure 3. ISI statistics for noisy input.Same as Fig.2, except that a noisy stimulus was used:I = I DC + 2I noise (t).

Table 1 .
Classification of channel noise models Noise Model Voltage dynamics Subunit dynamics Fraction open Na + channels Fraction open K + channels