## Abstract

Interspike interval (ISI) distributions of cortical neurons exhibit a range of different shapes. Wide ISI distributions are believed to stem from a balance of excitatory and inhibitory inputs that leads to a strongly fluctuating total drive. An important question is whether the full range of experimentally observed ISI distributions can be reproduced by modulating this balance. To address this issue, we investigate the shape of the ISI distributions of spiking neuron models receiving fluctuating inputs. Using analytical tools to describe the ISI distribution of a leaky integrate-and-fire (LIF) neuron, we identify three key features: *1*) the ISI distribution displays an exponential decay at long ISIs independently of the strength of the fluctuating input; *2*) as the amplitude of the input fluctuations is increased, the ISI distribution evolves progressively between three types, a narrow distribution (suprathreshold input), an exponential with an effective refractory period (subthreshold but suprareset input), and a bursting exponential (subreset input); *3*) the shape of the ISI distribution is approximately independent of the mean ISI and determined only by the coefficient of variation. Numerical simulations show that these features are not specific to the LIF model but are also present in the ISI distributions of the exponential integrate-and-fire model and a Hodgkin-Huxley-like model. Moreover, we observe that for a fixed mean and coefficient of variation of ISIs, the full ISI distributions of the three models are nearly identical. We conclude that the ISI distributions of spiking neurons in the presence of fluctuating inputs are well described by gamma distributions.

- excitation-inhibition balance
- integrate-and-fire model
- Hodgkin-Huxley-like model

neurons in vivo emit action potentials in an irregular manner (Softky and Koch 1993; Bair et al. 1994; Holt et al. 1996; Shinomoto et al. 1999; Compte et al. 2003; Maimon and Assad 2009). In contrast, neuronal activity recorded in cortical slices in vitro is regular and highly reproducible, suggesting that neuronal firing is inherently deterministic rather than stochastic (Mainen and Sejnowski 1995; Holt et al. 1996). Explaining the origin of irregular activity in vivo has been the focus of a number of studies in the last two decades, and theoretical (Shadlen and Newsome 1994, 1998; Tsodyks and Sejnowski 1995; Vreeswijk and Sompolinsky 1996; Amit and Brunel 1997; Troyer and Miller 1997) as well as experimental results (Destexhe et al. 2001, 2003; Shu et al. 2003; Haider et al. 2006) suggest that irregular firing arises from a highly fluctuating drive generated by a balance between excitatory and inhibitory synaptic inputs to the neurons.

Balanced synaptic inputs are a general mechanism to generate irregular firing, yet it is debated to what extent this mechanism reproduces the details of the statistics of cortical spike trains (Shinomoto et al. 1999; Sakai et al. 1999). Irregular firing in vivo is often thought of as a Poisson process, but a close examination of experimental recordings reveals a large variety of interspike interval (ISI) statistics. Here we focus on the shapes of ISI distributions, which have been found to range from narrow to bursting (Bair et al. 1994; Compte et al. 2003; Maimon and Assad 2009). It is not clear that fluctuating inputs to the neurons are sufficient to account for the full range of observed ISI distributions. The amount of fluctuations in the input can be modulated by the details of the balance between excitation and inhibition, and a range of ISI distributions can be produced in this manner. However, reproducing a particular shape might require specific neuronal properties in addition to the fluctuating input. Understanding what range of ISI distributions is spanned when input fluctuations are varied is therefore an important issue when interpreting experimental data.

In this study, we systematically examine the shape of ISI distributions generated by spiking neuron models receiving fluctuating inputs. Using analytical tools for the leaky integrate-and-fire (LIF) model, we first show that for long ISIs the distribution decays exponentially independently of the parameters of the input. As the amplitude of input fluctuations is increased at a given output firing rate, the corresponding exponential tail describes a larger and larger portion of the ISI distribution, up to a point where bursting appears and progressively takes over the ISI distribution. At a fixed firing rate, the ISI distribution of a LIF neuron progressively evolves between three regimes corresponding to an increasing value of the ISI coefficient of variation (CV) : *1*) a narrow distribution for low input fluctuations; *2*) an exponential distribution with an effective refractory period at intermediate input fluctuations; and *3*) an exponential distribution with bursts for high input fluctuations. In contrast, varying the mean ISI only marginally affects the shape of the ISI distribution, which thus approximately depends only on the CV. In view of these properties, we argue that the range of ISI distributions spanned by varying the fluctuating input can be described to a first approximation by the family of gamma distributions. Numerical simulations show that these findings hold also for the exponential integrate-and-fire (EIF) model and a Hodgkin-Huxley-type conductance-based model; the full ISI distributions of the three models lie very close to each other if the first two moments are fixed.

## MATERIALS AND METHODS

### Integrate-and-Fire Models

Integrate-and-fire models are simplified models of neurons in which action potentials (APs) are generated solely from the underlying dynamics of the membrane potential (Gerstner and Kistler 2002). We used generalized integrate-and-fire models, with dynamics given by (Fourcaud-Trocmé et al. 2003):
*V* is determined with respect to the resting potential of the cell, *c _{m}* = 1 μF/cm

^{2}is the membrane capacitance,

*g*= 50 μS/cm

_{L}^{2}is the leak conductance, Ψ(

*V*) is a spike generating current, and

*I*is the total current elicited by synaptic inputs to the neuron.

We studied two different versions of the integrate-and-fire model.

#### LIF model.

In this model, Ψ(*V*) = 0, there is no spike-generation current and an AP is emitted when the membrane potential crosses a fixed threshold value *V _{th}*. The membrane potential is subsequently reset to a value

*V*after a refractory period τ

_{r}*. The values used in the simulations were*

_{rp}*V*= 20 mV,

_{th}*V*= 10 mV, and τ

_{r}*= 2 ms.*

_{rp}#### EIF model.

In this model, the spike generation current is exponential:
*V _{th}*, it diverges to infinity in finite time. This divergence represents the firing of an AP. Following the divergence, the membrane potential is reset to a value

*V*after a refractory time τ

_{r}*. The parameter Δ*

_{rp}*quantifies the sharpness of the AP initiation. The parameter values used in most of this study were Δ*

_{T}*= 1 mV [a typical value for pyramidal cells (Badel et al. 2008)],*

_{T}*V*= 10 mV,

_{th}*V*= 3 mV and τ

_{r}*= 2 ms.*

_{rp}### Conductance-Based Model

The Wang-Buzsáki (WB) model is a modified Hodgkin-Huxley model (Wang and Buzsáki 1996). The dynamics of the membrane potential are given by
*c _{m}* is the membrane capacitance (

*c*= 1 μF/cm

_{m}^{2}),

