Patterns of interval correlations in neural oscillators with adaptation

Neural firing is often subject to negative feedback by adaptation currents. These currents can induce strong correlations among the time intervals between spikes. Here we study analytically the interval correlations of a broad class of noisy neural oscillators with spike-triggered adaptation of arbitrary strength and time scale. Our weak-noise theory provides a general relation between the correlations and the phase-response curve (PRC) of the oscillator, proves anti-correlations between neighboring intervals for adapting neurons with type I PRC and identifies a single order parameter that determines the qualitative pattern of correlations. Monotonically decaying or oscillating correlation structures can be related to qualitatively different voltage traces after spiking, which can be explained by the phase plane geometry. At high firing rates, the long-term variability of the spike train associated with the cumulative interval correlations becomes small, independent of model details. Our results are verified by comparison with stochastic simulations of the exponential, leaky, and generalized integrate-and-fire models with adaptation.


INTRODUCTION
The nerve cells of the brain are complex physical systems. They generate action potentials (spikes) by a nonlinear, adaptive, and noisy mechanism. In order to understand signal processing in single neurons, it is vital to analyze the sequence of the interspike intervals (ISIs) between adjacent action potentials. There is experimental evidence accumulating that the spiking in many cases is not a renewal process, i.e., a spike train with mutually independent ISIs, but that intervals are typically correlated over a few lags (Lowen and Teich, 1992;Ratnam and Nelson, 2000;Neiman and Russell, 2001;Nawrot et al., 2007;Engel et al., 2008) [further reports are reviewed in (Farkhooi et al., 2009;Avila-Akerberg and Chacron, 2011)]. These correlations are a basic statistics of any spike train with important implications for information transmission and signal detection in neural systems (Ratnam and Nelson, 2000;Chacron et al., 2001Chacron et al., , 2004Avila-Akerberg and Chacron, 2011) and man-made signal detectors (Nikitin et al., 2012). They are often characterized by the serial correlation coefficient (SCC) where T i and T i+k are two ISIs lagged by an integer k and · denotes ensemble averaging. ISI correlations can be induced via correlated input to the neural dynamics, e.g. in the form of external colored noise (Middleton et al., 2003;Lindner, 2004), intrinsic noise from ion channels with slow kinetics (Fisch et al., 2012), or stochastic narrow-band input Russell, 2001, 2005;Bauermeister et al., 2013).
Another ubiquitous mechanism for ISI correlations are slow feedback processes mediating spike-frequency adaptation (Chacron et al., 2000;Liu and Wang, 2001;Benda et al., 2005)a phenomenon describing the reduced neuronal response to slowly changing stimuli (Benda and Herz, 2003;Gabbiani and Krapp, 2006). In the stationary state, these adaptation mechanisms are typically associated with short-range correlations with a negative SCC at lag k = 1 and a reduced Fano factor as demonstrated by several numerical (Geisler and Goldberg, 1966;Wang, 1998;Liu and Wang, 2001;Benda et al., 2010) and analytical studies Farkhooi et al., 2011;Urdapilleta, 2011). The correlation structure of adapting neurons can show qualitatively different patterns, ranging from monotonically decaying correlations to damped oscillations when plotted as a function of the lag (Ratnam and Nelson, 2000). Because ISI correlations shape spectral measures (Chacron et al., 2004), they bear implications for neural computation in general. However, a simple theory that predicts and explains possible correlation patterns is still lacking.
In this article, we present a relation between the ISI correlation coefficient ρ k and a basic characteristics of nonlinear neural dynamics, the phase-response curve (PRC). The PRC quantifies the advance (or delay) of the next spike caused by a small depolarizing current applied at a certain time after the last spike (Ermentrout, 1996). For neurons which integrate up their input (integrator neurons), the PRC is positive at all times (type I PRC) whereas neurons, which show subthreshold resonances (resonator neurons), possess a PRC that is partly negative (type II PRC) (Ermentrout, 1996;Izhikevich, 2005;Ermentrout and Terman, 2010). Below we show that resonator neurons possess a richer repertoire of correlation patterns than integrator neurons do.

