Properties of networks with partially structured and partially random connectivity

Abstract

Networks studied in many disciplines, including neuroscience and mathematical biology, have connectivity that may be stochastic about some underlying mean connectivity represented by a non-normal matrix. Furthermore, the stochasticity may not be independent and identically distributed (iid) across elements of the connectivity matrix. More generally, the problem of understanding the behavior of stochastic matrices with nontrivial mean structure and correlations arises in many settings. We address this by characterizing large random N×N matrices of the form A=M+LJR, where M,L, and R are arbitrary deterministic matrices and J is a random matrix of zero-mean iid elements. M can be non-normal, and L and R allow correlations that have separable dependence on row and column indices. We first provide a general formula for the eigenvalue density of A. For A non-normal, the eigenvalues do not suffice to specify the dynamics induced by A, so we also provide general formulas for the transient evolution of the magnitude of activity and frequency power spectrum in an N-dimensional linear dynamical system with a coupling matrix given by A. These quantities can also be thought of as characterizing the stability and the magnitude of the linear response of a nonlinear network to small perturbations about a fixed point. We derive these formulas and work them out analytically for some examples of M,L, and R motivated by neurobiological models. We also argue that the persistence as N→∞ of a finite number of randomly distributed outlying eigenvalues outside the support of the eigenvalue density of A, as previously observed, arises in regions of the complex plane Ω where there are nonzero singular values of L(-1)(z1-M)R(-1) (for z∈Ω) that vanish as N→∞. When such singular values do not exist and L and R are equal to the identity, there is a correspondence in the normalized Frobenius norm (but not in the operator norm) between the support of the spectrum of A for J of norm σ and the σ pseudospectrum of M.

FIG. 1

(Color online) Top panel: the total power spectrum of steady state response $\overline{\parallel {\mathbf{x}}_{\omega }{\parallel }^2}$ as a function of input frequency ω, Eq. (2.33), for the system Eq. (2.2) with A = M + σJ, and M given by Eq. (2.35) with w = 1 and λ_n = ±i (with +i and −i alternating), respectively. Here, N = 700, σ = 0.5, and γ = 0.8. The input was fed into the last component of x (the beginning of the feedforward chain characterized by Eq. (2.35)), which for the matrix M has natural frequency -1. That is, the input was ${\mathbf{I}}_0\sqrt{2}\cos \omega t$ where I₀ was 1 for the last component and 0 for all other components. The green (thick dashed) curve is the ensemble average of the total power spectrum, ${\langle \overline{\parallel {\mathbf{x}}_w{\parallel }^2}\rangle }_J$, calculated numerically using the general formula Eq. (2.33), which is compared with an empirical average over 100 realizations of real Gaussian J (solid red line, mostly covered by the dashed green line). The pink (light gray) area shows the standard deviation among these 100 realizations around this average. The blue (thin) line shows the result when disorder, σJ, is ignored, i.e. A is replaced by its ensemble average M. Bottom panel: the eigenvalue spectrum of M + σJ (black dots). Red big dots at ±i show the eigenvalues of M. The red curve is the outer boundary of the eigenvalue spectrum of A as computed numerically using Eq. (2.5). The real and imaginary axes of the complex plane are interchanged, so that the frequency axis in the top panel can be matched with the imaginary part of the eigenvalues, i.e. the natural frequencies of Eq. (2.2).

FIG. 2

(Color online) The eigenvalue spectra of A = M + σJ for N = 2000 and M given by Eq. (2.35) with λ_n = 0, w = 1 for single realizations of real Gaussian J. σ = 0.95 and 0.5 in the left and rights panels, respectively. The red circles mark the circular boundaries of the spectral support given by Eq. (2.36). The insets show a comparison of the analytic formula Eq. (2.37) for the spectral density (black smooth trace) and histograms corresponding to the particular realization shown in the main plot (red jagged trace).

FIG. 3

