Impact of Network Structure and Cellular Response on Spike Time Correlations

Novel experimental techniques reveal the simultaneous activity of larger and larger numbers of neurons. As a result there is increasing interest in the structure of cooperative – or correlated – activity in neural populations, and in the possible impact of such correlations on the neural code. A fundamental theoretical challenge is to understand how the architecture of network connectivity along with the dynamical properties of single cells shape the magnitude and timescale of correlations. We provide a general approach to this problem by extending prior techniques based on linear response theory. We consider networks of general integrate-and-fire cells with arbitrary architecture, and provide explicit expressions for the approximate cross-correlation between constituent cells. These correlations depend strongly on the operating point (input mean and variance) of the neurons, even when connectivity is fixed. Moreover, the approximations admit an expansion in powers of the matrices that describe the network architecture. This expansion can be readily interpreted in terms of paths between different cells. We apply our results to large excitatory-inhibitory networks, and demonstrate first how precise balance – or lack thereof – between the strengths and timescales of excitatory and inhibitory synapses is reflected in the overall correlation structure of the network. We then derive explicit expressions for the average correlation structure in randomly connected networks. These expressions help to identify the important factors that shape coordinated neural activity in such networks.


Introduction
New multielectrode and imaging techniques are revealing the simultaneous activity of neural ensembles and, in some cases, entire neural populations [1][2][3][4]. This has thrust upon the computational biology community the challenge of characterizing a potentially complex set of interactions -or correlations -among pairs and groups of neurons.
Beyond important and rich challenges for statistical modeling [5], the emerging data promises new perspectives on the neural encoding of information [6]. The structure of correlations in the activity of neuronal populations is of central importance in understanding the neural code [7][8][9][10][11][12][13]. However, theoretical [9][10][11][14][15][16], and empirical studies [17][18][19] do not provide a consistent set of general principles about the impact of correlated activity. This is largely because the presence of correlations can either strongly increase or decrease the fidelity of encoded information depending on both the structure of correlations across a population and how their impact is assessed.
A basic mechanistic question underlies the investigation of the role of collective activity in coding and signal transmission: How do single-cell dynamics, connection architecture, and synaptic dynamics combine to determine patterns of network activity? Systematic answers to this question would allow us to predict how empirical data from one class of stimuli will generalize to other stimulus classes and recording sites. Moreover, a mechanistic understanding of the origin of correlations, and knowledge of the patterns we can expect to see under different assumptions about the underlying networks, will help resolve recent controversies about the strength and pattern of correlations in mammalian cortex [1,[20][21][22]. Finally, understanding the origin of correlations will inform the more ambitious aim of inferring properties of network architecture from observed patterns of activity [23][24][25].
Here, we examine the link between network properties and correlated activity. We develop a theoretical framework that accurately predicts the structure of correlated spiking that emerges in a widely used model -recurrent networks of general integrate and fire cells. The theory naturally captures the role of single cell and synaptic dynamics in shaping the magnitude and timescale of spiking correlations. We focus on the exponential integrate and fire model, which has been shown to capture membrane and spike responses of cortical neurons [26]; however, the general approach we take can be applied to a much broader class of neurons, a point we return to in the Discussion.
Our approach is based on an extension of linear response theory to networks [25,27]. We start with a linear approximation of a neuron's response to an input. This approximation can be obtained explicitly for many neuron models [28][29][30], and is directly related to the spike triggered average [31]. The correlation structure of the network is then estimated using an iterative approach. As in prior work [32][33][34], the resulting expressions admit an expansion in terms of paths through the network.
We apply this theory to networks with precisely balanced inhibition and excitation in the inputs to individual cells. In this state individual cells receive a combination of excitatory and inhibitory inputs with mean values that largely cancel. We show that, when timescales and strengths of excitatory and inhibitory connections are matched, only local interactions between cells contribute to correlations. Moreover, our theory allows us to explain how correlations are altered when precise tuning balance is broken. In particular, we show how strengthening inhibition may synchronize the spiking activity in the network. Finally, we derive results which allow us to gain an intuitive understanding of the factors shaping average correlation structure in randomly connected networks of neurons.

Network model and linear response theory
Our goal is to understand how the architecture of a network shapes the statistics of its activity. We show how correlations between spike trains of cells can be approximated using response characteristics of individual cells along with information about synaptic dynamics, and the structure of the network. We start by briefly reviewing linear response theory of neuronal responses [29,35,36], and then use it to approximate the correlation structure of a network.

Network model
To illustrate the results we consider a network of N nonlinear integrate-and-fire (IF) neurons with membrane potentials modeled by Here E i represents the mean synaptic input current from parts of the system not explicitly modeled.
A spike-generating current ψ(v i ) may be included to emulate the rapid onset of action potentials. Unless otherwise specified, we utilize the exponential IF model (EIF), so that [26]. Cells are subject to internally induced fluctuations due to channel noise [37], and externally induced fluctuations due to inputs not explicitly modelled [38]. We model both by independent, Gaussian, white noise processes, [39]. An external signal to cell i is represented by η i (t).
Upon reaching a threshold v th , an action potential is generated, and the membrane potential is reset to v r , where it is held constant for an absolute refractory period τ ref . The output of cell i is characterized by the times, t i,k , at which its membrane potential reaches threshold, resulting in an output spike train y i (t) = k δ(t − t i,k ). Synaptic interactions are modeled by delayed α-functions The N × N matrix J contains the synaptic kernels, while the matrix W contains the synaptic weights, and hence defines the network architecture. In particular, W ij = 0 represents the absence of a synaptic connection from cell j to cell i. Table 1 provides an overview of all parameters and variables.

Measures of spike time correlation
We quantify dependencies between the responses of cells in the network using the spike train autoand cross-correlation functions [40]. For a pair of spike trains, y i (t), y j (t), the cross-correlation function C ij (τ ) is defined as C ij (τ ) = cov (y i (t + τ ), y j (t)) .
The auto-correlation function C ii (t) is the cross-correlation between a spike train and itself, and C(t) is the matrix of cross-correlation functions. Denoting by N y i (t 1 , t 2 ) = t 2 t 1 y i (s)ds the number of spikes over a time window [t 1 , t 2 ], the spike count correlation, ρ ij (T ), over windows of length τ is defined as, .
We assume stationarity of the spiking processes (that is, the network has reached a steady state) so that ρ ij (T ) does not depend on t. We also use the total correlation coefficient ρ ij (∞) = lim T →∞ ρ ij (T ) to characterize dependencies between the processes y i and y j over arbitrarily long timescales.
The spike count covariance is related to the cross-correlation function by [7,41] cov We can interpret the cross-correlation as the conditional probability that cell j spikes at time t + τ given that cell i spiked at time t. The conditional firing rate, is the firing rate of cell j conditioned on a spike in cell i at τ units of time in the past, and

Linear response of individual cells
Neuronal network models are typically described by a complex system of coupled nonlinear stochastic differential equations. Their behavior is therefore difficult to analyze directly. We will use linear response theory [29,35,36,40] to approximate the cross-correlations between the outputs of neurons in a network. We first review the linear approximation to the response of a single cell. We illustrate the approach using current-based IF neurons, and explain how it can be generalized to other models in the Discussion.
The membrane potential of an IF neuron receiving input X(t), with vanishing temporal average, X(t) = 0, evolves according to The time-dependent firing rate, r(t), is determined by averaging the resulting spike train, y(t) = j δ(t − t j ), across different realizations of noise, ξ(t), for fixed X(t). Using linear response theory, we can approximate the firing rate by where r 0 is the (stationary) firing rate when = 0. The linear response kernel, A(t), characterizes the firing rate response to first order in . A rescaling of the function A(t) gives the spike-triggered average of the cell, to first order in input strength, and is hence equivalent to the optimal Weiner kernel in the presence of the signal ξ(t). [40,42]. In Fig. 1, we compare the approximate firing rate obtained from Eq. (4) to that obtained numerically from Monte Carlo simulations.
The linear response kernel A(t) depends implicitly on model parameters, but is independent of the input signal, X(t), when is small relative to the noise √ σ 2 τ ξ(t). In particular, A(t) is sensitive to the value of the mean input current, E. We emphasize that the presence of the background noise, ξ, in Eq. (3) is essential to the theory, as noise linearizes the transfer function that maps input to output.

Linear response in recurrent networks
The linear response kernel can be used to approximate the response of a cell to an external input. However, the situation is more complicated in a network where a neuron can affect its own activity through recurrent connections. To extend the linear response approximation to networks we follow the approach introduced by [27]. Instead of using the linear response kernel to approximate the firing rate of a cell, we use it to approximate a realization of its output Here y 0 (t) represents a realization of the spike train generated by an integrate-and-fire neuron obeying Eq. (3) with X(t) = 0.
The central assumption we make is that a cell acts approximately as a linear filter of its inputs. Note that Eq. (5) defines a mixed point and continuous process, but averaging y(t) in Eq. (5) over realizations of y 0 gives Eq. (4). Hence, Eq. (5) can be viewed as a natural generalization of Eq. (4) where the unperturbed output of the cell is represented as a point process, y 0 (t), instead of the firing rate, r 0 .
We first use Eq. (5) to describe spontaneously evolving networks where η i (t) = 0. Equation (1) can then be rewritten as where represents the temporal average.
As a first approximation of the spiking output of cells in the coupled network, we start with realizations of spike trains, y 0 i , generated by IF neurons obeying Eq.
This is equivalent to considering neurons isolated from the network, with adjusted DC inputs (due to mean network interactions). Following the approximation given by Eq. (5), we use a frozen realization of all y 0 i to find a correction to the output of each cell, with X(t) set to the mean-adjusted synaptic input, As noted previously, the linear response kernel is sensitive to changes in the mean input current. It is therefore important to include the average synaptic input E[f i ] in the definition of the effective mean input, E i .
The input from cell j to cell i is filtered by the synaptic kernel J ij (t). The linear response of cell i to a spike in cell j is therefore captured by the interaction kernel K ij defined by The output of cell i in response to mean-adjusted input, y 0 j (t) − r j , from cell j can be approximated to first order in input strength using the linear response correction We explain how to approximate the stationary rates, r j , in the Methods section.
The cross-correlation between the processes y 1 i (t) in Eq. (7) gives a first approximation to the cross-correlation function between the cells (See Methods), where we used f − (t) = f (−t). [25] obtained an approximation closely related to Eq. (8). They first obtained the cross-correlation between a pair of neurons which either receive a common input or share a monosynaptic connection. This can be done using Eq. (4), without the need to introduce the mixed process given in Eq. (5). [25] then implicitly assumed that the correlations not due to one of these two submotifs could be disregarded. The correlation between pairs of cells which were mutually coupled (or were unidirectionally coupled with common input) was approximated by the sum of correlations introduced by each submotif individually.
Equation (8) provides a first approximation to the joint spiking statistics of cells in a recurrent network. However, it captures only the effects of direct synaptic connections, represented by the second and third terms, and common input, represented by the last term in Eq. (8). The impact of larger network structures, such as loops and chains are not captured, although they may significantly impact cross-correlations [43][44][45]. Experimental studies have also shown that local cortical connectivity may not be fully random [46][47][48]. It is therefore important to understand the effects on network architecture on correlations.
To capture the impact of the full network structure, we propose an iterative approach which accounts for successively larger connectivity patterns in the network [33,34]. We again start with y 0 i (t), a realization of a single spike train in isolation. Successive approximations to the output of cells in a recurrent network are defined by To compute the correction to the output of a neuron, in the first iteration we assume that its inputs come from a collection of isolated cells: When n = 1, Eq. (9) takes into account only inputs from immediate neighbors, treating each as disconnected from the rest of the network. The corrections in the second iteration are computed using the approximate cell responses obtained from the first iteration. Thus, with n = 2, Eq. (9) also accounts for the impact of next nearest neighbors. Successive iterations include the impact of directed chains of increasing length: The isolated output from an independent collection of neurons is filtered through n stages to produce the corrected response (See Fig. 2

.)
Notation is simplified when this iterative construction is recast in matrix form 1 to obtain y n+1 (t) = y 0 (t) + (K * [y n − r])(t) where y n (t) = [y n i (t)] and r = [r i ] are length N column vectors, and K (k) represents a k-fold matrix convolution of K with itself.
and Y(t) = [Yij(t)] be n1 × n2 and n2 × n3 matrices of functions, respectively. We define the convolution of matrices (X * Y)(t) to be the n1 × n3 matrix of functions with entries defined by Expectations and convolutions commute for matrix convolutions as matrix expectations are taken entry-wise. Each entry of a matrix convolution is a linear combination of scalar convolutions which commute with expectations. Additionally, we adopt the convention that the zeroth power of the interaction matrix, K 0 ij (t), is the diagonal matrix with K 0 ij (t) = δ(t) when i = j. Hence K 0 ij (t) acts as the identity matrix under matrix convolution.

6
The n th approximation to the matrix of cross-correlations can be written in terms of the interaction kernels, K ij , and the autocorrelations of the base processes y 0 as (See Methods) where K − (t) = K(−t), X (kT ) = (X (k) ) T , and X (k) is the k-fold matrix convolution of X with itself.
If we apply the Fourier transform, , we find that for each ω, where X * denotes the conjugate transpose of X. The zero-mean Fourier transformsỹ n i of the spiking processes y n i is defined byỹ n i = F[y n i − r i ], andf = F(f ) for all other quantities. For a suitable matrix norm || · ||, when ||K|| < 1, we can take the limit n → ∞ in Eq. (12) [49], to obtain an approximation to the full array of cross-spectrã This equation can also be obtained by generalizing the approach of [27] (also see [13]). In the limit n → ∞, directed paths of arbitrary length contribute to the approximation. Equation (13) therefore takes into account the full recurrent structure of the network. We will use the spectral norm || · || 2 , and assume that in the networks we study ||K|| 2 < 1. This condition is confirmed numerically when we use Eq. (13).
Finally, consider the network response to external signals, η i (t), with zero mean and finite variance. The response of the neurons in the recurrent network can be approximated iteratively by . External signals and recurrent synaptic inputs are both linearly filtered to approximate a cell's response, consistent with a generalization of Eq. (4). As in Eq. (11), the n th approximation to the matrix of correlations is where C η (τ ) = E η(t + τ )η(t) T is the covariance matrix of the external signals. We can again take the Fourier transform and the limit n → ∞, and solve forC(ω). If ||K|| < 1, When the signals comprising η are white (and possibly correlated) corrections must be made to account for the change in spectrum and response properties of the isolated cells [27,50,51] (See Methods).
We note that Eq. (10), which is the basis of our iterative approach, provides an approximation to the network's output which is of higher than first order in connection strength. This may seem at odds with a theory that provides a linear correction to a cell's response, cf. Eq. (4). However, Eq. (10) does not capture nonlinear corrections to the response of individual cells, as the output of each cell is determined linearly from its input. It is the input that can contain terms of any order in connection strength stemming from directed paths of different lengths through the network.

Results
We use the theoretical framework developed above to analyze the statistical structure of the spiking activity in a network of IF neurons described by Eq. (1). We first show that the cross-correlation functions between cells in two small networks can be studied in terms of contributions from directed paths through the network. We use a similar approach to understand the structure of correlations in larger all-to-all and random networks. We show that in networks where inhibition and excitation are tuned for exact balance, only local interactions contribute to correlations. When such balance is broken by a relative elevation of inhibition, the result may be increased synchrony in the network. The theory also allows us to obtain averages of cross-correlation functions conditioned on connectivity between pairs of cells in random networks. Such averages can provide a tractable yet accurate description of the joint statistics of spiking in these networks.
The correlation structure is determined by the response properties of cells together with synaptic dynamics and network architecture. Network interactions are described by the matrix of synaptic filters, J, given in Eq. (2), while the response of cell i to an input is approximated using its linear response kernel A i . Synaptic dynamics, architecture, and cell responses are all combined in the matrix K, where K ij describes the response of cell i to an input from cell j (See Eq. (1)). The correlation structure of network activity is approximated in Eq. (13) using the Fourier transforms of the interaction matrix, K, and the matrix of unperturbed autocorrelations C 0 .

Statistics of the response of microcircuits
We first consider a pair of simple microcircuits to highlight some of the features of the theory. We start with the three cell model of feed-forward inhibition shown in Fig. 3A [52]. The interaction matrix,K(ω), has the formK where cells are indexed in the order E 1 , E 2 , I. To simplify notation, we omit the dependence of K(ω) and other spectral quantities on ω.
Note thatK is nilpotent, and the inverse of (I −K) may be expressed as Substituting Eq. (15) into Eq. (13) yields an approximation to the matrix of cross-spectra. For instance,C

III
. Figure 3B shows that these approximations closely match numerically obtained cross-correlations. C 0 X is the uncoupled power spectrum for cell X.
Equation (16) gives insight into how the joint response of cells in this circuit is shaped by the features of the network. The three terms in Eq. (16) are directly related to the architecture of the microcircuit: Term I represents the correlating effect of the direct input to cell E 2 from cell I. Term II captures the effect of the common input from cell E 1 . Finally, term III represents the interaction of the indirect input from E 1 to E 2 through I with the input from E 1 to I (See Fig. 3C). A change in any single parameter may affect multiple terms. However, the individual contributions of all three terms are apparent.
To illustrate the impact of synaptic properties on the cross-correlation between cells E 2 and I we varied the inhibitory time constant, τ I (See Fig. 3B and C). Such a change is primarily reflected in the shape of the first order term, I: Multiplication byJ E 2 I is equivalent to convolution with the inhibitory synaptic filter, J E 2 I . The shape of this filter is determined by τ I (See Eq. (2)), and a shorter time constant leads to a tighter timing dependency between the spikes of the two cells [25,[53][54][55][56]. In particular, Ostojic et al. made similar observations using a related approximation. In the FFI circuit, the first and second order terms, I and II, are dominant (red and dark orange, Fig. 3B). The relative magnitude of the third order term, III (light orange, Fig. 3B), is small. The next example shows that even in a simple recurrent circuit, terms of order higher than two may be significant.
More generally, the interaction matrices,K, of recurrent networks are not nilpotent. Consider two reciprocally coupled excitatory cells, E 1 and E 2 (See Fig. 4A, left). In this case, Equation (13) gives the following approximation to the matrix of cross-spectrã In contrast to the previous example, this approximation does not terminate at finite order in interaction strength. After expanding, the cross-spectrum between cells E 1 and E 2 is approximated byC Directed paths beginning at E 1 and ending at E 2 (or vice-versa) are of odd length. Hence, this approximation contains only odd powers of the kernelsK E i E j , each corresponding to a directed path from one cell to the other. Likewise, the approximate power spectra contain only even powers of the kernels corresponding to directed paths that connect a cell to itself (See Fig. 4A).
The contributions of different sub-motifs to the cross-and auto-correlations are shown in Figs. 4C, D when the isolated cells are in a near-threshold excitable state (CV ≈ 0.98). The auto-correlations are significantly affected by network interactions. We also note that chains of length two and three (the second and third submotifs in Fig. 4A) provide significant contributions. Earlier approximations do not capture such corrections [25].
The operating point of a cell is set by its parameters (τ i , E L,i , etc.) and the statistics of its input (E i , σ i ). A change in operating point can significantly change a cell's response to an input. Using linear response theory, these changes are reflected in the response functions A i , and the power spectra of the isolated cells,C 0 . To highlight the role that the operating point plays in the approximation of the correlation structure given by Eq. (13), we elevated the mean and decreased the variance of background noise by increasing E i and decreasing σ i in Eq. (1). With the chosen parameters the isolated cells are in a super-threshold, low noise regime and fire nearly periodically (CV ≈ 0.31). After the cells are coupled, this oscillatory behavior is reflected in the cross-and autocorrelations where the dominant contributions are due to first and zeroth order terms, respectively (See Figs. 4F,G).
Orders of coupling interactions: It is often useful to expand Eq. (13) in terms of powers of K [32]. The termK nC0 (K * ) m in the expansion is said to be of order n + m. Equivalently, in the expansion ofC ∞ ij , the order of a term refers to the sum of the powers of all constituent interaction kernelsK ab . We can also associate a particular connectivity submotif with each term. In particular, n th order terms of the formK ia n−1K a n−1 a n−2 · · ·K a 1 jC 0 jj are associated with a directed path j → a 1 → · · · → a n−2 → a n−1 → i from cell j to cell i. Similarly, the termC 0 iiK * ia 1 · · ·K * a n−2 a n−1K * a n−1 j corresponds to a n-step path from cell i to cell j. An (n + m) th order term of the form K ia n−1K a n−1 a n−2 · · ·K a 1 a 0C 0 a 0 a 0K * represents the effects of an indirect common input n steps removed from cell i and m steps removed from cell j. This corresponds to a submotif of the form i ← a n−1 ← · · · ← a 0 → b 1 → · · · → b n−1 → j consisting of two branches originating at cell a 0 . (See Fig. 5, and also Fig. 6A and the discussion around Eqs. (16,18).)

Statistics of the response of large networks
The full power of the present approach becomes evident when analyzing the activity of larger networks. We again illustrate the theory using several examples. In networks where inhibition and excitation are tuned to be precisely balanced, the theory shows that only local interactions contribute to correlations. When this balance is broken, terms corresponding to longer paths through the network shape the cross-correlation functions. One consequence is that a relative increase in inhibition can lead to elevated network synchrony. We also show how to obtain tractable and accurate approximation of the average correlation structure in random networks.
A symmetric, all-to-all network of excitatory and inhibitory neurons We begin with an all-to-all coupled network of N identical cells. Of these cells, N E make excitatory, and N I make inhibitory synaptic connections. The excitatory cells are assigned indices 1, . . . , N E , and the inhibitory cells indices N E + 1, . . . , N . All excitatory (inhibitory) synapses have weight The interaction matrixK may then be written in block form, Here 1 N 1 N 2 is the N 1 × N 2 matrix of ones,J X is the weighted synaptic kernel for cells of class X (assumed identical within each class), andÃ is the susceptibility function for each cell in the network. Although the effect of autaptic connections (those from a cell to itself) is negligible (See SI Fig. 2, their inclusion significantly simplifies the resulting expressions.
We defineμ E = N EJE ,μ I = N IJI , andμ =μ E +μ I . Using induction, we can show that Direct matrix multiplication yields which allows us to calculate the powersK kKl * when k, l = 0, An application of Eq. (13) then gives an approximation to the matrix of cross-spectra: The cross-spectrum between two cells in the network is therefore given by where X ∈ {E, I}. In Eq. (20) the first two terms represent the effects of all unidirectional chains originating at cell j and terminating at cell i, and vice versa. The third term represents the effects of direct and indirect common inputs to the two neurons. In Fig. 6A, we highlight a few of these contributing motifs.
Interestingly, when excitation and inhibition are tuned for precise balance (so thatμ =μ E +μ I = 0), Eq. (20) reduces to Effects of direct connections between the cells are captured by the first two terms, while those of direct common inputs to the pair are captured by the third term. Contributions from other paths do not appear. In other words, in the precisely balanced case only local interactions contribute to correlations.
To understand this cancelation intuitively, consider the contribution of directed chains originating at a given excitatory neuron, j. For τ > 0, the cross-correlation function, C ij (τ ), is determined by the change in firing rate of cell i at time τ given a spike in cell j at time 0. By the symmetry of the all-to-all connectivity and stationarity, the firing of cell j has an equal probability of eliciting a spike in any excitatory or inhibitory cell in the network. Due to the precise synaptic balance, the postsynaptic current generated by the elicited spikes in the excitatory population will cancel the postsynaptic current due to elicited spikes in the inhibitory population on average. The contribution of other motifs cancel in a similar way.
In Fig. 6B, we show the impact of breaking this excitatory-inhibitory balance on cross-correlation functions. We increased the strength and speed of the inhibitory synapses relative to excitatory synapses, while holding constant, for sake of comparison, the long window correlation coefficients, ρ(∞). at ≈ 0.05 Moreover, the degree of network synchrony, characterized by the short window correlation coefficients, is increased (See Fig. 6B inset). Intuitively, a spike in one of the excitatory cells transiently increases the likelihood of spiking in all other cells in the network. Since inhibition in the network is stronger and faster than excitation, these additional spikes will transiently decrease the likelihood of spiking in twice removed cells.
Linear response theory allows us to confirm this heuristic observation, and quantify the impact of the imbalance on second order statistics. Expanding Eq. (20) for two excitatory cells to second order in coupling strength, we find The complete cancellation between contributions of chains involving excitatory and inhibitory cells no longer takes place, and the two underlined terms appear as a consequence (compare with Eq. (21)). These underlined terms capture the effects of all length two chains between cells E i or E j , starting at one and terminating at the other. The relative strengthening of inhibition implies that chains of length two provide a negative contribution to the cross-correlation function at short times (cf. [57]). Additionally, the impact of direct common input to cells E i and E j on correlations is both larger in magnitude (because we increased the strength of both connection types) and sharper (the faster inhibitory time constant means common inhibitory inputs induces sharper correlations). These changes are reflected in the second order term |Ã| 2μ c in Eq. (22). In sum, unbalancing excitatory and inhibitory connections via stronger, faster inhibitory synapses enhances synchrony, moving a greater proportion of the covariance mass closer to τ = 0 (See Fig. 6B). To illustrate this effect in terms of underlying connectivity motifs, we show the contributions of length two chains and common input in both the precisely tuned and non-precisely tuned cases in Fig. 6C. A similar approach would allow us to understand the impact of a wide range of changes in cellular or synaptic dynamics on the structure of correlations across networks.
Random, fixed in-degree networks of homogeneous excitatory and inhibitory neurons Connectivity in cortical neuronal networks is typically sparse, and connection probabilities can follow distinct rules depending on area and layer [58]. The present theory allows us to consider arbitrary architectures, as we now illustrate.
We consider a randomly connected network of N E excitatory and N I inhibitory cells coupled with probability p. To simplify the analysis, every cell receives exactly pN E excitatory and pN I inhibitory inputs. Thus, having fixed in-degree, each cell receives an identical level of mean synaptic input. In addition, we continue to assume that cells are identical. Therefore, the response of each cell in the network is described by the same linear response kernel. The excitatory and inhibitory connection strengths are G E /(pN E ) and G I /(pN I ), respectively. The timescales of excitation and inhibition may differ, but are again identical for cells within each class.
The approximation of network correlations (Eq. (13)) depends on the realization of the connectivity matrix. For a fixed realization, the underlying equations can be solved numerically to approximate the correlation structure (See Fig. 7A). However, the cross-correlation between a pair of cells of given types has a form which is easy to analyze when only leading order terms in 1/N are retained.
Specifically, the average cross-spectrum for two cells of given types is (See SI Section 1) when i = j. This shows that, to leading order in 1/N , the mean cross-spectrum between two cells in given classes equals that in the all-to-all network (see Eq. (20)). Therefore our previous discussion relating network architecture to the shape of cross-correlations in the all-to-all network extends to the average correlation structure in the random network for large N .
[32] derived similar expressions for the correlation functions in networks of interacting Hawkes processes [59,60] by assuming either the network is regular (i.e., both in-and out-degrees are fixed) or has a sufficiently narrow degree distribution. Our analysis depends on having fixed in-degrees, and we do not assume that networks are fully regular. Both approaches lead to results that hold approximately (for large enough N ) when the in-degree is not fixed.
Average correlations between cells in the random network conditioned on first order connectivity As Fig. 7B shows there is large variability around the mean excitatory-inhibitory cross-correlation function given by the leading order term of Eq. (23). Therefore, understanding the average cross-correlation between cells of given types does not necessarily provide much insight into 13 the mechanisms that shape correlations on the level of individual cell pairs. Instead, we examine the average correlation between a pair of cells conditioned on their first order (direct) connectivity.
We derive expressions for first order conditional averages correct to O(1/N 2 ) (See SI Sec. 2). The average cross-spectrum for a pair of cells with indices i = j, conditioned on the value of the direct connections between them is Here we setJ ij = 0 if we condition on the absence of a connection j → i, andJ ij =J Y /p if we condition on its presence. The termJ ji is set similarly.
Although Eq. (24) appears significantly more complicated than the cell-type averages given in Eq. (23), they only differ in the underlined, first order terms. The magnitude of expected contributions from all higher order motifs is unchanged and coincides with those in the all-to-all network. Figure 7C shows the mean cross-correlation function for mutually coupled excitatory-inhibitory pairs. Taking into account the mutual coupling significantly reduces variability (Compare with Fig. 7B). To quantify this reduction, we calculate the mean reduction in variability when correlation functions are computed conditioned on the connectivity between the cells. For a single network, the relative decrease in variability can be quantified using 2 where T represents pairs of cells of a given type and connection (in the present example these are reciprocally coupled excitatory-inhibitory pairs), N T is the number of pairs of that type in the network, C CT T (τ ) is the leading order approximation of average correlations given only the type of cells in T (as in Eq. (23)), and C FOC T (τ ) the leading order approximation to average correlations conditioned on the first order connectivity of class T (as in Eq. (24)). Figure 7D shows µ error averaged over twenty networks. In particular, compare the reduction in variability when conditioning on bidirectional coupling between excitatory-inhibitory pairs shown in Figs. 7B,C, with the corresponding relative error in Fig. 7D (circled in red).

Discussion
We have developed a general theoretical framework that can be used to describe the correlation structure in a network of spiking cells. This theory allows us to find tractable approximations of cross-correlation functions in terms of the network architecture and single cell response properties. The approach was originally used to study the population response of the electrosensory lateral line lobe of weakly electric fish [27]. The key approximation relies on the assumption that the activity of cells in the network can be represented by a mixed point and continuous stochastic process, as given in Eq. (8). An iterative construction then leads to the expressions for approximate cross-correlations between pairs of cells given by Eq. (13).
[25] obtained formulas for cross-correlations that correspond to the first step in this iterative construction, given in Eq. (8). Their approach captures corrections due to direct coupling (first order terms) and direct common input (second order terms involving second powers of interaction kernels). Our approach can be viewed as a generalization that also accounts for length two directed chains, along with all higher order corrections. As Fig. 4 illustrates, these additional terms can yield significant contributions to the structure of cross-correlations. We note that another distinction between our approach and that of [25] is that the present approach also allows us to calculate corrected auto-correlations, whereas the framework of [25] did not allow for adjustments to autocorrelations.
[32] analyzed the correlation structure in networks of interacting Hawkes processes [59,60] using an approach similar to the one presented here. In both cases, we represented correlations between cell pairs in terms of contributions of different connectivity motifs. However, our methods differ in important ways: while their expressions are exact for Hawkes processes, [32] did not attempt to match their results to more physiological cell models, and did not account for the response properties of individual cells, although it may be possible to approximately fit the Hawkes models to such models. [32] examined only total spike count covariances, which equal the integrals of the cross-correlation functions. This leads to a loss of information about the temporal structure of correlations. However, as they note, their approach can be extended to obtain complete crosscorrelation functions for the Hawkes model.
To illustrate the power of the theory in analyzing the factors that shape correlations, we considered a number of simple examples for which the approximation Eq. (13) is tractable. We showed how the theory can be used both to obtain intuition for the effects that shape correlations, and to quantify their impact. In particular, we explained how only local connections affect correlations in a precisely tuned all-to-all network, and how strengthening inhibition may synchronize spiking activity.
Linear response methods are perturbative. For Eq. (5) to be valid, neurons need to respond approximately linearly to their inputs. This will only be true if inputs are "weak" relative to the dynamical operating point of the cell. We assume the presence of a white noise background which linearizes the transfer function, although it is possible to extend the present methods to colored background noise [26,61].
It may be surprising that an approach based on linear response theory can provide corrections to cross-correlations of arbitrary order in network connectivity. Corrections to firing activity which are higher order in network connectivity emerge from the linear correction in Eq. (10). A full expansion of firing activity would include terms arising from corrections to the input-output transfer of the individual cells beyond those captured by the linear response approximation. Formally, including such terms would be necessary to capture all contributions of a given order in network connectivity [33,34]. However, the high accuracy demonstrated by our method indicates that, at least in some cases, these additional terms are quite small. In short, our approximation neglects higher-order corrections to the input-output transfer of individual cells, which is compatible with the assumption that the presence of background noise causes this transfer to be close to linear.
As expected from the preceding discussion, simulations suggest that, for IF neurons, our approximations become less accurate as cells receive progressively stronger inputs. The "physical" reasons for this loss of accuracy could be related to interactions between the "hard threshold" and incoming synaptic inputs with short timescales. While the theory will work for short synaptic timescales, it will improve for slower synaptic dynamics, limiting towards being essentially exact in the limit of arbitrarily long synaptic time constants (Note the improvement in the approximation for the FFI circuit for the slower timescale as exhibited in Fig. 3). However, as we have shown, the theory remains applicable in a wide range of dynamical regimes, including relatively low noise, superthreshold regimes where cells exhibit strong oscillatory behavior. Moreover, the theory can yield accurate approximations of strong correlations due to coupling: for the bidirectionally coupled excitatory circuit of Fig. 4, the approximate cross-correlations match numerically obtained results even when correlation coefficients are large (ρ E 1 E 2 (∞) ≈ 0.8 in the excitable regime, ≈ 0.5 in the oscillatory regime).
Although we have demonstrated the theory using networks of integrate-and-fire neurons, the approach is widely applicable. The linear response kernel and power spectrum for a general integrate and fire neuron model can be easily obtained [30]. In addition, it is also possible to obtain the rate, spectrum, and rate response function for modulation of the mean conductance in the case of conductance-based (rather than current-based) synapses (See [62] and SI Sec. 3). As the linear response kernel is directly related to the spike triggered average [25,31], the proposed theoretical framework should be applicable even to actual neurons whose responses are characterized experimentally.
The possibilities for future applications are numerous. For example, one open question is how well the theory can predict correlations in the presence of adaptive currents [62]. In addition, the description of correlations in terms of architecture and response properties suggests the possibility of addressing the difficult inverse problem of inferring architectural properties from correlations [23][24][25]63]. Overall, it is our hope that the present approach will prove a valuable tool in moving the computational neuroscience community towards a more complete understanding of the origin and impact of correlated activity in neuronal populations.

Methods
Numerical methods Simulations were run in C++, and the stochastic differential equations were integrated with a standard Euler method with a time-step of 0.01ms. General parameter values were as follows: 5mV, ∆ T = 1.4mV, τ E = 10ms, τ I = 5ms, τ D,i = 1ms. Marginal statistics (firing rates, uncoupled power spectra and response functions) were obtained using the threshold integration method of [30] in MATLAB. All code is available upon request.
Calculation of stationary rates in a recurrent network The stationary firing rate of an IF neuron can be computed as a function of the mean and intensity of internal noise (E i , σ i ) and other cellular parameters (τ i , E L i , etc...) [64]. Denote the stationary firing rate of cell i in the network by r i , and by r i,0 (E, σ) the stationary firing rate in the presence of white noise with mean E and variance σ 2 . We keep the dependencies on other parameters are implicit. The stationary rates, r i , in the recurrent network without external input are determined self-consistently by where we used E[f i ] = j W ij E[y j ] = j W ij r j . This equality holds because the synaptic kernels, J ij , were normalized to have area W ij . These equations can typically be solved by fixed-point iteration.
Note that this provides an effective mean input, E i , to each cell, but does not give adjustments to the variance, σ i . We assume that the major impact of recurrent input is reflected in E i , and ignore corrections to the cell response involving higher order statistics of the input. This approach is valid as long as fluctuations in the recurrent input to each cell are small compared to σ i , and may break down otherwise [28].
Correction to statistics in the presence of an external white noise signals. Expression (14) can be used to compute the statistics of the network response to inputs η i (t) of finite variance. As noted by [27], when inputs have infinite variance additional corrections are necessary. As a particular example, consider the case where the processes are correlated white noise, i.e., when , where x c , x i are independent white noise processes with variance σ e . Then each η i is also a white noise process with intensity σ The firing rate of cell i in response to this input is r i = r 0 (E i , (σ i ) 2 + (σ e i ) 2 ), and the point around which the response of the cell is linearized needs to be adjusted.
Finally, we may apply an additional correction to the linear response approximation of autocorrelations. For simplicity, we ignore coupling in Eq. (14) (so thatK = 0). Linear response predicts thatC ii (ω) =C 0 ii (ω; σ 2 i ) + (σ e i ) 2 |Ã i (ω)| 2 , where we have introduced explicit dependence on σ 2 i , the variance of white noise being received by an IF neuron with power spectrumC 0 ii (ω; σ 2 i ), in the absence of the external signal. The approximation may be improved in this case by making the following substitution in Eq. (14) [27,51]: The response function A should be adjusted likewise.    (18)). Though we show only a few "chain" motifs in one direction, one should note that there will also be contributions to the cross-correlation from chain motifs in the reverse direction in addition to indirect common input motifs (See the discussion of  17) in the excitable (D) and oscillatory (G) regimes. In the oscillatory regime, higher order contributions were small relative to first order contributions and are therefore not shown. The network's symmetry implies that cross-correlations are symmetric, and we only show them for positive times. All cross-correlations are nearly indistinguishable from those obtained from simulations (See SI Fig. 3 for comparisons with simulations). Connection strengths were 40 mV · ms. The long window correlation coefficient ρ E 1 E 2 (∞) between the two cells was ≈ 0.8 in the excitable regime and ≈ 0.5 in the oscillatory regime. The ISI CV was approximately 0.98 for neurons in the excitable regime and 0.31 for neurons in the oscillatory regime.  (Solid line -mean cross-correlation, shaded area -one standard deviation from the mean, calculated using bootstrapping in a single network realization). (C) Mean and standard deviation for the distribution of cross-correlation functions conditioned on cell type and first order connectivity for a reciprocally coupled excitatory-inhibitory pair of cells. (Solid line -mean cross-correlation function, shaded area -one standard deviation from the mean found by bootstrapping). (D) Average reduction in L 2 error between cross-correlation functions and their respective first-order conditioned averages, relative to the error between the cross-correlations and their cell-type averages. Blue circles give results for a precisely tuned network, and red squares for a network with stronger, faster inhibition. Error bars indicate two standard errors above and below the mean. G E , G I , τ E , τ I for panels A-C are as in the precisely tuned network of Fig. 6, and the two networks of panel D are as in the networks of the same figure. v i , τ i , E L,i , σ i Membrane potential, membrane time constant, leak reversal potential, and noise intensity of cell i.
Mean and standard deviation of the background noise for cell i.
v th , v r , τ ref Membrane potential threshold, reset, and absolute refractory period for cells.
ψ(v), V T , ∆ T Spike generating current, soft threshold and spike shape parameters for the IF model [26].

W ij
The j → i synaptic weight, describes the area under a single post-synaptic current for current-based synapses.
J ij (t) The j → i synaptic kernel -equals the synaptic weight W ij multiplied times the synaptic filter for outputs of cell j.
N y i (t, t + T ), ρ ij (T ) Spike count for cell i, and spike count correlation coefficient for cells i, j over windows of length T .
ii Stationary rate, linear response kernel and uncoupled auto-correlation function for cell ij.
K ij (t) The j → i interaction kernel -describes how the firing activity of cell i is perturbed by an input spike from cell j. It is defined by K ij (t) = (A i * J ij )(t).
y n i (t), C n ij (t) The n th order approximation of the activity of cell i in a network which accounts for directed paths through the network graph up to length n ending at cell i, and the cross-correlation between the n th order approximations of the activity of cells i, j. g(t),g(ω)g(ω) is the Fourier transform of g(t) with the conventioñ g(ω) = F[g](ω) ≡ ∞ −∞ e −2πiωt g(t)dt