Echoes in correlated neural systems

Correlations are an established feature of neural activity [1] and evidence is increasing that they can provide insights into information processing in the brain [2]. The temporal relationship between the activity of pairs of neurons is described by correlation functions. Their shape has early been related to the direct coupling between neurons and to the common input shared by pairs of neurons. On the one hand, correlations may limit the signal-to-noise ratio of population rate signals [3], on the other hand they have shown to increase the amount of information available to unbiased observers [4]. Furthermore, synchronous neural activity has been proposed to bind elementary representations into more complex objects [5] and experimental evidence for such a correlation code is provided by task-related modulation of synchrony in the primary visual cortex [6] and the motor cortex [7].

The small magnitude [8] of pairwise correlations in the asynchronous irregular state [9] of the cortex has recently been related to the balance between excitation and inhibition in local networks [10, 11] and inhibitory feedback was identified as a general mechanism of decorrelation [12]. However, a quantitative theory explaining the temporal shape of correlation functions in recurrently impulse-coupled networks of excitatory and inhibitory cells remained elusive.

Assuming random connectivity with identical numbers and strengths of incoming synapses per neuron, as illustrated in figure 1, suggests by mean field arguments [13, 14] that the resulting activity of two arbitrarily selected neurons, and hence the power spectra of activities averaged over excitatory or inhibitory neurons, should be the same. Direct simulations, however, exhibit different power spectra for these sub-populations [15]. A similar argument holds for the covariance c_ff between the two neurons: if the covariance c between any pair of inputs is known, the covariance between their outgoing activity c_ff is fully determined [15–24]. By self-consistency, as both neurons belong to the same recurrent network, one concludes that c_ff = c. In particular the covariance averaged over excitatory pairs should be identical to the corresponding average over inhibitory pairs, which is in contrast to direct simulation (figure 1(b)). In this work, we elucidate why this mean field argument for covariances fails and derive a self-consistency equation for pairwise covariances in recurrent random networks which explains the differences in the power spectra and covariances.

Zoom In Zoom Out Reset image size

Theories for pairwise covariances have been derived for binary neuron models [11, 25] and for excitatory stochastic point process models [26]. However, the lack of either inhibition or delayed pulsed interaction limits the explanatory power of these models. A theory for networks of leaky integrate-and-fire (LIF) model neurons [27] is required, because this model has been shown to well approximate the properties of mammalian pyramidal neurons [28] and novel experimental techniques allow to reliably assess the temporal structure of correlations in the cortex [29]. Moreover, the relative timing of action potentials is the basis for models of synaptic plasticity [30], underlying learning in biological neural networks. Analytical methods to treat population fluctuations in spiking networks are well advanced [31] and efficient hybrid analytical–numerical schemes exist to describe pairwise covariances [32]. Here we present an analytically solvable theory of pairwise covariances in random networks of spiking leaky integrate-and-fire model neurons with delayed pulsed interaction in the asynchronous irregular regime.

We consider recurrent random networks of N excitatory and γN inhibitory leaky integrate-and-fire model neurons receiving pulsed input (spikes) from other neurons in the network. Each neuron has K = pN incoming excitatory synapses independently and randomly drawn from the pool of excitatory neurons, and γK = γpN inhibitory synapses (a homogeneous Erdős–Rényi random network with fixed in-degree). An impulse at time t arrives at the target neuron after the synaptic delay d and elicits a synaptic current I_i that decays with a time constant τ_s and causes a response in the membrane potential V_i (with a time constant τ_m) proportional to the synaptic efficacy J (excitatory) or −gJ (inhibitory), respectively. The coupled set of differential equations governing the subthreshold dynamics of a single neuron i is [33]

$$\begin{eqnarray} \tau_{\mathrm{m}}\frac{\mathrm{d}V_{i}}{\mathrm{d}t} & = & -V_{i}+I_{i}(t),\nonumber \\ \tau_{\mathrm{s}}\frac{\mathrm{d}I_{i}}{\mathrm{d}t} & = & -I_{i}+\tau_{\mathrm{m}}\sum_{j=1,j}^{N}J_{ij}s_{j}(t-d), \end{eqnarray} \tag{ 1 }$$

where the membrane resistance was absorbed into J_ij. If V_i reaches the threshold V_θ at time point tⁱ_k the neuron emits an action potential and the membrane potential is reset to V_r, where it is clamped for the refractory time τ_r. The spiking activity of neuron i is described by this sequence of action potentials, the spike train $s_{i}(t)=\sum _{k}\delta (t-t_{k}^{i})$.

The activity of a given neuron i depends on the history of the other neurons' activities s(t) = (s₁(t),...,s_N(t))^T in the network, so formally we can consider the spike train of neuron i as a functional of all other spike trains. The time-averaged covariance matrix expresses these interrelations and is defined as $\bar {\mathbf {c}}(\tau )=\langle \mathbf{s}(t+\tau )\mathbf {s}^{\mathrm {T}}(t)\rangle _{t}-\mathbf {r}\mathbf {r}^{\mathrm {T}}$, where r = 〈s〉_t is the vector of time-averaged firing rates. The diagonal contains the autocovariance functions (diagonal matrix a(τ)) which are dominated by a δ-peak at zero time lag and for τ ≠ 0 exhibit a continuous shape mostly determined by refractoriness, the inability of the neuron to fire spikes in short succession due to the voltage reset, as shown in figure 2(d). The off-diagonal elements contain the cross covariance functions c(τ) that originate from interactions. We therefore decompose the covariance matrix into $\bar {\mathbf {c}}=\mathbf {a}+\mathbf {c}$. A basic property of covariance matrices is the symmetry $\bar {\mathbf {c}}(\tau )=\bar {\mathbf {c}}^{\mathrm {T}}(-\tau )$, so we only need to consider τ > 0 and obtain the solution for τ < 0 by symmetry. Each spike at time t' can influence the other neurons at time t > t'. Formally we express this influence of the history up to time t on a particular neuron i in the network as s_i(t,{s(t')|t' < t}), which is a functional of all spike times until t. In the asynchronous state of the network [9] and for small synaptic amplitudes J_ij a single spike of neuron j causes only a small perturbation of s_i. The other inputs to neuron i effectively act as additional noise. We may therefore perform a linear approximation of j's direct influence on neuron i and average over the realizations of the remaining inputs s\s_j, as illustrated in figure 2(a). This linearization and the superposition principle lead to the convolution equation

$$\begin{eqnarray} \langle s_{i}(t)|s_{j}\rangle_{{\mathbf{s}{\backslash} s_{j}}} & = & r_{i}+\int_{-\infty}^{t}h_{ij}(t,t^{\prime})(s_{j}(t^{\prime})-r_{j})\, \mathrm{d}t^{\prime}\nonumber\\ &=& r_{i}+[h_{ij}\ast(s_{j}-r_{j})](t), \end{eqnarray} \tag{ 2 }$$

