When do correlations increase with firing rates in recurrent networks?

A central question in neuroscience is to understand how noisy firing patterns are used to transmit information. Because neural spiking is noisy, spiking patterns are often quantified via pairwise correlations, or the probability that two cells will spike coincidentally, above and beyond their baseline firing rate. One observation frequently made in experiments, is that correlations can increase systematically with firing rate. Theoretical studies have determined that stimulus-dependent correlations that increase with firing rate can have beneficial effects on information coding; however, we still have an incomplete understanding of what circuit mechanisms do, or do not, produce this correlation-firing rate relationship. Here, we studied the relationship between pairwise correlations and firing rates in recurrently coupled excitatory-inhibitory spiking networks with conductance-based synapses. We found that with stronger excitatory coupling, a positive relationship emerged between pairwise correlations and firing rates. To explain these findings, we used linear response theory to predict the full correlation matrix and to decompose correlations in terms of graph motifs. We then used this decomposition to explain why covariation of correlations with firing rate—a relationship previously explained in feedforward networks driven by correlated input—emerges in some recurrent networks but not in others. Furthermore, when correlations covary with firing rate, this relationship is reflected in low-rank structure in the correlation matrix.


Introduction
One prominent goal of modern theoretical neuroscience is to understand how the features of cortical neural networks lead to modulation of spiking statistics [1][2][3].This understanding is essential to the larger question of how sensory information is encoded and transmitted, because such statistics are known to impact population coding [4][5][6][7][8].Both experimental and theoretical inquiries are complicated by the fact that neurons are widely known to have heterogeneous attributes [9][10][11][12][13][14].
One family of statistics that is implicated in nearly all population coding studies is trial-to-trial variability (and co-variability) in spike counts; there is now a rich history of studying how these statistics arise, and how they effect coding of stimuli [15][16][17][18][19]. Recent work by numerous authors has demonstrated that the information content of spiking neural activity depends on spike count correlations and its relationship (if any) with stimulus tuning [15,17,[19][20][21].Since a population of sensory neurons might change their firing rates in different ways to stimuli, uncovering the general mechanisms for when spiking correlations increases with firing rate (or when they do not) is important in the context of neural coding.Thus, we study this question in a general recurrent neural network model.
One observation that has been made in some, but not all, experimental studies is that pairwise correlations increase with firing rates.This relationship has been observed in vitro [22] and in several visual areas: area MT [23], V4 [24], V1 [25,26], and notably, in ON-OFF directionally sensitive retinal ganglion cells [21,27].The retinal studies involved cells with a clearly identified function, and therefore allowed study of the coding consequences of this correlation/firing rate relationship.Both studies found that the stimulus-dependent correlation structure observed compared favorably to a structure in which stimulus-independent correlations were matched to their (stimulus-)averaged levels.This finding reflects a general principle articulated in other studies [17,19], that stimulus-dependent correlations are beneficial when they serve to spread the neural response in a direction orthogonal to the signal space.
While many studies have illustrated the connection between stimulus-dependent correlation structure and coding, these have (until recently: see [21,25,27]) largely taken the correlation structure as given, leaving open the question of how exactly a network might produce the hypothesized correlation structure [6,7] (see also the theoretical calculations in [21,27]).Theoretical studies of the mechanisms that contribute to correlation distributions have largely analyzed homogeneous networks (i.e. cells are identical, aside from E/I identity) [2,3,28,29], which does not allow an exploration of a correlation/firing rate relationship.Thus, how correlation coefficients can vary across a population of heterogeneously-tuned neurons is not yet well understood despite its possible implications for coding.
In this paper we investigated the relationship between correlations and firing rates in conductance-based leaky integrate-and-fire (LIF) neural network models, consisting of excitatory (E) and inhibitory (I) cells that are recurrently and randomly coupled.We introduced neural heterogeneity by allowing thresholds to vary across the population, which induced a wide range of firing rates, and induced different firing regimes by varying the strength of recurrent excitation.We found that with relatively strong excitation, pairwise correlations increased with firing rate.
In theoretical studies, this correlation-firing rate trend has been explained in feed-forward networks driven by common input [22,30,31].Here we investigated whether the correlation/firing relationship in recurrent networks can be explained by this theory, but where the source of input correlations is internally generated; i.e., from overlapping projections within the recurrent network.We first adapted a network linear response theory, to decompose predicted correlations into contributions from different graph motifs, which are subgraphs which form the building blocks of complex networks [28,32,33].We found distinct patterns in how motifs contribute to pairwise correlation, between the different spiking regimes: with weak excitation, negative third-order motifs partially cancel positive second-order motifs (as in Pernice et al. [28]), thus diluting a possible correlation/firing rate relationship, whereas third-order motifs reinforce second-order motifs with stronger excitation.
However, in both regimes second-order motifs -and specifically inhibitory common input -were still the dominant contributor to overall pairwise correlations.This allowed us to generalize theory from [22], and describe pairwise correlations in terms of a single-cell susceptibility function.Surprisingly, we found that correlations from inhibitory common input could either increase or decrease with firing rate, depending on how cells responded to fluctuations in inhibitory conductances.
We further show that a correlation-firing rate relationship has an important consequence for heterogeneous networks; it can shape low-dimensional structure in the correlation matrix.Low-dimensional structureoften modeled with a low-rank approximation to the correlation matrix -is important because it can be used to improve estimation [34] and even to reconstruct full correlation matrices from incomplete data [35][36][37]; such structure has been observed in experimental data [25,[38][39][40][41] but its origin is not always known.We demonstrate in our networks that when correlation co-varies with firing rate, the (E-E) correlation matrix could be accurately modeled with a low-rank approximation, and the low-rank projection in this approximation was strongly associated with firing rate.Thus we demonstrate that low-rank structure can result recurrent activity modulated by single-cell characteristics, as well as from a global input or a top-down signal [38].

Results
We studied asynchronous recurrent networks, where we varied the strength of excitation to get different firing behaviors.We find that the covariation of correlations with firing rates -a phenomenon observed in feed-forward networks -occurs here in one firing regime, but not the other.This could be explained in terms of how single cells responded to fluctuations in inhibitory conductance.Finally, we show that when correlations covary with firing rates, the correlation matrix admits a low-rank approximation.

Asynchronous firing in heterogeneous networks
We performed Monte Carlo simulations of recurrent, randomly connected E/I networks, as described in Methods: Neuron model and network setup.To connect to previous literature on asynchronous spiking, we compared networks with and without single-cell variability -referred to as heterogeneous and homogeneous respectively.Heterogeneity was introduced by allowing cell threshold to vary, which induced a corresponding range of firing rates (see Methods: Neuron model and network setup for details).We first chose parameters so that the networks exhibited the classical asynchronous irregular (Asyn) regime, in which each neuron has irregular Poisson-like spiking, correlations are low, and the population power spectra are flat [42].In Fig 1A we show raster plots from both the heterogeneous and homogeneous networks, in this regime.The heterogeneous network shows a gradient in its raster plot, because cells are ordered by decreasing firing rate.The population power spectra were flat, for both E and I cells and in both homogeneous and heterogeneous networks (Fig 1C).
When we increased excitation (by increasing both W EE and W IE , where W XY is the conductance strength from type Y to X; see Table 1 for parameter values), we observed occasional bursts of activity.However, the bursts do not occur at regular intervals and do not involve the entire population (we found excitatory bursts involved at most 25% of the population).The network is still moderately inhibition-dominated and neurons are spiking irregularly; example raster plots are shown in Fig 1B .The population power spectra (Fig 1D) are no longer flat (compare to the asynchronous regime, Fig 1C ); they show local maxima around 8 Hz, but it is not a pronounced peak.We will refer to this as the strong asynchronous (SA) regime [43].
In both Fig 1C and Fig 1D, we note that -despite the apparent differences in the distribution of spikes across the network, evident in the raster plots -both the autocorrelation functions (Fig 1C,D, insets) and the power spectra from the heterogeneous and homogeneous networks are very similar.Thus, we have a fair comparison to examine the role of heterogeneity, independent of other characteristics of the network.
The distribution of both excitatory and inhibitory firing rates are extremely narrow in the homogeneous network, but broad in the heterogeneous network (Fig 1E).This is expected, as each excitatory (inhibitory) cell in the homogenous network has the same uncoupled firing rate; because the number of synaptic inputs is likewise fixed, population variability in synaptic input is limited.Therefore the heterogeneous networks have a range of firing rates, which allows us to investigate the possibility of a relationship between (variable) firing rate and pairwise correlations.Population-averaged firing rates were very similar between the heterogeneous and homogeneous networks: in the asynchronous regime ⟪ν E ⟫ = 10.6 Hz (heterogeneous) and ⟪ν E ⟫ = 10.1 Hz (homogeneous), while ⟪ν I ⟫ = 44.3Hz (heterogeneous) and ⟪ν I ⟫ = 43.5 Hz (homogeneous).In both regimes Fano factors ranged between 0.9 and 1.1, consistent with Poisson-like spiking (more statistics are given in Tables S1 and S3).

Correlation increases with firing rate in the strong asynchronous regime
We next sought a possible relationship between pairwise correlations -quantified via the Pearson's correlation coefficient for spike counts, ρ ij ≡ Cov T (n i , n j ) Var T (n i ) Var T (n j ) -and single-cell firing rates.Such relationships have been found in feed-forward networks [22,30,31], and impact information transfer when considered in concert with stimulus selectivity (i.e.signal correlations) [7,8,15,19].In heterogeneous networks, the large range of firing rates -equivalently the large range of operating points -admits the possibility that cells at different operating points may differ in their ability to transfer correlations.
To investigate this we plotted pairwise correlations for each distinct excitatory pair ρ ij , versus the geometric mean of the firing rates √ ν i ν j , in both regimes (asynchronous and strong asynchronous), for a range of time scales (blue stars in Fig 2).We focus here on excitatory-excitatory (E-E) pairs, because excitatory synaptic connections provide the predominant means of propagating cortical sensory information to higher layers.Our results show a striking difference between the two spiking regimes; while there is no clear relationship with firing rate in the asynchronous regime (Fig In recurrent networks, the response of each cell is shaped by both direct and indirect connections through the network.We used the linear response theory described in Methods: Linear Response Theory and Methods: Computing statistics from linear response theory to predict the full correlation matrix C T at various time scales, including the limit of long time scales: C(0) = lim T →∞ theory successfully captured E-E correlations, both the full distribution of values and coefficients of individual cell pairs (details, including figures, can be found in: Supporting Information: S1 Text).
We then plotted the predicted correlation, Cij Cii Cjj , vs. geometric mean firing rate √ ν i ν j (magenta circles in Fig 2).The predicted correlations captured the same positive relationship observed in Monte Carlo results, with R 2 values of 0.47, 0.4, and 0.36.

Decomposition of correlation by graph motifs shows a distinct pattern for each spiking regime
Why does a correlation/firing rate relationship emerge in one spiking regime, but not the other?In feedforward networks, a positive correlation/firing rate relationship results from transferring common input through fluctuation-driven, asynchronously-firing cells [22,30].In contrast, the amount of shared input into two cells in a recurrent network is determined by both direct and indirect connections through the network.To separate the impact of different network pathways, we decomposed the linear response-predicted correlations at long time scales (i.e.C(0) = lim T →∞ 1 T C T ) into normalized contributions from n-th order motifs, as described in Methods: Quantifying the role of motifs in networks.Common input from a divergent connection, for example, results from the 2nd-order motif K * C 0 K.In Fig 3 In the asynchronous regime (top panel of Fig 3A ), first-order contributions ( R1 ) separate into three distinct "curves", reflecting a 1-1 relationship with firing rate conditioned on first-order connectivity (no connection between i and j; one connection between i and j; bidirectional connection between i and j).This might be expected if susceptibility to excitatory conductances, Ãg E (0), varies roughly like firing rate (as we will argue later); thus, the contribution from an i → j connection should approximately vary like: Second-order contributions are overall positive while third-order contributions are overall negative (consistent with [28]); neither appear to have a relationship with firing rate.Second-order contributions are conspicuously dominant; fifth and sixth order terms are near zero.
It is instructive to compare the results from the homogeneous network, shown in the bottom panel.Since firing rate is nearly uniform across the network, we would not expect to see a relationship; rather, we might use this to gauge the relative magnitude of contributions at different orders, independent of the complicating effect of firing rate.Indeed, we see similar pattens as in the heterogeneous network; R1 ij breaks into three distinct curves and is positive; R2 ij is generally positive and R3 ij generally negative; the magnitude of R2 ij overshadows both R1 ij and R3 ij .This qualitative picture changes when we consider the strong asynchronous regime, in Fig 3B .First-order contributions follow a similar pattern as in the asynchronous regime, and second-order contributions are likewise positive.However, third-order contributions are positive, and in the heterogeneous network they have a distinctly positive relationship with firing rate (top panel).Thus, in the asynchronous regime, negative third-order contributions partially cancel with positive second-order contributions; in the strong asynchronous regime, first, second, and third-order motifs reinforce each other, contributing to an overall positive relationship with firing rate (black dots).Despite these differences, second-order contributions appear to dominate both regimes.Therefore, we next analyze contributions from specific second-order motifs in Fig 4 .There are four distinct second-order motifs that can correlate two E cells.There are two types of chains, from K 2 C 0 and C 0 (K * ) 2 .An E → E → E chain tends to positively correlate; an E → I → E chain will negatively correlate; these are shown as blue and green respectively.There are two types of common input, from KC 0 (K * ); they correspond to common input from E and I cells, i.e.E ← E → E and E ← I → E. They both lead to positive correlations and are shown as red and magenta respectively.In the asynchronous regime (left panel of Fig 4A ), the dominant contributions are I common input (magenta) and negative (E → I → E) chains (green); correlating chains (blue) and excitatory common input  (red) are barely visible, as they are clustered near zero.This contrasts with the strong asynchronous case (right panel): blue and red dots are now visible, similar in magnitude, and show a clear 1-1 trend with firing rate.As a result, they prevent "dilution" of correlation from the decorrelating chains (green dots).In both regimes, inhibitory common input appears to be the dominant second-order motif.In Fig 4B we plot the contribution from different second-order motifs vs. the total contribution from second-order motifs, R2 ij (rather than geometric mean firing rate, √ ν i ν j ).In both panels, the inhibitory common input (magenta) clusters around the unity line, showing it is the best predictor of the total second-order contribution.
In conclusion, decomposition of pairwise correlations into graph motifs has shown us two important things: first, while third-order motifs contribute to the positive correlation/firing rate relationship observed in the SA regime, second-order motifs still dominate in both regimes.Second, inhibitory common input is the most important second-order motif in both regimes (Fig 4B).

Susceptibility to inhibition can either increase or decrease with firing rate
In feedforward networks -i.e. in the absence of a path between two cells -correlations in outputs (i.e.spike trains) must arise from correlations in inputs; for example, through shared or common inputs.We have found that inhibitory common input is the dominant contributor to pairwise correlations in both the asynchronous and strong asynchronous regimes; we now turn our attention to modeling this term (inhibitory common input) specifically.
Previous work that analyzed the relationship between the long-time correlation and firing rate in feedforward networks [22,30] quantified a susceptibility function that measures the ratio between output and input correlations: If both cells receive a large (but equal) number of uncorrelated inputs, c would be the fraction of inputs that are common to both i and j.
In the networks examined here, each cell had a fixed in-degree for both excitatory and inhibitory cells; however, for any given pair of cells i and j, the number of E and I inputs that synapsed onto both cells will vary from pair to pair.Thus, we next considered the possibility that our (negative) finding in the asynchronous network could be explained by accounting for variable c ij .
We focus on inhibitory common input, which is the dominant second-order contribution in the asynchronous network (Fig 4).We segregated pairs by whether they had 0, 1, 2, etc.. common inhibitory inputs; we then use this number as a proxy for c (recall that each excitatory cell had exactly 7 inhibitory inputs, so that this number divided by 7 approximates the common input fraction; two common inputs imply c ≈ 0.28 for example).We plot the results for the asynchronous network in Fig 5A, top panel (data for each distinct value of c is presented by color).As we might expect, correlation increases as c increases.However, for a fixed c, there is not an apparent relationship between firing rate and correlation; if anything, there appears to be a slight decrease.Correlation also increases with c in the strong asynchronous network (Fig 5A, bottom panel); however, here we also see a modest increase with geometric mean firing rate √ ν i ν j .
Previous theoretical work [22,30] identified an increase in susceptibility with firing rates in current-driven neurons; we next considered the possibility that this fails to hold for conductance-driven neurons.As described in Methods: Quantifying correlation susceptibility, we estimated correlation susceptibility for each pair of neurons, by using the susceptibility function for each neuron to conductance fluctuations (computed as part of the linear response theory), divided by a measure of the long-timescale spike count variance: We plotted the results for both networks in Fig 5B ; while susceptibility increases with firing rate in the strong asynchronous network (except for the largest firing rates), it actually decreases with firing rate in the asynchronous network.We can contrast with the estimated susceptibility to current fluctuations (i.e.A µ,i , with µ i , τ eff,i , and σ eff,i as in Eq 29) which we also computed for the same set of cell pairs, shown in Fig 5C.
Here, we see that S µ ij increases with firing rates, in both networks.To understand the difference between S g I ij and S µ ij , we consider how the cell responds differently to currents vs. conductances.Both networks are inhibition-dominated, with resting potentials in the range (−0.33, −0.28) (asynchronous) and (−0.183, −0.075) (strong asynchronous); since the cells are operating near the inhibitory reversal potential (E I = −0.5),fluctuations in the inhibitory conductance have limited effect on the firing rate of the cell.This effect is more pronounced in the asynchronous regime, where the small distance between the resting and synaptic potentials creates a saturating effect on the firing rate; previous authors found susceptibility decreases with firing rate, in a simple thresholding model in which the f-I curve saturates [22].In contrast, susceptibility to excitatory conductances, S g E ij , more closely resembled the current input susceptibility S µ ij .The cells are operating far from the excitatory reversal potential (E E = 6.5); therefore fluctuations in g E are multiplied by a relatively constant E E − V and thus have a "current-like" effect.
However, the difference between a threshold-linear vs. saturating f-I curve cannot be a perfect analogy, because for our neurons firing rate and susceptibility are not functions of a single parameter µ: instead, they can depend on all six parameters defining the cell; specifically, we define the single-cell susceptibility where This quantity is shown in Fig 6, where it is plotted vs. firing rate ν i (blue stars).Note that this is a negative quantity; since the susceptibility for a neuron pair S is the product (and therefore positive), an increase in S ⟨g I ⟩ i will result in a decrease in S ⟨g I ⟩ ij and vice versa.We have also used the asynchronous spiking assumption, that Cii (0) ≈ ν i .
By examining the dependence of susceptibility on its parameters, we hypothesized that it was most important to capture dependence on mean inhibitory conductance and threshold (see S1 Text: Approximating single-cell susceptibility in a heterogeneous network).We approximated S ⟨g I ⟩ i , by reevaluating the firing rate function in which σ I,i , ⟨g E,i ⟩, σ E,i and σ i have been replaced by their average values: i.e.
where In both the asynchronous and strong asynchronous regimes, the red stars form a scattered cloud around the average value ⟨ ⟨g I,i ⟩ ⟩ p , with no obvious relationship with θ i .This fact motivated a further simplification, i.e., we replaced ⟨g I,i ⟩ with its population average, ⟨ ⟨g .Here, we can see clearly that in the asynchronous regime, correlations should actually decrease with firing rate, for ν i > 5 Hz.In the strong asynchronous regime, correlations will increase with firing rate, saturating around 10-15 Hz.
Finally, recall that our actual network sampled a relatively small part of the (θ, ⟨g I ⟩) plane.This may be attributed to the fact that we generated firing rate diversity (and therefore heterogeneity), by modulating cell excitability through the cell threshold θ i .How might our results have changed, if we had generated firing rate diversity through some other mechanism?In both regimes, we can increase firing rates by either decreasing ⟨g I,i ⟩, or by decreasing θ (see Fig. S7).To explore this, we computed susceptibility values along another curve in the (θ, ⟨g I ⟩) plane; specifically, we held θ fixed and instead varied ⟨g I ⟩ (illustrated with black squares on Fig 7); i.e. where

Low-rank structure in neural correlations is mediated by correlation-firing rate relationship
Previous work has identified low-dimensional structure in neural correlation matrices [25,[38][39][40][41]; its origin is not always known [3].We next hypothesized that the positive correlation-firing rate relationship we observed in the strong asynchronous regime, might be reflected in low-dimensional structure in the correlation matrix.For simplicity, suppose that correlations were really represented by a function of firing rate (as in [22]): i.e. ρ ij = cS(ν i )S(ν j ).Then we could represent the off-diagonal part of the correlation matrix as C T = c S S T , where S is a length N vector such that S i = S(ν i ); that is, C T would be a rank-one matrix.We followed the procedure outlined in Methods: Low-rank approximation to the correlation matrix to approximate each correlation matrix, C T , as the sum of a diagonal matrix and low-rank matrix: where λ is given in closed form by the eigenvalues of C T : and σ 1 , u 1 are the first singular value and singular vector of C T .
In Fig 8, we show the results from heterogeneous networks in both the asynchronous (top panel in each subfigure) and strong asynchronous (bottom panel in each subfigure) regimes.We first show C T − λI, where   In the asynchronous network, the approximated correlations take on a narrow range (between 0 and 0.01, compared to between −0.015 and 0.03 for the measured coefficients) and do not show an obvious positive relationship.In the strong asynchronous regime, the range is more accurate ( between 0.02 and 0.1, vs. 0.01 and 0.15 for the measured coefficients) and the points cluster around the unity line.
In Fig 8D , we plot the weight of each cell in the first singular vector, (u 1 ) j vs. the firing rate ν j .We can clearly see a positive relationship in the strong asynchronous regime (bottom panel), suggesting that the positive relationship between correlation and firing rate is related to the success of the low-rank approximation.

Discussion
We simulated heterogeneous, asynchronous networks in order to investigate a possible relationship between firing rates and pairwise correlations in recurrent networks.We found that correlations can either increase or decrease with firing rates; this could be attributed to differences in how cells responded to fluctuations in inhibitory conductances.When correlations did increase with firing rates, this relationship was reflected in low-dimensional structure in the correlation matrix.This study offers an example of a practical consequence of the difference between treating synaptic inputs as conductances rather than currents; while most synaptic currents are more accurately modeled as conductances, current-based formulations are often used for analytical and computational simplicity.Although it is known that neural models responding to currents vs. conductances differ in their response dynamics [44][45][46], this approach is supported by findings that steady-state firing rates are qualitatively similar in both settings (e.g.[47]).Here, we found that refined features of the steady-state firing rate surface will govern susceptibility to common input in asynchronous networks; two "cuts" through this surface may yield divergent behavior with respect to correlation susceptibility, despite yielding similar firing rates.
In other words, the relationship between pairwise correlations and firing rate will depend on the means through which firing rate diversity is achieved.In our study, we created firing rate diversity by regulating cell excitability; if we had instead varied mean inhibitory input (by varying the number of inhibitory connections) or a background excitatory current (which would model diversity in stimulus tuning from feed-forward inputs), we would likely have seen a different pattern.Finally, the recurrent network will also shape the path the cells follow through the "firing rate surface"; to generalize [22] to recurrent networks, we need to identify both how firing rates are produced and how they are shaped by the recurrent network.
Low-dimensional structure has been a common finding in many large-scale neural recordings [25,[38][39][40][41]; while the origin is not always known, it is often interpreted as arising from a global input or top-down signal.This is an interpretation that arises naturally from the technique of factor analysis, in which one seeks to explain a data vector as the sum of a random vector and the linear combination of some number of latent factors [48] (for Gaussian random variables, each latent factor can literally be interpreted as a global input with a distinct pattern of projection onto the observed variables).In our network, we found that a single latent factor was effective at capturing correlations in the strong asynchronous regime; however, this latent factor did not reflect common input (there was no global external input into the network) but rather modulation from single-cell characteristics.Thus, we identify a novel mechanism that can contribute to low-dimensional structure in neural recordings.

Impact on stimulus coding
The networks studied here were not encoding a stimulus; correlations were generated by recurrent activity, given that each neuron had a baseline firing rate in the absence of recurrent input.However, we can readily connect this network to a stimulus coding task, in order to understand how the correlation-firing rate relationship can impact coding.
Consider a population of cells that is responsible for encoding a single scalar stimulus θ, such as movement direction or orientation of a visual stimulus, and that each cell has roughly a bell-shaped tuning curve.Furthermore, we model incoming stimulus by modulating a stimulus-dependent background current I i, (θ); i.e., cells which prefer the current stimulus have a higher level of current, and thus a higher firing rate, than cells which prefer an orthogonal or opposite stimulus.The network we studied here would model the response to a single stimulus θ 0 ; that is, the firing rate diversity we observe is present because some cells are strongly tuned to the current stimulus, while others are not.
We could extend this model, by resetting background currents to model a complete set of stimuli {θ 1 , θ 2 , ...θ n−1 }.For each stimulus θ j , correlations would show the rough firing rate dependence displayed in the strong asynchronous network, resulting in a stimulus-dependent correlation structure in which pairwise correlations vary like geometric mean firing rate.This is the structure analyzed in [21,27]: the authors found that such a stimulus-dependent correlation code enhances information, when compared to a stimulusindependent code with the same average correlation level.Intuitively, the mean population response lives on the surface of a (hyper-)sphere in neural response space; the population encodes location on this surface.Positive correlations between similarly tuned cells produce response distributions that are stretched in the radial direction, "orthogonal" to this sphere, and thus have a minimal impact on the encoded variable.
Moreover, the mechanism that produced stimulus-dependent correlations in [21,27] was similar to that shown here (see also [25]); common input modulated by stimulus-dependent gain factors.Here, we demonstrated how these stimulus-dependent gain factors might arise (or not) in a recurrent network.If excitation is tuned to put the network in the strong asynchronous regime, then the (stimulus-dependent) correlation structure that results will be favorable to coding.If excitation is tuned to put the network in the asynchronous regime, then correlations are overall low and not stimulus-dependent (although, given that average correlations are not matched, we do not here compare information contained within the two networks).

Future work
This work has, necessarily, focused only on a subset of network attributes that might affect firing statistics.One important feature is the frequency of higher-order graph motifs; experiments have shown that specific motifs will occur more frequently, than would be expected in an Erdős-Rényi network with fixed single-cell connection probability [49].Theoretical work has found that in networks of integrate-and-fire neurons, an overabundance of divergent and chain motifs will lead to enhanced correlation [33] (this finding does depend on the dynamical regime; different motifs impact correlations in networks of coupled oscillators [32]).In [33], the authors use the assumption of homogeneous single-cell characteristics to find parsimonious and instructive formulae for the average correlation, and give a roadmap for how this might be generalized to heterogeneous networks.We look forward to considering the combined effect of single-cell and network heterogeneity in future work.
Another source of cell-to-cell heterogeneity is how cells respond to stimuli, as emphasized in the previous discussion [17,20,21,27,50] (see [19] for a review).Here, we did not consider a specific sensory system with tuning but rather focus on the general question of how the distribution of correlation values arise in recurrent networks.Given the previous discussion, one next step will be to investigate how correlations covary with firing rates, when cell-to-cell heterogeneity is produced by stimulus tuning in a structured network responding to a single variable (such as direction or orientation).
Finally, for numerical tractability our simulations here were performed in relatively small networks.While high average correlations have been measured in experiments [51], theoretical models of asynchronous networks have found that correlations must go to zero as the system becomes large (N → ∞) [2].However, recent work has found that this does not have to be true, as long as spatial structure is introduced into the network [52].We anticipate that this may carry over to other forms of heterogeneity, such as single-cell variability, and that therefore the effect we observe here persists for larger networks.We look forward to reporting on this in future work.

Neuron model and network setup
We considered randomly connected networks of excitatory and inhibitory neurons.Each cell was a linear integrate-and-fire model with second-order alpha-conductances, i.e. membrane voltage v i was modeled with a stochastic differential equation, as long as it remained beneath a threshold θ i : When v i reaches θ i , it is reset to 0 following a refractory period: Each neuron was driven by a Gaussian, white background noise, with magnitude σ i depending only on the cell type; that is, ⟨ξ i (t)⟩ = 0 and ⟨ξ i (t)ξ i (t + s)⟩ = δ(s).The membrane time constant, τ m , and excitatory and inhibitory synaptic reversal potentials, E E and E I , are the same for every cell in the network.Each cell responded to synaptic input through conductance terms, g E,i and g I,i , which are each governed by a pair of differential equations: τ r,X dg where Y = {E, I} denotes the type of cell i and X = {E, I} denotes the type of the source neuron j.Each spike is modeled as a delta-function that impacts the auxiliary variable g X,i ; here t j,k is the k-th spike of Synaptic reversal potential 6.5 -0.5 cell j.The rise and decay time constants τ r,X and τ d,X and pulse amplitude α X depend only on the type of the source neuron; i.e. they are otherwise the same across the population.The parameter W Y X denotes the strength of X → Y synaptic connections, which are (once given the type of source and target neurons) identical across the population.The "raw" synaptic weight (listed in Table 1) is divided by N Y X , the total number of X → Y connections received by each Y -type cell.
We chose connections to be homogeneous and relatively dense, consistent with the local architecture of cortex.Connection probabilities ranged from 20 %-40 %, consistent with experimentally measured values [53][54][55].For our baseline network state, we then chose synaptic weights so the network is moderately inhibition-dominated (α E W IE < α I W II and α E W EE < α I W EI ); that is both E and I cells receive more inhibition than excitation) and shows noisy spiking consistent with the classical asynchronous state.Each neuron receives a fixed number of incoming connections, the identities of which are chosen randomly.(The specific cell ID numbers differ in the different simulations shown below.)For most of the networks we discuss here N = 100 with the 80/20 ratio typical of cortex (i.e n E = 80, n I = 20).Each excitatory cell receives N EE = 32 (40 %) excitatory and N EI = 7 (35 %) inhibitory connections; each inhibitory cell receives N IE = 16 (20 %) and N II = 8 (40 %) inhibitory projections.
In heterogeneous networks, the threshold θ i varied across the population.For both excitatory and inhibitory neurons, the thresholds θ i were chosen from a log-normal distribution between 0.7 and 1.4 (where the rest potential, V r = 0).To be precise, log θ i was chosen from a (truncated) normal distribution with mean −s 2 θ 2 and standard deviation s θ .With this choice, θ i has mean 1 and variance: e s 2 θ − 1.Thus we can view s θ as a measure of the level of threshold heterogeneity.
Throughout this paper, we set s θ = 0.2, which results in a wide range of firing rates compared to the homogeneous case.This was the only source of cell-to-cell heterogeneity; all other parameters were identical across the population, conditioned on neuron type.In homogeneous networks, the threshold was the same across the population: θ i = 1.
Monte Carlo simulations were performed using the stochastic forward-Euler method (Euler-Maruyama), with a time step much smaller than any time scale in the system (∆t = 0.01 ms).Each network was simulated for one second of simulation time, after an equilibration time.Then, a large number of realizations of this interval (n R = 10 5 ) were simulated.Spike counts were retained in each 1 ms window (for a total of 1000 windows) within a realization.With this large number of realizations/trials, the error bars on the resulting time-dependent firing rates were small.Therefore we emphasize that the firing rate pattern is largely driven by network connectivity; while firing is driven by random fluctuations in the background noise, any cell-to-cell variability in the trial-averaged firing rates are not an artifact of the finite number of trials.

Linear Response Theory
In general, computing the response of even a single neuron to an input requires solving a complicated, nonlinear stochastic process.However, it often happens that the presence of background noise linearizes the response of the neuron, so that we can describe this response as a perturbation from a background state.This response is furthermore linear in the perturbing input and thus referred to as linear response theory [56].The approach can be generalized to yield the dominant terms in the coupled network response, as well; we will use the theory to predict the covariance matrix of activity.
We first consider the case of a single cell: an LIF neuron responding to a mean zero current X i (t) (otherwise, the mean of X i can simply be absorbed into E i ).For a fixed input X i (t), the output spike train y i (t) will be slightly different for each realization of the noise ξ i (t) and initial condition v i (0).Therefore we try to work with the time-dependent firing rate, ν i (t) ≡ ⟨y i (t)⟩, which is obtained by averaging over all realizations and initial conditions.Linear response theory proposes the ansatz that the firing rate can be described as a perturbation from a baseline rate proportional to the input X i : ν i,0 is the baseline rate (when X = 0) and A i (t) is a susceptibility function that characterizes this firing rate response up to order [22,29,57].We now consider the theory for networks; here cell i responds to the spike train of cell j, y j (t), via the synaptic weight matrix W, after convolution with a synaptic filter F j (t): In order to consider joint statistics, we need the trial-by-trial response of the cell.We first propose to approximate the response of each neuron as: that is, each input X i has been replaced by the synaptic input, and J ij = W ij F j (t) includes both the i ← j synaptic weight W ij and synaptic kernel F j (normalized to have area 1); A i (t) is the susceptibility function from Eq 17.In the frequency domain this becomes where ỹi = F [y i − ν i ] is the Fourier transform of the mean-shifted process (ν i is the average firing rate of cell i) and f = F [f ] for all other quantities.In matrix form, this yields a self-consistent equation for ỹ in terms of ỹ0 : where Kij (ω) = Ãi (ω) Jij (ω) is the interaction matrix, in the frequency domain.The cross-spectrum is then computed To implement this calculation, we first solve for a self-consistent set of firing rates: that is, ν i is the average firing rate of where We must then compute the power spectrum ⟨ỹ 0 (ω)ỹ 0 * (ω)⟩ and the susceptibility A i (ω), which is the (first order in ) response in the firing rate r i (t) = r 0 i + A i (ω) exp(ıωt) in response to an input current perturbation X(t) = exp(ıωt) (here ı is used for √ −1, while i denotes an index).Both can be expressed as the solution to (different) first-order boundary value problems and solved via Richardson's threshold integration method [47,58].
In our simulations, we used conductance-based neurons; this requires modification, compared with the simpler current-based models.We first approximate each conductance-based neuron as an effective currentbased neuron with reduced time constant, following the discussion in [59].First, separate each conductance into mean and fluctuating parts; e.g.g E,i → ⟨g E,i ⟩ + (g E,i − ⟨g E,i ⟩).Then we identify an effective conductance g 0,i and potential µ i , and treat the fluctuating part of the conductances as noise, i.e. g E,i −⟨g E,i ⟩ → σ E,i ξ E,i (t): where We next simplify the noise terms by writing and assume that the fluctuating part of the voltage, v i − µ i , is mean-zero and uncorrelated with the noise terms ξ E,i (t) [59].That allows us to define an effective equation where and the fluctuating voltage, v i (t) − µ i , now makes no contribution to the effective noise variance.Finally, we consider how to model the conductance mean and variance, e.g.⟨g E,i ⟩ and σ 2 E,i .In our simulations, we used second order α-functions: each conductance g X,i is modeled by two equations that take the form τ r,X dg where X = E, I and the summation is over all type-X spikes incoming to cell i. (For notation purposes, αX,i includes all factors that contribute to the pulse size in Eq 16, including synapse strength and pulse amplitude.)The time constants τ r,X , τ d,X may depend on synapse type; the spike jumps αX,i may depend on synapse type and target cell identity.We assume that each spike train is Poisson, with a constant firing rate: i.e. each spike train is modeled as a stochastic process S(t) with Then a straightforward but lengthy calculation shows that where ν X,i is the total rate of type-X spikes incoming to cell i.
We now describe how these considerations modify the linear response calculation.First, for the selfconsistent firing rate calculation, Eq 22 is replaced by an equation with a modified time constant, conductance, and noise (Eq 29).
We next compute the susceptibility in response to parameters associated with the conductance, i.e. ⟨g E,i ⟩ and σ 2 E,i .This differs from the current-based case in two ways: first, there is voltage-dependence in the diffusion terms, which results in a different Fokker-Planck equation (and thus a different boundary value problem to be solved for the power spectrum ⟨ỹ 0 (ω)ỹ 0 * (ω)⟩).Second, modulating the rate of an incoming spike train will impact both the mean and variance of the input to the effective equation, Eq 23 (via µ i and σ X,i ).Furthermore, this impact may differ for excitatory and inhibitory neurons, giving us a total of four parameters that can be varied in the effective equation.However, neither consideration presents any essential difficulty [47].
Therefore we apply Richardson's threshold integration method directly to Eq 23: When we compute susceptibilities, the parameter to be varied is either a mean conductance -⟨g Thus we have a total of four susceptibility functions Ã⟨g E ⟩,i (ω), Ã⟨g I ⟩,i (ω), Ãσ 2 E ,i (ω), and Ãσ 2 I ,i (ω).Since the Fokker-Planck equation to be solved is linear, we can compute both susceptibilities separately and then add their effects.We now have the interaction matrix: where L(ω) plays a similar role as J, but for the effect of incoming spikes on the variance of conductance.Its relationship to J (either in the frequency or time domain) is given by the same simple scaling shown in Eq 35: i.e., for j excitatory, where the first factor comes from the effective spike amplitude αE,i (and is the scale factor proposed in [47], Eqn.(64)), and the second arises from using second-order (vs.first-order) alpha-functions.We use a modified version of the implementation given by [29] for Richardson's threshold integration algorithm [47,58] to compute rate ν i , power ⟨ỹ 0 i (ω)ỹ 0 * i (ω)⟩, and the various susceptibilities ( Ã⟨g E ⟩,i (ω), Ã⟨g I ⟩,i (ω), Ãσ 2 E ,i (ω), and Ãσ 2 I ,i (ω)) for an LIF neuron.We validated our code using exact formulas known for the LIF [60], and qualitative results from the literature [61].

Computing statistics from linear response theory
Linear response theory yields the cross spectrum of the spike train, ⟨ỹ i (ω)ỹ * j (ω)⟩, for each distinct pair of neurons i and j (see Eq 21).To recover a representative set of statistics, we rely on several standard formulae relating this function to other statistical quantities.
The cross correlation function, C ij (τ ), measures the similarity between two processes at time lag τ , while the cross spectrum measures the similarity between two processes at frequency ω: The Weiner-Khinchin theorem [56] implies that {C ij , Cij } are a Fourier transform pair: that is, In principle, the crosscorrelation C(t) and cross-spectrum C(ω) matrices are functions on the real line, reflecting the fact that correlation can be measured at different time scales.In particular, for a stationary point process the covariance of spike counts over a window of length T , n i and n j , can be related to the crosscorrelation function C ij by the following formula [4]: The variance of spike counts over a time window of length T , n i , is likewise given by integrating the autocorrelation function C ii : It can be helpful to normalize by the time window, i.e.
Cov T (n i , n j ) we can now see that for an integrable cross correlation function (and bearing in mind that the cross-spectrum is the Fourier transform of the cross correlation), that lim Thus, we can use Cij (0) and C ij (0) as measures of long and short time correlations respectively.Finally, the Pearson's correlation coefficient of the spike count defined as: is a common normalized measure of noise correlation, with ρ ∈ [−1, 1].While Cov T and Var T grow linearly with T (for a Poisson process, for example), ρ T,ij in general will not (although it may increase with T ).In general, ρ T,ij depends on the time window T ; however for readability we will often suppress the T -dependence in the notation (and use ρ ij instead).
Here, X µ (t) is a (possibly) time-dependent change in a parameter, such as input current or mean inhibitory conductance; y i,0 is the baseline spike train (when X = 0).A µ,i (t) is a susceptibility function that characterizes the cell's response (to the parameter variation) as long as X µ (t) is small [22,29,57].Following [22], the cross-spectrum of y can now be approximated as: where Cµ (ω) is the spectrum of the parameter variation.The susceptibility has an appealing interpretation in the limit ω → 0, as the derivative of the classical f-I curve: where ν i is the steady-state firing rate of cell i, assuming we can measure it for specific values of the parameter µ. lim ≈ Ãµ,i (0) Ãµ,j (0) This motivates the definition of a correlation susceptibility, which approximates the change in pairwise correlation induced by a parameter change experienced by both cells i and j: If this increases with firing rate -that is, if dν > 0 -then corrections will also increase with firing rate.We can further analyze this quantity by making an assumption for asynchronous spiking, that spike count variance is equal to spike count mean; i.e.Var T (n which motivates the definition of the single-cell quantity In general, the firing rate depends on all single cell parameters included in Eqn. ; i.e. there exists some function f such that (recall that the susceptibility for ω = 0 is the derivative of the firing rate with respect to the appropriate parameter (here, mean inhibitory conductance ⟨g I ⟩).We first consider the firing rate, shown in Fig. S1A.In the heterogeneous network, both excitatory and inhibitory firing rates have large ranges that span approximately an order of magnitude.The linear response theory accurately captures all aspects of the firing rate distributions.The inhibitory firing rates are higher than excitatory, consistent with this population receiving a stronger excitatory input (vs. the excitatory population; compare W IE N E and W EE N E , from Table 1 in the main text).In the homogeneous network, firing rates are strongly clustered around their mean values (shown in Table S1).They are also well-predicted, although the linear response theory does appear to slightly overestimate the inhibitory rates (see Table S2).Similarly, spike count variances match well (shown for long time windows (T = 100 ms) in Fig. S1B). of excitatory cells in a particular time window (Eqn.(40), main text).As in the Monte Carlo simulations, we have assumed spike count statistics to be stationary over time, so that for each T , spike counts n i and n j are treated as random variables sampled both over realizations (i.e trials) and time t.We computed these statistics for both short (T = 5 ms) and long (T = 100 ms) time windows, and illustrate them in Figure S1C; statistics from both heterogeneous (left panels) and homogenous (right panels) are show.E-E correlations are weakly positive, with a small fraction of pairs (∼ 5%) having values below zero.In all panels, the mean/median of the distribution are captured well by the linear response theory; however, the linear response calculation appears to slightly underestimate the simulated variance, as evidenced by the "taller and thinner" distribution shown in solid red (each histogram is computed by distributing 80 × 79 2 distinct coefficients over equally sized bins).The correlation values computed by linear response have comparable ranges in the heterogeneous and homogenous networks, similar to MC simulations and in contrast to first-order statistics.Linear response theory also predicts the distribution of spike count covariances (i.e. the numerator of Eqn. ( 40)): we show these in Fig. S1D.As for correlations, theory underestimates the observed variance of the distributions.However, it appears to capture the fat right tails in the heterogeneous network very well (Fig. S1D, top row).
We now turn our attention to the strong asynchronous (SA) regime, in which both types of excitatory connections were strengthened (see Table 1); the resulting network shows occasional, irregular bursts of concentrated activity (see Fig. 1B).Many of the overall trends are similar to the asynchronous case; we focus on the differences.
Excitatory firing rates were slightly under-predicted by linear response theory (Fig. S2A; see Table S3 for homogeneous rates).Similarly, spike count variances (Fig. S2B) were under predicted.Spike count correlations ρ EE are now positive, with few or no negative correlations (Fig. S2C).The mean is significantly under-predicted; the predicted distributions appear slightly narrower than the observed (Monte Carlo) distribution.Spike count covariances for long time windows are shown in Fig. S2D; the linear response theory appears to capture the qualitative shape of the distributions, particularly the fat right tail in the heterogeneous network (top panel).However, as for correlations (Fig. S2C), the mean is under-predicted in both networks (Table S4).

Linear response theory predicts the first-and second-order statistics of individual cells
We next investigate how well these statistics are predicted on a cell-to-cell basis.This is crucially important when individual correlation coefficients ρ ij within a simulation may vary over an order of magnitude or even in sign.For example, consider the heterogeneous network illustrated in Fig. 2C(bottom): E-E correlations were weakly positive (on average less than 0.01) but could range as high as 0.03 or as negative as -0.015 for some cell pairs.If I pick a specific cell pair i, j out of the population, can I predict where in this range ρ ij will fall?Predicting the correlation of specific cell pairs would be a valuable tool, as many models for heterogeneity -such as those based on population density methods (e.g.[?] ) -do not capture cell-to-cell variation.
We find that single-cell statistics are very accurately predicted for the heterogeneous network.In Figure S3, we show firing rate (Fig. S3A) and Fano factor (Fig. S3B) for three time windows: T = 5, 50, 100 ms.In each panel, both quantities from the Monte Carlo simulations and linear response theory are plotted, on a cell-by-cell basis.In Fig. S3A, the red stars give the firing rate of the uncoupled neurons (i.e.determined only by the threshold θ i and the level of background noise).The effect of coupling is to lower the firing rate of the E cells but to raise the firing rate of the I cells; this is captured very well by the fixed point iteration of Eqn.(22).There is still significant heterogeneity in the firing rates due to variable threshold, with high threshold neurons maintaining comparatively lower firing rates and low threshold neurons maintaining comparatively higher firing rates.
We now analyze the ability of linear response to predict two-cell statistics.In Fig. S3C we plot the spike count correlation ρ ij = Cov T (n i , n j ) Var[n i ] Var[n j ], for all possible E-E pairs in the heterogeneous network, at both T = 5 ms (top) and 100 ms (bottom).The values predicted by linear response theory matches well with the Monte Carlo simulations in both overall range and cell-to-cell; in both plots, the points cluster around the unity line.
We now consider how well linear response models the homogeneous network on a cell-to-cell basis.As in the heterogeneous network, single-cell statistics are accurately predicted (because both simulated and predicted single-cell statistics are nearly constant across the population, we report their values in Table

S2
). Firing rate is slightly overestimated, as is variance.Fano factor differs systemically with time interval: cell activities appear slightly "sub-Poisson" for T = 5 ms, but "super-Poisson" for T = 50, 100 ms.We then examined two-cell statistics: E-E correlations were weak and positive, and clustered in a cloud around the unity line (Fig. S3D), for both short (T =5 ms, top) and long (T =100 ms, bottom) time windows.
Although we mostly focus on E-E correlations here, we observed excellent results in predicting other two-cell statistics, for example excitatory-inhibitory (E-I) correlations.In Fig. S3E we show E-I correlations for T = 100 ms, for both the heterogeneous (top) and homogeneous (bottom) networks.E-I correlations took on a wider range of values; both positive and negative, with a range between [−0.15, 0.15] for T = 100 ms.In the homogeneous network they cluster in four distinct clouds (Fig. S3E, bottom): on closer inspection, these correspond to the presence or absence of direct connections between the pairs.For E-I pairs with no direct connection, correlations are weak and positive.Pairs with only a E → I connection are strongly positively correlated, while pairs with only an I → E connection are strongly negatively correlated.Pairs with BOTH connections are weakly negatively correlated, which may reflect the fact that W IE > W EI .We also find good results when we move to the strongly asynchronous case.This network has increased excitation (W EE = 9 and W EI = 8, vs. W EE = 5 and W EI = 0.5 in the asynchronous regime) and shows short bursts of activity (see Fig. 1); since this violates the assumption of constant firing rate, a priori we cannot be sure linear response theory will be successful.However, the theory is nonetheless successful at matching broad trends in firing rate, Fano factors, and cell-pair correlations (Fig. S4).There are differences between the simulations and linear response calculations.For excitatory neurons, firing rate is slightly overestimated (Fig. S4A), variance underestimated and Fano factor overestimated (Fig. S4B).For inhibitory neurons, firing rate appears to be very accurate; variance and Fano factor are slightly overestimated.We also see that ρ EE is systematically underestimated (Fig. S4C, heterogeneous; Fig. S4D, homogeneous;); ρ IE may also be slightly underestimated, but less so (see Fig. S4E).
Finally, the cell-by-cell accuracy of the linear response theory is reflected in the overall structure of the

Approximating single-cell susceptibility in a heterogeneous network
To understand how the single-cell susceptibility (Eqn.(4), main text) depends on the six parameters ⟨g I,i ⟩, etc., we plotted each parameter vs. firing rate (see Fig. S6 and Fig. S7 for the asynchronous and strong asynchronous regimes respectively).The panels in Fig. S6A show ⟨g I,i ⟩, σ I,i , ⟨g E,i ⟩, and σ E,i , which appear to be randomly scattered with no relationship to firing rate (there is also no apparent relationship in the derived relationships for effective time constant τ eff,i , effective potential µ i , and effective noise σ eff,i , Fig. S6B); σ i is constant for all cells.However, for both networks there is a clear relationship with θ i .Furthermore, of the four non-constant parameters with no discernible relationship, the values of ⟨g I ⟩ appear to have the greatest spread; we therefore hypothesized that we can approximate S ⟨g I ⟩ i , by reevaluating the firing rate function in which σ I,i , ⟨g E,i ⟩, σ E,i and σ i have been replaced by their average values: i.e. where and ⟨ ⋅ ⟩ p denotes the population average.To explore the role of how the cause of firing rate diversity might regulate correlations, we next plot firing rate as a function of threshold θ and mean inhibitory conductance ⟨g I ⟩; that is, we plot F (⟨g I ⟩, θ) (see Eqn. (74)).In both regimes, firing rate decreases from left to right and bottom to top (Fig. S8).Any curve transversal to the level curves of the firing rate will sample a wide range of firing rates, but possibly different susceptibilities.The points corresponding to the actual excitatory cells in our network are illustrated in red; black squares illustrate an alternate curve, where ⟨g I ⟩ is varied but θ = 1.

Fig 1 .Fig 2 .
Fig 1.Two firing regimes in heterogeneous networks.Monte Carlo simulations illustrating two firing regimes we consider in this paper.(A) Raster plots from the asynchronous (Asyn) regime.(B) Raster plots from the strong asynchronous (SA) regime, showing occasional bursts of activity.(C) Power spectra in the asynchronous regime.(D) Power spectra in the strong asynchronous regime.(E) Firing rates in the asynchronous (top panel) and SA (bottom panel) regimes.In (A-B), cells are ordered by increasing threshold value.Power spectra (C-D) are normalized to their maximum value and expressed in decibels/Hz.

Fig 3 .
Fig 3. Pairwise correlations are built from graph motifs.Contributions of different orders to prediction of E-E correlations in the asynchronous (A) and strong asynchronous (B) regimes with linear-response theory.Both heterogeneous (top panels) and homogenous (bottom panels) networks are shown.(A)First order contributions (blue dots) have a 1-1 correspondence (separated into three curves for three possible first order connectivity structures).Second order contributions (red) are positive, third order (green) are negative leading to cancellation; fourth (magenta) order contributions are negative while fifth and sixth are miniscule.(B) First and second order contributions are similar to (A) but with a stronger dependance on firing rate.Third order contributions (green) are distinctly positive so that there is not cancellation as in (A).See main text for further discussion.Bottom panels of (A,B) show results from homogeneous networks.The magnitude and range of each term is similar as in the corresponding heterogeneous network (top panel), but with less variation in firing rate.

Fig 4 .
Fig 4. Inhibitory common input is the dominant second-order motif in both asynchronous and strong asynchronous networks (A) Contributions of different 2nd-order motifs to prediction of E-E correlations in a heterogeneous network, in the asynchronous (left) and strong asynchronous (right) regimes.Left: inhibitory common input (magenta) is dominant but is partially cancelled by E → I → E chains (green); common E inputs (red) and E → E → E chains (blue) are near 0. Right: In contrast to the asynchronous regime, common E inputs (red) and E → E → E chains (blue) are positive and vary with firing rate, preventing dilution of correlation.(B) As in (A), but plotted vs. total contribution from second-order motifs R2 .The left panel shows that inhibitory common input (magenta) is balanced by E → I → E chains (green), because it is above the unity line.The right panel shows that the inhibitory common input (magenta) is instead reinforced by other motifs, because it is below the unity line.

Fig 5 .
Fig 5. Susceptibility to conductance fluctuations can explain correlation-firing rate relationships.In (A-C): heterogeneous asynchronous (top) and heterogeneous strong asynchronous (bottom).(A) Correlation (ρ) from I common inputs vs. firing rate, segregated based on the number of common inhibitory inputs.(B) Estimated correlation susceptibility to fluctuations in inhibitory conductances vs. firing rate (S g I ij ).(C) Correlation susceptibility to fluctuations in inhibitory currents vs. firing rate (S µ ij ).
) and ⟨ ⋅ ⟩ p denotes the population average.The results are also illustrated in Fig 6 (red triangles).In the asynchronous regime (Fig 6A), the results are remarkably close to the original quantities, indicating that using average parameter values has little effect; in the strong asynchronous regime (Fig 6B) the difference is larger, but the points appear to occupy the same "cloud".However, we can now visualize the susceptibility as a function of only two parameters, and we do so in Fig 7 by evaluating Ŝ⟨g I ⟩ i on a (θ, ⟨g I ⟩) grid; the points corresponding to the actual excitatory cells in our network are illustrated in red.
θ) Results are shown in Fig 6 (purple diamonds) and show a strikingly different relationship with firing rate; in the asynchronous regime, correlations should increase with firing rate for ν < 15 Hz; in the strong asynchronous regime correlations will increase with firing rate, saturating near 20 Hz.

Fig 6 .
Fig 6.How firing rate diversity is achieved in a heterogeneous network will affect susceptibility.Single-cell susceptibility function(s) for a conductance-based LIF neuron, as a function of firing rate ν.Successive approximations shown are: original single-cell susceptibility, S ⟨g I ⟩ i (Eq 4, blue stars); most parameters set to average value, Ŝ⟨g I ⟩ i (Eq 73, red triangles); all parameters but θ i set to average value, Ŝ⟨g I ⟩ i (Eq 9, gold squares); and θ fixed, Ŝ⟨g I ⟩ θ=1 (Eq 10, purple diamonds).(A) Asynchronous regime.(B) Strong asynchronous regime.

Fig 7 .
Fig 7. Susceptibility as a function of inhibitory conductance and threshold.Single-cell susceptibility function for a conductance-based LIF neuron, as a function of mean inhibitory conductance ⟨g I ⟩ and threshold θ: Ŝ⟨g I ⟩ (⟨g I ⟩, θ) (defined in Eq 73).Other parameters are set to the population average.Overlays show (⟨g I,i ⟩, θ i ) values of the actual cells in the network (red stars) and an alternative curve through the plane, (⟨g I ⟩, 1), along which comparable firing rate diversity can be observed (black squares).(A) Asynchronous regime.(B) Strong asynchronous regime.

Fig 8 .
Fig 8. Low-rank structure in correlation matrices.Approximating correlation matrices for the heterogeneous networks as a diagonal plus rank-one.In each column of (A-D), the asynchronous (top) and strong asynchronous (bottom) regimes are shown; T = 100 ms.(A)The shifted E-E correlation matrix, C T − λI, for an appropriately chosen λ. (B) A rank-one approximation to C T − λI.(C) True correlation coefficients vs. rank-one approximation, cell-by-cell.(D) Weight in the first singular vector, u 1 vs. geometric mean firing rate √ ν i ν j .

Fig
Fig S1.Theory predicts population statistics in the asynchronous regime.Distributions of spiking statistics, comparing the results of linear response theory with Monte Carlo simulations in the asynchronous regime.(A) Firing rates ν i , for the heterogeneous network.(B) Normalized spike count variances, Var T [n i ] T , heterogeneous network, T = 100 ms.(C) Spike count correlations, for 5 ms and 100 ms time windows: heterogeneous (left two panels) and homogeneous (right two panels).(D) Spike count covariance, 100 ms time windows: heterogeneous (top) and homogeneous (bottom).

Fig
Fig S2.Theory predicts population statistics in the strong asynchronous regime.Distributions of spiking statistics, comparing the results of linear response theory with Monte Carlo simulations in the strong asynchronous regime.(A) Firing rates ν i , for the heterogeneous network.(B) Normalized spike count variances, Var T [n i ] T , heterogeneous network, T = 100 ms.(C) Spike count correlations, for 5 ms and 100 ms time windows: heterogeneous (left two panels) and homogeneous (right two panels).(D) Spike count covariance, 100 ms time windows: heterogeneous (top) and homogeneous (bottom).

Fig
Fig S3.Theory predicts cell-by-cell statistics in the asynchronous regime.Distributions of spiking statistics, comparing the results of linear response theory with Monte Carlo simulations in the asynchronous regime, cell-by-cell.(A) Firing rates ν i , for the heterogeneous network.(B) Fano factor, heterogeneous network.Time windows (T ) shown are: 5 ms, 50 ms, 100 ms.(C) Spike count correlations, heterogeneous network.Time windows are: 5 ms (top) and 100 ms (bottom).(D) Spike count correlations, homogenous network.Time windows are: 5 ms (top) and 100 ms (bottom).(E) E-I spike count correlations, T = 100: heterogenous (top) and homogeneous (bottom).

Fig
Fig S4.Theory predicts cell-by-cell statistics in the strong asynchronous regime.Comparing the results of linear response theory with Monte Carlo simulations in the strong asynchronous regime, cell-by-cell.(A) Firing rates ν i , for the heterogeneous network.(B) Fano factor, heterogeneous network.Time windows (T ) shown are: 5 ms, 50 ms, 100 ms.(C) Spike count correlations, heterogeneous network.Time windows are: 5 ms (top) and 100 ms (bottom).(D) Spike count correlations, homogenous network.Time windows are: 5 ms (top) and 100 ms (bottom).(E) E-I spike count correlations, T = 100: heterogenous (top) and homogeneous (bottom).

Fig
Fig S5.Theory captures low-rank structure in correlation matrices.Approximating correlation matrices (for the heterogeneous networks) obtained from linear response theory, as a diagonal plus rank-one.In each column of (A-D), the asynchronous (top) and strong asynchronous (bottom) regimes are shown; T = 100 ms.(A) The shifted E-E correlation matrix, C T − λI, for an appropriately chosen λ. (B) A rank-one approximation to C T − λI.(C) True correlation coefficients vs. rank-one approximation, cell-by-cell.(D) Weight in the first singular vector, u 1 vs.geometric mean firing rate √ ν i ν j .

Fig
Fig S7.Effective parameters in the heterogeneous network: strong asynchronous regime Parameters used to estimate susceptibility, for all excitatory neurons in the network.Each parameter is plotted vs. firing rate.(A) Mean excitatory conductance ⟨g E,i ⟩ (top left), mean inhibitory conductance ⟨g I,i ⟩ (top right), excitatory conductance variability σ E,i (bottom left), and inhibitory conductance variability σ I,i (bottom right).(B) Effective time constant τ eff,i (top left), effective input current µ i (top right), effective current noise variability σ eff,i (bottom left), and threshold θ i (bottom right).

Fig S8 .
Fig S8.Firing rate as a function of inhibitory conductance and threshold.Firing rate of a conductance-based LIF neuron, as a function of mean inhibitory conductance ⟨g I ⟩ and threshold θ: Ŝ⟨g I ⟩ (⟨g I ⟩, θ) (defined in Eqn (19)).Other parameters are set to the population average.Overlays show (⟨g I,i ⟩, θ i ) values of the actual cells in the network (red stars) and an alternative curve through the plane, (⟨g I ⟩, 1), along which comparable firing rate diversity can be observed (black squares).(A) Asynchronous regime.(B) Strong asynchronous regime.
2, top row), the strong asynchronous regime shows a distinct positive trend with firing rate (Fig 2, bottom row).We can quantify a hypothesized relationship between ν and ρ with linear regression, and indeed find that geometric mean firing rate explains a substantial part of the variability of correlations in the strong asynchronous regime obtained from the Monte Carlo simulations, with R 2 values (i.e. percentage of variability explained) of 0.41, 0.37, and 0.34 for time windows of T = 5, 50, and 100 ms respectively (in contrast, R 2 values for the asynchronous network are below 0.005).
I,i ⟩ ⟩ p .The results are shown in Fig 6 (gold squares) and (as we should expect) allow us to discern a one-to-one relationship with firing rate ν i ; importantly, it appears to capture the average behavior of the actual susceptibility values S

Table 1 .
Excitatory connection strengths mediate between different firing regimes.Parameter W EI (I → E) W IE (E → I) W EE W II σ i (i ∈ E) σ i (i ∈ I)XY denotes X → Y connections; e.g.W EI denotes the strength of excitatory connections onto inhibitory neurons.The parameter σ i denotes the strength of background noise in units of (scaled) voltage, and depends only on cell type (E or I).

Table S1 .
Statistics from heterogeneous vs. homogeneous networks: asynchronous regimeComparing Monte Carlo simulations with predictions from linear response; firing statistics in the strong asynchronous regime.Statistics displayed here are: firing rates for both excitatory and inhibitory populations; Fano factor (FF) for both excitatory and inhibitory populations; spike count correlations for excitatory pairs only (ρ EE ).Standard deviations are reported across the population; i.e. across eighty (80) E cells, or twenty (20) I cells, or 3160 E-E pairs.