2.1 Neuronal networks with random recurrent connectivity

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

Here x_i(t) is the firing rate of neuron i. J_ij 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 ξ_i(t) with variance σ², that is, the expectation 〈ξ_i(t)ξ_j(t + τ)〉 = σ²δ_ijδ(τ).

We focus on the structure of long time scale covariation in the network, which are described by the long time window covariance . C_ij,ΔT is the covariance of the summed activity over a window of ΔT: C_ij,ΔT = 〈Δs_i(t)Δs_j(t)〉, . For biophysical neurons, C_ij,ΔT 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 (2)

Here I is the identity matrix, and A⁻¹, A^T are the matrix inverse and transpose (A^−T = (A⁻¹)^T). For simplicity we will set σ² = 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 (3) are drawn as independent and identically distributed (i.i.d.) Gaussian variables with zero mean and variance g²/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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. List of notations.

https://doi.org/10.1371/journal.pcbi.1010327.t001

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 λ_i 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_C(x) which is defined through the equality for all a, b. We also refer to p_C(x) as the spectrum (which should not be confused with the frequency spectrum obtained via Fourier transform). We will derive the limit of p_C(x) as N → ∞ and study how it depends on the connectivity parameters such as g = Nvar(J_ij).

The shape of p_C(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 (4) which can be calculated from the first two moments of p_C(x). We will also derive explicit expressions for D in random connectivity models.

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 λ_i 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), (5) where (6) and p_C(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).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

https://doi.org/10.1371/journal.pcbi.1010327.g001

The above result for the distribution p_C(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_C(x) does not have a simple relation to the spectrum of J because J is a non-normal matrix (i.e., J^TJ ≠ JJ^T). 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 − g²)⁻³ (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_C(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) (7)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

https://doi.org/10.1371/journal.pcbi.1010327.g002

As shown in Fig 2A, the power-law shape of p_C(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 − g² ≪ 1 and x ≫ 1. Define where x₊ ∝ (1 − g²)⁻³ is the upper edge of the support of p_C(x). Then, (8) 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_C(x) (Eq (5)) can be arbitrarily large as O((1 − g²)⁻²) (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 (9) which yields for the dimension (10)

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, (11) and in general, (12)

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_ij and J_ji are correlated, (13) 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_ij = J_ji, and κ = −1 corresponding to anti-symmetric matrix (or skew-symmetric J_ij = −J_ji). These cases are much simpler to analyze, because J is a normal matrix so p_C(x) can be derived from the well known eigenvalue distribution of J [31]. For symmetric random connectivity, (14)

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

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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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_ij = −J_ji. The pdf diverges at x = 1 as .

https://doi.org/10.1371/journal.pcbi.1010327.g003

Quantitatively using the equations above, we see that p_C(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), (16)

The p_C(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_ij, J_ji), −1 < κ < 1, we have derived an implicit equation for p_C,g,κ(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_c = (1 + κ)⁻¹ (based on the spectrum of J [31]). As κ increases, the relative coupling strengthg_r = g/g_c = 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_C,g,κ(x) with the same g_r to study effects due to κ beyond changing g_r. As shown in Fig 4D, when fixing the relative strength g_r the distribution narrows as κ increases.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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_r = g/g_c (Eq (18)). For small g_r, 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_r = 0.4 to control the main effect (see text). The changes in shape are now smaller and the support narrows with increasing κ.

https://doi.org/10.1371/journal.pcbi.1010327.g004

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)) (17)

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_C,g,κ(x) with reciprocal motifs are also reflected by the dimension measure, for which we derived a closed-form expression (Section E in S1 Text) (18)

Here μ₁ 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_c = (1 + κ)⁻¹, . The numerator of μ₁ is at least 2θ − 1 ≥ 2/(1 + κ) − 1 > 0. Therefore μ₁ diverges as . Since 1 + 4(g² − θ) ≥ 1 + 4(g² − g) ≥ (1 − 2/(1 + κ))² > 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_r (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_C,g,κ(x) will become similar to that of the anti-symmetric network p_{C,g,κ = −1}(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_C,g,κ(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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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^T. , 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^T 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).

https://doi.org/10.1371/journal.pcbi.1010327.g005

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_ij 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_αβ = k_αβK 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_ij. Importantly, the choice of connection probabilities and weights ensures that var(J_ij) 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_E = w₀ and w_I = −w₀. Thus, we can write J_ij 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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 = w₀ = 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.

https://doi.org/10.1371/journal.pcbi.1010327.g006

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) (19) where is the Fourier transform of the neural activity and z^† is the complex conjugate. C_ij(ω) can also be calculated by the Fourier transform of the time-lagged cross-correlation functions C_ij(τ) = 〈x_i(t)x_j(t − τ)〉 (Wiener-Khinchin theorem). Analogous to Eq (2)C(ω) obeys [25], (20)

Here z^† is the complex conjugate and |z| is the norm for a complex z. The transfer function a(ω) = (1 + iτω)⁻¹ summarizes the dynamics of single neurons in the network and corresponds to a response filter of e^−t/τ/τ, 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)) (21)

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(ω) ∝ ω². 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_C(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 μ_n is the n-th moments of the eigenvalue distribution, and the counterpart for the sampled distribution, (22)

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 (23)

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_C(x) (Eq 7) still holds under time sampling for any fixed α as g → 1 (Section H.3 in S1 Text).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

https://doi.org/10.1371/journal.pcbi.1010327.g007

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_s = 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) (24)

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 L² 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 σ² 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, σ² 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_R(x) = μp_C(μ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_R(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_C(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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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.

https://doi.org/10.1371/journal.pcbi.1010327.g008

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_i(t) to firing rates r_i(t) = ϕ(h_i(t)), and the recurrent interaction term is now (Eq (33) in Methods). When the connection strength g = Nvar(J_ij) 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_i(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].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Covariance spectrum in nonlinear dynamics.

A. Top: A histogram of covariance eigenvalues calculated from firing rate activities r_i(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_i(t) (aggregated across i and t) and the orange dashed curve is the histogram of r_i(t). The overlaying red curve shows the nonlinearity ϕ(x) as a reference. The 〈⋅〉 in the title of each plot represents averaging ϕ′(h_i(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.

https://doi.org/10.1371/journal.pcbi.1010327.g009

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).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

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_i−1 = 0.4, J_i+1,i = 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).

https://doi.org/10.1371/journal.pcbi.1010327.g010

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.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_ij, J_ji) (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_ii 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_ij. In the large-network limit, since individual connections are weak (e.g., ), allowing this self-coupling or setting J_ii = 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.

5.3 Power-law approximation of the eigenvalue distribution

The power-law property of p_C(x) for i.i.d. J_ij 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_C(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 − g²) or note that (1 − g²)³x₊x₋ = 1. Consider a x that is far away from the support edges as g → 1, given the above, this means, (34)

Note that since x₊/x₋ ∼ (1 − g²)⁻³, 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_g,α(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 (35)

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.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_i(t) (assuming bin size is ΔT large), (37)

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_i(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 L² norm cdf error using the Cramer-von Mises statistic (38)

Here n is the number of samples and x_i are the i-th empirical eigenvalues. Alternatively, the error can also be measured under L^∞ norm based on the Kolmogorov-Smirnov statistic (39) where x_i 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.