# NEURONAL STOCHASTIC VARIABILITY: INFLUENCES ON SPIKING DYNAMICS AND NETWORK ACTIVITY

EDITED BY: Mark D. McDonnell, Joshua H. Goldwyn and Benjamin Lindner PUBLISHED IN: Frontiers in Computational Neuroscience

#### *Frontiers Copyright Statement*

*© Copyright 2007-2016 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.*

*The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.*

*Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.*

*Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.*

*As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.*

> *All copyright, and all rights therein, are protected by national and international copyright laws.*

*The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use.*

ISSN 1664-8714 ISBN 978-2-88919-884-9 DOI 10.3389/978-2-88919-884-9

## About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

## Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

## Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

## What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office<EMAIL_ADDRESS>
## **NEURONAL STOCHASTIC VARIABILITY: INFLUENCES ON SPIKING DYNAMICS AND NETWORK ACTIVITY**

Topic Editors:

**Mark D. McDonnell,** University of South Australia, Australia **Joshua H. Goldwyn,** Ohio State University, USA **Benjamin Lindner,** Bernstein Center for Computational Neuroscience and Humboldt University Berlin, Germany

Stochastic fluctuations are intrinsic to and unavoidable at every stage of neural dynamics. For example, ion channels undergo random conformational changes, neurotransmitter release at synapses is discrete and probabilistic, and neural networks are embedded in spontaneous background activity.

The mathematical and computational tool sets contributing to our understanding of stochastic neural dynamics have expanded rapidly in recent years. New theories have emerged detailing the dynamics and computational power of the balanced state in recurrent networks. At the cellular level, novel stochastic extensions to the classical Hodgkin-Huxley model have enlarged our understanding of neuronal dynamics and action potential initiation. Analytical methods have been developed that allow for the calculation of the firing statistics of simplified phenomenological integrate-and-fire models, taking into account adaptation currents or temporal correlations of the noise.

This Research Topic is focused on identified physiological/internal noise sources and mechanisms. By "internal," we mean variability that is generated by intrinsic biophysical processes. This includes noise at a range of scales, from ion channels to synapses to neurons to networks.

The contributions in this Research Topic introduce innovative mathematical analysis and/or computational methods that relate to empirical measures of neural activity and illuminate the functional role of intrinsic noise in the brain.

**Citation:** McDonnell, M. D., Goldwyn, J. H., Lindner, B., eds. (2016). Neuronal Stochastic Variability: Influences on Spiking Dynamics and Network Activity. Lausanne: Frontiers Media. doi: 10.3389/978-2-88919-884-9

# Table of Contents

*05 Editorial: Neuronal Stochastic Variability: Influences on Spiking Dynamics and Network Activity* Mark D. McDonnell, Joshua H. Goldwyn and Benjamin Lindner

## **Theme 1: Stochastic Variability at the Single-Cell Level**

Cian O'Donnell and Mark C W van Rossum

*34 Diffusion approximation-based simulation of stochastic ion channels: which method to use?*

Danilo Pezo, Daniel Soudry and Patricio Orio

*49 Modeling the influence of short term depression in vesicle release and stochastic calcium channel gating on auditory nerve spontaneous firing statistics*

Bahar Moezzi, Nicolangelo Iannella and Mark D. McDonnell

Aurel A. Lazar and Yiyin Zhou

## **Theme 2: Influence of Variability on Network Dynamics**

*97 Self-consistent determination of the spike-train power spectrum in a neural network with sparse connectivity*

Benjamin Dummer, Stefan Wieland and Benjamin Lindner

*109 Structured chaos shapes spike-response noise entropy in balanced neural networks*

Guillaume Lajoie, Jean-Philippe Thivierge and Eric Shea-Brown

*132 Heterogeneity of heterogeneities in neuronal networks*

Fabiano Baroni and Alberto Mazzoni

*134 A Markov model for the temporal dynamics of balanced random networks of finite size*

Fereshteh Lagzi and Stefan Rotter

## Editorial: Neuronal Stochastic Variability: Influences on Spiking Dynamics and Network Activity

#### Mark D. McDonnell 1 \*, Joshua H. Goldwyn2 and Benjamin Lindner 3, 4

*1 Computational and Theoretical Neuroscience Laboratory, Institute for Telecommunications Research, University of South Australia, Mawson Lakes, SA, Australia, 2 Department of Mathematics, Ohio State University, Columbus, OH, USA, 3 Theory of Complex Systems and Neurophysics, Bernstein Center for Computational Neuroscience, Berlin, Germany, 4 Department of Physics, Humboldt University Berlin, Berlin, Germany*

Keywords: balanced network, channel noise, heterogeneity, Hodgkin-Huxley model, neural networks, neuronal variability, stochastic dynamics

## **The Editorial on the Research Topic**

## **Neuronal Stochastic Variability: Influences on Spiking Dynamics and Network Activity**

Stochastic variability is present across all scales of brain activity. At the single-cell level, for instance, synaptic transmission is mediated by stochastic release of neurotransmitter and membrane potentials fluctuate due to random conformational changes of ion channels. When these cell-level sources of stochastic variability emerge at the network level, they generate fluctuating currents that drive complex network dynamics. Even if intrinsic cellular noise sources are neglected, the interaction of many nonlinear units in recurrent networks typically leads to an effective network noise which is often mathematically tractable in a stochastic framework.

