How Noisy Adaptation of Neurons Shapes Interspike Interval Histograms and Correlations

Channel noise is the dominant intrinsic noise source of neurons causing variability in the timing of action potentials and interspike intervals (ISI). Slow adaptation currents are observed in many cells and strongly shape response properties of neurons. These currents are mediated by finite populations of ionic channels and may thus carry a substantial noise component. Here we study the effect of such adaptation noise on the ISI statistics of an integrate-and-fire model neuron by means of analytical techniques and extensive numerical simulations. We contrast this stochastic adaptation with the commonly studied case of a fast fluctuating current noise and a deterministic adaptation current (corresponding to an infinite population of adaptation channels). We derive analytical approximations for the ISI density and ISI serial correlation coefficient for both cases. For fast fluctuations and deterministic adaptation, the ISI density is well approximated by an inverse Gaussian (IG) and the ISI correlations are negative. In marked contrast, for stochastic adaptation, the density is more peaked and has a heavier tail than an IG density and the serial correlations are positive. A numerical study of the mixed case where both fast fluctuations and adaptation channel noise are present reveals a smooth transition between the analytically tractable limiting cases. Our conclusions are furthermore supported by numerical simulations of a biophysically more realistic Hodgkin-Huxley type model. Our results could be used to infer the dominant source of noise in neurons from their ISI statistics.