*I*=

_{L}*g*(

_{L}*V*−

*V*) is the leak current (

_{L}*g*= 0.1 mS/cm

_{L}^{2};

*V*= −65 mV),

_{L}*I*=

_{Na}*g*

_{Na}m^{3}

*h*(

*V*−

*VN*) is the sodium current,

_{a}*I*=

_{K}*g*

_{K}n^{4}(

*V*−

*V*) is the delayed rectifier potassium current, and

_{K}*I*(

*t*) is the total synaptic input current.

The activation of the sodium current is assumed to be instantaneous:
*h* and *n* are given by:
*x* = *h*, *n*.

The functions *α _{x}* and

*β*are given by:

_{x}The maximum conductance densities and reversal potentials are: *g _{Na}* = 35 mS/cm

^{2},

*V*= 55 mV,

_{Na}*g*= 9 mS/cm

_{K}^{2}, and

*V*= −90 mV.

_{K}### Current-Based Inputs to the Neurons

In most of this study, we used current-based inputs to the neurons. We assumed that the inputs to the neurons consist of a sum of large number of synaptic inputs, each individual synaptic input being of small amplitude. We therefore used the diffusion approximation (Tuckwell 1988) and represented the total input as a Gaussian random processes parametrized by its mean μ and SD σ, so that in *Eq. 1*
* _{m}* =

*c*/

_{m}*g*is the membrane time constant (τ

_{l}*= 10 ms) and η is a gaussian process of zero mean, unit variance, and correlation time τ*

_{m}*. Note that both μ and σ are in units of mV. In particular, for the LIF model μ corresponds to the mean membrane potential of the neuron in absence of threshold.*

_{n}The parameters μ and σ are effective parameters that can be related to the properties of synaptic inputs. For instance, in the case of Poisson presynaptic neurons and fast postsynaptic currents, we have (Brunel and Hakim 1999)
*N _{E}* and

*N*are the numbers of excitatory and inhibitory synapses,

_{I}*J*and

_{E}*J*(with

_{I}*J*< 0) are the peak values of excitatory and inhibitory postsynaptic currents, and

_{I}*v*and

_{E}*v*are the firing rates of excitatory and inhibitory presynaptic neurons. As μ scales linearly with the number of presynaptic inputs while σ scales as the square root of the number of inputs, fluctuations in the input are large only if inhibitory and excitatory inputs balance in such a way that μ and σ are of the same order of magnitude.

_{I}In most of this study we considered fast synaptic inputs, in which case η is a white-noise process, i.e., τ* _{n}* = 0. When τ

*≠ 0, η(*

_{n}*t*) was generated using

*t*) is a white noise process.

### Conductance-Based Inputs to the Neurons

Numerical simulations were performed to verify that the results obtained for current-based inputs under the diffusion approximation hold also in the more realistic case of conductance-based synapses of finite amplitude.

The inputs to the neurons consisted of a barrage of excitatory and inhibitory synaptic inputs, so that the total input current is given by:
*E _{e}* and

*E*are the reversal potentials of excitatory and inhibitory synapses and

_{i}*g*(

_{e}*t*) and

*g*(

_{i}*t*) are synaptic conductances given by

*t*

_{k}

^{(e)}and

*t*

_{k}

^{(i)}are the times of excitatory and inhibitory synaptic inputs, respectively. The synaptic strength is measured by the dimensionless parameters

*a*and

_{e}*a*, and the synaptic decay times are instantaneous.

_{i}The diffusion approximation holds also for conductance-based synaptic inputs (Richardson 2004). Within the diffusion approximation, the difference between a conductance-based input and a current-based input is, however, merely a rescaling of the membrane time constant (Richardson 2004).

In terms of synaptic parameters, the mean μ of the effective total input, its SD σ, and the effective membrane time constant τ* _{m}* are given by (Richardson 2004):

*v*and

_{e}*v*are the total rates of excitatory and inhibitory synaptic inputs.

_{i}The domain of validity of the diffusion approximation was examined by varying the strength of excitatory synapses over an order of magnitude while keeping the mean and CV of the firing constant. The mean firing rate and the CV are determined by the three parameters μ, σ, and τ* _{m}*, while the synaptic inputs are specified by four parameters

*a*,

_{e}*a*,

_{i}*v*, and

_{e}*v*. The excitatory synaptic weight can therefore be varied while keeping μ, σ, and τ

_{i}*fixed. However, the range within which*

_{m}*a*can be varied is in practice limited by the constraints that

_{e}*a*,

_{i}*v*, and

_{e}*v*be positive. To obtained the desired mean and CV of the firing, we therefore also varied the value of the reset potential. The parameters used in the simulations are displayed in Table 1.

_{i}### Description of the Dynamics

In most of this study, we used a mathematical description of the membrane potential dynamics to analyze the ISI distribution of the LIF neuron. As we are interested in the dynamics between two APs, we consider an initial condition where at time *t* = 0 the neuron has just emitted an AP and its membrane potential lies at the reset value *V _{r}*. We then follow the dynamics until the membrane potential crosses the threshold, at which point the membrane potential is reset to

*V*, and a new “trial” begins. The fluctuating input current is a stochastic process independent from trial to trial.

_{r}The ISI distribution is equivalent to the distribution of times when the membrane potential crosses the threshold. Instead of directly studying the first passage time distribution, here we adopt a complementary approach and focus on the stochastic dynamics of the membrane potential in presence of an absorbing boundary condition at the threshold *V _{th}*. We then introduce the probability density

*P*(

*V*,

*t*) for the corresponding stochastic variable to reach

*V*at time

*t*[see Risken (1984) and Tuckwell (1988)] starting from

*V*at

_{r}*time 0*. For simplicity, in results we refer to

*P*(

*V*,

*t*) as the membrane potential distribution.

The initial condition imposes
*T*, it does not contribute to *P*(*V*, *t*) for *t* > *T* so that *P*(*V*, *t*) = 0 for *V* ≥ *V _{th}* and in particular

The probability of emitting a spike at time *T*, i.e., the ISI distribution, is given by the probability current through the threshold at time *T* (Tuckwell 1988)
*V _{th}*, the probability density

*P*(

*V*,

*t*) obeys

*t*. In particular

*Eq. 21*implies that ∫

_{−∞}

^{Vth}

*P*(

*V*,

*t*)d

*V*< 1.

The dynamics of the distribution *P*(*V*, *t*) obey the Fokker-Planck equation (Risken 1984):
*L* is a linear partial differential operator that depends on the neuronal model, and τ* _{m}* is the membrane time constant of the neuron.

