The correlated state in balanced neuronal networks

Understanding the magnitude and structure of inter-neuronal correlations and their relationship to synaptic connectivity structure is an important and difficult problem in computational neuroscience. Early studies show that neuronal network models with excitatory-inhibitory balance naturally create very weak spike train correlations. Later work showed that, under some connectivity structures, balanced networks can produce larger correlations between some neuron pairs, even when the average correlation is very small. All of these previous studies assume that the local neuronal network receives feedforward synaptic input from a population of uncorrelated spike trains. We show that when spike trains providing feedforward input are correlated, the downstream recurrent neuronal network produces much larger correlations. We provide an in-depth analysis of the resulting “correlated state” in balanced networks and show that, unlike the asynchronous state of previous work, it produces a tight excitatory-inhibitory balance consistent with in vivo cortical recordings.


I. INTRODUCTION
Correlations between the spiking activity of cortical neurons have important consequences for neural dynamics and coding [1][2][3]. A quantitative understanding of how spike train correlations are generated and shaped by the connectivity structure of neural circuits is made difficult by the noisy and nonlinear dynamics of recurrent neuronal network models [4][5][6][7]. Linear response and related techniques have been developed to overcome some of these difficulties [8][9][10][11][12][13][14], but their accuracy typically require an assumption of sparse and/or weak connectivity and, in some models, an additional assumption that neurons receive uncorrelated, feedforward Gaussian white noise input. However, cortical circuits are densely connected and receive spatially and temporally correlated synaptic input from outside the local circuit [15][16][17][18].
An alternative approach to analyzing correlated variability in recurrent neuronal network models is motivated in part by the widely observed balance between excitatory and inhibitory synaptic inputs in cortex [19][20][21][22][23][24][25][26]. When synaptic weights are scaled like 1/ √ N where N is the size of a model network, a cortex-like balance between excitation and inhibition arises naturally at large network size, which defines the "balanced state" [27,28]. Early work on balanced networks assumed sparse connectivity to produce weak spike train correlations, but it was later shown that keeping connection probabilities O(1) naturally produces weak, O(1/N ), spike train correlations, defining the "asynchronous state" [29]. While these extremely weak spike train correlations are consistent with some cortical recordings [30], the magnitude of correlations in cortex can depend on stimulus, cortical area, layer, and behavioral or cognitive state, and can be much larger than predicted by the asynchronous state [6,[31][32][33][34][35]. This raises the question of how larger correlation magnitudes can arise in balanced cortical circuits. Later theoretical work showed that larger correlations can be obtained between some cell pairs in densely connected networks with specially constructed connectivity structure [36][37][38][39], offering a potential explanation of the larger correlations often observed in recordings. These previous theoretical studies of correlated variability in balanced networks assume that the recurrent network receives feedforward synaptic input from an external population of uncorrelated spike trains, so feedforward input correlations arise solely from overlapping feedforward synaptic projections. In reality, feedforward synaptic input to a cortical population comes from thalamic projections, other cortical areas, or other cortical layers in which spike trains could be correlated.
We extend the theory of densely connected balanced networks to account for correlations between the spike trains of neurons in an external, feedforward input layer. We show that correlations between the feedforward synaptic input to neurons in the recurrent network are O(N ) in this model, but cancel with O(N ) correlations between recurrent synaptic input to produce O(1) total input correlation and O(1) spike train correlations on average, defining what we refer to as the "correlated state" in densely connected balanced networks. This correlated state offers an alternative explanation for the presence of moderately large spike train correlations in cortical recordings. We derive a simple, closed form approximation for the average cross-spectral density between neurons' spike trains in the correlated state in term of synaptic parameters alone, without requiring the use of linear response theory or any other knowledge of neurons' transfer functions. We show that the tracking of excitatory synaptic input currents by inhibitory currents is more precise and more similar to in vivo recordings [22] in the correlated state than in the asynchronous state. We also investigate the applicability of linear response approximations to correlated variability in densely connected balanced networks. Taken together, our results extend the theory of correlated variability in balanced networks to the biologically realistic assumption that presynaptic neural populations are themselves correlated.

II. MODEL AND BACKGROUND
We consider recurrent networks of N integrate-andfire model neurons, N e of which are excitatory and N i inhibitory. Neurons are randomly and recurrently interconnected and also receive random feedforward synaptic input from an external population of N x neurons whose spike trains are homogeneous Poisson processes with rate r x (Fig. 1A).
The membrane potential of neuron j in population a = e, i obeys the exponential integrate-and-fire (EIF) dynamics with the added condition that each time V a j (t) exceeds V th , it is reset to V re and a spike is recorded. We additionally set a lower bound on the membrane potential at V lb = −100mV. Spike trains are represented as a sum of Dirac delta functions, where t a,j n is the nth spike time of neuron j in population a = e, i, x. The total synaptic input current to neuron j in population a = e, i is decomposed as T a j (t) = E a j (t) + I a j (t) + X a j (t) where U a j (t) = for U = E, I, X and b = e, i, x respectively where * denotes convolution, J ab jk is the synaptic weight from neuron k in population b to neuron j in population a, and α b (t) is a postsynaptic current (PSC) waveform. Without loss of generality, we assume that α b (t) = 1. We use where H(t) is the Heaviside step function, though our results do not depend sensitively on the precise neuron model or PSC kernel used. For calculations, it is useful to decompose the total synaptic input into its recurrent and external sources, where R a j (t) = E a j (t) + I a j (t) is the recurrent synaptic input from the local circuit.
Local cortical circuits contain a large number of neurons and individual cortical neurons receive synaptic input from thousands of other neurons within their local circuit and from other layers or areas. Densely connected balanced networks have been proposed to model such large and densely interconnected neuronal networks [29,38]. In such models, one considers the limit of large N (with N x , N e and N i scaled proportionally) with fixed connection probabilities and where synaptic weights are scaled like O(1/ √ N ) [28,29]. This scaling naturally captures the balance of mean excitatory and mean inhibitory synaptic input, as well as the tracking of excitation by inhibition, observed in cortical recordings [29]. In particular, we consider a random connectivity structure in which where connections are statistically independent and j ab , p ab ∼ O(1) for b = e, i, x and a = e, i. We furthermore define the proportions For all examples we consider, q e = 0.8 and q i = q x = 0.2. We next introduce notational conventions for quantifying the statistics of spike trains and synaptic inputs in the network. The mean firing rates of neurons in population a = e, i, x is defined by r a for a = e, i, x and it is useful to define the 2 × 1 vector, r = [r e r i ] T . The mean is technically interpreted as the expectation over realizations of the network connectivity, but for large N it is approximately equal to the sample mean over all neurons the network. Similarly, mean-field synaptic inputs to neurons in populations a = e, i are defined by for U = E, I, X, R, T and, in vector form, U = [U e U i ] T For quantifying correlated variability, we use the crossspectral density (CSD) is the cross-covariance function. The argument, f , is the frequency at which the CSD is evaluated. The CSD is a convenient measure of correlated variability because it simplifies mathematical calculations due to the fact that it is a Hermitian operator and because most commonly used measures of correlated variability can be written as a function of the CSD. For example, the cross-covariance is the inverse Fourier transform of the CSD. Spike count covariances over large time windows can be written in terms of the CSD by first noting that the spike count is an integral of the spike train [4], For large t 0 , the cross-spectrum between two integrals is related to the zero-frequency CSD, Hence, the spike count covariance over the interval [0, t 0 ] is given asymptotically by t 0 S a j , S b k (f = 0) for large t 0 . Following this result, we quantify covariability between spike trains and between synaptic currents using the zerofrequency CSD, which we estimate by taking the covariance between integrals as in Eq. (3) using t 0 = 250ms. This provides a simple, easily estimated quantity for quantifying covariance.
Most of our computations are performed at the level of population averages, so we define which is a scalar function of frequency, f , for each a, b = e, i, x and U, Z = E, I, X, S, R, T . It is also convenient to define the 2 × 2 mean-field matrix form, for U , Z = E, I, X, S, R, T . We also define the recurrent and feedforward mean-field connectivity matrices, For the exponential kernels we use, η b (f ) = 1/(1 + 2πif τ b ). The zero-frequency values, w ab = w ab (0) = p ab j ab q b , define time-averaged interactions and mean-field connection matrices, W = W (0) and W x = W x (0). This choice of notation allows us to perform computations on mean-field spike train and input statistics in a mathematically compact way. To demonstrate this, we first review the mean-field analysis of firing rates in the balanced state [27,28,[40][41][42]. Mean external input is given by X = √ N W x r x and mean recurrent input by R = √ N W r so that mean total synaptic input is given by In the balanced state, T , r ∼ O(1), which can only be obtained by a cancellation between external and recurrent synaptic inputs. This cancellation requires W r + W x r x ∼ O(1/ √ N ) so that [27,28,[40][41][42] lim N →∞ in the balanced state. Hence, the balanced state can only be realized when this solution has positive entries, r e , r i > 0, which requires that [27,28,40] X e /X i > w ei /w ii > w ee /w ie . Below, we perform analogous derivations of mean-field CSDs in balanced networks.

III. A REVIEW OF THE ASYNCHRONOUS BALANCED STATE
We first review prior theoretical work that derives meanfield CSDs when spike trains in the external population are uncorrelated Poisson processes (Fig. 1A,B), so Since the derivations of these results [38] and analogous derivations for networks of binary neuron models [29] are presented elsewhere and since the derivations are similar to those presented below, we only review the results here and give the details of the derivation in Appendix A.
Since spike trains in the external population are uncorrelated, correlations between the external input to neurons in the recurrent network arise solely from overlapping feedforward synaptic projections with [29,38] (see Appendix A) where W * x is the conjugate transpose of W x . If the spike trains in the external layer were uncorrelated, but not Poisson, then r x in this equation could be replaced by the power spectral density of the spike trains.
It would at first seem that this O(1) external input correlation would lead to O(1) correlations between neurons' spike trains, but this is prevented by a cancellation between positive and negative sources of input correlation. In particular, correlations between neurons' recurrent synaptic inputs, R, R , are also positive and O(1), but these positive sources of input correlations are canceled by negative correlations between neurons' recurrent and external inputs, X, R , in such a way that the total synaptic input correlation is weak, where R, X = X, R * . This cancellation is realized when the mean-field CSD between spike trains satisfies [38] S, S = 1 where N ×o(1/N ) → 0 as N → ∞ and W − * is the inverse of W * . The first term in this equation represents spike  Raster plot of 200 randomly selected neurons from population X and E respectively. D) Histogram of external (X, green) recurrent (R = E + I, purple) and total (T = X + E + I, black) to all excitatory neurons. Currents here and elsewhere are reported in units CmV /s where Cm is the arbitrary membrane capacitance. E) Mean external (green), recurrent (purple), and total (black) input to excitatory neurons for networks of different sizes, N . F) Mean excitatory (red) and inhibitory (blue) neuron firing rates for difference network sizes. Solid curves are from simulations and dashed curves are from Eq. (4). G) Mean covariance between pairs of excitatory neurons' external inputs (green), recurrent inputs (purple), total inputs (black), and mean covariance between the recurrent input to one excitatory neuron and external input to the other (yellow) for different network sizes. Covariances were computed by integrating the inputs over 250ms windows then computing covariances between the integrals, which is proportional to zero-frequency CSD and has a closer relationship with spike count covariance (see Eq. (3) and surrounding discussion). Integrated currents have units CmmV , so their covariances have units C 2 m mV 2 . H) Zoomed in view of black curve from E on a log-log axis (mean total input covariance, black) plotted alongside the function c/N (dashed gray) where c was chosen so that the two curves match at the largest N value. I) Mean spike count covariance between excitatory neuron spike trains (red), between inhibitory neuron spike trains (blue), and between excitatory-inhibitory pairs of spike trains (purple). Counts were computed over 250ms time windows. Solid curves are from simulations, dashed from Eq. (7) evaluated at zero frequency. Network size was N = 10 4 in B-D.
train correlations inherited from external inputs, namely X, X . The second term, C 0 , represents correlations generated intrinsically by chaotic or chaos-like dynamics in the network [28,38,40,43]. For example, a network with deterministic, constant external input, X a j (t) = X a , generates correlated variability [40] despite the fact that such networks are deterministic once an initial condition is specified and therefore X, X = 0 for such networks. While an exact expression for C 0 is unknown, we show empirically below that its effects are small compared to correlations inherited from external input, at least at low frequencies for the network parameters that we consider. Therefore, the following approximation is accurate To demonstrate these results, we first simulated a network of N = 10 4 randomly and recurrently connected neurons receiving feedforward input from a population of N x = 2000 uncorrelated Poisson-spiking neurons (Fig. 1A,B). As predicted, spiking activity in the recurrent network was asynchronous and irregular ( Fig. 1C; mean spike count correlation between neurons with rates at least 1 Hz was 5.2 × 10 −4 ) with approximate balance between external (X) and recurrent (R) synaptic input sources (Fig. 1D). Varying the network size, N , demonstrates the O( √ N ) growth of mean external (X) and recurrent (R) synaptic input currents that cancel to produce O(1) mean total input current (T ) (Fig. 1E), as predicted by the mean-field theory of balance. As a result, firing rates converge to the limiting values predicted by Eq. (4) (Fig. 1F).
As predicted by the analysis of the asynchronous state, the mean covariances between individual sources of synaptic inputs appear O(1) (Fig. 1G), but cancel to produce much smaller, O(1/N ), total input covariance (Fig. 1G,H). Mean spike count covariances also appear O(1/N ) and show good agreement with the closed form approximation from Eq. (7) (Fig. 1I).
The fact that the approximation in Eq. (7) accurately captures the scaling of spike count covariances from simulations implies that correlated variability inherited by overlapping synaptic inputs dominate intrinsic correlations, C 0 , which are ignored in Eq. (7). Similar results were found in previous work [38]. This contrasts with previous findings in networks of binary neurons, in which intrinsic variability appeared to dominate [44].

IV. THE CORRELATED STATE IN BALANCED NETWORKS
Above, we analyzed correlated variability when spike trains in the external population were uncorrelated, S x , S x = 0, which produced asynchronous spiking in the recurrent network, S, S ∼ O(1/N ). We now relax this assumption by considering moderate correlations between neurons in the external layer ( Fig. 2A), Most previous work analyzing spike train correlations in recurrent networks relies on knowledge of the "correlation susceptibility" or "transfer" function, which quantifies the mapping from synaptic input covariance to spike train covariance, e.g., the mapping from T e , T e to S e , S e . However, susceptibility functions depend sensitively on the neuron model being used and their derivation typically relies on diffusion approximations that assume neurons' input currents can be approximated by Gaussian white noise and are weakly correlated [8,11]. These assumptions are not valid for the densely connected, correlated networks with exponentially-decaying PSC kernels we consider here.
Instead of assuming knowledge of a covariance transfer or susceptibility function, our derivation relies only on an assumption that transfer of mean-field covariance is O(1), specifically that In other words, mean-field spiking statistics are not drastically different in magnitude than mean-field input statistics because neurons transfer inputs to outputs in an O(1) fashion. Importantly, we do not need to know the value of the fraction in Eq. (8).
In addition to this assumption, all of our derivations follow from a few simple arithmetical rules that rely on the bilinearity of the operator ·, · . Specifically (see Appendix A for derivations), for any U , Z = E, I, X, R, S, S x , T where A * is the conjugate-transpose of A and where we omit smaller order terms here and below (see Appendix A for details). Eq. (9) follows from the fact that total input is composed of recurrent and external sources, T = R+X. Eqs. (10) and (11) (12) is simply a property of the Hermitian crossspectral operator.
We first derive the CSD between external inputs to neurons in the recurrent network. Applying of Eqs. (11) and (12) gives Hence, O(1) covariance between the spike trains in the external population induces O(N ) covariance between the external input currents to neurons in the recurrent network. This is a result of the effects of "pooling" on correlations and covariances, namely that the covariance between two sums of N correlated random variables is typically O(N ) times larger than the covariances between the individual summed variables [29,45,46]. We next derive the CSD between spike trains and external inputs. First note that However, it follows from our assumption that neuronal which is only consistent if there is a cancellation between the two terms on the right hand side. Specifically, we must have that since X, X ∼ O(N ) from above where we have omitted terms smaller than O( √ N ). We can now calculate the CSD between spike trains in the recurrent network. First note that However, our assumption of O(1) transfer implies that S, T ∼ S, S so which is only consistent if there is cancellation between the terms on the right hand side. This cancellation can only be realized if which is O(1) and where we have omitted terms of smaller order. Notably, evaluating Eq. (14) does not depend on knowledge of neurons' correlation susceptibility functions or any other neuronal transfer properties, but only depends on synaptic parameters and input statistics. In summary, O(1) covariance between spike trains in the external population produces O(N ) covariance between neurons' external inputs, but O(1) covariance between spike trains in the recurrent network on average. We hereafter refer to this state as the "correlated state" since it produces moderately strong spike train correlations in contrast to the asynchronous state characterized by extremely weak spike train correlations. The reduction from O(N ) external input covariance to O(1) spike train covariance arises from a cancellation mechanism analogous to the one that reduces O(1) external input correlation to O(1/N ) spike train correlations in the asynchronous state (see above).
To demonstrate these results, we simulated a network of N = 10 4 neurons identical to the network from Fig. 1 except that spike trains in the external population were correlated Poisson processes ( Fig. 2A) with Here, r x = 10Hz is the same firing rate used in Fig. 1 and c = 0.1 quantifies the spike count correlation coefficient between the spike trains in the external population over large counting windows. See Appendix B for a description of the algorithm used to generate the spike trains. The recurrent network exhibited moderately correlated spike trains in contrast to spike trains in the asynchronous state (Fig. 2B, compare to Fig. 1C; mean spike count correlation between neurons with rates at least 1 Hz was 0.077) . The mean CSDs between spike trains in the recurrent network closely matched the theoretical predictions from Eq. (14) (Fig. 2C). As in the asynchronous state, external and recurrent synaptic input sources approximately canceled (Fig. 2D), as predicted by balanced network theory.
Varying N demonstrates that the network exhibits the same cancellation between O( √ N ) mean external and recurrent synaptic input sources and that Eq. (4) for the mean firing rates is accurate (Fig. 2E,F). As predicted by the analysis of the correlated state, the covariance between individual sources of input currents appear O(N ) (Fig. 2G), but cancel to produce much smaller, approximately O(1), total input covariance (Fig. 2G,H). Mean spike count covariances also appear O(1) and converge toward the limit predicted by Eq. (14) (Fig. 2I). Hence, despite the complexity of spike timing dynamics in densely connected balanced networks, mean-field spike train covariances is accurately predicted by a simple, linear equation in terms of synaptic parameters. The derivation of this equation does not require the use of linear response theory, which can be problematic for densely connected networks with synaptic kinetics and non-vanishing correlations.

V. THE CORRELATED STATE PRODUCES TIGHT BALANCE BETWEEN EXCITATORY AND INHIBITORY INPUT FLUCTUATIONS CONSISTENT WITH CORTICAL RECORDINGS
We have so far considered cancellation between positive and negative sources of input correlations at the mean-field level, i.e., averaged over pairs of postsynaptic neurons (Figs. 1G,H and 2G,H). In vivo cortical recordings reveal that this cancellation occurs even at the level of single postsynaptic neuron pairs. When one neuron was clamped near its inhibitory reversal potential and another neuron is clamped near its excitatory reversal potential (with spiking suppressed), recorded membrane potential fluctuations are approximately mirror images of one another (Fig. 3A, top). Similarly, if both neurons are held near their excitatory reversal potential (Fig. 3A, middle) or both near their inhibitory reversal potential (Fig. 3A, bottom), recorded membrane potential fluctuations are highly correlated. This implies that fluctuations in the excitatory and inhibitory synaptic input to one neuron are strongly correlated with fluctuations in the excitatory and inhibitory input to other nearby neurons (see [22] for more details and interpretation).
To test whether this phenomenon occurred in our simu-  lations, we randomly chose two neurons and decomposed their synaptic input into the total excitatory (E + X) and the inhibitory (I) components. In the asynchronous state, input current fluctuations were fast and largely unshared between neurons or between current sources in the same neuron (Fig. 3B), in contrast to evidence from in vivo recordings. Input current fluctuations in the correlated state were slower, larger, and most importantly largely synchronized between neurons (Fig. 3C), consistent with evidence from in vivo recordings. This precise tracking of fluctuations in excitatory and inhibitory synaptic currents is referred to as "tight balance" [47] (as opposed to the "loose balance" of the asynchronous state). The results would be similar if we decomposed inputs into their external (X) and recurrent (R = E + I) sources instead of excitatory (E + X) and inhibitory (I).
To better understand this striking difference between input currents in the asynchronous and correlated states, we first computed the average covariance between the excitatory and inhibitory input to pairs of (excitatory) neurons in the network. These averages have the same dependence on network size, N , as they do when input currents are broken into external and recurrent sources (compare Fig. 4A,B to Figs. 1G and 2G). Specifically, in the asynchronous state, covariances between individual current sources are O(1) on average, but cancel to produce weak O(1/N ) covariance between the total synaptic input to neurons on average (Fig. 4A). In the correlated state, the average covariance between individual input sources is O(N ) and cancellation produces O(1) average total input covariance (Fig. 4B).
Hence, despite the precise cancellation of positive and negative sources of input covariance at the mean-field level in the asynchronous state (Fig. 4A), this cancellation is apparently not observed at the level of individual neuron pairs (Fig. 3B). To see why this is the case, we computed the distribution of input current covariances across all pairs of excitatory neurons. We found that these distributions were broad and the distribution of total input covariance was highly overlapping with the distributions of individual input current sources (Fig. 4C, the black distribution overlaps with the others). This implies that cancellation does not reliably occur at the level of individual pairs since, for example, the total input covariance for a pair of neurons can be similar in magnitude or even larger than the covariance between those neurons' excitatory inputs. The distributions of input covariances were strikingly difference in the correlated state. The distribution of total input covariances was far narrower than the distributions of individual input current sources and the distributions were virtually non-overlapping (Fig. 4D). Hence, a precise cancellation between positive and negative sources of input covariance must occur for every neuron pair, leading to the tight balance observed in Fig. 3C.
These results are more precisely understood by computing the variance across neuron pairs of input covariances as N is varied. In the asynchronous state, the variance of input covariances from all sources appear to scale like O(1) (Fig. 4E). Since the mean input covariance between individual sources are also O(1) (Fig. 4A), the overlap between distributions in Fig. 4C is expected. In the correlated state, the variances of input covariances appear to scale like O(N ) except for the variance of the total input covariance, which appears to scale like O(1) (Fig. 4F). Since the variances scale like O(N ), the standard deviations scale like O( √ N ). This, combined with the fact that the mean input covariances between individual sources scale like O(N ), implies that the distributions in Fig. 4E will be non-overlapping when N is large. The same conclusions would be reached if we decomposed inputs into their external (X) and recurrent (R = E + I) sources instead of total excitatory (X + E) and inhibitory (I).

VI. TESTING LINEAR RESPONSE APPROXIMATIONS TO PAIRWISE SPIKE TRAIN COVARIANCE
As discussed above, Eqs. (7) and (14) are appealing because, unlike many previous approximations of spike train CSD and spike count covariance in recurrent networks [8][9][10][11]13], they do not require knowledge of neurons' correlation susceptibility functions or any other neuronal transfer properties. However, they are limited because they only give the population-averaged covariance or CSD.
Previous studies that provide approximations for the full N × N matrix of all spike train CSDs or spike count covariances use linear response approximations that rely on assumptions that the network is sparsely or weakly coupled and/or that external input is uncorrelated or weakly correlated Gaussian white noise [10,11]. These assumptions permit the application of Fokker-Planck and linear response techniques. External input in our models is not weakly correlated and is not approximated well by Gaussian white noise, so these approaches are not fully justified in out model, but might nevertheless yield a useful approximation. We next derive and test a linear response approximation to spike train covariability in our model.
First define to be the full N × 1 vector of spike trains in the recurrent network obtained by concatenating the excitatory and inhibitory spike train vectors. Define the N × 1 synaptic input vectors T (t), R(t), and X(t) similarly. The total input vector can be written as Here, is the N × N recurrent connectivity matrix and Note that the connection weights are time dependent and * denotes the matrix product with multiplication replaced by convolution [11].
The linearity of transfer from S to R and from S x to X implies that and Fourier transform of the matrix J(t) and similarly for J x (f ). This implies, for example, that the N × N matrix of external input CSDs is given by   and similarly for the CSDs between recurrent inputs, R, R = J S, S J * In summary, the transfer of spike train CSDs to input CSDs follows easily from the linearity of transfer from spike trains to inputs [4]. However, a closed expression for the matrix of spike train CSDs in the recurrent network, S, S , requires knowledge of the transfer from total input to spike train CSDs.
If spike trains were related linearly to total input, S = A * T for some diagonal matrix, A(t), then we could derive an exact closed equation for S, S . However, integrateand-fire neuron models transfer their input currents to spike trains nonlinearly. Nevertheless, we can obtain an approximation to spike train CSDs by assuming an approximate linear transfer of CSDs, specifically that [11] S, X ≈ A T , X , where A(f ) is a diagonal matrix with A jj (f ) the "susceptibility function" of neuron j [8,[48][49][50]. Note that these equations would be exactly true if neural transfer were linear, S = A * T . It can be derived from these assumptions that (see Appendix C for derivation and [11] for similar derivations) One problem with testing the approximation in Eq. (18) is that we do not have an estimate of A jj (f ). Spike count covariances over large windows are proportional to zerofrequency CSD (see Eq. (3) and surrounding discussion). Evaluated at zero frequency, a neuron's susceptibility function is equal to its gain, i.e. the derivative of the neuron's firing rate with respect to its mean input [49,50], We therefore tested Eq. (18) for spike count covariances by estimating the gain empirically from simulations. In doing so, we found that Eq. (18) is only moderately accurate at approximating spike count covariances from simulations, both in the asynchronous state (Fig. 5A) and in the correlated state (Fig. 5B). We suspected that some of the error was due to numerical inaccuracies introduced by inverting the large, ill-conditioned matrices in Eq. (18).
To test this hypothesis, we re-wrote Eq. (18) in a mathematically equivalent formulation that does not involve matrix inverses, where Id is the N × N identity matrix. We tested the accuracy of Eq. (19) at zero frequency using the same empirically estimated values of the gains and using the matrix, S, S (0), of spike count covariances estimated from simulations and found that it was very accurate, especially in the correlated state (Fig. 5C,D). Since Eqs. (18) and (19) are mathematically equivalent, this suggests that

much of the observed inaccuracy of Eq. (18) is due to the presence of inverses of large, ill-conditioned matrices.
One shortcoming of Eqs. (18) and (19) is that they require an estimate of A. We obtained this estimate empirically by simulating the entire network and numerically estimating the gains, which greatly limits the utility of the equations. We next wondered whether the accuracy of the equations is sensitive to the accuracy of the gain estimate, or if a rough approximation to the gains would be sufficient. To answer this question, we tested Eq. (19) again, but replaced the individual estimated gains on the diagonal of A with their mean. In other words, we used A(0) = gId where g is the average gain estimated from simulations and Id is the N × N identity matrix. This replacement greatly decreased the accuracy of Eq. (19) (Fig. 5E,F), suggesting that an accurate estimate of the individual gains is necessary for applying and interpreting Eqs. (18) and (19).
In summary, linear response approximations are fairly accurate in the densely connected balanced networks with spatially and temporally correlated noisy feedforward input studied here (Fig. 5C,D). However, their utility is limited by two factors: First, that using linear response theory to approximate spike train CSDs or spike count covariances requires the inversion of large, ill-conditioned matrices, which introduces substantial error (Fig. 5A,B). Secondly, that the application of linear response approximation requires estimates of neurons' individual susceptibility functions or gains. The mean-field equations (7) and (14) do not have these shortcomings, but only give the population-averaged CSDs and covariances. It should also be noted that correlations were weak in our simulations, even in the correlated state (mean spike count correlations between neurons with firing rate above 1 Hz was 0.077). Stronger correlations can be obtained with alternative parameters, which could potentially make Eqs. (18) and (19) less accurate.

VII. CORRELATED VARIABILITY FROM SINGULAR MEAN-FIELD CONNECTIVITY STRUCTURE
We have shown that O(1) spike train correlations can be obtained in balanced networks by including correlations between neurons in an external layer ( S x , S x ∼ O(1)), defining what we refer to as the "correlated state." Previous work shows that O(1) spike train correlations can be obtained in the recurrent network with uncorrelated external spike trains ( S x , S x = 0) when the mean-field connectivity matrix is constructed in such a way that the recurrent network cannot achieve the cancellation required for these states to be realized [37][38][39]. This is most easily achieved using networks with several discrete sub-populations or networks with distance-dependent connectivity. For simplicity, we restrict our analysis to discrete sub-populations. We first extend the mean-field theory from above to such networks, then generalize and extend the analysis from previous work to understand the emergence of O(1) correlations when the mean-field connectivity matrix is singular.
The recurrent networks considered above have two statistically homogeneous sub-populations: one excitatory and one inhibitory and the external population is a single homogeneous population. Suppose instead that there are K sub-populations in the recurrent network, with the kth popualtion containing N k = q k N neurons where k q k = 1. Connectivity is random with p jk denoting the connection probability from population k to j, and j jk / √ N denoting the strengths of the connections for j, k = 1, . . . K. All neurons in population k are assumed to have the PSC kernel η k (t) which is again assumed to have integral 1. Similarly, suppose that the external network contains K x sub-populations each with N x k = q x k N x neurons where q x k = k q x,k = 1. Feedforward connection probabilities and strengths are given by p x jk and j x jk / √ N for j = 1, . . . K and k = 1, . . . , K x . Assume that q k , p jk , j jk , q x k , p x jk , and j x jk are all O(1). We then define the K × K mean-field recurrent connectivity matrix by [W ] jk = p jk j jk q k η k and the mean-field feedforward connectivity matrix by [W x ] jk = p x jk j x jk q x k η x k . For all of the networks considered above, we had K = 2 and K x = 1.
When W is an invertible matrix, Eqs. (4), (7), and (14) are equally valid for networks with several subpopulations as they are for the simpler networks considered above. Hence, the mean-field theory of firing rates and correlations extends naturally to networks with several populations [38][39][40][41][42]51]. However, when W is singular, Eqs. (4), (7), and (14) cannot be evaluated. Instead, Eq. (4) can be re-written as When W is singular, this equation only has a solution, r, when X = −W x r x is in the range or "column space" of W . Otherwise, balance is broken. An in-depth analysis of firing rates in such networks is provided in previous work [41,42,51] (and extended to spatially continuous networks in [40,51]), so we hereafter assume that X is in the range of W and balance is achieved. A similar analysis may be applied to spike train CSDs. For simplicity, we assume here that spike trains in the external population are uncorrelated, S x , S x = 0, since this is the case considered in previous work and since this is the case in which a singular W breaks the asynchronous state. Eq. (7) can be re-written as where we have ignored smaller order terms. When W is singular, Eq. (21) is not guaranteed to have a solution, S, S . More precisely, a solution can only exist when the K × K matrix, X, X , is in the range of the linear operator L defined by In that case, Eq. (21) has a solution so that S, S ∼ O(1/N ) and the asynchronous state is still realized. However, if X, X is not in the range of L, the asynchronous state cannot be realized because Eq. (21) does not have a solution.
Darshan et al. [39] looked at a similar scenario except the singularity of their networks made them incapable of cancelling internally generated covariance, in contrast to the external input covariance considered here. Other work [37,38] analyzed the scenario with external input covariance and singular connectivity, as well as the extension to spatially extended networks. However, this previous work did not completely analyze the structure of correlations, but only showed that the asynchronous state was broken. We next show that correlations in the recurrent network are O(1) and derive their structure.
Let v 1 , v 2 , . . . , v n be a basis for the nullspace of W * and define which is a self-adjoint matrix that defines the orthogonal projection onto the nullspace of W * . Note that P is a Hermitian matrix (P = P * ) and P W = W * P = 0 (the zero matrix). Define the projection A 0 = P AP for any matrix A. Unless X, X is carefully constructed otherwise, we can expect that Then take the projection of both sides of Eq. (22) above to get Hence, the total input CSD is O(1) when X, X is not in the range of L, even though it is X, X /N when W is invertible (i.e., in the asynchronous state). Moreover, the structure of T , T is given to highest order in N by X, X 0 = P X, X P , which can be computed exactly from knowledge of the structure of X, X and W . When neural transfer from T to S is O(1) (see Eq. (8) and surrounding discussion), this implies that S, S ∼ O(1) so that the asynchronous state is broken when X, X is not in the range of L. While we cannot be certain that S, S has the same structure as T , T , it should have a similar structure as long as neural transfer of correlations is similar for each sub-population.
To demonstrate these results, we consider the same network from above with re-wired feedforward projections from the external population. Specifically, divide the excitatory, inhibitory, and external populations each into two equal-sized sub-populations, labeled e 1 , i 1 , x 1 , e 2 , i 2 , and x 2 where population a k contains N a /2 neurons. Hence the network has the same total number of neurons as before, but we have simply sub-divided the populations. To distinguish this network from the one considered in Figs. 1 and 2, we refer to the previous network as the 3-population network and to this modified network as the 6-population network. We re-wire the feedforward connections so that x 1 only connects to e 1 and i 1 , and x 2 only projects to e 2 and i 2 . Specifically, we set the connection probabilities to p aj x k = 2p ax if j = k and p aj x k = 0 if j = k for a, b = e, i and j, k = 1, 2, where p ab are the connection probabilities for the 3-population network and p aj b k for the 6population network. This re-wiring assures that neurons in the recurrent network receive the same number of feedforward connections on average from the external population. The recurrent connectivity structure is not changed at all. Specifically, we set p aj b k = p ab for a, b = e, i. All connection strengths are unchanged, j aj b k = j ab for a = e, i and b = e, i, x and all neurons in the external population have the same firing rate, r x , as before. See Fig. 6A for a schematic of this network.
The feedforward mean-field connectivity matrix can be written in block form as where 0 is the 2 × 1 zero-matrix and W 2×1 x = [w ex w ix ] T is the 2 × 1 feedforward connectivity matrix for the 3population network. Note that W x is 4 × 2 since there are 4 populations in the recurrent network and 2 populations in the external population. The recurrent mean-field connectivity matrix is is the 2 × 2 recurrent connectivity matrix for the 3population network. Note that W is 4 × 4. Here, w ab = p ab j ab q b η b are the same values used above for analyzing the 3-population network. Even though W is non-invertible, X = W x [r x r x ] T is in the range of W for this example, so firing rates in the balanced state can be computed using Eq. (20), and are identical to the firing rates for the 3-population networks considered above.
The nullspace of W * is spanned by the orthonormal vectors so the projection matrix is given in block form by P = 1 2 where I 2 is the 2 × 2 identity matrix.
The external input CSD is determined by the average number of overlapping feedforward projections to any pair of neurons in the recurrent network (multiplied by their connection strength and r x ), which gives (in block form) where 0 is the 2 × 2 zero matrix and X, X 2×2 is the external input CSD from the 3-population network, given by Eq. (5). Therefore, by Eq. (23), In other words, the mean total input CSD between excitatory neurons in the same subgroup (two neurons in e 1 or two neurons in e 2 ; diagonal blocks above) is positive and equal to half the mean external input between the same neurons. Hence, the cancellation by the recurrent network only reduces the external input CSD by a factor of 1/2, as opposed to the O(1/N ) reduction realized in the asynchronous state (when W is invertible). In contrast, the mean total input CSD between excitatory neurons in opposite subgroups (one neuron in e 1 and the other in e 2 ; off-diagonal blocks above) has the same magnitude as for same-subgroup pairs, but is negative. This represents a competitive dynamic between the two groups since they inhibit one another (recurrent connections are net-inhibitory in balanced networks [27,42]), but receive different feedforward input noise. Interestingly, the average CSD between all pairs of spike trains is still O(1/N ) in this example, but it is easy to design examples with singular W in which this is not true. A similar example was considered in previous work [38], but external input was generated artificially instead of coming from an external population.
Simulating this network for varying values of N shows that firing rates approach those predicted by the balance equation (20) (Fig. 6B), confirming that balance is realized. Pairs of excitatory neurons in the same group (both neurons in e 1 or both neurons in e 2 ) receive positively correlated external input and recurrent input (Fig. 6C, purple and green curves) that are partially canceled by negative correlations between their recurrent and excitatory input (Fig. 6C, yellow curve). Because the cancellation is only partial, the correlation between the neurons' total inputs is O(1) (Fig. 6C, black curve) in contrast to the asynchronous state (compare to Fig. 1G,H where cancellation is perfect at large N ). The total input covariance agrees well with the theoretical prediction from Eq. (24) (Fig. 6C, dashed gray line). As a result of this lack of cancellation between total input covariance, spike count covariances are also O(1) and positive between same-group pairs (Fig. 6D). For opposite group pairs (one neuron in e 1 and the other in e 2 ), cancellation is also imperfect, but this leads to negative total input covariance, in agreement with the theoretical prediction from Eq. (24) (Fig. 6E), and leads to negative spike count covariances between neurons in opposite populations (Fig. 6F).
In summary, we have analyzed two mechanisms to generate O(1) spike train correlations in balanced networks. For the first mechanism (Fig. 2), spike trains in the external population are correlated so that external input correlations are O(N ). Cancellation is achieved so that spike train correlations are reduced to O(1). For the other mechanism (Fig. 6), external input correlation is O(1), but precise cancellation cannot be achieved so that spike trains inherit the O(1) correlations from the input. How could these two mechanisms be distinguished in cortical recordings? Under the first mechanism, we showed that fluctuations of inhibitory input to individual neurons closely tracks fluctuations of other neurons' excitatory inputs (Fig. 3C). This should not be the case under the second mechanism because precise cancellation is not realized. Indeed, plotting the excitatory and inhibitory input to three excitatory neurons (two in e 1 and one in e 2 ) shows that input fluctuations are not closely tracked (Fig. 7). This provides a way to distinguish the two mechanisms from paired intracellular recordings. Indeed, the first mechanism (which we refer to as the "correlated state") appears more consistent with the cortical recordings considered here (compare Fig. 3A to Figs. 3C and 7).

