Statistical structure of neural spiking under non-Poissonian or other non-white stimulation

Nerve cells in the brain generate sequences of action potentials with a complex statistics. Theoretical attempts to understand this statistics were largely limited to the case of a temporally uncorrelated input (Poissonian shot noise) from the neurons in the surrounding network. However, the stimulation from thousands of other neurons has various sorts of temporal structure. Firstly, input spike trains are temporally correlated because their firing rates can carry complex signals and because of cell-intrinsic properties like neural refractoriness, bursting, or adaptation. Secondly, at the connections between neurons, the synapses, usage-dependent changes in the synaptic weight (short-term plasticity) further shape the correlation structure of the effective input to the cell. From the theoretical side, it is poorly understood how these correlated stimuli, so-called colored noise, affect the spike train statistics. In particular, no standard method exists to solve the associated first-passage-time problem for the interspike-interval statistics with an arbitrarily colored noise. Assuming that input fluctuations are weaker than the mean neuronal drive, we derive simple formulas for the essential interspike-interval statistics for a canonical model of a tonically firing neuron subjected to arbitrarily correlated input from the network. We verify our theory by numerical simulations for three paradigmatic situations that lead to input correlations: (i) rate-coded naturalistic stimuli in presynaptic spike trains; (ii) presynaptic refractoriness or bursting; (iii) synaptic short-term plasticity. In all cases, we find severe effects on interval statistics. Our results provide a framework for the interpretation of firing statistics measured in vivo in the brain.


Introduction
The biophysics of action potential generation in a single nerve cell is well understood in the framework of the celebrated equivalent circuit model by Hodgkin and Huxley (Koch 1999). However, there is a huge difference between the properties of neurons in isolation as described in the Hodgkin-Huxley model and in vivo, i.e. embedded in a network of other neurons. In particular, spontaneous neural activity is largely caused by the massive quasi-random input from surrounding cells (Destexhe et al 2003), which also strongly affects a cell's computational properties (Brunel et al 2001;London et al 2010). Hence, when modeling single neuron activity, it is vital to properly account for the characteristics of the spike trains that constitute this input.
Theoretical studies often assume that input spike trains have no temporal structure but are completely random in time with a spiking statistics given by a Poisson process. This assumption has been successfully used to explain and self-consistently determine the irregular, Poisson-like firing patterns of cortical neurons and the emergence of network oscillations (see, e.g. (Brunel 2000)). However, in biological neural networks the Poisson assumption is rarely strictly fulfilled (Baddeley et al 1997). Sources that induce temporal correlations are signal-related processes (Baddeley et al 1997), neuronal refractoriness (Câteau and Reyes 2006) and bursting (Bair et al 1994), adaptation (Wang 1998) and network-generated oscillations (Buzsáki and Draguhn 2004). Also the synapse with its conductance dynamics (Koch 1999) as well as short-term synaptic plasticity (Fortune and Rose 2001) contribute to correlate the effective input seen by the postsynaptic cell. As a consequence of both the nonlinear neural dynamics and the temporally structured driving, the interspike-interval (ISI) statistics of the driven cell is complex (Bair et al 1994;Baddeley et al 1997;Compte et al 2003;Nawrot et al 2007).
Mathematically, it is convenient to neglect any correlations in the input because it implies that the postsynaptic ISI can be modeled as the first-passage time of a Markov process. For the important class of integrateand-fire models with white noise (Burkitt 2006) efficient numerical schemes (Richardson 2008) and in simple cases even explicit analytical expressions exist (Gerstein and Mandelbrot 1964) that fit the ISI histograms of some neurons surprisingly well (Gerstein and Mandelbrot 1964;Fisch et al 2012). The challenge of analyzing neural activity driven by temporally correlated fluctuations, i.e. by a colored noise, is that we need to consider a non-Markovian process, for which no standard techniques exist to compute the first-passage-time (i.e. ISI) statistics. For the special case of a colored noise that can be represented by one or a few stochastic differential equations, researchers have used a Markovian embedding (Hänggi and Jung 1995), i.e. an extension of the phase space by the degrees of the noise dynamics ending with a first-passage-time problem in a higher-dimensional space. Analytical approaches using this trick have been largely limited to the case of exponentially correlated noise (for an exception, see (Bauermeister et al 2013)), that can be mimicked by an Ornstein-Uhlenbeck process (one additional degree of freedom). Using perturbation techniques (small or large correlation time or small noise intensity) approximations have been worked out for the firing rate (Brunel and Sergi 1998;Moreno-Bote and Parga 2004;Alijani and Richardson 2011), the spike train's auto-correlation (Brenner et al 2002;Moreno-Bote and Parga 2006), the ISI density (Lindner 2004;Schwalger and Schimansky-Geier 2008) and ISI correlations (Lindner 2004), and the dynamical response (Brunel et II  I  II III I I II   II II II I I  I I II I II II I II   I I I I I I I I  I   I I  I I I I I II   spike Fig. 1 Illustration of the colored-noise problem for neurons receiving massive spike train input through N synapses. Presynaptic spike trains have temporal structure, e.g. due to refractoriness, bursting, or rate modulations, which is expressed by non-flat power spectra S kk (f ), k = 1, . . . , N , for each single spike train (red lines). This is contrasted to a Poisson statistics corresponding to flat (constant) Poisson spectra (blue dashed line). Each spike train is filtered by the synaptic conductance dynamics and subsequently summed to produce a total input current that drives the neuron. The input current is approximately Gaussian with a power spectrum (or auto-correlation function) formed by single spike train statistics, cross-correlations and synaptic filter. The goal of this study is to link this input statistics to the spiking statistics of the postsynaptic neuron (output statistics).
Richardson 2011) and applied to experimental data (Fisch et al 2012). However, the diverse temporal correlation patterns of in-vivo input mentioned above cannot be captured by an exponentially correlated noise.
To account for the variety of correlation structures, a theory is needed that characterizes the firing statistics for an arbitrary input correlation function. In this paper, we put forward a set of explicit formulas that relate the correlation structure of the input noise to the ISI statistics of a perfect integrate-and-fire (PIF) neuron, which is a canonical model for the mean-driven firing regime (Sec. 2.1). The analysis is based on the assumption that the correlated noise input can be represented by a multi-dimensional Ornstein-Uhlenbeck process (Markovian embedding) and that the noise is weak compared to the mean driving current. This allows us to formulate the first-passage-time problem in terms of a multi-dimensional Fokker-Planck equation, which is solved for arbitrary dimensions of the Ornstein-Uhlenbeck process using weak noise approximation techniques. The details of the rather lengthy calculations are presented in the Appendix. In Sec. 2.2, we apply our general results to generic situations of correlated input and compare the ISI statistics to the commonly studied case of uncorrelated Poisson inputs. In Sec. 2.3, we investigate whether our findings for the PIF model also hold for more realistic nonlinear integrate-and-fire models. For the suprathreshold firing regime and weak noise, we use the phase response curve to derive expressions that generalize some of the results for the PIF model. For the exponential integrate-and-fire model, we use these analytical expressions as well as simulation results for the subthreshold regime to show that similar effects of non-Poissonian inputs also occur under more general conditions. We finally discuss further applications of the theory (Sec. 3).

