Fluctuations and information filtering in coupled populations of spiking neurons with adaptation

Finite-sized populations of spiking elements are fundamental to brain function, but also used in many areas of physics. Here we present a theory of the dynamics of finite-sized populations of spiking units, based on a quasi-renewal description of neurons with adaptation. We derive an integral equation with colored noise that governs the stochastic dynamics of the population activity in response to time-dependent stimulation and calculate the spectral density in the asynchronous state. We show that systems of coupled populations with adaptation can generate a frequency band in which sensory information is preferentially encoded. The theory is applicable to fully as well as randomly connected networks, and to leaky integrate-and-fire as well as to generalized spiking neurons with adaptation on multiple time scales.


I. INTRODUCTION
Multi-scale modeling of complex systems has led to important advances in fields as diverse as complex fluid dynamics, chemical biology, soft matter physics, meteorology, computer science, and neuroscience [1][2][3][4][5][6]. In these approaches, mathematical methods such as meanfield theories and coarse-graining provide the basis to link properties of microscopic elements to macroscopic variables. In many cases, macroscopic variables fluctuate due to a finite number of microscopic elements. For instance, in the brain, neurons can be grouped into populations of 50 to 1000 neurons [7] with similar properties [8,9]. Fluctuations of the global activity of such populations are not captured in classical mean-field theories [10,11], which assume infinite system size. Here we put forward a theory for the fluctuating macroscopic activity in networks of pulse-coupled elements occurring in neuronal networks [12], queuing theory [13] and synchronizing fireflies [14].
Adaptation characterized by a reduced response of a neuron to slow compared to fast inputs is a wide-spread phenomenon in the brain and has important implications for signal processing [31][32][33][34] and the spontaneous activity of single neurons [35][36][37]. On the population level, adaptation has been recently analyzed using a quasi-renewal (QR) theory [38]. The QR framework uses tools of renewal point process theory [39] to treat neurons with arbitrary refractoriness. In particular, the dynamics of the population activity is determined by an integral equation [11,40,41]. These studies, however, have assumed an infinitely large population.
Here we present a theory for the interaction of finitesized populations of adapting neurons. The theory is valid for the broad class of neuron models that can be approximated by a QR point process. This includes integrate-and-fire (IF) as well as GLM neurons, for which parameters can be reliably extracted from experimental data [9,12,42].
Based on this theory we analyze information filtering (noise shaping) in neuronal populations. We show that in a single population no noise shaping occurs, but that in coupled populations, band-pass-like noise shaping is possible due to adaptation and connectivity.
This article is structured as follows: First, we present the general dynamics of the population activity and its fluctuations. We then describe how randomly connected, adapting neurons can be treated in this framework. For fluctuations about a stationary state, we linearize the dynamics and compute the spectral density of the population activity. We then determine its coherence with external input signals and quantify information transmission.