MODEL
Spike frequency adaptation can be modeled by Hodgkin-Huxley type neurons with a depolarization-activated adaptation current (Wang, 1998;Ermentrout et al., 2001;Benda and Herz, 2003). However, the spiking of such conductance-based models can in many instances be approximated by simpler multi-dimensional integrate-fire (IF) models that are equipped with a spike-triggered adaptation current (Treves, 1993;Izhikevich, 2003;Brette and Gerstner, 2005); adapting IF models perform excellently in predicting spike times of real cells under noisy stimulation (Gerstner and Naud, 2009). Here, we consider a stochastic nonlinear multi-dimensional IF model for the membrane potential v, N auxiliary variables w j (j = 1, . . . , N) and a spike-triggered adaption current a(t): The membrane potential v(t) is subject to weak Gaussian noise ξ(t) with ξ(t)ξ(t ) = 2Dδ(t − t ) and noise intensity D. The dynamics is complemented by a spike-and-reset mechanism: whenever v (T) reaches a threshold υ(t), a spike is registered at time t i = t and v(t) and w(t) = [w 1 (t), . . . , w N (t)] T are reset to v(t + i ) = 0 and w(t + i ) = w r (where t + i denotes the right-sided limit t → t i + 0). At the same time, a(t) suffers a jump by ≥ 0 as seen from Equation (2c), which resembles high-threshold adaptation currents (Wang, 1998;Liu and Wang, 2001). The constant input current μ is assumed to be sufficiently large to ensure ongoing spiking even in the absence of noise. Note that the model is non-dimensionalized by measuring time in units of the membrane time constant τ m ∼ 10 ms and voltage in units of the distance between reset and spike-initiating potential (a typical value is 15 mV). In particular, the adaptation time constant τ a is measured relative to τ m and the unit of the firing rate is τ −1 m ∼ 100 Hz.
An important special case, the adaptive exponential integrateand-fire model (Brette and Gerstner, 2005) with purely spiketriggered adaptation and a white noise current with constant mean is illustrated in Figure 1. It assumes an exponential nonlin- -Trocmé et al., 2003;Badel et al., 2008) and corresponds to N = 0. Time courses of v(t) and a(t) are shown in Figures 1A1,B1 for two distinct correlation patterns possible in this model. The ISIs T i = t i − t i−1 are obtained as differences between subsequent spiking times t i . The sequence T i , T i+1 , T i+2 displays patterns of shortlong-long ( Figure 1A1) and short-long-short ( Figure 1B1), corresponding to a negative SCC, which decays monotonically with the lag k ( Figure 1A3) or to an SCC oscillating with k ( Figure 1B3). In the following, we develop a theory to analyze these and other values v = 0, a = a * , is indicated by a thick black line. For weak adaptation (A2) a short ISI T i causes positive deviations δa i = a i − a * and δa i+1 = a i+1 − a * of peak values leading to long ISIs T i+1 and T i+2 and, hence, to a negative ISI correlation at all lags (A3). Because of the qualitatively different limit cycle for strong adaptation (B2), deviations δa i and δa i+1 differ in sign, yielding an oscillatory correlation pattern (B3). correlation patterns possible in multi-dimensional adapting IF models.