For generalized integrate-and-fire models receiving a white noise current input, the Fokker-Planck operator *L* depends on a single variable, the membrane potential *V*, and reads (Fourcaud-Trocmé et al. 2003)
*V*) is the spike-generating function (see *Eq. 2*).

It can be shown (Risken 1984) that the operator *L* with the absorbing boundary condition *f*(*V _{th}*) = 0 (see

*Eq. 19*) is autoadjoint and negative definite, so that its eigenvalues are real and strictly negative, i.e., its eigenfunctions ϕ

*obey*

_{i}*> 0.*

_{i}The eigenfunctions ϕ* _{i}* form a complete set on which the membrane potential distribution

*P*(

*V*,

*t*) can be expanded. In this eigenbasis, the formal solution of the Fokker-Planck equation reads

*are set by the initial condition of*

_{i}*Eq. 18*and can be determined using the completeness equation for the eigenfunctions ϕ

*. The ISI distribution is then obtained by combining*

_{i}*Eqs. 20*and

*25*, and is given by a sum of decaying exponentials [for an alternative derivation, see Ricciardi and Sato (1988)]:

As the eigenvalues (λ* _{i}*) are strictly positive, the distribution

*P*(

*V*,

*t*) does not reach a steady state but instead decays asymptotically to zero for large

*t*. Moreover, if a finite gap is present between the smallest (dominant) eigenvalue λ

_{1}and the next (subdominant) eigenvalue λ

_{2}, the dynamics are dominated by a single exponential for large

*t*:

At long times the membrane potential distribution *P*(*V*, *t*) thus factorizes in a product of a time-independent function ϕ_{1}(*V*), and an exponential decay term *e*^{−λ1t/τm}. The time-independent function ϕ_{1}(*V*) is called a quasistationary distribution in the mathematical literature (Collet et al. 1995). In results, we refer to it as the resting distribution of membrane potentials and we denote it by *P*_{1}(*V*) = ϕ_{1}(*V*).

### Eigenvalue Decomposition of the Fokker-Planck Operator

#### LIF model.

Writing *y* = (*V* − μ)/σ, *y _{r}* = (

*V*− μ)/σ and

_{r}*y*= (

_{th}*V*− μ)/σ, the Fokker-Planck operator

_{th}*L*is given by

For a given λ, the equation
_{λ}(*y*) that is integrable on (−∞, *y*_{th}). It decreases exponentially fast [like *y* ^{λ} exp(−*y*^{2})] as *y* → ∞ and can be expressed as a linear combination of confluent hypergeometric functions (Brunel and Hakim 1999):
*M*[*a*, *b*, *c*] is the confluent hypergeometric function (Abramowitz and Stegun 1970), and *N*(λ) ensures 〈ϕ_{λ}|ϕ_{λ}〉 = 1.

The boundary condition ϕ_{λ}(*y _{th}*) = 0 determines the eigenvalues, which are strictly positive and form a discrete set {λ

_{i}} (Ricciardi and Sato 1988). In the following, we therefore write ϕ

*≡ ϕ*

_{i}_{λi}. The ϕ

*are normalized so that 〈ϕ*

_{i}*|ϕ*

_{i}*〉 = δ*

_{j}*.*

_{ij}Expanding the membrane potential distribution *P*(*y*, *t*) on the complete set formed by the functions {ϕ* _{i}*}, we write

_{i}are determined by the initial condition, which, from the completeness relation for the set {ϕ

*}, implies α*

_{i}*= e*

_{i}^{yr2}ϕ

*(*

_{i}*y*).

_{r}In the limit of vanishing firing rate obtained for μ → −∞, i.e., *y _{th}* → ∞, the LIF dynamics become equivalent to the harmonic oscillator, so that the eigenvalues take integer values λ

*→*

_{i}*i*− 1 (up to terms exponentially small in

*y*

_{th}

^{2}), and, ϕ

*→*

_{i}*H*

_{i−1}/

*N*(λ), where

*H*is the

_{i}*i*Hermite polynomial. In particular, since

_{th}*H*

_{1}(

*y*) = 2

*y*, ϕ

_{2}(

*y*) is proportional to

_{r}*y*, and changes sign for

_{r}*y*= 0, i.e.,

_{r}*V*= μ. This approximation is however, accurate only for very low firing rates (

_{r}*v*< 1 Hz).

#### EIF model.

For the EIF model, to our knowledge the eigenfunctions of the Fokker-Planck operator cannot be written explicitly. Instead, we integrate the eigenvalue equation numerically using the method introduced in (Richardson 2007).

The eigenvalue equation is

Introducing the probability density *J* defined by

This system of equation can be integrated by choosing *V _{min}* <

*V*and

_{r}*V*>

_{max}*V*, and integrating the system of equations from

_{th}*V*to

_{min}*V*. The initial condition is ϕ(

_{max}*V*) = 1 and

_{min}*J*(

*V*) = 0, corresponding to the constraint that

_{min}*J*(

*V*) decays exponentially for

*V*→ −∞. The result of this numerical integration is the value of ϕ(

*V*). The location of zeros of ϕ(

_{max}*V*) as a function of λ determines the eigenvalues we are looking for.

_{max}### Approximating the ISI Distribution

The full ISI distribution can be approximated by truncating the sum of exponentials in *Eq. 26* and keeping only the sum *P*_{ISI}^{(n)} of the first *n* terms (Alili et al. 2005). Increasing the number of terms improves the approximation at short intervals. However, for some parameters, the finite sum can take negative values for short intervals, therefore, we impose *P*_{ISI}^{(n)}(*T*) = 0 for *T* < *T*_{0}, where *T*_{0} is the largest zero of Σ_{i = 0}^{n} α* _{i}*∂

*ϕ*

_{V}*(*

_{i}*Vt*)

_{h}*e*

^{−λit/τm}.

The approximation we use is therefore

### Memory of the Previous Spike

The memory *t _{mem}* of the previous spike is defined as the time at which the dominant term in

*P*

_{ISI}

^{(n)}is larger by an order of magnitude than any other term:

### Spike-Train Autocorrelation Function

The autocorrelation function of a spike train is defined as

The Fourier transform *Ã*(ω) of the autocorrelation function, also called the power-spectrum, is related to the Fourier transform *P̃ _{ISI}*(ω) of the ISI distribution via (Gerstner and Kistler 2002)

*P̃*(ω) is known analytically (Tuckwell 1988), so the autocorrelation function

_{ISI}*A*(

*t*) in time can be determined with the help of the fast Fourier transform (Ostojic et al. 2009).

