Effect of Synaptic Heterogeneity on Neuronal Coordination

Recent advancements in measurement techniques have resulted in an increasing amount of data on neural activities recorded in parallel, revealing largely heterogeneous correlation patterns across neurons. Yet, the mechanistic origin of this heterogeneity is largely unknown because existing theoretical approaches linking structure and dynamics in neural circuits are restricted to population-averaged connectivity and activity. Here we present a systematic inclusion of heterogeneity in network connectivity to derive quantitative predictions for neuron-resolved covariances and their statistics in spiking neural networks. Our study shows that the heterogeneity in covariances is not a result of variability in single-neuron firing statistics but stems from the ubiquitously observed sparsity and variability of connections in brain networks. Linear-response theory maps these features to the effective connectivity between neurons, which in turn determines neuronal covariances. Beyond-mean-field tools reveal that synaptic heterogeneity modulates the variability of covariances and thus the complexity of neuronal coordination across many orders of magnitude.


I. INTRODUCTION
Neuronal networks in the brain display largely heterogeneous activity: common observables such as firing rates [1][2][3], coefficients of variation (CVs) [4], and pair-wise correlations [5][6][7][8] are widely distributed across neurons.This has important implications for coding and information processing in the brain, as the coordinated activity across the enormous number of units in neuronal circuits is thought to underlie all complex functions [9][10][11][12].The causes of heterogeneity in neuronal dynamics are diverse: intrinsic neuronal properties, external inputs, and the network connectivity itself are all sources of variability.While these structural and dynamic heterogeneities can be readily observed with modern experimental techniques [13][14][15], understanding their mechanistic relations requires theoretical tools that are currently still lacking.
In this study, we focus on the effects of connectivity and investigate the influence of heterogeneity in connections on the activity of networks of identical neurons receiving homogeneous external input.Previous work [16] has shown that a considerable fraction of the variance, in the distribution of firing rates across neurons and in the CV of individual neurons' spike trains, in such networks can already be explained by the distributed number of inputs the neurons in a network receive.In this study, we go beyond single neuron activities and focus on the statistics of pair-wise correlations and the related covariances, which measure how strongly the activities of pairs of neurons co-fluctuate.Such coordination builds the basis for collective network activity and function.
With the exception of small organisms such as Caenorhabditis elegans [17], the microconnectome of most biological neuronal networks is unknown.However, overall connectivity properties and statistics, like the connection probabilities between different cortical areas [18,19] and cell types [15], the distance dependence of connections [15,20], or the statistics of synaptic strengths [15,[21][22][23][24][25][26] are available nowadays.Hence, rather than a one-to-one relation between microconnectome and pairwise covariances [27][28][29][30][31][32], a relation between connectivity and covariance on a statistical level would readily allow the inclusion of this knowledge.To derive such a relation, common population-level theories [27,[33][34][35][36][37][38] cannot be used because they can only describe population-averaged observables and, in particular, do not capture heterogeneity in covariances within populations.Here, we instead employ mean-field theory on the single-neuron level [30], which we systematically compare to network simulations, and we go beyond mean-field theory by including nontrivial fluctuation terms to obtain the statistics of covariances between individual neuron pairs.The main difficulty of a single-neuron level approach is that the predictions of the theory for individual neurons strongly depend on the specific details of the connectivity.To get a description on the level of connectivity statistics, we perform a disorder average, a technique originally developed for spin-glass systems [39,40] that allows retaining information about the connectivity statistics while averaging over the realization randomness.As our main results, we show how to systematically calculate higher moments of neuronal activity averaged over the disorder in the connectivity using replica and beyond-mean-field theory, and we use this technique to derive a relation between the mean and variance of covariances and the mean and variance of the network connectivity.First results based on a similar but reduced theoretical approach have already been successfully applied in the neuroscientific context to infer the dynamical regime of cortical networks [7] and to explain spatial properties of coordi- nation structures [8] and dimensionality [41].
To summarize, we investigate the origin of neuronal coordination structures, as experimentally observed across various species and cortical areas, by analyzing covariances in a prototypical network model of cortical dynamics [42], namely sparsely connected excitatory and inhibitory neurons that operate in the balanced state [34].In this model, all neurons have identical parameters and receive homogeneous, uncorrelated external input.As in biological cortical networks, the sparsity in the connectivity between neurons [15] as well as the wide distribution in synaptic amplitudes [15,[21][22][23][24][25][26] constitute the source of variability in connections and thereby the dynamics: Rates, CVs, variances, covariances, and hence correlation coefficients are all described by distributions with sizable variance (see Fig. 1).
The following sections investigate the sources of the variance in these quantities.Sec.II introduces mean-field theory on the single neuron level.In Sec.III, we derive the main results on how to compute disorder-averaged moments of neuronal activity, and we calculate explicit expressions for the mean and variance of covariances.In Sec.IV, we discuss our findings and their limitations in the context of the existing literature.

II. BACKGROUND: LINEAR-RESPONSE THEORY OF SPIKING NEURONAL NETWORKS ON A SINGLE-NEURON LEVEL
To understand the origin of the distribution of covariances, we start with analyzing a simulated network on a single-neuron level.The example network throughout this study comprises 8000 excitatory (E) and 2000 inhibitory (I) leaky integrate-and-fire (LIF) neurons that make connections according to distinct populationspecific statistics.We mimic two abundant features of heterogeneity in connectivity of brain circuits, that are sparse connections and distributed synaptic weights.To do so we consider random sparse connectivity J , with connection probability 10 %, giving rise to an excitatory indegree K E = 800 and an inhibitory indegree K I = 200.To compensate for the imbalance in excitatory and inhibitory neuron count, we follow the work by Brunel [42] and scale the strengths of existing inhibitory connections with respect to excitatory ones by a factor g = −6 to obtain an asynchronous irregular dynamic regime.In addition, we distribute synaptic weights of existing connections according to population-specific normal distributions j E ∝ N (j, 0.2j) and j I ∝ N (gj, 0.2j), such that the overall heterogeneity in network connectivity is comprised of the random sparseness of connections and the variable strength of connections.For more details of the model, see Appendix IV A.
Working point Given the parameters of the simulated network of leaky integrate-and-fire neurons, especially the specific realization of the connectivity matrix J , we determine the stationary working point, comprising the input statistics (µ, σ) and the firing rates ν, as done by Brunel and Hakim [43] and Brunel [42].To this end, we first neglect correlations between the neurons and approximate the neurons' inputs as independent Gaussian white noise processes.In this diffusion approximation, the mean input µ i and input variance σ 2 i of neuron i are given by with membrane time constant τ m , membrane capacitance C, constant input current I ext , and excitatory and inhibitory external Poisson noise with rates ν ext,E and ν ext,I which are fed into the system via weights j and gj, respectively.The firing rates are given by the Siegert function [44] with refractory period τ r , and rescaled reset and threshold voltages These equations can be solved iteratively in a selfconsistent manner.Given the working point, we can determine the coefficients of variation using (see Appendix A.1 of Ref. [42]; note that they use different units) (4) Linearization The full dynamics of LIF neurons are nonlinear.However, as covariances measure cofluctuations of neurons around their working points, we can study covariances by analyzing linearized dynamics as long as the fluctuations are sufficiently small.Grytskyy et al. (see Sec. 5 of Ref. [31]) show that a network of LIF neurons can be mapped to a linear rate model with output noise with neuronal activity x i (t), normalized linear response kernel h(t), synaptic delay d, and uncorrelated Gaussian white noise ξ i (t), ⟨ξ i ⟩ = 0, ⟨ξ i (s) ξ j (t)⟩ = D ij δ (s − t), with diagonal noise strength matrix D ij = δ ij D ii .The matrix W , referred to as effective connectivity, combines the connectivity matrix J with the sensitivity of neurons to small fluctuations in their input.It is formally given by the derivative of the stationary firing rate of neuron i [Eq.(3)] with respect to the firing rate of neuron j evaluated at the stationary working point (see Appendix A of Ref. [45]) with Spike-count covariances In this study we are interested in spike-count covariances in spiking networks with spike counts n i occurring within bins of size T , where the average, indicated by the brackets ⟨⋅⟩, is taken across all bins that can be viewed as trials with different realizations of the external input.As shown in Dahmen et al. [7,Materials and Methods], for stationary processes and large bin sizes spike-count covariances C ij can be mapped to the time-lag integrated covariances c ij (τ ) between spike trains of neurons i and j (see also 46, 47, for more details on definitions of covariances see Appendix IV B) In the following the term covariance always refers to C ij .Making use of the Wiener-Khinchin theorem (IV C) allows expressing the time-lag integrated covariances in terms of the neuronal activities' Fourier components X i (ω) at frequency zero which can be evaluated by Fourier transforming Eq. ( 5), yielding For calculating the covariances, we therefore only need the effective connectivity W and the noise strength D.
The correlation coefficients follow as To estimate the noise strength D, we assume that the spike trains are described sufficiently well as renewal processes for which the variances are given by [48] Using that D is by definition diagonal (see Appendix IV D for limitations on exactly matching simulated spikecount covariances with a linear rate model with uncorrelated white noise input), we can solve Eq. ( 9) for D, which results in with The above expressions can be combined to compute theoretical estimates of the quantities measured in the simulation.To solve the self-consistency equations for the firing rates and to compute the covariances, we make use of the Python package NNMT [49], which includes optimized implementations of the equations introduced above.A comparison of theoretical and simulation results is shown in Fig. 2. For the chosen parameters, simulation and theory correlate strongly, and the theory appears to capture the primary sources of heterogeneity in the rates, covariances, and correlation coefficients.Note that such a good match between theory and simulation cannot be observed in all parameter regimes of the spiking network; the validity of the assumptions made and the resulting theoretical estimates depend on the network state (see Appendix IV D for further discussion on valid parameter regimes).Fig. 2 also reveals some unexplained variance, particularly pronounced in the covariances and correlations.This variance is the result of the finite simulation time and the associated uncertainty in the estimated covariances.As we show in Appendix IV K, the covariance estimate bias can be significant and it can only be corrected for on a statistical level rather than for individual covariances.Focusing on the statistics of covariances, however, has further advantages: For realistic network sizes, Eq. ( 9) is a high-dimensional equation that depends on each and every connection in the network.Understanding general mechanisms relating network structure and dynamics is therefore difficult.The covariance statistics instead summarize the most important aspects of covariances and, for large neuron populations, can be assumed to be self-averaging [40,50,51], which makes them less dependent on connectivity details.Second, Eq. ( 9) cannot be used for inference based on experimentally measured parameters because as of yet it is not possible to determine the effective connectivity or covariances of all neurons in a network.And last, as stated above, we demonstrate that covariance statistics are more robust measures than single-neuron covariances, both with respect to finite measurements as well as to the assumptions made in the derivation above.

III. STATISTICAL DESCRIPTION OF COVARIANCES
Expression (9) reveals that the statistics of the covariances C, in particular their heterogeneity, is determined by the statistics and heterogeneity of the effective connectivity matrix W and the external noise strength D. Our aim here is to derive a description of the cross-covariance statistics in terms of the statistics of W and D. To this end, we derive analytical expressions for the mean and the variance of the time-lag integrated cross-covariances averaged over the heterogeneities of the system.
To do this, simply averaging Eq. ( 9) is not feasible due to W appearing in the inverse matrix (1 − W ) −1 .Performing an average over a random connectivity is, however, a well known problem in the theory of disordered systems [40,50,52,53], where it is handled on the level of generating functions.To proceed analogously, we start with Eq. ( 8), which expresses the covariances in terms of the moments of the dynamic variables' Fourier components at frequency zero.This allows us to write the covariances in terms of the moment-generating function Z (J ) of the zero-frequency Fourier components X i of the dynamical equation ( 5) (see Appendix IV E for more details): and of the unnormalized moment-generating function Z (J ).
Here, X are auxiliary variables that can be used to calculate the response function ⟨X i Xj ⟩ of neuron i to a perturbation of neuron j by introducing additional sources J in the moment-generating function (see Appendix IV E).Equation (13) shows that calculating the disorder-average of the covariances boils down to calculating the disorder-average of the moment-generating function.In the following two sections, we use this approach to calculate the mean of the cross-covariances ⟨C⟩ W ,D and subsequently the variance of the crosscovariances ⟨δC 2 ⟩ W ,D , where ⟨⋅⟩ W ,D refers to the average over the randomness in W and D.

A. Mean of cross-covariances
Disorder average We begin with the mean crosscovariances, focusing first on the average over the ensemble of connectivities, indicated in the following by ⟨⋅⟩ W .In the moment-generating function [Eq.( 13)], W occurs linearly in the exponent of Z (J ), which is advantageous for performing the disorder average.However, the averaging procedure is complicated by two aspects: (1) W contributes to the noise strength D through the variance-rescaling matrix B −1 , and (2) the normalization Z (0) depends on W .However, as illustrated in Fig. 3, in practice the first point does not appear to be a problem: Fig. 3(a) indicates that the specifics of D are largely determined by the details of the variances CV 2 ν, because a different realization of W essentially yields a similar D, and Fig. 3(b) suggests that the effect of the disorder average on D is minimal.For these reasons, we treat D as though it was independent of the explicit realization of W .To address the second point, an alternative approach based on the moment-generating functional for the full time-dependent dynamics (see Appendix IV E) could be utilized.This moment-generating functional has a unit determinant normalization independent of W [54].The disorder average of its frequency space complement, however, introduces cross-frequency couplings that complicate the further analysis.Here, instead, we follow Dahmen et al. [7], and separate the averages over Z (J ) and Z (0), as we find that this factorization approach does yield accurate results.This leaves us with the task of calculating The disorder average only affects the coupling term and can be expressed using the moment-generating functions for independently drawn connections The moment-generating function can be written in terms of a cumulant expansion ϕ ij (X) = exp (∑ ∞ k=1 κ k,ij X k /k!), with k-th cumulants κ k,ij .For fixed connection probability, the number of inputs to a neuron scales with the network size N .To keep the input and its fluctuations finite when increasing the network size, we require synaptic weights to scale with 1/ √ N [34,55], such that the cumulant expansion is an expansion in 1/ √ N .A truncation at the second cumulant (∝ N −1 ) maps W to a Gaussian connectivity with distribution N (M , ∆/N ), such that × exp [S 0 (X, X) + J T X] 1 Note that a systematic approach to this factorization approximation would be to employ the replica trick Z (0) −1 = lim n→0 Z (0) n−1 and jointly average ⟨ Z (J) Z (0) n−1 ⟩ W in the limit n → 0, or to average the joint moment-generating functional for all time points or frequencies that by construction has a trivial normalization Z (0) = 1 (see IV E).  √ N ) as well as variances ∆ ij = O(1).The higher-order cumulants are suppressed by the large network size and therefore do not significantly affect the dynamics of the network.As we show below, our theory with the Gaussian approximation of the connectivity therefore faithfully recovers the correlation statistics in the sparse excitatory-inhibitory spiking network with population-specific connection statistics.
Auxiliary field formulation To deal with the fourpoint coupling term in Eq. ( 15), we define auxiliary variables Q i ∶= 1 N ∑ j ∆ ij X j X j , which we formally introduce by inserting an identity in the form of a Fouriertransformed delta distribution: .
The auxiliary variables Qi are introduced to express the delta distribution as an integral.This leads to Here diag (Q) ij refers to a diagonal matrix with diagonal elements Q i .As the action S Q, Q (X, X) at fixed auxiliary variables describes an auxiliary free theory, Eq. ( 17) describes the activity of linear rate neurons in a network with disorder-averaged connectivity M that interact with fluctuating external variables Q and Q. Inserting Eq. ( 17) into Eq.( 14) yields ∫ DP ∫ D P exp (−N P T P ) ZP , P (0 with joint probability distribution ∫ DP ∫ D P exp (−S (P , P )) and properly normalized moment generating function . These equations imply that the disorder-average of arbitrary moments ⟨X i1 ⋯X i k ⟩ can be calculated by determining the corresponding moments ⟨X i1 ⋯X i k ⟩ Q, Q with respect to the auxiliary free theory and averaging them over the auxiliary variables: Saddle-point approximation Due to the prefactor N in Eq. ( 19) and the scalar products in ZQ, Q (0) with N contributions, we expect p (Q, Q) to peak sharply for N → ∞, such that we can perform a saddle-point approximation.To lowest order, we expect which yields with second moments evaluated at the saddle point.The moments can be calculated explicitly by solving the Gaussian integrals (see Appendix IV F).Using the and solving for the saddle point yields Finally, making use of the Wiener-Khinchin theorem [Eq.(40)] and inserting the solution of the saddle-point equations into Eq.( 20) yields the mean covariances averaged across the disorder of the connectivity: Averaging over the disorder in D then yields Here D denotes the disorder-averaged noise strength [cf.Fig. 3(b) and discussion after Eq. ( 31)].Note that the saddle point Q * (D = D) together with D yields an effective noise strength which shifts average variances and covariances.Importantly, it is only the heterogeneity in the connectivity W that causes this shift.Average covariances are insensitive to heterogeneity in the noise strengths D; they only depend on the average D.

B. Variance of cross-covariances
Replica method Calculating the variances of covariances across the ensemble of possible network connectivities requires making use of the replica method [40,56] and deriving an expression for the disorder-averaged moment-generating function of the replicated system ⟨Z (J ) Z (K)⟩ W , as this allows calculating disorder averages of arbitrary squared moments , which occur in the first term in Eq. ( 23).The procedure is completely analogous to the previous section's derivations.However, the disorder average now affects the term , where X and Y refer to the activity in the first and second replicon, respectively.A cumulant expansion up to second order introduces -along four-point couplings separately in X and Y similar to the one in Eq. ( 15)a replica coupling term To deal with the four-point couplings, we again introduce auxiliary variables and obtain a relation similar to Eq. ( 20), ∫ DP ∫ D P exp (−S (P , P )) , where Q and Q are shorthand notations denoting all auxiliary variables, and (Y , Ỹ ) are given by Eq. ( 18).Saddle-point approximation As in Sec.III A, we approximate p (Q, Q) as a delta function at the saddlepoint Q * , Q * (for details see Appendix IV F), and with Eq. ( 24) to lowest order we get where we used Wick's theorem, which is allowed by the fact that, for Q and Q given and fixed, ZQ, Q (J , K) describes a Gaussian theory, and the fact that all crossreplica correlators ⟨X i Y j ⟩ Q * , Q * vanish at the saddle point (see Appendix IV F).
Fluctuations around the saddle-point Equation ( 27) implies that the variance of covariances is zero in the saddle-point approximation, and we need to account for Gaussian fluctuations of the auxiliary fields around their saddle points by making a Gaussian approximation of p (Q, Q).The crucial fluctuations are the ones of Q XY and QXY , as they can potentially preserve the replica coupling and thus lead to nonvanishing variance contributions of cross-replica correlators ⟨X i Y i ⟩ Q * , Q * .Away from the saddle points, the correlators in Eq. ( 26) depend on Q and Q in a complicated manner.To render the integrals in Eq. ( 26) solvable in the Gaussian approximation, we perform a Taylor expansion of the correlators around the saddle points Q * , Q * , which effectively is an expansion of ZQ, Q (J , K) (see Appendix IV G for more details).In the first term of Eq. ( 26), leading order fluctuations in Q XY and QXY depend on correlators with an odd number of variables of each replicon.Therefore, this term cannot yield a contribution to the variance due to fluctuations of Q XY and QXY .The major replica coupling arises from the second and third terms in Eq. ( 26).We note that the third term contains off-diagonal elements of correlators which are suppressed by a factor 1/N with respect to the diagonal ones.Therefore, we can neglect this term for cross-covariances as well and only keep the second term in Eq. ( 26) as the leading order contribution.For autocovariances the second and third terms in Eq. ( 26) are the same, yielding an additional factor 2.
where we used that cross-replica correlators vanish at the saddle point.Inserting the above fluctuation expansion result around Q * XY and Q * XY into Eq.( 26) leads to Next, we consider the Gaussian approximation of where S (2) contains the second derivatives with respect to the auxiliary fields , which allows evaluating the correlators of the auxiliary fields in Eq. ( 29) (see Appendix IV H for details).Inserting the results, to leading order we find (see Appendix IV I for details) To get the variances rather than the second moments, we subtract the squared mean covariances ⟨C ij ⟩ 2 W .However, for the setup that we study here the squared mean cross-covariances are of the order O (1/N2 ) and therefore negligible.Taking into account that R = (1 − M ) −1 ≈ 1, which holds as long as the network is inhibition dominated 2 , we find the following expression for the disorderaveraged variance of cross-covariances (see Appendix IV I for full expression) where we wrote S = ∆/N .However, if the noise strength D has to be estimated using Eq. ( 11), this expression is still dependent on the specific realization of W , both implicitly through the estimates of the single-neuron rates and CVs described in Sec.II and explicitly through the matrix B [Eq. (12)].Since the right hand side of Eq. ( 31) depends nonlinearly on D, averaging over the statistics of D introduces terms depending on the heterogeneity of D. However, Fig. 9 in the Appendix shows that heterogeneity in D -both via the explicit dependence on W and via the implicit dependence through distributed firing rates and CVs -is negligible for the statistics of cross-covariances.This can be understood by considering the structure of Eq. ( 31): The matrices (1 − S) −1 are multiplied with D, such that any heterogeneity in D is averaged out.An E-I network is an illustrative example, with (1 − S) = 1 + U with a 2 × 2 block matrix U whose entries are homogeneous in each population block, such that the matrix product effectively is an average over D.
To obtain an average D that is not depending on a specific realization of W , we follow Eq. ( 11) and set which inserted into the disorder-averaged expression for the autocovariances [Eq.( 22)] yields the correct autocovariances: Here we used (1 − M ) The realization-independent estimates ν i and CV 2 i of the rates and CVs, respectively, can be obtained using standard population-resolved mean-field theory [42,43], which only requires knowing the statistics of W .A procedure similar to the one described in Sec.II can be used: In the population view, however, the indices i, j no longer denote single neurons but rather populations of equal neurons.In Eq. (1) J ij is replaced by K ij J ij and J 2 ij in Eq. ( 2) is replaced by K ij J 2 ij , where K ij is the indegree from population j to population i, and J ij then is interpreted as the mean synaptic weight from population j to population i.
Replacing D in Eq. ( 31) by Eq. ( 32) yields a fully realization-independent disorder-averaged estimate of the variance of cross-covariances.

C. Singularities
Next, we discuss the interpretation of the derived formulas.Thereto, we need to have a closer look at the effective noise strength D + diag [Q * (D)], which occurs in both the mean [Eq.(22)] and the variances [Eq.( 31)] of covariances.Using Eq. ( 32), we find that the impact of heterogeneity on the effective noise cancels: where a i = CV 2 i ν i is the vector of estimated autocovariances.This is because we specifically choose the noise strength D such that autocovariances match those from the spiking networks: As heterogeneity is increased, external fluctuations get amplified by the factor (1 − S) −1 in Eq. (34).To achieve that autocovariances do not diverge, external inputs need to be scaled down according to Eq. (32).Hence, the mean and variance of crosscovariances are given by Note that any inverse matrix can be written as , where adj (A) denotes the adjugate matrix.As a result, the elements of an inverse matrix A −1 diverge if the determinant of the matrix A vanishes, which occurs when at least one eigenvalue of A is zero.Therefore, the divergence behavior of the mean and variance of covariances is determined by the eigenvalues of M and S with real parts close to 1.
Eq. ( 35) reveals that mean cross-covariances are determined by the mean connectivity M .By choosing D to match the autocovariances of the spiking model, they are, in particular, unaffected by network heterogeneity, represented by S. A range of important network properties, such as population structure determining E-I balance [34,36,37,45,57], spatial structure like distance-dependent connection probabilities [8,[58][59][60][61][62], or low-rank structures [63], can be encoded in M .Divergences in mean covariances, caused by eigenvalues of M close to 1, can thus be indicative of phenomena like loss of E-I balance with excessive excitation (cf.Fig. 8D of Ref. 57) or instability of the homogeneously active state in spatially organized networks [64].
Variances of cross-covariances are determined by network heterogeneity [Eq.(36)], encoded in the connectivity variance S and are to leading order independent of the mean connectivity M .Note that subleading terms nevertheless can become sizable if eigenvalues of M are close to the instability line at 1.As demonstrated by Aljadeff et al. [65], if S is a block structured matrix encoding different populations, its eigenvalue spectrum is circular, with a spectral radius r that is determined by the square root of the maximum eigenvalue of S.
For the E-I network studied here, the matrices M and S have nontrivial block structure with one excitatory and one inhibitory block (see Appendix IV M).In this case, the spectral radius is given by [66] r 2 = N E σ 2 E + N I σ 2 I .The spectral radius, a measure of network heterogeneity, increases when the variance of synaptic strength grows, which is controlled by an interplay between the connection probabilities of different populations and the variances of the associated synaptic weights.Intuitively, as explained in Dahmen et al. [8], multisynaptic signal transmission is very efficient in a network with a large spectral radius, such that pairs of neurons influence each other via a large number of neuronal pathways, possibly including differing numbers of excitatory and inhibitory neurons.The effects of these various pathways add up, and the large variety of potential pathways results in a broad distribution of covariances.
We see that the effects of M and S are mostly independent of one another, allowing the mean and variance of covariances to vary separately.This, however, applies only to synaptic weights that are identically and independently distributed.If the weights are correlated, such as through chain structures in the connectivity, the respective eigenvalues cannot be changed independently.A more detailed analysis of this behavior is to be published elsewhere.As a final remark, it is worth noting that the independence of the mean covariances of S confirms that previously employed population models [31,36,37,45,57], which neglect the variance of connectivity, are valid for computing mean covariances.
To illustrate how the mean and variance of covariances change as functions of the network heterogeneity, we plot Eq. ( 22) and Eq.(31) with Eq. (32) for spectral radii between 0 and 1 in Fig. 4 (predicted linear).We kept the working point roughly constant for the different spectral radii by maintaining the mean µ and variance σ 2 of the total input to each neuron while modifying the synaptic efficacy.To compensate for the increased intrinsic input and fluctuations at larger spectral radii, we reduced the mean and fluctuations of the external input.
Confirming the discussion of Eq. ( 35) and Eq. ( 36), when the spectral radius is modified, the variances of covariances vary by several orders of magnitude, whereas mean covariances remain in the same order of magnitude.A range of prior research [31,37,57] has shown that a divergence of mean covariances would be observed as a function of E-I balance, e.g. by altering g.Here we focus on network scenarios away from the excitatory instability (fixed g = −6) and therefore do not see a divergence of mean covariances.Nevertheless, we observe a change of mean covariances when changing the spectral radius.This is because in the sparse random network chosen here, the variance of the synaptic weights is not independent from the mean of the weights.Adjusting the spectral radius requires modifying the weights, resulting in the residual change in the mean covariances visible in Fig. 4(a-c).Note that by keeping the working point of the network constant across spectral radii, we also keep the noise strength factor in Eq. ( 22) and Eq. ( 31) constant [cf.Eq. ( 35) and Eq. ( 36)].If the external noise strength was instead determined by a fixed external process, i.e D independent of W , then mean covariances would also diverge as a function of the spectral radius due to the factor (1 − S) −1 , which enters the noise strength term via

D. Comparison of prediction and measurement of covariance statistics
To check how closely the predictions match the outcomes of spiking network simulations, we ran ten simulations for different spectral radii of the effective connectivity matrix W similar to the one shown in Fig. 1 using the parameters specified in Appendix IV A. The spectral radius was modified by changing the synaptic weight scale j (see Table II).In the sparse E-I networks that we consider here, this affects both the mean connectivity M and the variance of connections S [cf.Eq. ( 39)].We ensured that the spiking networks have roughly similar working points for the different spectral radii by adjusting the external input as explained in Appendix IV D. The theory of Eq. ( 22) and Eq. ( 31) yields, for any pair of neurons in the network, predictions for the mean and variance of covariances measured across different network realizations.In brain circuits, one finds statistical similarities between different neurons, for example, given by their population membership.These similarities result in symmetries in the matrices M and S. In the current example of the excitatory-inhibitory network, the statistical similarity among all excitatory neurons and the similarity among all inhibitory neurons is reflected in the block structure of the matrices M and S.This block structure results in a similar block structure for the matrices for the mean and variance of covariances.Any off-diagonal element of one of the blocks (EE, EI, or II) is representative of the statistics of cross-covariances for this type of neuron pair.If networks and groups of statistically similar neurons are sufficiently large, then -by the selfaveraging property [7,40,50,51] -an empirical average across these statistically similar neurons is insensitive to the particular network realization and can be compared to the results for the statistics across network realizations from the theory.We thus computed the empirical mean and variance of the measured covariances for each type of neuron pair (EE, EI, and II), corrected the variances for bias due to finite simulation time (see Appendix IV K for details), and compared the results to the predictions by Eq. ( 22) and Eq.(31).The results are displayed in Fig. 4: The top row shows changes in mean covariances with spectral radius of W [via according changes in M , see Eq. ( 22)], and the bottom row shows changes in the variance of covariances with spectral radius of W [via according changes in S; see Eq. ( 31)].
We observe that the order of magnitude of mean and variance are well predicted by Eq. ( 22) and Eq. ( 31), which is especially evident for the variances [Fig.4(d-f)], which span several orders of magnitude.However, there is some quantitative discrepancy between the predictions of the presented linear theory and the results of the simulated spiking network, which is visible in Fig. 4(a-c), indicating that a linear theory cannot fully capture the nonlinear spiking dynamics at high spectral radii, where potential nonrenewal effects of spiking arise [67].To verify that the discrepancy originates mostly from the linearresponse approximation rather than our disorder-average approximations, we plotted the predictions of the linear theory [Eq.( 9)] for 20 different network realizations: At small spectral radii, the predicted disorder-average-based mean is equal to the empirical mean of the linear networks, and for large spectral radii, the predicted mean appears to be within the range of two standard deviations around the empirical mean.This shows that the deviations to the spiking network results mostly stem from the linear-response approximation.The remaining difference between the predicted and the empirical mean in linear networks could be explained by the fact that, for high spectral radii, the effective connectivity matrix contributes much more strongly to the noise strength, such that we can no longer disregard its contribution to the noise strength (cf.Fig. 3) and averaging over W and D separately is no longer feasible.

IV. DISCUSSION
In this study, we introduce theoretical tools based on statistical physics of disordered systems to investigate the role of heterogeneous network connectivity in shaping the coordination structure in neural networks.While the presented methods are applicable to arbitrary independent connectivity statistics, for illustration we focus our analysis on the prototypical network model for cortical dynamics by Brunel [42], which is a spiking network of randomly connected excitatory and inhibitory leaky integrate-andfire neurons receiving uncorrelated external Poisson input.This model has been extensively studied before using mean-field and linear-response methods to understand neuronal spiking statistics such as average firing rates and CVs [30,68] as well as average cross-covariances between populations of neurons [27,37,45,57,58].In this study, we go beyond the population level and in-   (7)] measured in network simulations of spiking E-I networks.Solid lines show the theoretical predictions of the mean [Eq.( 22)] and variance [Eq.( 31)] of cross-covariances, respectively, using the noise strength estimate [Eq.(32)].The dashed lines show the population-resolved empirical mean and variance of cross-covariances from Eq. ( 9) averaged over 20 realizations of W , and the shaded area depicts a two-standard-deviation range around the mean.The variances computed from the simulation results have been corrected for bias due to finite simulation time (see Appendix IV K).Same excitatory-inhibitory network model as in previous figures.For model details and simulation parameters see Appendix IV A.
troduce tools from field theory of disordered systems to study the heterogeneity of activity across individual neurons.We show how to turn a linear-response result on the link between covariances and connectivity [27][28][29][30][31] into a field-theoretic problem using moment-generating functions.Then we apply disorder averages, replica theory, and beyond-mean-field approximations to obtain quantitative predictions for the mean and variance of crosscovariances that take into account the statistics of connectivity, but are independent of individual network realizations.We show that this theory can faithfully predict the statistics of cross-covariances of spiking leaky integrate-and-fire networks across the whole linearly sta-ble regime.In doing this, we fixed the statistics of individual neurons according to their theoretical prediction and showed that this one working point, defined by the firing rates of all neurons in the network, can correspond to very distinct correlations structures.Furthermore, we demonstrate that while the heterogeneity in singleneuron activities directly impacts the statistics of neuronal autocovariances, it does not have a sizable impact on the heterogeneity in cross-covariances (cf.Fig. 9).The latter heterogeneity is determined by the heterogeneity in neuronal couplings, quantified by the spectral radius of effective connectivity bulk eigenvalues.
Technically, by employing linear-response theory, we study two systems: the spiking leaky integrate-and-fire network and a network of linear rate neurons.We derive a procedure to set the external input noise of the linear model in such a way that the covariance statistics of the spiking network and the linear network match quantitatively.This way, the autocovariances are fixed to values determined by single-neuron firing rates and CVs, as predicted by renewal theory for spike trains.Consequently, autocovariances remain finite in the matched rate network even when approaching the point of linear instability.This is achieved by reducing external input fluctuations to account for the increased intrinsically generated fluctuations when increasing the heterogeneity in network connectivity.As a result, also neuronal cross-covariances remain finite close to linear instability.The variance of cross-covariances nevertheless displays a residual divergence, which is why, within the linear regime, mean cross-covariances only vary mildly, while the variance of cross-covariances spans many orders of magnitude when changing the spectral radius of bulk connectivity eigenvalues.
The methods presented here are restricted to the linearly stable network regime, usually referred to as the asynchronous irregular state of the Brunel model [42].We show that, while mean covariances are low in this state [5,34,36], individual cross-covariances between pairs of neurons can still be large, reflected by the large variance of cross-covariances in strongly heterogeneous network settings.Linear stability can, for example, be realized in excitatory-inhibitory networks if the overall recurrent feedback in the network is inhibition dominated or only marginally positive [37] and if synaptic amplitudes are not too strong.Previous work [67,69] has shown that the here-considered model transitions to a different asynchronous activity state if synaptic amplitudes become larger.This state, however, is not well described by linear-response theory, as slow network fluctuations and nontrivial spike-train autocorrelations emerge, causing deviations from the renewal assumptions on spike trains used here.Note that such slow network fluctuations have not been observed in previous studies on spontaneous activity in macaque motor cortex [7,8] and mouse visual cortex [41] that employed first results of the more general theoretical approach presented here to explain experimentally observed features, such as the large dispersion of covariances, long-range neuronal coordination, a rich repertoire of time-scales and low-dimensional activity.These studies relied on Wick's theorem to calculate the variance of covariances, which is, however, restricted to linear systems.Here we instead employ a more general replica approach that can be straightforwardly applied to nonlinear rate models [52], as extensively studied in the recent theoretical neuroscience literature [65,[70][71][72][73][74][75][76].Importantly, the replica theory reveals in a systematic manner that the variance of covariances is an observable that is O(1/N ) in the network size and requires beyond-mean-field methods to be computed.In mean-field or saddle-point approximation, the replica coupling term that yields the nontrivial variance of covariances vanishes.We here calculate the next-toleading-order Gaussian fluctuations around saddle points that yield good quantitative results across the whole linear regime.The fact that the linear rate model captures the covariance statistics of the spiking leaky integrateand-fire model further shows that the presented results on the link between connectivity and covariances do not depend on model details and are generally valid in the linear regime, which enables applications to experimental data [7,8,41].
In this paper, we focus on intrinsic mechanisms for heterogeneity and study the first-and second-order statistics of network connectivity.The formalism can be applied to any network topology, as arbitrary connectivity structures can be encoded in the mean and variance matrices that are the central objects of the theory.The results of the formalism are particularly useful for comparing covariance statistics within a single network with groups of statistically equivalent neurons sharing the same connectivity statistics, such as different cell types within neural circuits.Such statistical equivalence imposes symmetries on the structure of the matrices M and S that encode the first-and second-order connectivity statistics.These symmetries allow a dimensionality reduction of the problem and, by the self-averaging property [7,40,50,51], the comparison of theoretical results to empirical averages over covariances of different pairs of neurons within a single circuit.Here we study a network with two populations, excitatory and inhibitory, leading to a block structure in both M and S.This example has been chosen as the simplest but relevant setting that goes beyond the fully homogeneous random network in Dahmen et al. [7] with a correspondingly simpler homogeneous theory.Extensions to more populations or neuron clusters as well as more complex population-specific connectivity statistics are straight forward by adding more blocks to M and S. Another application of the here derived theory to a more complex scenario, including interneuron distance dependence of connection statistics, has been studied (without derivation) by Dahmen et al. [8].Notably, in our theory we assume that connection weights are independently drawn across different neuron pairs from an arbitrarily complex probability distribution with finite cumulants.The focus on mean M ij and variance S ij of this distribution is justified as long as connection weights scale at least as O (1/ √ N ), a scaling that is often employed to preserve fluctuations in the limit of large networks [34,55].In this scaling, effects of higher-order-connectivity cumulants are suppressed by the typically large network size, and the Gaussian approximation of the connectivity yields accurate results for network dynamics, as here demonstrated on the example of excitatory-inhibitory networks with population-specific connectivity statistics comprising sparseness (Bernoulli distribution) in addition to distributed (Gaussian) synaptic amplitudes of existing connections.Generalization of dynamic meanfield methods to heavy-tailed connectivity, which cannot be expanded in cumulants, has been proposed for studying single-neuron activity statistics [77,78].A similar approach may potentially be combined with the methods presented here to investigate cross-covariances.Furthermore, extensions to correlated connection weights, reflecting an over-or under-representation of reciprocal, convergent, divergent, and chain motifs, have been proposed in Ref. [41].
In addition to network connectivity, external inputs can be correlated and heterogeneous and thereby cause heterogeneity in covariances of local circuits.Previous works have shown that external inputs can have a strong impact on local covariances, especially in the limit of infinite network size [36,58,79,80].For biologically realistic network sizes of local circuits, intrinsically generated covariances via local recurrent activity reverberations, however, make up a substantial contribution to cross-neuronal coordination [7,57].This contribution is explainable with the here-presented methods.In general, more research is required to decipher the precise interplay between intrinsic heterogeneity and external inputs to arrive at a complete picture for the mechanistic origin of heterogeneous covariance structures in local circuits.
Many previous studies have linked connectivity and dynamics on an average level, taking into account particular connection pathways between neural populations [19], clustering [81,82], or the spatial dependence of connections [83][84][85].In contrast, we here focus on heterogeneity as a key feature of neural network connectivity and show that it yields a wealth of complex coordination patterns that are progressively becoming experimentally accessible via recent advances in measurement techniques [15].Our theoretical framework to systematically incorporate structural heterogeneity and predict dynamical heterogeneity in biologically plausible neural network models enables the use of this experimental knowledge about neural systems.Likewise, the current framework can be used for the inverse problem of inferring network properties from measured covariances, as we demonstrated in Ref. [7] for the spectral radius in homogeneous random networks (see Appendix IV M for an extension to E-I networks) and in Ref. [8] for long-range, multisynaptic interactions in networks with spatially organized connectivity.Our work thus opens new avenues for the interpretation of data on network structure and dynamics and proposes a change of focus from population-averaged observables to higher-order statistics that uncover the central role of heterogeneity in biological networks.
All code and data to reproduce the simulations and figures of this study are publicly available under the Ref. [86].

APPENDIX A. Network model and NEST simulation
We simulate networks of leaky integrate-and-fire neuron models, where the subthreshold dynamics of the membrane potential V i of neuron i is given by with total input current I i (t) that consists of recurrent input via connections with strength J ij and delay d as well as external input: Synaptic currents are instantaneous without synaptic filtering [88,89].The external input is decomposed into a constant current I ext and Poisson spike trains s ext,E (t) of rate ν ext,E and s ext,I (t) of rate ν ext,I that affect neurons with excitatory weight j and inhibitory weight gj, respectively.R and C denote the membrane resistance and capacitance, respectively.More information on the model parameters and their values can be found in Table I and Table II.We consider random sparse connectivity J , with connection probability 10 %.Sparse connections are realized with a fixed excitatory indegree K E = 800 and inhibitory indegree K I = 200, with potential self-connections and prohibiting multiple connections between the same pair of neurons.To compensate for the imbalance in excitatory and inhibitory neuron count, we scale the strengths of existing inhibitory connections with respect to excitatory ones by a factor g = −6 to obtain an asynchronous irregular dynamic regime [42].We distribute synaptic amplitudes of existing connections according to population-specific normal distributions j E ∝ N (j, 0.2j) and j I ∝ N (gj, 0.2j).The parameter j determining the overall scale of connection strength is varied to modify the overall heterogeneity of connections as measured by the variance and thereby the spectral radius of bulk connectivity eigenvalues (Table II).

B. Time-lag integrated covariances
The cross-covariance function of two stochastic zero-mean processes x i (t) and x j (t) is defined as where the average is over the ensemble of realizations of the processes.If the stochastic processes are stationary, the cross-covariance function solely depends on the time-lag Here we are considering the time-lag integrated covariances, as they can be linked to the experimentally accessible spike-count covariances [6,7], which can be interpreted as a zero-frequency Fourier transform.The Wiener-Khinchin theorem (IV C) allows expressing the time-lag integrated covariances in terms of the time series's Fourier components X i (ω) at frequency zero  II.Parameters adjusted for setting different spectral radii while keeping the firing rate constant: The spectral radius of the connectivity is set by choosing different synaptic strengths j of connections.The synaptic strength not only affects the mean connectivity but also the variance of connections, i.e., the heterogeneity of the network, that determines the spectral radius [see Eq. ( 39)].The parameters of the external inputs, which model the total excitatory and inhibitory input from external populations of neurons, are furthermore adjusted to maintain a constant firing rate across different spectral radii.This is done to isolate effects of changing spectral radii on correlations from effects of changing firing rates (see Appendix IV E and Fig. 6).For small spectral radii and thus weak recurrent input, strong external input is needed to drive the network to moderate firing rates, while for large spectral radii and strong recurrent input, only weak external input is needed for moderate firing rates.

C. Wiener-Khinchin theorem
Here in parts we follow the book by Gardiner [90].Let x(t) and y(t) be stochastic, stationary processes.Stationary means that for any n-tuple (t 1 , t 2 , . . ., t n ) of time points and any real number u the samples x (t 1 ) , . . ., x (t n ) follow the same distribution as the samples x (t 1 + u) , . . ., x (t n + u) [91].Consequently, we may define a raw correlation function as which, due to the assumption of stationarity, does not depend on the time t.The average is over the ensemble of realizations of the processes.If the Fourier transforms of x and y exist, we may calculate the ensemble average over X (ω) and Y (ω) as where we used the identity ∫ dt e −iωt = 2πδ(ω) which follows from 1 2π ∫ dω e iωt 2πδ(ω) = 1, so that 2πδ(ω) is the Fourier transform of the constant function and vice versa.Equation (41) states that the cross spectrum between two stationary processes vanishes except at those frequencies ω = −ω ′ , where it is proportional to a δ-distribution times the Fourier transform of the autocorrelation function.

D. Validity of theoretical predictions
In this section, we discuss the conditions under which the theory and simulation described in this paper yield the same results.There are several factors to consider: the limits of the theory we built upon, the limitations of the newly presented theory, and the simulation's constraints.
The estimation of covariances presented in this paper relies on the proper estimation of firing rates and CVs, for which we employ Eq. ( 3) and Eq. ( 4) [42], which rely on a diffusion approximation that substitutes spiking input by uncorrelated white noise input.Due to this approximation, these formulas have their own limitations, and they do not yield good estimates in all parameter regimes, as shown in Fig. 5 and Fig. 6(a) and (b) for a single neuron receiving Poisson input.Note that a correction due to the finite amplitude of spiking input could in principle be accounted for [92]; the diffusion approximation was shown to typically overestimate firing rates, as can be seen, for example, in Fig. 2(a).In a network context, the diffusion approximation furthermore neglects any nontrivial structure of the autocorrelation of inputs, as well as their cross-correlation structure.Apart from an offset, Fig. 2(a) also shows a clustering of inhibitory firing rates in theory.The difference between the two clusters can be traced back to the presence or absence of self-connections: Firing rates of neurons with self-connections are more in line with simulated rates.This fact, however, does not seem to generalize as networks with other spectral radii do not show any apparent clustering in firing rate predictions (Fig. 7).Because the estimates for firing rates and CVs are used to calculate the effective connectivity matrix and noise strength, a poor estimate has a direct impact on the covariance estimation.Furthermore, the quality of the rate estimates affects how closely the simulated network matches its analytical counterpart due to the way we set the parameters for the simulation: We fix the mean and variance of the single neuron input, and therefore their firing rates ν set , and adjust the external input to set the spectral radius r, which we estimate using the result of Rajan and Abbott [66] for random Bernoulli E-I networks with connection probability p.The effective weights w eff,E (ν), w eff,I (ν) are computed using Eq. ( 6).Once we simulate the network, we can measure the firing rates, extract the connectivity matrix, and compute the effective connectivity matrix realized in the simulation.Its largest eigenvalue determines the spectral radius r sim .A comparison of r set and r sim is shown in Fig. 6(c).They do not coincide perfectly, which is a direct result of the unreliable estimation of the firing rates, which are slightly overestimated by the theory [see Fig. 2(a), Fig. 6(a), and Fig. 7(a), (d), and (g)].
To make sure the simulated network is always in a linearly stable regime, we restrict our analysis to spectral radii r set ≤ 0.90, for which r sim ≤ 0.94.We estimate the noise strength D by computing the variances using Eq. ( 10), assuming that D is diagonal, and inverting Eq. ( 9) which yields Eq. (11).First of all, the equation for the variances Eq. ( 10) relies on the assumption that the spike trains are well described by renewal processes [48].Therefore, the noise strength estimate is reliable only if the spike trains are not too bursty.However, even for networks with CV ≈ 1 we observed that for large spectral radii this approach of estimating the noise strength can yield negative values for D, which has no physical interpretation.Measuring the covariances in a simulation and inverting Eq. ( 9) without restricting D to be diagonal, yields a matrix that seems to be almost diagonal, shown in Fig. 8(a).Setting the off-diagonal elements to zero and using the result to compute the covariances via Eq.( 9), however, reveals that the off-diagonal contribution cannot be neglected [Fig.8(b)], which means that the external noise sources do have to be correlated to explain the observed covariance.In cases in which the lowest eigenvalue of D full is negative, we conclude that it is not possible to find a physical linear system (positive-definite D) that explains the individual pair-wise covariances observed in the spiking network simulation with a large spectral radius.Our theoretical predictions for the mean and variance of cross-covariances, Eq. ( 22) and Eq. ( 31), based on D computed with Eq. ( 11) and its averaged analog, Eq. ( 32), nevertheless yield quantitatively matching results with respect to the spiking network simulations also in this regime (Fig. 4), because, as we show in Fig. 9, the results only depend on the average of D. The theory based on the statistics of connections is therefore found to be more robust than the theory based on individual connectivity realizations.
Finally, simulations have one major limitation: their finite simulation time, which results in a biased estimation of the covariances at the single-neuron level.As seen in Fig. 7(b), (c), (e), (f), (h), and (i), there is some variance in the simulations that is not explained by the theory.This variance is caused by the finite simulation time and vanishes for longer simulations.The relative unexplained variance is larger for small spectral radii, since the firing rates of the neurons are slightly smaller in these networks leading to poorer estimation of covariances, and overall the covariances are smaller for small spectral radii.

E. Derivation of moment generating function
As discussed in Sec.II, in absence of correlated external input and in the regime of low average covariances, covariances can be understood in linear response theory [31], where the dynamical equation of LIF neurons describes a model network of Ornstein-Uhlenbeck processes [Eq.( 5)].Grytskyy et al. [31] further showed that relation Eq. ( 9) between time-lag-integrated covariances and effective connections is independent of the particular filter kernel h(t) and whether noise is injected in the input or output of neurons.Therefore, we here for simplicity choose Gaussian white noise in the input and h(t) to be an exponential kernel with unit time constant.The stochastic differential equation becomes with generating functional [7] Z The latter can easily be interpreted in the Fourier domain due to the linearity of Eq. ( 42) and the invariance of scalar products under unitary transforms, with Fourier-transformed variables denoted by capital letters.The scalar product in the frequency domain reads XT X = ∑ j ∫ dω Xj (−ω)X j (ω).The generating functional factorizes into generating functions for each frequency ω.
As we use Eq. ( 8) to calculate the time-lag integrated covariances, we only require the zero-frequency components X (0).In the following, we therefore only discuss zero frequency and omit the frequency argument, i.e. we write X ≡ X (0) and correspondingly for sources J .After integrating over all nonzero frequencies, we obtain the generating function for zero frequency, with the single-frequency scalar product defined as XT X = ∑ j Xj Xj , integration measures ∫ DX = ∏ j ∫ ∞ −∞ dX j and ∫ D X = ∏ j 1 2πi ∫ i∞ −i∞ d Xj , and normalization prefactor λ.We introduce another source variable J so that later we can also compute correlators that include X The Gaussian integrals are solved as follows The identity matrix and the matrix of ones commute, therefore we can use det ( A B C D ) = det(AD − BC), and we get The normalization condition Z (J = 0) = 1 yields λ = |det (1 − W )|, and the generating function becomes = exp ( respectively.We obtain the time-lag-integrated covariances

F. Saddle points and correlators of activity fields
The saddle points including correlators evaluated at the saddle point Q * , Q * .To evaluate them, we need to solve the Gaussian integral in Eq. ( 17) where we added the additional source term J T X to allow for the calculation of correlators including X.We can rewrite the equation as , where the prefactor (2π) −N comes from the integration measure D X , such that Deriving the normalized moment-generating function Z (J , J ) = Z (J , J ) / Z (0, 0) twice with respect to J yields ) , which is solved by Q * = 0. Inserting this result, we find Inserting the correlators into the saddle-point equations and solving for Q * yields The saddle point of the auxiliary fields Q XX , Q XY , Q Y Y in the replica theory are determined by finding the zeros of the first derivative of the action [Eq.(25)].This yields and in a fashion analogous to the derivation above, we find Inserting the latter solution into the saddle-point equations again yields a linear self-consistency equation for Q * XY with the solution G. Fluctuations around saddle-points Here we showcase how to perform a fluctuation expansion of the correlator Other correlators follow analogously.Following the definition in Eq. ( 24), the correlator is given by . Now, we expand Z Q, Q around the saddle points Using the definition of ZQ, Q (J , K) in Eq. ( 25), its derivative is given by such that normalizing and evaluating at the saddle point and for zero sources yields The second term on the right-hand side of Eq. ( 46) vanishes at the saddle point and for J = K = 0 The derivative with respect to QXY,k can be computed analogously, with Xk Ỹk replaced by ∑ l ∆ kl X l Y l .Therefore, the first-order expansion in the replica coupling term reads which we use in Eq. ( 28).

H. Correlators of auxiliary fields
We consider the Gaussian approximation of p (Q, Q) with where S (2) contains the second derivatives with respect to the auxiliary fields Using S (2)  = ( 0 S (2)T 21

I. Disorder-averaged variance of covariances
Starting with Eq. ( 26), we find Inserting the derived expressions for the correlators into Eq.( 29) yields where we used .
Finally, we obtain the second moment to leading order and for the covariance we find However, a naive implementation of Eq. ( 49) is numerically unstable due to the diverging integrals.To proceed, we rewrite Eq. ( 49) using the following steps: We make a change of variables v ′ = s − v, w ′ = s − w, where we immediately drop the prime, yielding We switch the order of integration using ∫

M. Inference of connectivity features from covariances
Equations ( 35) and (36) relate the mean connectivity M and the connectivity's variance S to the mean covariances ⟨C⟩ W ,D and the variance of covariances ⟨δC 2 ⟩ W ,D , respectively.These relations can be inverted to infer features of M and S from the statistics of covariances.Here, we focus on S and, in particular, the spectral radius r = √ N E S E + N I S I , which substantially modulates the size of the variance of covariances, as shown in Fig. 4. In the E-I network considered in this work, S has the structure with {1} XY being the matrix of ones of dimension N X × N Y .The inverse matrix (1 − S) −1 then follows as with Focusing on the off-diagonal elements of each block of F yields quadratic equations that can be solved for SE and SI and subsequently, using Eq. ( 52) and Eq. ( 53), for the connectivity parameters S E and S I : with σ(x) = 1 for x ≥ 0 and σ(x) = −1 for x < 0, F EE = ⟨δC 2 EE ⟩ W ,D / (CV 2 ν)

2
. As shown in Fig. 12, this in particular allows the inference of the spectral radius of bulk connectivity eigenvalues, from the variance of covariances.In principle, a similar procedure could be followed to infer the mean connectivity M E and M I from the mean covariances.However, we performed a number of simplifications (1 − M ) −1 ≈ 1 in calculating the effective noise [see, e.g., Eq. (33)] to arrive at Eq. (35).Therefore, the inference of M E and M I based on measured mean covariances is expected to be less accurate and would potentially require a more careful mathematical treatment.

Figure 1 .
Figure 1.Simulation of excitatory-inhibitory (E-I) network of leaky integrate-and-fire (LIF) neurons.(a) Spike trains of a sample of 20 excitatory (E, blue) and 5 inhibitory (I, red) neurons.(b) Distributions of firing rates.(c) Distributions of coefficients of variation (CVs).(d) Distributions of variances of spike counts measured in time bins of 1s [cf.Eq. (7)].(e) Distributions of spike-count covariances for different pairings of neurons.(f) Distributions of spike-count correlation coefficients for different pairings of neurons.For model details and simulation parameters see Appendix IV A for spectral radius r = 0.49.

Figure 2 .
Figure 2. Simulation results and theoretical estimates for E-I network of LIF neurons.Here ρ denotes the Pearson correlation coefficient.(a) Firing rates ν.(b) Covariances C. (c) Correlation coefficients κ.For model details and simulation parameters see Appendix IV A for spectral radius r = 0.49.

Figure 3 . 2 XT
Figure 3.Effect of averaging noise strength over disorder in W and effect of applying B −1 on variances.(a) Noise strength computed using the procedure described above (noise strength in the following) vs. noise strength computed using a new realization of W (dark gray), and noise strength vs. variances (light gray).(b) Noise strength vs. noise strength computed using an average over 100 realizations of W (dark gray) and noise strength vs. variances (light gray).Same excitatoryinhibitory network model as in previous figures.For model details and simulation parameters see Appendix IV A for spectral radius r = 0.49.

Figure 4 .
Figure 4. Statistics of covariances at different spectral radii of the effective connectivity matrix W . Covariance statistics are resolved according to population membership of neuron pairs (EE, both neurons excitatory; EI, one neuron excitatory, one neuron inhibitory; II, both neurons inhibitory).Dots show the population-resolved [(a)-(c)] mean or [(d)-(f)] variance of the spike count cross-covariances [Eq.(7)] measured in network simulations of spiking E-I networks.Solid lines show the theoretical predictions of the mean [Eq.(22)] and variance [Eq.(31)] of cross-covariances, respectively, using the noise strength estimate [Eq.(32)].The dashed lines show the population-resolved empirical mean and variance of cross-covariances from Eq. (9) averaged over 20 realizations of W , and the shaded area depicts a two-standard-deviation range around the mean.The variances computed from the simulation results have been corrected for bias due to finite simulation time (see Appendix IV K).Same excitatory-inhibitory network model as in previous figures.For model details and simulation parameters see Appendix IV A.

Figure 5 .Figure 6 .
Figure 5. Validity of firing rate and CV prediction for single LIF neuron with instantaneous synapses provided with two independent Poisson inputs.(a) Simulated rates and (b) simulated CVs at given mean input µ and input variance σ 2 .(c) Relative error ϵ = |νsim − ν thy | /νsim of rate and (d) relative error of CV prediction using Eq.(3) and Eq.(4).Black pixels denote error values larger than 1.00.

Figure 8 .
Figure 8. Noise strength properties.(a) First 100 entries of D full computed from simulated covariances.(b) Comparison of covariances computed using D full and using D diag = diag (D full ).Same excitatory-inhibitory network model as in previous figures.For model details and simulation parameters see Appendix IV A for spectral radius r = 0.9.

ds e x 2 Figure 12 .
Figure 12.Set spectral radius rset (defined in Appendix IV D) vs inferred spectral radius r inf .Same excitatory-inhibitory network model as in previous figures.For model details and simulation parameters see Appendix IV A.