Theoretical framework and analytical results
We consider a perfect integrate-and-fire (PIF) neuron, which represents a reasonable approximation of neural spike generation in the mean-driven firing regime. The membrane potential V (t) obeys where C m is the membrane capacity and I 0 is a constant membrane current. Instances at which V (t) reaches the threshold V T = 15 mV and is reset to V = 0 mV define the spike times {t i } and thus a sequence of ISIs . . , N , is the input spike train of the th presynaptic neuron with spike times t ( ) j . The convolution (j * x in, )(t) = dt j (t )x in, (t − t ) of the spike train with the post-synaptic current j (t) approximates the synaptic filtering by the conductance dynamics of synapse . The statistics of stationary input spike trains can be characterized by the input firing rates ν = x in, (t) and the input covariance matrix C k (τ ) = x in,k (t)x in, (t+τ ) −ν k ν , or equivalently, by the spectral statistics S k (f ) = dτ e 2πif τ C k (τ ). For the case of the commonly used Poisson statistics, the spike train power spectrum is simply a constant equal to S (f ) = ν . In general, however, a spike train has a frequency-dependent power spectrum reflecting the presence of temporal correlations (Fig.1, left).
If the number of inputs is large, I syn (t) is approximately Gaussian with an auto-correlation function that is shaped by both the spatio-temporal correlations C k (τ ) of presynaptic spike trains (Câteau and Reyes 2006;Lindner 2006) and the synaptic filter ( Fig. 1, middle). In this Gaussian approximation, the dynamics Eq. (1) can be cast into the formV = µ+ση(t), where µ subsumes the constant parts of the input, σ 2 is the variance of the synaptic drive and η(t) is a zero-mean, unitvariance Gaussian process with correlation function Fig. 2 Schematic of the probability flow in the space (v, y). Initially, the probability density p(v, y, 0) is concentrated at v = 0 and proportional to the distribution upon firing p spike (y). Snapshots of p(v, y, t) are schematically depicted at three different times (0 < t (a) < t (b) < t (c) ). The density moves towards the threshold for t > 0 and is not reset upon crossing the hyperplane v = V T . The p.d.f. of the nth spike time, P n (t), is equal to the total flux d d y J v (nV T , y, t) through the hyperplane v = nV T . The boundary condition that prevents recrossings with negative velocity µ+σb T y < 0 (region below dashed line) can be neglected for σ/µ 1.
C in (τ ) and power spectrum S in (f ) = dt e 2πif t C in (t), cf. Appendix A.1. In the case of a flat power spectrum (white noise), the model generates a renewal spike train with an inverse Gaussian ISI statistics. In contrast, we are concerned here with the case of a correlated (colored) noise η(t) with an arbitrary correlation function and a non-flat power spectrum leading to a nonrenewal spike train, i.e. the ISI sequence displays serial correlations. One can show that regardless of the input's correlation structure, the output firing rate of our simple model is always r = µ/V T (Bauermeister et al 2013). However, the ISI probability density and its higher moments as well as the serial correlations among ISIs are strongly shaped by input correlations.

Markovian embedding
The voltage dynamics V (t) is non-Markovian because η(t) is correlated in time. Apart from a few special cases (e.g. (Salinas and Sejnowski 2002;Lindner 2004;Droste and Lindner 2014)), exact solutions of the associated FPT problem are not known for non-Markovian processes. To obtain approximations of the FPT statistics, we assume here that the noise is weak as expressed by the small parameter = σ/µ 1. In particular, this excludes the singular case of purely white noise, for which the variance diverges. Furthermore, we assume that there exists a (possibly high-dimensional) Markovian embedding (Hänggi and Jung 1995) such that η(t) can be represented as a projection of a multivariate Ornstein-Uhlenbeck process (OUP): Here, b is a constant projection vector, Y is a ddimensional OUP, ξ(t) is a vector of Gaussian white noise processes obeying ξ i (t)ξ j (t ) = δ i,j δ(t − t ), and A and B are constant matrices determining the drift and diffusion strength of the OUP, respectively. To ensure a finite variance of the OUP, we require that the eigenvalues of A have negative real parts. The Markovian embedding is possible if the correlation function of the noise can be approximated by a sum of exponential functions and damped harmonic oscillations with decay rates and oscillation frequencies given by the (complex) eigenvalues of A. In the frequency domain, this amounts to approximating the power spectrum of the noise by a sum of Lorentzian functions. Thus, by allowing arbitrarily large dimensions of the OUP, the class of attainable correlation functions or power spectra is large. In particular, many biologically relevant correlation functions may be approximated by a sufficiently large number of exponentials, possibly with complex-valued rates. In Sec. 2.2, we show that even power spectra with an algebraic decay proportional to 1/f γ are captured by our approach because such noise arises from infinitely many exponentially correlated processes with a continuum of correlation times (see e.g. Sobie et al (2011)).
Explicitely, the correlation function and the power spectrum corresponding to Eq. (2) are given by C in (τ ) = b T Σe |τ |A T b and S in (f ) = −2b T A 2 + (2πf ) 2 I d −1 AΣb, respectively (Risken 1984). Here, e τ A is the matrix exponential function, I d is the identity matrix and Σ is the covariance matrix of Y. Its elements Σ ij = Y i Y j can be computed from Eq. (34) given in the appendix A.2. Fortunately, the inverse problem, i.e. finding matrices A and B that yield a desired correlation function C in (τ ), does not need to be solved. Indeed, the Markovian representation (2) is only an intermediate step that allows us to formulate the evolution of the system by the multi-variable Fokker-Planck equation (FPE) (Risken 1984) Here, p(v, y, t) is the probability density for the joint process [V (t), Y(t)] (Fig. 2) and D = 1 2 B T B is the diffusion matrix. Remarkably, we will see that the final expressions for the ISI statistics can be re-expressed solely in terms of the correlation function C in (τ ). Thus, an explicit construction of the OUP, i.e. of the Markovian embedding Eq. (2), is not necessary for the application of our theory.
In order to calculate the statistics of the stationary ISI sequence {T i }, we shift the time origin to the zeroth spike, t 0 = 0. The statistics of t n = n i=1 T i is given by the n th -order interval density P n (t). Because the right hand side of Eq. (1) does not depend on V , the time of the n-th spike is equivalent to the first-passage time from V = 0 to the n-fold threshold V = nV T (Middleton et al 2003;Lindner 2004). For an ensemble of trajectories with V (0) = 0, the first-passage-time density is equal to the total probability flux at time t across the threshold (Fig. 2) if trajectories that have reached the threshold are removed from the ensemble. Thus, Besides the natural boundary conditions [p(v, y, t) = 0 for v → −∞ and y i → ±∞], the removal of trajectories can be realized by requiring that p(nV T , y, t) = 0 for all y for which µ + σb · y < 0. The latter ensures that there is no re-flux of probability into the domain v < nV T from above threshold (cf. (Brunel and Sergi 1998)). However, for weak noise, σ/µ 1, negative velocitiesV = µ + σb · Y are highly unlikely, so that this boundary condition can be safely ignored. The initial condition corresponds to the state of the ensemble conditioned on a spike at t = 0. Hence, p(v, y, 0) = δ(v)p spike (y), where p spike (y) = (1 + b · y) p s (y) is the stationary density of the OUP upon firing and p s (y) = exp − 1 2 y T Σ −1 y / (2π) d |Σ| is the stationary density of the OUP at an arbitrary time.