## RESULTS

### Stochastic Dynamics of the Membrane Potential

We consider a single neuron receiving a barrage of background synaptic inputs due to the spontaneous activity of the presynaptic network, as typically observed in vivo (Anderson et al. 2000; Destexhe et al. 2001, 2003). Instead of taking into account individual synaptic inputs, we represent their compound effect as an ongoing fluctuating input to the neuron. Indeed, if the amplitudes of individual synaptic inputs are small, the compound input current resulting from many synaptic inputs can be represented by a Gaussian stochastic process in time (Stein 1965, 1967; Tuckwell 1988). In most of this study, we adopt this diffusion approximation, which is a good description of the activity in many cortical areas (Anderson et al. 2000; Destexhe et al. 2001, 2003), although not all (DeWeese and Zador 2006). We moreover consider a stationary situation in which the mean and SD of the compound input do not vary in time. The values of this mean and SD depend on the conductances and reversal potentials of inhibitory and excitatory synapses as well as on the firing rates of presynaptic cells (Richardson 2004). If inhibitory and excitatory inputs are approximately balanced, the SD of the compound input can be large compared with its mean. Here we explore the effects of such a fluctuating input by varying directly the mean and SD of the effective compound input rather than the parameters of the presynaptic inputs.

The neuron receiving the fluctuating current is represented either as a generalized integrate-and-fire neuron or a WB conductance-based neuron (see materials and methods). Our aim is to study the distribution of ISIs produced by these models. We therefore examine the dynamics of the membrane potential between two APs, starting from the reset due to the previous AP, to the time when the membrane potential crosses the threshold for AP generation. We call such a period of time a trial. As the background synaptic inputs arrive at random times, the precise trajectory of the membrane potential is different on every trial. We therefore examine the statistics over many such trials.

Independently of the details of the neuronal model, if the mean input current μ is far above threshold, the membrane potential dynamics are close to deterministic. The membrane potential crosses the threshold at approximately the same time in every trial, and the firing is close to periodic (see Fig. 1*A*, *i*). In contrast, if the mean input current μ is below threshold, the membrane potential will follow stochastic dynamics, and the threshold-crossing times are widely distributed as shown in Fig. 1*B*, *ii–iii* (Troyer and Miller 1997).

To study the distribution of ISIs produced by different levels of fluctuating inputs, we use analytic tools that we combine with numerical simulations. Averaging over many trials, the distribution of ISIs is equivalent to the distribution of times at which the membrane potential crosses the threshold for the first time. In the mathematical literature, this is a well-studied quantity called the first-passage time distribution. In particular, the ISI distribution for the LIF neuron receiving a white-noise input is equivalent to the first-passage time distribution of the Ornstein-Uhlenbeck process and has been extensively investigated (Tuckwell 1988; Ricciardi 1977).

Instead of directly studying the first passage time distribution, here we adopt a complementary approach and focus on the membrane potential distribution *P*(*V*, *t*), which directly determines the ISI distribution (see materials and methods). The dynamics of this distribution obey the Fokker-Planck equation, which is essentially the counterpart of the membrane potential equation: the membrane potential equation describes the dynamics in the absence of background synaptic inputs, while the Fokker-Planck equation describes the stochastic dynamics in the presence of fluctuating inputs. The Fokker-Planck equation is very general and can be written for any neuronal model receiving a diffusive input. However, analyzing the Fokker-Planck equation is in practice difficult if the dynamics depend on variables other than the membrane potential, such as the conductances of various channels. For that reason, we first concentrate on integrate-and-fire models of a single variable, the LIF and EIF models, receiving an input uncorrelated in time. In a second stage, we use numerical simulations to verify that the findings obtained for these models are valid also for a conductance-based neuronal model with several dynamical variables, for fluctuating inputs correlated in time, and for full conductance-based synaptic inputs.

### Resting State and Exponential Tail of the ISI Distribution

A mathematical analysis of the membrane potential dynamics reveals that, under very general conditions, the ISI distribution of integrate-and-fire models can be written as sums of decaying exponentials (see material and methods for details):
* _{m}* is the membrane potential,

*v*are decay rates corresponding to eigenvalues of the Fokker-Planck equation, and β

_{i}*are constant coefficients.*

_{i}From *Eq. 43*, it is apparent that for large ISIs the distribution decays exponentially, a property generically observed in experimental data. Indeed if the smallest (dominant) decay rate *v*_{1} and the next (subdominant) decay rate *v*_{2} are separated by a finite difference, for *t* large the dynamics are dominated by a single exponential:

From the point of view of membrane potential dynamics, the origin of this exponential decay can be traced back to the fact that for sufficiently long times since the previous spike, the membrane potential converges to a resting distribution. The resting distribution is a direct counterpart of the resting potential reached by the membrane potential in absence of noise: for a subthreshold input, in absence of noise the membrane potential converges to a fixed resting potential μ, while in presence of noise the membrane potential explores a resting distribution around μ. Once that distribution is reached the probability that the membrane potential crosses the threshold is constant in time and the firing becomes a Poisson process. This mathematical description is fully consistent with the intuitive picture of Troyer and Miller (1997), who proposed that the membrane potential dynamics can be divided in two epochs, an initial transient where the dynamics are governed by postspike effects and a resting state at long times in which all memory from the previous AP has been lost.

For the LIF model, the decay rates and coefficients in *Eq. 43*, as well as the resting state, can be computed analytically (Alili et al. 2005). Details of the computation are provided in materials and methods, and Fig. 1 displays an illustration of the convergence to the resting state. The LIF model in the presence of fluctuating inputs possesses only two independent parameters, which we choose to be the mean μ and SD σ of the membrane potential. We examined the dominant decay rate and resting membrane potential distribution while systematically varying μ and σ. To allow for a direct comparison between models, we parametrized the results by the mean firing rate *v* and CV of the ISI distribution, rather than μ and σ, as the two sets of parameters (*v*, CV) and (μ, σ) are related by a one-to-one mapping (Vilela and Lindner 2009b). More precisely, we increase the amplitude σ of input fluctuations, while keeping *v* fixed by decreasing μ. This is equivalent to increasing the CV at fixed firing rate.

The value of the dominant decay rate *v*_{1} and the shape of the resting distribution are tightly related, as the dominant decay rate represents the rate at which the membrane potential crosses the threshold while exploring the resting distribution. Consequently, the more the resting distribution is concentrated near the threshold, the larger the dominant decay rate. Figure 2*A* displays *v*_{1} as a function of the CV of the firing for three different values of the mean firing rate *v*, while Fig. 2*B* displays the resting distribution for *v* = 30 Hz and three values of the CV. Both the dominant decay rate and the resting distribution depend on the level of input fluctuations. For weak input fluctuations (small values of the CV), the weight of the resting distribution is concentrated near the threshold *V _{th}* = 20 mV, and correspondingly the dominant decay rate is large. As the CV increases, the maximum of the resting distribution shifts away from the threshold and