II. RESULTS
A. Dynamics of globally coupled renewal models.
Our main quantity of interest is the population activity is the spike train of neuron i with spike times t k i and N denotes the number of neurons. In experiments or simulations, the measured activityĀ(t) would be determined by temporal filtering of the population activity, i.e.Ā(t) =´∞ −∞ A(t − s)f (s)ds with a normalized filter function f (s) with finite support. Below we will use a rectangular filter f (s) = θ(s)θ(∆t − s)/∆t, where θ(s) is the Heaviside step function.
To determine the fluctuation statistics of A(t), we generalize the integral equation of an infinite population [11] to large but finite N . Let us first consider a homogeneous population of all-to-all connected renewal neurons. In this case, the spikes of each neuron occur with an instantaneous rate or hazard function ρ H (t,t), which only depends on its last spike timet ≤ t and the synaptic input determined by the history H(t) = {A(t )} t <t of the population activity. Note that for uncoupled stationary networks, the hazard reduces to ρ H (t,t) = ρ(t −t) as it should be for a renewal model. The probability density of the next spike time t givent is given by Our approach is to use the Gaussian approximation for large N , i.e. we calculate the first-and second-order statistics of A(t) as a functional of its past activity (see Appendix B). However, the dynamics of A depends on its own history H(t) and on the occupation density of refractory states t−t across the population [43]. Thus, we have to average over the possible refractory states consistent with a given history H. To perform the Gaussian approximation, the dynamics of the full system is coarse-grained by discretizing time with a small time step ∆t which is still large enough to include many spikes of the population. For large N , the number of neurons that fire in the time bin at t and had their last spike in bint is a Gaussian random number with mean and variance n 0 (t)P H (t,t)∆t, where n 0 (t) = N´t +∆t t A(s) ds is the past spike count at t < t. Summing overt and treating ∆t as macroscopically infinitesimal, we find the conditional mean activity (see Appendix B) which is equal to the population integral [11] for the infinite system. For finite N , A(t) will be of the form where the deviation δA(t) has zero mean and a diverging standard deviation a(t)/(N ∆t) because the variance of the spike count in [t, t + ∆t] is given by N ∆t a(t). Importantly, δA(t) cannot be described by a white noise process but future values δA(t + τ ), τ > 0, are correlated with δA(t), because they share a common history H(t).
In fact, a neuron that fired its last spike att < t cannot have its next spike at both times t and t+τ , which induces a negative correlation for the deviations at t and t + τ . We find (see Appendix B) for τ ≥ 0 the conditional correlation function where · H(t) denotes the average conditioned on the history of A before t. Thus, the correlation function is in general explicitly time-dependent.
B. Adaptation and random connectivity.
In the presence of adaptation, the instantaneous rate of a neuron depends on all its previous spikes so that it can no longer be described by renewal theory. Here we describe how adapting neurons in networks may still be approximated by a quasi-renewal process. Specifically, we consider a homogeneous population of neurons modeled by the spike-response model with escape-noise [11], also known as GLM [9,12,42], with hazard function ( Fig. 1). That is, neuron i produces a spike in a small time interval [t, ∆t) with probability ρ i (t)∆t. This probability depends on the input potential h i (t) = κ * N j=1 w ij s j (t) + I (t), which is driven by presynaptic spike trains s j (t) (with synaptic weight w ij ) and external input I(t). The membrane filter kernel is given by κ(t) = θ(t − τ s ) exp(−(t − τ s )/τ m ), where τ s and τ m are the synaptic delay and the membrane time constant, respectively. The operation * denotes the convolution (f * g)(t) =´∞ −∞ f (t − s)g(s)ds and θ(t) is the Heaviside step function.
The variable ϑ i defined as ϑ i (t) = (s i * η)(t) can be interpreted as a dynamic firing threshold that is triggered by the neuron's own output spike train s i (t) [9]. Here, η(t) is a feedback kernel that consists of two parts, η = η a + η r : a short-range refractory kernel η r (t) = θ(t) [J r θ(t − τ abs ) exp(−(t − τ abs )/τ m ) + θ(τ abs − t)D] mainly affected by the last spike, and a long-range adaptation kernel η a (t) = J a θ(t) exp(−t/τ a ) that accumulates the spike history on a longer time scale τ a . An absolute refractory period is included in η r (t) by setting it to D = 10 12 for 0 < t < τ abs . Our choice of the kernels corresponds to a leaky IF model with dynamic threshold [44] and a reset by a constant amount −J r after each spike.
The parameter c in Eq. (4) sets a baseline firing rate and δu sets the strength of intrinsic noise ("softness" of threshold). Fits of this model to pyramidal neuron recordings yielded δu ≈ 4mV [45]. Our standard parameter set given below corresponds to an amplitude of single post-synaptic potentials of 0.25mV (excitatory, exc.) and −1.1mV (inhibitory, inh.). In the following, we measure voltage in units of δu, so that δu = 1 in dimensionless units. For the synaptic weights w ij , we use a homogeneous random network as specified in Appendix A, below.
The dependence of the term exp(−ϑ i ) in Eq. (4) which describes the feedback of the neuron's own spiking history s i can be approximated by the explicit contribution of the last spike of the neuron att and the average effect of previous spikes up tot [38]: Here, the average is taken over all previous spike times t k i <t. As shown in [38], this average can be approximated by exp ´t −∞ e −η(t−t ) − 1 s i (t ) dt . Replacing further the firing rate s i (t ) by the population activity A(t ), the threshold ϑ i (t) becomes for all neurons with last spike att. The kernel γ τ (s) = θ(s − τ )(1 − e −η(s) ) represents the effect of adaptation in the quasi-renewal (QR) approximation. Furthermore, for homogeneous random networks and large N , the local field h i caused by synaptic input to neuron i is determined by A and an effective weightw = N −2 i,j w ij [10], hence where J s = Nw. The above steps enable us to treat neuronal adaptation and network coupling in a quasirenewal framework with hazard function t + ∆t), with ∆t = 2ms, black) resulting from 500 randomly connected adapting neurons (4) receiving a common input current (I(t), blue). The theoretical expectation a(t) (green, Eq. 1) based on QR approximation (7), and expected fluctuations (red, one std., a(t) ± (N ∆t) − 1 2 a(t) 1 2 ). At time t, a(t) depends on the actual historyĀ(t ) (black) for t < t. Standard parameters (see Appendix A), except for I(t) as shown and c = 5s −1 . Inset shows mean (green) and std. (red) of deviations δA(t) =Ā(t) − a(t) as a function of √ a, averaged over 25 repetitions of the displayed I(t) (dots) vs. theory (lines).
Note that ρ H (t,t) is identical for all neurons which have fired their last spike att.
In Fig. 2 the population activity of a spiking neural network simulation is compared to the theoretical prediction (2). To evaluate (1) numerically, we iteratively compute S H (t + ∆t,t) = S H (t,t)(1 − ρ H (t,t)∆t), with S H (t, t) = 1, and use P H = ρ H S H . The QR population integral describes the response of the population activity and its fluctuations for stationary, as well as for slowly or rapidly varying inputs.