(Color online) The norm squared of the response to impulse, ∥x(t)∥², of the system Eq. (2.2), for A = M + σJ, with binary J, and M given by Eq. (2.35) (with λ_n = 0) describing a N-long feedforward chain with uniform weights w. Here, w = 1, σ = 0.5, $\gamma =1.005\sqrt{{\sigma }^2+w^2}\simeq 1.124$, and N = 700. The green (thick dashed) curve shows our result, Eq. (2.38), for the average squared impulse response, 〈∥x(t)∥²>_J, which lies on top of the red (thick solid) curve showing the empirical average of ∥x(t)∥² over 100 realizations of binary J. The five thin dashed black curves show the result for five particular realizations of J, and the pink (light gray) area shows the standard deviation among the 100 realizations. The standard deviation shrinks to zero as N → ∞, and ∥x(t)∥² for any realization lies close to its average for large N. For comparison the purple (thin, lowest) curve shows ∥x(t)∥² obtained by ignoring the effect of quenched disorder, i.e. by setting A = M.

FIG. 4

(Color online) The eigenvalue spectra of A = M + σJ for a binary J with σ = 0.1 and M given by Eq. (2.41) with K = 1 (corresponding to w_b = 1 for all the diagonal 2 × 2 blocks in Eq. (2.40)). The main panels show the eigenvalues for single realizations of J, with N = 600 (left) and N = 60 (right). The red circles mark the boundaries of the spectral support, Eq. (2.44). Since A is real in this case, its eigenvalues are either exactly real, or come in complex conjugate pairs; the spectrum is symmetric under reflections about the real axis. However, such signatures of the reality of the matrix appear only as subleading corrections to the spectral density ρ(z); they are finite size effects which vanish as N → ∞. The insets show a comparison of the analytic formula Eq. (2.45) (black curve) and the empirical result, based on the eigenvalues of the realizations in the main panels, for the proportion, n_< (r), of eigenvalues lying within a radius r of the origin (red dots). The random fluctuations and the average bias of the empirical n_< (r) are both already small for N = 60, and negligible for N = 600.

FIG. 5

(Color online) The eigenvalue spectra of A = M + σJ for the M given by Eq. (2.42) in the balanced case, v^T u = 0. Here, N = 800, σ = 1 and µ = 12 (see equation Eq. (2.43)). The black dots are the superimposed eigenvalues of A for 20 different realizations of complex Gaussian J. The small red circle enclosing the vast majority of the eigenvalues has radius σ = 1, corresponding to the standard circular law Eq. (2.46). A Θ(N) number of eigenvalues lie within this circle. $Aϴ\left(\sqrt{N}\right)$ number lie just outside of this circle in a thin boundary layer which shrinks to zero as N → ∞. Finally, a Θ(1) number of eigenvalues lie at macroscopic distances outside the unit circle. ^{The dashed blue circle shows radius r}0 ^{given by Eq. (2.44);} outliers can even lie outside this boundary.

FIG. 6

(Color online) The number of eigenvalues of M + σJ, for the M given by Eq. (2.42), lying outside the circle of radius σ vs. N (red line). Here, σ = 1, µ = 12 and v^T u = 0. The numbers (red points connected by solid red lines) are obtained by numerically calculating the eigenvalues and counting the outliers for 200 realizations of J, and taking the average of the counts over all realizations, for N = 100, 200, 400, 800, 1600 (error bars show standard error of mean). The black dashed line plots $\sqrt{N}$ for comparison with our theoretical result Eq. (2.47); the (dashed) blue line which includes subleading corrections to $\sqrt{N}$, is obtained by numerically solving Eq. (5.42) and substituting the result in Eq. (5.43) (these formulae are in turn obtained from Eqs. (2.8)–(2.9) in Sec. V B).

FIG. 7

(Color online) The number, N_> (r), of outlier eigenvalues of A = M + σJ, for the M given by Eq. (2.42), lying farther from the origin than r, as a function of r. Here, σ = 1, µ = 12 and v^T u = 0. The vertical line marks |z| = r₀ ≃ 3.54 where r₀ is given by Eq. (2.44). The colored (shades of gray) connected points are N_> (r) for realizations of A, based on 200 samples of J, each color for a different N, for N = 100, 200, 400, 800, 1600 and 3200 (error bars show standard error of sample mean). Note the lack of scaling of N_> (r) with N.