GENERAL THEORY
In our model Equation (2), a(t) is the only variable that keeps a memory of the previous spike times thereby inducing correlations between ISIs. Over one ISI the time course of adaptation is an exponential decay, relating two adjacent peak values a i = a(t + i ) and a i+1 = a(t + i+1 ) by ( Figures A1,B1). We assume that in the deterministic case (D = 0) our model has a finite period T * (i.e., the model operates in the tonically firing regime) and, hence, for D = 0 the map (3) has a stable fixed point The asymptotic deterministic dynamics can be interpreted as a limit-cycle like motion in the phase space from the reset point to the threshold and back by the instantaneous reset [cf. Figures 1A2,B2]. Weak noise will cause small deviations in the period δT i = T i − T * ≈ T i − T i that are mutually correlated with coefficient The peak adaptation values, however, also fluctuate, δa i = a i − a * , and both deviations are related by linearizing Equation (3): A second relation between the small deviations can be gained by considering how a small perturbation in the voltage dynamics affects the length of the period. This effect is captured by the infinitesimal phase response curve (PRC), Z(t), t ∈ (0, T * ) (Izhikevich, 2005;Ermentrout and Terman, 2010) (see Section 4 for the precise definition). During the interval T i + 1 , the voltage dynamics in Equation (2a) can be written asv Compared to the deterministic limit cycle, the dynamics is perturbed by the weak noise and the small deviation in the adaptation δa i e −(t − t i )/τ a yielding in linear response Combining Equations (5), (6) we obtain the stochastic map where Note that local stability of the fixed point a * requires that |αϑ| < 1. The covariance c k = δa i δa i + k of the auto-regressive process Equation (7) can be calculated by elementary means and using Equation (5) we obtain for k ≥ 1: In order to compute α and ϑ via Equation (8), we have to calculate T * and Z(t) (a * then follows from Equation (4)), which can be done for simple systems analytically. Our main result, Equations (8), (9), allows to draw a number of general conclusions. It shows that the SCC is always a geometric sequence with respect to the lag k that can generate qualitatively different correlation patterns depending on the value of ϑ and thus on PRC and adaptation current. Because |αϑ| < 1 and 0 < α < 1, the prefactor A in Equation (9) is always positive. Consequently, ρ 1 is negative for ϑ < 1 and positive for ϑ > 1. Looking at Equation (8), we find that a positive PRC inevitably yields ϑ < 1. This implies that adapting neurons with type I PRC possess negative correlations between adjacent ISIs. Intuitively, a short ISI causes in the following on average a higher inhibitory adaptation during the subsequent ISI. Such an inhibitory current always enlarges the ISI in type I neurons-hence, a short ISI is followed by a long ISI.
The sign of higher lags is determined by the base of the power: for ϑ > 0 correlations decay monotonically, whereas for ϑ < 0 the SCC oscillates. Two special cases are ϑ = 0 with a negative correlation at lag 1 and vanishing correlations at all higher lags and ϑ = 1 where all correlations vanish. Overall, we find five basic patterns corresponding to the cases −α −1 < ϑ < 0, ϑ = 0, 0 < ϑ < 1, ϑ = 1 and 1 < ϑ < α −1 . These basic patterns cover all interval correlations discussed in previous theoretical studies Urdapilleta, 2011). Our geometric formula generalizes the theory for the perfect IF model with adaptation  to more realistic, nonlinear multi-dimensional IF models with adaptation.
The cumulative effect of the correlations can be described by the sum over all ρ k , which determines the long-time limit of the Fano factor and the low-frequency limit of the spike train power spectrum (for a definition of these quantities, see Section 4.2). Evaluating the geometric series yields This shows that adaptation in neurons with type I resetting (ϑ < 1) leads to a negative summed correlation and hence a reduced long-term variability. Furthermore, at high firing rates achieved by a strong input current μ, the sum in Equation (10) can be approximated by In particular, for strong adaptation ( τ a v T ) the sum is only slightly larger than −1/2. Note that by virtue of the fundamental relation lim t→∞ F(t) = C 2 V 1 + 2 ∞ k = 1 ρ k (Cox and Lewis, 1966) (see Section 4.2), the smallest possible value for the sum is −1/2 in order to ensure the non-negativity of the Fano factor F(t). At this minimal value the long-term variability as expressed by the Fano factor vanishes even for a non-vanishing ISI variability as quantified by the coefficient of variation C V . The latter quantity can also be estimated using the weak-noise theory: From Equation (7) one can calculate the variance of a i and using Equation (5) an approximation for C 2 V ≈ δT 2 i /T * 2 can be obtained as follows:

ONE-DIMENSIONAL IF MODELS WITH ADAPTATION
In the simplest case ( where v 0 (t) is the limit cycle solution and Z(T * ) = [f (v T ) + μ − a * + ] −1 is the inverse of the velocityv 0 (T * ) at the threshold, which is always positive. Thus, the PRC is positive for all t ∈ (0, T * ), i.e., one-dimensional IF models show type I behavior. From our general considerations, this implies a negative SCC at lag k = 1. The sign of the correlations at higher lags can be inferred from the sign of ϑ, for which one can show (Section 4) that Because Z(0) > 0, the sign of ϑ is determined by the sign of f (0) + μ − a * . For weak adaptation such that a * < f (0) + μ (achieved by a sufficiently small value of or τ a , Figure 1A), we will have ϑ > 0 and a negative correlation at all lags ( Figure 1A3).
In this case, a short ISI occurring by fluctuation will cause a positive deviation δa i (Figure 1A2, green arrow). Geometrically, it is plausible that such a positive deviation causes a likewise positive deviation δa i+1 in the subsequent cycle ( Figure 1A2, red arrow). Because a positive deviation is associated with a long ISI, the initial short ISI is on average followed by longer ISIs.
In marked contrast, for strong adaptation such that a * > f (0) + μ (achieved by a sufficiently large value of or τ a ), ϑ becomes negative and hence the SCC's sign alternates with the lag. This alternation of the sign can be understood by means of the phase plane. Let us again consider a positive deviation δa i due to a short preceding ISI ( Figure 1B2, green arrow). Becausev 0 (0) = f (0) + μ − a * < 0, the neuron is reset above the v-nullcline and hence hyperpolarizes at the beginning of the interval, i.e., the trajectory makes a detour into the region of negative voltage (corresponding to a "broad reset" in Naud et al. (2008)). A positive deviation δa i leads to a larger detour (green trajectory) causing a sign inversion and hence a negative deviation δa i+1 (Figure 1B2, red arrow). Because a positive (negative) deviation corresponds on average to a long (short) ISI, the alternation of δa i also entails an alternation of the ISI correlations. Thus, the distinction between monotonic and alternating patterns relates to a qualitative distinction of the voltage trace after resetting [cf. "sharp" vs. "broad" resets in Naud et al. (2008)].
As demonstrated in Figures 1A3,B3, our theory works well for the adapting exponential integrate-and-fire model. We next demonstrate the validity of our approach over a broad range of firing rates (Figure 2) for another important 1D model, the adapting leaky integrate-and-fire model (Treves, 1993) for which f (v) = −γv and (here T * has still to be determined from a transcendental equation). Changing the firing rate by varying the input current μ, we find a good agreement for the first two correlation coefficients and the sum of all ρ k ; the approximation of the CV shows deviations from simulation results when the input current μ becomes small (approaching the fluctuation-driven regime). In accordance with previous findings (Wang, 1998;Liu and Wang, 2001;Benda et al., 2010;Nesse et al., 2010;Urdapilleta, 2011), the first correlation coefficient ρ 1 displays a minimum corresponding to strong anti-correlations between adjacent intervals. The correlations at lag 2 can be positive for a finite range of firing rates if the adaptation strength is sufficiently large (Figure 2B), whereas for moderate adaptation we find a negative ρ 2 at all firing rates ( Figure 2A). In both cases, however, the sum of SCCs approaches a value close to −1/2 for high firing rates as predicted by Equation (11) (Figure 2, bottom). This is strikingly similar to experimental data from weakly electric fish, in which some electro-receptors display a monotonically decaying SCC and some show an oscillatory SCC (Ratnam and Nelson, 2000) but all cells exhibit a sum close to −1/2 ( Ratnam and Goense, 2004). Finally, Figure 2 reveals a local maximum of the CV for some suprathreshold current μ-an effect that has been described by Nesse et al. (2008).

GENERALIZED INTEGRATE-AND-FIRE MODEL WITH ADAPTATION
Different correlation patterns become possible if we consider a type II PRC, which is by definition partly negative and can lead to a negative value of the integral in Equation (8), and hence to ϑ ≥ 1. This corresponds to a non-negative SCC at lag 1, which is infeasible in the one-dimensional case. To test the prediction ρ 1 ≥ 0, we study the generalized integrate-and-fire (GIF) model  with spike-triggered adaptation. This model is defined by f 0 (v, w) = −γv − βw and f 1 (v, w) = (v − w)/τ w . Using the method described in Section 4, the PRC is obtained as where ν = γ + 1/τ w , = β+γ τ w − ν 2 4 and w 0 (t) is one component of the deterministic limit-cycle solution [v 0 (t), w 0 (t), a 0 (t)] that we calculated numerically.
In Figure 3B we demonstrate that all possible correlation patterns can be realized in the GIF model and that the predicted SCCs agree quantitatively well in theory and model simulations (for comparison, see the SCC for the LIF in Figure 3A). To each distinct pattern belongs a range of ϑ (Figure 3, left), determined by the area under the weighted PRCZ(t) = a * τ a e − t τa Z(t). The Figures 3A,B) illustrates, why an adapting GIF neuron can show vanishing (Figure 3Biv) or even purely positive ISI correlations (Figure 3Bv). In case of type II resetting, inhibitory input can shorten the ISI because of the negative part in the PRC; here inhibition acts like an excitatory input. Consequently, a short ISI will induce a stronger inhibition (adaptation) that now causes a likewise short interval and results thus in a positive correlation between adjacent ISIs. Also, the shortening effect of the adaption current in the early negative phase of the PRC can be exactly balanced by the delaying effect of the late positive phase of the PRC (pseudo-renewal case, in which the area underZ is zero).