C. Linearized population dynamics.
The amount of information transmitted and processed in sensory areas of the brain is limited by the fluctuations of the population activities [46][47][48]. Likewise, in decision networks, finite-size induced fluctuations determine the reliability of decisions [49,50]. The spontaneous activity of cortical networks is typically asynchronous, and is believed to underlie cortical information processing [51]. In order to analytically determine the power spectrum of the spontaneous population activity, we linearize the dynamics around the large N limit. To this end, we assume that in the limit N → ∞ and for constant external input I(t) = I 0 the network dynamics has an equilibrium point with activity A 0 corresponding to an asynchronous firing state. With this equilibrium activity we can associate a renewal neuron model that is obtained from the original model, Eq. (4), by replacing h i (t) and ϑ i (t,t) by h 0 = κ * (JA 0 + I 0 ) and ϑ 0 (t −t) = η(t −t) + (γ t−t * A 0 ), respectively. In the following, we will use the subscript "0" to refer to quantities of the associated renewal model. For finite system size, N < ∞, the activity will deviate from A 0 . Through (5) and (6) the fluctuations ∆A(t) = A(t) − A 0 also lead to fluctuations ∆h(t) = h(t) − h 0 (t) and ∆ϑ(t,t) = ϑ(t,t) − ϑ 0 (t −t) , which in turn influence A(t). Our goal is to determine the spectral properties of ∆A(t).
To simplify the derivations, we approximate the QR kernel γ τ by its average over the inter-spike-interval den- Expanding (1)-(3) to first order in ∆A, ∆h and ∆ϑ, yields the linearized stochastic dynamics (see Appendix C) determines the linear response of the expected activity a to a perturbation ∆A. For our model (4), the kernel L in this expression is given by L(t) = θ(t)´∞ 0 ρ 0 (s)S 0 (s+t)ds but L can be derived for most common neuron models [11], or, alternatively, may be estimated from neural recordings. The noise term ξ(t) is stationary Gaussian noise with correlation function for all τ ; cf. (3). Eq. (9) shows that the population activity in the stationary state is a Gaussian process with memory, where finite-size fluctuations are described by the colored noise ξ(t).