FIG. 8

(Color online) The squared norm of response to impulse, x(t) ², of the system Eq. (2.2), for A = M + σJ, with log-normal J, and M given by Eq. (2.40) describing N/2 doublet feedforward chains weights w_b. Here, w_a = |w_b|^{2 b} = 3, σ = 0.4, γ = 1, and N = 1400. The green (thick dashed) curve shows our result, Eqs. (2.48)–(2.49), for the average norm squared which, except for a small window around its peak, lies on top of the red (thick solid) curve showing the empirical average of x(t) ² over 100 realizations of binary J. The five thin dashed black curves show the result for five particular realizations of J, and the pink (light gray) area shows the standard deviation among the 100 realizations. The standard deviation shrinks to zero as N → ∞ and x(t) ² for any realization lies close to its average for large N. For comparison the purple (thin, lowest) curve shows x(t) ² obtained by ignoring the effect of quenched disorder, i.e. by setting A = M.

FIG. 9

(Color online) The eigenvalue spectra of A = M + LJR with M, L and R given by Eqs. (2.52)–(2.55) with neurons belonging to one of three different types (C = 3). The main panels show the eigenvalues for two particular realizations of J. In both panels, N = 2000, f₁ = 0.6, f₂ = f₃ = 0.2, l_c = 1, σ₁ = r₁ = 0.76, σ₂ = r₂ = −0.57, σ₃ = r₃ = −1.71 (so 〈σ_c〉_c = 0 and $r_0^2={\langle {\sigma }_c^2\rangle }_c$ = 1), and J_ij had real entries with log-normal distribution; in the left (right) panel, the normally distributed log₁₀ J_ij had standard deviation 0.5 (0.75). The solid red circles mark the boundaries of the spectral support as given by Eq. (2.60), and the dashed blue circles show the radii given by Eq. (2.61). The insets compare n_>(r) based on the numerically calculated eigenvalues shown in the main panels (connected red dots), with that found by solving Eq. (2.62) (black curve). In the right panel’s inset we have also plotted (green connected dots lying slightly above the red connected circles) the empirically calculated n_>(r) for a single realization with the same ensemble parameters, but with N = 8000; the convergence to the universal limit at N → ∞ is significantly slower in the right panel in which the distribution of J_ij had a considerably heavier tail.

FIG. 10

(Color online) The eigenvalues (black dots) of A = M + J R, with M and R given by Eqs. (2.75)–(2.76) with g = 0.01, a = 1.02 and N = 2000. This matrix governs the dynamics of small perturbations away from a non-trivial random fixed point in a clustered network of neurons (see Eq. (2.74)), studied in Ref. [48]. The cyan dots on the real line are the eigenvalues of M, and the red curve is the boundary of support of the eigenvalue distribution, as calculated numerically from Eq. (2.5).

FIG. 11

The first two lines define different elements of Feynman diagrams: the Green’s function for J = 0 (zero disorder), ${\mathit{G}}_{ab}^{\alpha \beta }(\eta ,z;0)$, the covariance of two J elements, the ensemble averaged Green’s function, G(η, z) ≡ G(η, z; J) _J, and the self-energy Σ(η, z), Eq. (3.24) (the matrix indices for G(η, z) and Σ(η, z) are arranged as for ${\mathit{G}}_{ab}^{\alpha \beta }(\eta ,z;0)$). The third line is the diagrammatic representation of the expansion Eq. (3.22) of G(η, z; J) before averaging over J, where the J’s are represented by dashed lines. Averaging over Eq. (3.2) is performed by pairing all J ’s and connecting them with the wavy lines representing 〈J J〉. In the large N limit, the contribution of crossing pairings is suppressed by negative powers of N ; the sum of all non-crossing diagrams, shown on the fourth line, yields the leading contribution to G(η, z) for large N. The last line shows the diagrammatic representation of Eq. (3.23), which if iterated generates all the non-crossing diagrams. Alternatively, G(η, z) can be found by solving this self-consistent equation directly.

FIG. 12

