The spectrum of covariance matrices of randomly connected recurrent neuronal networks with linear dynamics.

A key question in theoretical neuroscience is the relation between the connectivity structure and the collective dynamics of a network of neurons. Here we study the connectivity-dynamics relation as reflected in the distribution of eigenvalues of the covariance matrix of the dynamic fluctuations of the neuronal activities, which is closely related to the network dynamics' Principal Component Analysis (PCA) and the associated effective dimensionality. We consider the spontaneous fluctuations around a steady state in a randomly connected recurrent network of stochastic neurons. An exact analytical expression for the covariance eigenvalue distribution in the large-network limit can be obtained using results from random matrices. The distribution has a finitely supported smooth bulk spectrum and exhibits an approximate power-law tail for coupling matrices near the critical edge. We generalize the results to include second-order connectivity motifs and discuss extensions to excitatory-inhibitory networks. The theoretical results are compared with those from finite-size networks and the effects of temporal and spatial sampling are studied. Preliminary application to whole-brain imaging data is presented. Using simple connectivity models, our work provides theoretical predictions for the covariance spectrum, a fundamental property of recurrent neuronal dynamics, that can be compared with experimental data.

1 Introduction

Collective dynamics in networked systems are of great interest, with numerous applications in many fields, including neuroscience, spin glasses, social and ecological networks [1]. Many studies of neuronal networks have focused on how certain statistics of dynamics depend on the network’s connectivity structure [2-4], including the population average [5] and variance [6] of pairwise correlations. These dynamic features can be estimated experimentally by repeatedly recording a small number of neurons at each time. In this sense, they may be regarded as local or marginal features of dynamics. On the other hand, certain global or joint dynamic features may only be possible or efficient to estimate by recording a population of neurons simultaneously. An important example is the eigenvalues of the covariance matrix, which are complicated nonlinear functions of all the matrix elements. These eigenvalues arise naturally when performing the widely used Principal Component Analysis (PCA) of population activity, where they correspond to the amount of variance contained in each principal component of the activity. Although one can in principle fill out the covariance matrix through repeated pairwise recordings, the matrix is much more efficiently calculated from simultaneously recorded data. Furthermore, a sample covariance matrix calculated from simultaneously recorded data also requires particular methods to address the effect due to the finite recording length (Section 3.7) which would be different from repeated recordings. Another example of global dynamic features that has recently received substantial interest [7-11] is the effective dimensionality of neural population activity. When describing the data distribution in terms of linear subspaces, the dimensionality can be defined based on the moments of the covariance eigenvalues. Many recent experimental studies have observed a low dimensional dynamics of neurons in the brain [12, 13], and theoretical investigations have illustrated the importance of having a low dimensionality for brain function and computation [14], such as when representing stimuli [15] and generating motor outputs [13].

As the experimental techniques of measuring the activity of large population of neurons in biological networks become increasingly available, new opportunities arise for studying how the network’s connectivity structure affects these global aspects of population dynamics. In this work, we study the eigenvalue distribution (i.e., spectrum) of the covariance matrix of spontaneous activity in a large recurrent network of stochastic neurons with random connectivity. We focus on several basic and widely used models of random connectivity, including independent and identical Gaussian distributed connectivity [2] (Section 3.1), networks with second-order connectivity motifs [5, 16–20] (Section 3.3), and random Excitation-Inhibition (EI) networks (Section 3.5). Random connectivity has been a fundamental model in theoretical studies of neuronal network dynamics [2, 8, 21]. It can be motivated as a minimal model to capture the complex, disordered connections observed in many neuronal circuits, such as in the cortex. Some aspects of these covariance spectra might be distinct from those under ordered, deterministic connectivity (Section 4.2).

The dynamics considered here is simple where the activity fluctuations around the steady-state are described by a linear response [22, 23], which experimentally is related to spontaneous or persistent neural activity in absence of structured spatial-temporal stimuli. Despite the simple dynamics and minimal connectivity model, we find the resulting spectrum has a continuous bulk of nontrivial shape exhibiting interesting features such as a power-law long tail of large eigenvalues (Section 3.2), and strong effects due to the non-normality of the connectivity matrix (Section E.2 in S1 Text). These covariance spectra highlight interesting population-level structures of neuronal variability shaped by recurrent interactions that were previously unexplored.

Using the theory of the covariance spectrum, we derive closed-form expressions for the effective dimensionality (previously known for the simple random i.i.d. Gaussian connectivity [6]) We show that the continuous bulk spectrum has the advantage over low-order statistics such as the dimensionality thanks to its robustness to low rank perturbations (Section 3.3 to 3.5 and 3.8). Our analytically derived eigenvalue distributions can be readily compared to real activity data of recurrent neural circuits or simulations of more sophisticated computational models. We provide ready-to-use code to facilitate such applications (see Data Availability Statement). An example of such an application for a whole-brain calcium imaging data is presented in Section 3.8.

2 Model

2.1 Neuronal networks with random recurrent connectivity

We consider a recurrent network of linear rate neurons driven by noise

Here x(t) is the firing rate of neuron i. J describes the recurrent interaction from neuron j to i. τ is a time constant describing how quickly the firing rates changes in response to inputs. The network is driven by independent Gaussian white noise ξ(t) with variance σ2, that is, the expectation 〈ξ(t)ξ(t + τ)〉 = σ2δδ(τ).

We focus on the structure of long time scale covariation in the network, which are described by the long time window covariance . C is the covariance of the summed activity over a window of ΔT: C = 〈Δs(t)Δs(t)〉, . For biophysical neurons, C typically settles to its limiting value when ΔT > 50ms [24]. It can be shown [25] that the long time window covariance C (also the zero-frequency covariance, see Section 3.6) is

Here I is the identity matrix, and A−1, A are the matrix inverse and transpose (A− = (A−1)). For simplicity we will set σ2 = 1 unless stated otherwise. The covariance matrix C can also be estimated from experimental data consisting of simultaneously recorded neurons (Methods). We consider generalizations beyond the long time window covariance in Section 3.6.

Our analysis and results start from the covariance-connectivity relation Eq (2), which also describes, or closely approximates, the network dynamics in other models (Section 5.2) including networks of integrate-and-fire or inhomogeneous Poisson neurons [23, 26–28], fixed point activity averaged over whitened inputs, and structural equation modeling in statistics [29]. These models provide additional motivations for the covariance (Eq (2)), which may allow interpreting our results in experiments where the neural activity is driven by stimuli [10].

For many biological neural networks, such as cortical local circuits, the recurrent connectivity is complex and disordered. Random connectivity is a widely used minimal model to gain theoretical insights on the dynamics of neuronal networks [2, 4]. We first consider a random connectivity where are drawn as independent and identically distributed (i.i.d.) Gaussian variables with zero mean and variance g2/N (referred subsequently as the i.i.d. Gaussian connectivity). The covariance spectrum follows directly from results in random matrices [30, 31]. We then show how to generalize to other types of random connectivity, including: those with connectivity motifs (Section 3.3), Erdős-Rényi random connectivity, networks with excitation and inhibition (Section 3.5). The theory we derived assumes the network is large and is exact as N → ∞, and we verify their applicability to finite-size networks numerically. A list of variable notations is given in Table 1 for ease of reference.

Table 1

List of notations.

Notation	Description
N	number of neurons
C	covariance matrix, Eq (2)
pC(x)	pdf of eigenvalues of C, Section 2.2
x ±	support edges of pC(x)
J	matrix of connection weights Section 2.1
σ 2	variance of white noise input, Section 2.1
g	connection strength var(Jij) = g2/N, (3)
g c	maximum g constrained by stability
g r	g/gc
κ	ρ(Jij, Jji) reciprocal motif cumulant, Eq (13)
μ	mean of eigenvalues 1N∑i=1Nλ_i
D	dimension, Eq (4)
M	number of time samples, Eq (37)
α	N/M, Section 3.7
N s	number of sampled neurons, Section 3.7
f	Ns/N, Section 3.7

2.2 Covariance eigenvalues and dimensionality

Principal Component Analysis (PCA) is a widely used analysis of population dynamics, where the activity is decomposed along orthogonal patterns or Principal Components (PCs). The PCs are the eigenvectors of the covariance matrix C (Eq (37)), and the associated eigenvalues λ are nonnegative and show the amount of activity or variance distributed along the modes. In this work, we focus on the distribution of these covariance eigenvalues, described by the (empirical) probability density function (pdf) p(x) which is defined through the equality for all a, b. We also refer to p(x) as the spectrum (which should not be confused with the frequency spectrum obtained via Fourier transform). We will derive the limit of p(x) as N → ∞ and study how it depends on the connectivity parameters such as g = Nvar(J).

The shape of p(x) can provide important theoretical insights on interpreting PCA. For example, it can be used to separate outlying eigenvalues corresponding to low dimensional externally driven signals from small eigenvalues corresponding to fluctuations amplified by recurrent connectivity interactions [32] (Section 3.8). the spectrum is also closely related to the effective dimension of the population activity. In many cases, the linear span of the activity fluctuations is full rank, N. Nevertheless, most of the variability is embedded in a much lower dimensional subspace. A useful measure of the effective dimension, known as the participation ratio [8, 33] is given by which can be calculated from the first two moments of p(x). We will also derive explicit expressions for D in random connectivity models.

3 Results

3.1 Continuous bulk spectrum with finite support

For networks with i.i.d. Gaussian connectivity (Section 2.1), there is one parameter g describing the overall connection strength. For stability of the fixed point and the validity of the linear response theory around it, g is required to be less than 1 [2]. The parameter σ in Eq (2) just scales all λ and thus is hereafter set to 1 for simplicity. Our main theoretical result is the following expression for the probability density function (pdf) of the covariance eigenvalues in the large N limit (Section A in S1 Text), where and p(x) = 0 for x > x+ and x < x−. The distribution has a smooth, unimodal shape and is skewed towards the left (Fig 1C). Near both support edges, the density scales as (Section A.1 in S1 Text).