D. Fluctuations in coupled populations.
Let us now turn to K populations consisting of N = (N 1 , . . . , N K ) neurons. Parameters of neurons and coupling are homogeneous within each population but may differ between one group and the next. To incorporate network coupling, the network input (6) where J is the coupling matrix and A(t) is a vector of population activities. For each population the dynamics are given by (9) but Q(t) is now a K × K matrix of coupling kernels. Using the Fourier transformf (ω) =´∞ −∞ f (t)e −iωt dt, this matrix can be written asQ(ω) =P 0 + iωA 0L (KJ −G), where the matrices P 0 , A 0 , L, K, G are defined as the diagonal matrices of the vectors P 0 , A 0 , L, κ, γ, respectively. The power spectrum, defined as the Fourier transform of the Here, † denotes the adjoint matrix (conjugate transpose). It is instructive to rewrite this expression in terms of the power spectrum of the associated renewal model [52] Here is the diagonal matrix containing the linear response functions of the associated renewal models with respect to current perturbations I(t) [11]. Eq. (10) shows that finite-size fluctuations are characterized by the renewal spectrumC 0 (ω), shaped by recurrent input (via J andK) and adaptation (viaG), and reduced by the factor N −1 . There are several known limit cases: first, for vanishing adaptation,G = 0, we recover the linear response result of [28][29][30] for networks of white-noise driven IF neurons. Our formula shows that adaptation appears as an additional diagonal term in the effective coupling matrix J −K −1G , and hence can be interpreted as an inhibitory self-coupling. Second, if both adaptation and recurrent connections vanish,G = 0 and J = 0, we arrive atC A (ω) = N −1C 0 (ω) because the superposition of independent spike trains does not change the shape of the power spectrum. Third, our result also includes the frequently employed Hawkes process [18][19][20], which is recovered for a constant single neuron spectrum C 0 = A 0 and vanishing adaptation,G = 0. For a comparison of our result to simulations, see Fig. 3, which will be discussed below. In section II F we make use of Eq. (10) to quantify information filtering in neural populations (Fig. 4). A comparison to the special cases of Eq. (10) described above is shown in Fig. 5.
In Fig. 3(a)-(b) the spectral density is shown compared to simulations, where the frequency equals ω/(2π). The spectra are well described by the novel theory, which captures refractoriness, recurrent feedback and the reduction of power at low frequencies due to adaptation. The latter arises from negative correlations between ISIs typical for adapting neurons [37]. Interestingly, this purely nonrenewal effect is well accounted for by our quasi-renewal theory. Since adaptation effects are most prominent at low frequency, we examined the dependence of the power on the model parameters for ω → 0 (Fig. 3(c)-(d)). Our theory describes the simulations well across the studied parameter range. Neuronal networks in the brain are subject to external influences, either due to sensory input or ongoing activity in other brain areas. How do neuronal populations respond to small, time-dependent input cur- rents I(t) = (I 1 , . . . , I L ) with spectral densityC I (ω)?
To answer this question, we proceed as before and linearize Eqs. (1)-(3) with respect to the small fluctuations ∆h k (t) = (κ k * (J∆ A + M I) k )(t) of the local field. Here, we restrict our analysis to independent inputs I, such thatC I is diagonal with entriesC I,i , but also included a K ×L mixing matrix M, which allows us to model shared input. The resulting spectral density is given by the sum Thus, additional fluctuations due to the stimulus I are shaped by both the single neuron filterR 0 M and the network and adaptation filterB (10), combining the effects of recurrent connectivity and adaptation.
F. Information transmission.
Our theory allows us to quantify the transmission of information from external input signals I(t) through a system of coupled neural populations. The coherence between the signal j and the activity of population i can be regarded as a frequency resolved measure of information transmission. Information theory [53][54][55] states that the mutual information rate is bounded from below by −´∞ 0 log 2 [1 − Γ ij (ω)] dω 2π . Since adaptation attenuates the response to slowly changing signals, one might expect that it also attenuates low frequency information content. For a single population, however, this is not the case, but instead the coherence is low-pass, i.e. it monotonically decreases for increasing frequency [29]. Here we show that in coupled populations of adapting neurons, coherences can be nonmonotonic allowing the neural circuit to preferentially encode information in certain frequency bands. Put differently, a multi-population setup can realize an information filter.
Using Eq. (12), we find the general form of the coherence matrix: In this expression, the numerator represents the contribution of the signal I j (t) to the power spectrum of population i. This effective signal power is divided by the total power spectrum of population i, which consists of direct (k = i) and indirect (k = i) sources of variability. Both sources contain internally generated noise due  to finite size N k as well as signal power. However, the diagonal elements (k = i) of the shaping matrices can be much stronger than the off-diagonal elements (k = i) depending on the coupling matrix J. Therefore, we expect that the direct source of variability dominates in the denominator.
If there is only one population and signal (K = 1, L = 1), the term |B 11 (ω)| 2 occurs in both numerator and denominator and cancels. Thus, coupling and adaptation do not shape the coherence in a single population. Furthermore, the signal term |R 0,1 (ω)| 2 M 11CI,1 (ω) is matched in both numerator and denominator, which leads to a flat coherence at frequencies where the signal dominates the finite N noise. At high frequencies, the neural response amplitude |R 0 | 2 decays due to the leaky membrane, but the spontaneous spectrumC 0 has a constant high-frequency limit equal to A 0 . One therefore typically observes a low-pass like information transfer characteristics of single neurons or populations [29,56,57].
For several populations (K > 1), however, we can distinguish two cases: If the signal is read out at a receiving population, M ij = 0, the signal power in the numerator is matched by the dominating direct signal power in the denominator, and hence the shaping of the signal power cancels. In contrast, if read out at a different population, M ij = 0, the signal j contributes only indirectly to the power spectrum of population i via synaptic connections. Thus, we expect that the shape of signal power and power spectrum (i.e. numerator and denominator, respectively) is generally different if the transmission path involves multiple populations.
As an example of this mechanism, we show a feedforward chain of excitatory and inhibitory populations (Fig. 4, K = 6, L = 1). In the first layer, the effective signal power is reduced at low frequencies because of adaptation and inhibitory feedback. However, the power spectrum of A exc,1 (t) is dominated by the same signal power and hence exhibits a very similar reduction of lowfrequency power. Consequently, the coherence (being the ratio of these two spectra) is rather flat at low frequencies and shows a decay at higher frequencies (Fig. 4,  red lines). This low-pass characteristics changes at later stages in the chain (green & blue): The signal term is increasingly more shaped by adaptation and coupling properties, whereas the noise spectrum changes less. As a result, the coherence shows a maximum at a finite frequency (Fig. 4, blue lines). This band-pass structure becomes more pronounced from layer to layer, representing a form of information filtering. Coherence functions with band-pass characteristics have been observed in neurons postsynaptic to electro-receptor afferents in electric fish [58,59].