Introduction
The firing of action potentials of a neuron in vivo is a genuine stochastic process due to the presence of several sources of noise [1]. The spontaneous neural activity (e.g. the firing of a sensory cell in absence of sensory stimuli) [2,3] as well as the response of neurons to stimuli cannot be understood without taking into account these fluctuations [4]. Moreover, noise can have a positive influence on neural function, e.g. by stochastic resonance [5,6], gain modulation [7], decorrelation of spiking activity [8], enhancement of signal detection [9], or fast transmission of noise coded signals [10,11]. For these reasons, reduced stochastic models of neural activity have been suggested [12][13][14] and analytical methods have been developed to calculate the statistics of spontaneous activity and the response to periodic stimuli [15][16][17]. Studying such reduced models allows to relate specific mechanisms with certain statistics of neural firing. Vice versa, analytical expressions for the firing statistics of model neurons may be used to infer unknown physiological details from experimental data.
Spike-frequency adaptation is another common feature of neural dynamics that, however, is still poorly understood in the context of stochastic spike generation. Associated adaptation currents which act on time scales ranging from 50ms to seconds are ubiquitous throughout the nervous system [18]. Prominent examples of adaptation mechanisms include M-type currents, calcium-gated potassium currents (I AHP ), and slow inactivation of sodium currents. Functional roles of spike-frequency adaptation include forward masking [19], high-pass filtering [20][21][22], and response selectivity [23][24][25]. If the neuron is driven by fast fluctuations, adaptation reveals itself in the interspike interval statistics of neural firing, most prominently in the occurrence of negative correlations among interspike intervals [26][27][28][29][30][31]. These features can be phenomenologically captured in generalized integrate-and-fire (IF) models via introduction of a slow inhibitory feedback variable, either acting as a dynamic threshold or as an inhibitory conductance or current [28,29,[32][33][34] or in even more simplified models [35][36][37].
In previous studies on stochastic models with adaptation, fluctuations were considered to be fast, e.g. Poissonian synaptic spike trains passing through fast synapses or a white Gaussian input current representing a mixture of intrinsic fluctuations and background synaptic input. In particular, the dominating intrinsic source of fluctuations is ion channel noise [1,[38][39][40]. This kind of noise is not only contributed by the fast ionic conductances, which establish the spike generating mechanism, but also by the channels that mediate adaptation currents. If the number of adaptation channels is not too large, the stochastic opening and closing of single channels will contribute a fluctuating component to the adaptation current. This noise contribution, which was so far ignored in the literature, and its impact on the ISI statistics is the subject of our study. Here, we only consider the simplest adaptation channel model which corresponds to an M-type adaptation current. Our results, however, also apply to other sources of noise emerging from a slow adaptation mechanisms as, for instance, slow Ca 2z fluctuations in the case of calcium-gated potassium currents (see discussion).
In this paper, we analyze a perfect integrate-and-fire (PIF) model in which a population of N a channels mediate a stochastic adaptation current. We approximate this model by simplified stochastic differential equations (diffusion approximation). For slow adaptation, we are able to show that the latter is equivalent to a PIF neuron driven by a slow external noise. As a surprising consequence, pure adaptation channel noise induces positive ISI correlations in marked contrast to negative ISI correlations evoked by the commonly studied combination of fast noise and deterministic adaptation [28]. Furthermore, the ISI histogram is more peaked and displays a heavier tail than expected for a PIF model with fast current noise. Our results for the PIF (positive ISI correlations, peaked histograms) are qualitatively confirmed by extensive simulations of a conductance-based model with N a stochastic adaptation channels supporting the generality of our findings.

Results
Our main concern in this paper was the effect of noise associated with the slow dynamics of adaptation on the interspike interval statistics. Specifically, the noise was regarded as the result of the stochastic activity of a finite number of slow ''adaptation channels'', e.g. M-type channels ( [41,42]). We contrasted this slow adaptation channel-noise with the opposite and commonly considered case of a deterministic adaptation mechanism together with fast current fluctuations [28,29]. In both cases, the spike frequency and the variability of the single ISI were similar, although the sources of spiking variability were governed by completely different mechanisms. In a real neuron this distinction would correspond to the case where one source of variability dominates over the other.
In order to demonstrate the results on two different levels of complexity, we conducted both analytical investigations of a tractable integrate-and-fire model and a simulation study of a biophysically more realistic Hodgkin-Huxley-type neuron model. For the first model we chose the perfect (non-leaky) integrate-andfire (PIF) model [2,8]. This model represents a reasonable description in the suprathreshold firing regime, in which a neuron exhibits a stable limit cycle (tonic firing). The model was augmented with an inhibitory adaptation current mediated by a population of N a adaptation channels (Fig. 1A). For simplicity, we assumed binary channels that switch stochastically between an open and a closed state. The transition rates depend on the presence or absence of an action potential. This can be approximated by passing the membrane potential V through a steady-state activation probability, w ? (V ), that attains values close to unity during action potentials, i.e. when the voltage exceeds the threshold, and is near zero for potentials below the firing threshold.
The PIF model, however, describes only the dynamics of the subthreshold voltage. In particular, it does not explicitely yield the suprathreshold time course of the action potential, but only the time instant of its onset (given by the threshold crossings). Nevertheless, the spike-induced activation of adaptation channels can be introduced in the PIF model. To this end, we approximated the activation function w ? (V ) by a piecewise constant function of time, w ? (t), which attains unity for t AP~1 ms after the onset of each action potential and is zero otherwise (Fig. 1B, middle). A weak subthreshold activation of the adaptation current as observed for the M-current does not change the qualitative results of the paper (see Discussion). The time constant of the first-order channel kinetics was set by the parameter t w .
Although our model aims at the stationary firing statistics, we would like to stress that it exhibits spike-frequency adaptation in the presence of time-varying stimuli. In particular, spike-frequency adaptation in response to a step stimulus is retained regardless of the considered noise source, channel numbers or approximations made during the theoretical analysis (Fig. 1C). This is a nice feature of the PIF model, for which the firing rate does not depend on the nature and magnitude of the noise. This allowed us to vary the noise properties without altering the adaptation properties.

Diffusion approximation of channel noise
For the theoretical analysis the channel model describing the dynamics of each individual channel could be considerably simplified by a diffusion approximation. As shown in Methods, the dynamics of the finite population of adaptation channels can be described by (i) the deterministic adaptation current and (ii) additional Gaussian fluctuations with the same filter time as the adaptation dynamics. Together with our integrate-and-fire dynamics for the membrane potential V (Methods, Eq. (35)), we obtained a multi-dimensional Langevin model that approximates the IF model with stochastic ion channels (Methods, Eq. (20), (35)): Put differently, a finite population of slow adaptation channels (instead of an infinite population and hence a deterministic

Author Summary
Neurons of sensory systems encode signals from the environment by sequences of electrical pulses -so-called spikes. This coding of information is fundamentally limited by the presence of intrinsic neural noise. One major noise source is ''channel noise'' that is generated by the random activity of various types of ion channels in the cell membrane. Slow adaptation currents can also be a source of channel noise. Adaptation currents profoundly shape the signal transmission properties of a neuron by emphasizing fast changes in the stimulus but adapting the spiking frequency to slow stimulus components. Here, we mathematically analyze the effects of both slow adaptation channel noise and fast channel noise on the statistics of spike times in adapting neuron models. Surprisingly, the two noise sources result in qualitatively different distributions and correlations of time intervals between spikes. Our findings add a novel aspect to the function of adaptation currents and can also be used to experimentally distinguish adaptation noise and fast channel noise on the basis of spike sequences.
adaption dynamics) entails the presence of an additional noise g(t) with a correlation time t w (time scale of the deterministic adaptation) and a variance which is inversely proportional to the number of channels. The membrane potential V of the PIF model is thus driven by four processes: (i) the base current m, (ii) the white current fluctuations j(t) of intensity D (representing an applied current stimulus, channel noise originating from fast sodium or delayed-rectifier potassium currents, or shot-noise synaptic background input), (iii) the slow noise g(t) due to stochasticity of the adaptation dynamics, and (iv) the deterministic feedback of the neuron's spike train w(t) due to the deterministic part of the adaptation (see Fig. 2). In Eq. (1), the parameter b determines the strength of adaptation and s 2 is set by the duration of the action potential t AP relative to the mean ISI (Methods, Eq. (41)).
To study the effect of the two different kinds of noise, we focused on two limit cases: In the limit of infinitely many channels, the adapting PIF model is only driven by white noise. In this case, Eq. (1) reads We call this case deterministic adaptation. A dimensionality analysis shows that the ISI statistics are completely determined by the quantities t w m and t w D if one assumes b and t AP as constants (see Methods). Thus, decreasing m and D by some factor and simultaneously increasing t w by the same factor (yielding a decreased firing rate r!m) would, for instance, not alter the statistical properties of the model.
In the opposite limit, we considered only the stochasticity of the adaptation current but not the white noise. Setting D~0 we obtain We call this case (and the corresponding model based on individual adaptation channels) stochastic adaptation. As shown in Methods, this case is determined by the quantities t w m and t w =N a (assuming b and t AP as constants). For instance, decreasing m and increasing t w and N a by the same factor (and thereby lowering the firing rate r) would again result in an equivalent model with the same statistical properties.
The fraction of open channels W~N op N a performs a random walk with discontinuous jumps. The direction of jumps depends on Figure 1. Integrate-and-fire dynamics with adaptation channels. A Channel model: a population of N a independent voltage-gated ion channels, which can be either in an open or a closed state, mediate an adaptation current through a neuron's membrane. B PIF model: Subthreshold dynamics of the membrane potential V (bottom). The variable V (measured in units of V th ) is reset to a value V reset~0 after crossing the threshold at V~V th . Action potentials are not generated explicitely. Instead, the effect of an action potential is captured by the activation function w ? (t), which is set to one in a short time window of 1ms following each threshold crossing of the IF model (middle panel). The adaptation current is proportional to the fraction of open channels W (top panel). The sample traces were obtained from a simulation of Eq. (35) with N a~1 000 channels, white noise intensity D~0:01V 2 th ms, adaptation time constant t w~1 00ms, base current m~0:4V th =ms and maximal adaptation current b~3V th =ms. C The time-dependent firing rate (top) in response to a step stimulus (bottom) is independent of the source of noise (stochastic adaptation -solid line, deterministic adaptation plus white noise -dashed line). The gray line shows the theory given by Eq. (55). doi:10.1371/journal.pcbi.1001026.g001 the presence of spikes, which in turn is affected by W ( Fig. 2A). The diffusion approximation of W and the separation of deterministic and stochastic components are illustrated in Fig. 2B and 2C, respectively. Although the increments of the continuous diffusion process have the same (Gaussian) statistics as the original discontinuous process on a time interval larger than the mean ISI, the short-time statistics is rather different ( Fig. 2A,B). Therefore, it is not obvious whether the diffusion approximation yields a good approximation to the ISI statistics, and in particular, how this approximation depends on the number of channels N a and the adaptation time constant t w . To clarify this issue, we performed both simulations based on individual channels (''channel model'') and simulations of Eq. (1) (''diffusion model''). It turned out, that the diffusion approximation yields a fairly accurate approximation for the shape of the ISI density, the coefficient of variation and the serial ISI correlations even for small channel populations. However, significant deviation were found for higher-order statistics like the skewness and kurtosis of the ISIH (see next section).

Interspike interval statistics of the adapting PIF model
The calculation of the ISI statistics (histogram and serial correlations) of the PIF model with noise and spike-frequency adaptation is generally a hard theoretical problem. Here we put forward several novel approximations for the simple limit cases Eq.
(2) and Eq. (3). For typical adaptation time constants that are much larger than the mean ISI we found the ISI histogram in the case of pure white noise (N a ??, Eq. (2)) mapping the model to one without adaptation and renormalized base currentm m (Methods, Eq. (52)). This corresponds to a mean-adaptation approximation [18,[43][44][45][46][47], because the adaptation variable w(t) is time-averaged by the linear filter dynamics in Eq. (2b) for t w much large than the mean ISI (rt w &1). However, this approximation cannot account for ISI correlations, because any correlations between ISIs are eliminated in the limit t w ?? -in fact, the reduced model is a renewal model. For this reason, we developed a novel technique to calculate serial correlations for a PIF neuron with adaptation and white noise driving, which is valid for any time constant t w (see Methods).
In the opposite limit of only adaptation fluctuations (D~0, Eq. (3)), we could calculate analytically the ISI histogram, the skewness and kurtosis of ISIs as well as the ISI serial correlations by mapping the problem to one without an adaptation variable but a colored noiseg g(t) with renormalized parameters. Specifically, the IF dynamics for only adaptation channel noise reduces to where the effective parameters are scaled by a common scaling factor:m As before, a spike is fired whenever V reaches V th~1 , whereupon the voltage is reset to V reset~0 . We call this model (Eq. (4)-(6)) the colored noise approximation. For the perfect integrate- and-fire model driven by a weak colored noise, i.e. for the model described by Eq. (4), analytical expressions for the ISI density and the serial correlation coefficient are known [48]. In addition to that, we derived novel analytical expressions for the skewness and kurtosis of the ISIs (see Methods).
Interestingly, the scaling factor in Eq. (6) has a concrete meaning in terms of spike-frequency adaptation: l coincides with the degree of adaptation in response to a step increase of the base current (see Methods, Eq. (56)).
ISI density. Fig. 3A shows ISI histograms (ISI densities) for the case of deterministic adaptation. We found, that the ISI densities can be well described by inverse Gaussian probability densities with mean V th =m m given by Eq. (64) (see Methods). In the opposite case of stochastic adaptation, the ISI variability solely depends on the number of slow adaptation channels (Fig. 3B). For a small channel population (N a~1 00) the discreteness of the adaptation W still appears in the ISIH as single peaks that cannot be averaged out. This is related to realizations of the channel noise for which the fraction of open channels does not change during the ISI; realizations for which the fraction changes at least once lead to the continuous part of the ISI density. In contrast, the diffusion model yields a purely continuous curve, that looks like a smoothed version of the ISIH of the model with channel noise. As N a increases, the discrete peaks in the latter become more and more dense and insignificant, and the ISIH of the channel model is well approximated by the diffusion model. Furthermore, the theory for the colored noise approximation, Eq. (4), coincides well with the diffusion model, Eq. (1), and hence for sufficiently large N a it also fits well the ISIH of the channel model.
One central claim of this paper is that ISI histograms of neurons, for which slow channel noise dominates the ISI variability, cannot be described by an inverse Gaussian (IG) distribution in contrast to cases where fast fluctuations dominate. We recall that the IG distribution yields the ISI histogram for a PIF model driven by white noise without any (deterministic or stochastic) adaptation, so a priori we cannot expect that this density fits any of the cases we consider here. However, as mentioned above, the ISI density can be captured by an IG for deterministic adaptation (Fig. 4A,C). In fact, the main effect of a slow adaptation is to reduce statically the mean input current which is reflected in our approximation by going from m tom m. Slight deviations of the simulated ISI histogram from the IG can be seen for large intervals where the simulated density displays a stronger decay than the IG (Fig. 4C). With adaptation, large intervals are prevented because for large times (after the last spike) the inhibitory effect of the adaptation current subsides -a feature that is not present in the static approximation for the reduced base current which was made above. Nevertheless, the deviations are small and will be hardly visible when comparing the IG density to the histogram of limited experimental data sets. By contrast, the ISI histogram in the case of stochastic adaptation as illustrated in Fig. 4B and D, possesses a much stronger peak and decays slower at large ISIs than the IG with the same mean and variance of the ISI; comparison to the IG density that has the same mean and mode or the same mode and CV gave comparably bad fits (data not shown). Instead, the colored noise approximation as outline above, describes the simulation data fairly well.
This suggests that both cases -deterministic and stochastic adaptation -might be distinguishable from the shape of the ISI histograms even if mean and CV of the ISIs are comparable as for the data in Fig. 4. To this end, we introduced new measures a s and a e based on the skewness and the kurtosis (excess) of the ISI distribution that are exactly unity for an IG distribution (see Methods, Eq. (61) and (62)). Indeed, Fig. 5 reveals that deterministic and stochastic adaptation are well separated with respect to the rescaled skewness and kurtosis a s and a e . In particular, these quantities are clearly larger than unity for stochastic adaptation meaning that the ISI density is more skewed and more peaked compared to an IG, which confirms our previous observations. On the other hand, for deterministic adaptation, a s and a e are smaller than unity in accordance with our previous observations that the tail of the ISI density decays slightly faster than an IG density.
The rescaled kurtosis reveals also differences between the channel and the diffusion model. In Fig. 6A, the CV still matches almost perfectly for both models even at extremely small channel numbers, where the Gaussian approximation is expected to fail. This is also remarkable in the light of the discrete structure of the ISIH for small channel numbers (cf. Fig. 3 for N a~1 00). However, in Fig. 6B it becomes apparent that the two models differ with respect to higher-order measures as a e ; for increasing numbers of channels the differences decrease. Fig. 5 and 6 also support the colored noise approximation, which describes the diffusion model quite accurately. This suggests, that the heavy-tailed and pronouncedly peaked ISIH in the case of stochastic adaptation can be simply explained by the effect of a long-correlated, colored noise. It is known that for the related leaky IF model such correlations result in ISIHs with a large kurtosis [49]. To examine the role of long-range temporal correlations in shaping the ISI density we analyzed the dependence of the ISI statistics on the time-scale separation between adaptation time constant and mean ISI. This can be quantified by the ratio t w =ST i T~rt w , where r and ST i T denote the stationary firing rate and the mean ISI, respectively.
In the case of stochastic adaptation, we obtained analytical expressions for a s and a e using the colored noise approximation for weak noise (see Methods). For the following discussion, it suffices to consider the zeroth order of the weak-noise expansion, which is given by and  These expressions only depend on the non-dimensional parameter d~ST i T=t t~(lt w r) {1 , i.e. on the product of rescaled adaptation time constant and firing rate. From Eq. (7) and (8), it can be shown that a s w1 and a e w1 for rt w w0 and that both quantities converge to unity in the limit t w r?0. For t w much larger than the mean ISI, i.e. when rt w is large, the leading orders saturate at a (0) s~2 and a (0) e~2 4=5. These predictions are confirmed by simulations using different t w at a fixed value of m~0:4 (r~100Hz, Fig. 7), thereby varying the time scale separation rt w (l remains constant). In particular, both rescaled skewness and kurtosis are larger than unity and increase strongly at moderate fluctuations (N a~2 00, Fig. 7A,B). The increase is more pronounced for the diffusion model compared to the channel model. At large rt w the simulation data deviate from the first-order approximation, because higher-order terms in the small noise expansion cannot be neglected. The agreement becomes better when the number of channels becomes larger (see Fig. 7C,D for N a~5 00). These observations were qualitatively confirmed by corresponding simulations at a smaller base current (m~0:04) leading to a lower firing rate of r~10Hz (data not shown). In particular, quantitatively similar curves were obtained when at r~10Hz the number of channels was increased to N a~2 000 in order to maintain the same effective noise level. It should be noted, that in the case of low firing rates the weak-noise expansion might become infeasible if the number of channels is too small. This is, for example, the case for N a~2 00 at r~10Hz, for which the small parameter (Eq. (68)) becomes larger than unity. For deterministic adaptation a s and a e approach unity for rt w &1 as predicted by the mean-adaptation approximation, for which the ISIH is given by the IG (Eq. (64)). In the opposite limit of small rt w the parameters a s and a e also approach unity. This is intuitively clear, because for t w %ST i T the adaptation w(t) decays quickly to zero after each spike. Hence, the base current is almost always equal to m except for a very short period after a spike where the driving is m{bw(t). Put differently, the dynamics can be approximated by a PIF model with a constant driving m and an effective reset value V reset &{bt AP v0. In this case, the ISIs are again distributed according to the IG statistics.
In the intermediate range, where the time scale of the adaptation is of the same order as the mean ISI, a pronounced minimum of a s and a e is observed in the case of deterministic adaptation. This is due to the decay of adaptation at such a rate that large ISIs are suppressed. As a consequence, the tail of the ISIH decays faster and the ISIH becomes less skewed compared to an IG. The same qualitative behavior was verified in simulations at a lower firing rate r~10Hz (data not shown).
ISI correlations. Another clear distinction between stochastic and deterministic adaptation is revealed by the correlations between ISIs. In several modeling studies it has been found that negative feedback mechanisms like adaptation currents in the presence of white noise give rise to negative correlations between adjacent ISIs [28,29]. However, a theoretical explanation of this effect has not been provided yet. Therefore we developed a theory based on the dynamics close to the deterministic limit cycle of the adaptation dynamics (see Methods). This dynamics can account well for the correlations between ISIs in the case of deterministic adaptation. Specifically, the serial correlation coefficient (SCC) for two ISIs with lag n (see Eq. (63)) is given by the geometric sequence where and Noting that 0vnv1 and qv1, we find that the prefactor in front of the term (nq) n{1 is negative. Thus, correlations at odd lags are always negative, whereas ISIs with even lag are anti-correlated only if qw0. If on the other hand qv0, ISIs with even lags exhibit a positive SCC, giving rise to oscillations of r n . Both cases, purely negative correlations with an exponential decay and oscillating correlations, are indeed observed for deterministic adaptation (Fig. 8A). The slight deviation of the theoretical prediction is due to the short rise phase of w following each spike in the simulations instead of the instantaneous increase of w assumed in the derivation of Eq. (9).
In striking contrast, the case of stochastic adaptation yields positive ISI correlations with a slow exponential decay (Fig. 8B) as predicted by the theory, Eq. (72). From the formula it becomes evident, that the decay constant is to first-order given by the ratio d {1~t t=ST i T of effective correlation timet t~lt and mean ISI. The good agreement of the colored noise theory suggests, that adaptation noise effectively acts as a colored noise that slowly modulates the ISIs. It is known, that systems with a slow stochastic driving exhibit positive ISI correlations [48][49][50]. In fact, in the absence of additional fast fluctuations the ISIs are strongly correlated with the slow noise, which retains a memory of previous ISIs. For instance, if T i vSTT due to a large, positive fluctuation of g, this will on average cause a likewise small T iz1 vSTT, because the slow dynamics of g tends to persist at positive values in the course of several subsequent ISIs.
ISI correlations are strongly governed by the time scales in the system. We therefore investigated the role of the time scale separation parameter rt w~tw =ST i T on the serial correlations of adjacent ISIs (Fig. 9). For deterministic adaptation we found that adjacent ISIs become most anti-correlated at a finite value of rt w close to unity. In the limits rt w ?0 and rt w ??, however, r 1 vanishes as predicted by the theory and as observed in previous studies [29,34]. This is intuitively clear, because in the first limit the adaptation variable w(t) cannot accumulate by subsequent spikes, and hence no memory of previous ISIs is retained; in the latter limit because w(t) converges to its (constant) mean value for rt w ??. Interestingly, the ISI correlations seem to be almost independent from the noise intensity D if D is not too large (compare Fig. 9A and B for two different values of D). This insensitivity of the correlation coefficient to the noise intensity could be anticipated from the analytical theory (see Methods), in which the noise dependent term cancels out in the ratio of ISI covariance and ISI variance.
For stochastic adaptation, the positive correlations become strongest for t w much larger than the mean ISI ( Fig. 9), i.e. rt w &1. The channel and diffusion model agree generally quite well, except for very small and very large rt w . In the latter case, the Gaussian approximation becomes worse, because the opening and closing events of the channels are extremely rare. The expected number of channel transitions in a time window of length t is N(t)~2(1{SwT)SwTN a t=t w with SwT~t AP =ST i T. For instance, taking the extreme case t w~1 0 4 ms at the standard parameters N a~2 00, SwT~0:1 one would on average observe only 0:036 transitions on the time scale of a single ISI (t*10ms). On this time scale the fraction of open channels can hardly be approximated by a Gaussian process. As expected, we obtained a better agreement between channel and diffusion model at large rt w by increasing the number of channels (Fig. 9B). The decrease of r 1 for the diffusion model and the colored noise approximation at very large rt w might be due to the fact that the ISI variance grows faster with rt w than the covariance ST i T iz1 T{ST i T 2 , thus the correlation coefficient is suppressed by the variance. A similar effect has been observed for the LIF model [49].
Mixed case of fast and slow noises. So far, we found that the two limit cases of the adapting PIF model can be well distinguished by the values of the shape parameters a s and a e relative to unity and the correlation coefficient r 1 relative to zero. Do these quantities also allow for an unambiguous distinction of the dominating source of noise in the more realistic case where both kinds of noise are present? To answer this question, we performed simulations of the adapting PIF model for a fixed intensity of the white noise (''fast fluctuations'') but different sizes of the population of adaptation channels. Thereby, we could vary the ratio of the two different types of noise. Note, that the mean adaptation current is kept constant in our setting. This can be realized by scaling the single channel conductance or the membrane area with 1=N a (see Methods).
For small channel numbers, i.e. large channel noise, we observed both large values of the rescaled kurtosis a e w1 and a positive serial correlation coefficient of adjacent ISIs (Figs. 10A,B, white region) indicating the strong impact of the colored noise effect. As expected, at the other end of large channel population sizes the pure white noise case could be recovered. In between, we found a critical channel number at which both the rescaled kurtosis crossed the line a e~1 and the serial correlation coefficient changed its sign. This simultaneous change suggests, that below the critical channel number the ISI statistics was dominated by slow adaptation channel noise, whereas above this critical size it was dominated by the white noise input (gray-shaded region).

Effects of a stochastic adaptation current on the ISI statistics of a Hodgkin-Huxley type model
We investigated whether our theoretical predictions based on a simple integrate-and-fire model are robust with respect to a more detailed model of the Hodgkin-Huxley type. To this end, we performed simulations of the conductance-based Traub-Miles model with a M-type adaptation current [51]. As in the previous model we separately considered the two cases of white noise input and a slow M-type channel noise to get an intuition of the individual effects on the ISI statistics. Fig. 11 demonstrates that the ISI histograms show essentially the same features as in the PIF model: in the case of white noise input the shape of the ISIH could be well approximated by an inverse Gaussian distribution which was uniquely determined by the firing rate and the CV. In the case of a stochastic M-type current there is a strong disagreement between the ISIH and an inverse Gaussian with the same rate and CV. In particular, ISIHs exhibited again a sharper peak compared to the relatively broad inverse Gaussian.
The different ISI statistics for the case of deterministic and stochastic adaptation are analyzed more closely in Fig. 12. As in the PIF model (cf. Fig. 5) the rescaled skewness and kurtosis are significantly smaller for white noise than for adaptation noise in a wide range of CVs (Fig. 12A,B). This is in accordance with the pronounced peak of the ISIH in the case of stochastic adaptation (Fig. 11B). However, the values are not strictly separated by a s=e~1 as in the PIF model. This discrepancy is not surprising,   given that the Traub-Miles dynamics with constant input and white noise driving does not exactly yield an inverse Gaussian ISI density but only an approximate one. Importantly, however, the rescaled kurtosis a e quickly saturates at a finite value in the large t w limit (albeit not at unity, Fig. 12C). This is markedly different from the case of stochastic adaptation. In this case, the rescaled kurtosis increases strongly as it was observed for the PIF model. In a similar manner, the rescaled skewness also showed this distinct behavior for stochastic vs. deterministic adaptation, although the increase of the rescaled skewness was not as strong as for the rescaled kurtosis (data not shown).
A clear distinction between both cases appears in the serial correlations of ISIs (Fig. 12D). Similar as in the PIF model, the case of deterministic adaptation is characterized by negative ISI correlations at lag one, which are strongest at an intermediate time scale t w . Furthermore, the case of stochastic adaptation exhibits positive correlation coefficients r 1 , which show a maximum at an intermediate value of t w . This is also in line with the PIF model. The correlations decay rapidly with the lag for deterministic adaptation (Fig. 12E) and decay exponentially for stochastic adaptation (Fig. 12F). As in the PIF model, the exponential decay is slower for large time constants t w .
Finally, we inspected the case in which both white noise and slow adaptation noise is present (Fig. 12G,H). As in Fig. 10 for the PIF model, we fixed the noise intensity of the white noise and varied the number of adaptation channels N a . In the Traub-Miles model one finds qualitatively similar curves as in the PIF model. In particular, the serial correlation coefficient at lag one, shows a transition from positive to negative ISI correlations at a certain number of adaptation channels (Fig. 12H). As for the PIF model, this value can be used to define two regimes -one dominated by adaptation noise (white region) and another one dominated by white noise (gray-shaded region). In the adaptation-noise dominated regime the parameter a e is larger than in the white-noise dominated regime (Fig. 12G).
The observation that key features of the ISI statistics in the presence of a stochastic adaptation current seem to be conserved across different models suggests a common mechanism underlying these features. As we saw, this mechanism is based upon the fact that a stochastic adaptation current can be effectively described by an independent colored noise. The long-range temporal correlations of this noise naturally yield positive ISI correlations and a slow modulation of the instantaneous spiking frequency. The latter typically involves a large kurtosis due to the increased accumulation of both short and long ISIs. A significant amount of colored noise can effect the kurtosis and the ISI correlations so strongly, that details of the spike generation seem to be of minor importance. Thus, it becomes plausible that the spiking statistics of a rather complex neuron model could be explained by a simple integrate-and-fire model including a stochastic adaptation current.

Discussion
In this paper, we have studied how a noisy adaptation current shapes the ISI histogram and the correlations between ISIs. In particular, we have compared the case of pure stochastic adaptation with the case of a deterministic adaptation current and an additional white noise current. Using both a perfect IF model that is amenable to analytical calculations and a more detailed Hodgkin-Huxley type model, we found large differences in the ISI statistics depending on whether noise was mediated by the adaptation current or originated from other noise sources with fast dynamics. As regards the ISI histogram, stochasticity in the adaptation leads to pronounced peaks and a heavy tail compared to the case of deterministic adaptation, for which the ISI density is close to an inverse Gaussian. To quantify the shape of ISI histograms we proposed two measures that allow for a simple comparison with an inverse Gaussian probability density that has the same mean and variance. The first one is a rescaled skewness (involving the third ISI cumulant); the second is a rescaled kurtosis (involving the fourth ISI cumulant). Both quantities possess the property that they assume unity for an inverse Gaussian distribution. If they are larger than unity as in the case of stochastic adaptation the ISI density is more skewed or respectively has a sharper peak and a heavier tail than an inverse Gaussian density with the same variance. If these measures are smaller than one, the ISI histogram tends to be more Gaussian like. Most strikingly, we found that for a stochastic adaptation current the rescaled skewness and kurtosis strongly increase when the time scale separation of adaptation and spiking becomes large (t w =ST i T&1). By contrast, for a deterministic adaptation current the rescaled kurtosis saturates close to one in this limit.
Another pronounced difference arises in the ISI correlations. For a deterministic adaptation current and a white noise driving one observes short-range anti-correlations between ISIs as reported previously (e.g. [29]). In contrast, with slow adaptation noise ISIs exhibit long-range positive correlations. In the presence of both types of noise, the serial correlation coefficient changes continuously from positive to negative values when the ratio of white noise to adaptation noise is increased. The two domains might be useful in determining the dominating source of noise from a neural spike train.
Interestingly, the perfect integrate-and-fire model augmented with an adaptation mechanism predicted all the features seen in the spiking statistics of the Traub-Miles model with stochastic adaptation and/or white noise input. This indicates the generality and robustness of our findings. It also justified the use of the adapting PIF model as a minimal model for a repetitive firing neuron with spike-frequency adaptation. It seems, that in the suprathreshold regime the details of the spike generator are of minor importance compared to the influence of adaptation and slow noise. By means of the PIF model one can theoretically understand the underlying mechanism leading to the large kurtosis and the positive ISI correlations in the case of stochastic adaptation. This rests upon the fact that slow adaptation noise effectively acts as an independent colored noise with a large correlation time. One can think of the colored noise as a slow external process that slowly modulates the instantaneous firing rate or, equivalently, slowly changes the ISIs in the sequence. Such a sequence of many short ISIs in a row and a few long ISIs gives rise to a large skewness and kurtosis and positive serial correlations. In previous works, slow processes which cause positive ISI correlations were often assumed to originate in the external stimulus [49,50,52]. Here, we have shown that an intrinsic process, i.e. the fluctuations associated with the stochasticity of adaptation, yields likewise positive ISI correlations. Our finding also provides an alternative explanation of positive ISI correlations in experimental studies [30,53]. Moreover, in vivo recordings from a looming-sensitive interneuron in the locust optic lobe have revealed both positive correlations at large firing rates and negative correlations at low firing rates [23]. Because this neuron exhibits pronounced spike-frequency adaptation an intriguingly simple explanation for these observations would be the presence of both fast noise and stochastic adaptation (corresponding to our mixed case). In this case, a large firing rate could indeed lead to a large effective correlation time of the noise associated to the adaptation mechanism and thus to positive ISI correlations.
Spike-frequency adaptation has been ascribed to different mechanisms (see e.g. [18]), involving for instance, calciumdependent potassium currents I AHP [42], slow voltage-dependent M-type currents I M [41,42] and slow recovery from inactivation of sodium currents [54]. Here, we chose the M-current as an example to illustrate the emergence of noise in the adaptation mechanism. In this specific case, it was the finite number of M-type potassium channels that gave rise to slow channel noise. For the other commonly studied adaptation mechanism, the I AHP [18,23,24,26,29,51], we have to deal with two possible sources of noise: the finite number of potassium channels N K and fluctuations of the local Ca 2z concentration c. Proceeding in a similar fashion as for I M , we would obtain I AHP~ g g AHP K(V {E K ), with the fraction K of open potassium channels, obeying Here, the Gaussian white noises j k and j c approximately represent the channel noise and the concentration fluctuations due to stochastic removal of calcium, respectively. The calcium gating is characterized by the steady-state activation k ? (c). For simplicity, the increase of calcium D caused by an action potential is assumed to be deterministic. Importantly, however, the channel dynamics is fast compared to the slow removal of calcium, i.e. t k %t c . Following [18] the open probability of the potassium channels k:SKT adiabatically adjusts to c (i.e. k&k ? (c)) and the relationship is roughly linear (i.e. k ? (c)!c). Thus, we have I AHP !(czg k )(V {E K ), where the ''channel noise'' g k possesses a correlation time t k . If this correlation time is much smaller than the mean ISI, the channel noise can be approximately treated as a white noise. But this means, that a PIF neuron with a calcium-dependent I AHP instead of I M can likewise be approximated by Eq. (1): the fast channel noise can be included into the white noise term ffiffiffiffiffiffi ffi 2D p j(t) and the slow fluctuations of the calcium concentration assume the role of the slow adaptation noise g(t). Approximating I AHP again by a voltage-independent current, the PIF model with I AHP would read These equations can indeed be put into the form of Eq. (1) by splitting the deterministic and the noise part of c. This illustrates that the main results derived in this paper are not specific to a certain adaptation current, but apply quite generally to any noise associated to the slow dynamics of adaptation.
The adaptation currents I M and I AHP have been distinguished with respect to their ability to synchronize coupled neurons [51] and regarding the influence on neural coding [55]. The difference consists in whether the current is activated solely by spikes as in the case of I AHP or whether it is also activated by subthreshold voltages as for I M . For the sake of clarity, we have set the activation function w ? (t) of the M-type adaptation current in the PIF model equal to zero at subthreshold voltages, i.e. between spikes (Eq. (25)). Thus, the adaptation current in the PIF model, unlike the M-current, was only activated during action potentials. It is, however, easy to show that the results of this paper are unchanged if subthreshold activation is allowed. For simplicity, let us consider the extension that in-between spikes the steady-state activation function w ? is equal to the value v1, i.e. instead of Eq. (25), (26) (see Methods) we have This only increases the mean adaptation to SwT~Sw ? (t)T~z(1{ )rt AP (cf. Eq. (37)). Similarly, the variance changes according to s 2~S wT{SwT 2 (cf. Eq. (40) in Methods). As a result, the effective base current is now given bỹ with the new scaling factor The colored noise approximation can be carried out in an analogous manner yielding the same result Eq. (4) (again with t t~lt w ,s s 2~l s 2 , but the new scaling factor l, Eq. (19)). Thus, it can be expected that in the presence of subthreshold activation of w(t) the colored noise effect (i.e. pronounced peak of ISIH, positive ISI correlations) in the case of stochastic adaptation is preserved. Furthermore, l still serves as the degree of adaptation, i.e. the ratio of steady-state to initial gain when a step current is applied.
The analytical calculation of higher-order statistics in the presence of adaptation is a fundamental theoretical problem, which has been largely ignored so far (for a recent exception see [37]). Here, we succeeded to provide explicit expressions for the ISI histogram and their serial correlations for both white noise driving and noise in the adaptation dynamics. This was achieved by analyzing a spike generator and a channel model that are as simple as possible. There are certainly a lot of details that can be modeled in a more realistic way. For instance, it is known that the M-channel kinetics is governed by several time scales and more than two internal states [56]. Furthermore, channels might not be strictly independent, but channel clusters might exhibit cooperative behavior [57]. The latter case, would actually increase the level of channel noise compared to the case of independent channels, i.e. cooperativity would contribute to stochastic adaptation.
For many neurons physiological details like the number of ion channels are hard to obtain directly from experiments. Instead given the spike train statistics of a neuron, our study could be useful to judge whether M-channels or other adaptation mechanisms could potentially contribute to the neuronal variability. Furthermore, it is not impossible to think of experiments, in which the number of adaptation channels is reduced (e.g. by the mild application of a channel blocker) and thus the effects of stochastic adaptation is affected in a controlled way. Another possibility to test our predictions would be to vary the firing rate of the neuron by increasing or decreasing the input current. In this way, the time scale of spiking would change relative to the time scale of adaptation and, thus, the colored noise effect of adaptation noise could be enhanced or attenuated, respectively.
Channel noise can crucially influence neural firing especially in the absence of synaptic input [58,59]. This could be particularly relevant for the irregular discharge patterns of certain receptor cells. So far, channel noise has been studied mostly in the context of stochastic Na z and K z channel gating involved in the spike generation itself. These channels are considered to be fast. Because we were mainly interested in the effect of slow adaptation channels compared to fast fluctuations resulting from fast ion channels or synaptic activity, we lumped all fast noise sources into an unspecified additive white noise source. This is certainly a simplification; e.g. it has been shown in experiments that voltage noise due to Na z channels depends on the mean voltage itself [39,40]. More detailed models of the various sources of noise are worth the efforts in future investigations. However, we do not expect that such sophisticated models would change our results qualitatively, because they mainly hinge upon the presence or absence of long time scales.
Realistic numbers of M-type channels per neuron are difficult to estimate and the numbers used in this paper must be seen as a tuning parameter for the channel noise intensity. Channel densities of the M-type have been estimated to be of the order of one functional channel per 4mm 2 [60]. Assuming a spherical cell with a diameter of 10mm one obtains of the order of 100 channels. Thus, the channel numbers used in this study (N a~2 0-1000) seem to be reasonable; and hence the M-current could be a potential source of fluctuations.
The diffusion approximation for the stochastic dynamics of ion channel populations (also known as Langevin or Gaussian approximation) has been studied by several authors [38,[61][62][63][64] (see also [65] in the context of chemically reacting particles). Here, we have shown how one can map the stochastic dynamics of a population of ion channels with negative feedback to the macroscopic current dynamics plus an additive colored noise (see [38] for a related treatment). In other words, the dynamics could be reduced to an analytically accessible Langevin equation for voltage and adaptation. In particular, we investigated the effect of the diffusion approximation on the statistics of interspike intervals and found a fairly good agreement with the channel model, despite the small number of channels. This seems surprising, given that for a typical parameter set -N a~2 00, open probability SwT*0:1 and t w~1 00ms -one expects only 1:8 closing transitions (between spikes) and 0:2 open transitions (during action potentials) per millisecond. Apparently, on the time scale of 1 ms the number of transition events is not Gaussian distributed. However, we found that the main effect of the channel noise consists in a slow modulation of the instantaneous firing rate on the large time scale oft t~lt, whereas high-frequency components of the noise are of minor importance. Thus, on relevant time scales of the order of 10-100 ms the average number of transitions is much larger and a Gaussian approximation seems to be reasonable.
Spike-frequency adaptation has been commonly studied with regard to its mean effect on the firing rate [18,[43][44][45][46]. It has been shown that these effects can be exhaustingly analyzed using a universal firing rate model [18]. In this paper, however, it became evident that higher-order statistics and fluctuation effects may differ and may be used to distinguish different kinds of noise sources.

Model for the stochastic adaptation current
To analyze the slow, voltage-dependent adaptation channels in a simple setup we consider a population of N a independent ion channels that reside in an open or a closed state. For each channel, we thus have the simple reaction kinetics shown in Fig. 13A.
For N a channels, one can either perform N a independent simulations of one two-state process or one simulation with N a states where the number of open channels N op can be increased or decreased by one: Here, N cl (t)~N a {N op (t) denotes the number of closed channels. The rates a w and b w for the transitions between the closed and open state can be related to the (voltage-dependent) kinetics of a typical gating variable by choosing Therein, w ? (V ) is the steady-state open probability of a single channel when the membrane potential is clamped at V and t w sets the time scale of the channel kinetics. Note, that both w ? (V ) and t w are accessible from experiments. The master equation for the open probability of the two-state model reads or after insertion of Eq. (21) which follows the common scheme for gating variables with voltage-dependent kinetics in Hodgkin-Huxley type models. We will use the function w ? in two different versions.
In the conductance-based Traub-Miles model, we will use the kinetic scheme with the rate functions given by Inserting this relation into Eq. (21), both rate functions show a quasi sigmoidal dependence on the voltage (see Fig. 13B), such that a w is appreciably different from zero only during the action potential (i.e. for V w{50mV) whereas b w is only finite in the opposite range of a subthreshold membrane voltages.
For the integrate-and-fire (IF) model we employ a simplified variant that distinguishes only between two values of w ? : one attained when the voltage is subthreshold and another one for the duration of the action potential t AP . Specifically, at a spike event of the IF model, w ? is set to one for a duration of t AP~1 ms and otherwise it is set to zero (Fig. 1B). Thus, w ? can be expressed as a function of time t and the last spike time t i Ã ƒt: This function and the set of spike times ft i g are related by where h( : ) denotes the Heaviside function. Hence, in the simplified model, the dependence of the rates on the membrane voltage is substituted by an explicit dependence on time t and the most recent spike time t i Ã ƒt: In the definition of w ? (t), Eq. (25), we require that the duration of the pulse t AP is much smaller than the mean ISI, so that overlaps of two subsequent pulses are unlikely. Note, that for simplicity channels can be activated only during spikes (a w~1 =t w ), but not in between spikes (a w~0 ). However, one can show that the following results are not changed qualitatively if one allows for subthreshold activation (w ? (t)w0), as observed for M-type potassium currents (see Discussion).
Taking the open probability exactly as the fraction of open channels amounts to taking the limit of an infinite population of channels. In the Traub-Miles model or the integrate-and-fire model with deterministic adaptation, this corresponds to adding an adaptation current of the form [18] In this equation, g g a denotes the maximal conductance (per unit membrane area) and E a constitutes the reversal potential of the adaptation current.
For a finite channel population, the fraction of open channels is given by the stochastic quantity W (t):N op (t) N a , where N op (t) evolves according to the kinetic scheme Eq. (20). In contrast to Eq. (28), this gives rise to a stochastic adaptation current When the channel number is varied, we assume that the maximal conductance per unit membrane area g g a remains constant. Thus, a change of the channel number can be realized either by a variation of the total membrane area or by a change of the channel density in conjunction with a simultaneous scaling of the single channel conductance. Such procedures enable a change of the amount of channel noise, without changing the mean current per unit membrane area.
Diffusion approximation. The channel model, Eq. (20), can be approximated by a single Langevin equation for N op (t) if the number of transitions in a sufficiently large time interval, which is still smaller than the decay time t w , is treated as a Gaussian random variable [65]. Under this condition N op obeys where s 2 (W ,w ? (t)) is given by Furthermore, a separation of the adaptation into a deterministic and a stochastic part, W~wzg, will be useful for the interpretation of our results. In these new variables Eq. (31) can be rewritten as two equations: Note that Eq. (1), which will be our final diffusion model, also involves the additive-noise approximation presented below.

Perfect integrate-and-fire model with adaptation
The perfect integrate-and-fire (PIF) model [2] constitutes a minimal model for a neuron possessing a stable limit cycle. In this model the subthreshold voltage is determined by the equation where m~I 0 =C m and {bW are proportional to the base current I 0 and adaptation current I a (t), respectively; C m is the membrane capacitance and the scaling factor for the adaptation current reads b~ g g a (E 0 {E a )=C m . Here, we used an effective-time-constant approximation [66], where we substituted in Eq. (29) V by the average voltage E 0~S VT to obtain a voltage-independent adaptation current [18]. The last term in Eq. (35) represents fast Gaussian input fluctuations of intensity D and correlation function Sj(t)j(t')T~d(t{t') (here and in the following, the angular brackets denote an ensemble average). The model Eq. (35) is complemented by the fire-and-reset rule: upon reaching the threshold V th a spike is elicited and V is reset to V reset , V reset vV th . Because Eq. (35) is invariant with respect to a constant shift in V , we can choose the reset value as the origin, i.e. V reset~0 . The threshold crossing events define the spike times t i , i~0,1,…of the PIF model.
It is a feature of the PIF model that the firing rate is directly proportional to a constant driving current and independent of the noise. However, even for the slowly varying driving current m{bW (t) in Eq. (35), one can define an instantaneous firing rate which can be seen as the slowly varying firing rate that is obtained by averaging the spike train over the time scale t w of the adaptation. Averaging over time scales much larger than t w yields a constant firing rate r, because the process is stationary. Thus, the (stationary) firing rate can be obtained by averaging Eq. (36), which gives r~m{bSwT ð Þ =V th . On the other hand, the firing rate is related to the mean fraction of open channels by averaging Eq. (33): Solving for r yields where l~(1zbt AP =V th ) {1 .
Additive-noise approximation. The noise term in Eq. (31) or (34) is multiplicative, because both arguments of the noise strength s 2 W ,w ? (t) ð Þdepend on the dynamical variables. For slow adaptation, rt w &1, and large N a , however, the relative fluctuations of s 2 (W ,w ? (t)) are small and, hence, s 2 can be treated as constant (additive-noise approximation). This can be seen as follows: The rapid changes of s 2 (W ,w ? (t)) due to the switchings of w ? (t) between 0 and 1 on the time scale of the mean ISI are averaged out by the linear filter dynamics Eq. (34) if rt w &1. Thus, in Eq. (32) w ? (t) can be replaced by its local average t AP r r(t). Because of Eq. (36), s 2 can therefore be written as a quadratic function of W . As explained below (Sec. ''Colored noise approximation''), the variance of W is ls 2 N a . Thus, using SW T~SwT and SW 2 T&lSs 2 T=N a zSwT 2 the mean of s 2 can be self-consistently determined. It can be shown that &SwT{SwT 2 ð40Þ where D~bt AP =V th and r is given by Eq. (38). In the second line we neglected the term that depends on N {1 a , because for large N a it is much smaller than unity. Similarly, the variance of s 2 can be calculated neglecting terms that are third and forth order in W . As a result, we obtain Var s 2 Â Ã~( 1{D) 2 4(1zD)N a SwTzO SwT 2 À Á .
The relative fluctuations of the noise intensity are given by ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Var s 2 ½ p . Ss 2 T. Inserting the expressions for mean and variance of s 2 to first order in SwT, we find that the additive-noise approximation s(W ,w ? (t))&Ss 2 T is roughly valid if For the standard parameter set used in this paper, D~3, m~0:4, t AP~1 , V th~1 , we find the condition 1 ffiffiffiffiffiffiffiffiffiffiffiffi 0:4N a p %1. For instance, for N a~2 00 the relative fluctuations of s 2 are 11%.
Non-dimensional form. It is useful to regard a nondimensional form of Eq. (35). To this end, we measure henceforth voltages in units of V th and time in units of t w , i.e. we introduce the variablesV V~V =V th andt t~t=t w . Furthermore, we will use the adaptation variable A~t w b V th W . Then Eq. (31) and (35) can be rewritten as where the firing threshold isV V th~1 and the adaptation feedback is given by the function The dynamics of the PIF model is completely determined by the four non-dimensional parameterŝ In the last step, Eqs. (38) and (41) were used. The second term of s s 2 can be neglected, because it arises from the small term (rt AP ) 2 in Eq. (41). In fact, rt AP is the ratio between the spike duration and the mean ISI, which can be assumed to be small. Thus, the channel noise intensity is approximatelŷ Again, the adaptation variable A can be written as the sum of a deterministic part a(t) and a stochastic part g a (t), satisfying the equations In accordance with Eq. (37), the mean adaptation is proportional to the firing rate [18], wherer r is the firing rate in units of t {1 w . Mean adaptation approximation. In the case of deterministic adaptation (N a ??), the membrane potential is driven by the currentm m{a and by white noise. If the firing rate of the neuron is large compared to the decay rate of the adaptation current, i.e. ifr r&1, the trace a(t) is effectively smoothed by the linear filter equation Eq. (49). As a consequence, the relative deviations of a(t) from its mean value SaT are small, so that the replacement a(t)?SaT yields a good approximation for calculating the ISI density. In the PIF model, SaT can be selfconsistently determined from Eq. (50), because for the constant input currentm m{SaT the firing rate of the PIF model is simply given byr r~m m{SaT. Solving for the firing rate yieldŝ with l~(1zD) {1 . Thus, the mean adaptation approximation of the PIF model reads Colored noise approximation. In the case of stochastic adaptation (D~0) one can approximate the adaptation variable A~azg a by an effective colored noise. If the adaptation is slow, i.e.r r&1, the spike train looks almost periodic at small time scales with a slowly varying frequency r r(t t)~m m{A (in units of 1=t w ), because of the absence of fast fluctuations. As a result, the instantaneous firing frequency r r(t t) can be well estimated already by averaging over a few spikes. Because the linear filter in Eq. (49) averages a ? (t t) on time scales of the order of 1, the dynamics of A is retained if in Eq. (44) a ? (t) is replaced by its slowly varying average D : r r(t): Multiplication with l~(1zD) {1 yields l _ A A~{AzDlm mz ffiffiffiffiffiffiffiffiffiffiffi ffi 2l 2ŝ s 2 p j a (t t): Finally, we perform the variable transformationĝ g~A{Dlm m, so that we obtain a colored noise driven PIF neuron model l _ g ĝ g g~{ĝ gz ffiffiffiffiffiffiffiffiffiffiffi ffi 2l 2ŝ s 2 p j a (t t): As can be seen from Eq. (53a), the stationary firing rate for stochastic adaptation is again given by Eq. (51).
Nonstationary firing rate. If the stimulus depends on time, i.e. ifm m~m m(t t), the firing rate is a function of time as well. Here, the firing rater r(t t) is understood in the sense of an ensemble average or peri-stimulus histogram. In the PIF model, such a firing rate is given by the ensemble average of the instantaneous rate r r~m m{A, i.e.r r(t t)~m m(t t){SA(t t)T. The ensemble average of A(t t) can be obtained by averaging Eq. (44). Solving for the firing rate yieldsr As an example, let us consider a step of the base current at timê t t Ã such thatm m(t t)~m m 0 zh(t t{t t Ã )m m 1 . In this case, the timedependent firing rate reads This represents exactly the exponential decay of the firing rate in response to a step stimulus, which is typical for spike-frequency adaptation. Eq. (55) has been used in Fig. 1C.
The dynamical response to a step stimulus yields an interesting interpretation of the scaling factor l: from Eq. (55) one finds that which is the ratio of firing rate increases after adaptation to the increase directly after stimulus onset. In Eq. (56),r r 0~lm m 0 , r r Ã~r r(t Ã z) andr r ?~l (m m 0 zm m 1 ) are the firing rates before stimulus onset, directly after the step and att t~?, respectively. Thus, l quantifies the ''percentage adaptation'' or the ''degree of adaptation'' [26].

Statistical measures and known expressions for the PIF model
In the following, T i~ti {t i{1 denotes the interspike intervals (ISIs), i.e. the intervals between adjacent spikes, andT T i~Ti =t w are the corresponding ISIs in units of the adaptation time constant. The statistics of ISIs can be characterized by different measures. A single intervalT T i is distributed according to the probability density P(T T)~Sd(T T{T T i )T, i.e. the normalized ISI histogram. The shape of the ISI density can be characterized using the cumulants. The first cumulant equals the mean ISI and the inverse of the firing rater r, The second cumulant equals the variance and is related to the coefficient of variation (CV), which is a measure of ISI variability: We further consider the third cumulant, which is related to the skewness c s of the ISI distribution, defined as and the fourth cumulant, which determines the kurtosis (or excess) c e of the ISI distribution. It is defined as Roughly speaking, the kurtosis indicates how much of the variability is due to events that strongly deviate from the mean value. For instance, a unimodal ISI density with a heavier tail compared to another ISI density with the same CV, tends to exhibit a larger kurtosis. This is typically accompanied by a more pronounced peak close to the mean value to balance the heavy tail.
In this paper, we want to compare the ISI density with an inverse Gaussian probability density serving as a reference statistic. For an inverse Gaussian ISI distribution (see below) one observes that the skewness is proportional to the CV and the kurtosis scales like the squared CV. This suggests to introduce rescaled versions of the skewness and the kurtosis as follows: ð62Þ By defining the rescaled skewness and kurtosis in this manner, we obtain measures that are identical to unity for an inverse Gaussian ISI density. For larger (smaller) values, the ISI density is respectively more (less) skewed and more (less) peaked compared to an inverse Gaussian density. This procedure is somewhat analogous to the definition of the CV, for which the Poisson process serves as a reference for the ISI variability with C V~1 .
Furthermore, we consider the correlations among ISIs as quantified by the serial correlation coefficient [67]: where due to stationarity the expression does not depend on the integer i, i.e. on the position in the sequence of ISIs. Alternatively, the averages involved in Eq. (63) can also be taken along the sequence (. . . ,T T i{1 ,T T i ,T T iz1 , . . .).
For the adapting PIF model, the two limit cases that are considered in this paper could be reduced to simplified models, for which analytical results are partly known. Firstly, in the case of deterministic adaptation, i.e. N a ??, the ISI density can be approximately described by the first-passage-time density of a biased Brownian motion described by Eq. (52) (mean adaptation approximation). A classical result for this purely white noise driven PIF neuron with a constant drift lm m is that the ISI density is given by the so-called inverse Gaussian [2] This distribution has a mean and a CV Furthermore, by construction we have The mean adaptation approximation would wrongly predict that the ISIs are uncorrelated. The reason is that in the PIF model driven by only white Gaussian noise any memory of the ISI history is erased upon reset. A better account of ISI correlations is given below.
Secondly, the case of stochastic adaptation can be reduced to a PIF neuron driven by a reduced base current lm m and a colored noise with correlation time l and variance lŝ s 2 (Ornstein-Uhlenbeck process, see Eq. (53)). In the case of a weak colored noise, it is useful to define the small parameterŝ For %1 an approximation for the ISI density is given by [48] with c 1 (T T)~T T=lze {T T=l {1, c 2 (T T)~1{e {T T=l and mean ISI ST T i T. ForT T%l, i.e. for ISIs much smaller than the correlation time of the colored noise, the expression for P OU (T T) reduces to [50]: Although this formula captures the ISI density at small ISIs, it is of limited use, because the second and higher ISI moments diverge. Throughout the paper we have therefore used the full expression Eq. (69). The mean ISI and the firing rate do not depend on the noise statistics, in fact they are equal to the white noise case, Eq. (65), as shown below (derivation of the ISI cumulants). The squared CV can be obtained to second order in [48] using the methods presented below: with d~(l 2m m) {1~ST T i T l. Similarly, the rescaled skewness and kurtosis are derived for weak colored noise below. Finally, the serial correlation coefficient can be computed analytically for weak noise [48] (nw0): Serial correlation coefficient of adaptive PIF neuron driven by white noise Here we derive an expression for the serial correlation coefficient of a PIF neuron with deterministic adaptation current and white noise driving. We consider the following subthreshold augmented with the usual reset rule forV V at the spike threshold V V th~1 . The adaptation variable a jumps by an amount D at each spiking event.
Deterministic limit cycle. If we setD D~0 a trajectory in the plane spanned by the variablesV V and a will converge to a stable limit cycle as sketched in Fig. 14. It consists of three segments: (i) the drift ofV V from reset to threshold which lasts for a timeT T Ã and the simultaneous decrease of a from a Ã to a Ã e {T T Ã , (ii) the instantaneous jump of a by an amount D following the threshold crossing, and (iii) the instantaneous reset ofV V . The time of the drift partT T Ã amounts to the whole period of the limit cycle, i.e. the interspike interval of the deterministic system. From this it is clear, that the limit cycle is determined by the two conditions Figure 14. Illustration of the deterministic dynamics of the adaptive PIF neuron. The flow in the phase space spanned by (V,a) possesses a stable limit cycle, which consists of the voltage drift from V reset at a~a Ã to the threshold (A), the following increase of a by an amount D (B) and the subsequent voltage reset (C). The period of the limit cycleT T Ã is solely due to the time of the process A, whereas the processes B and C occur instantaneously. A trajectory starting off the limit cycle at (V reset ,a Ã zda i ) reaches the threshold at the timeT T Ã zdT T i . The deviation from a Ã after one period has an absolute magnitude jda iz1 jvjda i j. doi:10.1371/journal.pcbi.1001026.g014 Solving for the limit cycle parameters a Ã andT T Ã yieldŝ where we introduced the abbreviation Deviations from the limit cycle due to a small perturbation. The serial correlation coefficient for the case of weak noise can be understood by considering small deviations from the limit cycle. Weak noise leads to a distribution of a values upon reset, which is centered about a Ã . Hence, the value of a immediately after reset can be represented as a i~a Ã zda i , where da i is a random number with Sda i T~0. The index i denotes the index in the sequence of interspike intervals, i.e. a i is the initial value of a at the beginning of the ith ISI.
Similarly, noise leads to a distribution of ISIs centered aboutT T Ã and for small noise intensities the mean ISI coincides withT T Ã . Thus, the ith ISI can be represented asT T i~T T Ã zdT T i , where dT T i is a random number with SdT T i T~0. The deviation from the limit cycle after one period, da iz1 , is a function of da i and dT T i . A simple calculation reveals that which after expanding e {dT T i to first order in dT T i and omitting the higher-order term da i dT T i =a Ã reduces to or solving for dT T i : The last equation allows to express the interval correlations SdT T i dT T izn T in terms of the correlations c n :Sda i da izn T: Note, that c n is symmetric, in particular, c {1~c1 .
Correlations of subsequent deviations from the limit cycle. In order to calculate the correlations c n~S da i da izn T the statistics of dT T i for a given initial perturbation da i is needed. To this end, we assume that dT T i can be represented in the following form: where X i is an independent random number with SX i T~0 and SX i X j T~d ij . The functions m(da i ) and s T (da i ) denote for a fixed da i the mean deviation ofT T i fromT T Ã and the standard deviation ofT T i , respectively. The representation Eq. (83) is indeed justified, because in the weak-noise limit the distribution of dT T i for fixed da i is close to a Gaussian distribution. Since any Gaussian random variable is completely determined by its first two moments, it can always be put in the form Eq. (83). Furthermore, for small perturbations da i one can expand the mean deviation up to first order yielding m(da i )&m'(0)da i (the zeroth-order term must vanish, since ST T i T~T T Ã for da i~0 ). Inserting Eq. (83) into Eq.(80) yields a stochastic map for da i : where we defined Multiplying Eq. (84) by da iz1{n and averaging over the full ensemble gives where we have taken into account that for n §1 the perturbation da iz1{n is uncorrelated with X i . Mean ISI for a given perturbation da i . It remains to determine m(da i )~m'(0)da i zO(da 2 i ), i.e. how much the mean ISI is changed by a small perturbation da i . To this end, we integrate Eq. (73) from reset to threshold and obtain Linearizing the exponential e {dT Ti and averaging Eq. (87) over an ensemble with fixed da i (neglecting the higher-order noise term) yieldsV which after subtraction of the unperturbed equation (Eq. (76)) leads to We mention, that this expression is consistent with the result in [68] in the limit of vanishing noise and D%m m. Inserting the linear coefficient m'(0) into Eq. (85) gives therefore an expression for q: Using Eq. (82) and Eq. (86) in the expression for the serial correlation coefficient r n~S dT T i dT T izn T=SdT T 2 i T, n §1, leads to Eq. (9). Our result is the limiting expression for vanishing noise (in fact, Eq. (9) does not depend on the noise intensity at all). Corrections for moderate noise intensity may be feasible with the results in [68].

ISI cumulants of a PIF neuron driven by colored noise
Following [48], we derive here the ISI cumulants of the colorednoise driven PIF model, Eq. (53), in the weak-noise approximation. The ISI cumulants are required to compute the rescaled skewness and kurtosis. The Fokker-Planck equation associated to Eq. (53) for the time-dependent joint probability density p(V V ,ĝ g,t t) reads For the description of the ISI density it is necessary to use the initial condition that corresponds to the distribution ofĝ g upon spiking. The initial condition is for weak noise,~ŝ s 2 (lm m 2 )%1, well approximated by (see [48]) The ISI density is equal to the time-dependent probability flux across the threshold lineV V~1 if threshold crossings of trajectories with negative _ V V V V , i.e. crossing from above the threshold, are prohibited. This is achieved by imposing a reflecting boundary on the half lineV V~1,ĝ gv{lm m. For %1, however, negative _ V V V V are highly unlikely, so that the free process without reflecting boundary generates a flux that is a reasonable approximation of the ISI density.
To carry out a weak-noise expansion of Eq. (91) we change to the variables In these variables the Fokker-Planck equation forp p(x,y,t') reads which arises from the subsequent application of a Laplace and a Fourier transformation to the probability densityp p(x,y,t'). The cumulants k n of the ISI density can be obtained from the characteristic function, see e.g. [69]: k n~( {l) n d n ds n lnG(s) From Eq. (97) and (99) it is clear, that all cumulants can be generated from the function Q(x,v,s). Applying the Laplace and Fourier transformation (as in Eq. (98)) to Eq. (94) we arrive at an equation for Q, which reads with the boundary condition Q(x,v,s)~0: To solve Eq. (100) for weak noise, we expand Q in powers of the small parameter ffiffi p : Inserting into Eq. (100) gives a hierarchy of first-order differential equations for the coefficients Q n : L x Q n zvL v Q n~{ (szv 2 )Q n ziL x L v Q n{1 , n §2: ð105Þ The solutions can be obtained order-by-order using the method of characteristics. Here, we report only the first three coefficients Q n :

Traub-Miles model with adaptation
To test the generality of our findings, we simulated the conductance-based Traub & Miles model modified by B. Ermentrout [51]. It is a single compartment model with an additional M-type current, i.e. a slow voltage-dependent potassium current, inducing spike-frequency adaptation. In order to contrast the effects of deterministic versus stochastic adaptation on the firing statistics of the conductance-based model, we simulated two versions with either additive white Gaussian noise or adaptation channel noise. For the first model with fast fluctuating current noise and deterministic adaptation, the membrane potential V measured in mV is determined by where C m~1 mF cm 2 denotes the membrane capacitance, I is the base current, and D indicates the intensity of Gaussian white noise with correlation function Sj(t)j(t')T~d(t{t'). The deterministic ionic currents are given by the following equations [51]: Sodium current: I Na~ g g Na m 3 h(V {E Na ) g g Na~1 00mS cm 2 , E Na~z 50mV In this model, the adaptation time constant t w (V ) is a voltagedependent function that we reparameterized such that t w roughly corresponds to the time constant governing the exponential buildup of w during periodic firing at 100 Hz in a simulation of equation (123) without the I M current and with D~0.
For the second model with adaptation channel noise, the voltage is described by The currents I Na , I K and I L are the same as in the first model. The M-type adaptation current, however, is modeled as a stochastic current I s M~ g g M W (V {E M ) where W~N op N a indicates the fraction of open channels.
As in the PIF model, we assumed the adaptation channels to be two-state ion channels with the transition rates a w (V )~w ? (V )=t w (V ) and b w (V )~(1{w ? (V ))=t w (V ). The gating of the adaptation channels was simulated using the Gillespie algorithm [70,71]. This algorithm calculates the time interval until the next state transition, determines the reaction type, here channel opening or closing, and updates the number of channels in each possible state accordingly. For a given time step, the number of channels N op in the open state is then used to calculate the fraction of open channels W as well as the stochastic adaptation current I s M . Furthermore, since the transition rates depend on V , we restricted the maximal transition time to 1ms.
In the model with stochastic adaptation current, the maximal channel conductance g g M and the constant base current (I~18mA cm 2 ) were chosen to result in a CV & 0.6 and a firing rate of 100 Hz for a simulation of N a~1 00 ion channels carrying the adaptation current. For the simulation with deterministic adaptation current and additive white noise the base current was adjusted (see corresponding figure captions) to yield the same rate while keeping the conductance g g M the same.