DISCUSSION
We have found a general relation between two experimentally accessible characteristics: the serial interval correlations and the phase response curve of a noisy neuron with spike-triggered adaptation. The theory predicts distinct correlation patterns like short-range negative and oscillatory correlations that have been observed in experiments (Ratnam and Nelson, 2000;Nawrot et al., 2007) and in simulation studies of adapting neurons (Chacron et al., 2000;Liu and Wang, 2001). Beyond negative and oscillatory correlations, we have found, however, that resonator neurons with spike-frequency adaptation can exhibit purely positive ISI correlations or a pseudo-renewal process with uncorrelated intervals. Adaptation currents that are commonly associated with negative ISI correlations (Wang, 1998;Chacron et al., 2001;Liu and Wang, 2001;Chacron et al., 2003;Benda et al., 2010;Nesse et al., 2010) can thus induce a rich repertoire of correlation patterns. Despite the multitude of patterns, there is a universal limit for the cumulative correlations at high firing rates [cf. Equation (11)], which shows that the longterm variability of the spike train is in this limit always reduced in agreement with experimental studies (Ratnam and Goense, 2004).
Our analytical results apply to arbitrary adaptation strength and time scale but require that (1) the noise is weak and white, (2) the deterministic dynamics shows periodic firing with equal ISIs (i.e., a limit-cycle exists) and (3) the adaptation current is purely spike-triggered with (4) a single exponential decay time.
Regarding the weak-noise assumption, we found from numerical simulations quantitative agreement with our theory for values of the coefficient of variation (CV) up to 0.4, which is, for instance, typical for neurons in the sensory periphery (Ratnam and Nelson, 2000;Neiman and Russell, 2004;Vogel et al., 2005). This holds even in the subthreshold regime at low CVs, where the deterministic system does not follow a limit cycle. In this case, T * has to be replaced by the mean ISI. Moreover, we found qualitative agreement even for moderately strong noise with values of the CV up to 0.8, which is typical for cortical non-bursting neurons in vivo (e.g. Figure 3 in Softky and Koch (1993)).
In the absence of a deterministic limit-cycle, i.e., in the fluctuation-driven regime at high CVs, different mathematical approaches have to be employed, such as those based on a hazard-function formalism (Muller et al., 2007;Nesse et al., 2010;Farkhooi et al., 2011). Furthermore, for some parameter sets, we also observed repeat periods of the deterministic system that involved multiple ISIs corresponding to a periodic ISI sequence with T i = T i + n , where the smallest period is n ≥ 2. Such cases can realize bursting (Naud et al., 2008), which we did not consider in the present study. However, we expect that these parameter regimes yield interesting correlation patterns because already in the noiseless case a periodic ISI sequence exhibits correlations between ISIs.
Regarding the last two assumptions, it seems that the analytical derivation cannot be easily extended to the cases of adaptation currents activated by the subthreshold membrane potential ("subthreshold adaptation" Ermentrout et al., 2001;Brette and Gerstner, 2005;Prescott and Sejnowski, 2008;Deemyad et al., 2012) and multiple-time-scale adaptation (Pozzorini et al., 2013). Ermentrout et al. (2001) have shown that the inclusion of subthreshold adaptation can lead to type II PRCs, which according to our theory could qualitatively change the correlation patterns. An adaptation dynamics depending on the subthreshold membrane potential also involves a fluctuating component because v is noisy. According to , this stochasticity could contribute positive correlations. The combined effect of spike-triggered, subthreshold and stochastic adaptation currents on the sign of the SCC is not clear.
The important cases of the fluctuation-driven regime and multiple-time-scale adaptation have been recently analyzed with respect to the first-order spiking statistics including the stationary firing rate as well as the mean response to time-dependent stimuli (Richardson, 2009;Naud and Gerstner, 2012). The secondorder statistics, which describes the fluctuations of the spike train ("neural variability," cf. Section 4.2) and which limits the information transmission capabilities of neurons, is however still poorly understood theoretically in these cases. How adaptation shapes second-order statistics in the cases of multiple adaptation time scales, fluctuation-driven spiking and sub-threshold adaptation is an interesting topic for future investigations.
As an outlook we sketch, how our theory could be used to constrain unknown physiological parameters by measured SCCs and PRCs. For instance, from the mean ISI we can estimate T * = T . Furthermore, knowing ρ 1 = −A(α, ϑ)(1 − ϑ) as well as the ratio ρ 2 /ρ 1 = αϑ one can eliminate ϑ and solve for α. This allows to estimate the unknown adaptation time constant τ a = −T * / ln α and the amplitude of the adaptation current Although experimental PRCs are notoriously noisy (Izhikevich, 2005), the integral over Z(t) determining our estimate of a * is less error-prone. Combining our approach with advanced estimation methods for the PRC (Galán et al., 2005), may thus provide an alternative access to hidden physiological parameters using only spike time statistics.