III. CONCLUSIONS
We have shown that fluctuations in finite-sized networks of spiking neurons are captured by a colored noise term added to the population integral equation of the infinite system. Our approach yields spectral densities of the population activity in randomly or fully connected multi-population networks which are in excellent agreement with simulation results. Our quasi-renewal theory includes refractory effects and adaptation on multiple time scales. In contrast to earlier treatments of neuronal refractory effects in population dynamics [17,19,20] or linear response formulas for adaptive neurons [60], the QR population integral, derived directly from the neuron model definition, captures the time-dependent, nonlinear dynamics and adaptation of neural population activity.
We applied our theory to information filtering by coupled populations of spiking neurons with adaptation. We showed that, although impossible in single populations due to a cancellation of signal and noise terms of the coherence, coupled populations can filter information through adaptation mechanisms and neuronal interactions. This mechanism might be exploited in the layered structure of cortical circuits, or in sensory systems of insects where signals traverse a sequence of nuclei.
In this paper we treated populations of point neurons with static synapses and applied linear response theory. How to generalize our theory to incorporate effects of nonlinear dendritic integration, spike-synchrony detection and short-term synaptic plasticity, which all contribute to information filtering [61][62][63], is an important question that merits further investigation. Nonetheless, due to the versatility of GLM models, our theory already provides a useful tool for interpreting neural data at the population level. For example, our theory suggests that in in-vitro experiments with optogenetically evoked input currents and simultaneous measurements of neural activity [64], system parameters may be identified based on the relation of the spectra, Eq. (11). Moreover, largescale neural systems can now be analyzed as coupled populations of model neurons with single-cell parameters extracted from experiments, and simulated using a muti-scale approach. M. D. and T. S. contributed equally to this work.

