Investigating the Correlation–Firing Rate Relationship in Heterogeneous Recurrent Networks

The structure of spiking activity in cortical networks has important implications for how the brain ultimately codes sensory signals. However, our understanding of how network and intrinsic cellular mechanisms affect spiking is still incomplete. In particular, whether cell pairs in a neural network show a positive (or no) relationship between pairwise spike count correlation and average firing rate is generally unknown. This relationship is important because it has been observed experimentally in some sensory systems, and it can enhance information in a common population code. Here we extend our prior work in developing mathematical tools to succinctly characterize the correlation and firing rate relationship in heterogeneous coupled networks. We find that very modest changes in how heterogeneous networks occupy parameter space can dramatically alter the correlation–firing rate relationship.


Introduction
One prominent goal of theoretical neuroscience is to understand how spiking statistics of cortical networks are modulated by network attributes [9,28,42]. This understanding is essential to the larger question of how sensory information is encoded and transmitted, because the statistics of neural activity impact population coding [7,[15][16][17]37]. 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 [1,10,18,25,33]. Recent work has demonstrated that the information content of spiking neural activity depends on spike count correlations and its relationship (if any) with stimulus tuning [1,6,18,25,44].
An important relationship observed in many experimental studies is that pairwise correlations on average increase with firing rates. This has been observed in vitro [8] and in several visual areas: area MT [2], V4 [5] (when measured between cells in the same attentional state), V1 [19,36], and notably, in ON-OFF directionally sensitive retinal ganglion cells [11,44]. The retinal studies involved cells with a clearly identified function, and therefore allowed investigation of the coding consequences of the observed correlation-firing rate relationship. These studies found that the stimulusdependent 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 [18,25], that stimulus-dependent correlations are beneficial when they serve to spread the neural response in a direction orthogonal to the signal space.
These findings thus provide strong motivation for understanding what network mechanisms can produce this positive (and perhaps beneficial) correlation-firing rate relationship. This correlation-firing rate trend has been explained theoretically in feedforward networks driven by common input [8,26,38]; however, many cortical networks are known to be dominated by strong recurrent activity [24,34,39]. On the other hand, theoretical studies of the mechanisms for correlations in recurrent networks have largely analyzed homogeneous networks (i.e., identical cells, aside from excitatory/inhibitory cell identity) [9,13,27,28,40,41], and have not considered how correlations vary with firing rates with realistic intrinsic heterogeneity. Thus, how spike count correlations can vary across a population of heterogeneouslytuned, recurrently connected neurons is not yet well understood despite its possible implications for coding.
In a previous paper, we addressed this gap by developing a mathematical method to compactly describe how second-order spike count statistics vary with both intrinsic and network attributes [4]. Specifically, we adapted network linear response theory [14,27,43] to account for heterogeneous and recurrent networks, which in turn allows us to identify important network connections that contribute to correlations via a single-cell susceptibility function [8]. Here, we will use this method to survey a broad family of recurrent networks to understand how three factors influence the relationship between correlations and firing rates; how the neurons occupy parameter space, the strength of recurrent excitation, and the strength of background noise. This work thus provides a more complete picture of how even modest changes in important circuit parameters can alter the correlation-firing rate relationship.

Results
First, we summarize how a network linear response theory allows us to use the singleneuron firing rate function to approximate correlations. We then apply this theory to examine three factors that can modulate the correlation-firing rate relationship: direction in effective parameter space, strength of recurrent excitation, and strength of background noise.