Contribtutions to Eq. (4.17) in the non-crossing approximation. The first line shows Eq. (4.17) written using the expansion Eq. (3.22). The diagram shows the contribution of the m-th and n-th terms in the expansion for two Green’s functions, respectively. Thus the top (bottom) solid line contains m (n) factors of J, shown by dashed lines. In the large N limit, averaging each summand over J boils down to summing all non-crossing pairings (NCP) of the dashed lines. The second row shows a specific non-crossing pairing for the diagram shown in the first line. Finally, summing over all m, and n and all NCP’s, is equivalent to replacing all solid lines (representing G(η_i, z_i; J = 0)) with thick solid lines representing the non-crossing average Green’s function, G(η_i, z_i) (calculated according to Eqs. (3.28)–(3.26)), and summing over all NCP’s with every pairing connecting the straight lines on top and bottom (and not each to itself). This procedure yields the ladder diagrams, the sum over which is shown in the third line.

FIG. 13

The first row is the diagrammatic representation of Eqs. (4.19)–(4.21). In the last term, ρ’s and λ’s are summed over. It shows the sum of all ladder diagram contributing to Eq. (4.17) (i.e. the last line of Fig. 12) in terms of D, which is defined in the second row. The first term on the right side of the first row equation (the ladder with zero rungs) is the disconnected average Eq. (4.20); it corresponds to taking the average of each Green’s function in Eq. (4.17) separately and then multiplying. The last row shows an iterative form of the equation in the second row, which can be solved to give the expression Eqs. (4.23) and (4.26) for D.

FIG. 14

The contribution to $D_{ad;cb}^{\alpha \delta ;\gamma \beta }$ (1; 2) from the second term in the series shown in the second row of Fig. 13, in more detail. The covariance of J in the form Eq. (4.22) is used to write this expression in a more manageable form. The repeated indices, r, t, u, s, are summed over 1 and 2. The matrices inside the loop multiply each other in cyclic order, giving rise to the trace Tr (G(2)π^tG(1)π^u). The whole diagram gives $\frac{1}{N}{\sum }_{rs}{({\pi }^r\otimes 1)}_{ad}^{\alpha \delta }{\left[{\sigma }^1{\Pi }^D{\sigma }^1\right]}_{rs}{({\pi }^s\otimes 1)}_{cb}^{\gamma \beta }$ where the “polarization matrix” ${\Pi }_{tu}^D$ was defined in Eq. (4.25).

FIG. 15

(Color online) The orbits (shown by thin red paths) for two diagrams for the spectral density in a complex J ensemble. The non-crossing diagram on top has three orbits: orbit (1) is the external orbit connecting the two ends of the Green’s function, while orbits (2) and (3) are the internal orbits. As in Eqs. (A1) and (A2), they contribute $\mathrm{tr}\left({\sigma }^{+}\mathit{G}(\eta ,z;0){\pi }^{r_1}\mathit{G}(\eta ,z;0)\right),\mathrm{Tr}\left({\pi }^{3-r_1}G(\eta ,z;0){\pi }^{r_2}\mathit{G}(\eta ,z;0)\right)$ and $\mathrm{Tr}\left({\pi }^{3-r_2}\mathit{G}(\eta ,z;0)\right)$ respectively, with r₁ and r₂ summed over 1 and 2 (cf. Eq. (3.21)). The trace contributed by each of the three orbits is O(N ), which when combined with the three factors of 1/N accounting for the two wavy lines and the normalization of the external orbit’s trace, yield an O(1) expression for this diagram. By contrast, the crossing diagram on the right has no internal orbits. Its only external orbit contributes $\mathrm{Tr}\left({\sigma }^{+}\mathit{G}(\eta ,z;0){\pi }^{r_2}\mathit{G}(\eta ,z;0){\pi }^{3-r_1}\mathit{G}(\eta ,z;0){\pi }^{3-r_2}\mathit{G}(\eta ,z;0){\pi }^{r_1}\mathit{G}(\eta ,z;0)\right)$ which after normalization is O(1). Accounting for two factors of 1/N coming from the wavy lines, we then see that this crossing diagram is O(N ⁻²) and hence is suppressed as N → ∞.