This Research Topic brings together works that address the pressing challenges of developing computational tools and mathematical theories that advance our understanding of stochastic neural dynamics. Six contributions cover stochastic variability at the single-cell level. Moezzi et al. study synaptic coupling between inner hair cells and auditory nerve fibers. Three works update our understanding of ion channel noise in stochastic versions of the Hodgkin-Huxley equations (O'Donnell and Van Rossum; Pezo et al.; Rowat and Greenwood). Puzerey and Galán quantify information transmission in a stochastic Hodgkin-Huxley neuron model that receives barrages of balanced excitatory and inhibitory inputs. Lazar and Zhou communicate a modeling framework that includes dendritic processing of noisy inputs and channel-noise influenced spike generation.

Edited and reviewed by: *Si Wu, Beijing Normal University, China*

#### \*Correspondence:

*Mark D. McDonnell<EMAIL_ADDRESS>
Received: *27 February 2016* Accepted: *01 April 2016* Published: *21 April 2016*

#### Citation:

*McDonnell MD, Goldwyn JH and Lindner B (2016) Editorial: Neuronal Stochastic Variability: Influences on Spiking Dynamics and Network Activity. Front. Comput. Neurosci. 10:38.*

*doi: 10.3389/fncom.2016.00038*

The remaining four studies offer new perspectives on network dynamics. Dummer et al. works out the requirements for self-consistent input/output statistics for neurons embedded in recurrent networks. Lagzi and Rotter develop a Markov chain model that clarifies the stochastic dynamics of balanced networks. Mejias and Longtin explore effects of neural heterogeneity on network response properties. Lajoie et al. make elegant use of random dynamical systems theory to analyse stimulus encoding in in chaotic networks.

Two commentary articles are also part of this research topic: the commentary of Thomas on Lajoie et al. and the commentary of Baroni and Mazzoni on Mejias and Longtin.

## STOCHASTIC VARIABILITY IN SINGLE NEURON DYNAMICS

Moezzi et al. study the sub-cellular origins of long and short-term correlations in inter-spikeintervals (ISIs) in spontaneously firing individual fibers of the auditory nerve, and the form of the distribution of ISIs. They hypothesize the existence of a pre-synaptic mechanism in auditory inner hair cells that causes randomly delays in the availability of synaptic vesicles immediately following prior release of a vesicle. They propose a model for this that produces a simulated ISI distribution equivalent to those in standard existing models. Building on this, the authors tackle the subcellular source of ISI correlations. They report that introducing a model of ion channel noise in the calcium channels of pre-synaptic inner hair cells, when combined with the random vesicle availability model, enables simulation results from the overall model to match the qualitiative nature of ISIs in auditory nerve fiber spontaneous spiking.

Rowat and Greenwood analyse the distributions of interspike intervals generated from simulations of stochastic Hodgkin-Huxley models using several available algorithms. They conclude that appropriate use of the diffusion approximation (e.g., Goldwyn and Shea-Brown, 2011) can yield accurate estimations of these distributions when compared to "microscale" simulations of Markov chain models of ion channel kinetics. The authors draw on their previous studies of Hodgkin-Huxley and related models (Baxendale and Greenwood, 2011; Rowat and Greenwood, 2011) to explain features of the computed interspike interval distributions including local maxima ("bumps") and exponential shape in the tails of of these distributions.

Pezo et al. concur that diffusion approximations can provide accurate and efficient approximations to Markov chain models of channel noise. They caution, however, against using these methods for simulations of fewer than ∼1000 channels. In this low channel number regime, Markov chain methods can be used without sacrificing speed or accuracy. Pezo et al. make the important point that membrane area is subdivided into compartments when simulating dynamics of spatially-extended neurons. The local channel count within each membrane compartment can be small, and thus Markov chain methods may be preferred in multi-compartment computations. Taken together, the Pezo et al. and Rowat and Greenwood articles provide a set of benchmark simulations that will guide researchers interested in studying channel noise in neurons.

O'Donnell and Van Rossum remind us that channel noise depends on the channel type. In stochastic versions of the Hodgkin-Huxley model, for instance, membrane fluctuations are driven by Na+ and K+ channels. Differences in the properties of these ion channels include the probabilities of being in the open states, the time scales of channels opening and closing, and the maximal conductances per channel. O'Donnell and Van Rossum introduce a method to quantify how each channel type contributes to membrane potential fluctuations and how these fluctuations trigger spontaneous action potentials. Importantly, their method can be applied to any channel type in a conductance-based model (as they illustrate with an analysis of a model of CA1 hippocampal neuron). As such, the methods presented by O'Donnell and Van Rossum should be essential to researchers seeking to estimate membrane potential fluctuations induced by diverse channel types.

Puzerey and Galán study responses of the stochastic Hodgkin-Huxley model to barrages of excitatory and inhibitory inputs. They quantify information transmission (measured from estimates of spike train entropy) for a range of synaptic time scales and delays between excitation and lagging inhibition. They find that synaptic kinetics modulate information transmission the most when synaptic currents are balanced and delays in inhibition are small (0.8 ms). The framework of these simulations (neuron model with channel noise, driven by balanced and noisy synaptic inputs) lays a foundation for future studies that will integrate sources of stochastic variability at the neuron and network levels.

Lazar and Zhou seek to identify the theoretical limits of precision with which noisy sensory neurons can encode and decode stimuli. To this end, they introduce an innovative abstract neuron model, comprised from a dendritic module and a spikegenerating module, and study small circuits of such neuron models. The model exhibits a strong grounding in physiological detail; it includes active dendrites; it can use biophysical models such as the Hodgkin-Huxley model for spike generation; and most pertinently to this research topic, includes two intrinsic noise sources. The first source of noise is dendritic variability; the second is ion channel noise due to a finite number of channels (also the focus of O'Donnell and Van Rossum; Pezo et al.; Rowat and Greenwood in this topic). The authors provide extensive mathematical analysis and simulation results, and based on this argue that a duality between stimulus decoding and functional identification holds.

## STOCHASTIC VARIABILITY IN NEURON NETWORK DYNAMICS

The study by Dummer et al. explores an important condition for the temporal correlations of spiking neurons within homogeneous recurrent networks. In such networks, the mean activity of every unit is determined in a self-consistent manner, i.e., its mean input (coming from similar neurons) is simply related to its mean output. Likewise, temporal input correlations are proportional to temporal output correlations and this fact can be used to determine the spike-train correlation function or, equivalently, the spike-train power spectrum, from iterative simulations of a single neuron. Dummer et al. compare two such iterative schemes to simulations of a sparse recurrent network and find excellent agreement. Moreover, their study proves that the emergent network noise can be strongly colored and that the shape of the power spectrum (the "noise color") depends in a nontrivial way on cellular and network parameters.

Although single-cell spike train variability is an important proxy of network stochasticity, it does not necessarily give reliable estimates of the multi-cell spike-train variability that is important for encoding of time-dependent stimuli at the population level. Lajoie et al. study the multi-cell variability of a recurrent network of deterministic quadratic integrate-andfire neurons (a highly chaotic system) under the influence of a set of frozen-noise stimuli. Results from the theory of random dynamical systems are used to estimate an upper bound of the noise entropy. This entropy quantifies the variability caused by different initial conditions of the chaotic system. Lajoie et al. show that, surprisingly, this upper bound of the multi-cell noise

entropy is an order of magnitude lower than a naïve estimate that is based on the single cell's noise entropy. These results may pave the way for the estimation of the mutual information between multicellular spike patterns and complex spatio-temporal stimuli (see also the commentary article of Thomas).

A significant open question in understanding information processing and learning in mammalian cortex is whether the ubiquitous heterogeneity in the anatomy and physiology of cortical neurons and their synaptic connections are necessary for function, or makes little difference. Mejias and Longtin systematically investigate this question for a physiological parameter. They introduce heterogeneity into a standard population model of sparsely-connected excitatory and inhibitory cortical neurons that traditionally assumes homogeneity in all parameters. The heterogeneity in the model takes the form of random distributions for the membranepotential threshold for action potential generation in individual neurons. The authors use simulations and mathematical analysis to explore the effects of heterogeneity, and report intriguing differences that result from whether it is excitatory neurons that are modeled as heterogeneous, or inhibitory neurons. Following publication, a commentary from Baroni and Mazzoni summarized the paper of Mejias and Longtin and discuss the many ways in which the work could be extended. They point out that different heterogeneities could interact in complex unpredictable ways, hence coining the phrase heterogeneity of heterogeneities.

Last but not least, Lagzi and Rotter also consider recurrently connected networks of excitatory and inhibitory neurons. They specifically consider the balanced-state of such networks, where population activity fluctuates but has a steady-state stable mean. Like Dummer et al. and Lajoie et al. the model considered has no external noise, and the fluctuations are a result of the underlying chaotic dynamics of the system. The central contribution of Lagzi and Rotter is to propose a two-state Markovian model which they show exhibits the same phenomenology and statistical properties

## REFERENCES

as simulations of recurrently connected cortical networks in the balanced regime. The authors additionally derive corresponding mean-field equations and provide associated analysis.

## OUTLOOK

The work collected in this Research Topic represents important advances in our understanding of stochastic variability of dynamics in neurons and networks. Methodological issues—such as preferred numerical methods for simulating channel noise in stochastic versions of Hodgkin-Huxley-type models—have been clarified. Moving forward, we expect to see more examples of models that are tightly constrained by known physiological data (as in Moezzi et al.). We have highlighted how variability, as interpreted in these works, can arise from stochastic cellular mechanisms (synapses, ion channels) and from nonlinear and chaotic network dynamics. Of great interest will be future studies that bridge these scales by exploring, for instance, how celllevel noise filters through populations of neurons to shape the dynamics and signal processing capabilities of neural networks.

## AUTHOR CONTRIBUTIONS

All authors listed have made substantial, direct and intellectual contribution to the work, and approved it for publication.

## FUNDING

MM's contribution was supported by an Australian Research Fellowship from the Australian Research Council (project number DP1093425).

## ACKNOWLEDGMENT

We thank all authors and reviewers of manuscripts submitted to this research topic.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 McDonnell, Goldwyn and Lindner. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

## The ISI distribution of the stochastic Hodgkin-Huxley neuron

#### *Peter F. Rowat 1 \* and Priscilla E. Greenwood2*

*1 Institute for Neural Computation, University of California, San Diego, La Jolla, CA, USA 2 Mathematics Department, University of British Columbia, Vancouver, BC, Canada*

#### *Edited by:*

*Benjamin Lindner, Bernstein Center for Computational Neuroscience, Germany*

#### *Reviewed by:*

*Benjamin Lindner, Bernstein Center for Computational Neuroscience, Germany Sten Rüdiger, Humboldt University*

*at Berlin, Germany*

#### *\*Correspondence:*

*Peter F. Rowat, Institute for Neural Computation, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA e-mail<EMAIL_ADDRESS>
## The simulation of ion-channel noise has an important role in computational neuroscience. In recent years several approximate methods of carrying out this simulation have been published, based on stochastic differential equations, and all giving slightly different results. The obvious, and essential, question is: which method is the most accurate and which is most computationally efficient? Here we make a contribution to the answer. We compare interspike interval histograms from simulated data using four different approximate stochastic differential equation (SDE) models of the stochastic Hodgkin-Huxley neuron, as well as the exact Markov chain model simulated by the Gillespie algorithm. One of the recent SDE models is the same as the Kurtz approximation first published in 1978. All the models considered give similar ISI histograms over a wide range of deterministic and stochastic input. Three features of these histograms are an initial peak, followed by one or more bumps, and then an exponential tail. We explore how these features depend on deterministic input and on level of channel noise, and explain the results using the stochastic dynamics of the model. We conclude with a rough ranking of the four SDE models with respect to the similarity of their ISI histograms to the histogram of the exact Markov chain model.

**Keywords: ISI distribution, Hodgkin-Huxley, stochastic dynamics, stochastic differential equation, ISI histogram, Gillespie algorithm, Kurtz approximation**

## **1. INTRODUCTION**

Channel noise is important because it contributes to spike-time variability (Sigworth, 1980; White et al., 2000) and has been shown to be essential for subthreshold oscillations in stellate cells (White et al., 1998; Dorval and White, 2005), and for response variability in sensory cells (Fisch et al., 2012). In addition it contributes importantly to intrinsic irregular firing in cortical interneurons (Englitz et al., 2008; Stiefel et al., 2013), while in certain small neurons a single channel opening can initiate a spike (Lynch and Barry, 1989).

In this paper we compare published SDE approximation methods that simulate the stochastic Hodgkin-Huxley (HH) neuron model, by comparing the inter-spike-interval (ISI) distributions produced when driven by a constant DC current I. Theoretical work on the ISI distributions of stochastic neuron models was carried out by Chow and White (1996); Gerstein and Mandelbrot (1964); Gutkin and Ermentrout (1998); Tuckwell (2005), and Wilbur and Rinzel (1983).

In all cases the deterministic model used as a basis for the various stochastic schemes is the classical model of Hodgkin and Huxley (1952). This model was introduced to describe action potentials in the squid giant axon, and remains a foundation of modern neuroscience. Its dynamics comprise a subcritical Hopf bifurcation together with a switching region in phase space where a fixed point is near to a limit cycle, the two being separated by an unstable limit cycle (**Figures 1A,B**). Thus, the deterministic

HH model has a bistable range: when the input current, *I*, lies between 6.2 and 9.8µA/cm2 (approximately) it is either spiking tonically—represented by the system traversing the locally stable limit cycle—or is quiescent—represented by the system spiraling inside the unstable limit cycle in toward the fixed point. When noise is present, and a trajectory traverses the switching region where the fixed point is close to the stable limit cycle, the system can switch between limit cycle behavior and quiescence. Thus, its overall behavior exhibits irregular switching between bursts of tonic spiking and periods of quiescence. This stochastic behavior continues to occur for a considerable range of the input, I, both below and above the deterministic bistable region *I* = [6.2, 9.8] (Yu and Lewis, 1989; Rowat, 2007).

Here we study the dependence of firing and quiescence patterns on the way noise is modeled, as reflected in the resulting distribution of inter-spike intervals (ISIs). The models we investigate are presented and studied in the papers by Fox and Lu (1994), Fox (1997), Goldwyn et al. (2011), Linaro et al. (2011), Orio and Soudry (2012), and Güler (2013).

In the standard stochastic model for the HH neuron each potassium channel has four binary gates, all of which must be open for potassium to be conducted. Each sodium channel has three activation gates and one inactivation gate, producing eight states, the channel being open only when it is in a particular one of these states [for complete details see

Rowat and Greenwood ISI distribution of the stochastic Hodgkin-Huxley neuron

**FIGURE 1 | (A)** Bifurcation diagram of the deterministic Hodgkin-Huxley model. Solid line, stable fixed point (SFP). Long-dashed line, unstable fixed point (UFP). Dot-dashed lines, extreme values of voltage on the stable limit cycle (SLC). Short-dashed lines, extreme voltage values on the unstable limit cycle (ULC). **(B)** The switching region in the two-dimensional Morris-Lecar model neuron. A similar, but four-dimensional, region is present in the Hodgin-Huxley phase-space. **(B1)** The fixed point is very close to the stable limit cycle with an unstable limit cycle (ULC) between them. Inset: the switching region, enlarged. **(B2)** This shows a noisy trajectory that emits one spike followed by sub-threshold oscillation inside the ULC and then another spike.

Rowat (2007)]. The voltage dependent rates of moving between states have been established from data by Hodgkin and Huxley (1952).

A Markov chain algorithm for keeping track of the number of channels in each state was developed by Chow and White (1996), Gillespie (1977), and Skaugen and Walløe (1979). The algorithm was used by Rowat (2007) to compute several aspects of HH stochastic dynamics. This "exact" method of simulation of the stochastic HH we call the Micro model.

Because of both computation speed and ease of analysis, it is useful to replace the Micro Markov chain model with a system of stochastic differential equations, an SDE model. In fact we were shown how to do this already in a paper of Kurtz (1978), where a system of SDE's is constructed that approximates a density dependent Markov chain at a rate depending on population size N, with error of order log (*N*)/*N*.

Without knowing about the results of Kurtz (1978), authors of a number of papers, Fox (1997); Güler (2013); Linaro et al. (2011), and Orio and Soudry (2012), have devised systems of SDEs to approximate the Micro model. In fact Orio and Soudry (2012), using heuristics, derived the same set of approximating SDEs which a theorem of Kurtz (1978) defines for the Micro model. Here we call this model the Orio-Kurtz or simply Orio model, and sometimes the Kurtz approximation. However, it should be kept in mind that Kurtz proved the approximation for general Markov chain models in 1978. The complete Langevin system of SDEs proposed by Fox (1997), which requires taking two matrix square-roots at every time-step, was implemented for the first time by Goldwyn et al. (2011) so we sometimes refer to this as the Fox-Goldwyn model.

In neuroscience it is widely accepted that the distribution of spike timing, not simply mean spike frequency, is important. While the Micro, Fox-Goldwyn, Güler, Linaro, and Orio-Kurtz models all produce nearly the same mean spike frequency, it is not known how well these models capture the inter-spike-interval (ISI) distributions of the Micro model. Here we generate and compare the ISI distributions of the four SDE models with the ISI distributions of the Micro model, for a range of input current, I, that includes the region of bistability in the deterministic Hodgkin-Huxley model.

We find that, in fact, the ISI distribution is quite similar for all these models over a large range of (constant) deterministic inputs, I, and over a large range of channel numbers, which are proportional to A, the membrane area used. This result is, in fact, expected from analytic considerations, since, arguably, all the models are at least fairly well-approximated by the same diffusion process (Kurtz, 1978), which is based on a approximation theorem with a known error rate and can be regarded as a "gold standard." We will see that both the Fox-Goldwyn and the Orio-Kurtz approximations give ISI distributions which are both quite close to that of the approximated Markov chain model. The Güler model compares well with the Fox-Goldwyn and Orio-Kurtz approximations while the Linaro model is somewhat less successful.

Our second main result is the first detailed description of the form of the ISI distribution of the stochastic HH, which appears in Section 3. In addition we explore how each of the features of the ISI distribution depends on I, the input to the model neuron, and on the number of channels in play, which is functionally related to the standard deviation of the noise in the system.

#### **1.1. FOUR SDE MODELS OF HODGKIN-HUXLEY NOISE**

The current conservation equation for voltage V (mV) and applied current *I* (µA/cm2) in the deterministic HH model is

$$C\dot{V} = I - \left[\mathbf{g}\_{Na}(V - V\_{Na}) + \mathbf{g}\_{K}(V - V\_{K}) + \mathbf{g}\_{L}(V - V\_{L})\right] \quad (1)$$

where the constants are given in **Table 1**.

Equation (1) is the deterministic basis of the four stochastic differential equation (SDE) models we study: the Fox approximation, the Orio-Kurtz approximation, and the models of Güler

**Table 1 | Hodgkin-Huxley parameters for model simulations.**

**Table 2 | Correspondence between si variables and states** *m***jh***k***.**

**Table 3 | Changes in ISI distribution parameters with changes in noise level and applied current.**

(2013) and Linaro et al. (2011). Simulation of the Markov chain model called " Micro" is detailed in Rowat (2007). A potassium channel has four activation *n*-gates, where each gate has the (opening, closing) rates (α*n*(*V*), β*n*(*V*)). The corresponding Markov network is

where a channel in state *ni*, *i* = 0, 1,..., 4, has *i* open *n*-gates. The channel is closed except when in state *n*4.

A sodium channel has 3 activation *m*-gates and one inactivation *h*-gate, where the *m*-gates have (opening, closing) rates (α*m*(*V*), β*m*(*V*)), and the *h*-gates have rates (α*h*(*V*), β*h*(*V*)). The corresponding Markov network is

and the channel is open only in state m3h1. In state m*j*h*k*, there are *j* open *m*-gates and *k* open *h*-gates. The gate transition rates are given by the following functions:

$$\begin{aligned} \alpha\_m(V) &= \frac{0.1(V+40)}{1 - e^{-(V+40)/10}}, & \beta\_m(V) &= 4e^{-(V+65)/18},\\ \alpha\_h(V) &= 0.07 \, e^{-(V+65)/20}, & \beta\_h(V) &= \frac{1}{1 + e^{-(V+35)/10}},\\ \alpha\_n(V) &= \frac{0.01(V+55)}{1 - e^{-(V+55)/10}}, & \beta\_n(V) &= 0.125 \, e^{-(V+65)/80}, \end{aligned} (4)$$

In the Orio and Linaro models one writes an SDE for the proportion of channels in each of the states shown above. Because the K+ channels have 5 states and the Na+ channels have 8 states, these approximations consist of a system of 13 SDEs for Orio or 11 SDEs for Linaro, plus Equation (1), which has no explicit noise term.

In the K+-channel equations that follow, each variable *ni* represents the proportion of channels in state *ni*, for *i* = 0, 1,..., 4. In the subsequent Na+-channel equations, for ease of notation we use variables *s*0, *s*1,..., *s*7 to stand for the proportions of channels in states m0h0, m0h1, ..., m3h1. The correspondence is given in the **Table 2**.

When Equation (1) is integrated, the values of the K+- and Na+-conductances are given by the definitions:

$$\mathbf{g}\_{\mathbf{K}} = n\_4 \mathbf{\bar{g}}\_{\mathbf{K}}, \quad \mathbf{g}\_{\mathbf{Na}} = s\_7 \mathbf{\bar{g}}\_{\mathbf{Na}} \tag{5}$$

Here we write, in algorithmic form, how the SDEs which, in fact, follow from the Kurtz approximation (Kurtz, 1978), are formulated by Orio and Soudry (2012). A few words about the proof of the Kurtz approximation are in Section 1.2. To obtain each of the 13 equations for proportions of channels in each state, we first write the equation as an ordinary differential equation (ODE), thinking of the dynamics for a particular state as deterministic, i.e., the rates are deterministic input and outputs. For each state that is directly linked to the current state, we add to the right hand side a 2-component deterministic term with a positive input component and a negative output component. Then, for each of these deterministic terms on the right hand side, we add a noise term which has the form √*x dW* where *x*, in each case, is the deterministic term with any '−' signs changed to '+,' and *dW* is a Brownian increment. This gives the effective variance when the rates are considered as Poisson rates instead of deterministic rates. In each pair of directly linked states, the stochastic mass going out of one state is the same as the stochastic mass going into the other state. In these cases the terms √*xdW* are kept separate and have opposite signs in the two equations. We see examples in Equations (6) and (7) below. This simple description of the procedure for obtaining the Kurtz approximation (Kurtz, 1978) together with its justification is sketched in Greenwood and Gordillo (2009). Another version is in Orio and Soudry (2012), and its supplement S1.

The full system of SDEs for the potassium channel is given by:

$$
\dot{n}\_0 = (-4\alpha\_n n\_0 + \beta\_n n\_1) + \frac{\dot{\xi}\_1}{\sqrt{N\_K}} \sqrt{4\alpha\_n n\_0 + \beta\_n n\_1}
$$

$$
\dot{n}\_1 = (4\alpha\_n n\_0 - \beta\_n n\_1) + (2\beta\_n n\_2 - 3\alpha\_n n\_1)
$$

$$
$$

$$
\dot{n}\_2 = (3\alpha\_n n\_1 - 2\beta\_n n\_2) + (3\beta\_n n\_3 - 2\alpha\_n n\_2)
$$

$$
$$

$$
\dot{n}\_3 = (2\alpha\_n n\_2 - 3\beta\_n n\_3) + (4\beta\_n n\_4 - \alpha\_n n\_3)
$$

$$
$$

$$
\dot{n}\_4 = (\alpha\_n n\_3 - 4\beta\_n n\_4) - \frac{\dot{\xi}\_4}{\sqrt{N\_K}} \sqrt{\alpha\_n n\_3 + 4\beta\_n n\_4}
\tag{6}
$$

Here ξ*i*, *i* = 1,..., 4 are Gaussian noise terms with mean 0 and standard deviation 1. Note that there are 5 SDEs but only 4 noise terms.

The SDEs for the sodium Markov network (3) can be read off directly from the network using the recipe described above between Equations (5) and (6):

˙*s*0 = (β*ms*2 − 3α*ms*0) + (β*hs*1 − α*hs*0) + ξ20 √*NNa* β*ms*2 + 3α*ms*0 + ξ10 √*NNa* β*hs*1 + α*hs*0 ˙*s*1 = (α*hs*0 − β*hs*1) + (β*ms*3 − 3α*ms*1) − ξ10 √*NNa* β*hs*1 + α*hs*0 + ξ31 √*NNa* β*ms*3 + 3α*ms*1 ˙*s*2 = (3α*ms*0 − β*ms*2) + (2β*ms*4 − 2α*ms*2) + (β*hs*3 − α*hs*2) − ξ20 √*NNa* 3α*ms*0 + β*ms*2 + ξ42 √*NNa* 2β*ms*4 + 2α*ms*2 + ξ23 √*NNa* α*hs*2 + β*hs*3 ˙*s*3 = (3α*ms*1 − β*ms*3) + (2β*ms*5 − 2α*ms*3) − (α*hs*2 − β*hs*3) − ξ31 √*NNa* 3α*ms*1 + β*ms*3 + ξ53 √*NNa* 2β*ms*5 + 2α*ms*3 − ξ23 √*NNa* α*hs*2 + β*hs*3 ˙*s*4 = (2α*ms*2 − 2β*ms*4) + (3β*ms*6 − α*ms*4) + (β*hs*5 − α*hs*4) − ξ42 √*NNa* 2α*ms*2 + 2β*ms*4 + ξ64 √*NNa* 3β*ms*6 + α*ms*4 + ξ54 √*NNa* β*hs*5 + α*hs*4 ˙*s*5 = (2α*ms*3 − 2β*ms*5) + (3β*ms*7 − α*ms*5) + (α*hs*4 − β*hs*5) − ξ35 √*NNa* 2α*ms*3 + 2β*ms*5 + ξ75 √*NNa* 3β*ms*7 + α*ms*5

$$\begin{aligned} -\frac{\dot{\xi}\_{54}}{\sqrt{N\_{\text{Na}}}} & \sqrt{\alpha\_{\text{h}}s\_{4} + \beta\_{\text{h}}s\_{5}} \\\\ \dot{s}\_{6} &= (\alpha\_{\text{m}}s\_{4} - 3\beta\_{\text{m}}s\_{6}) + (\beta\_{\text{h}}s\_{7} - \alpha\_{\text{h}}s\_{6}) - \frac{\dot{\xi}\_{64}}{\sqrt{N\_{\text{Na}}}} \sqrt{\alpha\_{\text{m}}s\_{4} + 3\beta\_{\text{m}}s\_{6}} \\ & + \frac{\dot{\xi}\_{76}}{\sqrt{N\_{\text{Na}}}} \sqrt{\beta\_{\text{h}}s\_{7} + \alpha\_{\text{h}}s\_{6}} \\\\ \dot{s}\_{7} &= (\alpha\_{\text{m}}s\_{5} - 3\beta\_{\text{m}}s\_{7}) - (\alpha\_{\text{h}}s\_{6} - \beta\_{\text{h}}s\_{7}) - \frac{\dot{\xi}\_{75}}{\sqrt{N\_{\text{Na}}}} \sqrt{\alpha\_{\text{m}}s\_{5} + 3\beta\_{\text{m}}s\_{7}} \\ & - \frac{\dot{\xi}\_{76}}{\sqrt{N\_{\text{Na}}}} \sqrt{\alpha\_{\text{h}}s\_{6} + \beta\_{\text{h}}s\_{7}} \end{aligned} \tag{7}$$

Güler (2013) presented a different stochastic Hodgkin-Huxley model. This model also approximates the stochastic dynamics of the membrane potential, arising from random opening and closing of sodium and potassium channels, by a system of seven differential equations, five of them stochastic, together with a modified version of Equation (1), appearing here as Equation (8). The stochastic dynamics, which follow, in a sense, more directly from the approach pioneered by Fox and Lu (1994) than the others, are approximated using carefully constructed diffusion coefficients. In addition, the drift coefficients contain stochastic components, *qK* and *qNa*, designed to capture "non-trivial cross-correlation persistence" (NCCP) effects, namely correlations between transmembrane voltage fluctuations and the component of open channel fluctuations due to gate multiplicities (Güler, 2011). Since properties of the NCCP effects are similar to those of a harmonic Brownian oscillator, the equations that describe *qK* and *qNa* are written as those of a Brownian oscillator. Güler argues that NCCP effects have a major influence on excitability, spontaneous firing, and spike coherence. Güler reports that his model captures very accurately the functional correspondence between input current and mean spike frequency as obtained from the Micro structure (Markov network) model, as well as the mean spike frequency obtained from the Linaro model.

In the Güler SDE model, the current conservation Equation (1) is modified to read:

$$C\dot{V} = -\overline{\mathfrak{g}}\_K \psi\_K (V - V\_K) - \overline{\mathfrak{g}}\_{\text{Na}} \psi\_{\text{Na}} (V - V\_{\text{Na}})$$

$$-\mathfrak{g}\_L (V - V\_L) + I \tag{8a}$$

$$\text{where} \quad \psi\_K = n^4 + \sqrt{\frac{n^4(1-n^4)}{N\_K}} \, q\_{Na} \tag{8b}$$

$$\text{and} \quad \psi\_{\text{Na}} = m^3 h + \sqrt{\frac{m^3 (1 - m^3)}{N\_{\text{Na}}}} \, h \, q\_{\text{Na}} \tag{8c}$$

and the periodic stochastic variables *qK* and *qNa* satisfy two second-order linear SDEs written as four first-order SDEs:

$$
\pi \dot{q}\_{\mathbb{K}} = p\_{\mathbb{K}} \tag{8d}
$$

$$\tau \dot{p} \mathbf{x} = -\varkappa \mathbf{x} p \mathbf{x} - \omega\_K^2 \left[ \alpha\_n (1 - n) + \beta\_n n \right] q \mathbf{x} + \xi \mathbf{x} \tag{8e}$$

$$
\pi \dot{q}\_{\rm Na} = p\_{\rm Na} \tag{8f}
$$

$$\tau \dot{\rho}\_{\rm Na} = -\varphi\_{\rm Na} p\_{\rm Na} - \omega\_{\rm Na}^2 \left[ \alpha\_m (1 - m) + \beta\_m m \right] q\_{\rm Na} + \xi\_{\rm Na} \tag{8g}$$

The gating variables n,m, and h are given by three more SDEs as in Fox and Lu's (1994) paper:

$$
\dot{n} = \alpha\_n (1 - n) - \beta\_n n + \eta\_n \tag{8h}
$$

$$
\dot{m} = \alpha\_m (1 - m) - \beta\_m m + \eta\_m \tag{8i}
$$

$$
\dot{h} = \alpha\_h (1 - h) - \beta\_h h + \eta\_h \tag{8\text{j}}
$$

The Gaussian noise terms have zero means, with variances given by

$$\text{Var}(\xi\_K) = \gamma\_K T\_K \left[ \alpha\_n (1 - n) + \beta\_n n \right] \tag{9a}$$

$$\text{Var}(\xi\_{\text{Na}}) = \wp\_{\text{Na}} T\_{\text{Na}} \left[ \alpha\_m (1 - m) + \beta\_m m \right] \tag{9b}$$

$$Var(\eta\_n) = \frac{\alpha\_n(1-n) + \beta\_n n}{4N\_K},\tag{9c}$$

$$\text{Var}(\eta\_m) = \frac{\alpha\_m (1 - m) + \beta\_m m}{\Re \text{N}\_{\text{Na}}},\tag{9d}$$

$$\text{Var}(\eta\_h) = \frac{\alpha\_h (1 - h) + \beta\_h h}{N\_{\text{Na}}} \tag{9e}$$

where the values of the fixed parameters in the Equations (8e,g) and (9a,b) are:

$$
\gamma\_K = \gamma\_{Na} = 10, \,\omega\_K^2 = 150, \,\omega\_{Na}^2 = 200, \, T\_K = 400, \, T\_{Na} = 800.
$$

The functions α*x* and β*x*, *x* = *n*, *m*, *h*, were given earlier in Equation (4). There are similarities between the Kurtz approximation and the Güler model, e.g., the Güler Equations (9c–e) specify that the diffusion terms of the SDEs (8h,i,j) are similar to the drift terms with '−' changed to '+' just as in the Kurtz approximation. There are significant differences seen in Güler's (8b,c) and in the fact that his SDEs (8d,e) form a second order SDE with a single noise term, and similarly for his SDEs (8f,g). Still it may be that Güler's model is an approximation to the Micro model in the same sense as the Kurtz approximation, or nearly so.

The Linaro model (Linaro et al., 2011) starts from the same current conservation Equation (1), appearing as Linaro et al. (2011; Equation 18). As in the Kurtz approximation, 11 SDEs are introduced Linaro Equation (19), but were obtained through the introduction of Orstein-Uhlenbeck processes for M-1 of the M elements of an M-state Markov process and applying this to the K- and Na- Markov processes. Hence both the drift terms and the diffusion coefficients take a different form from the Kurtz approximation Linaro Equation (19). In view of these differences it is perhaps surprising that the Linaro model produces ISI distributions which are rather close to those produced by the Micro model, the Fox-Goldwyn model, the Orio-Kurtz approximation, and the Güler model.

Diffusion approximations for this stochastic HH Markov chain model have also been studied by Bruce (2009); Goldwyn et al. (2011), and Huang et al. (2013). Engel et al. (2008) and Verechtchaguina et al. (2007) also study ISI histograms for a different modeled neuron and by a different approach.

## **1.2. KURTZ'S STRONG APPROXIMATION THEOREM FOR MARKOV CHAINS**

Here we describe briefly a theorem of Kurtz (1978) and how it applies to approximate the stochastic HH model by the system of SDE's consisting of Equations (1) and (4–7). A more complete version, including an alternate approach using a van Kampen expansion, is described by Baxendale and Greenwood (2011).

In fact one can approximate any normed density dependent Markov process, *XN*(*t*) = *X*(*t*)/*N*, with values in Z*d*, for large population size *N*, by a diffusion process with small error. The method of Kurtz (1978) represents a Z*d* -valued Markov process as a sum of Poisson processes. The essential step is replacing each normed compensated, or conditionally centered, Poisson process with a scalar Brownian motion, where an error of order log (*N*)/*N* is introduced. The resulting stochastic system can be written as

$$d\tilde{X}^N(t) = F(\tilde{X}^N(t))dt + \frac{1}{\sqrt{N}}C(\tilde{X}^N(t))dW(t),\qquad(10)$$

where *F* is the vector field of conditional means of the terms in Kurtz's sum, and the *d* × *d* diffusion coefficient matrix function *C*(*z*) is chosen so that *C*(*z*)*C*(*z*)∗ = *B*(*z*), the covariance function arising from interactions of the terms. One avoids computing the square root of the matrix *B* by retaining the conditional centerings as separate Brownian increments as in Orio and Soudry (2012), Equation (13). We see these terms written out in Equations (5) and (6). This produces a sum of noise terms in each equation so that in distribution the system is the same as (9). The paper (Allen et al., 2008) gives the details of this process.

#### **1.3. THE FORM OF THE ISI DISTRIBUTION**

Suppose we have a recording of membrane potential from a neuron firing in response to a fixed input current I, or we are looking at the output of a simulation of a neuron firing model such as one of those we are considering. An inter-spike interval (ISI) is the time between two successive downward crossings of the recording across a potential level chosen to be well above the range of sub-threshold oscillations. In general the successive ISIs of a simulation of a stochastic model are regarded as independent whereas those of a real neuron are not necessarily so. However, we do not pursue this question here. We are interested only in the distribution of the random ISIs.

The mean spiking frequencies of three models, Micro, Fox-Lu, and Güler are compared for a range of input currents in Figures 6–8 by Güler (2013). These means are not exactly the same but are rather similar. Here we look instead at the entire distribution of ISIs. In Figure 1A of Rowat (2007) we find, already, with Area = 100µm2 and *I* = 0µA/cm2, ISI histograms for the stochastic HH model considered by Chow and White (1996). Figure 16 of Rowat (2007) shows that the histogram of **Figure 1A** is nearly identical to that obtained when Gaussian noise is added to the HH current balance equation for a particular level of noise and a particular constant deterministic input. The effects of carefully modeled channel noise and an equivalent level of Gaussian noise added to the current balance Equation (1) are found to be nearly indistinguishable on the basis of the resulting ISI histograms. This observation motivates the present study where we compare ISI distributions more systematically and for additional recently studied SDE models of the stochastic HH equation.

The form of the ISI histograms indicates that the mean of the ISI distribution is in fact an inadequate parameter to use for comparison of stochastic models. The distribution is not unimodal but instead has the following characteristic form (see **Figures 2**, **3**). For short time intervals there is a tall, narrow peak even on a log scale which represents the distribution of times taken by those individual spike firings which are preceded by one or more spikes, i.e., the times taken by the simulated stochastic path to traverse the locally stable limit cycle of the dynamics when there was a preceding spike. The fact that this first peak is narrow indicates that the variance of the time taken by a stochastic firing is small. The

**FIGURE 2 | Each panel compares the ISI distributions generated by the Micro, Güler, Linaro, Orio, and Fox methods, for a particular area, applied current combination.** Here, the current is constant at 6.0µA/cm2 while the area takes values 100µm2 **(A)**, 500 **(B)**, 900 **(C)**. Equivalently, the noise amplitude decreases from **(A to C)**.

area under this first peak indicates the proportion of ISIs in runs, or "bursts," of two or more spikes. As was found by Rowat (2007), **Figure 4**, the height of this spike increases with the deterministic input, I.

The second obvious feature of the histogram, plotted on a log scale, as in Figure 16B of Rowat (2007) and **Figures 2**, **3**, is that the right tail of the distribution is exponential, as indicated in our linear-log plots by a straight line, starting soon, though not immediately, after the initial peak. The histogram appears increasingly noisy as the length of the ISI time interval increases since the amount of data for estimation decreases, and because short histogram bars are enlarged by the log-scale.

The fact that the right tail of the ISI histogram is exponentially decreasing for the stochastic Morris Lecar model is studied in detail and explained by Rowat and Greenwood (2011).

**FIGURE 3 | Each panel compares the ISI distributions generated by the algorithmic methods Micro, Güler, Linaro, Orio, and Fox, for a particular area, applied current combination.** Here, the area is constant at 400µm2 while the applied current takes values 10.0µA/cm2 **(A)**, 6.0µA/cm2 **(B)**, 2.0µA/cm2 **(C)**.

The explanation is based on the dynamics of the sub-threshold stochastic process which is in a conditional equilibrium beginning soon after a sub-threshold interval begins. The exit time from such an equilibrium must be exponentially distributed. This argument will apply to all stochastic HH models.

A third feature, less obvious but quite distinct, is that there are one or more small bumps or local maxima in the ISI histogram just after the initial peak, and before the exponential tail begins.

The one or more small local maxima in the ISI histogram are explained by considering the dynamics of the stochastic HH model just at the end of firing, i.e., as the stochastic path crosses into the basin of attraction of the locally stable fixed point of the deterministic HH model. The first orbit begins near the outer edge of the basin of attraction—the unstable limit cycle—and so the probability that the next firing (traverse of the locally stable limit cycle) comes after just one sub-threshold orbit is relatively high: see **Figure 1B**. Hence the probability that the ISI ends at the time taken by one such orbit from the end of the previous firing is relatively high. This produces the first small local maximum of the ISI distribution (see **Figures 2**, **3**). Given that the next firing does not occur in this first small orbit of the fixed point, the next subthreshold orbit takes again a similar length of time, and again the probability of firing near the end of this second orbit is somewhat increased, producing the second local maximum, and so on. When there is less noise the pattern is more distinct. After one or two, or at most three such random orbits, the stochastic path is very nearly in its temporary, conditional stationary distribution concentrated close around the fixed point, so the remainder of the ISI distribution is exponentially distributed as discussed above. This description is made more explicit in the Morris Lecar twodimensional system since the "inside" and "outside" of a limit cycle are well defined. In the four-dimensional HH system, the argument is made plausible by examining the relationships in 4 dimensions between the stable limit cycle, the unstable limit cycle, and the fixed point, and by projections onto a 2D plane (e.g., Rowat (2007), Figure 10, and the discussion).

The ISI distribution of any stochastic HH neuron model can, therefore, be resolved into three sections, an initial peak representing the distribution of times taken up by contiguous spike firings, followed by one or more local maxima representing the additional times taken up by each of a few subthreshold orbits of the fixed point immediately after firing, followed by an exponential tail representing time until escape from the subthreshold state once the process is in its conditionally stationary distribution. Thus, a complete comparison of ISI histograms for simulations of stochastic HH models built from different noise models can be made by comparing the defining parameters of these three components: the center, height, and width of the initial peak, the shapes and placement of the local maxima, and the parameters of the exponential tails. We can use these criteria for comparing the ISI histograms produced by the four stochastic HH models described in Section 1.1.

## **2. METHODS AND IMPLEMENTATION DETAILS**

All model computer runs used the standard Hodgkin-Huxley parameter values, as in **Table 1**, and all data sets used for the ISI histograms had 105 elements. All histograms were normalized so that their bars sum to 1, and all are displayed with a log scale on the y-axis because the first peak is often an order of magnitude higher than the second and third peaks. The integration time-step was 0.005 ms. An SDE was used for each state. If a potassium variable became negative the random number generator was called again. If a sodium variable became negative it was immediately reset to zero. At the end of each integration step, but before integrating Equation (1), the potassium variables *ni*, *i* = 0,..., 4 were normalized to satisfy *ni* = 1, and the sodium variables *sj*, *j* = 0,..., 7 were normalized by *sj* = 1. This seems more correct than defining the last variable (*n*0 or *s*0) in terms of the others, since it preserves the relative values of the variables, but has the disadvantage of using two extra SDEs. However, for the Orio method and any Kurtz-type approximation it does not increase the number of random number generator calls.

## **3. RESULTS**

In **Figures 2**, **3** we see ISI histograms from simulations of the Markov chain model and the four SDE models, which are labeled Micro, Güler, Linaro, Orio, and Fox, to refer to the five ways of modeling HH channel noise described in Section 1.1. The noise level is proportional to Area−1/2 since the standard deviation of the Na+-channel noise (K+-channel noise) is ∝ *N*−1/2 *Na* ∝ *N*−1/2 *K* and the number of channels is a constant times the area. In **Figure 5** where the area A = 400µm2, *NNa* = 60 × *A* = 24000 and *NK* = 18 × *A* = 7200 Separate sets of plots show the results for fixed applied current *I* = 6µA/cm2 and area *A* = 100, 200, ..., 1000 µm2, and for fixed area A = 400µm2 and applied currents *I* = 2, 3,..., 12µA/cm2. See **Figures 4**, **5**. All the histograms show the features detailed in Section 3: an initial peak, followed by local maxima, followed by an exponential tail, that appears linear on a log scale. The Güler, Micro, Orio, and Fox histograms are nearly identical in all respects and the Linaro plots are very similar.

The histograms computed by the Güler method all have a slightly higher proportion of area under the initial peak, corresponding to runs of successive spikes, than Micro-generated histograms, while Orio-generated histograms always have slightly

lower proportion of ISIs in runs than Micro (differing by no more than 4%), while the Linaro histograms have proportions of ISIs in runs that are 8% lower than in Micro histograms.

Since the plots for the different SDE models are so similar it may not be worth dwelling on the differences. The exponential parameters of the tails are very close as evidenced by the nearly parallel plots. Also the timing of all features is nearly identical, showing that the different ways of modeling noise by SDEs have had little effect on the pattern of firing of the simulated neuron.

In **Figures 2**, **3** we see what happens to the three parts of the ISI histogram as the area A increases and the applied current *I* increases. Since A is proportional to the numbers of channels, in fact the standard deviation of the noise per unit time is proportional to A−1/2. We see in **Figure 2** that the slope of the exponential part of the histogram, i.e., the slope of the last part, becomes less negative as A increases, equivalently, as the noise decreases. The height of the initial peak barely changes, but we see that the bumps in the middle part of the histogram become more prominent as noise decreases. When A is fixed and I, the deterministic input, decreases, in **Figure 3** we see that the negative slope, i.e., the negative exponent of the exponential tail also decreases but, opposite to the case when noise decreases, the bumps (only one can be seen) become less prominent. The height of the initial peak decreases considerably, and its position moves right as *I* decreases. These effects are studied in more detail for the Orio model in **Figures 4**, **5**. We discuss their interpretation in the next section, but give a summary in **Table 3**.

## **4. DISCUSSION**

The models we have simulated produce similar histograms with the same basic features in good alignment. Here we discuss further how these features depend on two important parameters of the stochastic Hodgkin-Huxley model, the deterministic input, I, and the strength of the stochastic input which is proportional to A−1/2. Notice that these two parameters can be regarded as measures of deterministic and stochastic input, respectively.

First let us focus on how these two inputs affect the parameter of the ISI tail distribution. We see from the log-linear plots

**FIGURE 6 | Comparison of exponential tail exponents generated by the five methods, for fixed current** *I* **= 6***.***0 µA/cm2, plotted as a function of Area−1***/***2 where Area takes the values 1000, 900, . . . , 100** µ**m2.**

in **Figures 2**–**5** that the negative slope of the final segment of the ISI distribution, which is the negative exponent of the exponential tail of the distribution, increases with increasing deterministic I, as well as with increasing stochastic input, A−1/2, i.e., these become steeper with decreasing A and with increasing I. In **Figures 6**, **7** the tail exponents are plotted as functions of A−1/2 and of I, respectively.

To understand this result we return to the state space picture of the stochastic dynamics of the neuron model, represented by **Figure 1B2**, which shows the dynamics for the analogous Morris-Lecar model. Between firings the state of the neuron is in the subthreshold region centered on the fixed point but well inside the unstable limit cycle except while traversing the switching region. The probability that it moves out of this region and fires is greater if either the noise has greater standard deviation or if the deterministic input to which the noise is added is greater. Hence an increase of either *I* or A−1/2 should have a similar effect on the parameter of the ISI exponential tail distribution. Furthermore, as *I* changes the configuration pictured in **Figure 1B** changes. As *I* moves toward the bifurcation at *I* ≈ 9.8, the bifurcation diagram in **Figure 1A** shows us that the stability of the fixed point decreases, becoming zero at *I* = 9.8, while the unstable limit cycle shrinks and disappears. Correspondingly, the subthreshold region shrinks in size, but does not disappear since one sees short intervals of subthreshold behavior for values of *I* ranging at least as high as 12.0. In **Figure 8** one might note that when *I* = 12 none of the curves have reached 1. Both these effects cause the probability of firing to increase. Reduction in stability means it is easier to escape from the fixed point, while reduction in size of the unstable limit cycle means the size of the subthreshold regime is smaller, thus reducing the expected time to reach the deterministic basin of attraction of the stable limit cycle. The combination of these effects seem to cause the relation between *I* and the exponential tail to be concave, as in **Figure 7**, instead of nearly linear as in the case of noise, as in **Figure 6**. Switching to spiking and maintenance of spiking become more probable, and the exponential tail of the ISI distribution becomes steeper.

Next we consider the effect of increasing *I* or A−1/2 on the bumps in the middle part of the ISI histograms. To understand this we need to review some theory about the behavior of a similar neuron model during its subthreshold phase. This was studied in detail for the Morris-Lecar model by Ditlevsen and Greenwood (2013). The same analysis applies here with some alterations because the deterministic HH model is 4-dimensional. If we linearize the stochastic model at the fixed point we obtain a linear stochastic system of the form

$$dX = -AX\,dt + C\,dW.\tag{11}$$

The deterministic matrix C is obtained by evaluating the stochastic diffusion coefficient matrix at the fixed point. The matrix A will have a pair of complex eigenvalues, −λ ± *i*ω with negative real part λ ω and other negative eigenvalues −β, −γ where β, γ are greater than λ. This means the system, started near the fixed point, moves rapidly toward the plane defined by the eigenvectors corresponding to the complex eigenvalues. Thus, the 4-dimensional system can be studied in terms of a 2-dimensional system. Other examples are found in Baxendale and Greenwood (2011).

When the neuron fires and then becomes subthreshold, the stochastic path enters the region "inside" the unstable limit cycle at its edge and proceeds to roughly circle the fixed point at a frequency ω, the orbit being damped at a rate λ while also being restrained from damping by the stochastic aspect of the

model. The path of this process can be approximated in terms of a fixed rotation multiplied by a 2-dimensional Ornstein-Uhlenbeck process as shown by Baxendale and Greenwood (2011). The sample path of the process orbits the fixed point for some time, with frequency ω, until it arrives at its stationary distribution. Before stationarity sets in we see one or more decreasing "bumps" in the ISI histogram, with frequency ω, and after stationarity sets in we have the exponential tail, being the escape distribution from a stationary distribution by a standard argument.

As examples, according to computation by Hassard (1978): for *I* = 5, λ ± *i*ω = −0.097 ± 0.521*i*, −β = −0.129, −γ = −4.60; for *I* = 9, λ ± *i*ω = −0.015 ± 0.578*i*, −β = −0.137, −γ = −4.73. We find that the spacing between the second and third bumps in **Figure 4**, where *I* = 6, and also for larger areas (not shown), is approximately 12 ms, which is in rough agreement with 2π/ω ≈ 12.05 ms for *I* = 5 above. It is notable that the eigenvalue frequency ω, and thus the bump spacing, bears no particular relationship to the frequency of the unstable limit cycle (ULC), as computed by Rinzel and Miller (1980). Let *I*1 ≈ 9.8 be the subcritical Hopf bifurcation. For *I* close to *I*1, *I* < *I*1, the eigenvalue frequency and the ULC frequency are the same at approximately 90 Hz, but as *I* decreases toward *Iv* ≈ 6.26, ω decreases by only 9% while the ULC frequency decreases steeply from 55% from 90 Hz at *I*1 to approximately 40 Hz at *I* = 7.5 then smoothly reverses direction and increases back up to about 50 Hz at *I*υ. One might also note the reduction in λ as *I* increases from 5 to 9, while β and γ are both much larger than λ, as predicted by Equation (10).

Here we make a comment on the existence of exponential tails for *I* > *I*1. The underlying mechanism of this has been discussed in Rowat and Greenwood (2011). Numerically, it has been shown by Rowat (2007) and Tateno and Pakdaman (2004) that the probability *p*(*I*) that a spike is followed by a non-spike is continuous across *I*1. Note that 1 − *p*(*I*) is the proportion of ISIs in runs of two or more spikes (see **Figures 8**, **9**).

When *I* > *I*1, say with *I* = 11 or *I* = 12, the equilibrium point is now unstable with dominant eigenvalues λ ± *i*ω, where λ is small, 0 < λ ω. Since λ is small one sees numerically that a deterministic trajectory started very close to the unstable

**FIGURE 9 | Proportion of ISIs in runs of two or more spikes, compared by method, for fixed current** *I* **= 6***.***0 µA/cm2, as a function of Area−1***/***2.** Area takes values 1000, 900, . . . , 100µm2.

equilibrium makes several very tight, small, slowly expanding spirals around the equilibrium before switching out to the stable limit cycle (SLC)—i.e., the spiking cycle.

When *I* is below the Hopf bifurcation, with dominant eigenvalues λ ± *i*ω where λ < 0 and −λ ω, Baxendale and Greenwood (2011) identify the stochastic process whereby deterministic damped oscillations, with the addition of noise, show sustained oscillations at an amplitude well above the expected noise level. The exit time from this stochastic equilibrium process is what creates the exponential tail when *I* < *I*1.

In view of the observation above, that when *I* > *I*1, several small tight spirals may occur before a deterministic trajectory, started close to the unstable FP, switches out to the SLC, it seems reasonable to propose that in the presence of noise there is a short-term stochastic process that tends to contract the slowly expanding deterministic spirals, thus creating a conditional equilibrium for a short time before the trajectory switches out to the SLC. Thus, the exit time from this conditional stochastic equilibrium has an exponential distribution that creates the exponential tail when *I* > *I*1.

The number and pronounced definition of the bumps become less as the noise increases because the onset of stationarity is hastened by more noise. We see this effect in **Figure 4**. In **Figure 5** we observe that changing *I* with the noise level fixed has much less effect on the size and definition of the "bumps."

Finally, how do the parameters *I* and *A*−1/2 affect the height, width, and location of the main peak of the ISI histogram? In **Figure 8** we have plotted the proportion of ISI mass in runs of two or more spikes as a function of input *I* for fixed A = 400µm2. This is represented in the histogram as the area under the first peak. It increases roughly linearly with *I* except for saturation at 0 and 1. The main peak moves right as *I* decreases and as A increases.

Note that any occurrence of a run of two or more spikes corresponds to the occurrence of a spike immediately followed by another spike, hence the proportion of ISI histogram mass under the first peak is in fact the probability that a spike is followed by another spike. Equivalently, the histogram mass or area under the tail (including any "bumps") is the probability that a spike is followed by a period of quiescent behavior.

**Figures 8**, **9** show that the proportion of ISIs occurring in runs of two or more spikes increase roughly linearly with *I* and with *A*−1/2, respectively, when the other variable is fixed. This is reflected in the ISI distribution as an increase in the height and width of the initial peak. The reasons for increasing steepness of the exponential tail apply equally to the increase we observe in **Figures 8**, **9**. A large negative tail exponent implies less area, or "probability mass," under the tail and more mass in the main peak. Thus, the increases seen in **Figures 8**, **9** correlate well with the increases in negative tail exponent in **Figures 6**, **7**.

In **Table 4**, we give numerical values for the exponential tail exponent, the proportion of ISIs in runs of two or more spikes, and the running times, across all the models, for one specific (Area, I) combination, namely Area = 400µm2, *I* = 6.0µA/cm2. These computer simulations were all run consecutively on the same hardware. We see that for these parameters, the Fox-Goldwyn and Orio-Kurtz methods are equally close (within a few percent), the Güler method a little further away, and the Linaro method further away again.

ISI densities were also computed by Verechtchaguina et al. (2007) and Engel et al. (2008) by a different method and for a different neuron. An electrical circuit was used to capture the frequency-dependent subthreshold dynamics in stellate and pyramidal cells of the entorhinal cortex, which was converted to a noise-driven harmonic oscillator; from this they analytically computed ISI densities.

## **5. CONCLUSION**

**Figures 6**–**9** and **Table 4** together show that the Fox-Goldwyn, and Orio-Kurtz methods both generate ISI histograms very close to those of Micro. The Güler histograms are not quite as close and the Linaro histograms are only a little further off.

According to Kurtz' theorem the Orio method gives an error of at most log (*N*)/*N* which is 0.0001 for the data sets computed here (*N* = 105). Hence it should be regarded as a "gold standard" for producing a good approximation to the ISI distribution of the Markov chain model. However, when computation time is an issue, one might well prefer to use the Güler model which runs about three times as fast as the Linaro and Orio models. This was true for our Python implementations on a 2.5 GHz Intel Core i5 and will no doubt generalize to other languages and systems. We used the same basic code framework for the Güler, Linaro, and Orio methods. The main reason for the increased speed is that

*Column 2 lists the exponential exponent of the histogram tail, while column 3 lists the probability that a spike is immediately followed by another spike. This is equivalent to the proportion of spikes that occur in a run of two or more spikes. In the lower four rows of columns 2 and 3, the first figure is the actual parameter value while the second figure is its difference from the corresponding value for the Micro simulation. The fourth column gives the running time on our hardware and the fifth the computer language used. Both the Fox-Goldwyn and Linaro code were retrieved from the ModelDB repository.*

the Güler simulation calls the random number generator much less often than the others. In addition, the Güler method uses considerably fewer algebraic operations. Unfortunately the Fox-Goldwyn model was implemented in Fortran so its computation time cannot reasonably be compared with the other three SDE models.

Although the Güler method generates histogram parameters further away from the Micro histogram parameters than either the Fox-Goldwyn or Orio-Kurtz histogram parameters, one must bear in mind that when introducing harmonic Brownian oscillator-type SDEs, there are six phenomenological parameters in the Güler method that were carefully chosen by examination of simulations of Micro voltage data in a subthreshold regime, with *I* = −4µA/cm2 (to avoid spikes). It may be that if these parameters were chosen with reference to Micro simulation voltage data generated with another *I*-value, e.g., in the middle of the bistability interval [6.2, 9.8], the parameters of the Güler histograms could be much closer to the parameters of the Micro histograms.

## **REFERENCES**

Yu, X., and Lewis, E. R. (1989). Studies with spike initiators: linearization by noise allows continuous signal modulation in neural networks. *IEEE Trans. Biomed. Eng.* 36, 36–43. doi: 10.1109/10. 16447

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 26 May 2014; accepted: 25 August 2014; published online: 08 October 2014.*

*Citation: Rowat PF and Greenwood PE (2014) The ISI distribution of the stochastic Hodgkin-Huxley neuron. Front. Comput. Neurosci. 8:111. doi: 10.3389/fncom. 2014.00111*

*This article was submitted to the journal Frontiers in Computational Neuroscience. Copyright © 2014 Rowat and Greenwood. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Systematic analysis of the contributions of stochastic voltage gated channels to neuronal noise

#### *Cian O'Donnell 1,2\* and Mark C. W. van Rossum2 \**

*1 Computational Neurobiology Laboratory, Salk Institute for Biological Studies, La Jolla, CA, USA*

*2 School of Informatics, Institute for Adaptive and Neural Computation, University of Edinburgh, Edinburgh, UK*

#### *Edited by:*

*Benjamin Lindner, Bernstein Center for Computational Neuroscience, Germany*

#### *Reviewed by:*

*Andre Longtin, University of Ottawa, Canada Eric Shea-Brown, University of Washington, USA Richard Naud, University of Ottawa, Canada*

#### *\*Correspondence:*

*Cian O'Donnell, Computational Neurobiology Laboratory, Salk Institute for Biological Studies, 10010 N Torrey Pines Road, La Jolla, CA 92037, USA e-mail<EMAIL_ADDRESS>Mark C. W. van Rossum, School of Informatics, 10 Crichton Street, Edinburgh EH8 9AB, UK e-mail<EMAIL_ADDRESS>
## **1. INTRODUCTION**

An obvious characteristic of behavior is the variability that one observes from trial to trial in even the most controlled settings. This behavioral variability is reflected at the neural level in the noisy character of spike trains. Various hypotheses have been put forward for a potential functional role of neural variability, such as stochastic resonance (McDonnell and Abbott, 2009), prevention of synchrony (van Rossum et al., 2002), and probabilistic sampling (Buesing et al., 2011). A number of factors can contribute to trial-to-trial variability: non-stationarity and unobserved modulation of the nervous system; chaotic network dynamics resulting from deterministic single neuron dynamics (van Vreeswijk and Sompolinsky, 1996); and biophysical noise. In this paper we concentrate on the latter, and in particular on the noise from voltage-gated ion channels.

Ion channels are pore-forming macromolecular proteins that allow the selective passage of ionic currents in and/or out of the cell (Hille, 2001). Each ion channel can, at any given moment, occupy only one of multiple discrete conformational states; at least one of which is an open/conducting state, and at least one of which is a closed/non-conducting state. Transitions between states are exceedingly rapid (<1µs) and, like all molecular reactions, stochastic in nature—they are driven by thermal agitation. In the case of voltage-gated channels (VGCs) considered here, the transition probabilities depend on the cell's membrane potential. Channels are commonly modeled as Markov processes, which

Electrical signaling in neurons is mediated by the opening and closing of large numbers of individual ion channels. The ion channels' state transitions are stochastic and introduce fluctuations in the macroscopic current through ion channel populations. This creates an unavoidable source of intrinsic electrical noise for the neuron, leading to fluctuations in the membrane potential and spontaneous spikes. While this effect is well known, the impact of channel noise on single neuron dynamics remains poorly understood. Most results are based on numerical simulations. There is no agreement, even in theoretical studies, on which ion channel type is the dominant noise source, nor how inclusion of additional ion channel types affects voltage noise. Here we describe a framework to calculate voltage noise directly from an arbitrary set of ion channel models, and discuss how this can be use to estimate spontaneous spike rates.

**Keywords: channel noise, voltage-gated ion channels, Hodgkin–Huxley, spontaneous firing, simulation**

lead to accurate predictions of the noise in macroscopic currents recorded from neurons (Hille, 2001).

Because spike generation appears reliable during somatic current injection (Calvin and Stevens, 1968; Bryant and Segundo, 1976; Mainen and Sejnowski, 1995) and the number of VGCs is large, it is typically believed that the stochastic gating of VGC contributes little to the total observed variability in neuronal spiking. However, such a conclusion might be premature. First, during somatic current injection the dendrites are typically more hyperpolarized compared to the realistic case where the neuron receives synaptic input, and the noise typically decreases strongly with hyper-polarization (see below). Moreover, the total number of channels in a neuron might be large, but in spatial compartments such as narrow axons or dendrites, the number of channels is typically small.

Several experimental studies have focused on the physiological consequences of ion channel noise. Sigworth (1980) used fluctuation analysis to estimate the number of Na+ channels at a single frog node of Ranvier ∼30000, and subsequently used formulae from Lecar and Nossal (1971a,b) to estimate fluctuations in the current threshold of action potential generation due to channel noise. Johansson and Arhem (1994) found that the stochastic opening of a small number of channels in cultured hippocampal neurons were sufficient to trigger spontaneous action potentials. White et al. (1998) recorded subthreshold membrane potential oscillations in stellate cells of layer II entorhinal cortex (EC) and found that they could only reproduce the co-existence of both oscillations and spiking in a computational model if they included stochastic gating of Na+ channels, suggesting a form of periodic stochastic resonance. Subsequently, Dorval and White (2005) used the dynamic clamp technique to inject a "virtual" Na+ conductance which was either deterministic or stochastic to EC stellate cells *in vitro*. Only the stochastic conductance could reproduce the observed membrane potential oscillations. Similarly, Dudman and Nolan (2009) used computational models of the same cell type to demonstrate that stochastic channel gating can also account for the clustered firing patterns displayed by these cells when stimulated by steady current pulses *in vitro*. Diba et al. (2004) characterized somatic subthreshold voltage noise in cultured hippocampal neurons due to stochastic ion channel gating. Voltage fluctuations were small, with a standard deviation <0.3 mV and based on pharmacological experiments appeared to arise primarily from K+ channels. Jacobson et al. (2005) reported similar results from neocortical pyramidal cells from layer IV/V of rat somatosensory cortex brain slices, with similar amplitude (submillivolt) voltage fluctuations. Finally, Kole et al. (2006) used fluctuation analysis to measure the properties and distribution of hyperpolarization-activated cation (Ih) channels in LV neocortical pyramidal cells *in vitro*. They found that although the Ih single-channel conductance was exceedingly small (<1 pS), Ih channels contribute substantially to voltage noise in the distal dendrites of these cells.

A great deal of theoretical and numerical studies have looked at membrane noise from stochastic ion channels, beginning with Lecar and Nossal, who used stochastic differential equations and a reduced dynamical system model of the action potential to attempt to quantify the magnitude of membrane noise on action potential threshold fluctuations (Lecar and Nossal, 1971a,b). Skaugen and Walløe (1979) were the first to examine the consequences of stochastic gating of ion channels through numerical simulation. They found that in the stochastic version of Hodgkin– Huxley (HH) squid giant axon model the current threshold was lowered compared to deterministic models, the membrane could spike spontaneously, and that the frequency-current curve was smeared around the threshold. Subsequent simulation work by DeFelice and colleagues (Clay and DeFelice, 1983; Strassberg and DeFelice, 1993) further elaborated on the direct link between the microscopic (stochastic) and macroscopic (deterministic) versions of the HH model. Rubinstein (1995) simulated a model of the frog node of Ranvier and reproduced the spread in action potential firing threshold due to stochastic channel gating predicted by Lecar and Nossal, in agreement with earlier experiments (Verveen, 1960). Chow and White (1996) examined the dependence of spontaneous firing rate in the stochastic HH model on membrane patch area and found it to decrease exponentially with area. They approximated the system as a boundary escape problem, with stochastic gating of the activation subunit of the Na+ channel as the noise source. They found that the mean escape time as a function of area agreed well with numerical simulation results (we will comment on this finding below). Manwani and Koch (1999) used a perturbative approach to compare the contributions of thermal noise, channel noise and synaptic "noise" (from Poissonian inputs) to total membrane voltage noise in a single compartment. Steinmetz et al. (2000) used similar methods to demonstrate the voltage and channel type dependence of ion channel noise spectra for both the HH model and a commonly used neocortical pyramidal cell model (Mainen et al., 1995). In the present study we employ methods very similar to both of these works, but toward a different goal: we aim to systematically separate all of the contributing factors that determine the contribution of an ion channel type to voltage noise and spontaneous firing.

In general detailed simulation of stochastic channels will give the most accurate answer regarding the noise and the contribution of the different channels. But as stochastic simulation of the full channel kinetics is very involved, several recent studies have developed approximate stochastic-differential equation models that efficiently capture the essence of the noise statistics of discrete ion channels (Goldwyn and Shea-Brown, 2011; Goldwyn et al., 2011; Linaro et al., 2011; Orio and Soudry, 2012). Our objective here is complementary but different: rather than developing a precise model for the noise we seek to *estimate* the contribution of the various channel types. Intuitively it is not clear what properties of a given channel type are important determinants for noise. This is relevant when a full state diagram of a certain channel type is not available, but nevertheless a coarse estimate of its contribution to noise is desired. At the same time by breaking down the various factors that determine the magnitude of the noise of a certain channel type, a deeper insight in the results from simulations and experiments can be obtained.

Our study is split in four parts: We use simulations to demonstrate that in the well-studied stochastic Hodgkin–Huxley model, most spikes are due to stochastic K+ channel gating and not the Na+ channel. This is of interest as in the literature conflicting findings can be found (see Discussion). Next, we review the different factors that explain how stochastic channel noise leads to noise in membrane voltage, using a linear, weak noise analysis and explain why the K+ noise is dominant. While these individual contributing factors are well-known, a concise account was in our opinion lacking. Third, we examine the relation between voltage noise and spontaneous spikes using an approach recently introduced for integrate-and-fire neurons. We show that this relation is complex, but that nevertheless rough estimates are possible. Finally, we apply the same methods to analyze a CA1 pyramidal neuron model to show that the approach is easily generalizable to other neuron models.

## **2. MATERIALS AND METHODS**

All stochastic channel simulations were implemented using the Parallel Stochastic Ion Channel Simulator (PSICS) which is specifically designed for efficient simulations of stochastic ion channel gating in single neuron models (see Cannon et al., 2010 and http://psics.org/ for algorithmic details). Current noise injection in deterministic HH models was done using NEURON (Carnevale and Hines, 2006). Analysis was done using MATLAB (The Mathworks).

## **2.1. HODGKIN–HUXLEY MODELS**

All simulations of the Hodgkin–Huxley model used the standard voltage-dependent equations for Na+ and K+ gating schemes, at 6.3◦C (Hodgkin and Huxley, 1952). For completeness the standard parameters are given in **Table 1**.

In the HH squid axon model, the sodium conductance obeys

$$\mathcal{g}\_{\rm Na}(V, t) = \mathcal{y}\_{\rm Na} \rho\_{\rm Na} A m^3(V, t) h(V, t)$$

where γ*Na* is the single channel sodium conductance, ρ*Na* is the density of channels per area, and *A* is the area. The gating variables *m* and *h* move between off and on-states with a voltage-dependent rate α*m*, *h*(*V*) from the off-state to the on-state, and back with a rate β*m*, *h*(*V*). These rates have been empirically established as

$$\alpha\_m(V) = \frac{0.1(V+40)}{1 - e^{-(V+40)/10}} \quad \beta\_m(V) = 4e^{-(V+65)/18}$$

$$\alpha\_h(V) = 0.07e^{-(V+65)/20} \quad \beta\_h(V) = \frac{1}{1 + e^{-(V+35)/10}}$$

where *V* is membrane voltage in mV and the transition rates have units 1/ms.

In the limit of very many channels the gating variables are a continuous quantity, namely the probability to find them in the on-state. The dynamics of gating variable *m*(*V*, *t*) can be written as

$$\frac{dm(V,t)}{dt} = \frac{m\_{\infty}(V) - m(V,t)}{\mathfrak{r}\_{m}(V)}\tag{1}$$

where *m*∞(*V*) is the steady-state value for the activation variable,

$$m\_{\infty}(V) = \frac{\alpha\_m(V)}{\alpha\_m(V) + \beta\_m(V)}\tag{2}$$

and τ*m*(*V*) is its time constant

$$\pi\_m(V) = \frac{1}{\alpha\_m(V) + \beta\_m(V)}\tag{3}$$

and analogous for the inactivation variable *h*(*V*, *t*). The dynamics in the continuum limit is fully deterministic.

When the system is described stochastically, the gating variables of each channel are binary variables that switch between off-state and on-state. For the Na+ channel to be open, all its 3 *m*'s and the *h* switch need to be in the on-state. To (naively) simulate

**Table 1 | The Hodgkin–Huxley parameters for model simulations.**

this case the transitions are drawn stochastically using a random number generator, using a time-step δ*t* such that αδ*t*,βδ*t* 1.

Likewise, the K+ conductance is given by

$$\mathbf{g}\_{\mathcal{K}}(V, t) = \mathcal{Y}\_{\mathcal{K}} \rho\_{\mathcal{K}} \mathcal{A} n^4(V, t)$$

It has four identical activation variables, labeled *n*, with rates

$$\alpha\_n(V) = \frac{0.01(V+55)}{1 - e^{-(V+55)/10}} \quad \beta\_n(V) = 0.125e^{-(V+65)/80}$$

In all simulations the single channel conductance for both Na+ and K+ was 20 pS. Although this value is close to that reported experimentally for the squid giant axon K+ conductance (Llano et al., 1988), it is slightly larger than experimental estimates for the Na+ conductance (Bezanilla, 1987). These values were chosen for simplicity (it removes one confounding factor when comparing channel type noise contributions) and to enable comparison with the literature (Strassberg and DeFelice, 1993; Chow and White, 1996; Schneidman et al., 1998). In line with the literature, leak channels were modeled deterministically, although in more realistic models they should be made stochastic as well.

For the simulations and analysis of the hippocampal pyramidal cell model (**Figure 6**), we use the channel models for active Na+, delayed rectifier K+ (Kdr), and A-type K+ channel (Ka) exactly as previously published by Migliore et al. (1999), Jarsky et al. (2005). However, our model was single-compartment while these previous studies looked at multi-compartment model neurons. For consistency with the HH simulations we also choose a single-channel conductance of 20 pS. The channel densities were matched to the macroscopic conductance densities of the soma in the model of Jarsky et al. (2005), implying Na+: 20 channels /µm2, Kdr: 20/µm2, and Ka: 24/µm2. In addition to these active channels, we added two voltage-independent leak channels, one permeable to Na+ and one permeable to K+, which we simulated deterministically. We chose their densities 0.0065 mS/cm2 (Na+) and 0.0185 mS/cm2 (K+), to fit two constraints: a total leak conductance of 0.025 mS/cm2 (Jarsky et al., 2005), and a resting voltage of −65 mV. Finally, as according to Migliore et al. (1999) we set reversal potentials *ENa* = +55 mV and *EK* = −90 mV. When attempting to analytically calculate the membrane impedance for this model we unfortunately found that it diverged upward at around −60 mV. This singularity is problematic because it would break our small voltage noise assumption. Hence we instead estimated the impedance for this model empirically as is done in experiments. We injected sine wave currents over a large frequency range to the deterministic version of the model, and measured the amplitude of the resulting voltage responses.

#### **2.2. CALCULATION OF POWERSPECTRA OF K+ AND NA+ NOISE**

As is well known, the current noise power spectrum from a population of ion channels follows directly from the channel kinetic scheme and its autocovariance function (DeFelice, 1981). For ease of presentation we first summarize this calculation for a two-state channel. The conditional probability that a two-state channel is open at time *t* given that it was open at time 0, *po*|*o*(*t*) is

$$p\_{o|o}(t) = p\_{\infty} + (1 - p\_{\infty})e^{-t/\tau}$$

where *p*∞ is the steady-state open probability and τ is the correlation time. The autocorrelation *po*(*t*0)*po*(*t*) = *p*∞*po*|*o*(*t*). The autocovariance of a single channel *Co*(*t*) open probability can then be written as

$$\begin{aligned} C\_o(t) &= \langle \left[ p\_o(t\_0) - p\_\infty \right] \left[ p\_o(t) - p\_\infty \right] \rangle \\ &= p\_\infty p\_{o|o}(t) - p\_\infty^2 \\ &= p\_\infty (1 - p\_\infty) e^{-t/\tau} \end{aligned}$$

The autocovariance of the current through a population of *N* such channels, *CI*(*t*), is simply given by

$$C\_I(t) = \text{Ni}^2 C\_o(t)$$

where *i* is the single-channel current. Note that the autocovariance function at *t* = 0 is equal to the variance, *CI*(0) = *Ni*2*p*∞(1 − *p*∞) = σ2, and decays exponentially with time constant τ , so that when *t* τ , *CI*(*t*) → 0.

The Wiener-Khinchin theorem states that the power spectrum is equal to the real part of the Fourier transform of the autocovariance function

$$S\_I(\omega) = 4 \int\_0^\infty \mathcal{C}\_I(t) e^{-2\pi i \theta t} dt$$

$$= S\_I(0) \frac{1}{1 + \left(2\pi f \tau\right)^2} \tag{4}$$

where *SI*(0) = 4*Ni*2*p*∞(1 − *p*∞)τ and the pre-factor 4 arises from our definition of spectral density. Thus for the two-state channel population, the power spectrum is a single Lorenztian function with a corner frequency *fc* = 1/(2πτ ). Above this corner frequency, the power of the signal falls off ∝ 1/*f* 2.

In an analogous way we can calculate the power spectra of the HH Na+ and K+ channel populations. For the HH K+ channel composed of four identical and independent sub-units, the conditional probability that the channel is open at time *t* given that it was open at time 0 is

$$p\_{K,o|o}(t) = \begin{pmatrix} n\_{\infty} + (1 \ - \ n\_{\infty})e^{-t/\tau\_n} \end{pmatrix}^4 \tag{5}$$

Hence, the autocovariance of the current noise from the HH K+ channel population is a sum of exponentials,

$$\begin{split} C\_{IK}(t) &= N\_K i\_K^2 \left( n\_{\infty}^4 p\_{K, o \mid o}(t) - \left( n\_{\infty}^4 \right)^2 \right) \\ &= N\_K i\_K^2 n\_{\infty}^4 \left[ (1 - n\_{\infty})^4 e^{-4t/\tau\_{\pi}} + 4n\_{\infty} (1 - n\_{\infty})^3 e^{-3t/\tau\_{\pi}} \right] \\ &+ 6n\_{\infty}^2 (1 - n\_{\infty})^2 e^{-2t/\tau\_{\pi}} + 4n\_{\infty}^3 (1 - n\_{\infty}) e^{-t/\tau\_{\pi}} \right] \end{split}$$

The corresponding power spectrum of the current noise is

$$\begin{split} S\_{IK}(\omega) &= 4N\_k n\_{\infty}^4 i\_K^2 \tau\_n \left[ (1 - n\_{\infty})^4 \frac{4}{16 + \omega^2 \tau\_n^2} \right. \\ &\left. + n\_{\infty} (1 - n\_{\infty})^3 \frac{12}{9 + \omega^2 \tau\_n^2} + n\_{\infty}^2 (1 - n\_{\infty})^2 \frac{12}{4 + \omega^2 \tau\_n^2} \right. \end{split}$$

$$\left[ + n\_{\infty}^{3} (1 - n\_{\infty}) \frac{4}{1 + \omega^{2} \tau\_{n}^{2}} \right]$$

This is the sum of four Lorenztians with corner frequencies equal to 4/(2πτ*n*), 3/(2πτ*n*), 2/(2πτ*n*), and 1/(2πτ*n*). Because at the resting potential *n*∞ is close to zero, the first term in the square brackets with correlation time τ*n*/4 will dominate the power spectrum.

Similarly, for the Na+ current noise power spectrum one has

$$p\_{\rm Na, o|o} = \left(m\_{\infty} + (1 - m\_{\infty})e^{-t/\tau\_m}\right)^3 \left(h\_{\infty} + (1 - h\_{\infty})e^{-t/\tau\_h}\right)^3$$

and the Na+ current noise autocovariance is

$$\begin{split} C\_{\text{Na}}(t) &= N\_{\text{Na}}i\_{\text{Na}}^{2} \left( m\_{\infty}^{3}h\_{\infty}p\_{\text{Na},o|o} - \left( m\_{\infty}^{3}h\_{\infty} \right)^{2} \right) \\ &= N\_{\text{Na}}i\_{\text{Na}}^{2}m\_{\infty}^{3}h\_{\infty} \left[ 3m\_{\infty}^{2}h\_{\infty}(1 - m\_{\infty})e^{-t/\tau\_{m}} \right. \\ &\left. \left. + 3m\_{\infty}h\_{\infty}(1 - m\_{\infty})^{2}e^{-2t/\tau\_{m}} \right. \right. \\ &\left. + h\_{\infty}(1 - m\_{\infty})^{3}e^{-3t/\tau\_{m}} \right. \\ &\left. + m\_{\infty}^{3}(1 - h\_{\infty})e^{-t/\tau\_{h}} \right. \\ &\left. + 3m\_{\infty}^{2}(1 - m\_{\infty})(1 - h\_{\infty})e^{-t/\tau\_{m} - t/\tau\_{h}} \right. \\ &\left. + 3m\_{\infty}(1 - m\_{\infty})^{2}(1 - h\_{\infty})e^{-2t/\tau\_{m} - t/\tau\_{h}} \right. \\ &\left. + (1 - m\_{\infty})^{3}(1 - h\_{\infty})e^{-3t/\tau\_{m} - t/\tau\_{h}} \right] \end{split}$$

The corresponding power spectrum of the Na+ current noise is

*SINa*(ω) = 4*NNai* 2 *Na*(*m*3 ∞*h*∞) 2 1 − *m*∞ *m*∞  3τ*m* 1 + ω2τ 2 *m* + 1 − *m*∞ *m*∞ 2 6τ*m* 4 + ω2τ 2 *m* + 1 − *m*∞ *m*∞ 3 3τ*m* 9 + ω2τ 2 *m* + 1 − *h*∞ *h*∞  τ*h* 1 + ω2τ 2 *h* + 3 1 − *h*∞ *h*∞  1 − *m*∞ *m*∞  τ*m*τ*h* τ*m* + τ*h* 1 1 + (ωτ*m*τ*h*/(τ*m* + τ*h*)) 2 + 3 1 − *h*∞ *h*∞  1 − *m*∞ *m*∞ 2 τ*m*τ*h* τ*m* + 2τ*h* 1 1 + (ωτ*m*τ*h*/(τ*m* + 2τ*h*)) 2 + 1 − *h*∞ *h*∞  1 − *m*∞ *m*∞ 3 τ*m*τ*h* τ*m* + 3τ*h* 1 1 + (ωτ*m*τ*h*/(τ*m* + 3τ*h*)) 2

Near rest *m*∞ 1 and *h*∞ ≈ 0.6, so that the third and last Lorentzians dominate the powerspectrum. The corresponding corner-frequencies are 3/(2πτ*m*) and (τ*m* + 3τ*h*)/(2πτ*m*τ*h*), which are virtually identical.

We calculate the quasi-active (linearized) membrane impedance using standard methods (Mauro et al., 1970; Koch, 1999).

## **3. RESULTS**

## **3.1. STOCHASTIC POTASSIUM CHANNELS CAN TRIGGER SPONTANEOUS ACTION POTENTIALS**

Previously, it has been demonstrated that a Hodgkin–Huxley (HH) type neural model with discrete Markovian stochastic ion channels instead of the classic continuous deterministic rate equations can fire spontaneous action potentials if the membrane patch is small (Skaugen and Walløe, 1979; Clay and DeFelice, 1983; Strassberg and DeFelice, 1993; Chow and White, 1996; Schneidman et al., 1998). However, the relative contributions of the Na+ and K+ channel populations to spontaneous activity are less well understood. To investigate, we simulate the HH squid axon model using the PSICS simulator (Cannon et al., 2010) with stochastic Markovian ion channels while varying the membrane patch area under three different conditions: first, both sodium (Na+) and potassium (K+) channels stochastic ("all stochastic"), second, Na+ stochastic but K+ deterministic, and third, Na+ deterministic but K+ stochastic. Comparing the spontaneous firing rate between the three conditions allows us to find whether Na+ or K+ channels contribute most to spontaneous activity.

As observed previously (Chow and White, 1996), if a fixed density of ion channels is assumed, then the firing rate decreases approximately exponentially with increasing membrane patch area (**Figure 1**) such that membrane areas greater than ∼400µm2 produce almost no spontaneous action potentials approximating the deterministic model. This exponential dependence of spontaneous rate with membrane area is consistent with a stochastic barrier-escape problem (Chow and White, 1996). As the channels are independent, the voltage variance is proportional to the number of channels *N*, but inversely proportional to the square of the area *A*, since the input impedance decreases linearly with area. Combining these two opposing factors, the spontaneous rate scales as exp ( − *A*2/*N*) ∝ exp ( − *A*).

When either Na+ or K+ channels are switched to deterministic mode, spontaneous firing rate is reduced compared to the fully stochastic mode. Surprisingly, however, stochastic K+ channel gating alone triggers greater spontaneous firing rates than stochastic Na+-channel gating alone (**Figure 1B**). At first impression, this result might be counter-intuitive because the opening of Na+ channels is necessary for the initiation of an action potential, while the much slower K+ channels are conventionally considered responsible for the re-polarizing phase. A simple conceptual model for spontaneous spike generation might therefore be that the chance opening of a few Na+ channels depolarizes the membrane and activates the runaway Na+ channel opening underlying the action potential. However, stochastic *closure* of K+ channels can also depolarize the membrane, similarly activating Na+ channels to trigger an action potential. We test this possibility by examining the dynamics of Na+ and K+ currents. We adapt the 'spike-triggered average' (STA) measure from the neural coding literature. Here we determine the average total current of a given ion channel population *x* at time interval *t* prior to a spontaneous action potential at time *ti*, averaged over *n* such events,

$$I\_{\mathbf{x}}^{\mathrm{STA}}(t) = \frac{1}{n} \left[ \sum\_{i=1}^{n} I\_{\mathbf{x}}(t\_i - t) \right]$$

In the "all stochastic" mode (**Figure 2A**, solid curves), we find that the STA potassium current *ISTA K* (*t*), drops between 8 and 2

ms before the spike, while there is a simultaneous increase in the STA sodium current *ISTA Na* (*t*). Nearer to the spike the Na+ and K+ currents grow rapidly but in opposite directions as the action potential forms. A positive current corresponds to depolarization. This depolarizing action of the K+ current change be clearly seen in **Figure 2B** where we plot *IK*(*t*) and *INa*(*t*), the change in Na+ and K+ relative to resting current. Importantly, the change in *ISTA K* precedes the increase in *ISTA Na* (**Figure 2B**), suggesting that spontaneous action potential firing in this model is primarily driven by K+ channel fluctuations, not Na+ noise. We test this explanation by simulating K+ channel conductance in deterministic mode and repeating the STA measurement. As expected, in this case spontaneous spikes are not preceded by a drop in K+ conductance, but instead driven by an elevated Na+ conductance fluctuation (**Figure 2A**, dotted curves).

Further examination of the Na+ and K+ current dynamics confirm these findings. In **Figure 2C** we plot the timecourse of *INa*(*t*) (y-axis) vs. *IK*(*t*) (x-axis) from 10 to 1 ms prior to each recorded action potential together on the same plot. Regions to the lower-right of the identity diagonal (dashed black line) indicate timepoints where K+ current fluctuations are contributing more to voltage depolarization than Na+ current fluctuations. The mean STA curve for the "all stochastic" model (solid curve) initially takes off into this right-hand region. In contrast, the Na+ only stochastic (dotted curve) moves only slightly to the right of the origin before taking off in the vertical (Na+-driven) direction.

**Figures 2A–C** plot the average behavior. As the system is stochastic, we expect that occasionally Na+ fluctuations trigger spikes as well. To get a sense of the spike-to-spike variability in **Figure 2D** we plot the percentage of cases where the potassium current exceeds the sodium current, *IK*(*t*) > *INa*(*t*), as a function of time before spike. In general, this quantity is timedependent because the depolarization due to K+ fluctuations will most likely be maximal at some time between 8 and 2 ms before the spike. Once the full spike upswing begins (∼1 ms before *t* = 0) Na+ always dominates. In the all-stochastic case, the earliest phase of most action potentials are K+-driven. In contrast, in the Na+ stochastic case (dotted line) the majority of spikes are driven primarily by Na+ fluctuations.

A slightly different picture appears for the K+ stochastic simulations. In this mode, Na+ fluctuations are removed so all spikes must by initially triggered by K+ fluctuations. Consequently, the spike-triggered average K+ current fluctuation amplitude is even greater than in the all-stochastic model (**Figures 2A–C**) and an even larger percentage of spikes are driven by momentary fluctuations in K+ currents (**Figure 2D**). In summary, these simulations show that K+ channel noise is the dominant driver of spontaneous spiking in the stochastic Hodgkin–Huxley model.

Schneidman et al. (1998) looked at stochastic Na+ and K+ channel trajectories during spike initiation to address whether a sufficiently strong stimulus can override the intrinsic channel noise. However, in contrast to our STA analysis for determining the contributions of Na+ vs. K+ noise to spontaneous spiking, they examined stimulus-driven firing by injecting a fluctuating current into the model neuron. In this stimulus-driven case, the observed trajectories of the Na+ vs. K+ currents combine the effect of the stimulus current dynamics and the effect of the channel noise. Our analysis however shows that the drop in noisy K+ current occurs naturally before the spike.

#### **3.2. THE FACTORS DETERMINING A CONDUCTANCE'S CONTRIBUTION TO MEMBRANE NOISE**

These results lead to the questions: What properties of the HH K+ conductance cause it to trigger more spontaneous action potentials than the Na+ conductance? And how can the contribution of an arbitrary channel type be estimated? We proceed by first calculating the resulting noise in the membrane voltage. In the limit of large areas (small noise), this can be calculated exactly. In the subsequent section we relate the voltage noise to spontaneous spike rates. This will turn out to be only approximately possible.

There are at least five possible factors that determine a channel population's contribution to membrane noise:

## *3.2.1. Open probability (11* **×** *noisier for K***+***)*

First, Na+ and K+ have different steady-state open probabilities at resting membrane potential. The steady-state probability of a single ion channel being open is identical to the steady-state permeability fraction of the corresponding macroscopic conductance in the classic HH formalism. The steady-state Na+ conductance is obtained from the product of the steady-state values of the *m*∞ and *h*∞ gating variables:

$$\operatorname{g\_{Na\infty}}(V) = \bar{\operatorname{g}}\_{\mathrm{Na}}\left[m\_{\infty}(V)\right]^{\mathfrak{F}}h\_{\infty}(V)$$

where *g*¯*Na* is the maximal conductance through the Na+ channel population. Hence the open probability *po* = (*m*∞)3*h*∞. The gating variables *m*∞ and *h*∞ can in turn be expressed in terms of the forward and backward gating rates α*m* and β*m*, and α*h* and β*h*, see Equation 2.The steady-state K+ open probability equals *po* = [*n*∞(*V*)] 4. At the resting potential of −65 mV in the HH squid axon model, the steady-state open probabilities are ∼0.000089 for the Na+ and and ∼0.010 for the K+ channels. At any instant the open probability follows a binomial distribution so that the variance in the single channel current σ2 *i* = *i* 2*po*(1 − *po*), where *i* is the single-channel current. The variance is parabolic in *po*: zero when *po* = 0 or 1, and maximal at *po* = 0.5. As below spike threshold, most ion channels have very low open probabilities, the standard deviation can be approximated by

$$
\sigma\_i = i\sqrt{p\_o}
$$

Therefore, ion channel populations with greater *po* at resting membrane potential tend to have larger current fluctuations than populations with lesser *po*. This effect predicts that the standard deviation of the K+ channel noise is 10.7 times larger than Na+ channel noise.

## *3.2.2. Number of channels (1.8* **×** *noisier for Na***+***)*

Second, because the channels act independently, the standard deviation of the number of open channels grows proportional to √*N*. Thus channel populations with greater *N* have greater fluctuations in their absolute number of open channels. The Na+ population has 3.33× more channels in the standard HH model than the K+ population (**Table 1**), yielding Na standard deviation larger by a factor √10/3 ≈ 1.8.

## *3.2.3. Driving force (10* **×** *noisier for Na***+***)*

The third factor is the difference in driving force for each conductance. As the HH model assumes that these ion channel current-conductance relationships are Ohmic (linear), the current through an open channel is proportional to the difference between the membrane potential and the channel's driving force,

$$i\_{\mathbf{x}} = \gamma\_{\mathbf{x}} (V\_{\mathbf{m}} - E\_{\mathbf{x}})$$

where *ix* is the single-channel current, γ*x* is the single-channel conductance, *Vm* is the membrane potential and *Ex* is the conductance's reversal potential, given by the Nernst equation. In the HH model, *ENa* = +50 mV, *EK* = −77 mV, and *Vrest* = −65 mV, giving Na+ a driving force of +115 mV and K+ a driving force of −12 mV. This means that the driving force for the Na+ current is 9.6× greater than the K+ current at *Vrest*.

## *3.2.4. Single-channel conductance (identical for Na***+** *and K***+***)*

Fourth, the single-channel conductance γ*x* is another important factor determining a channel's contribution to membrane noise. For the same channel population current per unit squared cell membrane, a larger γ*x* implies smaller *N*, and a larger *ix*. Hence channels with a larger single-channel conductance will have greater current fluctuations. In our implementation of the stochastic HH model, however, we assume, like most other stochastic models, the same single-channel conductance for both Na+ and K+ (20 pS). The value of 20 pS is close to experimentally measured estimates for Na+ (14 pS) (Bezanilla, 1987), while the K+ conductance in the squid axon is probably made up of multiple different channel types, with single-channel conductances estimated at 10, 20, and 40 pS (Llano et al., 1988).

These four factors can be put together to construct a binomial model of the amplitude of channel noise at steady state. This model does not have any notion of dynamics or channel kinetics. We calculate the steady-state open probabilities at rest directly from the Hodgkin–Huxley equations, and test the binomial model's ability to reproduce simulated voltage-clamp data and probe its predictions on the relative magnitudes of Na+ and K+ channel noise (**Figure 3**).

We find that, as expected, the binomial model exactly predicts the conductance and current fluctuations from voltage-clamp simulation data at resting potential of −65 mV (**Figures 3C,F**). We use the binomial model to estimate the steady-state standard deviation in open channel numbers and total current from the Na+ and K+ populations at a range of membrane potentials (**Figures 3A–B,D–E**). As expected from the above analysis the Na+ current standard deviation is about 1.7× that of the K+ current (matched in simulations, see circle symbols in **Figure 3E**).

## *3.2.5. Channel gating kinetics and membrane filtering (3* **×** *noisier for K***+***)*

The fifth factor is that the Na+ and K+ conductances have different gating kinetics. These differences are important because the current fluctuations from each ion channel populations are filtered differentially by the membrane impedance, hence altering each channel's contributions to membrane voltage noise.

In the Methods we calculate the power-spectra of the Na+ and K+ current noise assuming a constant membrane potential and small noise. Both powerspectra are sums of multiple Lorentzians,

$$S\_I(f) = \sum\_k \frac{a\_k}{1 + \left(f/f\_c^k\right)^2}$$

where the *f k c* are the corner frequencies of the Lorentzians (the frequency at which the powerspectrum is half of the zero frequency magnitude), and *ak* are (voltage-dependent) coefficients. The full expressions are given in the Methods, but the K+ noise spectrum is dominated by a Lorentzian with corner frequency *fc* = 4/(2πτ*n*). At the resting potential τ*n* ∼ 5.5 ms (Equation 3), so that the dominant Lorentzian has a corner frequency of 115 Hz. The Na+ spectrum is dominated by a Lorentzian with *fc* = 3/(2πτ*m*). At rest τ*m* ∼ 0.24 ms, τ*h* ∼ 8.5 ms, so that the

dominant Lorentzian for Na+ has a much higher corner frequency of 1980 Hz—note however that we include all Lorentzian terms in the presented results, not just the dominant one. The analytically calculated spectra match well the estimated spectra of the simulated stochastic currents, **Figure 4A**.

To calculate the voltage response to the current noise, we approximate the active membrane by a linear impedance. The voltage noise power spectrum *SV*(*f*) is given by the generalized Ohm's law

$$S\_V(f) = S\_I(f)|Z(f)|^2$$

where *SI*(*f*) is the power spectrum of a current noise source and *Z*(*f*) is the membrane impedance. The impedance relates

colored curves are theoretical, derived from channel kinetic schemes. **(B)** Variance of Na+ and K+ current noise from both simulation and theory. Note Na+ channel noise variance is greater than K+ channel noise. **(C)** Total HH membrane impedance at −65 mV as a function of signal frequency. **(D)** Voltage noise of Na+ and K+ channels calculated from current noise spectra and membrane impedance. Gray curve is sum of Na+ and K+ noise, while

K+ channels contribute more to voltage noise than Na+ channels. **(F)** Comparison of voltage noise variance from theory and simulation as a function of membrane area. Note curves substantially diverge only for areas < 500µm. For those small areas spontaneous spiking occurs and the associated large voltage fluctuations are not part of the theory. All other panels use a membrane area of 1000µm2.

changes in voltage to changes in currents. In active membranes the input impedance is not given by just the capacitance and leak conductance, but also by any other channels open at rest and their reaction to small changes in the voltage. The impedance can be found by linearizing the four-dimensional (*V*, *m*, *h*, and *n*) HH equations around the resting state (e.g., Mauro et al., 1970; Koch, 1984; Carnevale and Hines, 2006). In the HH model the presence of Na+ and K+ conductances introduce a resonance in the impedance at ∼100 Hz, but the 1/*f* behavior still dominates at higher frequencies (**Figure 4C**). The resonant peak in the impedance is the electrical signature of an inductor. Although such an inductor has no physical counterpart in the biological cell membrane, the delayed-rectifier K+ conductance opposes changes in membrane potential and for small currents behaves as a phenomenological inductance (Mauro et al., 1970; Koch, 1984, 1999).

Now we combine the membrane impedance with the current noise spectra to calculate each channel population's contribution to voltage noise. In **Figure 4D** we plot the theoretical power spectra of the voltage noise from the HH Na+ and K+ channel populations, calculated at *Vm* = *Vrest* = −65 mV. The sum of the Na+ and K+ power spectra give the total voltage noise power spectrum (gray line in **Figure 4D**). This predicts almost exactly the power spectrum measured from simulation (dashed line in **Figure 4D**). The voltage noise variance from each channel population equals the integrated power spectrum:

$$
\sigma\_\mathbf{x}^2 = \int\_0^\infty \mathcal{S}\_\mathbf{x}(f) df
$$

where subscript *x* indicates the relevant channel population. These are graphed in **Figure 4E**. It is clear the K+ channel fluctuations contributes ∼ 4× more *voltage* noise variance than Na+. While the Na+ *current* noise has a greater variance than the K+ current noise, it is filtered more strongly by the membrane impedance.

Can we summarize the total effects of membrane filtering on noise from the two channel types? One way to quantify this effect is to take the ratio of voltage and current noise standard deviations for each channel type, *rx* = σ*Vx*/σ*Ix*. Doing so we find that *rNa* = 44.5 M and *rK* = 141.7 M, implying that, after all other factors are accounted for, membrane filtering attenuates Na+ noise ∼ 3× more than K+ noise.

As both the noise spectra and the impedance are voltagedependent, these calculations assume that the voltage remains at a fixed potential, which holds if the fluctuations are small. Because here we simulate a large membrane area (1000µm2) with low membrane resistance, voltage changes are small and there is only a small discrepancy between the voltage noise calculated analytically and the estimate from simulation (**Figure 4E**). In particular for small membrane areas the voltage fluctuations will be large, nevertheless the approximation remains good down until areas where spontaneous spikes appear, **Figure 4F**. At this point, currents associated to the spike will dominate the measured current.

In summary, the contribution of each channel type to membrane noise is determined by their number, single-channel conductance, voltage dependencies and gating kinetics. This is true for any neuron model. In the case of the HH model, the sum properties of K+ channels at subthreshold voltages make their contributions to membrane noise greater than that from Na+ channels.

#### **3.3. EFFECT ON SPONTANEOUS SPIKE RATES**

So far we have seen that the channels in the HH model contribute differentially to the noise in the membrane voltage and that K+ channels have the largest contribution. One would expect that therefore K+ channels are the most important contributor of noise-driven spontaneous spike activity as well. This is indeed the case as we have seen in **Figure 1** but quantitatively the link between the spontaneous rate and the subthreshold membrane voltage fluctuations is not trivial.

The analysis of spontaneous spiking rate in the HH model to correlated noise is a complicated stochastic differential equation problem. Reducing the spiking mechanism to a one dimensional escape problem, Chow and White (1996) derived the spontaneous rate using multiplicative, white (uncorrelated) current noise to approximate the Na+ noise. But it is not obvious how such an analysis can be extended to colored (correlated) noise. The timeconstants of the K+ noise, Na+ noise and the membrane are all of similar magnitude (**Figure 4**), complicating any perturbative expansion. Even in much simpler integrate-and-fire neuron models, the treatment of correlated noise is complicated, resulting in a two-dimensional Fokker-Planck equation that can only be solved in certain limits (Brunel and Sergi, 1998; Moreno-Bote and Parga, 2004; Alijani and Richardson, 2011).

We first examine how the correlation time in the Gaussian noise model affects spontaneous firing rates in a HH neuron. Traditionally, studies have kept the variance of the injected noise fixed while varying the correlation time. However, as shown above, this can lead to widely different voltage fluctuations due to the differential membrane filtering. The idea we examine here is that instead of the current variance, the voltage variance is a better predictor of the spontaneous rate. This was recently shown in exponential integrate-and-fire neurons: the noise driven firing rate was relatively independent of the noise correlation time when the voltage variance was kept constant (Alijani and Richardson, 2011).

In order to research this in the HH model, we simulated a deterministic HH neuron and injected a Gaussian noise current with a correlation time τ and with a correlation function

$$
\langle I(T)I(T+t) \rangle = \frac{\sigma\_V^2}{z^2(\pi)} e^{-|t|/\pi}
$$

The function *z*(τ ) is an impedance that relates the voltage variance to the current variance of injected colored noise with timeconstant τ . It is given by *z*2(τ ) = *SI*(*f*)|*Z*(*f*)| 2*df* , where *SI*(*f*) = 4τ 1 + (2π*f* τ )2 is the power-spectrum of the noise current, and *Z*(*f*) the linearized impedance of the HH model. Its shape reflects the resonance in the impedance, **Figure 5A**. As a result of this impedance correction the membrane voltage variance in the limit of small fluctuations equals σ2 *V* and is thus independent of τ .

When the amplitude of the noise is scaled up, spontaneous spikes appear. The resulting spontaneous rate is shown in **Figures 5B** (black solid curve),**C**. Remarkably, across six orders of magnitude in the noise time-scale, ranging from noise much faster than any channel or membrane to extremely slow noise, the spontaneous rate varies only one order of magnitude (from 40 to 5 Hz for the noise in **Figure 5B**). With less noise this ratio can be larger as slow noise might then be unable to evoke the spontaneous spikes, **Figure 5C**. Furthermore, the relation is monotonic, which is useful when comparing two noise sources. The shape

**FIGURE 5 | Spontanous action potentials in the HH model driven with Gaussian noise with varying amplitude and correlation time. (A)** The impedance factor *z*(τ ) gives the voltage fluctuations resulting from a correlated noise current as a function of the noise correlation time. It is used to ensure the voltage variance is identical as the time constant of the noise is varied. Membrane area of 1000µm2. **(B)** Spontaneous activity as a result of a Gaussian current noise, as a function of the correlation time of the noise. When the variance of the current is fixed, the resulting rates vary strongly depending on the correlation time (dashed curve). However, when the noise is calibrated to yield the same variance in the

membrane voltage irrespective of the noise correlation time, the firing rate is much less variable even across six orders of magnitude (solid curve). **(C)** The firing rate vs. the impedance corrected current noise (expressed in σ*V* ) for various values of the noise time-constant. **(D)** The spontaneous firing rate vs. membrane area for K+ (red) and Na+ (blue) noise and their various approximations. Solid curve: full simulation, redrawn from **Figure 1B**. Thick dashed curve: approximation using Gaussian noise with identical variance and time-constant. Thin dashed curve: binomial approximation of the noise. In particular for Na+ noise the approximations yield rates that are substantially too low.

of the curve highlights that fast fluctuating noise is typically an effective driver of neurons, while slow varying noise will tend to inactivate the Na channel close to threshold and is less effective.

Instead, the *current* variance is a much worse predictor of the firing rate (dashed line). Note that in the limit of very slow noise, the dynamics decouple and the spontaneous rate equals *P*(*I*)*f*(*I*)*dI*, where *P*(*I*) is the distribution of currents and *f*(*I*) is the neuron's deterministic f-I curve.

Finally, we ask if we can use these results to estimate the spontaneous rates caused by K+, Na+, or in fact any arbitrary channel noise. In the previous section the noise currents were approximated by colored Gaussian noise (an OU process) with variance and correlation time derived from the channel kinetics at the resting voltage, and filtered by the membrane linearized around rest. These approximations hold very well in the subthreshold regime, i.e., for small noise—equivalent to large membrane areas — see **Figure 4E**, but it is a priori unclear whether they also hold for larger noise amplitudes when spontaneous spikes appear.

We used a colored Gaussian current noise to model the K+ and Na+ noise, and injected this into a deterministic HH model. For example, in case of the K+ channel the variance according to the above sections is

$$
\langle \delta I^2 \rangle = \rho\_K A [(V\_{\text{rest}} - E\_K) \gamma\_K n\_{\infty}^4 (V\_{\text{rest}})]^2
$$

and its correlation time is τ*n*/4. The spontaneous rate of the neuron driven by this noise was examined as a function of the membrane area. Although this could be accidental, for K+ noise, the noise model gave an almost perfect fit to the fully stochastic simulations, **Figure 5D** compare solid and thick dashed curves. However, the approximated Na noise lead to far too few spontaneous spikes, **Figure 5D**. Its standard deviation had to be increased by some 50% to fit the simulated spontaneous rates. This need for a fudge factor shows that for smaller areas the Na noise is not well described by additive, colored Gaussian noise. There are many possible cause for this mismatch: the binomial instead of Gaussian current distribution (which additional simulations showed to be a small effect, **Figure 5D** thin dashed curve), the voltage dependence of the impedance, the dynamics of the full HH system, and likely most important, the strong dependence on the noise current on the membrane voltage, **Figure 3E**.

In summary then, while we find that sub-threshold noise can be estimated accurately, caution is needed when extrapolating to the spiking regime.

#### **3.4. APPLICATION TO CA1 PYRAMIDAL NEURON MODEL**

Above we showed how to break down the factors determining a given ion channel type's contribution to voltage noise and spontaneous spiking, using the Hodgkin–Huxley squid axon model as an example. However, our approach is completely general and could in principle be applied to any neuron model. To demonstrate its straightforward application, we now perform the same analysis on a well-studied mammalian cell type: the rodent hippocampal CA1 pyramidal neuron. We studied a single-compartment model of this cell type using a well-validated model from the literature (Migliore et al., 1999; Jarsky et al., 2005). The original model contained three active channel types: an Na+ channel, an A-type K+ channel, Ka, and a delayed rectifier K+ channel, Kdr. We built stochastic versions of these active conductances with parameters exactly as previously used (Jarsky et al., 2005).