Using the Single-Neuron Firing Rate Function to Characterize Correlation Susceptibility
The response of a neuron is generally a nonlinear function of model parameters. However, background noise linearizes this response, so that we can use a linear theory to describe the change in response to small changes in parameters. We assume we have some way to approximate the change in firing rate which occurs as a result of a change in parameter X: ν i,0 is the baseline rate (when X = 0) and A X,i (t) is a susceptibility function that characterizes the firing rate response [8,20,41]. The parameter which is modulated has often been chosen to be a current bias μ [8,41]; however, it could also be the mean or variance of a conductance, a time constant, or a spike generation parameter [29,30]. In a coupled network, the parameter change X i arises from inter-neuron coupling; substituting X i (t) → j (J X,ij * ν j ) and moving to the spectral domain, we find whereỹ i = F[y i − ν i ] is the Fourier transform of the mean-shifted process (ν i is the average firing rate of cell i) andf = F[f ] for all other quantities;K ij (ω) = A X,i (ω)J X,ij (ω) is the interaction matrix in the frequency domain (which may depend on the parameter being varied, i.e. X); ỹ 0 (ω)ỹ 0 * (ω) is the power spectrum of the uncoupled neural response. Using the usual series expansion for (I−K(ω)) −1 , we see that the covarianceC(ω) ≡ ỹ(ω)ỹ * (ω) naturally decomposes into contributions from different graph motifs: Each instance ofK includes the susceptibility function in the spectral domain, A X (ω), which modulates the effect of each directed connection by the responsiveness of the target neuron to the parameter of interest. As noted by previous authors [27], the validity of the expansion in Eq. (3) relies on the spectral radius ofK, ρ(K), remaining less than one. We checked that this remains true for all networks we examined in this paper, with a maximum of ρ(K) = 0.9564.
We next justify why-at least for long-time correlations-we can estimate each of these terms using the single-neuron firing rate function. First, in the small frequency limit, A X,i (ω) will coincide with the derivative of the firing rate with respect to the parameter: For the common input motif,K(ω)C 0 (ω)K * (ω), we can write where we have separated excitatory (E) and inhibitory (I) contributions, using g E and g I to denote the mean conductance of each type, and assumed that the synaptic coupling kernels,J jk , depend only on E/I cell identity (and have thus dropped the first subscript, which adds no additional information). This provides a formula for the long-time covariance, but we are typically interested in the correlation; fortunately, the small frequency limit also allows us to obtain a normalized correlation measure from the cross-spectrum, because where Cov T (n i , n j ) and Var T (n i ) denote covariance and variance of spike counts in a time window T ; that is, ρ T ,ij is the Pearson correlation coefficient (which varies between −1 and 1). Finally, letting ω → 0 and normalizing with the assumption that spiking is close to a Poisson process, as expected for an asynchronously firing network, so that C kk ≈ ν k : The above expression summarizes how the impact of excitatory and inhibitory common input are modulated by the single-neuron firing rate function, ν, and its derivatives.
Thus far, we have presented results previously obtained by others [27,40,41]. We now depart from these authors by focusing specifically on the behavior of the term in Eq. (9); and for simplicity, the behavior of this modulating factor for two identical neurons; e.g.
In principle, analogous expressions govern larger motifs, such as chains, that involve additional evaluations of ν and its derivatives; for example, excitatory length-3 chains arising fromK 3C0 would be: However, we found that, for a wide range of networks, direct common input-and inhibitory common input in particular-was the dominant contributor to pairwise correlations ( Fig. 6(A)). We now examine how different network mechanisms modulate the correlationfiring rate relationship, focusing on three factors: direction in effective parameter space, strength of recurrent excitation, and strength of background noise.