PHASE-RESPONSE CURVES OF ADAPTING IF MODELS
We use the phase-response curve Z(t ) to characterize the shift of the next spike following a small current pulse applied at a given "phase" t ∈ [0, T * ] of an ISI. More precisely, let us assume that the last spike occurred at time t 0 = 0. Then, the next spike time t 1 of the perturbed limit cycle dynamicsv = f 0 (v, w) + μ − a + δ(t − t ), v(0) = 0, w(0) = w r , a(0) = a * , 0 < t ≤ T * will be shifted by some amount δT(t , ) = t 1 − T * . The infinitesimal PRC can be defined as the limit where the sign has been chosen such that a spike advance (δT < 0) due to a positive stimulation ( > 0) leads to a positive PRC. The definition of Z(t) by the shift of the next spike differs from the PRC that describes the asymptotic spike shift but is equivalent to the so-called "first-order PRC," which is often measured in experiments (Netoff et al., 2012).

Adjoint equation and boundary conditions
The PRC can be computed using the adjoint method (see e.g. Ermentrout and Terman (2010)). To this end, the dynamics is linearized about the T * -periodic limit cycle solution y 0 (t) = [v 0 (t), w 0 (t), a 0 (t)]. The linearized limit-cycle dynamics y(t) = y 0 (t) + δy(t) corresponding to Equation (2) is given bẏ with the Jacobian matrix The linear response of the ISI to perturbations of the limit-cycle dynamics in an arbitrary direction is given by the vector The remaining N + 1 boundary conditions are obtained by the following consideration: On the limit cycle , a phase φ : → [0, T * ] can be introduced in the usual way by inverting the map t → y 0 (t) and setting φ = t. Because we are interested in the shift of the next spike, it is useful to define the isochrons (sets of equal phase) as the sets of all points in phase space that will lead to the same first spike time. Put differently, phase points belonging to the same isochron will have their first threshold crossing in synchrony. As a consequence, the threshold hyperplane defined by the condition v = v T is a special isochron corresponding to the phase φ = T * . Note that this definition of the phase implies that the reset line defined by the condition v = 0, w = w r does generally not correspond to φ = 0 but to positive phases if a < a * and negative phases if a > a * . Thus, off-limit-cycle trajectories suffer a phase jump upon reset. Close to the threshold, the isochrons are parallel to the threshold, and thus, a perturbation perpendicular to the vdirection does not change the phase. This insensitivity implies the boundary conditions Z w 1 (T * ) = . . . = Z w N (T * ) = Z a (T * ) = 0. Note that a definition of the PRC based on the asymptotic spike shift would require periodic boundary conditions (Ladenbauer et al., 2012 subject to the boundary conditions Z w k (T * ) = 0, k = 1, . . . , N.
The PRC with respect to a is determined bẏ The matrix in Equation (21) is again evaluated on the limit cycle at v = v 0 (t), w = w 0 (t) and is therefore time-dependent. An analytic solution of Equation (21) is possible for one-dimensional models with adaptation (N = 0) or general linear IF models although in most cases the deterministic period T * still has to be computed numerically.

RELATION BETWEEN SECOND-ORDER STATISTICS OF SPIKE COUNT, SPIKE TRAIN AND INTERSPIKE INTERVALS
A stationary sequence of spike times {. . . , t i−1 , t i , t i+1 , . . . } is often characterized by the statistics of the spike train x(t) = i δ(t − t i ), the spike count N(t) = t 0 dt x(t ) or the sequence of ISIs {T i = t i − t i−1 }. In particular, neural variability can be quantified by the second-order statistics of these different descriptions as, for instance, the spike train power spectrum the Fano factor and the coefficient of variation C V = (T i − T i ) 2 / T i and SCC ρ k as defined in Equation (1). These statistics are connected by the fundamental relationship (Cox and Lewis, 1966) (see also (van Vreeswijk, 2010)) It shows that the summed SCC has a strong impact on the long-term variability of the spike train. In particular, a negative sum yields a more regular spike train on long time scales than a renewal process with the same C V .