where we define as the linear response kernel $h_{ij}(t,t^{\prime })=\langle \frac {\delta s_{i}(t)}{\delta s_{j}(t^{\prime })} \rangle _{\mathbf {s}{\backslash } s_{j}}$ the functional derivative of s_i(t) with respect to s_j(t'), formally defined in appendix A. The kernel h_ij quantifies the effect of a single spike at time t' of neuron j on the expected density of spikes of neuron i at time t by a direct synaptic connection from neuron j to neuron i, earlier introduced as the 'synaptic transmission curve' [34]. This density vanishes for t < t' due to causality. For the stationary network state studied here, the kernel further only depends on the time difference τ = t − t'. In a consistent linear approximation there are no higher order terms, so the effects of two inputs s_j and s_k superimpose linearly. In general, the response of a cell is typically supra-linear in the number of synchronously arriving excitatory spikes [34–36]. If, however, the network state to be described is sufficiently asynchronous, as is the case here, a linear approximation is adequate. Using this linear expansion and the definition of the covariance matrix, the off-diagonal elements of the covariance matrix fulfil a linear convolution equation

$$\begin{eqnarray} \mathbf{c}(\tau) & = & [\mathbf{h}\ast\left(\mathbf{a}+\mathbf{c}\right)](\tau)\quad\textrm{for }\tau>0. \end{eqnarray} \tag{ 3 }$$

The equation is by construction valid for τ > 0 and needs to be solved simultaneously obeying the symmetry condition c(t) = c(−t)^T. Of up to linear order in the interaction, equation (3) is the correct replacement of the intuitive self-consistency argument sketched in figure 1.

Zoom In Zoom Out Reset image size

In order to relate the kernel h_ij to the leaky integrate-and-fire model, we employ earlier results based on Fokker–Planck theory [33]. For small synaptic amplitudes J ≪ V_θ − V_r and weak pairwise covariances the summed synaptic input $\tau _{\mathrm {m}}\sum _{j=1}^{N}J_{ij}s_{j}(t-d)$ can be approximated as a Gaussian white noise with mean $\mu _{i}=\tau _{\mathrm {m}}\sum _{j}J_{ij}r_{j}$ and variance $\sigma _{i}^{2}=\tau _{\mathrm {m}}\sum _{j}J_{ij}^{2}r_{j}$ for neurons firing with rates r_j and Poisson statistics. For short synaptic time constants τ_s ≪ τ_m the stationary firing rate r_i(μ_i,σ²_i) (A.1) depends on these two moments [33]. In a homogeneous random recurrent network the input to each neuron is statistically the same, so the stationary rate r_i = r is identical for all neurons. It is determined by the self-consistent solution of r(μ_i,σ²_i), taking into account the dependence of μ_i and σ²_i on the rate r itself [9].

The integral $w_{ij}=\int _{0}^{\infty }h_{ij}(t)\,\mathrm {d}t$ of the response kernel is equivalent to the dc susceptibility $w_{ij}=\frac {\partial r_{i}}{\partial r_{j}}$ [37] and has earlier been termed 'asynchronous gain' [34]. The approximation is second order in the synaptic amplitude w_ij = αJ_ij + βJ²_ij. The first order term originates from the dependence of μ_i on r_j, the second order term stems from the dependence of σ²_i on r_j (A.3).

Throughout this work we choose the working point of the neurons in the network such that the firing of the cells is driven by strong fluctuations and the mean membrane potential is close to threshold, in order to have irregular firing. This is achieved by appropriate choices of the external driving Poisson sources, as described in appendix E. For the small networks considered here, it is realistic to assume that about 50% of the inputs to a neuron come from outside the local network [38]. Figure 2(b) shows the deflection of the firing rate from the baseline caused by an impulse in the input averaged over many repetitions. For sufficiently strong fluctuations σ_i, the quadratic term in J_ij is negligible, as seen from figure 2(c). For smaller membrane potential fluctuations σ_i we expect the linear approximation to be less accurate.

The kernel shows exponential relaxation with an effective time constant τ_e that depends on the working point (μ_i,σ_i) and the parameters of the neuron, which is obtained by a least squares fit of a single exponential to the simulated response in figure 2(b) for one particular amplitude J_ij. We therefore approximate the response h_ij(t) as

$$\begin{eqnarray} w_{ij}h(t) & = & \Theta(t-d)\;\frac{w_{ij}}{\tau_{\mathrm{e}}}\,\mathrm{e}^{-\frac{t-d}{\tau_{\mathrm{e}}}}, \end{eqnarray} \tag{ 4 }$$

where d is the synaptic delay and Θ the Heaviside step function.

In experiments covariance functions are typically averaged over statistically equivalent pairs of neurons. Such averages are important, because they determine the behavior of the network on the macroscopic scale of populations of neurons. We therefore aim at a corresponding effective theory. Assuming identical dynamics (1), all neurons have, to good approximation, the same autocovariance function a(t) and response kernel h(t). So each incoming excitatory impulse causes a response w h(t), an inhibitory impulse −gw h(t). We define the covariance function averaged over all pairs of excitatory neurons as $c_{\mathrm {e}\mathrm {e}}(\tau )=\frac {1}{N^{2}}\sum _{i,j\in \mathcal {E},i\neq j}c_{ij}(\tau )$ (setting N(N − 1) ≃ N² for N ≫ 1), where $\mathcal {E}$ denotes the set of all excitatory neurons. The pairings c_ei,c_ie, and c_ii are defined analogously. Inserting equation (3) into the average c_ee(τ), the first term proportional to the autocovariance a(t) only contributes if neuron j projects to neuron i. For fixed i, there are K such indices j, so the first term yields $\sum _{i,j\in \mathcal {E},i\neq j}h_{ij}\ast a_{j}=NKw\, h\ast a$. The second sum $\sum _{i,j\in \mathcal {E},i\neq j}h_{ik}\ast c_{kj}$ can be decomposed into a sum over all intermediate excitatory neurons $k\in \mathcal {E}$ and over all inhibitory neurons $k\in \mathcal {I}$ projecting to neuron i. Replacing the individual covariances by their population average, c_ee and c_ie, respectively, and considering the number of connections and their amplitude we obtain $wNK\, h\ast \left (c_{\mathrm {e}\mathrm {e}}-\gamma g\, c_{\mathrm {i}\mathrm {e}}\right )$, with γ = N_i/N_e and g = w_i/w_e the relative number of inhibitory neurons and the relative linearized inhibitory synaptic amplitude. Similar relations hold for the remaining three averages, so we arrive at a two-by-two convolution matrix equation for the pairwise-averaged covariances for τ > 0

$$\begin{eqnarray} \fl \mathbf{c}(\tau) &=& \left[h\ast\mathbf{M}\mathbf{c}\right](\tau)+\mathbf{Q}\left[h\ast a\right](\tau)\nonumber\\ \textrm{with }\mathbf{M} &=& Kw\left(\begin{array}{@{}cc@{}} 1 & -\gamma g\\ 1 & -\gamma g \end{array}\right),\quad\mathbf{Q}=\frac{Kw}{N}\left(\begin{array}{@{}cc@{}} 1 & -g\\ 1 & -g \end{array}\right),\nonumber\\ \textrm{and }\mathbf{c}(\tau) &=& \left(\begin{array}{@{}cc@{}} c_{\mathrm{e}\mathrm{e}}(\tau) & c_{\mathrm{e}\mathrm{i}}(\tau)\\ c_{\mathrm{i}\mathrm{e}}(\tau) & c_{\mathrm{i}\mathrm{i}}(\tau) \end{array}\right). \end{eqnarray} \tag{ 5 }$$

The convolution equation only holds for positive time lags τ. For negative time lags it is determined by the symmetry c(−τ) = c^T(τ). The solution of this equation can be obtained by an extension of the method used in [26] employing Wiener–Hopf theory [39] to the cross spectrum $\mathbf {C}(\omega )=\int _{-\infty }^{\infty }\mathbf {c}(t)\mathrm {e}^{-\mathrm {i}\,\omega t}\, \mathrm {d}t$ in the frequency domain, as shown in appendix B (here capital letters denote the Fourier transform of the respective lower case letters). With the definition of the propagator $\mathbf {P}(\omega )=\left (\mathbf {1}-\mathbf {M} H(\omega )\right )^{-1}$ the cross spectrum takes the form $\mathbf {C}(\omega )=\mathbf {P}(\omega )\left [D_{+}(\omega )\mathbf {Q}+D_{+}(-\omega )\mathbf {Q}^{\mathrm {T}}-A(\omega )|H(\omega )|^{2}\mathbf {Q}\mathbf {M}^{\mathrm {T}}\right ]\mathbf {P}^{\mathrm {T}}(-\omega ),$ where we split the term H(ω)A(ω) = D₊(ω) + D₋(ω) so that d₊(τ) and d₋(τ) vanish for times τ < 0 and τ > 0, respectively. For the averaged cross spectrum, the matrix Q is defined in (5), the non-averaged cross spectrum can be recovered as a special case setting Q = M = W, because the convolution equations (3) and (5) have the same structure and symmetries. If all eigenvalues of MH(ω) have an absolute value smaller than unity, the propagator P(ω) can be expanded into a geometric series in which the nth term contains only interactions via n steps in the connectivity graph. This expansion has been used to obtain the contribution of different motifs to the integral of covariance functions [40] and their temporal shape [41].

In the following, we neglect the continuous part of the autocovariance function, setting a(t) = rδ(t), because the δ-peak is typically dominant. With this replacement we mainly neglect the trough around 0 due to the relative refractoriness. An estimate of the error can be obtained considering the respective weights of the delta peak (r) and the trough. From the relation $\int _{-\infty }^{\infty }a(t)\, \mathrm {d}t=r\mathrm {CV^{2}}$ [42, 43] the integral weight of the trough follows as r(CV² − 1), which is small for irregular spike trains with a coefficient of variation CV close to unity.

For an arbitrary causal kernel h it follows that D₊(ω) = r H(ω), so the cross spectrum takes the form

$$\begin{eqnarray} \fl \mathbf{C}(\omega) &=& r\frac{Kw}{N}\left(\begin{array}{@{}cc@{}} 1 & -g\\ 1 & -g \end{array}\right)U(\mathrm{i}\,\omega)+\textrm{c.c. trans. } + r(1+g^{2}\gamma)\frac{(Kw)^{2}}{N}\left(\begin{array}{@{}cc@{}} 1 & 1\\ 1 & 1 \end{array}\right)|U(\mathrm{i}\,\omega)|^{2}\nonumber \\ \textrm{with } U(z) &=&\frac{1}{H^{-1}(-\mathrm{i}\,z)-L}\quad {\rm and}\quad L =Kw(1-\gamma g).\nonumber \end{eqnarray}\noindent \tag{ 6 }$$

The limit ω → 0 corresponds to the time integral of the cross covariance function, approximating the count covariance for long time bins [2, 8]. With A(0) = r, the integral correlation coefficient averaged over neuron pairs fulfils the equation

$$\begin{eqnarray} \fl \frac{\mathbf{C}(0)}{A(0)} = \frac{Kw}{N}\frac{1}{1-L}\left(\begin{array}{@{}cc@{}} 2 & 1-g\\ 1-g & -2g \end{array}\right) +\frac{(Kw)^{2}}{N}\frac{1+g^{2}\gamma}{(1-L)^{2}}\left(\begin{array}{@{}cc@{}} 1 & 1\\ 1 & 1 \end{array}\right), \end{eqnarray}\noindent \tag{ 7 }$$

which has previously been derived from noise-driven linear rate dynamics [12]. The quantity L plays a key role here: it determines the feedback magnitude of in-phase fluctuations of the excitatory and inhibitory populations. The stability of the average firing rate requires this feedback to be sufficiently small [9], i.e. L < 1, indicated by the pole in equation (7) at L = 1. Typically cortical networks are in the balanced regime [9, 13, 14], i.e. L < 0. For such inhibition-dominated networks, the denominator in equation (7) is larger than unity, indicating a suppression of covariances [12]. As shown in figure 3(a), the prediction (7) agrees well with the results of simulations of leaky integrate-and-fire networks for a broad range of network sizes N.

Zoom In Zoom Out Reset image size

Previous works have investigated neural networks in the thermodynamic limit N → ∞ [11, 13], scaling the synaptic amplitudes $J\propto 1/\sqrt {N}$ in order to arrive at analytical results. Such a scaling increases the feedback on the network level $L\propto \sqrt {N}$ and therefore changes the collective network state. Equation (7) provides an alternative criterion to scale the synapses while keeping the dynamics comparable: The two terms in equation (7) depend differently on the feedback L. In order to maintain their ratio, we need to keep the population feedback L constant. The synaptic amplitude J (approximately ∝w (A.3)) hence needs to scale as J∝1/N. In addition, the response kernel h of each single neuron must remain unchanged, requiring the same working point, characterized by the mean μ_i and fluctuations σ_i in the input into each cell. A constant mean directly follows from L = const, but the variance due to local input from other neurons in the network decreases as 1/N. To compensate, we supply each neuron with an additional external uncorrelated balanced noise whose variance appropriately increases with N (as described in detail in appendix E). Figure 3(b) shows that the shape of the covariance functions is invariant over a large range of network sizes N; in particular the apparent time lag of inhibition behind excitation observed as the asymmetry in figure 3(b) does not vanish in the limit N → ∞. The magnitude of the covariance decreases as 1/N as expected from equation (7), because Kw = const.

Global properties of the network dynamics can be inferred by considering the spectrum of equation (6), those complex frequencies z_k at which the expression has a pole due to the function U(z). These poles are resonant modes of the network, where the real part ℜ(z_k) denotes the damping of the mode and the imaginary part ℑ(z_k) is the oscillation frequency. A pole appears whenever z_k is a single root of U⁻¹(z_k) = H⁻¹(−iz_k) − L = 0. With the Fourier representation $H(\omega )=\frac {\mathrm {e}^{-\mathrm {i}\omega d}}{1+\mathrm {i}\omega \tau _{\mathrm {e}}}$ of the response kernel (4), the poles correspond to the spectrum of the delay differential equation $\tau _{\mathrm {e}}\frac {\mathrm {d}y}{\mathrm {d}t}(t)=-y(t)+Ly(t-d)$ (cf [44]) which describes the evolution of the population-averaged activity. As shown in appendix C, the location of the poles can be expressed by the branches k of the Lambert W function, the solution of W_ke^W_k = x [45], as

$$\begin{eqnarray} z_{k} & = & -\frac{1}{\tau_{\mathrm{e}}}+\frac{1}{d}W_{k}\left(L\frac{d}{\tau_{\mathrm{e}}}\,\mathrm{e}^{\frac{d}{\tau_{\mathrm{e}}}}\right),\quad k\in\mathbb{N}_{0}. \end{eqnarray} \tag{ 8 }$$

The spectrum only depends on the population feedback L, the delay d and the effective time constant τ_e of the neural response kernel. This explains why keeping L constant while scaling the network in figure 3 yields shape-invariant covariance functions. A typical spectrum is shown in figure 4(a) as an inset, where each dot marks one of the poles z_k (8). The two principal branches of W_k are the modes with the largest real part ℜ(z_k), and hence with the least damping, dominating the network dynamics. The remaining branches appear as conjugate pairs and their real parts are more negative, corresponding to stronger damping. Investigating the location of the principal branches therefore enables us to classify the dynamics in the network. Their dependence on the delay is shown in figure 4(a) as a parametric plot in d. The point at which the two real principal solutions turn into a complex conjugate pair marks a transition from purely exponentially decaying dynamics to damped oscillations. This happens at a sufficiently strong negative coupling L or a sufficiently long delay d, precisely when the argument of W_k is smaller than −e⁻¹ [45], leading to the condition

$$\begin{eqnarray} L & < & -\frac{\tau_{\mathrm{e}}}{d}\,\mathrm{e}^{-\frac{d}{\tau_{\mathrm{e}}}-1}. \end{eqnarray} \tag{ 9 }$$

The gray cross marks this point in figure 4(a); the gray curve shows the corresponding relation of feedback and delay in the phase diagram of figure 4(b). In the region below the curve, the dominant mode of fluctuations in the network is thus damped oscillatory, whereas above the curve fluctuations are relaxing exponentially in time.

Zoom In Zoom Out Reset image size

For a sufficiently long delay d the principal poles may assume positive real values, leading to ongoing oscillations—a Hopf bifurcation. The condition under which this happens can be derived from H⁻¹(ω_crit) = L, as detailed in appendix C. Equating the absolute values on both sides leads to the condition $\omega _{\mathrm {crit}}\tau _{\mathrm {e}}=\sqrt {L^{2}-1}$: oscillations can only be sustained if the negative population feedback is sufficiently strong, L < −1. The oscillation frequency increases the stronger the negative feedback. The condition for the phases leads to the critical delay required for the onset of oscillations (see appendix C for details)

$$\begin{eqnarray} \frac{d_{\mathrm{crit}}}{\tau_{\mathrm{e}}} & = & \frac{\pi-\arctan(\sqrt{L^{2}-1})}{\sqrt{L^{2}-1}}. \end{eqnarray}\noindent \tag{ 10 }$$

This relation is shown as the black curve in the phase diagram of figure 4(b). The oscillatory frequency on the bifurcation line, at the onset of oscillations, can be expressed as

$$\begin{eqnarray} 2\pi f_{\mathrm{crit}}d=\omega_{\mathrm{crit}}d & = & \pi-\arctan(\sqrt{L^{2}-1}), \end{eqnarray}\noindent \tag{ 11 }$$

which is shown in the left sub-panel of the phase diagram of figure 4(b). Consequently, the oscillation frequency f_crit at the onset is between (4d_crit)⁻¹ and (2d_crit)⁻¹, depending on the strength of the feedback L, approaching f_crit = (4d_crit)⁻¹ at the onset with increasing negative feedback.

Changing the synaptic delay homogeneously for all synapses in the network allows us to observe the transition of the network from exponentially damped to oscillatory damped and finally to oscillatory dynamics. For a short delay of d = 0.5 ms the dynamics is dominated by the single real pole near −1/τ_e (brown dot in figure 4(a)) and the covariance function is exponentially decaying (figure 4(c)). Increasing the delay to d = 1 ms the principal poles split into a complex conjugate pair as the delay crosses the gray curve in figure 4(b) so that side troughs become visible in the covariance function in figure 4(c). Further increasing the delay, the network approaches the point of oscillatory instability, where a Hopf bifurcation occurs, marked by black crosses in figure 4(a) and the black curve in figure 4(b). The damping of oscillations decreases as the system approaches the bifurcation (figure 4(c)). The structure of the autocovariance function of single spike trains (figure 4(d)) is dominated by the dip due to refractoriness of the neuron after reset. At pronounced network oscillations, the autocorrelation function shows a corresponding modulation. The neglect of the oscillating continuous part of the autocorrelation function in the theory does apparently not have a pronounced effect on the obtained cross correlation functions, evidenced by the good agreement between direct simulation and theory for d = 5 ms. For weaker oscillations, e.g. at d = 3 ms, the coherence time of the oscillations is typically shorter than the width of the dip in the autocorrelation. The phase coherence of the oscillation is then shorter than the typical inter-spike-interval and single neuron spike trains are irregular. This state is known as synchronous irregular activity [9], where the activity of a single neuron is irregular, but collectively the neurons participate in a global oscillation.

If the network is not in the oscillatory state, all modes are damped in time, i.e. all poles z_k of the function U(z) appearing in equation (6) lie in the left complex half-plane, ℜ(z_k) < 0. Hence, the dynamics is stable and we can expect to obtain a unique solution for the covariance as an observable of this stable dynamics. We perform the Fourier transform to the time domain using the residue theorem $u(t)=\frac {1}{2\pi \mathrm {i}}\oint _{E_{t}}U(z)\, \mathrm {e}^{zt}\, \mathrm {d}z=\sum _{z_{k}\in E_{t}}\mathrm {Res}(U,z_{k})\, \mathrm {e}^{z_{k}t}$, where the integration path E_t proceeds along the imaginary axis from −i∞ to i∞ and is then closed in infinity in the left half-plane (for t > d) or the right half-plane (for 0 < t < d) to ensure convergence, resulting in (see appendix D for the detailed calculation)

$$\begin{eqnarray} u(t) & = & \sum_{k=0}^{\infty}\frac{1}{\left(1+z_{k}\tau_{\mathrm{e}}\right)d+\tau_{\mathrm{e}}}\Theta(t-d)\,\mathrm{e}^{z_{k}(t-d)}. \end{eqnarray} \tag{ 12 }$$

The back transform of $V(\omega){\stackrel{\mathrm{def}}{\,=\,}}|U(\omega)|^{2}$ proceeds along similar lines and results in

$$\begin{eqnarray} v(t) & = & \sum_{k=0}^{\infty}\frac{1}{\left(1+z_{k}\tau_{\mathrm{e}}\right)d+\tau_{\mathrm{e}}}\,\frac{\mathrm{e}^{z_{k}|t|}}{(1-z_{k}\tau_{\mathrm{e}})-L\,\mathrm{e}^{z_{k}d}}. \end{eqnarray} \tag{ 13 }$$

The population-averaged covariance functions in the time domain then follow from equation (6) for t > 0 as

$$\begin{eqnarray} \mathbf{c}(t) & = & r\,\frac{Kw}{N}\left(\begin{array}{@{}cc@{}} 1 & -g\\ 1 & -g \end{array}\right)\; u(t) + r\,\frac{(Kw)^{2}}{N}(1+g^{2}\gamma)\left(\begin{array}{@{}cc@{}} 1 & 1\\ 1 & 1 \end{array}\right)\; v(t), \end{eqnarray} \tag{ 14 }$$

which is the central result of our work. Figure 5 shows the comparison of the theory (14) with the covariance functions obtained by direct simulation. The analytical expression reveals that the covariance functions are composed of two components: the first term in equation (14) has heterogeneous matrix elements and hence depends on the neuron types under consideration. Its origin is illustrated in figure 5(a): if one of the neurons emits a spike, as indicated, this impulse travels along its axon and reaches the target neurons after one synaptic delay d. Depending on the type of the source neuron, the impulse excites (synaptic amplitude J) or inhibits (−gJ) its targets. Its effect is therefore visible in the pairwise covariance function as a positive (blue) or negative (red) deflection, respectively, in figure 5(c). This deflection not only contains the direct synaptic effect, but also infinitely many reverberations of the network, seen formally in equation (12). This expression is not proportional to the kernel h(t) directly, but is rather a series including the whole spectrum of the network. The shape of the spike echo consequently shows onsets of reverberations at integer multiples of the synaptic delay (figure 5(c)), being transmitted over multiple synaptic connections. The contribution of the second term in equation (14) follows the intuitive argument illustrated in figure 5(b). The incoming activity from the network to a pair of neurons is correlated. As the input statistics is the same for each neuron, this contribution is identical for any pair of neurons (green curve in figure 5(c)). The sum of both components results in the covariance functions shown in figures 5(d)–(f). The same analytical solution is shown for different delays in figure 4(c) showing good agreement with direct simulation. For different sizes of simulated networks in figures 3(b)–(d) the analytical expression (14) explains why the spike echo does not become negligible in the thermodynamic limit N → ∞: for fixed population feedback L, both contributions in equation (14) scale as 1/N, so the relative contribution of the echo stays the same. This also explains the apparent paradox (see figure 1), that covariance functions in recurrent networks not only depend on the input statistics, but in addition the spike feedback causes a reverberating echo. The power spectrum of population averages is dominated by pairwise covariances, explaining the different spectra observed in the excitatory and inhibitory population activity [15]. Scaling the network so as to keep the marginal statistics of single neurons constant, $J\propto w\propto 1/\sqrt {N}$ [11, 13] changes the spectrum (8), because the feedback increases as $L\propto \sqrt {N}$ which can ultimately lead to oscillations as shown in figure 4(c).

Zoom In Zoom Out Reset image size

The present work qualitatively explains certain features of the correlation structure of simultaneously recorded synaptic currents of two cells in vivo. Novel experimental techniques are able to separate contributions of excitatory and inhibitory inputs [29]. We calculate such covariances in a random network and show that the covariance between synaptic impulses decomposes into a linear combination of the covariance of the spiking activity and the autocovariance functions (see the caption of figure 6). Each synaptic impulse has a certain time course, here modeled as a first-order low-pass filter with time constant τ_s = 2 ms (see (1)). The covariances between these filtered currents are shown in figure 6. Their temporal structure resembles those measured in the cortex in vivo [29, their figures 1(e) and (f)]: covariances between afferents of the same type are monophasic and positive, while the covariances between excitatory and inhibitory afferents are biphasic and mostly negative. The lag reported between inhibitory and excitatory activity [29, their figure 2(b)], which was also observed in binary random networks [11, 25], is explained by the echo of the spike contributing to the covariance function. In contrast to previous work, we take the delayed and pulsed synaptic interaction into account. Without delays and with binary neurons [11, 25] the echo appears as a time lag of inhibition with respect to excitation.

Zoom In Zoom Out Reset image size

Measurements of membrane potential fluctuations during quiet wakefulness in the barrel cortex of mice [46] showed that correlations between inhibitory neurons are typically narrower than those between two excitatory neurons [46, their figures 4(A), 5(B) and 5(C) and (E)]. These results qualitatively agree with our theory for covariances between the spiking activity, because fluctuations of the membrane potential are uniformly transferred to fluctuations of the instantaneous firing intensity. The direct measures of spiking activity [46, their figure 6] confirm the asymmetric correlation between excitation and inhibition. The low correlation between excitatory neurons reported in that study may partly be due to the unnormalized, firing rate-dependent measure and the low rates of excitatory neurons.

The qualitative features of the cross correlation functions, namely their different widths and their asymmetry for excitatory–inhibitory pairs, are generic and agree with experimental results. They are fully explained by the decomposition into an echo term and a term corresponding to the feed-forward transmission of correlation. This decomposition merely relies on the fact that the autocorrelation of a spike train has a delta peak and that a spike triggers a response in the target cell with positive or negative sign for excitatory and inhibitory connections, respectively. Hence, we expect that these features are robust and survive also for more realistic network models with many heterogeneous subpopulations. For weakly correlated fluctuations, if the input to each cell is sufficiently noisy, a linear approximation of the neuronal response provides a viable first order approximation also for nonlinear neuron dynamics, as shown here. We suspect that a deeper reason why such a linearization is possible is the inherent decorrelation [12] by negative feedback in networks in the inhibition dominated regime. The decorrelation keeps population fluctuations small and hence prevents strong excursions that would exceed the validity of the linear approximation.

Oscillations in the γ range (25–100 Hz) are ubiquitous in population measures of neural activity in humans, and have earlier been explained in networks of leaky integrate-and-fire model neurons [9] by the Hopf bifurcation induced by delayed negative feedback. For the regime of high noise we here uncover a simpler analytical condition for the onset (10) and frequency (11) of fast global oscillations. For lower noise, deviations of the nonlinear leaky integrate-and-fire dynamics from the linear theory presented here are expected.

A traditional motivation to scale the network size to infinity is to obtain an analytic solution for an otherwise difficult problem. Biologically realistic networks have a finite size. In the present work we have determined the correlation structure for such finite-sized networks analytically. We present the scaling of the correlation functions in figure 3 to relate our work to previous results that applied scaling arguments to obtain the correlation structure [11]. The latter work investigated networks of binary neurons without conduction delays and assumed a scaling of the synaptic amplitudes $J\propto 1/\sqrt {N}$. Such a scaling increases the overall feedback strength $L\propto pNw\propto \sqrt {N}$. This has two consequences: firstly, as seen from (7), the relative contributions of the spike echo and the feed forward term change with L and hence with network size. Therefore the shape of correlation functions depends on the network size. Secondly, the overall dynamic regime of the network is affected by the scaling. This can be seen from figure 4(b). For non-zero synaptic conduction delays, the network eventually becomes oscillatory if the network size exceeds a certain value. This happens precisely at the point where L crosses the bifurcation line shown in figure 4(b). We therefore propose an alternative scaling J∝1/N here for which we show that it preserves the shape of correlation functions and the overall network state. Only the magnitude of the cross correlation functions decreases ∝1/N. However, a caveat of this scaling is that while it preserves the global network properties, it affects the working point of each individual neuron, because the fluctuations due to the local synaptic input decrease $\propto 1/\sqrt {N}$. We alleviated this shortcoming by supplying each neuron with additional, uncorrelated noise.

Our results are based on a simplified network model composed of two homogeneous (excitatory and inhibitory) subpopulations of leaky integrate-and-fire neurons. The real cortex, in contrast, is highly heterogeneous in several respects: Its layered structure with layer-specific neuron and connection properties (e.g. time constants, spike thresholds, synaptic weights, in-degrees) requires the distinction of more than two subpopulations. Even within each subpopulation and for each combination of subpopulations, the neuron and connection parameters are not constant but broadly distributed. Further, a plethora of cortical neurons, in particular various types of interneurons, exhibit a much richer dynamical repertoire (e.g. resonating behavior) than the leaky integrate-and-fire neuron model (regular spiking integrator). In principle, the mathematical framework presented in this paper can be extended to networks composed of n (n > 2) heterogeneous subpopulations of different neuron types. This would require to account for subpopulation-specific working points (e.g. due to layer-specific firing rates; see [47]) and for the effect of parameter distributions [48, 49] and the single-neuron dynamics on the effective linearized subpopulation responses. This extended theory would result in an n-dimensional linear algebraic equation for the subpopulation-averaged cross spectra, similar to (6) where n = 2. For n > 2, this equation most likely needs to be solved numerically.

A fundamental assumption of the presented theory is that the network states show irregular single neuron dynamics. This requirement arises from the analytical description replacing spike trains by spike densities and a stochastic realization of spikes. Regular spike trains are outside the scope of such a description. Moreover, the approximation of the neuronal response to linear order is only a viable approach in sufficiently asynchronous network states with low correlations. States with stronger correlations, such as observed in convergent–divergent feed-forward structures [34, 50], require an explicit treatment of the nonlinear response [36].

From a physics viewpoint, neuronal networks unite several interesting properties. They do not reach thermodynamic equilibrium even in the stationary state, as detailed balance does not hold for all pairs of states of the system. In detailed balance, the rate of transition from one state to another is equal to the rate of the reverse transition. For a leaky integrate-and-fire neuron the state of the neuron is uniquely determined by its membrane voltage. In the stationary state neurons fire with a constant rate, so there is a continuous flux of the neurons' voltage from reset up to threshold. Imagining a discretization of the voltage axis we see that a pair of adjacent voltage-intervals between reset and threshold is more often traversed from lower to higher voltage than in the reverse direction, so obviously detailed balance does not hold even for a single neuron in the stationary state. Moreover, the interaction between pairs of neurons is directed, delayed, pulsed and depends on the flux of the sending neuron's state variable at threshold. In contrast, pairwise interactions frequently studied in physics, like the Coulomb interaction or exchange interaction, can be expressed by a pair potential and are thus symmetric (undirected), instantaneous and depend directly on the state variables (e.g. spatial coordinates or spins) of the pair of interacting particles. Non-equilibrium systems are at the heart of ubiquitous transport phenomena, like heat or electric conduction. Understanding fluctuations in such a system marks the starting point to infer macroscopic properties by the assertion of the fluctuation–dissipation theorem that connects microscopic fluctuations to macroscopic transport properties. Despite the non-equilibrium dynamics and the non-conservative pairwise interaction, in this paper we develop a simple analytical framework merely based on linear perturbation theory that explains the time-dependent covariance functions of the activity of pairs of integrate-and-fire model neurons in a recurrent random network. Formally our approach resembles the step from the kinetic Ising model near equilibrium to its non-equilibrium counterpart, the network of binary neurons [25]. A difference is the spiking interaction considered in our work, which led us to describe each neuron in terms of the flux over threshold (spike train) rather than by its state variables (membrane voltage and synaptic current). In this respect, we follow the established mean-field approach for spiking neuronal systems [9, 14]. However, while this mean field approach proceeds by assuming vanishing correlations to obtain the dynamics of the ensemble-averaged activity, we here derive and solve the self-consistency equation for the pairwise-averaged covariances of the microscopic system.

The typical time scale of covariance functions found here coincides with the time window of biological synaptic plasticity rules [30], so that non-trivial interactions of dynamics and structure are expected. It is our hope that the novel capability to resolve the temporal structure of covariances in spiking networks presented here proves useful as a formal framework to further advance the theory of these correlated non-equilibrium systems and in particular serves as a further stepping stone in the endeavor to understand how learning on the system level is implemented by the interplay of neuronal dynamics and synaptic plasticity.

We are grateful to Birgit Kriener, Sonja Grün and George Gerstein for discussions and exploratory simulations in fall 2005 that led to the observation of asymmetric covariance functions. This work is partially supported by the Helmholtz Association: HASB and portfolio theme SMHB, the Next-Generation Supercomputer Project of MEXT, EU grant 15879 (FACETS), EU grant 269921 (BrainScaleS). All network simulations were carried out with NEST (http://www.nest-initiative.org).

The response kernel h_ij needs to be related to the dynamics of the neuron model (1). Here we present an approximation of this kernel which is sufficiently accurate to allow quantitative predictions, but yet simple enough to enable an analytical solution for the correlation structure. If the synaptic time constant is short, τ_s ≪ τ_m, the synaptic amplitude J can be thought of as the amplitude of the jump in the membrane potential V caused upon arrival of an incoming impulse. If correlations between incoming spike trains are sufficiently small, the first and second moments of the summed impulses $\tau _{\mathrm {m}}\sum _{j}J_{ij}s_{j}(t-d)$ are $\mu _{i}=\tau _{\mathrm {m}}\sum _{j}J_{ij}r_{j}$ and $\sigma _{i}^{2}=\tau _{\mathrm {m}}\sum _{j}J_{ij}^{2}r_{j}$, respectively, if the inputs' statistics can be approximated by Poisson processes of rate r_j each. For small J and a high total rate, the system of differential equations (1) is hence approximated by a stochastic differential equation driven by a unit variance Gaussian white noise ξ

$$\begin{eqnarray*} \tau_{\mathrm{m}}\frac{\mathrm{d}V_{i}}{\mathrm{d}t} & = & -V_{i}+I_{s,i}(t),\\ \tau_{\mathrm{s}}\frac{\mathrm{d}I_{i}}{\mathrm{d}t} & = & -I_{i}+\mu_{i}+\sigma_{i}\sqrt{\tau_{\mathrm{m}}}\xi(t). \end{eqnarray*}$$

The stationary firing rate in this limit is given by [33]

$$\begin{eqnarray} &&\fl r_{i}^{-1} = \tau_{\mathrm{r}}+\tau_{\mathrm{m}}\sqrt{\pi}\left(F(y_{\theta})-F(y_{r})\right),\nonumber\\ &&f(y) = \mathrm{e}^{y^{2}}(1+\mathrm{erf}(y))\quad F(y)=\int^{y}f(y)\, \mathrm{d}y\\ &&\textrm{with }\,y_{\theta,r} = \frac{V_{\theta,r}-\mu_{i}}{\sigma_{i}}+\frac{\alpha}{2}\sqrt{\frac{\tau_{\mathrm{s}}}{\tau_{\mathrm{m}}}}\quad\alpha=\sqrt{2}\left|\zeta\left(\frac{1}{2}\right)\right|, \nonumber \end{eqnarray} \tag{ A.1 }$$

with Riemann's zeta function ζ. The rate r_i is the density of action potentials over time. The response of the firing density of the neuron i at time point t with respect to a point-like deflection of the afferent input s_j at time point t' defines the response kernel as the functional derivative

$$\begin{eqnarray*} &&\fl \left\langle \frac{\delta s_{i}(t)}{\delta s_{j}(t^{\prime})}\right\rangle _{\mathbf{s}{\backslash} s_{j}} = \lim_{\epsilon\rightarrow0}\frac{1}{\epsilon}\langle s_{i}(t,\{\mathbf{s}(\tau)+\epsilon\delta(\tau-t^{\prime})\mathbf{e}_{j}|\tau<t\}) -s_{i}(t,\{\mathbf{s}(\tau)|\tau<t\})\rangle _{\mathbf{s}{\backslash} s_{j}}\\ &&\fl \hphantom{\left\langle \frac{\delta s_{i}(t)}{\delta s_{j}(t^{\prime})}\right\rangle _{\mathbf{s}{\backslash} s_{j}}} \stackrel{\mathrm{def}}{=} h_{ij}(t-t^{\prime})=w_{ij}\; h(t-t^{\prime}). \end{eqnarray*}$$

Here we used the homogeneity, namely the identical input statistics of each neuron i, leading to the same temporal shape h(t) independent of i and the stationarity, so that the kernel only depends on the time difference t − t'. We choose h(t) to have unit integral and define w_ij as the integral of the kernel. We determine the temporal integral of the kernel as

$$\begin{eqnarray} w_{ij} & = & \int_{-\infty}^{\infty}h_{ij}(t)\; \mathrm{d}t\nonumber\\ & = & \frac{\partial r_{i}}{\partial r_{j}}. \end{eqnarray} \tag{ A.2 }$$

The second equality holds because the integral of the impulse response equals the step response [51]. Further, a step in the density s_j corresponds to a step of r_j. Up to linear approximation the effect of the step in the rate r_j on the rate r_i can be expressed by the derivative considering the perturbation of the mean μ_i and the variance σ²_i upon change of r_j. Using equation (A.1) we note the chain rule $\frac {\partial r_{i}}{\partial r_{j}}=-r_{i}^{2}\frac {\partial r_{i}^{-1}}{\partial r_{j}}$. The latter derivative follows as $\frac {\partial r_{i}^{-1}}{\partial r_{j}}=\frac {\partial r_{i}^{-1}}{\partial y_{\theta }}\frac {\partial y_{\theta }}{\partial r_{j}}+\frac {\partial r_{i}^{-1}}{\partial y_{r}}\frac {\partial y_{r}}{\partial r_{j}}$. The first derivative in both terms yields $\frac {\partial r_{i}^{-1}}{\partial y_{A}}=\sqrt {\pi }\tau _{\mathrm {m}} f(y_{A})$ with y_A∈{y_θ,y_r}. The second derivative evaluates with $y_{A}=\frac {A-\mu _{i}}{\sigma _{i}}+\frac {\alpha }{2}\sqrt {\frac {\tau _{\mathrm {s}}}{\tau _{\mathrm {m}}}}$ to $\frac {\partial y_{A}}{\partial r_{i}}=-\frac {1}{\sigma _{i}}\tau _{\mathrm {m}} J_{ij}-\frac {A-\mu _{i}}{\sigma _{i}^{2}}\frac {\tau _{\mathrm {m}} J_{ij}^{2}}{2\sigma _{i}}=-\tau _{\mathrm {m}}\frac {J_{ij}}{\sigma _{i}}(1+\frac {A-\mu _{i}}{\sigma _{i}}\frac {J_{ij}}{2\sigma _{i}})$. So together we obtain

$$\begin{eqnarray} &&\fl w_{ij} = \alpha J_{ij}+\beta J_{ij}^{2},\nonumber\\ &&\textrm{with }\alpha = \sqrt{\pi}(\tau_{\mathrm{m}} r_{i})^{2}\frac{1}{\sigma_{i}}\left(f(y_{\theta})-f(y_{r})\right),\\ &&\textrm{and }\beta = \sqrt{\pi}(\tau_{\mathrm{m}} r_{i})^{2}\frac{1}{\sigma_{i}}\left(f(y_{\theta})\, J_{ij}\frac{V_{\theta}-\mu_{i}}{2\sigma_{i}^{2}}-f(y_{r})\, J_{ij}\frac{V_r-\mu_{i}}{2\sigma_{i}^{2}}\right).\nonumber \end{eqnarray} \tag{ A.3 }$$

The autocovariance A(ω) in the frequency domain ($F(\omega )=\mathcal {F}[f](\omega )=\int _{-\infty }^{\infty }f(t)\,\mathrm {e}^{-\mathrm {i}\,\omega t}\, \mathrm {d}t$) has two different terms. The first term is a constant r due to the spiking with rate r resulting from the delta peak rδ(t) in the time domain. The second term is the continuous function a_c(t), for example due to the refractoriness of the neuron. It further follows from a(t) = a(−t) that A(ω) = A(−ω). For τ > 0 the covariance matrix fulfils the linear convolution equation (5). As this equation only holds for the positive half of the time axis, we cannot just apply the Fourier transform to obtain the solution. For negative time lags τ < 0 the covariance matrix is determined by the symmetry c(τ) = c^T(−τ). Here we closely follow [26] and employ the Wiener–Hopf theory [39] to derive an equation for the cross spectral matrix in the frequency domain that has the desired symmetry and solves (5) simultaneously. To this end we introduce the auxiliary matrix $\mathbf {b}(\tau )=\left (h\ast (\mathbf {M}\mathbf {c})\right )(\tau )+\mathbf {Q}(h\ast a)(\tau )-\mathbf {c}(\tau )$ for −∞ < τ < ∞. Obviously, b(τ) = 0 for τ > 0. Since the defining equation for b holds on the whole time axis, we may apply the Fourier transform to obtain $\mathbf {B}(\omega )=H(\omega )\left (\mathbf {M}\mathbf {C}(\omega )+\mathbf {Q} A(\omega )\right )-\mathbf {C}(\omega )$. Solving for C

$$\begin{eqnarray} \mathbf{C}(\omega) & = & (\mathbf{1}-\mathbf{M} H(\omega))^{-1}\left(H(\omega)A(\omega)\mathbf{Q}-\mathbf{B}(\omega)\right) \end{eqnarray} \tag{ B.1 }$$

and using the symmetry C(ω) = C^T(−ω) we obtain the equation

$$\begin{eqnarray*} &&\fl \left(H(\omega)A(\omega)\mathbf{Q}-\mathbf{B}(\omega)\right)\left(\mathbf{1}-\mathbf{M}^{\mathrm{T}}H(-\omega)\right) = \left(\mathbf{1}-\mathbf{M} H(\omega)\right)\left(H(-\omega)A(-\omega)\mathbf{Q}^{\mathrm{T}}-\mathbf{B}^{\mathrm{T}}(-\omega)\right). \end{eqnarray*}\noindent$$

We observe that QM^T = MQ^T is symmetric and with A(ω) = A(−ω) the term proportional to |H(ω)|² cancels on both sides, resulting in

$$\begin{eqnarray} &&\fl H(\omega)A(\omega)\mathbf{Q}+\left(\mathbf{1}-\mathbf{M} H(\omega)\right)\mathbf{B}^{\mathrm{T}}(-\omega) = H(-\omega)A(\omega)\mathbf{Q}^{\mathrm{T}}+\mathbf{B}(\omega)\left(\mathbf{1}-\mathbf{M}^{\mathrm{T}}H(-\omega)\right). \end{eqnarray} \tag{ B.2 }$$

We next introduce D(ω) = H(ω)A(ω) which we split into D(ω) = D₊(ω) + D₋(ω), chosen such that d₊(t) (in the time domain) vanishes for t < 0 and d₋(t) vanishes for t > 0. Consequently the Fourier transforms of both terms may have poles in distinct complex half-planes: D₊(ω) may only have poles in the upper half-plane ℑ(ω) > 0 and the function vanishes for lim_{|ω|→∞,ℑ(ω)<0}D₊(ω) = 0, following from the definition of the Fourier integral. For D₋(ω) the half-planes are reversed. The analytical properties of H(ω) are thus similar to those of D₊(ω); those of B(ω) are similar to D₋(ω). We sort the terms in (B.2) such that the left hand side only contains terms that vanish at infinity in the lower half-plane ℑ(ω) < 0, the right hand side those that vanish at infinity in the upper half-plane ℑ(ω) > 0

$$\begin{eqnarray} &&\fl D_{+}(\omega)\mathbf{Q}-D_{-}(-\omega)\mathbf{Q}^{\mathrm{T}}+\left(\mathbf{1}-\mathbf{M} H(\omega)\right)\mathbf{B}^{\mathrm{T}}(-\omega)\nonumber\\ &&= D_{+}(-\omega)\mathbf{Q}^{\mathrm{T}}-D_{-}(\omega)\mathbf{Q} +\mathbf{B}(\omega)\left(\mathbf{1}-\mathbf{M}^{\mathrm{T}}H(-\omega)\right). \end{eqnarray} \tag{ B.3 }$$

The left hand side consequently is analytic for ℑ(ω) < 0, the right hand side is analytic for ℑ(ω) > 0, so (B.3) defines a function that is analytic on the whole complex plane and that vanishes at the border for |ω| → ∞. Hence by Liouville's theorem it is 0 and we can solve the right hand side of (B.3) for B

$$\begin{eqnarray*} \mathbf{B}(\omega) & = & \left(D_{-}(\omega)\mathbf{Q}-D_{+}(-\omega)\mathbf{Q}^{\mathrm{T}}\right)\left(\mathbf{1}-\mathbf{M}^{\mathrm{T}}H(-\omega)\right)^{-1}. \end{eqnarray*}$$

Inserted into (B.1) this yields the definition $\mathbf {P}(\omega )=\left (\mathbf {1}-\mathbf {M} H(\omega )\right )^{-1}$

$$\begin{eqnarray*} &&\fl \mathbf{C}(\omega) = \mathbf{P}(\omega)\left(H(\omega)A(\omega)\mathbf{Q}-\left(D_{-}(\omega)\mathbf{Q}-D_{+}(-\omega)\mathbf{Q}^{\mathrm{T}}\right)\mathbf{P}^{\mathrm{T}}(-\omega)\right)\\ &&\fl \hphantom{\mathbf{C}(\omega)} = \mathbf{P}(\omega)\left(\left(D_{+}(\omega)+D_{-}(\omega)\right)\mathbf{Q}(\mathbf{1}{-}\mathbf{M}^{\mathrm{T}}H({-}\omega)) -D_{-}(\omega)\mathbf{Q}+D_{+}({-}\omega)\mathbf{Q}^{\mathrm{T}}\right)\mathbf{P}^{\mathrm{T}}({-}\omega)\\ &&\fl \hphantom{\mathbf{C}(\omega)} = \mathbf{P}(\omega)\left(D_{+}(\omega)\mathbf{Q}+D_{+}(-\omega)\mathbf{Q}^{\mathrm{T}}-A(\omega)|H(\omega)|^{2}\mathbf{Q}\mathbf{M}^{\mathrm{T}}\right)\mathbf{P}^{\mathrm{T}}(-\omega). \end{eqnarray*}$$

The latter expression can be brought to the form

$$\begin{eqnarray} &&\fl \mathbf{C}(\omega) = D_{+}(\omega)\mathbf{P}(\omega)\mathbf{Q}+D_{+}({-}\omega)\mathbf{Q}^{\mathrm{T}}\mathbf{P}^{\mathrm{T}}({-}\omega) {+}\left(D_{+}(\omega)H(-\omega)+D_{+}({-}\omega)H(\omega){-}A(\omega)|H(\omega)|^{2}\right)\nonumber\\ &&\times\mathbf{P}(\omega)\mathbf{M}\mathbf{Q}^{\mathrm{T}}\mathbf{P}^{\mathrm{T}}(-\omega), \end{eqnarray} \tag{ B.4 }$$

which has the advantage that the first two terms have poles in distinct half-planes ℑ(ω) > 0 and ℑ(ω) < 0, respectively. This means these terms contribute for either positive and negative times, respectively, and the last term contributes for positive and negative times.

With the Fourier representation $H(\omega )=\frac {\mathrm {e}^{-\mathrm {i}\,\omega d}}{1+\mathrm {i}\,\omega \tau _{\mathrm {e}}}$ of the delayed exponential kernel (4) the averaged cross spectrum (6) contains the two functions U(z) and V (z) = |U(z)|² defined on the complex frequency plane z = i ω. These functions may exhibit poles. The function U has a pole z_k whenever the denominator has a single root H⁻¹(−iz) − L = 0, which amounts to the condition (1 + z_kτ_e) e^z_kd = L. These complex frequencies can be expressed by the Lambert W function, the solution of We^W = x [45], by

$$\begin{eqnarray} &&(1+z_{k}\tau_{\mathrm{e}})\,\mathrm{e}^{z_{k}d} = L,\nonumber\\ &&\left(\frac{d}{\tau_{\mathrm{e}}}+z_{k}d\right)\,\mathrm{e}^{\frac{d}{\tau_{\mathrm{e}}}+z_{k}d} = L\frac{d}{\tau_{\mathrm{e}}}\,\mathrm{e}^{\frac{d}{\tau_{\mathrm{e}}}} \end{eqnarray} \tag{ C.1 }$$

as

$$\begin{eqnarray*} z_{k} & = & -\frac{1}{\tau_{\mathrm{e}}}+\frac{1}{d}W_{k}\left(L\frac{d}{\tau_{\mathrm{e}}}\,\mathrm{e}^{\frac{d}{\tau_{\mathrm{e}}}}\right), \end{eqnarray*}$$

leading to (8). The Lambert W_k(x) function has infinitely many branches k [45]. The principal branch has two real solutions, if x > −e⁻¹. The remaining branches appear in conjugate pairs. For x < −e⁻¹ the principal solutions turn into a complex conjugate pair. This happens at a sufficiently strong negative coupling L or long delays d

$$\begin{eqnarray*} &&L\frac{d}{\tau_{\mathrm{e}}}\,\mathrm{e}^{\frac{d}{\tau_{\mathrm{e}}}} < -\mathrm{e}^{-1},\\ &&L < -\frac{\tau_{\mathrm{e}}}{d}\,\mathrm{e}^{-\frac{d}{\tau_{\mathrm{e}}}-1} \end{eqnarray*}$$

or

$$\begin{eqnarray*} &&\frac{d}{\tau_{\mathrm{e}}} >W_{0}\left(-\frac{1}{Le}\right). \end{eqnarray*}$$

The principal poles may assume positive real values, leading to oscillations. The condition under which this happens can be derived from (C.1). At the point of transition the pole can be written as z = iω_crit; it is a solution to (1 + iω_critτ_e) e^i ω_critd = L. In order for this equation to be fulfilled, the absolute value and the phase must be identical on both sides. The equation for the absolute value requires 1 + (ω_critτ_e)² = L². This means there are only oscillatory solutions if the magnitude of the feedback exceeds unity L < −1. Since the poles come in conjugate pairs, we can assume w.l.o.g. that ω_crit > 0. The condition for the absolute value hence reads

$$\begin{eqnarray} \omega_{\mathrm{crit}}\tau_{\mathrm{e}} & = & \sqrt{L^{2}-1}. \end{eqnarray} \tag{ C.2 }$$

This is the frequency of oscillation at the onset of the Hopf bifurcation. For strong feedback |L| ≫ 1 the frequency increases linearly with the magnitude of the feedback. The condition for the agreement of the phase angles reads ∠(1 + iω_critτ_e) + ω_critd = 0, so $\frac {\Im (1+\mathrm {i}\,\omega _{\mathrm {crit}}\tau _{\mathrm {e}})}{\Re (1+\mathrm {i}\,\omega _{\mathrm {crit}}\tau _{\mathrm {e}})}=-\frac {\Im (\mathrm {e}^{\mathrm {i}\,\omega _{\mathrm {crit}}d})}{\Re (\mathrm {e}^{\mathrm {i}\,\omega _{\mathrm {crit}}d})}=\tan \omega _{\mathrm {crit}}d$, which leads to tan ω_critd = −ω_critτ_e. This equation has a solution in $\frac {\pi }{2}\leqslant \omega _{\mathrm {crit}}d\leqslant \pi$. In the limit of the vanishing delay d → 0 the frequency goes to infinity, as the solution converges to $\omega _{\mathrm {crit}}d=\frac {\pi }{2}$. This corresponds to the frequency $f_{\mathrm {crit}}=\frac {1}{4d}$. Inserting (C.2) leads to $\tan (\frac {d}{\tau _{\mathrm {e}}}\sqrt {L^{2}-1})=-\sqrt {L^{2}-1}$, which can be solved for the critical delay

$$\begin{eqnarray} \frac{d}{\tau_{\mathrm{e}}} & = & \frac{\pi-\arctan(\sqrt{L^{2}-1})}{\sqrt{L^{2}-1}}, \end{eqnarray} \tag{ C.3 }$$

where we took care that the argument of the tangent is in $[\frac {\pi }{2},\pi ]$. So with (C.2) and (C.3) the oscillatory frequency at the transition can be related to the synaptic delay as

$$\begin{eqnarray*} 2\pi f_{\mathrm{crit}}d=\omega d & = & \pi-\arctan(\sqrt{L^{2}-1}). \end{eqnarray*}$$

In the non-oscillatory state all poles z_k (8) have a negative real part. The function U(z) = ((1 + zτ_e) e^zd − L)⁻¹ in (6) then has all poles in the left complex half-plane, ℜ(z_k) < 0 ∀k. We perform the Fourier back transform

$$\begin{eqnarray} u(t) & = & \frac{1}{2\pi}\int_{-\infty}^{\infty}U(\mathrm{i}\,\omega)\, \mathrm{e}^{\mathrm{i}\,\omega t}\, \mathrm{d}\omega\nonumber\\ & = & \frac{1}{2\pi \mathrm{i}}\oint_{E_{t}}U(z)\, \mathrm{e}^{zt}\, \mathrm{d}z, \end{eqnarray}\noindent \tag{ D.1 }$$

replacing the integration path by a closed contour E_t following the imaginary axis from −i∞ to i∞. In order to ensure convergence of the integral, for t < d we need ℜ(z) > 0, so we close E_t<d in infinity within the right half-plane, where the integrand vanishes. Since there are no poles in the right half-plane, for t < d the path E_t<d does not enclose any poles, so u(t) = 0. For t ⩾ d the path E_t⩾d must be closed in the left half-plane to ensure convergence of (D.1), so the residue theorem yields

$$\begin{eqnarray} &&u(t) = \Theta(t-d)\sum_{z_{k}\in E_{t>d}}\mathrm{Res}(U,z_{k})\, \mathrm{e}^{z_{k}t}. \end{eqnarray}\noindent \tag{ D.2 }$$

The residue can be calculated by linearizing the denominator of U(z_k + z) around z_k

$$\begin{eqnarray*} (1+(z_{k}+z)\tau_{\mathrm{e}})\,\mathrm{e}^{(z_{k}+z)d}-L & = & (1+z_{k}\tau_{\mathrm{e}})\,\mathrm{e}^{z_{k}d}\,\mathrm{e}^{zd}-L+z\tau_{\mathrm{e}}\,\mathrm{e}^{z_{k}d}\,\mathrm{e}^{zd}\\ & = & L(\mathrm{e}^{zd}-1)+z\tau_{\mathrm{e}}\,\mathrm{e}^{z_{k}d}(1+zd)+O(z^{2})\\ & = & z(Ld+\tau_{\mathrm{e}}\,\mathrm{e}^{z_{k}d})+O(z^{2}), \end{eqnarray*}\noindent$$

which yields

$$\begin{eqnarray*} \mathrm{Res}(U,z_{k}) & = & \lim_{z\rightarrow0}zU(z_{k}+z)\\ & = & \lim_{z\rightarrow0}\frac{z}{z(Ld+\tau_{\mathrm{e}}\,\mathrm{e}^{z_{k}d})}=\frac{1}{Ld+\tau_{\mathrm{e}}\,\mathrm{e}^{z_{k}d}}\\ & = & \frac{\mathrm{e}^{-z_{k}d}}{(1+z_{k}\tau_{\mathrm{e}})d+\tau_{\mathrm{e}}}, \end{eqnarray*}\noindent$$

where we used (C.1) in the last step. The poles of V (z) = U(z)U(−z) are located in both half-planes, consequently v(t) is non-zero on the whole time axis. Here we only calculate v(t) for positive times t > 0, because it follows for negative times by symmetry v(−t) = v(t). The path has to be closed in the left half-plane, where the poles z_k have the residues

$$\begin{eqnarray*} \mathrm{Res}(V,z_{k})=\mathrm{Res}(U,z_{k})U(-z_{k}) & = & \frac{1}{\left(1+z_{k}\tau_{\mathrm{e}}\right)d+\tau_{\mathrm{e}}}\;\frac{1}{(1-z_{k}\tau_{\mathrm{e}})-L\,\mathrm{e}^{z_{k}d}}. \end{eqnarray*}\noindent$$

So by applying (D.2) the functions u and v are

$$\begin{eqnarray} u(t) & = & \sum_{_{k=0}}^{\infty}\frac{1}{\left(1+z_{k}\tau_{\mathrm{e}}\right)d+\tau_{\mathrm{e}}}\Theta(t-d)\,\mathrm{e}^{z_{k}(t-d)},\nonumber\\ v(t) & = & \sum_{k=0}^{\infty}\mathrm{Res}(V,z_{k})\left(\Theta(t)\,\mathrm{e}^{z_{k}t}+\Theta(-t)\,\mathrm{e}^{-z_{k}t}\right)\nonumber \\ & = & \sum_{k=0}^{\infty}\frac{1}{\left(1+z_{k}\tau_{\mathrm{e}}\right)d+\tau_{\mathrm{e}}}\,\frac{1}{(1-z_{k}\tau_{\mathrm{e}})-L\,\mathrm{e}^{z_{k}d}}\,\mathrm{e}^{z_{k}|t|}. \end{eqnarray}\noindent \tag{ D.3 }$$

The amplitude of the modes decreases with k. For all figures in the paper we truncated the series after k = 30.

All network simulations were performed using NEST [52]. The parameters of the leaky integrate-and-fire neuron model (1) throughout this work are τ_m = 20 ms, τ_s = 2 ms, τ_r = 2 ms, V_θ = 15 mV, V_r = 0 mV. All simulations are performed with precise spike timing and time stepping of 0.1 ms [53]. Figures 1 and 3–6 of the main text all consider recurrent random networks of N excitatory and γN inhibitory leaky integrate-and-fire model neurons receiving input from randomly drawn neurons in the network and external excitatory and inhibitory Poisson input, so that the first and second moments are μ_i = 15 mV and σ²_i = 10 mV, respectively. Unless stated explicitly, we use N(1 + γ) = 10 000 neurons except in figure 3 where the number of neurons is given in the legend. Each neuron has K = pN incoming excitatory synapses with synaptic amplitude J independently and randomly drawn from the pool of excitatory neurons, and γK = γpN inhibitory synapses with amplitude −gJ (a homogeneous Erdős–Rényi random network with fixed in-degree), realizing a connection probability p = 0.1. Cross covariance functions are measured throughout as the covariance between two disjoint populations of 1000 neurons each taken from the indicated populations in the network. Correlation functions are evaluated with a time resolution of 0.1 ms. In figure 1 we use a synaptic delay d = 1 ms.

In figure 3 we keep the feedback of the population rate constant L = Kw(1 − γg) = const. Increasing the size of the network N the synaptic amplitude J (which is proportional to w in the linear approximation) needs to scale as J = J₀N₀/N, where we chose J₀ = 0.1 mV and N₀ = 10 000 here. The variance caused by local input from the network then decreases with increasing network size ∝1/N, while the local mean is constant because L = const. Each cell receives in addition uncorrelated external balanced Poisson input, adjusted to keep the mean μ_i = 15 mV, and fluctuations σ_i = 10 mV constant. This is achieved by choosing the rates of the external excitatory (r_e,ext, amplitude J_ext = 0.1 mV) and inhibitory (r_i,ext, amplitude −gJ_ext) inputs as

$$\begin{eqnarray} &&\fl r_{\mathrm{e},{\rm ext}} = r_{\mathrm{e},0}+r_{\mathrm{bal}},\quad r_{\mathrm{i},{\rm ext}}=r_{\mathrm{bal}}/g,\nonumber\\ &&\textrm{with }r_{\mathrm{e},0} = \frac{\mu_{i}-\mu_{\mathrm{loc}}}{J_{{\rm ext}}\tau_{\mathrm{m}}}\quad\textrm{and }\quad r_{\mathrm{bal}}=\frac{\sigma_{i}^{2}-\sigma_{\mathrm{loc}}^{2}-\tau_{\mathrm{m}} r_{\mathrm{e,0}}J_{{\rm ext}}^{2}}{\tau_{\mathrm{m}} J_{{\rm ext}}^{2}(1+g^{2})}, \end{eqnarray}\noindent \tag{ E.1 }$$

where μ_loc = τ_mrKJ(1 − γg) and σ²_loc = τ_mrKJ²(1 + γg²) are the mean and variance due to local input from other neurons of the network firing with rate r. From L = const. follows that also Kw = const., so that (7) predicts a scaling of the magnitude of the covariance functions in proportion to 1/N. Other network parameters are d = 3 ms and g = 5. The firing rate in the network is r = 23.6 Hz.

In figures 4 and 5 we use N = 10⁴, g = 6, J = 0.1 mV and the delay d as described in the captions; the remaining parameters are as in figure 3.

In figure 6 we use N = 10⁴, J = 0.1 mV, g = 5 and d = 2 ms and the remaining parameters as in figure 3. We obtain the filtered synaptic currents by filtering the spike trains with an exponential filter of time constant τ_s = 2 ms. This results in an effective filter for the cross covariances of $q(t)=\frac {\tau _{\mathrm {m}}^{2}}{2\tau _{\mathrm {s}}}\,\mathrm {e}^{-|t|/\tau _{\mathrm {s}}}$. The different contributions shown are c_IeIe = q*(J²pKa_e + J²K²c_ee), c_IiIi = q*(J²pKa_i + J²K²c_ii), c_IeIi = −q*(gJ²γK²c_ei) and c_IiIe = −q*(gJ²γK²c_ie), where * denotes the convolution.

For figure 2, we simulate two populations of N = 1000 neurons each. Each neuron receives independent background activity from Poisson processes and in addition input from a common Poisson process with rate r_c = 25 Hz causing in population 1 a positive synaptic amplitude of J and for population 2 a negative synaptic amplitude J (J is given on the x-axis). The synaptic amplitude of the background inputs is J_e = 0.1 mV for an excitatory impulse and J_i = −0.5 mV for an inhibitory impulse. The rates of the excitatory and inhibitory background inputs are chosen so that the first and second moments μ_i = τ_m(J_er_e + J_ir_i + J r_c) = 15 mV and σ²_i = τ_m(J²_er_e + J²_ir_i + J r²_c) = 10 mV are independent of J. The spikes produced by each population are triggered to the arrival of an impulse in the common input and averaged over a duration of 10 s to obtain the impulse response.