Direction in Parameter Space Modulates the Correlation-Firing Rate Relationship
Suppose we have a firing rate function of one or more intrinsic parameters (for exposition purposes, assume a function of two parameters (x 1 , x 2 )), i.e.
The parameters x j might include input bias, membrane time constant, ionic/synaptic reversal potentials, or a spiking threshold. Based on our arguments from the previous section, we will approximate correlation susceptibility by the quantitŷ where x 1 is an appropriately chosen parameter. Specifically, we will find, empirically, that the inhibitory common input is the dominant term, and therefore will use x 1 = g I , the mean inhibitory conductance. Thus, the units ofŜ in all figures are Hz/[g 2 I ], where g I is dimensionless (see Eq. (17)).
Heterogeneous firing rates can arise when each neuron occupies a different location in parameter space (i.e. a different (x 1 , x 2 ) point), thus causing both firing rate F and susceptibilityŜ to vary from neuron to neuron. We now ask: how doesŜ change with firing rate?
Note that this question, as stated, is ill-posed because F andŜ are both functions of two parameters (x 1 and x 2 ): there is not necessarily a one-to-one or even a functional relationship between these quantities. Suppose that, locally, a population of neurons can be described as lying along a parameterized path in the (x 1 , x 2 ) plane: i.e., (x 1 (s), x 2 (s)), for s min < s < s max . Then we can compute the directional derivative: where we have expressed the directional derivatives in terms of the local direction of the path: i.e. dx = ( dx 1 ds , dx 2 ds ) and the gradients of F andŜ. However, this depends not only on the functions F andŜ, but also on the direction dx.
To gain some intuition for why (and when) direction in (x 1 , x 2 ) space matters, we consider the networks studied in [4]. Previously, we simplified the firing rate function by setting all but two parameters (inhibitory conductance, g I,i , and threshold, θ i ) to their population average; i.e. F g I,i , θ i ≡ f g I,i , σ g I ,i p , g E,i p , σ g E ,i p , σ i p , θ i , (13) and · p denotes the population average. In Figs. 1(A) and (B), we show F andŜ ≡ ( ∂F ∂x 1 ) 2 /F thus computed, for the asynchronous network studied in that paper. We then surveyed a sequence of diagonal paths through the center (i.e., midpoint of the ranges of threshold and inhibitory conductance), with each path plotted in a different color. In Fig. 1(C) we plot firing rate (solid lines) and susceptibility (dashed-dotted lines) along each curve. Finally, in Fig. 1(D) we plot the susceptibility versus the firing rate, along each path.
This last panel shows that there is striking variability in the susceptibility-firing rate relationship, depending on the direction the path takes through the center of the (θ, g I ) plane. Any given firing rate (below ∼ 15 Hz) is consistent with either increase or decrease of susceptibility. All curves go through a single point in the (θ, g I ) plane, and therefore a single point in the (F,Ŝ) plane; here the direction Fig. 1 Firing rate and susceptibility (Ŝ), computed for the asynchronous (Asyn) network studied in [4], with all other parameters besides threshold θ and mean inhibitory conductance g I set to their average values (thereby leaving a two-dimensional parameter space). Here, we traverse the plane on straight-line paths defined by their angle through the origin. Although the units of g I are dimensionless, they are shown as the units forŜ for completeness. The units of θ (i.e., voltage) are also scaled to be dimensionless of theŜ-F relationship (i.e., whetherŜ increases or decreases with F ) can change rapidly with angle, as for the red, orange, and yellow curves.
We then extended these observations by traversing the phase space with two additional families of straight-line paths (Fig. 2); the radial paths are repeated in Figs. 2(A) and (B). For horizontal (Figs. 2(C) and (D)) and vertical (Figs. 2(E) and (F)) families, paths no longer intersect at a single point; nevertheless, a given firing rate is consistent with both increases and decreases in susceptibility. This is pronounced in Fig. 2(F), where at 15 Hz susceptibility decreases with firing rate in the orange, yellow and light green paths, but increases for the light blue, medium blue, royal blue, and indigo paths.
We performed the same computations on the strong asynchronous network studied in [4] that has larger excitatory coupling strength: results are shown in Fig. 3. A given firing rate could be consistent with either increase or decrease of susceptibility; we see this in the radial paths (Figs. 3(A) and (B)) and horizontal paths (Figs. 3(C) and (D)) for rates between 5-10 Hz. However, vertical paths (Figs. 3(E) and (F)) always have susceptibility increasing with firing rate (except perhaps at the highest firing rates). As in the asynchronous network, direction of susceptibility (increase vs. decrease) can change rapidly with angle, for paths that intersect a single point; see Figs. 3(A)-(B), red to orange to yellow.  [4], with all other parameters besides threshold θ and mean inhibitory conductance g I set to their average values (thereby leaving a two-dimensional parameter space: the other (averaged) parameters are σ g I ,i p = 0.6602 (see Eq. (30)), σ g E ,i p = 0.0026 (see Eq. (29)), σ i p = 6.3246, g E , i p = 0.0053 (see Eq. (17)). Here, we traverse the plane on three different families of straight-line paths. The dashed lines show paths visualized in [4]: θ = 1 (vertical, aqua blue) and g I = 1.83 (horizontal, orange/yellow)

Quantifying the Likelihood of a Positive Correlation-Firing Rate Relationship
In the previous section, we saw that the path a network occupies in effective parameter space can have a dramatic effect on the correlation-firing rate relationship: we next seek to formalize these observations. Let dx be the local direction in which we want to  [4], with all other parameters besides threshold θ and mean inhibitory conductance g I set to their average values (thereby leaving a two-dimensional parameter space: the other (averaged) parameters are σ g I ,i p = 0.5884 (see Eq. (30)), σ g E ,i p = 0.0378 (see Eq. (29)), σ i p = 4.7434, g E , i p = 0.0611, see Eq. (17)). Here, we traverse the plane on three different families of straight-line paths. Dashed lines show paths visualized in [4]: θ = 1 (vertical, aqua blue) and g I = 1.46 (horizontal, yellow/green) consider the behavior of F andŜ. If ∇Ŝ · dx and ∇F · dx are of the same sign, thenŜ increases with F ; if ∇Ŝ ·dx and ∇F ·dx have opposite signs, thenŜ decreases with F .
The more aligned ∇Ŝ and ∇F , the more paths lead to dŜ dF > 0; the more anti-aligned ∇Ŝ and ∇F , the more paths lead to dŜ dF < 0. For example, consider the limiting cases where: (i) if ∇Ŝ and ∇F point exactly in the same direction, then ∇Ŝ · dx and Where ∇Ŝ and ∇F are similarly aligned,Ŝ and F will both increase along most paths through that point. In each panel, gray shows the part of the x-plane where dŜ dF = ∇Ŝ·dx ∇F ·dx > 0, black where dŜ dF < 0. From left to right: ∇Ŝ and ∇F positively aligned; ∇Ŝ and ∇F orthogonal; ∇Ŝ and ∇F negatively aligned ∇F · dx are always same-signed for any dx; (ii) if ∇Ŝ and ∇F point in opposite directions, then ∇Ŝ · dx and ∇F · dx are never same-signed. Figure 4 illustrates how the alignment of these two quantities alters the region where correlation increases with firing rate.
To examine the utility of this idea, we reconsider the radial paths along which we previously computed susceptibility. The paths all go through a single point, so we can check the sign((∇Ŝ · x)(∇F · x)) for all directions through this point. In Figs. 5(A) and (C), white indicates positive and black negative. Figures 5(B) and (D) repeat the susceptibility-firing rate curves from Fig. 2(B) and Fig. 3(B). For the asynchronous network (Fig. 5(A)), the red, indigo, and royal blue paths are predicted to have negative dŜ/dF , as we can confirm in Fig. 5(B). Yellow, green, and light blue curves are predicted to have positive dŜ/dF . The orange curve is close to dF = 0; the true blue curve is close to dS = 0. For the strong asynchronous network, only the red curve is in the negative region of Fig. 5(C); this is also the only path with dŜ/dF < 0 in Fig. 5(D).
Of course, this prediction only applies to portions of the path near the point at which we computed the gradients; away from this point, gradients will change along with the direction of theŜ vs. F curve. For example, the royal blue curve in Fig. 5(B) increases with firing rate for small firing rates, and the light blue, true blue, and royal blue curves in Fig. 5(D) decrease with firing rate, (for large firing rates).
This motivates a direction-independent measure to quantify the fraction of paths that lead to an increase of correlation with firing rate. This is given exactly in terms of the angle between ∇Ŝ and ∇F : and in particular the fraction of x directions in whichŜ increases with F is given by Because cos −1 has range [0, π], this varies between 0 and 1. The more aligned ∇Ŝ and ∇F , the more paths lead to dŜ dF > 0; the more anti-aligned ∇Ŝ and ∇F , the more paths lead to dŜ dF < 0.

Strength of Recurrent Excitation Modulates the Correlation-Firing Rate Relationship
Our use of inhibitory susceptibility (i.e. Eq. (10)) to characterize the relationship between correlations and firing rates relied on intermediate assumptions, specifically: • Second-order motifs dominate pairwise correlations.
• Inhibitory common input is the dominant second-order motif.
Here, we check that these conditions are still satisfied for a range of neural network models.
In [4], we examined two spiking regimes achieved by varying the strength of excitation: both recurrent excitation W EE and excitatory input into the inhibitory population W I E . We next applied our theory to a dense grid of parameters (different networks), each identified by its location on the (W EE , W I E ) plane. Both excitatory strengths were varied from a minimum of their values for the asynchronous network (W EE = 0.5 and W I E = 5) to a maximum of 1.4 times their value in the strong asynchronous network (to W EE = 12.6 and W I E = 11.2). The parameter W XY is a dimensionless scale factor (see Eqs. (17)-(20)); for example in an all-to-all homogeneous network, W XY = 1 is when the presynaptic input is exactly the average population firing rate (filtered by the synapse).
To quantify the importance of paths of different length through the network, we can define the contributions at any specific order k by using the interaction net- Then we regressed the total correlation (C ij / C iiCjj ) against the contributions at each specific order (R k ij ); the corresponding fraction of variance explained (R 2 value) gives a quantitative measure of how well the total correlation can be predicted from each term.
Similarly, we quantity the importance of specific second-order motif types, by regressing the total contribution from second-order motifs (R 2 ij ) against the contribution from specific motifs. We performed this computation for each network (a total of 225 networks), and summarize the results in Fig. 6; to present the data compactly, we collapse R 2 values (all values are ∈ [0, 1]) for a fixed W EE and varying W I E by showing mean and standard deviation only (standard deviation as error bars). Secondorder contributions dominate for small to moderate W EE (Fig. 6(A)); other orders only compete when W EE has already exceeded the level of the strong asynchronous network (where the network is close to the edge of instability, and at the limit of validity for mean-field, asynchronous assumptions).
To illustrate the meaning of this statistic, in Fig. 6(B) we show contributions up to fourth order (R k ij , for k = 1, . . . , 4) vs. total correlation (C ij / C iiCjj ) for all E-E cell pairs in a network, for two individual networks included in the summary panel. Note that the second-order contributions cluster near the unity line in both cases, indicating that second-order contributions are the best predictor of total correlations, consistent with the R 2 values stated.
Of the second-order motifs, inhibitory common input is the dominant contribution at any value of W EE , except perhaps the last, at which point the excitatory population has unrealistically high firing rates (Fig. 6(C)). To summarize, we have confirmed that far from being limited to a few examples, the conditions we identified in [4] as ij ) against contributions from the distinct types of second-order motifs: inhibitory common input (magenta), excitatory common input (red), decorrelating chains (green), and correlating chains (blue). Both (A,C): Each data point represents the mean value from 15 networks with W I E between 5 and 11.2; error bars show standard deviation across these values. W XY = 1 is when the presynaptic input is exactly the average population firing rate (filtered by the synapse) in an all-to-all homogeneous network allowing us to focus on susceptibility to inhibition to explain pairwise correlations, appear to hold up for a variety of networks.
We note that our findings about the relative magnitudes of terms of different orders are purely empirical; that is, they are based on concrete numerical observations, rather than a priori estimates. Thus, it should be reassessed if anything about the parameters or network connectivity is changed. In particular, a likely reason for the relative weakness of first-order terms is that in these networks excitation is almost always weaker than inhibition; while first-order terms are always excitatory, second-order terms can involve excitation and/or (comparatively strong) inhibition.
Having confirmed the validity of our approach, we computed the susceptibility with respect to inhibition, for each of the networks examined in the previous section (because of instability, we restricted these computations to excitatory strengths ×1.2 the values used in [4]). Because background noise values varied slightly, we actually examined two network families; one in which we chose σ values as in the asynchronous network, another in which we chose σ values as in the strong asynchronous network. Also as in [4], we estimated susceptibility using network-averaged values of g E , g I , σ g E , and σ g I . Surprisingly, the difference in background noise dwarfed the recurrent excitation strengths, at least in accessing the relationship betweenŜ and firing rate. In Fig. 7, we showŜ vs.   [4]), and 8.6. Again, W XY = 1 is when the presynaptic input is exactly the average population firing rate (filtered by the synapse) in an all-to-all homogeneous network, so the coupling strengths vary significantly W XY = 1 is when the presynaptic input is exactly the average population firing rate (filtered by the synapse) in an all-to-tall homogeneous network, so these coupling strengths vary significantly. We see that, for the full range of recurrent excitation values,Ŝ vs. F curves in Fig. 7(A) are mostly decreasing;Ŝ vs. F curves in Fig. 7(B) are mostly increasing for low F , and saturating for high F .