ISI density and spike auto-correlations
As shown in Appendix A.2, the n th -order interval density P n (t) can be re-expressed in terms of the characteristic function q( , k, t) = dv e i v d d y e ik·y p(v, y, t) instead of p(v, y, t). This function obeys a tractable first-order differential equation (Fourier-transformed FPE) (see Appendix A.2). In the solution, the correlation function C in (τ ) = b T Σe |τ |A T b can be identified and re-substituted. The result is For n = 1 this formula yields the ISI density P 1 (t). Our result includes the known special cases for exponentially correlated noise (Lindner 2004) and narrow-band noise (Bauermeister et al 2013). It also allows to relate the output correlation function C out (τ ) = S(t)S(t+τ ) −r 2 to the input correlation function C in (τ ) by the formula (Cox and Lewis 1966a) In Fig. 4D and 5C, below, we truncated the sum at n = 4.

Cumulants and correlations of ISIs
Using the moment-generating functionP n (s) = ∞ 0 dt e −st P n (t), the k-th cumulant of t n is given by κ k,n = (−1) k d k lnP n /ds k | s=0 (Risken 1984). The moment generating function can be rewritten asP obeys a transformed FPE, which still contains mixed derivatives. To solve this equation for weak noise, we use the perturbation expansion ϕ = ∞ k=0 ϕ (k) k . This yields a hierarchy of first-order equations for the coefficients ϕ (k) , which can be solved iteratively (see Appendix A.2). Up to a correction of order 6 , the second and third cumulant of the nth-order interval are obtained as κ 2,n ≈κ 2,n =2r −2 2 h n + 2 g 2 n + C in (n/r)h n (6) κ 3,n ≈ 12r −3 4 g n h n , using the shorthand g n = g(n/r) and h n = h(n/r).
To quantify the variability of ISIs, we consider the (squared) coefficient of variation where sinc(x) = sin(πx)/(πx). This shows that the variance of ISIs is determined by the total power of the sincfiltered input noise. Hence the color of the noise (i.e. the form of S in (f )) and the output firing rate r = 1/ T (affecting the width of the filter but also the prefactor) directly influence the CV. Particular effects of input correlations become apparent in higher-order cumulants of the ISI density. Convenient is the comparison to the cumulants of the inverse Gaussian density, which is the exact solution for pure uncorrelated (white noise) input (Gerstein and Mandelbrot 1964). Here we consider a rescaled version of the skewness involving the third cumulant: α s = T κ 3,1 /(3κ 2 2,1 ). This quantity is equal to unity for any inverse Gaussian ISI density and is smaller or larger unity for a less or more skewed distribution, respectively . For weak noise, we find from Eqs. (6) and (7) That is, the skewness is the ratio of the areas under the functions 2S in (f )sinc(2f /r) and S in (f )sinc 2 (f /r). For a constant power spectrum, these areas are equal yielding α s = 1 as expected for a white-noise-driven neuron. A less or more skewed ISI density can be obtained by weighting the minima or maxima of sinc(2f /r) by the input spectrum, respectively. Correlations between ISIs can be quantified by the serial correlation coefficient (SCC where n is the lag between ISIs. The SCC is given by the second-order cumulants (Lindner 2004), yielding ρ n = κ 2,n+1 − 2κ 2,n + κ 2,n−1 2κ 2,1 The leading-order approximation (second line) shows that the SCC is proportional to the inverse Fourier transform of the filtered spectrum S in (f )sinc 2 (f T ) evaluated at discrete time lags t = n T . This is equivalent to the correlation function of the input noise smoothed by a moving average of window size T . As shown in the Appendix A.2, the leading-order approximations for the CV and the SCC, Eqs. (9) and (13), are consistent with exact analytical results for the variability of the spike count on large time scales, which have been obtained previously (Middleton et al 2003;Sobie et al 2011;Moreno-Bote et al 2014). Indeed, the Fano factor F (t), defined as the variance to mean ratio of the spike count in a time window of length t, can be approximated for general t > 0 as Here, {t} = t − T t/ T denotes the difference between t and the largest multiple of T equal or less than t. The long-time limit of Eq. (14) is exact and yields lim t→∞ F (t) = C 2 V 1 + 2 ∞ n=1 ρ n , which recovers a well-known identity for spike trains (Cox and Lewis 1966b).

Applications
By means of our theory, we now discuss the spiking statistics for three important mechanisms of temporal input correlations. We consider an excitatory and an inhibitory population of N E = 800 and N I = 200 presynaptic neurons, respectively. Each input spike generates an exponential postsynaptic current of the form j k (t) = J k e −t/τ k θ(t), k = E, I with synaptic time constants τ E = 4 ms and τ I = 8 ms (θ(t) denotes the Heaviside step function). In the case of dynamic synapses, the th synapse has its own time-dependent coefficient J k, (t) =J k F (t)D (t) with a constant prefactorJ k and facilitation and depression variables (see Appendix A.6); for static synapses J k, (t) =J k . If not stated otherwise, we use the base current I 0 = 0.02 (here and in the following, we state all current amplitudes and synaptic weights J k in units of C m V T /ms). Note that the synaptic filter alone leads to a colored noise input even if input spike trains are uncorrelated (Poissonian) and short-term plasticity is absent. Most previous studies of colored noise in neurons (Brunel and Sergi 1998;Fourcaud and Brunel 2002;Moreno-Bote and Parga 2006;Lindner 2004;Schwalger and Schimansky-Geier 2008) have investigated exactly this case of Poissonian input spike trains filtered by the conductance dynamics of static synapses. In contrast, we focus in the following on additional sources of input correlations that have a more drastic effect on postsynaptic ISI statistics.

Gaussian stimuli with complex temporal statistics
A number of studies have reported membrane potential fluctuations that indicate synaptic input with a power law spectrum (Pozzorini et al 2013;Destexhe et al 2003). Direct injection of power-law currents have also been used to mimic the statistics of natural stimuli . How does such a powerlaw statistics shape the spiking statistics of a neuron? To answer this question, we consider Poissonian input with a common rate modulation (signal). Specifically, the firing rate of the excitatory input streams is r Here, N is a normalization constant (see Appendix A.4). The power spectrum exhibits a power-law decay with exponent γ in the range between f 1 and f 2 , which implies long-range correlations that decay like t γ−1 . A standard calculation yields the total input spectrum S in (f ) (see Appendix A.4), which is illustrated in Fig. 3A. Although, strictly speaking, such a scale-free spectrum cannot be achieved by a finite-dimensional OUP, the input spectrum can nevertheless be used in our final formulas, provided the variance is sufficiently small. Figs.3B-D show good agreement between theory (red lines) and stochastic simulations (circles). The ISI density is significantly more skewed than an inverse Gaussian density with the same mean and variance (corresponding to white noise input), as observed previously for positive input correlations . The signature of the powerlaw correlations in the input is apparent in the serial correlations of ISIs that decay like n γ−1 . Likewise, the Fano factor increases algebraically like t γ , resembling "fractal" behavior found in auditory neurons (Lowen and Teich 1992; Peterson , the CV decreases, the SCC ρ 1 changes sign and the ISI density changes from less (α s < 1) to more skewed (α s > 1).
et al 2014). This scaling behavior can be understood by an asymptotic analysis of Eqs. (13) and (14) (see Appendix A.4).

Presynaptic refractoriness
Even in the absence of correlated stimuli, neuronal spiking is correlated in time. In contrast to Poisson statistics, neurons cannot exhibit arbitrarily small ISIs due to a refractory period after each spike. This causes correlations between subsequent spikes. A simple description of refractoriness is provided by a renewal point process with reduced probability for short ISIs (Gerstner and Kistler 2002). Such processes display often a reduced variability (C V,in < 1), indicating a spike train that is more regular than a Poisson process (C V,in = 1). In the following, we model presynaptic input spike trains by renewal processes with an inverse Gaussian ISI density (Gerstein and Mandelbrot 1964;Cox and Miller 1965), for which the exact form of the input power spectrum can be inferred (see Appendix A.5). This process can be parametrized by input firing rate and input CV allowing us to continuously vary the intensity and regularity of the input. For regular presynaptic spike trains (C V,in < 1), the spectrum of the total synaptic input current exhibits a dip at low frequencies and a maximum near the firing rate of presynaptic neurons (Fig. 4A) -effects of the refractory period that are seen in a large fraction of cortical cells (Bair et al 1994;Franklin and Bair 1995;Compte et al 2003). In the corresponding input correlation function the refractory period becomes manifest by a negative trough (Fig. 4A). Such input, as predicted by our theory, yields a less skewed ISI density than a Poisson-driven neuron (Fig. 4B). Although the spike correlations in the output does not differ drastically from the case of Poisson driving (Fig. 4D), correlations among intervals display pronounced short-term negative correlations (Fig. 4C) -a feature commonly observed for neurons with adaptation currents driven by fast (uncorrelated) noise (Chacron et al 2000;Liu and Wang 2001;Schwalger et al 2010;. However, here it is entirely due to an ubiquitous property of presynaptic spike trains, namely, a refractory period. Effects of refractoriness on the input side are restricted to small to moderate postsynaptic firing rates. At high output rates the synaptic low-pass filter overshadows refractory effects and induces positive correlations over multiple ISIs. Accordingly, our theory predicts a transition from negative to positive serial correlations and from less to more skewed ISI distributions (Fig. 4E).

Presynaptic burstiness
Another important source of temporal correlations is bursting activity, which is a widespread phenomenon in the nervous system. Bursting leads to spike patterns that are more irregular than a Poisson process because of the presence of short intraburst intervals and long interburst intervals. We mimic burstiness by a superposition of renewal processes with C V,in > 1.
For irregular spike input (here C V,in = 2.5), the input spectrum has increased low-frequency power as compared to the Poisson case (Fig. 5A, red vs. blue line). Such power spectra are observed in cortical cells (Bair et al 1994;Baddeley et al 1997;Compte et al 2003). The burst-like input leads to extended positive input correlations, which, similar to the case of positive power-law correlations, result in a more skewed  ISI distribution and pronounced positive serial correlations (Fig. 5B,D); effects of this are not pronounced at the level of the spike train output correlation function (Fig. 5C). In contrast, ISI correlations in the presence of Poisson input are very small -they are still present due to synaptic filtering.
Our approximations for the moments of the ISI statistics are based on a Gaussian approximation of the input shot-noise and on a perturbation expansion for weak noise ( 1). To test the range of validity of the perturbation theory, we compared the theory for CV and ISI correlations with simulations for different values of the perturbation parameter = σ/µ (Fig. 6). Remarkably, we find excellent agreement even beyond = 1. As expected, the perturbation theory becomes worse for even larger input variability (  2).

Neurons with synaptic short-term plasticity
It is known that input spectra can be changed by synaptic short-term dynamics (Lindner et al 2009;Merkel and Lindner 2010;Rosenbaum et al 2012;Droste et al 2013).
In particular, facilitating (depressing) synapses lead to a low-pass (high-pass) filtering of the input power. Even if we assume completely uncorrelated presynaptic spike trains, we thus expect synaptic dynamics to have similar effects as the presynaptic refractoriness or bursting discussed above. As an example for the difference in postsynaptic spiking statistics, we inspect the interval correlations that are induced by dynamic synapses. To this end, we use in our standard setup synapses endowed with a facilitation (depression) dynamics, using an established model of short-term plasticity (Dittman et al 2000;Merkel and Lindner 2010). Doing so, we find that for Poissonian spike trains impinging on facilitating synapses, the postsynaptic cell shows marked positive ISI correlations (cf. Fig. 7B, red line and circles) due to the low-pass filtering of the input. Conversely, depressing synapses cause pronounced negative correlations (Fig. 7B, green line and squares) because depression reduces input power at low frequencies. Our theory (lines in Fig.7B) is in excellent agreement with the simulations (symbols).

Beyond the perfect integrate-and-fire model
A crucial assumption that permitted a systematic perturbation solution was the voltage-independence of the membrane current I 0 . This is justified if the mean drive is so strong that it dominates any variations of the membrane current due to voltage changes. Here, we briefly discuss the significance of our results if this assumption is relaxed; for a similar analysis for whitenoise-driven IF neurons and neurons endowed with an adaptation current, see . We extend the subthreshold dynamics by a voltage-dependent term: Spikes are registered when V (t) reaches the value V spike ≥ V T , whereupon V is reset to zero. The nonlinear IF model Eq. (16) is a more general description of neural dynamics that can accurately model passive membrane properties (leak) and spike initiation dynamics. In particular, the exponential integrate-and-fire (EIF) model (Fourcaud-Trocmé et al 2003) given by the fits neural data of many neurons remarkably well (Badel et al 2008). The model Eq. (16) can operate in two different dynamical regimes: Firstly, if γ L Ψ (V ) + µ > 0 for all V < V spike , the deterministic dynamics (σ = 0) generates periodic spiking (supra-threshold or meandriven regime). For the EIF model, this condition reads γ L < µ/(V T −∆ T ). Secondly, for γ L > µ/(V T −∆ T ), the noiseless dynamics reaches a steady state below threshold, i.e. the neuron is quiescient unless it is excited by the noise (sub-threshold or fluctuation-driven regime). We will separately consider these two cases.

Mean-driven regime
If γ L < µ/(V T − ∆ T ), the mean ISI is approximately given by the deterministic ISI Furthermore, the dynamics is associated with a phaseresponse curve (PRC) Z(t ) that describes the change of the ISI δT = −αZ(t ) caused by a weak current pulse αδ(t − t ) at a given time t ("phase") in the interval (0, T ) (Ermentrout and Terman 2010). For the nonlinear IF model the PRC reads  Fig. 8 ISI statistics for the exponential integrate-and-fire (EIF) model, Eq. (16), in the mean-driven regime. The neuron is driven by inverse Gaussian renewal processes each with r = 5Hz and C V = 0.5 (cf. Fig. 4). Solid red lines depict the theory, Eq. (19), (20). The black dashed line shows the PIF-theory using an effective base current µ eff = V T / T , where T is given by Eq. (17). A The leak γ L is varied from γ L = 0 (limit of the PIF model) to the maximal value, at which the dynamics becomes subthreshold.
(J E = 0.005, J I = −0.01, µ = 0.1 V T /ms; corresponding to = 0.24) B The output firing rate is varied by the parameter µ (J E = 0.004, J I = −0.008, γ L = 0.01 ms −1 ). Other parameters: (see, e.g. ). Here, V 0 (t) is the deterministic solution of Eq. (16) with V 0 (0) = 0. As shown in Appendix A.7, the CV and the SCC can be related to the power spectrum S in (f ) of a weak colored noise perturbation ση(t) by Here,Z(f ) = 1 T T 0 dt Z(t)e 2πif t . Note that for the PIF model (i.e. γ L = 0), we obtain |Z PIF (f )| 2 = sinc 2 (f /r)/µ 2 . This exactly recovers the leading-order terms of CV and SCC, Eqs. (9) and (13). Figure 8 shows the estimation of the CV and ρ 1 using the PRC Z(t), which agrees well with simulation data. We have also tried an even simpler approximation in this case by using the formulas for an effective PIF model with the same mean ISI given by Eq. (17) (see also ). This approximation (dashed lines in Fig.7) already yields a good agreement. Most importantly, the effect of a presynaptic refractory period on postsynaptic variability and ISI correlations remains qualitatively the same as for the PIF neuron (cf. Fig. 4E and Fig. 8B): with increasing output firing rate, ISI correlations turn from negative to positive and output variability decreases.

Fluctuation-driven regime
To see whether similar effects of non-Poissonian stimuli also hold in the fluctuation-driven (excitable) regime, we performed simulations of the EIF model with sub-threshold currents µ ≤ µ * ≡ γ L (V T − ∆ T ) (Fig. 9). In contrast to the mean-driven regime, where the firing rate is mainly determined by the deterministic ISI, Eq. (17), the firing rate in the fluctuation-driven regime may also depend on the properties of the noise. To study the sensitivity of the spiking statistics to regularity of the input spike trains, we again model input spike trains by inverse Gaussian renewal processes with varying C V,in and fixed firing rates. Figure 9A (left panel) shows that the firing rate in this setup depends only weakly on the regularity of the input, whereas the base current µ has a much stronger effect. The weak increase of the firing rate with C V,in is caused by an increasing low frequency power due to higher input irregularity (cf. Fig. 4A and Fig. 5A). According to our main focus, we further investigated effects on higher-order ISI statistics at different values of µ ≤ µ * . To avoid trivial effects of the firing rate, we chose for each value of µ synaptic weights that yielded comparable firing rates ( Fig. 9A (right panel)). Directly at the bifurcation point µ = µ * from sub-to suprathreshold firing, the ISI densities were unimodal and more (less) skewed than the case of Poisson input for irregular (regular) input (Fig. 9B)similar to the findings for the mean-driven firing regime. Furthermore, the ISI sequence showed significant correlations that were positive over several lags for irregular input (C V,in = 2.5) and was negative at lag one for regular input (C V,in = 0.5). Poisson input did not result in any significant serial correlations. Qualitatively, these results persist in the subthreshold case, µ < µ * .
In addition, we also observed that the ISI density can become bimodal for regular input (Fig. 9C,D). The ISI densities become more skewed with increasing input irregularity, which is confirmed by the corresponding values of α s (Tab. 1). Furthermore, the skewness of the ISI density under Poisson stimulation approaches the value α s = 2/3 for decreasing µ. This corresponds to a gamma distribution, which is expected far below the bifurcation point, µ µ * , where the ISI statistics becomes Poissonian.

Discussion
In this paper, we put forward a theory for the spiking statistics of a mean-driven neuron stimulated by an arbitrary colored Gaussian noise. Within the developed framework, all statistics (ISI density, auto-correlogram, CV, skewness and SCC) can be expressed as nonlinear functions of two simple integrals of the input correlation function or, equivalently, of the input power spectrum. The leading order of the theory permits to understand the interaction between input time scales (e.g. the correlation time of the input stimulus) and the postsynap- Table 1 Values of the modified skewness α s corresponding to the ISI densities shown in Fig. 9B-D. Note that α s increases with input CV. Decreasing µ leads to a transition from an inverse Gaussian distribution (α s = 1) to a gamma distribution (α s = 2/3). Applying our theory, we found that distinct sources of input correlations can lead to similar ISI statistics, all showing pronounced differences to the case of a purely Poissonian drive. For instance, both presynaptic refractoriness and synaptic short-term depression lead to negative interval correlations and less skewed ISI distributions than with Poisson input. Conversely, ratecoded stimuli, presynaptic bursting, and synaptic facilitation cause positive ISI correlations and more skewed ISI densities. Our theory thus shows clearly the limitations for inferring neural dynamics from spike statistics only. Negative correlations, for instance, have often been found for neurons driven by fast fluctuations and adaptation currents, but, as we have shown, can be caused already by a much more basic feature, namely, presynaptic refractoriness. Despite this ambiguity, our formulas may still be applicable to draw conclusions about stimuli and neural parameters (see (Bauermeister  In a recent study (Peterson et al 2014), negative ISI correlations observed in auditory nerve fibres have been attributed to synaptic depression of the ribbon synapse that generates the input to a nerve fibre. In that case, the spike train was modeled as a thinned-out version of stochastic transmitter release events of the synapse. The transmitter release events already carried negative inter-event interval correlations, which were inherited by the thinned-out spike train. Our theory for negative correlations due to synaptic depression provides an al-ternative explanation based on the temporal structure of synaptic input and its effect on a dynamic spike generator.
Relations between input and output spike statistics as developed in this paper may be also instrumental for novel approaches to recurrent neural networks. In the simple situation of a completely homogeneous and sparse network, the statistics of the postsynaptic spike train (forming the input to other cells) has to coincide with the statistics of presynaptic spike trains (Lerchner et al 2006;Dummer et al 2014). Specifically, in our theory one should be able to find a solution for which C in (t) = C out (t). For heterogeneous and non-sparse networks, more self-consistency conditions should be met, not only for the auto-but also for the cross-correlation functions among neurons. In this sense, our paper constitutes only a first step in a theory of neural networks that incorporates the correlation statistics more accurately than the classical approach that uses the diffusion approximation.

A.1 Gaussian approximation of the synaptic input current
We would like to approximate the synaptic input current by a Gaussian process with the same mean and correlation function as I syn (t). To this end, let us rewrite Eq. (21) in terms of the delta-spike trains denotes the i-th spike time at the k-th synapse): In the case of stationary input, mean and correlation function are given by and Here, ν k = x k (t) is the stationary firing rate of the k-th input spike train and C kl (τ ) is the correlation function between spike trains defined by It is convenient to formulate Eq. (24) in the Fourier domain.
Defining the power spectrum as the Fourier transform of the correlation function we find Here and in the following,˜and * denote Fourier transform and complex conjugation, respectively. Using these approximations, we can write the current balance equation (1) aṡ where and ζ(t) being a zero mean Gaussian process with power spectrum If S ζ (f ) is not constant, ζ(t) will be referred to as colored noise (in contrast to white noise, which has a flat power spectrum).

A.2 Derivation of the interspike interval statistics
Our approach assumes that the colored noise can be represented as Specifically, the neuron model can be written aṡ with the additional fire-and-reset rule that V is reset to V = 0 whenever V reaches the threshold V T . A and B are drift and diffusion matrices, respectively. The eigenvalues of A must have negative real part and ξ(t) is a vector of Gaussian white noise processes that obey ξ i (t)ξ j (t ) = δ i,j δ(t − t ). Furthermore, we assume that the elements of b have the same physical dimensions as µ and that Y has non-dimensional units.
The variance of ζ(t) is given by σ 2 = b T Σb, where Σ is the covariance matrix of the OUP with elements The covariance matrix is related to the drift matrix A and the diffusion matrix by the Lyapunov equation (see e.g. (Risken 1984;van Kampen 1992)) This system of linear equations can be used to solve for the elements of the covariance matrix Σ ij .

Non-dimensional model
It is useful to non-dimensionalize the model Eq. (31) by measuring time in units of the mean ISI V T /µ and voltage in units V T . This corresponds to the introduction of dimensionless variableŝ This yields the non-dimensionalized Langevin equatioṅ The non-dimensional firing threshold isV T = 1. To make the small parameter in our perturbation calculation explicit, we have introduced the dimensionless parameter = σ/µ 1. In the following, the non-dimensional model Eq. (37) will be considered. For notational simplicity, we will omit in the following the hats on top of the non-dimensional quantities. The stationary auto-correlation function of the noise η(t) = b · Y(t) is for τ ≥0 given by (see e.g. (Risken 1984)) where e τ A is the matrix exponential defined by its series representation e τ A = ∞ k=0 1 k! (τ A) k . Note that by our rescaling of b, Eq. (36), the variance of η(t) is C in (0) = η 2 = 1. Furthermore, taking the Fourier transform of the correlation function yields the input power spectrum where I d is the d × d identity matrix.

Fokker-Planck equation
The stochastic dynamics (37) with fire-and-reset rule can be reformulated in terms of probability currents and densities. The probability current that describes the drift of V (t) is given by is the probability density of V (t) and Y(t) at time t. The d-dimensional OUP corresponds to d probability currents Here, i = 1, . . . , d and the non-dimensionalized diffusion matrix is defined by D = 1 2B TB . Because trajectories cannot leave the system, the probability is conserved, which can be expressed by the continuity equation The right hand side represents the reinjection and outflux of probability at the reset v = 0 (source) and the threshold v = 1 (sink), respectively, corresponding to the reset of trajectories.
With the definition of the currents, Eq. (40) and (41), the last equation can be written as a partial differential equation for p(v, y, t). Because the source and sink terms are sharply localized, we can omit these terms in the regions v < 0, 0 < v < 1 and v > 0 and we arrive at the Fokker-Planck equation (FPE), Eq.
(3) with non-dimensionalized parameters µ = 1 and σ = . The source and sink terms can be realized by imposing appropriate boundary conditions at v = 0 and v = 1. The sink term at v = 1 ensures that no probability can be found beyond the threshold and thus there is no backflux of probability from above the threshold, i.e. J v (1, y, t) = 0 for 1 + b · y < 0. With Eq. (40) this implies the condition (cf. (Brunel and Sergi 1998)) p(1, y, t) = 0 ∀y : 1 + b · y < 0.
Furthermore, integration of Eq. (42) from v = −ε to v = ε (and letting ε → 0), reveals that J v (v, y, t) suffers a jump at v = 0 by an amount J v (1, y, t). With Eq. (40) this implies the condition where v = 0− and v = 0+ denote the left-and right-sided limits. Finally, probability should not leak across infinite boundaries and the probability density must be normalized: dv d d y p(v, y, t) = 1.
In addition to the boundary conditions, the FPE must be supplied with an initial condition. Because we aim at the statistics of the interspike intervals, it is advantageous to consider an ensemble of neurons that have just fired a spike at time t = 0. This corresponds to preparing the system in a state, in which the membrane potential has been reset to zero, i.e. V (0) = 0, and the OUP is distributed according to the distribution of Y sampled at the moments of firing. If we denote this distribution by p spike (y), the initial condition for the FPE must be of the form p(v, y, 0) = δ(v)p spike (y).

Distribution upon firing
In order to determine p spike (y), we first seek the stationary solution p s (v, y) of the FPE. Setting ∂p/∂t = 0, we obtain the stationary FPE which is still difficult to solve analytically due to the discontinuous boundary condition (43) at v = 1. However, for weak noise, 1, the velocity 1 + b · y, is practically always positive, so that the boundary condition (43) can be safely ignored. A positive velocity also precludes trajectories from entering the region v < 0. As a consequence, the system is restricted to the region 0 ≤ v ≤ 1 and has the periodic boundary condition A solution of Eq. (47) is given by the stationary distribution of the OUP, i.e.
which clearly fulfills the boundary conditions (48). The distribution upon firing p spike (y) is proportional to the stationary probability current J (s) v (1, y) = (1 + b · y) p s (1, y) across the threshold (Lindner 2004). In fact, the probability per unit time that a trajectory crosses the threshold through the surface element dS(y) at the point y on the hyperplane defined by v = 1 is given by J (s) v (1, y)dy 1 · · · dy d . Thus, the probability density upon firing must be proportional to J (s) v (1, y). The factor of proportionality can be found from the normalization condition, which yields p spike (y) = (1 + b · y) p s (1, y) and hence the initial distribution p(v, y, 0) = δ(v) (1 + b · y) p s (y). (50)

Interspike interval density
An interspike interval is equivalent to the first-passage-time (FPT) with respect to the threshold V T = 1, i.e. the time a system that has been prepared according to a spike at t = 0 needs to reach the threshold for the first time. A standard approach to obtain the probability density of the FPT (i.e. the ISI density) is to consider an ensemble in which trajectories are not reset upon threshold crossing but are taken out of the system. As a consequence, the probability to find a trajectory with v < 1, is the probability that a trajectory survives until time t. Here, p(v, y, t) is the time-dependent solution of the FPE, for which the reset condition (44) has been dropped now. The FPT density is given by P (t) = −dG/dt, or Using the continuity equation (42) and the natural boundary conditions (45) yields where is the FPT density with respect to a general threshold V T = v. Equation (53) simply expresses that the FPT is equal to the total probability current through the threshold at time t. The ISI density can be rewritten as whereP (t, ) = dv e i v P(t, v) is the Fourier transform with respect to v. Using Eq. (54), we thus havẽ This function can be further related to the characteristic function where k is a d-dimensional vector with elements k i , i = 1, . . . , d. Indeed, comparing (56) with (57), we find the relatioñ where ∇ k q denotes the gradient of q with respect to k, i.e ∇ k q is a vector with elements ∂ i q ≡ ∂q/∂k i . Thus, if the characteristic function q( , k, t) is known, one can derive the ISI density via (55) and (58). In order to compute q( , k, t), it is not necessary to calculate the probability density p(v, y, t) explicitly from the FPE, which does not permit a closed form solution for the initial conditions (50). However, applying the Fourier transform to the FPE directly yields a first-order differential equation for q( , k, t): with initial condition This equation can be solved by the method of characteristics (see e.g. (Kamke 1965)). The characteristic equations arė The solution of the first equation is which can be solved for the initial vector c: The solution of Eq. (62) is given by where the Lyapunov equation (34) has been used. The prefactor Q 0 (c) can be determined from the initial condition (60) as Q 0 = 1 + i b T Σc. The integration in Eq. (65) has to be performed using the expression (63), followed by a resubstitution of the vector c = c(t, k) according to Eq. (64). As a result we obtain the characteristic function Using Eq. (58) leads tõ and which constitutes the final result for the ISI density. Let us provide two alternative expressions for h(t) and g(t). First, if A is invertible we can write these functions in terms of the drift matrix: Second, using the power spectrum of the noise we find where sinc(x) = sin(πx)/(πx).

Cumulants of the n-th order intervals
Apart from the ISI density the ISI statistics can be characterized by its cumulants (semi-invariants). We will calculate these cumulants also for the nth-order intervals because they will be used to find an expression for the SCC and the calculation does not pose any extra difficulties. In fact, as argued in the main text, the time t n of the n-th threshold crossing, i.e. the nth-order interval, is equal to the first-passage time of the free process with respect to the nth-fold threshold V T = n (given the same noise realization). The probability density of the nth-order interval is thus given by P(t, n) (cf. Eq. (54)). The k-th cumulants of the nth-order intervals is given by (van Kampen 1992) wherē is the Laplace transform of the n-th-order interval density, which plays the role of a moment generating function. Introducing the Laplace-Fourier transform of the probability density the moment generating function can be rewritten as The function ϕ satisfies the transformed FPE The mixed derivative ∂ v ∇ k ϕ makes an exact solution infeasible. However, using the perturbation expansion yields a hierarchy of first-order partial differential equations for the coefficients ϕ (k) (v, k, s). Specifically, we obtain where the the inhomogeneities I m are given by If we insert the ansatz (76) into the formula (74), we obtain a perturbation series for the moment generating function: whereP (0) n (s) = ϕ (0) (n, 0, s) and for m ≥ 1. Equation (77) The integration of the first two equations results in To solve the third equation, we make the ansatz Using this ansatz in Eq. (85) yields whereÎ n (τ, c) = I n exp sτ + 1 2 k T Σk . In terms of the lower-order functions ψ n−1 , the inhomogeneous parts can be written aŝ In order to regard these expressions as functions of τ and c, the vector k has to be understood as k(τ, c) = e −τ A T c. It turns out that in Eq. (88), the prefactor in front of ψ (m) vanishes because in light of Eqs. (34) and (84) Thus, ψ (m) can be easily integrated: For the evaluation of this formula, the derivatives of the lowerorder functions ψ (m−1) (τ, c) with respect to v and k are needed. To this end, it is useful to apply the chain rule: for a given scalar function f (τ, c) we can write and In Eq. (94), we made use of ∇ k τ = 0 and ∇ k c = e τ A T T = For m ≥ 1, the perturbation coefficientP (m) n (s) can be expressed in terms of ψ (m−1) (τ, c) by using Eqs. (82), (93) and (95): (96) Zeroth-order. For m = 0, we obtain from Eq. (93) Equation (87) then yields the leading-order of the momentgenerating function The derivatives of ψ (0) simply read First-order. Because the gradient in Eq. (96) equals zero for m = 1 (Eq. (99)), the first-order contribution toP n (s) vanishes: From Eqs. (90) and (93) we obtain Using Eq. (94), the gradient of ψ (1) reads Furthermore, the derivative with respect to v is and the mixed derivative takes the form Second-order. For m = 2, we find the first non-vanishing correction of the moment-generating function due to noise: Here, C in (t) is the noise correlation function as defined in Eq. (38). Furthermore, we obtain and the gradient The derivative with respect to v is more involved: using Eq. (95) first leads to Using the identity e −τ A T A T = − d dτ e τ A T and integrating by parts simplifies the expression to Taking the derivative of the last equation with respect to k results in where I d denotes the d × d identity matrix.
Third-order. Because the gradient ∇ k ψ (2) vanishes for c = 0, we find Using the derivatives in Eqs. (106), (108) and (109) we obtain the third-order perturbation solution Differentiating this expression with respect to k, setting c = 0, and taking the scalar product with the vector b leads to Here, C in (τ ) denotes the derivative of the correlation function with respect to τ .
Fourth-order. Finally, using Eqs. (96) and (112) yields the second non-vanishing correction: To summarize, the moment-generating function reads up to fourth order in : P n (s) = e −sn 1 + 2 s 2 n 0 dτ 2 This formula satisfies the general propertyP n (0) = 1 of a moment-generation function, which expresses the normalization ∞ 0 dt P n (t) = 1 of the nth-order interval density.
Cumulants. Knowing the moment-generating function, Eq. (114), we can compute the cumulants of the nth-order intervals by means of Eq. (71) up to fourth-order in . As expected, the first cumulant, or mean nth-order interval, is given by which is consistent with the fact that t n = T 1 + · · · + T n = n T i = n (recall that T i = 1 in the non-dimensionalized system).
The second cumulant is equal to the variance ∆t 2 n of the nth-order intervals, where ∆t n = t n − T i . From the formula for the cumulants, Eq. (71), it can be seen that only the coefficients of the quadratic terms s 2 in Eqs. (113) and (114) enter the expression for the second cumulant (P n (s) has no linear terms). Thus, we find κ 2,n ≡ ∆t 2 n = 2 2 n 0 dτ 2 The third cumulant only involves the cubic terms in Eq. (113), hence

Coefficient of variation and serial correlation coefficient
Formulas for a given noise correlation function. The coefficient of variation is given by the ratio of standard deviation and mean of the 1st-order interval (i.e. ISI). If we take into account that we measure time in units of the mean ISI, we find where κ 2,1 is given by Eq. (116) with n = 1. Furthermore, we can compute the SCC by using formula (116) for arbitrary n: If we further restrict ourselves to the leading order, we can derive a much simplified expression for the SCC in terms of the noise correlation function. In fact, writing ∆t 2 n = 2 2 n 0 dτ 2 τ 2 0 dτ 1 C in (τ 1 ) + O( 4 ) and accounting for the symmetry of C in (τ ), we obtain from Eq. (119) In terms of the noise power spectrum, Eq. (70), this expression can be rewritten as Formulas for a given matrix representation of the Ornstein-Uhlenbeck process. To elucidate the general structure of serial correlations, it is instructive to express the SCC in terms of the drift matrix A and covariance matrix Σ of the OUP. In leading order of , the variance of the nthorder intervals reads where we assume that the drift matrix A is invertible. Setting n = 1, we obtain the squared coefficient of variation: Using Eq. (119), the leading order of the SCC can now be written as (123) Apparently, the SCC has the simple exponential structure ρ n = c T 1 e (n−1)A T c 2 with constant vectors c 1 and c 2 . That means that the SCC is a linear combination of the entries in the matrix exponential function e (n−1)A T . In fact, let Q be the constant matrix made up by the (generalized) eigenvectors of A T . Then, A T can be diagonalized (or, or more generally, transformed into its Jordan normal form) through the representation A T = QJQ −1 . Thus, the SCC can be rewritten as where the vectors d 1 = Qc 1 and d 2 = Q −1 c 2 are constant. In the special case, where A T is diagonalizable (i.e.A T possesses d linearly-independent eigenvectors corresponding to the eigenvalues λ i , i = 1, . . . , d), the matrix exponential e (n−1)J is diagonal with the diagonal elements e λ i (n−1) . Thus, the SCC is simply a superposition of the exponentials e λ i (n−1) . Real eigenvalues contribute to exponential decays with "rate" 1/λ i . On the other hand, complex eigenvalues lead to damped oscillations, which, however, might appear rather irregular as e λ i (n−1) is sampled at discrete integer values (Bauermeister et al 2013). In any case, the SCC decays to zero for n → ∞ because, by assumption, the real parts of the eigenvalues of the drift matrix are negative.

Cumulative correlations and Fano factor
The exponential structure of the SCC, Eq. (123), also allows us to calculate the cumulative effect of serial correlation as expressed by the sum of all SCCs. Applying the formula for the geometric series, (125) This sum also determines the long-term variability of the spike train. Indeed, the fundamental relation (Cox and Lewis 1966a) lim t→∞ relates the Fano factor F (t) of the spike count or the power spectrum S out (f ) of the output spike train to the second order ISI statistics as expressed by the CV and the sum of SCCs. The Fano factor is defined as the variance to mean ratio of the spike count N (t) in the time window (0, t), i.e. F (t) = ∆N (t) 2 / N (t) (note that, in this definition, the origin of the time window is arbitrary; it does not have to coincide with a spike time). The power spectrum of the output spike train x out (t) = i δ(t − t i ) is defined as the Fourier transform of the spike train auto-correlation function, i.e. S out (f ) = dτ e 2πif τ x out (t)x out (t + τ ) . Using the expression for the squared CV, Eq. (122), and the summed SCC, Eq. (125), the long-time limit of the Fano factor is In order to obtain an approximation for the Fano factor for arbitrary length of the time window t, we follow (Middleton et al 2003): If t is large and 1, the spike count N (t) can be approximated by the free membrane potential V (t) with V (0) ∈ [0, 1). Indeed, for weak noise we can as-sumeV (t) > 0, so that the n-th spike time t n is equivalent to the crossing of the n-fold threshold, i.e. V (t n ) = n. This implies that the deviation δV (t) = V (t) − N (t) is bounded (in fact, 0 ≤ δV (t) < 1) and hence, N (t) = V (t) + O(1). This relation still holds if we assume the specific initial condition V (0) = 0. The trajectory of V (t) can then be obtained by integrating Eq. (37): Thus, the mean and variance are where h(t) has been defined in Eq. (66). For large times, the Fano factor is hence given by In the exact limit t → ∞ we re-obtain the leading-order term of Eq. (127), which shows the consistency of the two approaches.
For small time windows, the spike train can be regarded as periodic given the weakness of the noise. This argument holds irrespective of the specific correlation function of the noise, which allows us to apply the same approximation as in (Middleton et al 2003): where {t} = t − t is the difference between t and the largest integer t not exceeding t (sawtooth function). Following (Middleton et al 2003), a uniform approximation on the whole range of t can be obtained by adding the small and large time approximations of the Fano factor:

A.3 Model of the synaptic input current
To test our theory (Figs. 3-9), we assume that the pool of presynaptic neurons can be split up into N E excitatory and N I = N − N E inhibitory neurons. Furthermore, the postsynaptic synaptic currents induced by each spike are modeled by the exponential impulse response functions j E (t) = J E e −t/τ E θ(t) and j I (t) = J I e −t/τ I θ(t) with synaptic time constants τ E and τ I (synaptic filtering). This corresponds to the simple first order kinetics where x E/I,k (t) = j δ(t − t (E/I,k) j ) is the spike train arriving at the k-th excitatory/inhibitory synapse. J E and J I denote the synaptic efficacies. The mean synaptic current is thus I syn (t) = N E J E τ E ν E + N I J I τ I ν I , where we assumed that excitatory and inhibitory neurons fire with rate ν E and ν I , respectively.
In frequency domain, the synaptic filter becomes for i = E, I. In the special case, where the single presynaptic spike trains are independent and have power spectra S E (f ) and S I (f ), respectively, we find the power spectrum of the synaptic current from Eq. (27):

A.4 Signal encoded in a common rate modulation of Poissonian inputs
To model the encoding of a continuous signal s(t), we assume that s(t) represents the rate modulation of an inhomogeneous Poisson process. Specifically, we assume that the k-th presynaptic spike train is a Poisson process with rate r k (t) ≡ x k (t) = ν k (1 + ε k s(t)).
Here, s(t) is common to all inputs and is modeled as a meanzero Gaussian process with power spectrum S s (f ). For the correlation matrix, Eq. (25), we find S kl (f ) = ν k δ kl + ε k ε l ν k ν l S s (f ).
The power spectrum S I (f ) of the total input current is then given by Eq. (27), and the normalized spectrum is S in (f ) = S I (f )/ df S I (f ). Specifically, we find for homogeneous excitatory and inhibitory neurons, where only the excitatory population carries a signal (ε E = ε > 0, ε I = 0): n γ−1 , hence ρ n ∝ n γ−1 , 1 n f −1 1 .

A.5 Inverse Gaussian input
To model non-Poissonian input spike trains with a CV smaller or larger than unity, we assume presynaptic spike trains to be renewal processes with an inverse Gaussian ISI density: This distribution is parametrized by the rate ν and the coefficient of variation C V . The spike train power spectrum is given by (Stratonovich 1967) S IG (f ; ν, C V ) = 1 − |P (f ; ν, C V )| 2 |1 −P (f ; ν, C V )| 2 , whereP (f ; ν, C V ) is the Fourier transform of the inverse Gaussian, Eq. (140), given bỹ P (f ; ν, C V ) = exp 1 − 1 − 4πiC 2 V f /ν C 2 V (Holden 1976;Bulsara et al 1994;Chacron et al 2005). The power spectrum of the synaptic current is then, according to Eq. (27): A.6 Synapses with short-term plasticity Following (Dittman et al 2000;Merkel and Lindner 2010), we use a deterministic model for the synaptic dynamics. A presynaptic spike at the th synapse causes a postynaptic response with the history-dependent weight A (t) = F (t)D (t), where facilitation is governed by F (t) = F 0 + ((1 − F 0 ) −1 + F −1 C, (t)) −1 , F C, (t) = −F C, (t)/τ F + ∆ · x (t), and where depression dynamics obeẏ Here, x (t) is the presynaptic spike train, F 0 is the intrinsic release probability, ∆ quantifies the strength of facilitation and τ F (τ D ) is the timescale of facilitation (depression). Above, we have compared the limit cases of pure facilitation and pure depression, which corresponds to setting D (t) ≡ 1 (pure facilitation) or F (t) ≡ F 0 (pure depression).

A.7 Nonlinear integrate-and-fire models
Using the PRC, Eq. (18), the deviation of the ith ISI from its mean is given by This allows us to compute the serial correlations: If C in (t) does not change rapidly, we can approximate t k ≈ k T in the last equation. Hence dt Z(t)e 2πif t . From this follow Eqs. (19) and (20).
The sum of all SCCs evaluates to ∞ n=1 ρ n = 1 2 whereZ (0) is equal to the average PRC obtained by averaging over the interval [0, T ]. Furthermore, using Eq. (126) we obtain the long-time limit of the Fano factor or, equivalently, the spike train power spectrum in the zero-frequency limit as follows lim t→∞ F (t) = T lim f →0 S out (f ) = σ 2Z (0) 2 T S(0).
This equation implies that for any neural oscillator driven by weak colored noise, the spike count diffusion coefficient D eff = S out (0)/2 is proportional to the intensity D η = S(0)/2 of the input noise, irrespective of the nonlinearity.