*v*

_{1}decreases. The position of the maximum is approximately given by the value of μ, the effective resting potential that is determined by the mean of the fluctuating input (see

*Eq. 9*in materials and methods). Interestingly, the dominant decay rate rescaled by the mean firing rate appears to depend only on the value of the CV and not on the value of the mean firing rate

*v*itself (Fig. 2

*A*). Moreover,

*v*

_{1}/

*v*is approximately a power-law function of CV with an exponent close to −1.2.

### Shape of the Full ISI Distribution

Having described the ISI distribution for long ISIs, we now turn to the shape of the full distribution. *Equation 43* shows that the full ISI distribution can be approximated by a sum of exponentials, the coefficients and decay rates of which can be calculated without any fitting parameters. Here we describe the results obtained from this approximation for the LIF model.

For weak input fluctuations (Fig. 1*Bi*, CV = 0.22), the ISI distribution is narrowly distributed around its mean 〈*T*〉 =1/*v*, and the dominant exponential decay describes only the steep tail of the distribution. Including additional exponentials improves the approximation at short intervals. However, as the firing is close to periodic, the actual ISI distribution vanishes at short intervals, a behavior that cannot be reproduced with a finite sum of exponentials. We therefore set the analytical approximation of the ISI distribution to zero at short intervals (see materials and methods).

As the amplitude of input fluctuations is increased (while keeping the mean firing rate constant by decreasing the mean input), the ISI distribution becomes broader and the dominant exponential decay describes an increasingly large portion of the full ISI distribution. For intermediate input fluctuations (Fig. 1*Bii*, CV = 0.75), *P _{ISI}* (

*T*) displays a broad maximum and is well described by a sum of just two exponentials. The behavior at long intervals is well captured by the dominant exponential decay. At short intervals, the dominant approximation overestimates the ISI distribution, due to refractory effects. These refractory effects are purely due to the reset, as the effective refractory period is much larger than the absolute refractory period of the neuron, its actual value being dependent on the level of noise. This refractory period is taken into account by the subdominant exponential, the coefficient of which is negative.

For strong input fluctuations (Fig. 1*Biii*, CV = 1.2), the qualitative behavior of the ISI distribution changes. The ISI distribution becomes monotonically decaying, and at short times the single exponential approximation underestimates the true distribution. This observation indicates that the probability of firing per unit time is larger right after the previous spike than at longer times since the previous spike. In other words, the neuron tends to fire in burst due to the strongly fluctuating input (Barbieri and Brunel 2007; Vilela and Lindner 2009a). Including a second exponential takes into account this bursting, and the coefficient of this subdominant term is now positive. The transition from refractory firing to bursting thus takes place when the coefficient of the subdominant term changes signs, which approximately happens when the effective resting potential μ crosses below the reset potential *V _{r}*. The emergence of bursting can therefore be explained by the observation that for strong input fluctuations; the maximum of the resting distribution lies below the reset potential (cf. Fig. 2

*B*). Intuitively, if the reset is above the mean of the resting distribution of membrane potential, then immediately following a spike the membrane potential becomes closer to threshold than it would be on average, enhancing the likelihood of a short ISI than one would expect from the average firing rate.

In the intuitive picture of stochastic membrane potential dynamics proposed by Troyer and Miller (1997), the degree of irregularity of the firing is determined by the relative amount of time spent in the postspike transient regime and the resting distribution: the faster the resting state is reached, the more variable the firing. Our approximation of the ISI distribution allows us to actually compute the time needed to reach the resting distribution. We call this time the memory *t _{mem}* of the previous AP. Plotting this time as a function of the input fluctuations complements the intuitive picture and provides an additional characterization of the ISI distribution as it quantifies the range of ISIs in which the distribution is dominated by a single exponential.

Figure 3 represents the memory time *t _{mem}* in units of the mean ISI 〈

*T*〉 = 1/

*v*, as a function of the CV of the ISI distribution, for different values of the mean firing rate

*v*. While one might intuitively expect that

*t*decreases monotonically as the SD of input fluctuations is increased (as does the dominant decay rate, see Fig. 1

_{mem}*A*), this is not the case. As the CV is progressively increased,

*t*first decreases as expected but eventually reaches a minimum. This first branch corresponds to periodic and refractory firing. In that parameter range

_{mem}*t*lies close to the maximum of the ISI distribution (Fig. 1

_{mem}*B*), and as the amplitude of input fluctuation is increased (while decreasing the mean input μ to keep the firing rate constant) the ISI distribution is increasingly dominated by a single exponential. The location of the minimum of

*t*is given by the point where the coefficient β

_{mem}_{2}of the subdominant exponential vanishes, which happens approximately when the effective resting potential μ crosses the reset potential

*V*. An intuitive explanation is that when the reset potential is close to the peak of the resting distribution (as is the case when the reset potential is close to the effective resting potential) one expects the time to relax to the resting distribution to be minimal. As the CV is further increased beyond the location of the minimum,

_{r}*t*increases, reflecting the larger contribution of the subdominant term, and the increasing tendency of the neuron to fire in bursts. Note that the regime of increasing

_{mem}*t*and burst firing was not anticipated from the intuitive picture of Troyer and Miller. Interestingly, the minimum of

_{mem}*t*lies at a CV <1, so that for unit CV the firing of the neuron is not Poisson but rather bursting.

_{mem}So far we have examined the shape of the ISI distribution as the CV is increased at fixed firing rate *v*. The effect of the firing rate on the shape of the ISI distribution is displayed in Fig. 4, where we show the distribution of ISIs rescaled by the mean ISI 〈*T*〉 = 1/*v*, for several values of *v*. For a fixed CV, varying the mean firing rate affects only weakly the shape of the distribution: for large firing rates, the rescaled distributions are almost indistinguishable, while at small *v* the deviations become more noticeable. While the dominant decay rate appears to be essentially independent of the firing rate (see Fig. 2), subdominant terms depend more strongly on *v*.

Our analysis shows that for a given set of parameters for the fluctuating input, the ISI distribution is never exactly exponential. However, as the mean firing rate *v* is decreased the memory time decreases, so that the full ISI distribution becomes closer and closer to an exponential, but it becomes rigorously exponential only in the limit *v* → 0 (Nobile et al. 1985), and the convergence to this limit is rather slow. Indeed Fig. 3 shows that for a mean firing rate of *v* = 10 Hz the convergence time can still be significant while for a Poisson distribution it is vanishingly small.

### Comparison with the Gamma Distribution

We have identified three main features of the family of ISI distributions of spiking neurons in presence of fluctuating inputs: *1*) the distribution decays exponentially for large intervals; *2*) as the CV is increased, the distribution progressively evolves from narrow to refractory-exponential to bursting-exponential; and *3*) the shape of the ISI distribution is approximately independent of the firing rate. These properties are very similar to those of gamma distributions, a two-parameter family of distributions with probability densities given by

For the ISI distribution of spiking neurons, we have observed that the dominant decay rate *v*_{1}, the mean firing rate *v* and the coefficient of variation are related approximately via *v*_{1}/*v* = CV^{−1.2} (cf. Fig. 2*A*). For the gamma distributions, *v*_{1} = β, *v* = β/α, and CV = 1/α; hence, *v*_{1}/*v* = CV^{−2}. The two scaling relationships are comparable, but the exponents are clearly different. While the two families of distributions are qualitatively similar, our analysis shows that they are mathematically distinct. A direct comparison between the distributions however, shows that at fixed mean and CV, ISI distributions of spiking neurons deviate only weakly from gamma distributions (see Fig. 4), as the two families of distributions have similar asymptotics. We therefore conclude that gamma distributions provide a natural first approximation to ISI distributions of spiking neurons.

### ISI distributions Are Approximately Independent of the Neuronal Model

The analysis of the LIF model identified a progression between three types of ISI distributions as the input fluctuations are increased at a fixed firing rate: *1*) a narrow distribution, *2*) an exponential with an effective refractory period, and *3*) an exponential with bursts. Here we examine whether these three types of distributions are also present in more biophysically realistic spiking models. One possible concern is that some aspects of the firing statistics of the LIF neuron could originate from the hard threshold used to generate spikes in that model and that more realistic models might lead to ISI distributions of different shapes.

We first consider the EIF model, in which APs are initiated by a sodium-like, exponential conductance (Fourcaud-Trocmé et al. 2003). Compared with the LIF model, the EIF model includes one additional parameter, the spike-sharpness Δ* _{T}*, which allows it to produce APs with a range of different shapes. Fits to experimental data (Badel et al. 2008) indicate that for pyramidal cells Δ

*≈ 1 mV, while the Hodgkin-Huxley model corresponds to Δ*

_{T}*≈ 3.5 mV (Fourcaud-Trocmé et al. 2003), and the LIF model is recovered for Δ*

_{T}*= 0.*

_{T}The EIF model is a one-dimensional model, so its dynamics can be analyzed with the help of the Fokker-Planck equation. Similarly to the LIF model, the ISI distribution can be written as a sum of exponentials, and the corresponding decay rates can be computed numerically using the method introduced by Richardson (2004) (see materials and methods). This computation confirms that a gap is present between the dominant and the subdominant decay rate, implying that at long times the dynamics converge to a resting state and the ISI distribution decays exponentially with a rate *c*_{1}. For a given mean firing rate *v* and a given CV, the values of the dominant decay rate *v*_{1} for the EIF model are approximately equal to those of the LIF model, as shown in Fig. 2*A*. Numerical simulations in fact show that for a given *v* and CV the full ISI distributions of the LIF and EIF model are very similar to each other (Fig. 5). These observations hold independently of the value of the spike sharpness parameter Δ*t*.

We next turn to the WB model, a Hodgkin-Huxley like conductance-based model (Wang and Buzsáki 1996). It has been previously found that the f-I curve and the linear response function of the WB model can be reproduced by an EIF model with appropriately chosen values of the spike initiation parameter Δ*t*, reset potential *V _{r}*, and refractory time

*τ*(Fourcaud-Trocmé et al. 2003). Here we examined the ISI distribution of the WB model receiving a white-noise input. For given values of the mean firing rate

_{rp}*v*and CV, the ISI distribution of the WB model closely overlaps those of the LIF and EIF model (see Fig. 5). The most apparent difference is seen at short ISIs in the bursting regime (see Fig. 5

*Aiii*): due to a different postspike reset mechanism in the WB model compared with the integrate-and-fire models, the distributions are different for ISIs shorter than 5 ms. While the ISI distributions of the WB model and the fitted EIF model could have been expected to be similar, it is somewhat surprising that the ISI distribution of the WB model is so close to the LIF model and EIF model with sharper spikes, as the firing rate dynamics of these models differ significantly (Fourcaud-Trocmé et al. 2003; Ostojic and Brunel 2011).

Our results strongly suggest that the first two moments essentially determine the full ISI distribution independently of the model considered. This implies that the analytic results obtained for the LIF model extend directly to the EIF and WB model. In particular, these models also exhibit an exponential decay of the ISI distribution at long times due to the loss of memory of the previous spike, and the shape of the full ISI distributions depends approximately only on the value of the CV.

It should be noted that the finding that the ISI distributions depend weakly on the model implies that the spike-train autocorrelation functions are also very similar for the three models. The ISI distributions and the autocorrelation functions are indeed directly related (see materials and methods). A direct comparison between the autocorrelation functions of the LIF, EIF, and WB models indeed confirms that they are very close to each other (see Fig. 5). This observation complements the earlier finding that the power spectra of the LIF model, the non-leaky IF model, and the quadratic integrate-and-fire model are very similar for given values of the mean firing rate *v* and CV (Vilela and Lindner 2009a). Note also that the autocorrelation functions display qualitatively different features in the three firing regimes that we have distinguished: in the regular firing regime the autocorrelation is oscillatory (Fig. 5*Bi*); in the refractory-exponential regime the autocorrelation has a dip at short intervals (Fig. 5*Bii*); and in the bursting-exponential regime the autocorrelation has a bump at short intervals (Fig. 5*Biii*).

### Effects of Correlated Noise

So far we have examined the ISI distributions only for white noise inputs, i.e., in the limit of very fast synaptic timescales. Finite synaptic timescales lead to fluctuating inputs that are correlated in time, and previous studies have found that such correlations can qualitatively change the dynamics (Brunel et al. 2001; Fourcaud and Brunel 2002; Moreno-Bote and Parga 2004, 2006; Schwalger and Schimansky-Geier 2008). In particular, because of temporal correlations in the inputs, the firing is no longer a renewal process and consecutive ISIs can be correlated (Schwalger and Schimansky-Geier 2008), an effect we do not investigate here.

Here we examine to what extent the results we obtained in the case of uncorrelated inputs extend to the more realistic case of a Gaussian stochastic process with a finite correlation time τ* _{n}*. In that situation, the Fokker-Planck equation depends on two variables, and the boundary conditions become complex (Brunel et al. 2001; Fourcaud and Brunel 2002). The limit of large correlation times has been studied analytically (Schwalger and Schimansky-Geier 2008), and it has been found that in this case the ISI distribution decays as a power-law for large ISIs. That behavior, however, dominates only for very large correlation times (τ

*> 100 τ*

_{n}*, where τ*

_{m}*is the membrane time constant). Here we focus on shorter correlation times corresponding to synaptic filtering. We determined the ISI distribution for different models using numerical simulations.*

_{m}Figure 6 displays the ISI distributions for the LIF, EIF, and WB models for two different values of the input correlation time, τ* _{n}* = 6 ms and τ

*= 12 ms. Correlated noise modifies the ISI distributions with respect to the case of white-noise inputs, and the precise shape of the ISI distributions depends on the value of the input correlation time τ*

_{n}*. Yet, for a given τ*

_{n}*, the main results obtained in the uncorrelated situation hold:*

_{n}*1*) for long ISIs, the ISI distribution becomes exponential;

*2*) as the CV of the firing is increased at fixed firing rate, the ISI distribution progressively evolves from a narrow distribution to a refractory exponential and finally to a bursting exponential; and

*3*) for given firing rate

*v*and CV, the ISI distributions of the LIF, EIF, and WB models lie very close to each other.

### Effects of Conductance-Based Synaptic Inputs

The results presented in this study are based on the assumption that the compound effect of many synaptic inputs can be represented as a fluctuating input current to the neuron. This assumption is expected to be correct if the amplitude of individual postsynaptic currents is small with respect to the threshold of the neuron. In addition to the influx of current, synaptic inputs also lead to an increase of the membrane conductance. In the case of a constant barrage of weak synaptic inputs, it has been argued that the effect of this conductance increase can be correctly taken into account simply by decreasing the effective membrane time constant of the neuron (Richardson 2004), and this is the approximation we used.

To examine the validity of these approximations, we performed numerical simulations of a LIF neuron receiving conductance-based excitatory and inhibitory synaptic inputs (Richardson 2004; Burkitt 2001). We varied the peak conductances of the synapses and the firing rates of the presynaptic neurons in such a way as to obtain a postsynaptic firing rate of 30 Hz and three different levels of postsynaptic firing variability corresponding to CV = 0.2, 0.8, and 1.2 as in Figs. 1, 5, and 6. This essentially amounts to using four synaptic parameters to fix the three parameters corresponding to the mean and SD of the equivalent fluctuating input and the effective time constant, leaving one degree of freedom, which we chose to be the conductance of excitatory synapses. We therefore determined the ISI distributions corresponding to fixed mean and CV but three different strengths of excitatory synapses varying over an order of magnitude (denoted by 1×, 5×, and 10×) and corresponding to unitary excitatory postsynaptic potentials varying approximately from 0.01 to 0.1 mV. For some of the parameter values, the diffusion approximation is not quantitatively correct, in that the values of the mean and SD of the fluctuating input computed from the synaptic parameters do not predict the correct firing rate and CV. Nevertheless, for fixed mean firing rate and CV, the three parameter sets lead to indistinguishable ISI distributions that are perfectly reproduced with a model receiving a fluctuating input current (see Fig. 7).

An interesting implication of conductance-based synapses is that the mean and SD of the equivalent fluctuating current are constrained by the requirement that synaptic conductances must be positive. For fixed neural parameters, this constraint induces a limit on how irregular the firing of the postsynaptic neuron can be (Chance et al. 2002). However, the maximal attainable value of the CV depends on the neural parameters such as the reset and the threshold, as these parameters effectively set the relevant scale of fluctuation and determine whether an input is subreset [see also, Troyer and Miller (1997)]. In conclusion, while for neurons receiving current-based inputs the values of the reset and threshold potential play no role for the range of possible ISI distributions, this is not true for neurons with conductance-based inputs.

## DISCUSSION

In this study, we systematically examined how fluctuating inputs resulting from a balance between excitation and inhibition shape the distribution of ISIs produced by spiking neuron models. Using analytical tools for the LIF model, we showed that the ISI distribution can be accurately approximated by a sum of exponentials. For long ISIs, the distribution generically displays an exponential decay, corresponding to the fact that the membrane potential distribution converges to a resting distribution. The behavior at shorter ISIs depends on the magnitude of the fluctuating input and progressively evolves between three classes of features. If the mean input is suprathreshold and the fluctuations are weak, the ISI distribution is dominated by the effects of periodicity in the firing and is in consequence narrow and close to a Gaussian. If the mean input is between threshold and reset values, the ISI distribution is dominantly exponential but exhibits refractory effects at short intervals. If the mean input is below reset, the ISI distribution is dominantly exponential, with an increased tendency to fire at short intervals since the previous spike. Fluctuating inputs lead only to these three possible shapes of ISI distributions, the actual shape being determined by the CV of the ISIs and independent of the mean firing rate. The ISI distributions of LIF neurons are therefore well approximated by gamma distributions.

Numerical simulations show that the properties of the ISI distributions identified for the LIF model also hold for other spiking neuron models such as the EIF model and a Hodgkin-Huxley type conductance-based model. In addition, we observed that the full ISI distributions were essentially identical for the three models at equal mean firing rate and ISI CV, so that the first two moments seem to determine the full ISI distribution for a large class of models. These findings are moreover independent of whether the fluctuating input is modeled as white or correlated noise and of whether synapses are current or conductance based.

### Comparison with Previous Studies

The ISI distributions of model neurons have been the subject of a large number of previous works. In a pioneering study, Gerstein and Mandelbrot (1964) modeled the membrane potential dynamics as a random walk and found that the ISI distribution exhibits a power-law behavior. The model they examined corresponds to an integrate-and-fire model without leak. This model can be studied with the mathematical tools we used for the LIF model. However, in absence of a leak current, the ISI distribution cannot be written as a discrete sum of exponentials because the spectrum of the Fokker-Planck operator is continuous and no gap is present around the dominant decay rate. In contrast, as soon as a confining potential is present, such as a leak or even simply a lower bound on the membrane potential (Fusi and Mattia 1999), the spectrum becomes discrete, and our results hold.

In the mathematical literature, the study of ISI distributions is part of the large class of first-passage problems. In particular, the ISI distribution of the LIF model is equivalent to the first-passage time distribution of the Ornstein-Uhlenbeck process and has been extensively studied (Tuckwell 1988). Most of these studies directly examined the first-passage time distribution, the Laplace transform of which is known analytically (Siegert 1951), and focused on the moments of the ISI distribution. Ricciardi and Sato (1988) noted that the first-passage time distribution can be written as a sum of exponentials, the exponents being the poles of the Laplace transform of the distribution. Studying the ISI distribution with the help of the Fokker-Planck equation for the membrane potential dynamics allowed us to appreciate that this property is valid not only for the LIF model but for a much larger class of models. To our knowledge, only one study related the ISI distribution to the eigenvalues of the Fokker-Planck equation for the membrane potential dynamics (Alili et al. 2005), but the authors did not systematically examine the effects of the fluctuating input on the shape of the ISI distribution.

In recent years, several studies have examined dynamics of integrate-and-fire neurons using a similar mathematical method to ours, eigen-decompositions of the Fokker-Planck equation (Knight et al. 2000; Knight 2000; Mattia and Giudice 2002, 2004). While these studies exploited the same tool as we did, they addressed questions related to rate dynamics rather than ISI statistics. When studying rate dynamics, the boundary conditions in the Fokker-Planck equation need to take into account the reset following the firing of an action-potential [see e.g., Brunel and Hakim (1999)], while this is not the case for first-passage statistics. Because the boundary conditions are different, the spectrum of the Fokker-Planck operator is different in the cases of rate dynamics and first-passage distributions. In particular, in the case of first-passage distribution, the spectrum is always real and positive, while for rate dynamics eigenvalues can take on complex values.

### Interpreting Experimental Data

We have shown that the family of gamma distributions is a natural first approximation to the range of ISI distributions spanned by varying the parameters of the fluctuating inputs to the neurons. Gamma distributions have been used to fit ISI distributions measured from cortical neurons (Barbieri et al. 2001; Miura et al. 2007; Maimon and Assad 2009), which clearly indicates that the ISI distributions produced by spiking neuron models correspond closely to experimental findings. In particular, Miura et al. (2007) used dynamic clamp to inject fluctuating conductance-based inputs to neurons in cortical slices and found that the ISI distributions were well described by gamma distributions in direct agreement with our predictions. However, in that study the authors found only gamma distributions with α > 1 (see *Eq. 45*), i.e., they did not observe bursting exponential distributions, while in vivo recordings lead to distributions with both α > 1 and α < 1 (Maimon and Assad 2009). We predict that gamma distributions with α < 1 should be observed for subreset inputs, a range presumably not explored by Miura et al. (2007).

In many experimental recordings, the number of ISIs recorded for a given cell is not sufficient to construct a full histogram. Several previous studies have therefore focused on the moments of the ISI distribution rather than the full distribution and attempted to fit various models to experimental data (Tuckwell 1988; Lansky and Radil 1987; Inoue, Sato, and Ricciardi 1995; Shinomoto, Sakai, and Funahashi 1999; Sakai, Funahashi, and Shinomoto 1999). In a particularly extensive study, Shinomoto et al. (1999) computed the CV and skewness coefficient (equivalent to the third moment of the ISI distribution) for >600 spike trains and compared them with the values that can be obtained from the LIF model receiving an uncorrelated fluctuating input current. They numerically found that the CV and skewness coefficient (SK) of the LIF model cannot be varied independently but are instead constrained to a very narrow region of the CV − SK plane. This observation is in agreement with our finding that the value of CV essentially determines the shape of the ISI distribution. Moreover, although not explicitly shown by Shinomoto et al. (1999), the region in the CV − SK plane spanned by the LIF model is very close to the line SK = 2CV that corresponds to the gamma distribution. Examining the experimental values of the CV and skewness coefficient, the authors conclude that the LIF model with white noise inputs does not fully account for experimental data, as a number of experimental points deviate from the region in the CV − SK plane spanned by the LIF model. It should, however, be noted that >90% of experimental points lie close to the LIF region. Moreover, in a later study the authors showed that including a correlation time in the fluctuating input allows the LIF model to account for all experimental data points (Sakai, Funahashi, and Shinomoto 1999).

Recent experimental reports have suggested that the typical shape of ISI distributions vary systematically between different cortical areas (Shinomoto et al. 2005, 2009; Maimon and Assad 2009). In particular, Maimon and Assad (2009) used gamma distributions to fit the ISI distributions recorded in a primate neocortex in vivo and found that the average value of the parameter α of the distribution (see *Eq. 45*) appears to vary systematically between cortical areas, the firing being more regular (larger α) going from the visual to motor areas. An important question is whether the different types of ISI distributions may reflect different levels of intrinsic fluctuations in different cortical areas, corresponding to different details of excitation-inhibition balance in the network activity. The three types of ISI distributions we have found as the amount of noise are varied in spiking models resemble strikingly the cortical distributions found in (Maimon and Assad 2009). The observation that the magnitude of input fluctuations, rather than the details of the spike generation mechanism, determine the shape of the ISI distribution suggest that different amounts of network-generated noise in the different areas may be sufficient to explain the differences between ISI distributions in different cortical areas (Maimon and Assad 2009). In our models, the scale of intracellular fluctuations is set by the distance between threshold and reset, a parameter that might vary between different types of neurons. Teasing apart the effects of input fluctuations and cellular properties may therefore be difficult in extracellular recordings.

In this study, we have examined spiking models that fire only sodium-type spikes. While these models fire in bursts when the input fluctuations are large and the mean input is subreset, we found that they lead only to unimodal ISI distributions. In contrast, bimodal ISI distributions have been observed in some experimental studies (Bair et al. 1994; Compte et al. 2003) [but very rarely in Maimon and Assad (2009)]. Preliminary results indicate that bimodal distributions can be obtained by introducing a calcium-like conductance, the simplest model with such a conductance being the integrate-and-fire-or burst neuron (Smith et al. 2000). A calcium-like conductance leads to broad calcium spikes which trigger bursts of sodium APs, similar to the Poisson bursting model (Bair et al. 1994). Other types of mechanisms, such as subthreshold resonance, might also lead to bimodal ISI distributions (Engel et al. 2008). An additional possibility is that some aspects of the recorded ISI distributions are due to nonstationarities, an important experimental confound.

## GRANTS

This work was supported by funding from the European Community's Seventh Framework Programme through a Marie Curie International Outgoing Fellowship for Career Development.

## ACKNOWLEDGMENTS

I am deeply indebted to Vincent Hakim for valuable advice at the start of this project. I am grateful to Dan Rubin for a careful reading of the manuscript.

- Copyright © 2011 the American Physiological Society