Background Noise Modulates the Correlation-Firing Rate Relationship
To further explore the role of background noise, we repeated the susceptibility calculation on additional families of networks, now allowing background noise strengths σ E and σ I (i.e. to the excitatory and inhibitory populations) to vary separately. Input to excitatory cells was varied between σ E = 1.5 and 2.5; input to inhibitory cells was varied between σ I = 1.5 and 3. These noise values are relatively large; see Eq. (17) and note that voltage is of order 1. In Fig. 8(A) we display susceptibility vs. firing rate curves for 12 (σ E , σ I ) pairs; asterisks indicate σ E and σ I by color (green intensity for σ E and blue intensity for σ I ). Within each panel curves are colored as in Fig. 7: red intensity for W I E and blue intensity for W EE .
Surprisingly, most network families (i.e. (σ E , σ I )) were associated with a decrease in correlation with firing rate. The exceptions are (0.15, 0.25) (as in the strong asynchronous network in [4]) and (0.15, 0.3). This latter was most robustly associated with a positive correlation-firing rate relationship. Furthermore, the shape of susceptibility-firing curves did not appear to vary much with the strength of recurrent excitation (i.e., curves within each panel are similar). We next sought to investigate possible mechanisms for a positive correlation-firing rate relationship, by examining the effective parameters that govern the neural response: in essence, the network's "operating point" (see Eq. (26)). Possible choices include g I , g E , σ g E , σ g I , and the effective reversal potential E rev ; we found σ g E and σ g I to be largely functions of g E and g I , while E rev has a (nonlinear) functional relationship with g E and g I . Thus two parameters suffice, and here we chose to use g I and E rev . In Fig. 8(B), we plot average parameter values for each curve, colorcoded by (σ E , σ I ). Any given color is consistent with a relatively tight range of g I and (comparatively) broad range of E rev . As σ I increases (increasing blue intensity), inhibitory conductance g I increases and reversal potential E rev decreases. However, it was not apparent that any particular region in (g I , E rev ) parameter space was associated with a positive correlation-firing rate relationship, in that the values of g I and E rev that supported a positive relationship supported negative relationships as well.