Appendix A: Network simulations
We compare our theoretical results to simulations of networks of excitatory and inhibitory neurons defined by (4). In the case K = 1 (Fig. 3(a)), we use N = N and J = J s , in the case K = 2 ( Fig. 3 Generally, neurons of population i receive synapses from a random subset of pN j neurons of population j, each with synaptic weight w ij = J ij /(pN j ) and delay τ s . Unless the connection probability p is 1, self-connections are excluded. Networks were simulated for 2 · 10 4 s using NEST [65] (neuron model pp_psc_delta, temporal resolution 2ms). Standard parameters, unless indicated otherwise: N = 500, c = 10s −1 , τ m = 0.01s, τ a = 0.3s, ∆t = τ s = τ abs = 2ms, J r = 3, J a = 1, J s = 5, p = 0.2, I 0 = 0,C I = 0. In the exc.-inh. network, for the inh. neurons which typically show little adaptation [9] we deactivated adaptation by setting J a = 0 and c = 5s −1 . While it is possible to theoretically approximate the stationary interval distribution P 0 by searching for a selfconsistent rate A 0 as described in [38], here we use P 0 from simulated inter-spike-intervals of each population. From the measured P 0 (t) we derive A 0 = 1/´∞ 0 tP 0 (t)dt, S 0 (t) =´∞ t P 0 (t )dt and ρ 0 (t) = P 0 (t)/S 0 (t).

Appendix B: Detailed derivation of Eq. (3)
The aim is to find a dynamical equation for the population activity where n 0 (t) is the total number of spikes in the interval [t, t + ∆t]. More generally, we define n k (t), k ∈ Z, as the total number of spikes in [t − k∆t, t − (k − 1)∆t], i.e. the activity k time bins in the past. It is useful to consider

of neurons with last spike in
[t − k∆t, t − (k − 1)∆t] and next spike in [t, t + ∆t] ρ(t,t) hazard function: rate at t given last spike att S(t,t) survivor function: probability of no spike in [t, t] P (t,t) inter-spike-interval density: probability density of next spike at t given last spike att; P = ρ · S furthermore the total number of neurons that spiked in the time bin [t − k∆t, t − (k − 1)∆t], k = 1, 2, . . . , but had no further spike until time t. Let us denote this number by m k (t) (Fig. 6 and Table I). The number of neurons that spike in [t, t + ∆t] and had their last spike in the bin [t − k∆t, t − (k − 1)∆t] shall be denoted by δm k (t). These neurons decrease the number m k (t) of neurons from group k that had survived until time t + ∆t in the next time step, i.e.
The total number of spikes at time t, n 0 (t), is the sum over all possible last spike times, hence We will now express the activity n 0 (t) in terms of the past activity H t = {n k (t)} k=1,2,... , using the Gaussian approximation. This requires to compute the mean and correlation function of n 0 (t) given the past values n k (t), k = 1, 2, . . . . In the following, the averaging bracket · has to be understood as the conditional average · Ht = · {n k (t)} k=1,2,... , i.e. we will omit the conditioning subscript for simplicity. Although n k (t), the total number of spikes in bin t−k∆t, is fixed, the number m k (t) of neurons that had their last spike in bin t − k∆t is variable. It is this variability that we will average over (This corresponds to a statistical ensemble of populations that all have an identical history of population activity n k (t), k = 1, 2, . . . .).
Suppose we know the value m k (t) of the group of neurons with their last spike in [t − k∆t, t − (k − 1)∆t]. Then the expected number of spikes from that group in the next interval is where ρ(t,t) is the hazard function of the neurons (instantaneous rate at time t given last spike att), and x y denotes the expectation of x conditioned on y (in addition to the overall condition of a fixed history H t = where we have used (B4). We now use this result to calculate the expected number of spikes in the interval [t, t + ∆t]. Averaging over (B3) yields The average number of neurons that fired their last spike in [t − k∆t, t − (k − 1)∆t] and survived up to t can be expressed using the survival probability S(t,t) = exp(−´t t ρ(t ,t)dt ) as follows: We can now take the limit ∆t → 0 in (B6) and find where P (t,t) = ρ(t,t)S(t,t) is the inter-spike-interval density. Eq. (B8) is equivalent to Eq. (1).
To obtain the correlation function we can write for q ∈ Z n 0 (t)n 0 (t + q∆t) = ∞ k,l=1 Here, the spike numbers δm k (t) and δm l (t + q∆t) that refer to different groups k and l − q are uncorrelated. Correlations only arise for δm k (t) and δm k+q (t + q∆t), i.e. spikes that refer to the same group in the past. Thus, n 0 (t)n 0 (t + q∆t) = ∞ k,l=1 δm k (t) δm l (t + q∆t) δm k (t) δm k+q (t + q∆t) , (B10) so that the covariance is where ∆n 0 (t) = n 0 (t) − n 0 (t) . Therefore we need to compute δm k (t)δm k+q (t + q∆t) . (B12) To this end, let us consider the cases q = 0 and q > 0 separately. For q = 0 and large N , the number of neurons that spike in [t, t + ∆t] and had their last spike at t − k∆t is a Poisson variable with mean and variance ∆t P (t, t − k∆t) n k (t). Thus (B12) becomes For q > 0, we employ (B2) twice and obtain In order to evaluate each of these four correlators, we note that the probability that a neuron from group k "survives" until time t + q∆t given that it survived until time t is S(t + q∆t, t − k∆t)/S(t, t − k∆t) according to Bayes law. Thus, out of the m k (t) neurons that survived until time t, on average also survive until t + q∆t. Therefore, the correlator for 0 ≤ l < q can be written as Applying this result to (B14), we obtain How can we calculate the second moment of m k (t)? Recall that m k (t) is the part of the n k (t) neurons firing in bin t − k∆t that survived until time t. Thus m k (t) can be regarded as a binomially distributed random number with n = n k (t) trials and survival probability p = S(t, t− k∆t). This random number has mean np and variance np(1 − p). Hence, the second moment reads Likewise, because n k+1 (t+∆t) = n k (t). Inserting (B17) into (B16) we find Here, we have identified the derivative d /dtS(t,t) = −P (t,t). Note that this expression is of order O(∆t 3 ), whereas δm k (t) δm k+q (t+q∆t) is of order O(∆t 4 ). So we can neglect the second term on the right-hand side of Eq. (B11).
Putting all together, we find 1 ∆t 2 ∆n 0 (t)∆n 0 (t = t + q∆t) δm k (t)δm k+q (t + q∆t) Thus, using (B8), we arrive at the final result where δA(t) is Gaussian with conditional correlation function δA(t)δA(t + τ ) = N −1 δ(τ )ˆ∞ −∞ P (t, t )A(t ) dt − N −1ˆ∞ −∞ P (t + τ, t )P (t, t )A(t ) dt (B21) for τ ≥ 0. In the latter expression we extended the limits of integration to infinity, which is possible if we assume P (t, t ) = 0 for t < t . Equation (B21) tells us that the noise correlation function consists of two parts: a white (δ-correlated) part and a negative correlation due to neural refractoriness. Intuitively, since δm k (t) and δm k+1 (t + ∆t) share the same number of available neurons m k (t), a positive fluctuation of the number of spikes in the bin t of the neurons of group k reduces the number of neurons with last spike in bin k more than on average. Thus, the number of neurons of group k that can still fire in time bin t + ∆t is smaller than on average, which explains the negative correlations (Fig. 6).
Appendix C: Detailed derivation of Eq. (9) We aim to linearize the population integral a(t) = t −∞ P (t,t)A(t) dt, Eq. (1), around an equilibrium point A 0 , with small fluctuations ∆A(t) = A(t) − A 0 with a mean of zero. To this end, let us first note that the hazard function Eq. (7) then can be written as ρ(t,t) = ce h0−ϑ0(t−t) e ∆h(t)−∆ϑ(t,t) = ρ 0 (t −t)e ∆g(t,t) .