VIII. SUMMARY AND DISCUSSION
We analyzed correlated variability in recurrent, balanced networks of integrate-and-fire neurons receiving correlated feedforward input from an external population. We showed that correlations between spike trains in the recurrent network are small (O(1/N )) when spike trains in the external population are uncorrelated, consistent with previous work on the asynchronous state [29,38], but much larger (O(1)) when spike trains in the external population are correlated, giving rise to a "correlated state." In both states, strong correlations in the feedforward input are canceled by recurrent synaptic input due to the excitatory-inhibitory tracking that arises naturally in densely connected balanced networks. This cancellation allows for the derivation of a concise and accurate closed form expression for spike train CSDs in terms of synaptic parameters alone. Hence correlations in balanced networks are determined predominately by synaptic connectivity structure, not neuronal dynamics. The tracking of excitatory synaptic input by inhibition was observable on a pair-by-pair basis in the correlated state, but not the asynchronous state, suggesting that the correlated state is more consistent with in vivo recordings.
We only considered recurrent networks with two, statistically homogeneous neural populations: one excitatory and one inhibitory. Our analysis can be extended to multiple subpopulations as long as each sub-population contains O(N ) neurons, and also extends to networks with connection probabilities that depend on distance, orientation tuning, or other continuous quantities. This analysis has been developed for the asynchronous state in previous work [38] and is easily extended to the correlated state as well. The primary difference is that X, X is Previous work has shown that networks with multiple sub-populations and networks with distance-dependent connectivity can break the asynchronous state in balanced networks when the network connectivity structure is constructed in such a way that the recurrent network cannot achieve the cancellation required for the asynchronous state [37][38][39], leading to O(1) correlations between some cell pairs. We showed that the precise tracking of excitation by inhibition provides an experimentally testable prediction for distinguishing this mechanism from the one underlying the correlated state (see Fig. 7 and surrounding discussion).
Another alternative mechanism for achieving larger correlations in balanced networks is through instabilities of the balanced state. Such instabilities can create patternforming dynamics that produce correlated spiking [40,[52][53][54][55][56][57][58][59]. Future studies should work towards experimentally testable predictions that distinguish correlations that arise from instabilities from those that arise through the mechanisms considered here. For example, since instabilities generate correlations internally, they should produce weak correlations between activity in the recurrent network and activity in the external population(s) providing input to that network [57], in contrast to the mechanisms we consider here. Indeed, some recordings show that local circuit connectivity can increase correlations [60], which is consistent with internally generated correlations, but inconsistent with the mechanisms that we consider here.
In the correlated state, spike train correlations in the recurrent network are essentially inherited from correlations between spike trains in the external population. Hence, the O(1) correlations realized by this mechanism require the presence of another local network with O(1) correlations. This raises the question of where the O(1) correlations are originally generated. One possibility is that they could be generated in a presynaptic cortical area or layer through the alternative mechanisms discussed in the previous paragraph. Another possibility is that they originate from a network that is not in the balanced state at all. Non-balanced networks can easily achieve O(1) spike train correlations simply from overlapping synaptic projections. While cortical circuits are commonly believed to operate in a balanced state, correlations could originate in thalamus, retina, or other sub-cortical neural populations then eventually propagate to cortex.
The cancellation between variances of covariances observed empirically in Fig. 4F is, to our knowledge, a novel observation, but we were unable to derive it analytically. Path integral approaches have recently been applied to compute variances of covariances in linear network models with uncorrelated external input [61] ( X, X diagonal so X, X = 0). This previous work derives a simple relationship between a network's criticality and a parameter, ∆, that represents the ratio between the standard deviation of spike count covariances and the mean spike count variance. Specifically, they showed that in their network model, the eigenvalue spectrum of the network dynamics is given to leading order in N by λ max = 1 − 1/(1 + N ∆ 2 ). In our model, where external input is correlated, the variance of covariances and spike count variances both appear to be O(1), so that ∆ ∼ O(1) and their expression reduces to λ max = 1 − O(1/ √ N ) regardless of the spectral radius of W . We interpret this to mean that their formula is not accurate for networks with correlated external input. Specifically, their formula would erroneously predict nearcritical network dynamics if applied to a network that is far from criticality, but receives correlated external inputs. Future work should consider the possibility of extending their analysis of criticality to networks with external input correlation.
Our mean-field analysis applies to networks of integrateand-fire neuron models, which are arguably more biologically realistic than networks of binary model neurons that are often used for mean-field analysis of neuronal networks. Binary neuron networks are appealing due to the mathematical tractability of their mean-field analysis, but our work demonstrates that integrate-and-fire networks are similarly tractable, calling into question the utility and appeal of binary network models.
Two unproven assumptions underly our mean-field analysis of the correlated state. The first assumption is that neural transfer is O(1) (Eq. (8) and surrounding discussion). The second assumption is that individual connection strengths are not strongly correlated with individual CSD values so that the step from Eq. (A2) to (A3) is valid when ignoring smaller order terms. These assumptions are made in other work, even if not stated explicitly. We have been unable to prove these assumptions rigorously for the model studied here, leaving an open problem for future work.
We showed that linear approximations to spike train covariance developed for small, sparsely coupled net-works [10][11][12] can also be accurate for large, densely connected balanced networks (Fig. 5). However, their usefulness is limited by the need to invert large, ill-conditioned matrices and to approximate the susceptibility functions of individual neurons. The simpler equations we derived for mean-field spike train CSDs (Eqs. (7) and (14)) do not have these problems. Moreover, while linear response approximations require that neural transfer of input is approximately linear, our mean-field derivations did not depend on this assumption. Recent work has called for looking beyond linear analysis of neuronal networks [62]. Our analysis shows that, even in networks where neural transfer of inputs is nonlinear, linear mean-field analysis could still be accurate and useful.
In summary, we showed that correlations in balanced networks can be caused by feedforward input from a population of neurons with correlated spike trains, defining the "correlated state" that is quantitatively captured by a linear mean-field theory. In contrast to other mechanisms of correlation in balanced networks, the correlated state predicts a precise balance between the fluctuations in excitatory and inhibitory synaptic input to individual neuron pairs, consistent with some in vivo recordings [22].
where o( X, X /N ) scales smaller than X, X /N as N → ∞ and where C 0 ∼ O (1). Note that in the correlated state, X, X ∼ O(N ) so that the C 0 /N term can be absorbed into the o( X, X /N ) term. A sketch of the derivation for the correlated state is given in Results, and the derivation is similar in the asynchronous state. Here, we give the details of this derivation.
We first derive Eq. (5) for X, X in the asynchronous state, i.e. when spike trains in the external population are uncorrelated Poisson processes, so S x , S x = 0. In that case, the CSD between two neurons' external input is given by J ax jm J bx kn η x η * x S x m , S x n for a, b = e, i where we have used the fact that the crossspectral operator, ·, · is a Hermitian operator. Since external spike trains are uncorrelated Poisson processes, S x m , S x n = 0 when m = n and S x m , S x m = r x . Therefore, we can rewrite the equation above as From Eq. (2), the expectation of the summand in this equation is Hence, taking expectations across the N x elements of the sum and the coefficient in front of the sum gives X a , X b = q x j ax j bx p ax p bx | η x | 2 r x which, in matrix form, is equivalent to Eq. (5). We next derive Eq. (10). Noting that recurrent input to neuron j in population a = e, i is composed of excitatory and inhibitory components, R a j (t) = E a j (t) + I a j (t), we have where O(avg k,b S e k , U b k / √ N ) accounts for the diagonal terms not counted in the definition of S e , U b . Note that this step requires us to assume that individual values of the random variable, J ae jm , are not strongly correlated with individual values of S e m , U b k , so that the expectation of their product can be replaced by the product of their expectations. This assumption is implicit in derivations in other studies [29,38,61], even though it is never proven and often not made explicit.
Repeating this calculation for I a j , U b k and putting them together gives the average In matrix form, this becomes An identical calculation, replacing R with X and S with S x , gives For the correction terms, O(avg k,b,c S c k , U b k / √ N ) and O(avg k,b S x k , U b k / √ N ), to contribute at largest order in N , it needs to be true that for c = e, i, or x. In the correlated state, this is never the case, so the correction term can be ignored, giving Eqs. (10) and (11). In the asynchronous state, avg k,b S b k , S b k ∼ O(1) since spike train power spectral densities are O(1) due to intrinsically generated variability, but S, S ∼ O(1/N ). This causes the power spectral densities, i.e. auto-correlations, to contribute to correlated variability in the network at the largest order in N and ultimately leads to the presence of the C 0 term in Eq. (6) where (1/N )C 0 represents mean-field CSDs that would be obtained in the absence of external input correlations, X, X = 0 (see other work [38,39,44] for an in-depth treatment of intrinsically generated correlations).
Appendix B: Generation of correlated spike trains for external inputs.
To generate correlated, Poisson spike trains for the external population in the correlated state we used the multiple interaction process (MIP) method [63] with jittering. Specifically, we generated one shared "mother" process with firing rate r m = r x /c. Then, for each of the N x "daughter" processes, we randomly kept each spike in the mother process with probability c. As a result, each daughter process is a Poisson process with firing rate cr m = r x and a proportion of c of the spikes are shared between any two daughter processes. To get rid of perfect synchrony between the daughter processes, we jittered each spike time in each daughter process by a normally distributed random variable with mean zero and standard deviation τ c = 5ms. Upon jittering, the daughter processes remain Poisson and the resulting CSD between daughter processes is given by Eq. (15). Spike count correlations between the daughter processes over large time windows are exactly c. The daughter processes were used as the spike trains, S x j (t) in the external population in the correlated state. See [63] for a deeper analysis of this algorithm.
V re = −75mV, V lb = −100mV, ∆ T = 1mV, and V T = −55mV. Synaptic currents in figures are reported in units C m V /s. Covariances between synaptic currents are computed between integrals of the currents (see Eq. (3) and surrounding discussion), so the covariances have units C 2 m mV 2 .
All simulations and numerical computations were performed on a MacBook Pro running OS X 10.9.5 with a 2.3 GHz Intel Core i7 processor. All simulations were written in Matlab (Matlab R 2018a, MathWorks). The differential equations defining the neuron model were solved using a forward Euler method with time step 0.1ms. Statistics in Figs. 1D, 2C,D, and 4C,D were computed from a simulation of duration 50s. Statistics in Figs. 1E-I, 2E-I, and 4A,B,E,F were computed by repeating a simulation of duration 50s over ten trials for each value of N , then averaging over trials. For each trial, network connectivity was generated with a different random seed, so the statistics are averaged over time and over realizations of the "quenched" variability of network connectivity. Statistics in Fig. 5 were computed by repeating a simulation of duration 100s for 50 trials, with network connectivity the same for each trial. Statistics were then averaged over trials. Gains were estimated by fitting a rectified quadratic function to the relationship between all neurons' firing rates and mean total inputs (r j and T j ), then computing the derivative of the fitted quadratic at the input value for each neuron. The same approach was used in previous work [38,51] to estimate a mean-field gain.