Discussion
In this paper, we showed that using a single-cell firing rate function to examine the relationship between correlations and firing rates is feasible for a wide range of heterogeneous, recurrent networks. We focused on three factors that can modulate the correlation-firing rate relationship: how the network occupies effective parameter space, strength of recurrent excitation, and strength of background noise. Although there are many sets of parameters known to vary within a heterogeneous network of neurons, we have already demonstrated vastly different correlation-firing rate relationships with our methods, with a theory that can be readily applied to other model networks.
One possible application of this work is in designing neural networks for computational experimentation; just as modelers now modify cortical networks to obey experimental constraints on firing rates [3,35], we could also include a constraint on the desired correlation-firing rate relationship. Here we showed that we can quickly assess a wide range of possible network configurations for a positive correlationfiring rate relationship: in Sect. 2.5, for example, we performed calculations on 15 × 15 × 12 = 2700 heterogeneous networks, with a nominal amount of computing time.
One surprising finding in our computations was the relative insensitivity of the slope of the correlation-firing rate relationship to recurrent excitation (W EE , W I E ), as demonstrated in Figs. 7 and 8. This is striking in contrast to the strong sensitivity on display in Figs. 2(B) and 3(B). This difference is explained as follows: in every case where we computed the susceptibility for a self-consistent network (i.e. a solution of Eqs. (26)-(30) and (32)-(33)), the source of heterogeneous firing rates was neural excitability, set via a spiking threshold θ . The resulting effective parameters, such as inhibitory conductance g I , did not deviate strongly from their population mean values. In essence, all of these networks took a horizontal path through (θ, g I ) parameter space, as in Figs. 2(C), (D) and Figs. 3(C), (D). If we were to generate networks where heterogeneity arises from another source-such as the strength or frequency of inhibitory connections [23]-we might see different results. We look forward to exploring this in future work.
A priori, there is no reason to expect that the correlation-firing rate relationship in these recurrent networks can be simplified to a feedforward motif based on shared inhibitory input; this was purely an empirical observation (see Fig. 6(B)). We remark that others have shown that the effective input correlation can be canceled to have near zero input (and thus output) correlation on average in balanced networks [28,40], in contrast to some of the models considered here (i.e., strong asynchronous regime with more net excitation). The conditions for correlation cancellation in these model networks is beyond the scope of this study, but note that others have shown correlation cancellation does not always hold ( [21,22] via altering connection probabilities). The studies that demonstrate correlation cancellation often have faster (or equal) inhibitory synaptic time scales than excitatory: τ I ≤ τ E [21,28,32] ( [40] used current-based instantaneous synapses or τ I = τ E = 10 ms) while in our networks the inhibitory synapses have longer time scales (Table 1). Note that Fig. S4 of [28] shows that having effectively zero input correlation does not hold as the inhibitory time scales increase beyond the excitatory time scales. Finally, system size is another key parameter that could certainly affect the magnitude of the recurrent feedback. In contrast to [28,40], we did not account for how system size would affect correlation cancellation in these heterogeneous networks.
Although affirmative answers to whether correlations increase with firing rate in experiments were cited in the Introduction [2,5,8,11,19,36,44] we also note that many experiments have shown that the average correlation across a population can decrease with firing rate when a global state changes or a stimulus is presented. A recent review [9] shows that stimulus-induced decorrelation (with increased firing rate) occurs in a variety of brain regions and animals. This is slightly different from the situation we examine here, where we consider the relationship between correlations and firing rates within a stimulus condition. Regardless, the fact that the relationship between correlation and firing rate is not obvious points to the continued need for theoretical studies into the mechanisms of spike statistics modulation.
Finally, our finding that correlations often decrease, rather than increase, with firing rate stands in apparent contradiction to earlier work on feedforward networks [8,38]. On closer inspection, we can identify several reasons why our results differ; with conductance inputs (rather than currents) we have a quantitatively different relationship between input parameters and firing rates; furthermore with more adjustable single-neuron parameters, the same sets of firing rates may be observed with single-neuron parameters set in different ways. While the current clamp experiments described in [8] found a consistent increase of correlations with firing rates, we can hypothesize that the parallel dynamic clamp experiments in which pairwise correlations arise from common inhibitory input, would in fact show a decrease with firing rate. However, we also predict that whether an increase or decrease with firing rate is observed would depend on whether firing rates are modulated by varying the level of inhibitory input, or by otherwise varying the excitability of the cells (perhaps pharmacologically).