Fig 1

Covariance spectrum under random Gaussian connectivity.

A. Compare theory (Eq (5)) with finite-size network covariance using Eq (2) at N = 100, g = 0.5. The histogram of eigenvalues is a single realization of the random connectivity. B. Same as A. at N = 400. C. Covariance eigenvalue distribution at various value of g. As g increases the distribution develops a long tail of large eigenvalues. D. Dimension (normalized by network size) vs g. The dots and error bars are mean and sd over repeated trials from finite-size networks (Eq (2) and use Eq (4)). Note some error bars are smaller than the dots E. Covariance eigenvalues vs. their rank (in descending order). The circles are covariance eigenvalues from a single realization of the random connectivity with N = 100 (Eq (2)). The crosses are predictions based on the theoretical pdf (Eq (5)). F. Same as E. but for g = 0.9 and on the log-log scale. The red dashed line is the power law with exponent −3/2 derived from Eq (5), see Section 3.2.

The above result for the distribution p(x) follows from the derivation of the circular law distribution of the eigenvalues of the random matrix J [30, 31, 34, 35]. However, to the best of our knowledge, this is the first exposition of the explicit expression for the spectrum of C, (Eq (5), which is essential for fitting to empirical data Section 3.8) and for the study of network dynamics. We emphasize that p(x) does not have a simple relation to the spectrum of J because J is a non-normal matrix (i.e., JJ ≠ JJ). This point is further elaborated in Section 3.3.2. Although the above result is derived in the large N limit, it matches accurately the spectrum of C in networks of sizes of several hundred, as demonstrated in our numerical results, Fig 1A and 1B (see also Fig A in S1 Text for additional realizations and trial-averages). In PCA and other analyses, the covariance eigenvalues are plotted in descending order vs. their rank [10, 36]. We can use the theoretical pdf Eq (5) to predict this rank plot by numerically solving the inverse cumulative distribution function (cdf), i.e., quantile function, at probability . The closed form pdf (Eq (5)) allows for using the highly efficient Newton’s method to compute the quantiles. Fig 1E and 1F show a good agreement between the theory Eq (5) and a single realization of a N = 100 random network.

3.2 Long tail of large eigenvalues near the critical coupling

As g approaches the critical value of 1, the upper limit of the support x+ diverges as (1 − g2)−3 (Section 5.3 in Methods). This corresponds to an activity PC with diverging variance and is consistent with the stability requirement of g < 1. Note that the lower edge x− is always bounded away from 0 and has a limit of as g → 1. Analyzing the shape of p(x) for large x in the critical regime g → 1 yields a long tail of large eigenvalues, following a power law (Fig 2A and 2B, Methods)

Fig 2

Approximate power-law tail.