Neuron Model and Network Setup
We considered randomly connected networks of excitatory and inhibitory neurons. Each cell is 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 the voltage is below threshold θ i : When v i reaches θ i , a spike is recorded and voltage is reset to 0 following a refractory period: Each neuron receives 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 (see Table 1). The thresholds θ i are a significant source of heterogeneity, and they are selected from a log-normal distribution with mean 1 and variance e (0.2) 2 − 1; since the system size is moderate, the θ i 's were set to have C.D.F. (cumulative distribution function) values equally spaced from 0.05 to 0.95 for both E and I cells.
Each cell responds to synaptic input through conductance terms, g E,i and g I,i , which are each governed by a pair of differential equations: 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 (1) X,i ; here t j,k is the kth spike of 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, that is 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 2) is divided by N Y X , the total number of X → Y connections received by each Y -type cell.  Fig. 8, W EE was varied between 0.5 and 10.8 and W I E between 5 and 9.6; σ E was varied between 1.5 and 2.5 and σ I between 1.5 and 3.

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 [31]. 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 spiking activity. The derivation is presented in full in [20,29,30]; here, we present only the main points. We assume we have some way to approximate the change in firing rate which occurs as a result of a change in parameter: ν i,0 is the baseline rate (when X = 0) and A X,i (t) is a susceptibility function that characterizes this firing rate response up to order [8,20,41]. In order to consider joint statistics, we need the trial-by-trial response of the cell. First, we propose to approximate the response of each neuron by that is, each input X i has been replaced by a filtered version of the presynaptic firing rates y j .
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) andf = F[f ] for all other quantities. In matrix form, this yields a self-consistent equation forỹ in terms ofỹ 0 : whereK ij (ω) =Ã X,i (ω)J X,ij (ω) is the interaction matrix in the frequency domain. The cross-spectrum is then computed via To compute the interaction matrix for a network of conductance-based neurons, we use the effective time constant approximation (as in the supplemental for [41]). We first separate each conductance into mean and fluctuating parts, e.g., g E,i → g E,i + (g E,i − g E,i ) (see the discussion in [12]). Next we identify an effective conductance g 0,i and potential E rev,i , and treat the fluctuating part of the conductances as noise, i.e. g E,i − g E,i → σ g E ,i ξ E,i (t), so that Eq. (17) becomes where The parameters which govern the firing rate response will now be the conductance mean and variance, e.g. g E,i and σ 2 g E ,i ; we next compute how these depend on incoming firing rates for second-order α-function synapses (Eqs. (19) and (20)). We first simplify the equation for the auxiliary variable (Eq. (20)): so thatα X,i includes all factors that contribute to the pulse size in Eq. (20), including synapse strength and pulse amplitude. The time constants τ r,X , τ d,X and synapse jump sizesα X,i generally depend on cell type. Then assuming that each spike train is a Poisson process with a constant mean firing rate: i.e., each spike train is modeled as a stochastic process S(t) with a straightforward but lengthy calculation shows that g X,i (t) =α X,i ν X,i τ r,X , Var g X,i (t) = 1 2α 2 X,i ν X,i τ r,X τ r,X τ r,X + τ d,X where ν X,i is the total rate of type-X spikes incoming to cell i. Notice that modulating the rate of an incoming spike train will impact both the mean and variance of the input to the effective equation, Eq. (26) (via E rev,i and σ g 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. Therefore, we have four susceptibility functions to compute,Ã g E ,i (ω), A g I ,i (ω),Ã σ 2 g E ,i (ω), andÃ σ 2 g I ,i (ω). The first two capture the change in firing rate as a result of a change in mean conductanceg E,i → g E,i 0 + g E,i 1 exp(ıωt) or g I,i → g I,i 0 + g I,i 1 exp(ıωt)-while the final two address a change in variance-σ 2 g E ,i → (σ 2 g E ,i ) 0 + (σ 2 g E ,i ) 1 exp(ıωt) or σ 2 g I ,i → (σ 2 g I ,i ) 0 + (σ 2 g I ,i ) 1 exp(ıωt). Since the corresponding Fokker-Planck equation required to obtained these entities is linear, we can compute both susceptibilities separately and combine them to get the net effect. With these pieces, we now have the interaction matrix:K ij (ω) = Ã g E ,i (ω)J ij (ω) +Ã σ 2 g E ,i (ω)L ij (ω), j excitatory, A g I ,i (ω)J ij (ω) +Ã σ 2 g I ,i (ω)L ij (ω), j inhibitory, (34) whereL(ω) plays a similar role asJ, but for the effect of incoming spikes on the variance of conductance. Its relationship toJ (either in the frequency or time domain) is given by the same simple scaling shown in Eq. (33): i.e., for j excitatory, where the first factor comes from the effective spike amplitudeα E,i (and is the scale factor proposed in [29], Eq. (64)), and the second arises from using second-order (vs. first-order) alpha-functions.

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. (25)). The cross-correlation function, C ij (τ ), measures the similarity between two processes at time lag τ , while the crossspectrum measures the similarity between two processes at frequency ω: C ij (ω) ≡ ỹ i (ω)ỹ j (ω) .
The Weiner-Khinchin theorem [31] implies that {C ij ,C ij } are a Fourier transform pair: that is,C In principle, the cross-correlation C(t) and cross-spectrumC(ω) 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 cross-correlation function C ij by the following formula [17]: The variance of spike counts over a time window of length T , n i , is likewise given by integrating the autocorrelation function C ii : By normalizing by the time window and taking the limit as T → ∞, we can see that, for an integrable cross-correlation function, we can useC ij (0) as a measure of long-time covariance. Similarly, the long-time limit of the Pearson correlation coefficient of the spike counts, gives us a normalized measure of long-time correlation.