A. The exact pdf (solid line) of the covariance spectrum compared with the power-law approximation (dashed line, Eq (7)) at g = 0.7. Inset shows the log-log scale. B. Same as A. for g = 0.8. The approximation improves as g approaches the critical value 1. C. The log error between the exact pdf and approximation as a function of g and “distance” from the support edges. We quantify this “distance” as the minimum ratio of x/x− and (more details and motivations in Section A.2 in S1 Text. The plot shows the log error is small when this ratio is large, which means x being far away from the edges. The dashed line shows the attainable region of the ratio which increases with g.

As shown in Fig 2A, the power-law shape of p(x) is apparent for g = 0.7 and does not require a g very close to 1.

To better elucidate the range of validity of the above power law, we consider the regime where 1 − g2 ≪ 1 and x ≫ 1. Define where x+ ∝ (1 − g2)−3 is the upper edge of the support of p(x). Then, where F(z) = (1 + z)1/3 − (1 − z)1/3. Thus, far from the spectrum upper edge, z → 1 and we obtain Eq (7), whereas near the upper edge z → 0 and , which is the expected square-root singularity near the edge.

The power-law approximation of the probability density function Eq (5) translates to an approximation for the cumulative distribution function This also means a power law in the rank plot has an exponent of −3/2 when connection strength g is close to the critical value (Fig 1F), providing an alternative mechanism based on recurrent circuits for the experimental observations of power-law distributed covariance eigenvalues in the literature [10, 36] (see also the model Eq 31 and discussions in Methods).

Because the probability density is small in the power-law tail, large eigenvalues can appear to be sparsely located (Fig 1A) and potentially mistaken for statistical outliers. This underscores the importance of knowing the exact distribution and support edges for interpreting PCA results of population activity, topics which we revisit later (Section 3.8). Note that a long tail in the spectrum is a distinct feature of correlations arising from the recurrent network dynamics (see also a heuristic explanation in Section E.3 in S1 Text). For example, for the Marchenko–Pastur law that is often used for modeling empirical covariance spectra, the upper edge of its support relative to the mean is bounded by 4 (Methods). In contrast, the same ratio for p(x) (Eq (5)) can be arbitrarily large as O((1 − g2)−2) (see below for calculating the mean). This highlights the difference between covariance generated by finite samples of noise and correlations generated by the recurrent dynamics.

The long tail of the eigenvalue distribution is also reflected by a low effective dimension (Eq 4). In Section B in S1 Text we show that the mean and second moment of the eigenvalue distribution above, are given by which yields for the dimension

In particular, the relative dimension with respect to the network size D/N vanishes as g approaches 1 (Fig 1D). In comparison, D/N for the Marchenko–Pastur law (Eq 35) is at least .

While these low-order moments can be derived from previous methods (see [6] and Section B in S1 Text, and also a parallel work [37]), our method allows for the derivation of higher-order moments, such as, and in general,

3.3 Impact of the asymmetry in connectivity

We next consider generalizations of random connectivity beyond the i.i.d. Gaussian model (Sction 2.1). An important feature of biological neural networks is the presence of motif structures at various scales [16–18, 38], which correspond to an overabundance of certain subgraphs, relative to their frequency in an edge-shuffled network (i.e., an i.i.d. random graph with matching connection probability). Using a random connectivity model with Gaussian distributed entries [38] Section C.1 in S1 Text, we can study the effects of the four types of second order motifs (i.e., consisting of two edges). Among them, the diverging, converging, and chain motifs correspond to a low-rank component L of the connectivity , where the remaining part contains only reciprocal motifs. Based on the results in Section 3.4, we can focus below on studying the covariance spectrum under the simpler connectivity with only reciprocal motifs, because the bulk spectrum of is remains the same when adding the low-rank L (see examples in Section 3.4 and Section C in S1 Text). Note that the magnitude of diverging, converging, and chain motifs in the original connectivity still indirectly affect the bulk covariance spectrum by affecting the parameters g and κ (see below) of (see Section C.1 in S1 Text for the details of this relation).

For ease of notation, we still use J to represent a random connectivity matrix with potentially reciprocal motifs ( described above). This is equivalent to varying the degree of asymmetry of J [31]. In this case, each component of J is but J and J are correlated, with −1 ≤ κ ≤ 1. All other correlations are zero.

3.3.1 Symmetric and anti-symmetric random networks

First, we consider two extreme cases for the reciprocal motifs: κ = 1 corresponding to J = J, and κ = −1 corresponding to anti-symmetric matrix (or skew-symmetric J = −J). These cases are much simpler to analyze, because J is a normal matrix so p(x) can be derived from the well known eigenvalue distribution of J [31]. For symmetric random connectivity,

Here stability requires that . For anti-symmetric random connectivity,

Here the network is stable for all g. The derivations are given in the Section D in S1 Text.

Since the general shape and trend depending on g of the spectrum in the simpler symmetric case (Fig 3A) is qualitatively similar to the i.i.d. case (Fig 1C), one can use it to gain intuition, for example, of the thin tail of large eigenvalues as g approaches its critical value.

Fig 3

Covariance spectrum for the symmetric and anti-symmetric random connectivity.

A. The pdf of a covariance spectrum with random symmetric J with different g (note for stability). B. Same as A., but for random anti-symmetric J = −J. The pdf diverges at x = 1 as .

Quantitatively using the equations above, we see that p(x) of the symmetric random network (Fig 3A) has a power-law tail analogous to Eq (7) as g → 1/2 (i.e., large x) but with a different exponent from the i.i.d. case (Eq 7),

The p(x) of the anti-symmetric random network (Fig 3B) does not have a long tail as the upper limit of the support is always 1. Since J here is a normal matrix, this qualitative difference from the i.i.d. random connectivity can be understood considering the eigenvalues of J [31]. The eigenvalues of J all lie on the imaginary axis and therefore never approach 1. Thus, the eigenvalues of the covariance matrix do not develop a long tail when increasing g.

3.3.2 Connectivity with general asymmetry

For the Gaussian random connectivity with κ = ρ(J, J), −1 < κ < 1, we have derived an implicit equation for p(x) in the large N limit based on the results in [31] (Section E in S1 Text). Although a closed-form expression can be derived using the root formula for quartic equations, it seems quite cumbersome, hence we show here the numerical solutions of this equation. For a fixed g, as κ increases, the distribution broadens on both sides (Fig 4C). Intuitively, these effects (also Fig 4B) can be understood due to the change of the critical g for stability, which is now given by g = (1 + κ)−1 (based on the spectrum of J [31]). As κ increases, the relative coupling strengthg = g/g = g(1 + κ), which is also the maximum real part of J’s eigenvalues [31] increases, and the shape of the spectrum changes similar to increasing g in the case of i.i.d. connectivity (Fig 1C). This motivates us to further compare the distributions p(x) with the same g to study effects due to κ beyond changing g. As shown in Fig 4D, when fixing the relative strength g the distribution narrows as κ increases.

Fig 4

Impact of reciprocal motifs.

A. Compare theoretical covariance spectrum for random connectivity with reciprocal motifs and a finite-size network covariance using Eq (2)(g = 0.4, κ = 0.4, N = 400). B. The impact of reciprocal motifs on dimension for various g = g/g (Eq (18)). For small g, the dimension increases sharply with κ. C. The spectra at various κ while fixing g = 0.4. The black dashed line is the i.i.d. random connectivity (κ = 0). D. Same as C. but fixing relative g = 0.4 to control the main effect (see text). The changes in shape are now smaller and the support narrows with increasing κ.

Consistent with the above results, we have shown (Section E.1 in S1 Text) that for all intermediate values of −1 < κ < 1, the critical covariance spectrum has an asymptotic power-law tail with the same exponent as the i.i.d. random case (Eq (7))

In comparison, the κ = ±1 are singular cases in Eq 17 and have a different limiting power-law behavior (an exponent of and no long tail). This is intuitively consistent with the spectrum of J which becomes confined on a line for κ = ±1 rather than in an ellipses for −1 < κ < 1 [31] (see Section E.1 in S1 Text for more discussion).

The shape changes of p(x) with reciprocal motifs are also reflected by the dimension measure, for which we derived a closed-form expression (Section E in S1 Text)

Here μ1 is the mean of the distribution. Comparing with Eq 10, this shows the nontrivial dependence of dimension on the reciprocal motif strength κ. As g → g = (1 + κ)−1, . The numerator of μ1 is at least 2θ − 1 ≥ 2/(1 + κ) − 1 > 0. Therefore μ1 diverges as . Since 1 + 4(g2 − θ) ≥ 1 + 4(g2 − g) ≥ (1 − 2/(1 + κ))2 > 0, we have D/N vanishes as . The above limits under the critical g are similar to the κ = 0 case, Eq (10). Consistent with the shape changes, the dimension decreases with reciprocal motifs when fixing g and increases when fixing g (Fig 4B).

The general asymmetric random connectivity also provides an example of the strong effect of J being a non-normal matrix on the covariance spectrum. By continuity, one may expect that as κ decreases towards −1, the shape of p(x) will become similar to that of the anti-symmetric network p(x), which is bimodal for sufficiently large g (i.e., has another peak in addition to the divergence at 1, Fig 3B). Indeed, assuming a normal J predicts a covariance spectrum that is bimodal with a non-smooth peak in a large region of −1 < κ < 0 and g (Fig H in S1 Text). Intriguingly, the actual spectrum p(x) is unimodal for all but a minuscule region of (κ, g) where κ < −0.95, indicating a suppression of a bimodal spectrum under negative κ due to the non-normal J (Section E.2 in S1 Text).

3.4 Adding low rank connectivity structure

An important property of the spectrum of C is the robustness of its bulk component to the addition of low rank structured connectivity. Many connectivity structures that are important to the dynamics and function of a recurrent neuronal network can be described by a full rank random component plus a low rank component [39, 40]. For example, such components may arise from Hebbian learning [41] and by training neural networks by gradient decent [42]. A simple case is where we add a rank k structured matrix that is deterministic or independent to the random component [43, 44]. As shown in Section C in S1 Text, in large networks, the bulk covariance spectrum remains unchanged, but the low rank component may give rise to at most 2k outlying eigenvalues. This is illustrated by the example of rank-1 perturbation to J with i.i.d. Gaussian entries in Fig 5C and 5D, where the expected location of the outliers in the covariance spectrum can be predicted analytically (Fig 5E and 5F and Section C.3 in S1 Text). This is in contrast to the spectrum of J, where the same perturbations can lead to an unbounded number of randomly located eigenvalues [43, 45] (Fig 5A and 5B). In sum, the bulk spectrum of covariance is robust against low rank perturbations to the connectivity. Note, however, the relevance of the bulk spectrum for the network dynamics depends on the location of outliers. Outliers to the right of the bulk spectrum may indicate potential instability of the dynamics even for g < 1, as discussed in the example below.

Fig 5

Robustness of the covariance spectrum to low rank perturbations of the connectivity.

A. Eigenvalues of a Gaussian random connectivity J (Eq (3)), g = 0.4, N = 400. As N → ∞, the limiting distribution of eigenvalues is uniform in the circle with radius g ([34] red solid line). The black dashed line is the 0.995 quantile of the eigenvalue radius calculated from 1000 realizations. B. Same as A. but for the rank-1 perturbed J + xuv. , and x = 4.03. This example also corresponds to adding diverging motifs (Section 3.3 and Section C.1 in S1 Text). C. The histogram of covariance eigenvalues (Eq (2)) under the J in A. D. The bulk histogram of eigenvalues with J + xuv has little change and remains well described by the Gaussian connectivity theory (red line, Eq (5)). Besides the bulk, there are two outlier eigenvalues to the left and right (inset, arrows) E,F, Analytical predictions (solid and dashed lines) of the outlier locations given g and |x| when u, v are (asymptotically) orthogonal unit vectors that are independent of J (see other cases in Section C.3 in S1 Text). The y-axis is the outlier location subtracting the corresponding edge x±, Eq (6), so it is zero for small |x| before the outlier emerges (dashed line). The dots are the mean of the smallest (for the left outlier) or largest (right outlier) eigenvalues averaged across 100 realizations of the random J, N = 4000. The errorbars are the standard error of the mean (SEM, many are smaller than the dots).

3.5 Sparse excitatory-inhibitory networks

The Gaussian random connectivity has a non-zero connection weight for all pairs of neurons with probability 1, where many biological networks are sparsely connected. In addition, each neuron has both excitatory (positive) and inhibitory (negative) weights, in contrast to many neuronal networks that obey Dale’s Law, namely all neurons are either excitatory (with all outgoing weights positive) or inhibitory (with negative outgoing weights). We consider here a simple model of E-I network, consisting of N/2 excitatory and N/2 inhibitory neurons. The probability of each connection J to be nonzero, which may depend on the types of neurons i and j, is K/N, α, β = E, I. Thus, the mean number of inputs to a neuron of type α from a population of type β is K/2. All excitatory non-zero connections are of strength and the inhibitory connections are . We assume that K = kK where k = O(1) and K ≪ N.

To map this architecture on to the one studied above, we adopt the framework of [21] and consider the equivalent Gaussian connectivity with matching variance for each J. Importantly, the choice of connection probabilities and weights ensures that var(J) is (to the leading order for N ≫ 1) regardless of the cell type of neuron i, j. This allows us to define the effective synaptic gain as for all neurons. The mean of the connections between a presynaptic neuron of type β and postsynaptic α is where w = w0 and w = −w0. Thus, we can write J as a zero-mean i.i.d. Gaussian matrix with uniform variance and a rank-2 matrix of the means. As stated in Section 3.4, in such a case the bulk spectrum of the neurons’ covariance matrix is the same as Eq (5). In addition there are at most 4 outlier eigenvalues. For K ≫ 1, from the analysis of [21], the stability of the recurrent dynamics of a linear network with the above connectivity amounts to the requirement that all eigenvalues of the 2 × 2 matrix have negative real parts. Fulfilling this condition by choosing appropriate values for k (see example in Fig 6A and Section C.3 in S1 Text) guarantees that the outlier(s) due to the nonzero means are to the left of the bulk covariance spectrum so that the largest eigenvalue is x+(g), Eq (6). For K = O(1), the results in [43] show that the above condition is sufficient but not necessary for stability. For example, when all k are equal to k, which corresponds to a balance of excitation and inhibition [45], all eigenvalues of M are 0 and the dynamics is stable for g < 1 for large N. At the same time, there can be two outlying eigenvalues on the two sides of the bulk covariance spectrum (Fig 6B), whose expected location can be predicted (Section 3.4 and Section C.3 in S1 Text). Several additional examples including all inhibitory networks are considered in the Section F in S1 Text.

Fig 6

EI networks.

A. One realization of the covariance eigenvalues by Eq (2) with an EI network satisfying the stability condition (see text). The bulk spectrum is well described by the Gaussian random connectivity theory (solid line, Eq (5)). There is one small outlier to the left of the bulk (arrow). The parameters are g = w0 = 0.4, , , , , K = 60, N = 1000. To improve the accuracy of the theory to finite K, N, here we use a slightly modified connection weight, , for all excitatory non-zero connections, and similarly for inhibitory connections, such that holds exactly for finite N. B. Similar as A but for balanced EI network (see text) with k = k = 1, g = 0.4, K = 40, N = 400. Note there are two outliers on both sides of the bulk.

3.6 Frequency dependent covariance

We have so far focused on the long time window covariance matrix. This would be especially suitable for neural activity recordings with limited temporal resolution such as calcium imaging [46]. Temporal structures of correlation beyond the slow time scale can be described by the frequency covariance matrix (or coherence matrix) where is the Fourier transform of the neural activity and z† is the complex conjugate. C(ω) can also be calculated by the Fourier transform of the time-lagged cross-correlation functions C(τ) = 〈x(t)x(t − τ)〉 (Wiener-Khinchin theorem). Analogous to Eq (2)C(ω) obeys [25],

Here z† is the complex conjugate and |z| is the norm for a complex z. The transfer function a(ω) = (1 + iτω)−1 summarizes the dynamics of single neurons in the network and corresponds to a response filter of e−/τ, t > 0 for the model of Eq (1) (see also Section 5.2). The long time window covariance we have studied corresponds to C(ω = 0) (Wiener-Khinchin theorem).

For the i.i.d. Gaussian random connectivity J, we can show that the spectrum of C(ω) is given by the same Eq (5) for C(0) (up to a constant scaling) by replacing g with a frequency dependent g(ω) (compare with Eq (3))

We can use Eq (21) to study the scaling of frequency as g approaches the critical value of 1. In many cases, we can expect that the neuronal and synaptic dynamics lead to a smooth effective low-pass filtering of the recurrent input, such that for small frequency g(0) − g(ω) ∝ ω2. For the specific g(ω) in Eq (21), this can be directly verified. The low-pass filtering implies that the frequencies showing a critical covariance spectrum are those with .

Note, however, the simple replacement by an effective g may not apply to a connectivity matrix that does not have i.i.d. entries. For example, for networks with non-zero reciprocal motifs, the covariance spectrum changes qualitatively with frequency (Section D.1 in S1 Text).

3.7 Sampling in time and space

The theoretical spectra we have discussed are based on the exact covariance matrix (Eq (2)). For neural data, this is equivalent to the limit of the sample covariance (Eq (37)) when the number of time samples M is much larger than the number of neurons N. Note that if the activity data is first averaged or summed over a time window/bin (ΔT in Eq (37)) before calculating the sample covariance, then M is the number of bins. However, many large-scale neural recordings are in the so-called high dimensional regime, where N and M are comparable, that is, the ratio α = N/M remains finite or even greater than 1 for large N, M. It is thus important to study this effect of temporal sampling on the covariance eigenvalues to better relate to experimental data [7].

We refer to and as the time-sampled covariance and spectrum. The relation between the original spectrum p(x) and the time-sampled spectrum for a finite α ≥ 0 has been studied in [47]. The authors derived a general relation between the generating function of the eigenvalue distribution , where μ is the n-th moments of the eigenvalue distribution, and the counterpart for the sampled distribution,

We give an alternative derivation of this result using free probability [48-50] (Section H in S1 Text), which allows us to also generalize to the spatial sampling case (see below). For simplicity, here we describe the results for 0 ≤ α ≤ 1. For α > 1 where time samples are severely limited, the spectrum of the N − M nonzero eigenvalues can be calculated with small modifications (Section H.2 in S1 Text).

One corollary of Eq (22) is a simple formula for how the (relative) dimension changes under time sampling

These formulas show that both D and decrease with α (fewer time samples).

The relations Eqs (22) to (23) apply to any covariance matrix spectrum. For example, it reproduces the time-sampled dimension derived in [7] for a different model of covariance C. We now apply Eq 22 to the case of i.i.d. Gaussian connectivity to derive specific results of the time-sampled spectrum . Here Eq (22) becomes a cubic equation and can be solved analytically (Section H.3 in S1 Text). Consistent with the dimension, when α increases from 0 to 1, the support of the time-sampled distribution expands from both sides (Fig 7A). In particular, for any fixed α < 1 (so is positive definite), the left edge of the support x− decreases with g but is always bounded away from 0 even as g → 1, where (see also Fig L in S1 Text). Interestingly, the approximate power law of p(x) (Eq 7) still holds under time sampling for any fixed α as g → 1 (Section H.3 in S1 Text).

Fig 7

Effects of sampling in time and space on the covariance spectrum.

A. For the i.i.d. Gaussian random connectivity, how different levels of time samples α change the spectrum (Eq (22)). The non-sampled case corresponds to α = 0. g is fixed at 0.4. Inset: The relative dimension vs. α (Eq 23). The dots correspond to the pdfs with matched colors. B. Same as A. but for the spatial subsampling (Section I.2 in S1 Text), at g = 0.5. The non-sampled case corresponds to f = 1. The relative dimension in the inset is based on Section I.1 in S1 Text.

Another challenge for fitting to neural data is that often only a subset of neurons are observed instead of the entire local recurrent network. The unobserved neurons have an impact on the dynamics and affect the eigenvalues of the observed covariance matrix. We study this by considering randomly selecting N = fN, 0 < f ≤ 1 neurons and define their covariance as the space-sampled covariance. Using the free probability approach, we derive similar results as Eqs (22) and (23) (Section I in S1 Text) and apply them to derive the spectrum and dimension for the i.i.d. Gaussian connectivity under spatial sampling. In particular, the relative dimension increases with spatial sampling (i.e., decrease f) which is consistent with the shape of the spectrum where its support narrows in (Fig 7). Lastly, the power-law feature is also preserved under spatial sampling. For any fixed 0 < f ≤ 1, we show that as g → 1 and x → ∞ (see examples with g close to 1 in Fig L in S1 Text)

3.8 Fitting the theoretical spectrum to data

Our theory for the bulk covariance spectrum can be fitted to neural activity whenever the covariance eigenvalues can be calculated. The best-fitting theoretical spectrum can be found by minimizing the L2 or L∞ error between the empirical and theoretical cumulative distributions (Methods) with respect to parameters such as g. We note that the availability of closed-form or analytic solutions of the theoretical distributions makes this optimization highly efficient.

In many settings, the value of the baseline neuronal variance σ2 in Eq (2) is not known. But this can be easily addressed by scaling both the observed and the predicted eigenvalues to have a mean equal to 1. After fitting the connectivity parameter g for the normalized eigenvalues, σ2 can then be estimated using the original means of data and theory. For our theoretical spectra, the mean μ of covariance eigenvalues is available in closed-form (Eqs (9) and (18)), and the scaled pdf is easily found as p(x) = μp(μx).

Furthermore, the recorded neural activity is sometimes normalized for each neuron (e.g., by converting activity to z-scores). In this case, we need to analyze the eigenvalues of a correlation matrix whose entries are normalized as . Interestingly, we found that the correlation eigenvalue distribution for our random connectivity models in the large N limit is the same as the rescaled p(x) above. This is because the diagonal entries of C become uniform (thus converge to μ) for large N (Section J in S1 Text).

The fitting of the spectrum is also robust to outliers in the covariance eigenvalues (Section 3.4). In Section K in S1 Text we demonstrate an example where a rank-2 component is added to the covariance C. Since in practice the rank of the perturbation is unknown a priori, we use all eigenvalues in the fitting, and the fitted g is highly accurate despite the presence of outliers. We can also use the fitted g to help identify the outliers by separating them out based on the upper edge of p(x) support [32, 51].

We conclude with a preliminary application to whole-brain calcium imaging data in larval zebrafish. In [52], activities of the majority of the neurons in the larval zebrafish brain were imaged simultaneously during presentations of various visual stimuli and grouped into functional clusters based on their response similarity. These clusters reveal potential neural circuits and, in some cases, reveal a good match with known brain nuclei. Here we select a few clusters that contain a large number of neurons and are anatomically localized (Fig 8B). For each cluster, we calculate its sample correlation matrix during the spontaneous condition (no stimulus was presented) and then fitted the eigenvalues to the time-sampled spectrum with i.i.d. Gaussian connectivity (Section 3.7). Here the calcium activity is used but the correlation matrix is expected to be approximately the same if firing rates or long time window spike counts were used (Section K in S1 Text). Despite the simplicity of the model with only one parameter to tune, the results show good agreement with data and is significantly better than fitting using the Marchenko–Pastur law (Fig 8C), which models spatially independent noise (Section 5.4). Therefore, our theory provides a quantitative mechanistic explanation of how a long tail of covariance eigenvalues or equivalently low dimensional activity occurs in recurrent neural circuits.

Fig 8

Fitting the theoretical spectrum to data.

A. The anatomical map of neurons (dots) in the example functional clusters (different colors) across a larval zebrafish brain (scale bar is 50 μm, see text and [52]). B. Comparing the fitting error of the time-sampled random connectivity theory (Section 3.7) and the Marchenko–Pastur law. The errors are measured by Eq (38). The red dashed line is the diagonal. labeled on each plot is the relative dimensionality (Eq 4). The calcium activity is recorded at a frame rate of 2 Hz and a total of 600 frames of spontaneous activity [52] are used in here. See more details in Methods. Fitting results for all other clusters are in Fig O in S1 Text.

4 Discussion

In this work, we studied the eigenvalue density distribution of the covariance matrix of a randomly connected neuronal network, whose activity is approximated as noise driven linear fluctuations around a steady state. We derived an explicit expression for eigenvalue distribution in the large-network limit analytically in terms of the statistics of the connectivity such as coupling strength and second-order motifs. Our results also include closed-form expressions for the dimension measure generalizing known results [6]. Some of these dimensionality results are also derived in a recent parallel work [37], whereas the shape of the covariance spectrum was not studied in [37]. Knowing the exact shape and support of the bulk spectrum can facilitate separating outlying eigenvalues corresponding to low dimensional structure (coming from other unmodeled effects such as external input) from variability due to noise [51] (see an example in Fig N in S1 Text). The shape of the bulk spectrum reflects structured amplification of the neuronal noise by the random recurrent interactions. Since the bulk spectrum is not altered by low rank perturbations to the connectivity or to the activity (Section C in S1 Text), this allows for distinguishing different sources of variability in neural data. As the connection strength increases towards the critical edge of stability, the spectrum exhibits a power-law tail of large eigenvalues, with exponent −5/3 in pdf (or −3/2 in eigenvalue vs. rank plot). This power-law shape near the critical g provides concise theoretical characterizations of the spectra under various connectivities. Intriguingly, the same power law persists even when the shape of the spectrum is modified by connectivity motifs or due to finite temporal and spatial sample size. In contrast, when we move away from the asymmetric, random connectivity model, the exponent of the power law (if any) becomes different: −7/4 for symmetric random connectivity (Eq (16)), −2 for a normal connectivity J with matching eigenvalue distribution as i.i.d. Gaussian J (Section E.3 in S1 Text), and −d/4 − 1 for a d-dimensional ring network (see below). Based on these results, we conjecture that a power-law tail, whenever present for any covariance spectrum, reflects the qualitative nature of the connectivity and is a robust feature that will survive both temporal and spatial sampling with the same exponent (precise statement in Section I.3 in S1 Text).

Unlike the large eigenvalues [53], the interpretation of the bulk spectrum of PCA of neural activity data has received little attention. A notable exception is a recent work [54] which studied the power law of covariance spectrum of data near criticality based on the renormalization group method. By fitting experimental data to the theoretical spectrum, information from eigenvalues of all sizes can be used to estimate the effective connection strength (Fig 8). Our theory thus provides an important benchmark to compare with experimental data and advocates the bulk covariance spectrum as a powerful global description of collective dynamics in neuronal networks.

4.1 Nonlinear dynamics

One limitation of the work is the assumed dynamic regime where fluctuations of the neuronal activity are described by the linear response theory [22, 23] around a fixed point. While extending the theory to highly nonlinear activity such as chaotic dynamics [2] is left for future research, we provide some numerical examples of networks with nonlinear neurons to illustrate the applicability of our results.

We consider a model with nonlinear rate neurons driven by external noise by introducing an activation function ϕ(x) = tanh(x) in Eq (1) [2] [55]. The nonlinearity transforms the currents h(t) to firing rates r(t) = ϕ(h(t)), and the recurrent interaction term is now (Eq (33) in Methods). When the connection strength g = Nvar(J) and the noise magnitude σ are small, the neural activity will mainly fluctuate near 0 where ϕ(x) can be approximated by linearizing around 0 (Fig 9A, bottom). Indeed, the spectrum based on the linear theory (green dashed line in Fig 9A, top) approximates the numerical spectrum. For larger g and σ the deviation from the linear theory becomes significant as the firing rates are now strongly shaped by the nonlinearity (Fig 9B, 9C and 9D, bottom). Interestingly, these spectra can be well fitted by the linear-theory spectrum if g is replaced by a smaller effective . This reduction in the effective g can be qualitatively understood as the average slope of ϕ(x) over the distribution of h(t) (Fig 9 bottom). Note that, unlike the noiseless model in [2], the cases studied here are not chaotic even with g > 1. As shown in [55], the presence of external noise shifts the transition to chaotic dynamics from g = 1 to larger values. A theoretical characterization of effective g and other strongly nonlinear dynamics such as chaotic dynamics is left for future research. Future work could also consider cases with transient activity which are common in nonlinear systems and more general network architectures, such as multiple populations of EI networks [21] and incorporating distant dependent connectivity patterns based on known cortical microcircuit architectures [56-58].

Fig 9

Covariance spectrum in nonlinear dynamics.

A. Top: A histogram of covariance eigenvalues calculated from firing rate activities r(t) of simulating a N = 400 network model according to Eq 33. The eigenvalues are normalized to have a mean equal to 1 for easy comparison of the shape (or equivalently the eigenvalues of the correlation matrix Section 3.8). Here g = 0.4 and σ = 0.5 (see Methods for additional numerical details). The green dashed line is the time-sampled theoretical spectrum (Section 3.7) using actual g and α, only shown in A for clarity of the plots. The length of the simulated data corresponds to α = N/T = 0.1. The orange curve is also the time-sampled theoretical spectrum except for using an effective that is fitted numerically to best match the simulated eigenvalues. Bottom: The blue curve is the histogram of h(t) (aggregated across i and t) and the orange dashed curve is the histogram of r(t). The overlaying red curve shows the nonlinearity ϕ(x) as a reference. The 〈⋅〉 in the title of each plot represents averaging ϕ′(h(t)) over all i and t. B-D. Same as A except for σ = 1 and g = 0.6, 0.8, 1.2, respectively. Only the fitted theoretical spectrum (orange curve) is shown for clarity.

4.2 Ordered vs. disorder connectivity

We have studied the covariance spectrum under random connectivity, which is used as a model for complex recurrent networks. Here we ask whether features of these spectra are distinct results of the connectivity being random. To address this, we briefly explored the covariance spectra from several widely used examples of ordered connectivity for comparison.

First, consider a ring network [56] with translation invariant long-range connections, where the connection strength between neurons depends smoothly on their distance (Fig 10A inset, see Methods). In the large-network limit, the covariance spectrum becomes a delta distribution at 1 with a few discretely located large eigenvalues (Fig 10A). Next, we consider the ring network with short-range, in particular, Nearest-Neighbor (NN) connections. The covariance spectrum is now continuous (no outliers) and supported on an interval, but the pdf diverges at both edges as (Fig 10B).

Fig 10

Covariance spectra under some deterministic connectivity models.

A. Histogram of the covariance eigenvalues of a ring network with a long-range connection profile (inset, N = 100). Most eigenvalues are close to 1 and the rest of eigenvalues converge to discrete locations predicted by top Fourier coefficients (crosses) of the connection profile (Eq (36)). B. Same as A. but for a ring network with Nearest-Neighbor connections: J = 0.4, J = 0.2. The solid line is theoretical spectrum in large N limit which has two diverging singularities at both support edges. The effect of such singularities is also evident in the finite-size network at N = 400 (a single realization). C-F. Higher dimensional Nearest-Neighbor ring network (ad = 0.6, see Methods). As the dimension increases, the singularities in the pdf become milder and less evident, and the overall shape becomes qualitatively similar to the random connectivity case (Fig 1).

To seek further examples of ordered connectivity leading to a qualitatively similar covariance spectrum as the random connectivity, we consider the d-dimensional generalization of the NN ring (Methods). As dimension d increases, the smoothness of the pdf within and at the edges of the support increases, and the covariance spectrum becomes qualitatively similar to the random case [59] (Fig 10). Interestingly, as the connection strength approaches its critical value for stability, the covariance spectra also exhibit a power-law tail (Section L.3.1 in S1 Text; is also possible under other cases using [60]). To match the exponent of the random network d would be 8/3 ≈ 2.67. These comparisons suggest that the covariance spectrum’s overall smooth density and long tail shape is a shared property in highly connected networks with high rank connectivity matrices, including random networks and high dimensional short-range spatially invariant networks.

5 Methods

5.1 Models of random connectivity

Here is a summary of results on various random connectivity models.

i.i.d. Gaussian random connectivity : closed-form pdf and endpoints (Eq (5)), including the frequency-dependent covariance spectrum (Section 3.6), and a power-law tail approximation (Eq (7)).

Gaussian random connectivity with reciprocal motifs/asymmetry κ = ρ(J, J) (Section 3.3): analytic solution and endpoints (quartic root, Section E in S1 Text) and a power-law tail approximation (Eq 17). For special case of symmetric and ant-symmetric connectivity, closed-form pdf Eqs (14) and (15), including a frequency dependent covariance spectrum (Section G in S1 Text).

Erdős-Rényi and certain EI network Section 3.5: same bulk spectrum as the i.i.d. Gaussian case.

For all cases, the mean μ and the dimension D are derived in closed-form (Eqs 9, (10) and (18).

For simplicity, we do not require J to be zero (i.e., no self-coupling), but allow it, for example in the i.i.d. Gaussian model, to be distributed in the same way as other entries J. In the large-network limit, since individual connections are weak (e.g., ), allowing this self-coupling or setting J = 0 does not affect the covariance spectrum (Section A.3 in S1 Text).

5.2 Applications to alternative neuronal models and signal covariance

Although the relation between C and J (Eq (2)) is derived here in a linear rate neuron network Eq (1), it also arises in other models of networked systems.

Linearly interacting Poisson neurons

This is also called a multivariate Hawkes model [27]. This is a simple model for spiking neuron networks, but is versatile enough to capture for example the temporal spiking correlations seen in other more sophisticated nonlinear spiking neuron networks [26, 28]. A time-dependent Poisson firing rate is calculated as a filtered input spike train s(t) (sum of delta functions), and spikes are then drawn as a Poisson process given y(t),

Here we consider a homogeneous network where the baseline firing rate y0 and response filter A(t) is the same for all neurons.

The exact long time window spike count covariance matrix of this network can be shown to be [27] which is valid if the time varying y(t) does not often become negative (for example when any negative connections W are small compare to y0). Here , Y0 and Y are vectors of baseline and perturbed (with recurrent connections) firing rates of the neurons respectively. If we assume that the effective connection strength aW is weak so that we can approximate Y with Y0, then (26) becomes the same as Eq (2) with J = aW (note that for Poisson process ).

Another condition that ensures a uniform Y and does not restrict connections to be weak is a row balance condition of W sometimes assumed in EI networks [45],

This is not unreasonable to assume, for example, considering the homeostatic mechanisms of neurons.

Integrate-and-Fire neurons

As shown in [23, 26] the linear response theory [22] can be used to approximately describe the covariance structure of a network of generalized integrate-and-fire (IF) neurons

Here V is the membrane potential and a spike is generated when V reaches a threshold. y(t) = ∑ δ(t − t) is the spike train, and in Eq (29) is the “unperturbed” spike train in absence of recurrent connections W. Different choices of ψ(V) realize types of IF neurons, such as ψ(V) = Δ exp((V − V)/Δ) for the exponential IF neurons. During the asynchronous firing of neurons (no strong synchronized firing across the network), Eq (28) can be well approximated by

Here denotes convolution. W = {W} is the matrix of recurrent connection weights. A(t) is the linear response kernel for neuron i (e.g., an exponential decay) F(t) is the temporal kernel of the synapse. For simplicity, we assume that both A and F are uniform across the network.

It it shown [23, 26] that the long time window spike count covariance matrix C (in fact also the frequency covariance Eq (19)) is well approximated by

Here the scalar summarizes the cellular and synaptic dynamics. can be thought of as the baseline neuronal variance in the absence of recurrent connections (A = 0 in Eq (29). This expression of the covariance matrix again matches with Eq 2.

Fixed points over whitened input

The covariance we considered so far describes the structure of fluctuations of spontaneous dynamics without or under fixed external input, often referred as the noise covariance [61]. We can also consider the spectrum for a signal covariance. This perspective is needed to use our results to interpret experimental data where the neural activity is largely driven by stimuli for example [10].

Consider a network of linear firing rate neurons,

Here u is the external input to neuron i. Assume the network settles to a steady state, so the fixed point firing rates are given by

Now consider the network activity across an ensemble of input patterns, which has whitened statistics [62],

It is easy to see that the covariance of firing rates is given by σ2(I − J)−1(I − J)−, which is the same as Eq (2).

We note that Eq (31) or equivalently appears in broader contexts beyond neuroscience and is studied in the field of linear structural equation modeling (SEM) [29].

Nonlinear rate neurons

The following is a simple and classic nonlinear rate neuron model similar to that used in [2],

Here ϕ(x) = tanh(x) converts currents h(t) into firing rates r(t) = ϕ(h(t)). The white noise ξ and i.i.d. Gaussian random J are the same as in Eqs (1) and (3). The presence of noise allows for nontrivial activity for g < 1. Due to symmetry of ϕ(x) and J distribution, the average activity over time 〈h(t)〉 = 0 for large N. Note that the largest slope of ϕ(x) is 1 at x = 0, our main model Eq (1) can be viewed as a linear approximation to Eq (33) and the definition of connection strength g is consistent.

In Fig 9, N = 400, τ = 1 and Eq (33) is simulated using the Euler-Maruyama method with a time step of 0.01 for a duration of T0 = 40000. To get the long time window covariance, the firing rates r(t) are binned by a time window of 10, which is sufficient based on examining the decay of its autocorrelation function. After binning, the simulated data correspond to a time-sample parameter α = N/T = 0.1 (3.7) and this α is used in calculating the theoretical spectra. The empirical covariance matrix is then calculated from the binned firing rates according to Eq (37). As shown in Fig 9, the theoretical spectrum based on the linear dynamics can still be used to fit the simulated spectrum in this nonlinear network when g is replaced by a smaller effective .

5.3 Power-law approximation of the eigenvalue distribution

The power-law property of p(x) for i.i.d. J under critical g is probably known in random matrix theory (private communication), by results from the equivalent distribution studied in [30], although we do not know of a specific reference. We include a derivation based on the explicit expression of p(x) (Eq (5)) that is outlined below.

First note the limits of the support edges. As g → 1−, . For the lower edge, can be found by the Taylor expansion of (1 − g2) or note that (1 − g2)3x+x− = 1. Consider a x that is far away from the support edges as g → 1, given the above, this means,

Note that since x+/x− ∼ (1 − g2)−3, there is plenty range of x to satisfy the above for strong connections when g is close to 1. Under these limits, Eq (5) greatly simplifies as various terms vanish leading to

This explains the validity of the power-law approximation away from support edges. If we are only interested in the leading-order power-law tail in the critical distribution (i.e., g → 1− and then x → ∞), then there is a simpler alternative derivation that can we also apply to other connectivity models (see Section A.2 in S1 Text).

5.4 Comparison with the Marchenko–Pastur distribution

The Marchenko–Pastur distribution is widely used for modeling covariance eigenvalues arising from noise [32, 51, 63]. It is also the limit of the time-sampled spectrum p(x) (Fig 7 and Section 3.7) at weak connections g = 0. The Marchenko–Pastur law has one shape parameter α. We focus on the case when the covariance is positive definite which restricts 0 < α < 1 (otherwise there is a delta distribution at 0) and the pdf is

The first two moments are 1 and 1 + α, from which we know the dimension is 1/(1 + α) has a lower limit 1/2. The upper edge α+ is bounded by 4.

5.5 Deterministic connectivity

5.5.1 Ring network with short- and long-range connections

In a ring network, neurons are equally spaced on a circle (can be physical or functional space) and neuron i is associated with a location x = i/N, i = 0, …, N − 1. The connection between two neurons j and i only depends on the location difference x − x is therefore translation invariant.

For long-range connections, the connectivity has a shape determined by a fixed smooth periodic function f(x) on [0, 1),

In the large-network limit, the covariance eigenvalues have an approximate delta distribution at 1 except for a finite number of discretely located larger eigenvalues (Fig 10A). A precise statement of this result is described in Section L.1 in S1 Text. The outlying eigenvalues correspond to the leading Fourier coefficients of f(x).

For the Nearest-Neighbor (NN) connectivity, only J and J are non-zero and remain fixed as N → ∞.

5.5.2 Multi-dimensional ring network

For a d-dimensional ring, the neurons are equally spaced on a d-dimensional lattice which is periodic along each coordinate. We focus on the NN connectivity where each neuron is connected to 2d neighboring neurons with strength and along direction k. We show that the probability density function at both support edges scales as (for comparison, the random network edges scale as ). This means for dimension d ≥ 2, there is no singularity at the support edges (Fig 10).

To characterize the shape of the covariance spectrum (Fig 10), we further simplify by setting (see also Section L.3 in S1 Text for motivations based on 1D ring) and analytically derived p(x) (Section L.3.2 in S1 Text). For small dimensions d ≤ 3, there are distinct “inflection points” within the support. As d increases, this non-smooth feature becomes less evident and becomes hard to identify in empirical eigenvalue distributions from a finite-size network (not shown).

5.6 Fitting the theoretical spectrum to data

For neural activity data, C can be calculated from a large number of time samples of binned spike count s(t) (assuming bin size is ΔT large),

For calcium imaging data, the fluorescence is approximately integrating the spikes over the indicator time constant. So we can still apply Eq (37) by plugging in the fluorescence signal in place of s(t) to calculate the covariance C (omit the constant factor ΔT which does not affect fitting to the theory, Section 3.8).

We fit the theoretical spectrum to empirical eigenvalues by finding the connectivity parameter g that minimizes the error between the cumulative distribution functions (cdf) . This avoids issues such as binning when estimating the probability density function from empirical eigenvalues. We numerically integrate the theoretical pdf (Eq 5) to get its cdf. As seen below, the theoretical cdf only needs to be calculated at the empirical eigenvalues.

Motivated by methods of hypothesis testing on distributions, we measure the L2 norm cdf error using the Cramer-von Mises statistic

Here n is the number of samples and x are the i-th empirical eigenvalues. Alternatively, the error can also be measured under L∞ norm based on the Kolmogorov-Smirnov statistic where x are samples. Our code implements both measures.

In Fig 8B–8E, we fit the time-sampled theoretical spectrum with i.i.d. Gaussian connectivity (Section 3.7) to calcium imaging data in larval zebrafish [52]. The theoretical spectrum (once normalized by the mean, see Section 3.8) depends on two parameters g and α, but the latter is fixed to be N/M based on the data. Here N is the number of neurons in a cluster, and M is the number of time frames used in calculating the sample correlation matrix (Eq (37). The calcium fluorescence ΔF/F traces of each neuron are normalized to z-scores [52], which is consistent with calculating the eigenvalues of the correlation matrix (Section 3.8). g is then optimized to minimize the Cramer-von Mises error (Eq 37) between the data. The largest eigenvalue for each cluster is often much larger than the rest and is thus removed before the fitting. For comparison, we fit the same data to the Marchenko–Pastur law (Section 5.4) whose shape depends on the parameter α. Here we allow α to vary so that both models (random connectivity and MP law) have one parameter to be optimized to fit data.

Supplementary material.

Analytical derivations and additional figures.

(PDF)

Click here for additional data file.

16 Feb 2022

Dear Dr. Hu,

Thank you very much for submitting your manuscript "The spectrum of covariance matrices of randomly connected recurrent neuronal networks" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Tianming Yang

Associate Editor

PLOS Computational Biology

Thomas Serre

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors calculate the eigenvalue distribution of the covariance matrix of recurrent rate networks linearized around a fixed point driven by additive Gaussian white noise both for i.i.d. Gaussian connectivity and for various non-random connectivities, e.g. symmetric, antisymmetric, low-rank perturbations to i.i.d. Gaussian, ring network and describe the effects of partial temporal or spatial sampling.

The authors fit the eigenvalue distribution to whole-brain Calcium imaging of larval zebrafish and find a better match of the i.i.d. Gaussian connectivity compared to a Marchenko-Pastur distribution.

The contribution of the authors include:

· The analytic PDF of covariance matrix eigenvalues for linear rate networks with Gaussian connectivity.

· A detailed picture of how these distributions are shaped by structural perturbations such as low-rank updates, partial sampling, and partial (anti)symmetry

· A comparison of the eigenvalue distribution obtained from experimental recordings in larval zebrafish to both Marchenko-Pastur distribution and the novel analytically obtained eigenvalue spectrum

· A statement of common characteristics of the covariance matrices and covariance matrix eigenvalue densities, between the analytically solved random linear case and deterministic connectivity cases.

We wholeheartedly recommend this to be published in PLOS CB. The mathematical results are impressive, of high relevance to the audience of PLOS CB, manifold network structures are studied and the results seem mathematically rigorous (especially the Supplement to the extent that we studied it.) However, we feel the utility for the readership of PLOS CB would increase by considering the following major and minor comments:

Major Comments:

· While the mathematical results are very impressive, a detailed discussion of the limitations is largely missing:

- Under what conditions yields the PCA the dimensionality of a data set? (When can we assume the data follows a multivariate Gaussian distribution?). Is that the case for neural data? Why are your findings relevant for neuroscience?

- What are the limitations of a linear response theory? When would it break down?

- What is the relationship between ‘linear rate units’ and ‘neurons’. (Often ‘rate units’ are considered as populations of neurons. But how does this fit together with the reference to single neuron motifs (e.g. Song et al. 2005?))

- The relevance of the investigation of g approaching 1 from below is not clear. Isn’t at g=1 the network dynamics turning unstable and the linear response picture breaking down? Already close to g=1, the fluctuations are getting large, so one might expect that the linear response description is getting increasingly inaccurate because \\phi(x)≠x. So while the power-law scaling is an interesting observation, I am not sure it has any relevance for neuroscience (but of course, I’d love to be proven wrong on this).

* While ‘the equations speak for themselves’, an intuitive explanation of the observed results and the underlying mechanism is largely missing:

- Why does g->1 lead to a few dominant eigenvalues in the covariance matrix?

- Why do antisymmetry networks have a qualitatively different shape of the distribution of eigenvalues?

- Why do (anti)symmetry networks results in lower (higher) relative dimension D/N? Why does dependence of D(\\kappa) change when fixing g_r?

-

The Discussion section seems to add more new results (e.g. reference to ring network and other results in the Supplementary Material) instead of discussing the interpretation of the results and implications for theory and experiments.

Minor Comments:

* Consider changing the title to “The spectrum of covariance matrices of neuronal networks with linear dynamics” or “The covariance spectrum of recurrent neural networks with linear dynamics” or another title that contains the word 'linear' or 'linearly' would make it more clear that here a linear response framework is being used. This also leaves space for future studies that explore nonlinear regimes.

* Figure labels are so small that they are too small. Please increase the label and legend font size to make it readable.

* Make clear that PCA is based on pairwise statistics, so it doesn’t require simultaneous recording of all N neurons, it would be enough (assuming stationarity) for all pairs to be recorded one after the other. Therefore, I would like to politely disagree with the notion of 'local' and 'global' features, they might mislead some readers as only two-point interactions are considered here.

* The comparison of the experimentally obtained covariance spectra of zebrafish data to analytically obtained ones includes assumptions about the temporal correlations of the data (Marchenko-Pastur assumes that samples are independent, i.e. sampled at time frames that are sufficiently far apart such that temporal correlations can be ignored, to what extend is that assumption justified in your data? (It would be helpful to report the frame rate in the caption of the figure 8 and the decay time constant of GCamp6f and the number of time frames used to calculate the empirical covariance matrix. I think your reference 10 has a frame rate of 2.11 per s and GCamp6f has a decay constant of 1796±73ms according to (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3777791/#SD1) but please check.

* The experimentally recorded zebrafish was exposed to visual stimuli and the data was z-scored according to your methods, please comment on how both spatiotemporally structured external input and z-scoring affects your results.

* Line 126 should read semi-positive instead of positive.

* Line 139: Cite earlier papers on participation ratio.

* Line 160: “Matches pretty well” use precise language.

* Line 167/Figure 1A,B: showing more than one network realization would be more convincing.

* Figure 1: labels and legend too small

* Figure 1A,B: consider doing a stair step graph for the empirical histogram across multiple (e.g. 10) network realizations

* Line 187: c_0 seems not to be introduced in the main text.

* When you refer to Supplement, it would help if you refer to the section of the supplement because it is 60 pages long so things would be easier to find.

* Figure 5: inset too small, axes not legible

* Line 296: only K<>1 for the Gaussian assumption.

* It would be nice to plot both D(f) and D(\\alpha) (equation 23) next to figure 7AB so the effect of sparse time sampling becomes even visually more obvious.

* A closer look at the methods section (5.6) indicates that the authors incorporated temporal sampling corrections in their calculation of the theoretical spectral CDFs used to fit the neural data. However, it is not clear in the main text (3.8) that they did this.

* In section 3.8, where the authors fit distributions to data, a discussion of the data is missing. Would you expect the same PCA distribution for spikes instead of calcium? What is the effect of the Calcium response? It would be good to comment on the fact that the eigenvalues of a covariance matrix change after a change of variable.

* In figure 8, important details are missing: How many data points are used to estimate the covariance matrix? Why don’t you fit the model including symmetry/antisymmetry? What D/N do you numerically estimate for the different clusters? Do you have full access to all neurons Calcium activity? If not, do you apply your spatial sparsity model?

* Instead of motifs-> second-order motifs. (there are also higher order motifs)

* Please go carefully over all references to correct typos and misspelled names.

Stylistic:

* “iid” -> “i.i.d.”

* Highly complex->complex

* “Frequency spectrum in Fourier transform”-> “Frequency spectrum obtained via Fourier transform”

* “co-fluctuations” doesn’t to my knowledge exist, maybe use covariation instead?

* Instead of “the network’s Principal Component Analysis, we would propose the network dynamics Principal Component Analysis.

* Line 237: For denoting the imaginary i, I would not use a double-struck i.

Typos:

* Whitespace line 71 after dynamics.

* Excitatory-Inhibitory-> excitatory-inhibitory

* Line 559: “random matrix” -> “random matrices”

* In abstract: “The theoretical results are compared with those from finite-size networks and and the effects” should read “The theoretical results are compared with those from finite-size networks and the effects”.

* “the rank plot with exponent”

* “supported by a NIH grant” -> “supported by an NIH”

* “which correspond to overabundance of certain subgraphs” -> “which correspond to an overabundance of certain subgraphs”

* “as g approach the critical”-> “as g approaches the critical”

* “In large-network limit”-> “In the large-network limit”

* “the error can also be measure under” -> “the error can also be measured under”

* “Storing Infinite Numbers of Patterns in a Spin-Galss Model of Neural Networks”-> “Storing Infinite Numbers of Patterns in a Spin-Glass Model of Neural Networks”

* “Multiplciation of certain noncommuting random variables”. -> “Multiplication of certain noncommuting random variables”.

Rainer Engelken

Reviewer #2: Summary:

--------

In this manuscript, Hu and Sompolinsky derive analytical expressions

for the eigenvalue spectrum of covariance matrices and the related

dimensionality measure (participation ratio) for random networks (iid

weights / motifs / excitatory-inhibitory, and also some structured

networks for comparison). Their theory is based on linear(ized)

dynamics, where there is a simple relation between covariances and

connections. Yet, for the general case of non-normal connectivity

matrices, the derivation of the covariance spectrum from connectivity

parameters is involved. Their theoretical derivations, that for better

readability are mostly presented in the extensive supplement, are

based on the seminal work by Sommers, Crisanti, Sompolinsky and Stein

in 1988.

Given the large interest in dimensionality of neural activity in the

computational neuroscience community over the recent years, this work

is highly relevant and timely and an important contribution to the

understanding of the relation between structure and dynamics in neural

networks. The presentation is mostly clear, the theoretical

derivations are solid and the authors also show an application to

experimental data (whole-brain calcium imaging data from larval

zebrafish). We therefore believe this paper should eventually be

published in PLOS CB. There are, however, a few major and minor issues

(see below) that need to be resolved in a revised version before the

manuscript is ready for publication.

Major points:

-------------

-------------

Effect of other motifs:

-----------------------

The authors study in depth the effect of reciprocal connections /

symmetry of connections on the covariance spectrum. They also discuss

the other motifs (divergent, convergent, chain), but state that they

do not "affect the bulk spectrum of C" (main text line 220). In the

supplement below Eq. S26 the authors, however, state that also

divergent, convergent and chain motifs affect the bulk spectrum

indirectly by channeling their effect through changing the normalized

reciprocal motifs. Can the authors clarify this and elaborate in more

detail in the main text, which motifs affect the bulk spectrum in

which way? Also the critical coupling strength / stability should

depend on all motifs (see Eq. S26), so one would expect an influence

of all motifs on the bulk covariance spectrum.

Clarification of limitations:

-----------------------------

In Section 3.5 the authors should stress more that their theoretical

predictions are limited to networks where all connections have identical

variance. Their choice of E-I network is a special case that is

explicitly constructed such that this requirement holds. In typical

E-I networks the variance is not the same due to different population

sizes, connection probabilities, and synaptic strengths (that are not

precisely tuned to achieve the same variance of E and I connections).

Relation to Sommers et al. 1988 (S1988):

----------------------------------------

In supplement S1 the authors show how to derive the probability

density p_c for the covariance eigenvalues from the potential Phi that

has been computed in S1988. What did not become clear to us is how the

derivation deviates from the calculations in S1988. Because S1988

calculated the probability density p_J of the connectivity spectrum

rather than p_c. Both calculations seem to rely on the precision

matrix, but the detailed differences in calculations did not become

clear from reading the supplement. Can the authors please provide more

details in which steps of the calculations results from S1988 are

exactly used and in which steps results needed to be adapted to treat

the covariance spectrum rather than the connectivity spectrum?

Relation to Stringer, Pachitariu et al. Nature 2019 (SP2019)

------------------------------------------------------------

In line 188, the authors present their results as an alternative

mechanism for the experimentally observed power laws of covariance

spectra in SP2019. The authors, however, study the covariance

structure of spontaneous network dynamics while SP2019 studies

stimulus-evoked activity (responses in V1 due to visual stimuli).

Would one thus not expect differences in covariance spectra? Can the

authors briefly elaborate on this?

Relation to Dahmen, Recanatesi et al. biorxiv 2020 (DR2020):

------------------------------------------------------------

In line 208 the authors state that the dimensionality result Eq. 10

can also be derived based on results of their reference [11]. This has

already been done in DR2020

(https://www.biorxiv.org/content/10.1101/2020.11.02.365072v1.full.pdf).

It would be helpful to explain the relation of the current work to

these parallel developments. This seems particularly important given

the high similarities between figures of the current manuscript and

figures in DR2020, in particular

-Fig 1D <-> DR2020:Fig 2B

-Fig 4B <-> DR2020:Fig 4B

The approximation in l.257 for D/N as a function of g_r corresponds to

the approximate result in DR2020.

Also in the supplement "First two moments as a corollary of results in

[3]" it would be helpful of the authors cited DR2020 to clarify

the close relation of these parallel works.

Motivation of the manuscript:

------------------------------------------------

To motivate their work, the authors dinstinguish in the first

paragraph of the introduction between local and global features of

dynamics. They cite a few works falling in the local feature category

and argue that their work targets the global scale.

We, however, disagree with the authors' categorization.

While their refs [55,48] really discuss local features,

refs [27,11,18] study correlations (mean and variance across populations).

The latter emerge from collective network effects and should thereby be

attributed to the global feature category. In fact, since the dimensionality

and covariance spectrum are based on the very same covariances, the

authors' work falls in the same category as [27,11,18].

Minor points:

-------------

-------------

Power law exponent:

-------------------

Why is the power law shape of the covariance spectrum not continuous

in the parameter kappa? I.e. why is there a single power law exponent

for -1kappa=-1? Are these results an artefact of the calculations that are

different for kappa \\in {-1,1} and -1numerical comparisons for the power law exponents for various kappa?

The authors should extend their discussion on these points.

Wrong cross-references in the supplement:

-----------------------------------------

Please check again the cross-references in the supplement. Some of

them textually refer to Figures in the main text, but the hyperlink

takes us to figures in the supplement, e.g. in l.474: Ref to Fig 6

(in the main text) wrongly links to Fig. S6 in the supplement.

Figures:

--------

Fig 1: Colors for theory and simulation not consistent across panels

Fig 4A, Fig 5C&D, Fig 6A&B: legend labels for lines not consistent:

decide for "theory" or "Gaussian theory"

Fig. 6: Better explain the correction with the modified connection weight.

Fig. S6: y-axis labels are cut off

Strong non-normal effects:

--------------------------

Sections S5.2 and S5.3: I am not sure which point the authors want to

make here. The bimodality of the spectrum is an interesting

observation, but the comparison to the covariance spectrum of the

matching normal connectivity shows that the bimodality is not an

effect of the non-normality. In contrast the bimodality is more severe

for the matched normal. So the non-normality seems to suppress the

bimodality. If this is the main message then the authors should write

this more clearly in the supplement and also in Sec 3.3.2 of the main

text. How do the covariance spectra and matched normal spectra compare

for kappa_re>0?

Time and space sampling:

------------------------

We found supplemental sections S8 and S9 hard to follow, especially

for readers that are not familiar with free probability theory. We are

aware that the authors cannot go into all details on this issue in

this manuscript, but it would be helpful for the average reader of

PLOS CB if the authors could provide even more intuition on the

various steps of the calculations and possibly single out a good

references for an overview article that introduces the concept

of free prbability theory.

References:

-----------

l.272: We think it would be suitable to also cite Schuessler et al.

NeurIPS 2020, which showed that training random networks creates

low-rank connectivity components.

More detailed minor points (Supplement):

----------------------------------------

----------------------------------------

Eq.S3: indicate that -z is the argument corresponding to epsilon in

Eq.S2, i.e. d Phi/d epsilon(-z,\\eta)

In Eq. S2 epsilon shall be larger than zero in, but z can take on all

values in C excluding the real line. Could you please explain this a

bit more in detail?

l.165-166: Please define F_C also here in supplement.

Eq.S26: It it worth noting that \\hat{kappa}_*(J) is normalized by the

variance of J and \\hat{kappa}_*(\\tilde{J}) is normalized by the

variance of \\tilde{J}

Eq.S27: It is not clear where Sylvester's identity has been applied.

What is A and B in eq. S27?

l.340-341: z=eta? Please remind the reader here on the relation

between p_C and p_eta (shown in l.81).

Eq.S66: How to get from mu_1 to mu_1(C)?

l.426 and 436: Ref to Fig 3B correct?

Eq.S75: dx missing on the rhs?

l.513-514: what are overlines over omega and C referring to?

Typos (Main text):

------------------

l.22: and and -> and

l.125: the -> The

l.140: two moments -> two moments of

l.272: .[33,44]. -> [33,44].

l.528: broken sentence

l.536: weights) -> weights

l.559: random matrix -> random matrix theory

l.619: pdf -> cdf

l.623: measure -> measured

Typos (Supplement):

-------------------

l.121: E(\\lambda)^n -> E(\\lambda^n)

l.129: the its -> its

l.208: comma missing in between inline equations

l.354 and in following sections: p_P(x)->p_eta(x)

l.476: w0 ->w_0

l.523: C\\omega)->C(\\omega)

l.533-534: g missing in equation?

Reviewer #3: The manuscript analyzed theoretically the eigenvalue distribution of covariance matrix of the spontaneous neural activities in the randomly connected recurrent neuronal networks, and found the distribution has a finitely supported smooth bulk spectrum and exhibits an approximated power-law tail for coupling matrices near the critical edge. The simulation results consist with the theoretical prediction. I believe the mathematical derivations are rigorous and the results are solid.

However, I have some major concerns.

1. To analyze the eigenvalue distribution of covariance matrix, the authors made some important simplifications, for example, both the interactions between neurons and local dynamics of single neurons are linear. However, the local dynamics and interactions of our neural systems are all nonlinear, and lots of computing advantages are emerged based on these nonlinear characteristics. So could the authors provide some more simulation results to support that the main results can be extended to the nonlinear neural networks.

2. I am not clear what the eigenvalue distributions could give some new insights to the neural information processing mechanism? For example, the first few eigenvalues and eigenvectors of the covariance matrix can reflect the main computing characteristic or features of neural activities, so could the authors address the new computing advantages or functional meaning reflected by the eigenvalue distributions to the dynamics of neural activities?

3. The manuscript analyzed the covariance matrix of spontaneous activity, however, the neural systems usually show persistent activities when our brain performs some cognitive functions. Could the results extend to the study of persistent activities?

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Rainer Engelken

Reviewer #2: No

Reviewer #3: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

19 Apr 2022

Submitted filename: response to reviewers comments.pdf

Click here for additional data file.

1 Jun 2022

Dear Dr. Hu,

Thank you very much for submitting your manuscript "The spectrum of covariance matrices of randomly connected recurrent neuronal networks with linear dynamics" for consideration at PLOS Computational Biology. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. The reviewers appreciated the attention to an important topic. Based on the reviews, we are likely to accept this manuscript for publication, providing that you modify the manuscript according to the review recommendations.

Please prepare and submit your revised manuscript within 30 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to all review comments, and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Thank you again for your submission to our journal. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Tianming Yang

Associate Editor

PLOS Computational Biology

Thomas Serre

Deputy Editor

PLOS Computational Biology

***********************

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: General evaluation:

I thank the authors for addressing all major comments and questions. The manuscript is now much more clear and understandable. The additional figures greatly improve the presentation. I agree with the change of the title and consider the finding sufficient for publication in this journal.

There are just one minor issue remaining:

Minor comment:

Line 264: Instead of "Since J here is a normal matrix, this qualitative difference from the i.i.d. random connectivity can be understood considering the eigenvalues of J [49], which all lie on the imaginary axis and never approach 1 to cause large eigenvalues in the covariance when increasing g."

I would suggest writing instead for enhanced clarity: "Since J here is a normal matrix, this qualitative difference from the i.i.d. random connectivity can be understood considering the eigenvalues of J[49]. The eigenvalues of J all lie on the imaginary axis and therefore never approach 1. Thus, the eigenvalues of the corresponding covariance matrix don't develop a long tail."

Just for clarification: Our previous comment about line 296 was erroneous.

Reviewer #2: We thank the authors for addressing most of our comments in detail in

the revision. The presentation of results has much improved. There

are, however, a few points that are still unclear and need to be

clarified before the manuscript is ready for publication in PLOS CB.

Power law approximation:

------------------------

We apologize that parts of our previous question on the power-law

approximation seem to have been cut during the submission of the

comments. The missing part concerned the non-continuousness of the

power law exponent with kappa. We find unintuitive that there seems to

be an abrupt change in the distribution of bulk eigenvalues from e.g.

kappa=1- to kappa=1 (likewise for the other extreme case kappa=-1).

Because the approximated theory predicts different power law

exponents. Can the authors show simulations of p_c in logarithmic

yscale for readers to see what happens to the power law for fixed g_r

and varying kappa over the whole range including the edge cases

kappa=-1 and kappa=1 (Similar to Fig 4D but with log scale on y axis

and including cases closer and including kappa=+-1)? We understand

that the topology of connectivity eigenvalues changes from 2D to 1D at

the limiting kappas, but this happens in a smooth manner, so we would

have thought a smooth change in the covariance spectra.

Nonlinear network:

------------------

The analysis of the nonlinear network is a nice and convincing

addition to the manuscript. We, however, do find the authors'

description in lines 514ff potentially misleading. There the authors

state:

"Note that even for g > 1, the dynamics examined here is largely

driven by the noise rather than the intrinsic chaotic dynamics studied

in [50]."

It is not surprising that the linear approximation of the dynamics and

spectrum work well even for the case g>1 studied in Fig. 9. The

dynamics of all networks shown in Fig. 9 are not chaotic. In fact, the

dynamics are even linearly stable. As shown in Molgedey et al. PRL

1992 for discrete-time networks and in Schuecker et al. PRX 2018 for

continuous dynamics (same model as studied in this manuscript) the

transition to linear instability and the transition to chaos are both

shifted to larger values of g in the presence of external input (cf.

e.g. Fig. 3 in Schuecker et al.). Maybe the authors had this in mind,

but they should make the statement more precise to not suggest that

the approximation works well for networks with chaotic dynamics.

Minor points:

-------------

The xlabel in Fig 1D is missing.

Typos:

------

l479: a concise theoretical characterizations -> concise theoretical characterizations

l619: firing rates r_i(t) is binned -> firing rates r_i(t) are binned

Supp l.215f: This does not contradicting -> This does not contradict

Supp l.217: one introduce -> one introduces

Reviewer #3: I thank the authors for a thorough revision and comprehensively answering all of my questions. I have no further concerns.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Rainer Engelken

Reviewer #2: No

Reviewer #3: Yes: Yuanyuan Mi

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

References:

Review your reference list to ensure that it is complete and correct. If you have cited papers that have been retracted, please include the rationale for doing so in the manuscript text, or remove these references and replace them with relevant current references. Any changes to the reference list should be mentioned in the rebuttal letter that accompanies your revised manuscript.

If you need to cite a retracted article, indicate the article’s retracted status in the References list and also include a citation and full reference for the retraction notice.

14 Jun 2022

Submitted filename: response to reviewers comments_2.pdf

Click here for additional data file.

24 Jun 2022

Dear Dr. Hu,

We are pleased to inform you that your manuscript 'The spectrum of covariance matrices of randomly connected recurrent neuronal networks with linear dynamics' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Tianming Yang

Associate Editor

PLOS Computational Biology

Thomas Serre

Deputy Editor

PLOS Computational Biology

***********************************************************

15 Jul 2022

PCOMPBIOL-D-21-02214R2

The spectrum of covariance matrices of randomly connected recurrent neuronal networks with linear dynamics

Dear Dr Hu,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Zsofia Freund

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol