Bridging Scales: Spectral Theory Reveals How Local Connectivity Rules Sculpt Global Neural Dynamics in Spatially Extended Networks

Bridging Scales: Spectral Theory Reveals How Local Connectivity Rules Sculpt Global Neural Dynamics in Spatially Extended Networks | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Bridging Scales: Spectral Theory Reveals How Local Connectivity Rules Sculpt Global Neural Dynamics in Spatially Extended Networks View ORCID Profile Yuhan Huang , View ORCID Profile Keren Gao , View ORCID Profile Dongping Yang , View ORCID Profile Sen Song , View ORCID Profile Guozhang Chen doi: https://doi.org/10.1101/2025.05.23.655871 Yuhan Huang 1 School of Computer Science, Peking University , China 2 School of Information Science and Technology, University of Tokyo , Japan Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Yuhan Huang Keren Gao 1 School of Computer Science, Peking University , China Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Keren Gao Dongping Yang 3 Zhejiang Lab , China Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Dongping Yang Sen Song 4 Department of Biomedical Engineering, Tsinghua University , China Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Sen Song Guozhang Chen 1 School of Computer Science, Peking University , China Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Guozhang Chen For correspondence: guozhang.chen{at}pku.edu.cn Abstract Full Text Info/History Metrics Preview PDF Abstract The brain’s diverse spatiotemporal activity patterns are fundamental to cognition and consciousness, yet how these macroscopic dynamics emerge from microscopic neural circuitry remains a critical challenge. We address this by developing a spatially extended neural network model integrated with a spectral theory of its connectivity matrix. Our theory quantitatively demonstrates how local structural parameters, such as E/I neuron projection ranges, connection strengths, and density determine distinct features of the eigenvalue spectrum, specifically outlier eigenvalues and a bulk disk. These spectral signatures, in turn, precisely predict the network’s emergent global dynamical regime, encompassing asynchronous states, synchronous states, oscillations, localized activity bumps, traveling waves, and chaos. Motivated by observations of shifting cortical dynamics in mice across arousal states, our framework not only explains this repertoire of behaviors but also offers a principled approach to inferring underlying effective connectivity changes from macroscopic brain activity. By mechanistically linking neural structure to dynamics, this work provides a powerful tool for understanding brain function and paves the way for identifying potential biomarkers for neurological disorders. 1 Introduction The brain’s activity is remarkably diverse, forming complex spatiotemporal patterns that vary with an organism’s cognitive state and level of consciousness [ 1 – 4 ]. Propagating waves of neural activity, for example, are observed across numerous brain regions and species, from insects to humans [ 5 , 6 ]. Recent high-resolution imaging in mice, as they transition from anesthesia to wakefulness, has highlighted this complexity, revealing large-scale cortical waves during anesthesia and a shift towards more localized, intricate spatiotemporal patterns upon awakening [ 4 ]. Such observations emphasize that the spatial organization of neural activity is a fundamental aspect of brain function [ 7 – 10 ]. Consequently, a key challenge is to understand how these macroscopic dynamical states, which distinguish brain states, arise from the underlying microscopic neural architecture and effective connectivity rules. Since many critical brain dynamics are inherently spatial, involving coordinated activity across neural tissue, models must incorporate spatial extent to capture these phenomena faithfully [ 11 ]. While various computational models have shed light on specific neural phenomena [ 12 – 20 ], a comprehensive framework that quantitatively links a broad range of local network structures to the full spectrum of emergent global dynamics in spatially extended systems is still developing. Current theories often explain particular aspects, but a general theory predicting the emergence of, and transitions between, diverse spatiotemporal patterns from fundamental structural properties across a wide parameter space remains a critical need. To address this, we present a unifying theoretical framework centered on a spatially extended recurrent neural network model with excitatory (E) and inhibitory (I) populations. The core of our approach is a spectral analysis of the network’s connectivity matrix. Building upon theories of random matrices [ 21 – 25 ], we provide an analytical formulation for the spectral bulk arising from connection heterogeneity and, critically, demonstrate that the inherent spatial organization of the network itself constitutes a low-rank structure. This allows us to characterize a rich set of outlier eigenvalues that extend beyond a simple dependence on neuron types, reflecting specific spatiotemporal modes determined by the network’s geometry. We show that this comprehensive “spectral blueprint” - encompassing both the derived bulk and the spatially determined outliers - quantitatively links key local structural parameters (spatial reach of E/I connections, relative strengths, local density, weight variability) to a full repertoire of global dynamical phases. This framework offers predictive power, enabling us to anticipate the network’s spatiotemporal behavior from its effective connectivity structure. These predicted phases include stable asynchronous states, global synchrony, oscillations, localized bumps, traveling waves, and chaotic dynamics. Our work thus provides a principled understanding of how network structure dictates emergent function. 2 Spatial Extended Neural Networks 2.1 Model Definition To explore the link between structure, dynamics, and the eigenvalue spectrum, we constructed a minimal rate-based neural network model with biologically realistic features. The network consists of excitatory and inhibitory neurons, sparse connectivity that decays with distance, and synaptic weights drawn from a Gaussian distribution ( Fig. 1 ). Download figure Open in new tab Figure 1. (a) and (b) Mesoscopic optical imaging of the anesthetized and wake mouse cortex [ 4 ]. Model schematic. (c) Spatially extended neural network with excitatory (E) and inhibitory (I) neurons embedded in 2D space. (d) Synaptic weights follow a distance-dependent Gaussian distribution. (e) Connection probability follow a distance-dependent wrapped Gaussian distribution. We consider a network consisting of N E excitatory neurons and N I inhibitory neurons, where N E : N I = 4 : 1. Both excitatory neurons and inhibitory neurons are evenly distributed in the region [0, 1) × [0, 1). The dynamics of the neurons are described by the following equation: Here, h indicates the membrane potential of neurons, ϕ (·) is the activation function of neurons, α, β ∈ { E, I } denotes the neuronal populations, and ξ α represents the external input received by the population α . The external input is modeled as independent white noise with intensity ξ 0 = 0.1. The connections between neurons are sparse. The probability of connection between a α neuron located at x i and a β neuron located at x j is given by, where This indicates that the probability of connection between neurons decays according to a wrapped Gaussian profile with a characteristic decay length of d αβ . Here, represents the average out-degree, namely the number of α neurons to which a single β neuron projects. Two neurons are connected with a probability of . If two neurons are indeed connected, the connection weight follows an independent Gaussian distribution with a mean of and a variance of , namely, . The activation function ϕ (·) can be arbitrary. In our experiment, we chose the activation function ϕ E ( x ) = tanh( x ) for excitatory neurons and ϕ I ( x ) = 5 tanh( x/ 5) for inhibitory neurons. However, the choice of activation function has little impact on the conclusions, and we discuss this in detail in the Appendix A.5. In the Appendix A.5, we also show the results with the activation function as ϕ E ( x ) = 10[ x ] + and ϕ I ( x ) = 2[ x ] 2 . 2.2 A Repertoire of Emergent Dynamical Phases Asynchronous State The balance between excitatory and inhibitory interactions shapes the network’s synchronization and stability. When inhibition dominates or excitation and inhibition are balanced ( Fig. 2(a) ), the network exhibits low synchronization and small deviations from the steady firing rate. We define this regime as the asynchronous state ( Fig. 2(c) ). In the asynchronous state, neuronal membrane potentials fluctuate around the fixed point (zero) due to external input. Neural activity remains weakly correlated across neurons, and no spatial patterns emerge. This regime aligns with the classical asynchronous state described by [ 13 ]. Download figure Open in new tab Figure 2. Dynamical regimes in spatially extended networks. (a,b) Phase diagrams under different parameters (see Appendix A.7). Wave-chaos phase boundaries determined as in Appendix A.11. (c) Representative neural activity showing temporal evolution (left) and spatial patterns (right) with shared colorbar; enhanced versions in Appendix A.9 show the asynchronous state with individual scaling. Download figure Open in new tab Figure 3. Spectral signatures and spatial patterns of network dynamics. (a) Eigenvalue spectra and corresponding eigenvectors across spatial modes. (b) Example eigenvalue distributions (blue triangles: eigenvalues with the largest real part), leading eigenvectors, neural activity, and order parameters (blue point: data points) for distinct dynamical states. The four order parameters (rightmost column) indicate, from left to right: neural activity fluctuations, local synchrony, spatial patterning, and oscillations. Detailed definitions are given in the Appendix A.3. Synchronous State When excitation dominates ( Fig. 2(a) ), the network can enter a highly synchronized regime with large deviations from the steady-state firing rate, which we term the synchronous state phase. In this phase, the firing rates of all neurons shift away from the fixed point, and strong global synchronization emerges across the network ( Fig. 2(c) ). This globally synchronized activity resembles pathological brain states such as epileptic seizures [ 26 ]. Oscillatory Phase Neural oscillations emerge when excitation-inhibition coupling strengthens in highly synchronized networks ( Fig. 2(a) ). We term this dynamical regime the oscillatory phase, characterized by periodic membrane potential changes and synchronized neural activity ( Fig. 2(c) ). Our identified excitation-inhibition loop mechanism aligns with established experimental and theoretical work. This consistency supports the biological plausibility of our model, as evidenced by prior studies on neural oscillations [ 1 , 27 ]. Localized Bump Phase Mismatched projection ranges between neuron types generate structured spatial activity patterns. When excitatory neurons project locally while inhibitory neurons extend farther, the system exhibits bump phases and wave phases ( Fig. 2(a) ). The bump phase produces localized spatial patterns dependent on activation functions. Neural activity synchronizes within discrete regions, forming stripes ( Fig. 2(c) ). Rectified linear and threshold power-law functions yield spot-like patterns (see Appendix A.5), aligning with previous simulation results of spiking neural networks[ 14 , 16 ]. Similar wavelength-specific patterns occur in juvenile visual cortex [ 28 ], with persistent long-range correlations in adulthood [ 29 ]. Wave Phase The wave phase features propagating oscillations with location-dependent phase shifts ( Fig. 2(c) ). This mirrors biological observations of traveling alpha/gamma waves [ 6 , 5 ], suggesting our model captures essential mechanisms of spatial-temporal dynamics. Chaotic Phase High connection sparsity and weight variance drive neural activity into a chaotic regime ( Fig. 2(b) and Fig. 7(a) ). We term this disordered state the chaos phase, where spatial patterns collapse due to uncorrelated neural firing. Neurons exhibit large-amplitude, weakly correlated fluctuations in the chaos phase ( Fig. 2(c) ). Membrane potentials vary erratically, matching the second type of asynchronous state described in spiking networks [ 13 ]. This aligns with dynamic mean-field theory predictions of chaotic dynamics in rate-based neural networks [ 30 ]. 3 The Spectral Blueprint: Decoding Dynamics from Connectivity Structure 3.1 Effective Connectivity and its Eigenvalues Building on the neuroscience perspective of effective connectivity (EC) as a model-based measure of directed interactions [ 31 ], we propose that the linearized Jacobian matrix − I + J ϕ ′ ( x ∗ )–derived from the recurrent neural network’s fixed-point dynamics–serves as an analogous EC matrix. Here, structural connectivity J is dynamically modulated by nonlinear gains ϕ ′ ( x ∗ ), mirroring how anatomical constraints and state-dependent plasticity jointly shape brain network interactions. To unravel how such connectivity shapes collective behavior, we analyze its eigenvalue spectrum, which governs stability and activity patterns. In the following sections, we employ random matrix theory to characterize this spectrum, revealing universal dynamical regimes emergent from its structure. In the main text, we use the tanh activation function, so the connectivity matrix J can be regarded as the effective interaction matrix. We also present results with alternative activation functions in the appendix A.5 and A.8. 3.2 Eigenvalue Spectrum of Spatially Extended Networks The eigenvalue spectrum of spatially distributed neural networks’ connectivity matrix J comprises two distinct components: a bulk disk region and a set of spectral outliers, which respectively reflect heterogeneous neuron interactions and population-averaged connectivity patterns. Following the matrix decomposition , where denotes the expectation matrix and δ J represents zeromean random fluctuations, we observe that δ J governs the bulk spectrum through its stochastic components while generates spectral outliers through its low-rank structure. This aligns with the perturbation framework established by [ 22 ], wherein low-rank modifications to random matrices predominantly affect outlier positioning while preserving the original bulk spectral radius. The dichotomy between these spectral components provides a mathematical characterization of neural network dynamics - the deterministic outliers capture macroscopic interaction features, whereas the bulk spectrum encodes microscopic connection variability. 3.2.1 Outliers The matrix determines the outlier part of the eigenvalue spectrum [ 22 ]. The expected connectivity between neurons exhibits spatial translational invariance. We can decompose the activities of neurons into different spatial Fourier modes (See Appendix A.10). For each spatial Fourier mode, the effective connectivity matrix is where . The eigenvalue spectrum of the connectivity matrix reveals distinct outliers generated by spatial Fourier modes at different wave vectors k . Each outlier corresponds to collective population dynamics mediated by k -specific interaction submatrices. These submatrices encode spatially modulated couplings between excitatory and inhibitory populations, with effective weights determined by three components: (1) baseline inter-population connectivity ( between populations α and β ), (2) k -dependent modulation of spatial interaction range exp . 3.2.2 Bulk Disk The matrix δ J governs the bulk disk part of the eigenvalue spectrum, of which the radius is determined by the local sparsity of neuronal connections and the variance of the weights. Analytical tools developed in [ 23 , 24 ] characterize the spectral distribution of random matrices with independent, zero-mean, finite-variance entries. Applying these results, we can determine that the eigenvalues are extended within a circle of a specific radius , where the elements of the matrix M represent the variances of the elements of the connectivity matrix J . The variance of the elements of the connectivity matrix J is given by, The eigenvalues of this matrix are equivalent to those of the reduced matrix, where the elements represent the heterogeneity of connections between different types of neurons. For the one-dimensional case, this reduced matrix is Based on the above equations, the radius of the circular part of the eigenvalue spectrum is given by, We can observe that the heterogeneous matrix M αβ is composed of two parts: sparsity and variability in synaptic weights. The left half of Equation 5 represents the heterogeneity brought about by the sparse connections between neurons, while the right half of Equation 5 represents the heterogeneity due to the variability in the weights of the neuronal connections. As the local sparsity of neuronal connections and the variance of weights increase, the radius of the bulk disk part of the eigenvalue spectrum also increases. 3.3 Linking Spectral Features to Dynamical Phases The eigenvalue spectrum of the effective connectivity matrix J ϕ ′ ( x ∗ ) governs network dynamics: the dominant eigenvalue (largest real part) determines stability near the fixed point. A real part exceeding the critical threshold ( Re ( λ dom ) ≥ 1) quantifies deviation magnitude, while the dominant eigenvector specifies the spatial activity pattern. We thus classify neural networks into distinct dynamical phases based on their eigenvalue spectrum, as summarized in the table 1 . (See Appendix A.8 for details), View this table: View inline View popup Download powerpoint Table 1: Dynamical Phases and Spectral Features 3.4 The Role of Key Structural Parameters in Shaping the Spectrum Excitaion-Inhibition Balance The magnitude of excitatory and inhibitory interaction affects the magnitude of the real parts of the outlier eigenvalues, which in turn influences the degree of deviation of neural activity from the fixed point and the level of synchronization. As shown in Eq. 3 , the magnitude of excitatory and inhibitory interactions influences the elements of the effective interaction matrix. The more excitatory interaction there is, the larger the real part of the outlier eigenvalues, the greater the degree to which neural activity deviates from the fixed point, and the more synchronized the neural activity becomes; the opposite is also true. Excitation-Inhibition Loop The Excitation-Inhibition Loop is considered a crucial component for the emergence of neural oscillations [ 1 , 27 ]. Our theory explains this from the perspective of the eigenvalue spectrum. The greater the magnitude of the interaction between excitatory and inhibitory neurons, the more likely it is for the outlier eigenvalues to have an imaginary part, causing neural oscillations to emerge( Eq. 3 ). This is because and are located on the off-diagonal elements of the effective interaction matrix, and an increase in and can lead to the appearance of an imaginary part in the outliers. The mismatch of Excitation/Inhibition projection range The mismatch of projection range d αβ of types of neurons causes the emergence of spatial patterns. As shown in Eq. 3 , the elements of the effective matrix contain decay factors exp . If all the projection ranges d αβ are the same and let’s denote d αβ = d , the eigenvalues would follow the relationship as λ ( k ) = λ (0) exp (−|| k || 2 d 2 /2). In this case, only the eigenvalues corresponding to the wave vector k = 0 can have the largest real part. However, if the projection range d αβ of types of neurons is mismatched, eigenvalues corresponding to wave vector k ≠ 0 can have the largest real part. A common way for a spatial pattern to emerge is when the projection range of inhibitory neurons is greater than that of excitatory neurons, which is known as lateral inhibition [ 29 , 28 ]. Our theory explains this from the perspective of the eigenvalue spectrum. This is because the elements of the effective matrix contain decay factors , and the elements decay faster with k if d αβ is larger. Therefore, the inhibitory elements of the effective matrix may decay faster than excitatory ones. In some wave vectors k , the excitation exceeds the inhibition and the real part of these eigenvalues may exceed 1 and spatial patterns of neural activity emerge. Local sparsity Local sparsity, rather than the overall sparsity, plays a significant role in chaotic neural activity. In a spatially extended neural network, the number of neurons connected to a given neuron is certainly sparse compared to the total number of neurons. However, Eq. 5 indicates that what truly determines the dynamics is the ratio of the number of connections a neuron has to the number of neurons within its projection range , namely the “local sparsity” that plays a role in the radius of the bulk disk part. This is different from the situation without spatial distribution. This suggests that even under globally sparse conditions, the relative density of local connections between neurons in the brain may be the reason it can generate synchronized activities such as traveling waves. Besides, the concept of “local sparsity” also demonstrates many phenomena in numerical simulation, including both spike-based and rate-based models. In most numerical simulations that produce neural activity with spatial patterns, the connection between neurons is often very dense, sometimes even with full connection locally [ 8 , 9 , 19 ], while numerical simulations that produce neural activity without spatial patterns are often sparse [ 32 ]. 3.5 Analyzing Dynamical Transitions To elucidate phase transitions between distinct dynamical phases, we analyze the spatiotemporal organization of emergent activity patterns. We derive phase velocity fields from network activity to capture the local direction and speed of patterns like propagating waves. By applying Singular Value Decomposition (SVD) to these velocity fields, we identify dominant spatiotemporal modes of activity flow. This allows us to systematically study how these modes reconfigure as the system traverses different dynamical regimes and critical boundaries between them, offering a quantitative window into the nature of these transitions. Our analysis reveals a key phenomenon at phase boundaries: “mode mixing”, where SVD modes characteristic of both adjacent pure phases significantly contribute, indicating dynamically hybrid states consistent with underlying spectral properties near instability ( Fig. 4(a-b) ). Crucially, the dominant spatial SVD modes identified in our model (e.g., plane waves, spirals) exhibit compelling qualitative similarities to patterns observed in mesoscopic optical imaging of the mouse cortex across different arousal states ( Fig. 4(c-d) ) [ 4 ]. This correspondence suggests our model captures salient principles of spatiotemporal pattern formation and transition relevant to real brain dynamics. See Appendix A.2 for methods and supplementary results. Download figure Open in new tab Figure 4. (a) Phase diagram with SVD analysis on combined time series from marked parameter sets. (b) Projections onto SVD modes for: asynchrony center (blue), wave center (pink), and boundaries (circles/triangles). (c) SVD modes 1,3,4 of combined series. (d) Selected experimental SVD modes. Colored boxes highlight matching phase velocity patterns: green (planar waves), orange and red (clock wise & counter clock wise spirals, respectively). Download figure Open in new tab Figure 5. Analysis of all the phase boundaries and triple points. (a) Asynchrony-bump-wave phase diagram, and projection curves onto SVD modes of the combined time series of the parameter sets marked by crossings in the phase diagram. Sky blue line refers to asynchrony center, green line refers to oscillation center, dark blue line refers to synchrony center. Circles, triangle and squares refers to the corresponding parameter sets marked by numbers in the phase diagram. The inset shows the projection curves onto the first ten modes. (b) Same as (a), pink line refers to wave center, dark red line refers to bump center. (c) Same as (a), orange line refers to chaos center. Download figure Open in new tab Figure 6. Dynamical regimes in spatially extended networks with alternative activation functions. (a) The Phase diagram under alternative activation functions. (b) Representative neural activity and spatial patterns of excitatory neurons in different dynamical phases over time. Download figure Open in new tab Figure 7. Phase diagrams and order parameters under different parameters. 4 Conclusions and Discussions We introduced a spectral theory for spatially extended neural networks, quantitatively linking local connectivity (E/I projection ranges, strengths, and local density) to global dynamics via the eigenvalue spectrum of the connectivity matrix. This spectral blueprint, characterized by outlier modes and a bulk disk, accurately predicts a rich repertoire of emergent behaviors including asynchronous states, oscillations, bumps, waves, and chaos, providing a mechanistic bridge from circuit structure to dynamic function. This framework contextualizes how observed brain dynamics, such as state-dependent patterns [ 4 ], arise. It aligns with neural field theories [ 33 ] regarding pattern formation and offers refined insights into chaos generation compared to globally coupled models [ 30 ], highlighting the role of local density. Crucially, we posit that our model’s static connectivity represents time-varying effective connectivity in the brain, which is constantly reshaped by neuromodulation, stimuli, and attention [ 34 , 35 ]. Thus, the identified dynamical phases are attractor states within a larger phase space through which the brain navigates, enabling cognitive flexibility. This perspective also informs how pathological dynamics in disorders like epilepsy [ 36 , 26 ] or Parkinson’s [ 37 ] might stem from altered effective connectivity, offering avenues for spectral or dynamic biomarkers. Key limitations guide future work: extending the theory to nonlinear neural/node’s dynamics [ 13 ] and incorporating synaptic plasticity to model adaptive spectral changes and learning [ 38 ] are paramount. Addressing structural complexities beyond isotropic connectivity, modeling the explicit time-variance of effective connectivity, and robustly inferring spectral features from empirical data are also critical. Also, a more direct characterization of chaos needs to be done in future works. Furthermore, elucidating the direct computational roles of these spectrally-defined dynamical regimes remains a vital pursuit [ 40 ]. In conclusion, our spectral theory provides an analytically tractable lens for understanding selforganization in spatially extended neural systems. It explains how effective local architecture shapes global function and provides a foundation for investigating neural computation, brain state transitions, and the circuit bases of neurological disorders where dynamic adaptability is compromised. A. Appendixes A.1 Supplementary discussion Our calculation of the bulk disk radius builds upon foundational random matrix theory concerning the circular law, applicable even to non-Gaussian weight distributions [ 21 ], and methods for structured random networks [ 25 , 24 , 23 ]. The emergence of outlier eigenvalues from low-rank perturbations to such random matrices is also well-established [ 22 ]. Our novelty lies in systematically deriving both these spectral components for spatially extended E/I networks and explicitly linking them to distinct dynamical phases. While prior work has explored outlier eigenvalues in non-spatial networks, often focusing on single global outliers [ 41 ], population modes [ 42 ] without considering spatial, or randomly distributed local outliers that can be removed by zero-sum constrain[ 25 , 43 ], and while some studies noted the presence of spatially-organized eigenvalues as a secondary observation without direct calculation or dynamic linkage [ 44 ], our framework uniquely connects the full set of spatially-indexed outlier eigenvalues (for wave vectors) and the bulk disk radius to the emergence of diverse spatiotemporal patterns and chaos, respectively. The explicit calculation of eigenvalues associated with specific wave vectors is crucial for understanding spatially patterned activity. While some models of cortical networks have implicitly or explicitly involved such wave numbers, they often relied on simplifications from neural field theory for analytical tractability [ 18 , 16 , 19 ] or focused on conditions for specific instabilities, such as how E/I imbalance can lead to pattern formation at particular wave numbers [ 14 ]. Our approach provides a more general matrix-based spectral method that directly ties the parameters of the spatially explicit network (projection ranges, strengths) to the entire outlier spectrum without necessarily reducing to a continuum field limit, thereby offering a direct bridge from discrete network structure to emergent spatial dynamics like bumps and waves. A.2 Phase Velocity Field analysis and Wave Pattern Detection Signal Pre-processing for Phase Extraction To analyze wave patterns from time-series data recorded at multiple spatial sites (e.g., from optical imaging or electrophysiological arrays), we first extract the instantaneous phase from each recording site. This process involves two main steps: Band-pass Filtering: The raw signal from each site, denoted as s raw ( t ), is first band-pass filtered to isolate activity within a frequency range of interest. This step serves to remove high-frequency noise and focus the analysis on specific neural oscillations (e.g., delta: 0.5-4 Hz; theta: 4-8 Hz [ 45 ]; alpha: 8-13 Hz) or on frequency bands containing significant signal power as determined by power spectral analysis of the recordings. The choice of frequency band is also constrained by the temporal resolution of the recording technique. Let the filtered signal be x ( t ). Phase Extraction via Hilbert Transform: The instantaneous phase ϕ ( t ) is extracted from the filtered signal x ( t ) at each site using the Hilbert transform. The analytic signal z ( t ) is constructed as, where ℋ{ x ( t )} is the Hilbert transform of x ( t ), A ( t ) is the instantaneous amplitude, and i is the imaginary unit. The instantaneous phase ϕ ( t ) is then obtained as, where arctan2( y, x ) is the two-argument arctangent function that correctly resolves the phase into all four quadrants. The result of this pre-processing is a phase time series ϕ j ( t ) for each spatial recording site j . Optical Flow Method for Phase Velocity Field Estimation To quantify the propagation of phase patterns across the spatial recording array, we employ an optical flow method, specifically the Horn-Schunck algorithm, adapted for phase data. This method estimates a 2D velocity field ( v x ( x, y, t ), v y ( x, y, t )) that describes the motion of surfaces of constant phase. Let I ( x, y, t ) = ϕ ( x, y, t ) represent the instantaneous phase at spatial location ( x, y ) and time t. The core assumption of optical flow is brightness constancy, which for phase translates to phase constancy along a trajectory: I( x + v x dt, y + v y dt, t + dt ) = I ( x, y, t ). A first-order Taylor expansion yields the optical flow constraint equation: where , and are the spatial and temporal partial derivatives of the phase field. To solve for the two unknown velocity components ( v x , v y ) from this single equation, the Horn-Schunck method introduces a global smoothness constraint, minimizing an energy functional E : where α 2 is a regularization parameter that weights the smoothness term. Minimization of this functional leads to a system of coupled partial differential equations for v x and v y . Discretizing these equations (e.g., using finite differences for derivatives and the Laplacian ∇ 2 ) results in a large system of linear equations: where v x and v y now represent values at discrete grid points (x i , y j ). The Laplacian terms and are typically approximated using a five-point stencil, e.g., where is the average of v x at the four cardinal neighbors of (x, y), and Δx = Δy is grid spacing. This system is solved iteratively for each time step t : where (k) denotes the iteration number, and (related to α 2 and grid spacing encapsulates the smoothness constraint. Iterations proceed until convergence or for a fixed number of steps. Implementation Details Partial Derivatives: Spatial derivatives I x , I y were computed from the phase maps ϕ ( x, y , t) at each time t using a Sobel filter or a five-point central difference scheme, averaged between two consecutive time frames t and t + dt . The temporal derivative I t was computed using a forward difference between ϕ ( x, y, t + dt ) and ϕ ( x, y, t ), potentially after local spatial averaging to reduce noise. Boundary Conditions: For calculating spatial derivatives near boundaries, appropriate schemes (e.g. Neumann boundary conditions where derivatives are zero) were applied. For the averaging terms in the iterative update, Neumann or zero-padding (Dirichlet-like) boundary conditions were typically used for the velocity components outside the defined spatial grid. The specific choice of derivative computation and boundary conditions was validated on test data to ensure reasonable velocity fields. Regularization Parameter α: The value of α (or λ s ) was chosen empirically to balance adherence to the optical flow constraint with the smoothness of the resulting velocity field, often by visual inspection of results on synthetic or sample data. This procedure yields a phase velocity field ( v x ( x, y, t ), v y ( x, y, t )) for each time point, which can then be further analyzed to characterize wave properties like direction, speed, and coherence. Supplementary Results of Phase Velocity Field Analysis Our spectral understanding of how network structure dictates distinct dynamical phases can also illuminate the nature of transitions between these phases. To quantitatively characterize these transitions, particularly at phase bound-aries and triple points where multiple dynamical tendencies may coexist or compete, we analyze the spatiotemporal structure of phase velocity fields. Analyzing the instantaneous phase ϕ ( x, y, t ) offers several advantages for understanding organized spatiotemporal patterns. The phase captures the relative timing of oscillatory activity across different spatial locations, making it particularly sensitive to propagating waves, synchronized domains, and their complex interactions. By then computing the optical flow of these phase maps, we obtain a velocity field ( v x ( x, y, t ), v y ( x, y, t )) that directly quantifies the local direction and speed of these emergent patterns. This provides a rich, time-varying representation of the network’s collective spatiotemporal organization, which is amenable to techniques like Singular Value Decomposition (SVD) for identifying dominant modes of activity flow. The methodology for obtaining phase velocity fields from neural activity (either model-generated or experimental) involves band-pass filtering (e.g., 0-40 Hz, based on signal power spectrum), Hilbert transform to extract instantaneous phase ϕ ( x, y, t ), and an optical flow algorithm to compute the velocity vectors of phase propagation. Full details of this pre-processing and optical flow computation are provided in Appendix A.2. To analyze the structure of these time-varying phase velocity fields V field ( x, y, t ), we employ SVD. The velocity field data across all spatial sites (both v x and v y components) and time points is arranged into a matrix X , where rows typically represent time and columns represent flattened spatial components. SVD decomposes this matrix as X = UΣV † , where columns of U are temporal modes, columns of V are spatial modes, and Σ contains the singular values indicating the contribution of each mode. To study phase transitions, we first identify a set of common spatial modes by performing SVD on a combined dataset of phase velocity fields from parameter sets spanning different phases, phase boundaries, and triple points (see markers in Fig. 4a ). Then, for each individual parameter set (numbered markers at a phase boundary in Fig. 4a ), its phase velocity field time series is projected onto these common spatial modes. The resulting projection weights quantify how much each common spatial mode contributes to the dynamics of that specific parameter set. By examining the profile of these projection weights (variance explained by each mode) for parameter sets systematically chosen along a path crossing a phase boundary, we can characterize the transition. Fig. 4b illustrates this for points near the asynchrony-wave phase boundary. Parameter sets located on a phase boundary exhibit projection profiles that are hybrid, sharing features with the profiles of the pure phases they separate. SVD modes characteristic of both pure asynchrony and pure wave states contribute significantly. At a triple point, the dynamics reflect a richer mixture, with contributions from modes associated with all three converging phases (see Fig.5 ). This suggests that at these critical junctures in parameter space, the system’s dynamics are not committed to a single attractor but can explore or blend features of multiple underlying dynamical regimes. This “mode mixing” at boundaries and triple points can be intuitively understood from our spectral theory. Near these critical regions, the eigenvalue spectrum of the connectivity matrix may exhibit near-degeneracies, where multiple eigenvalues (corresponding to different potential dynamical patterns, e.g., a k = 0 oscillatory mode and a k ≠ 0 wave mode) have comparable real parts close to the instability threshold. Small amounts of perturbation can then cause the system to fluctuate between, or simultaneously express aspects of, these competing dynamical modes. This is consistent with experimental observations where brain activity can show transient or mixed features, especially during state transitions [ 46 , 47 ]. Encouragingly, the dominant SVD spatial modes extracted from our model’s phase velocity fields, such as plane waves or spiral patterns ( Fig. 4c ), show qualitative similarities to modes extracted from mesoscopic optical imaging of mouse cortex [ 4 ] across different arousal states ( Fig. 4d ). This suggests that the underlying principles governing the formation and transition of spatiotemporal patterns in our model may capture salient aspects of real brain dynamics. A.3 Order Parameters In order to numerically validate the correctness of the prediction of phase. We introduce order parameters CV . (Coefficient of Variation) and Mean Acti . (Mean Activity) to detect the magnitude of the neural activity fluctuation, introduce Osc. index (Oscillation Index) to detect the neural oscillation, introduce Local Sync . (Local Synchronization) and Moran’s Index to detect the spatial patterns of neural activity. Mean Acti (Mean Activity) This order parameter calculates the magnitude of neurons’ firing rate. This order parameter is useful for neural networks with hyperbolic tangent activation functions, given by CV (Coefficient of Variation) This order parameter describes the magnitude of neural activity fluctuation. It calculates the ratio between the standard deviation and mean of the neural firing rate. This order parameter is useful for neural networks with rectified linear, supra-linear and etc. activation functions: Osc. Index (Oscillation Index) This parameter describes the magnitude of neural oscillation. Similar to [ 48 ], we define this order parameter as the fraction of the Fourier spectrum energy of the peak concentrated at oscillation frequency. The Fourier spectrum is calculated by averaging the Fourier spectrum of all the excitatory neurons. Local Sync (Local Synchronization) This parameter describes the degree of local synchronization of neural activity. It calculates the ratio between the mean of the firing rate and the mean of the absolute value of the firing rate of local neurons. We define local neurons as the neurons within a square with a side length of 10 times the inter-neuron spacing. We chose this side length so that it’s similar to the detection range of local field potential. This order parameter is useful for neural networks with hyperbolic tangent activation functions: Moran’s Inde This parameter detects the existence of spatial patterns. It calculates the ratio between the mean of correlation between local neurons and the mean of correlation between all neurons. The definition of “local neurons” is the same to Local Sync: A.4 Numerical Simulation For the theoretical prediction of phase, we predict 61 × 61 gird points for a single phase diagram. For numerical simulation of phase prediction, the grid points is 21 × 21. For each grid point, we independently initialize the connectivity matrix and conduct 5 times simulations with a total simulation time of 100 times of time constant of membrane potential τ and a step length of 0.01τ. We only begin to calculate the order parameters after 25τ of simulation in order to avoid the influence of the transient process. Our numerical experiments were conducted on a computing cluster consisting of 16 nodes, each equipped with an Intel(R) Xeon(R) CPU E5-2407 0 @ 2.20GHz. A.5 Numerical Experiments with Alternative Activation Functions A.6 Numerical Results of Order Parameters and Phase Diagrams A.7 Parameters of Phase Diagrams View this table: View inline View popup Download powerpoint A.8 Relation between Eigenvalues and Dynamics We want to understand the relationship between connectivity structure and dynamics of spatial distributed neural networks. To begin with, we first consider a simple linear neural network, of which dynamics follows: Because it’s a linear dynamical system, the dynamics can be decomposed into different independent modes. Let’s assume the connectivity matrix J is diagonalizable. J = A Λ A −1 , where A is composed of eigenvectors of connectivity matrix, A = [ ν 1 , ν 2 ,…, ν N ], and Λ = [λ 1 , λ 2 ,…, λ N ] is composed of eigenvalues of connectivity matrix. The dynamical equation can be rewritten as, Therefore, we can consider the dynamics of each eigenvector component to be independent. The activity of neurons is a superposition of components in different directions. Let h(t) = ∑ i c i (t)ν i . c i is the magnitude of independent components, which satisfies, For a component with Re(λ i ) < 1, the magnitude is bounded and fluctuates around 0. For a component with Re(λ i ) ≥ 1, the magnitude increases over time as a rate of e Re ( λ i ) −1 . Therefore, the component with an eigenvalue with the largest real part dominates the dynamics of neural networks. For a non-linear neural network, we can also use eigenvalues to understand the relationship between connectivity structure and eigenvalues. We only need to perform a linear expansion around the fixed point of neural activity. Let’s denote the fixed point of membrane potential as h ∗ , and the deviation from the fixed point as δh. The neural activity follows, We can perform a linear expansion around the fixed point h ∗ . The deviation δh from fixed point follows, Therefore, we can consider a nonlinear neural network equivalent to a linear neural network with an effective connectivity matrix . The dynamical regime of neural networks can be characterized through spectral analysis of the effective connectivity matrix . When all eigenvalues satisfy Re(λ i ) < 1, neural activity remains bounded near the fixed point, exhibiting small-amplitude fluctuations around the steady-state firing rate. This regime corresponds to asynchronous irregular activity due to the absence of dominant eigenmodes. When spectral outliers emerge with Re(λ i ) ≥ 1, the corresponding eigenmodes dominate network dynamics. These regimes can be systematically classified ( Table 1 ) based on two spectral properties of the dominant eigenvalues: 1) temporal frequency (real vs. complex eigenvalues) and 2) spatial frequency (wave vector k of eigenvectors). Spatial organization in eigenvectors generates distinct spatiotemporal patterns, enabling classification into four phases: synchronized state, oscillatory phase, bump attractor, and traveling wave. In contrast, when dominant eigenvalues reside within the bulk spectral disk, the system enters a chaotic phase characterized by: (i) spatially unstructured eigenvectors, (ii) large-amplitude fluctuations, and (iii) weak inter-neuronal correlations. Intuitively, the absence of spatial patterning in neural activity stems from the structural homogeneity of dominant spectral bulk eigenmodes. Weak inter-neuronal correlations emerge from high-dimensional superposition of components satisfying Re(λ i ≥ 1). While spectral analysis provides initial insights, analytical determination of phase boundaries requires dynamical mean-field theory, as detailed in Appendix A.11. This method is not fully mathematically rigorous. First, it cannot fully deal with the case of multiple fixed points. Second, it cannot fully predict the dynamic behavior when the network activity goes far from a fixed point. Besides, it cannot fully characterize the dynamic behavior where the fixed point is heterogeneous between different neurons. Therefore, to verify the validity of our theory under nonlinear conditions, we performed many numerical experiments as aforementioned to demonstrate that our theory is indeed correct under nonlinear conditions and is useful and informative in understanding the relationship between neural networks’ dynamics and their connectivity structure. A.9 Neural Activity of the Asynchronous State As shown in Fig. 8 , in the asynchronous state, neurons are weakly correlated, therefore, there are no spatial patterns. However, [ 12 ] reveals low-rank structure of the connectivity matrix can produce outliers in the eigenspectrum of the long-time window covariance matrix. It’s valuable to study the eigenspectrum of the spatially extended neural networks’ covariance matrix in the future. Download figure Open in new tab Figure 8. Neural activity of asynchronous state. (a) temporal evolution of membrane potentials of excitatory neurons. (b) spatial patterns of excitatory neurons’ firing rate. Download figure Open in new tab Figure 9. The eigenspectrum of a random matrix with elements from i.i.d Gaussian distribution with unit variance and zero mean. Besides, the neural activity deviates from the fixed point due to the external input. In our experiments, the external input is independent white noise with small magnitude, thus the magnitude of deviation of neural activity is small. However, with other large input, the neural activity can have different behavior. As illustrated in Fig. 8 , neural populations in the asynchronous state exhibit weak pairwise correlations, resulting in the absence of emergent spatial patterns. Recent theoretical advances by [ 12 ], however, demonstrate that low-rank connectivity structures can induce spectral outliers in the longtime window covariance matrix. This finding motivates systematic investigation of spectral properties of covariance matrices in spatially extended neural networks a promising direction for future research. Notably, weak deviations from the fixed point emerge through external input modulation. In our experimental paradigm, network input originates from independent white noise sources with low amplitude, resulting in small deviations from the fixed point. The neural activity under other types of external input needs to be further explored. A.10 Eigenvalues and Eigenvectors of Spatially Distributed Neural Networks A.10.1 Circular Law A classical result of the random matrix is the circular law. If you take an n × n matrix with independent and identically distributed (i.i.d.) entries , then as n grows large, the eigenvalues of the matrix become uniformly distributed inside the disk with radius σ in the complex plane [ 49 ]. This result can be further generalized to more general random matrices where the i.i.d distribution is not Gaussian but has a finite variance. Theorem 1 Let A n be the n × n random matrix whose entries are i . i . d. complex random variables with mean 0 and variance 1. The empirical spectral distribution of then converges (both in probability and in the almost sure sense) to the uniform distribution on the unit disk [ 21 ] . This theorem means that the circular law can not only be used in ideal Gaussian distribution but also can be used in the case of other distributions like sparse connection and biologically plausible log-normal distribution, etc. However, it still requires the distribution to be identical independent distribution, while there are multiple types of neurons and connection probabilities that decay with the distance between neurons. Therefore, we still can not use this theorem in biologically plausible spatially distributed neural networks. Thanks to techniques of free probability theory [ 23 ], we can deal with the case that the distribution is independent, zero-mean but not identical. Besides, using dynamical mean-field theory can also derive similar results [ 24 ]. Theorem 2 Let A n = (σ ij ) be an n × n deterministic matrix, and X n = (X ij ) be an n × n random matrix with i . i . d. centered entries of unit variance. Define the rescaled matrix: where ° denotes the Hadamard product. Let denote the empirical spectral distribution (ESD) of Y n . For variance profiles has a positive density on the centered disc of radius , where and ρ(V n ) is its spectral radius. [23] Using this theorem, we can deal with biologically plausible neural networks with multiple types of neurons and spatial distribution. Although the connection strength distribution between different types of neurons are different, and the distribution also varies with distance between neurons, we only need to calculate a profile matrix σ ij above to determine the radius of a eigenspectrum disk. A.10.2 Outliers of Eigenspectrum The theorems mentioned above are all required of a zero-mean condition. However, in biologically plausible case, due to Dale’s Law, the distribution of connection strength between certain types of neurons cannot have a zero-mean. For example, the connection strength between excitatory neurons must have a positive mean. [ 41 ] considers a special case of non-zero-mean distributions. In their model, neurons are sparsely connected with a certain probability p for both excitatory and inhibitory neurons. The connection strength are the same among excitatory synapses and inhibitory synapses separately. They showed that the eigenspectrum of the connectivity matrix is composed of a bulk disk part and an outlier. This outlier lies in the position of the mean connection strength of all synapses. Namely, the non-zero mean of a distribution creates outliers. The example above is actually a special case of low-rank perturbation on a zero-mean random matrix. We can gain the following intuition. The connectivity matrix can be decomposed into two parts: a determined mean part and a zero-mean random matrix part. The zero-mean random matrix part creates a bulk disk eigenspectrum. And the determined mean part is actually a low-rank matrix because its rank is less than the number of types of neurons. This low-rank perturbation creates outliers in eigenspectrums [ 22 ] provided a rigorous mathematical theorem on outliers and low-rank perturbation. Theorem 3 Let X n be an iid random matrix with finite fourth moment, and for each n, let C n be a deterministic matrix with rank O(1) and operator norm O(1). Assume that for large n, C n has no eigenvalues in {z ∈ ℂ : 1 + ε < |z| < 1 + 3ε} and has j = O(1) eigenvalues in {z ∈ ℂ : |z| ≥ 1 + 3ε}. Then, almost surely, for large n , has exactly j eigenvalues in {z ∈ ℂ : |z| ≥ 1 + 2ε}, and these eigenvalues satisfy as n → ∞ for each 1 ≤ i ≤ j [22] . With this theorem, we can finally characterize the eigenspectrum of biologically plausible neural networks. The sparsity and variance of connection strength both contributed to the zero-mean random part, which creates a bulk disk part of eigenspectrum. And the Dale’s law force the connectivity matrix have a determined part. This part is often low-rank, thus creating outliers of the eigensepctrum. A.10.3 Eigenspectrum of connectivity matrices of spatially extended neural networks As mentioned above, the eigenvalues of the connectivity matrix J consist of a circular bulk part and a set of outliers. [ 22 ] noted that for a random matrix subjected to a low-rank perturbation, the eigenvalues of the new matrix largely remain within the original circle, with any outliers located where the eigenvalues of the perturbation matrix lie. We can conceptualize the connectivity matrix J as comprising the expected values of each element, , and the deviations from this expectation, δ J . The matrix δ J corresponds to the aforementioned random matrix, while corresponds to the low-rank perturbation matrix. Therefore, δ J dictates the bulk part of the eigenvalue spectrum, while determines the outlier part. A.10.4 Outliers Eigenvectors and the Spatial Translation Invariance of Neural Networks The matrix determines the outlier part of the eigenvalue spectrum [ 22 ]. In order to calculate the eigenvalues of the matrix , we need to utilize its property of spatial translation invariance Let’s start with a relatively simple one-dimension case. The expected component of the connectivity matrix, , can be expressed as a block matrix: where represents the expected connectivity from neurons of type β to neurons of type α, forming a matrix of size N α × N β . It satisfies the following relation, The element depends only on the distance between neurons i and j and their respective types. In our model, the proportion of excitatory to inhibitory neurons satisfies N E : N I = 4 : 1. We can consider four excitatory neurons and one inhibitory neuron as forming a small unit, and the neural network consists of repeated instances of this unit. When the neural network is collectively translated by several unit distances in physical space, the matrix remains unchanged, indicating its translational invariance. We define the following block matrix P , which satisfies: where P α is a matrix of size N α × N α : The matrix P represents the translation of the neural network by one small unit and has the property: This implies that the matrices P and share common eigenvectors. The eigenvalues of P satisfy , where . Each eigenvalue λ n is fivefold degenerate, with the eigenspace spanned by the vectors , which satisfy the following properties. These vectors can be expressed as , where is a vector of length N E and is a vector of length N I . For l = 0, 1, …, 3: For l = 4: Thus, in the basis: , the matrix can be written as a block diagonal matrix, where each block corresponds to an effective connectivity matrix for a specific Fourier mode. The effective connectivity matrix for a Fourier mode with wave vector k is given by: This expression is accurate for long-wavelength modes. the matrix takes the form: where For the k = 0 Fourier mode, the outliers are precisely the two eigenvalues given by , each of which is non-degenerate. For k ≠ 0 modes, since the eigenvalues corresponding to ± k coincide, they exhibit twofold degeneracy. For the two dimensions case, the mathematical derivation is almost the same. In the case of two dimensions, the spatial translation can be done in both directions. Therefore, the eigenvectors are plane waves. Effective Connectivity Matrix and Eigenvalues For each spatial Fourier mode, we can further simplify the effective connectivity matrix as follows: where , n x , n y = 0, ±1, ±2, … As k takes on different values, the effective connectivity matrices for different spatial Fourier modes yield different eigenvalues. These eigenvalues constitute the outliers in the eigenvalue spectrum of the connectivity matrix. In the case without spatial distribution, the number of outliers corresponds to the number of neuron types; in the presence of spatial distribution, as the spatial scale of the network increases, the number of outliers also increases and is ordered according to their corresponding wave vectors k. When the spatial scale approaches infinity, these outliers will be arranged along a continuous curve (here, we reference the figure from spatial effect). The two-dimensional case is similar, except that the wave vector k is a vector that takes values , n x , n y = 0, ± 1, ± 2, The eigenvalues in the two-dimensional case are also given by the effective connectivity matrix eigenvalues for different k values. The eigenvectors corresponding to the outlier eigenvalues are provided by the spatial Fourier modes indicated by k. Bulk disk part The matrix δ J governs the circular part of the eigenvalue spectrum. [ 23 ] and [ 24 ] provide formulations for the eigenvalue distribution of random matrices with independent entries, mean zero, and finite variance. The eigenvalues are distributed within a circle of a specific radius , where the elements of the matrix M represent the variances of the elements of the connectivity matrix J . The variance of the elements of the connectivity matrix J is given by: Similar to the mean part, the matrix is also a spatial translation invariant matrix, because both the variance and connection probability of the synapses between two neurons are only related to the relative distance between these two neurons. Therefore, the variance matrix M can also be diagonalized into a series of small block matrices. Let’s start with the one dimension case, where is the submatrix under a wave vector k . And its elements satisfy . Because the radius of the bulk disk of an eigenspectrum only depends on the largest eigenvalue of , and all the elements of are positive, we only need to consider the zero wave vector submatrix . Zero wave vectors are spatially uniform, thus the elements of submatrix are the average of original elements of . For simplicity, we denote the matrix as M αβ . For the one-dimensional case, this reduced matrix M αβ is: For the two-dimensional case, the reduced matrix M αβ is: A.11 Phase Boundary of Chaos Phase As mentioned above, the theoretical prediction based on linearization and eigenspectrum cannot fully characterize the dynamical behavior far away from a fixed point. The chaos phase is one of these cases. The bulk disk of an eigenspectrum corresponds to neural activity without spatial patterns, while outliers are related to spatially ordered neural activity. If both the radius of the bulk disk and the real part of outliers are larger than 1, they will compete against each other and we cannot directly determine whether the neural activity is spatially ordered. However, with the tool of dynamical mean-field theory (DMFT) [ 30 ], we can mathematically rigorously derive the condition for spatial patterns in the chaos phase. For a neural network with the radius of a bulk disk part greater than 1, if the deviation of the mean of the neurons’ membrane potential gradually amplifies after a perturbation, eventually the neural activities of different neurons will become synchronized, and the neural network will no longer be in the chaos phase. [ 50 ] calculated how the mean and variance of the neuronal membrane potential evolve after being subjected to perturbation, where, Here, taking the average ⟨·⟩ refers to averaging over the distribution of the neuronal membrane potentials. The distribution of the neuronal membrane potentials follows a Gaussian distribution, the mean and variance of which can be theoretically calculated based on dynamic mean-field theory [ 50 ]. Let’s denote the response matrix as follows, For a neural network with the radius of bulk disk part larger than 1, if all the eigenvalues of a response matrix are less than 0, the chaotic neural activity is stable and the neural activity is within a chaos phase. I further extended this result to spatially extended neural networks. We can regard a neural network with spatial distribution to be composed of lots of populations with area Δ S . Therefore, a neural network with spatial distribution is equivalent to a neural network with infinite populations. In our model, In the model presented in the main text, the activation function for neurons is the hyperbolic tangent function, tanh. The distribution of membrane potential is a Gaussian distribution with zero mean. Therefore, the matrix B and E are zero matrix. The response matrix R is a block diagonal matrix. The submatrix in the upper left corner represents the mean of membrane potential response to external inputs, and the submatrix I − D in the lower right corner represents the variance of membrane potential response to external inputs. The stability of the mean of membrane potential determines the phase boundary between the chaos phase and other phases with locally synchronized neural activity. Therefore, the largest real part of eigenvalues of the matrix determines the phase boundary of the chaos phase. Similar to the calculation above, the matrix is spatially translation invariant. Therefore, we can diagonalize this matrix into a series of submatrices corresponding to different wave vectors k. The elements of submatrices are as follows, where , n x , n y = 0, ± 1, ± 2, Averaging is performed over the distribution of neuronal membrane potentials. Using the tool of dynamical mean-field theory [ 50 ], the distribution of neuronal membrane potentials follows a Gaussian distribution. Let’s denote the expectation of membrane potential as u α , the expectation of firing rate as m α , the variance of membrane potential as Δ α , and the autocorrelation of firing rate as C α . They satisfy, where h β is the direct current (DC) input to the β type of neurons. z follows a standard Gaussian distribution with unit variance. where h β is the direct current (DC) input to the β type of neurons. z follows a standard Gaussian distribution with unit variance. If the eigenvalues of all the submatrices are less than 1, the chaotic neural activity is stable and the neural network is in the chaos phase. Otherwise, the chaotic neural activity is unstable. We can assign its phase based on the wave vector k and whether the eigenvalues of instability are complex to determine which phase it belongs to. References [1]. ↵ Gyorgy Buzsaki and Andreas Draguhn . Neuronal oscillations in cortical networks . science , 304 ( 5679 ): 1926 – 1929 , 2004 . OpenUrl Abstract / FREE Full Text [2]. Rory G Townsend and Pulin Gong . Detection and analysis of spatiotemporal patterns in brain activity . PLoS computational biology , 14 ( 12 ): e1006643 , 2018 . OpenUrl CrossRef [3]. Michael Schirner , Xiaolu Kong , BT Thomas Yeo , Gustavo Deco , and Petra Ritter . Dynamic primitives of brain network interaction . NeuroImage , 250 : 118928 , 2022 . OpenUrl CrossRef PubMed [4]. ↵ Yuqi Liang , Junhao Liang , Chenchen Song , Mianxin Liu , Thomas Knöpfel , Pulin Gong , and Changsong Zhou . Complexity of cortical wave patterns of the wake mouse cortex . Nature Communications , 14 ( 1 ): 1434 , 2023 . OpenUrl CrossRef PubMed [5]. ↵ Lyle Muller , Frédéric Chavane , John Reynolds , and Terrence J Sejnowski . Cortical travelling waves: mechanisms and computational principles . Nature Reviews Neuroscience , 19 ( 5 ): 255 – 268 , 2018 . OpenUrl CrossRef PubMed [6]. ↵ Honghui Zhang , Andrew J Watrous , Ansh Patel , and Joshua Jacobs . Theta and alpha oscillations are traveling waves in the human neocortex . Neuron , 98 ( 6 ): 1269 – 1281 , 2018 . OpenUrl CrossRef PubMed [7]. ↵ Yifan Gu , Yang Qi , and Pulin Gong . Rich-club connectivity, diverse population coupling, and dynamical activity patterns emerging from local cortical circuits . PLoS computational biology , 15 ( 4 ): e1006902 , 2019 . OpenUrl CrossRef PubMed [8]. ↵ Guozhang Chen and Pulin Gong . Computing by modulating spontaneous cortical activity patterns as a mechanism of active visual processing . Nature communications , 10 ( 1 ): 4915 , 2019 . OpenUrl CrossRef PubMed [9]. ↵ Guozhang Chen and Pulin Gong . A spatiotemporal mechanism of visual attention: Superdiffusive motion and theta oscillations of neural population activity patterns . Science advances , 8 ( 16 ): eabl4995 , 2022 . OpenUrl CrossRef PubMed [10]. ↵ Zhuda Yang , Junhao Liang , and Changsong Zhou . Critical avalanches in excitation-inhibition balanced networks reconcile response reliability with sensitivity for optimal neural representation . Physical Review Letters , 134 ( 2 ): 028401 , 2025 . OpenUrl CrossRef PubMed [11]. ↵ Stephen Coombes , Peter beim Graben , Roland Potthast , and James Wright . Neural fields: theory and applications . Springer , 2014 . [12]. ↵ Yu Hu and Haim Sompolinsky . The spectrum of covariance matrices of randomly connected recurrent neuronal networks with linear dynamics . PLoS computational biology , 18 ( 7 ): e1010327 , 2022 . OpenUrl CrossRef PubMed [13]. ↵ Srdjan Ostojic . Two types of asynchronous activity in networks of excitatory and inhibitory spiking neurons . Nature neuroscience , 17 ( 4 ): 594 – 600 , 2014 . OpenUrl CrossRef PubMed [14]. ↵ Robert Rosenbaum and Brent Doiron . Balanced networks of spiking neurons with spatially dependent recurrent connections . Physical Review X , 4 ( 2 ): 021039 , 2014 . OpenUrl [15]. Yang Qi and Pulin Gong . Dynamic patterns in a two-dimensional neural field with refractoriness . Physical Review E , 92 ( 2 ): 022702 , 2015 . OpenUrl [16]. ↵ Ryan Pyle and Robert Rosenbaum . Spatiotemporal dynamics and reliable computations in recurrent spiking neural networks . Physical review letters , 118 ( 1 ): 018103 , 2017 . OpenUrl CrossRef PubMed [17]. Francesca Mastrogiuseppe and Srdjan Ostojic . Linking connectivity, dynamics, and computations in low-rank recurrent neural networks . Neuron , 99 ( 3 ): 609 – 623 , 2018 . OpenUrl CrossRef PubMed [18]. ↵ Chengcheng Huang , Douglas A Ruff , Ryan Pyle , Robert Rosenbaum , Marlene R Cohen , and Brent Doiron . Circuit models of low-dimensional shared variability in cortical networks . Neuron , 101 ( 2 ): 337 – 348 , 2019 . OpenUrl CrossRef PubMed [19]. ↵ Noga Mosheiff , Bard Ermentrout , and Chengcheng Huang . Chaotic dynamics in spatially distributed neuronal networks generate population-wide shared variability . PLOS Computational Biology , 19 ( 1 ): e1010843 , 2023 . OpenUrl CrossRef PubMed [20]. ↵ Srdjan Ostojic and Stefano Fusi . Computational role of structure in neural activity and connectivity . Trends in Cognitive Sciences , 2024 . [21]. ↵ Terence Tao , Van Vu , and Manjunath Krishnapur . Random matrices: Universality of esds and the circular law . 2010 . [22]. ↵ Terence Tao . Outliers in the spectrum of iid matrices with bounded rank perturbations . Probability Theory and Related Fields , 155 ( 1 ): 231 – 263 , 2013 . OpenUrl CrossRef [23]. ↵ Nicholas Cook , Walid Hachem , Jamal Najim , and David Renfrew . Non-hermitian random matrices with a variance profile (ii): properties and examples . Journal of Theoretical Probability , 35 ( 4 ): 2343 – 2382 , 2022 . OpenUrl CrossRef [24]. ↵ Johnatan Aljadeff , David Renfrew , Marina Vegué , and Tatyana O Sharpee . Low-dimensional dynamics of structured random networks . Physical Review E , 93 ( 2 ): 022302 , 2016 . OpenUrl CrossRef [25]. ↵ Kanaka Rajan and Larry F Abbott . Eigenvalue spectra of random matrices for neural networks . Physical review letters , 97 ( 18 ): 188104 , 2006 . OpenUrl CrossRef PubMed [26]. ↵ Helen E Scharfman . The neurobiology of epilepsy . Current neurology and neuroscience reports , 7 ( 4 ): 348 – 354 , 2007 . OpenUrl CrossRef PubMed Web of Science [27]. ↵ György Buzsáki and Xiao-Jing Wang . Mechanisms of gamma oscillations . Annual review of neuroscience , 35 ( 1 ): 203 – 225 , 2012 . OpenUrl CrossRef PubMed Web of Science [28]. ↵ Haleigh N Mulholland , Matthias Kaschube , and Gordon B Smith . Self-organization of modular activity in immature cortical networks . Nature communications , 15 ( 1 ): 4145 , 2024 . OpenUrl CrossRef PubMed [29]. ↵ Gordon B Smith , Bettina Hein , David E Whitney , David Fitzpatrick , and Matthias Kaschube . Distributed network interactions and their emergence in developing neocortex . Nature neuroscience , 21 ( 11 ): 1600 – 1608 , 2018 . OpenUrl CrossRef PubMed [30]. ↵ Haim Sompolinsky , Andrea Crisanti , and Hans-Jurgen Sommers . Chaos in random neural networks . Physical review letters , 61 ( 3 ): 259 , 1988 . OpenUrl CrossRef PubMed Web of Science [31]. ↵ Matthew D Greaves , Leonardo Novelli , Sina Mansour L , Andrew Zalesky , and Adeel Razi . Structurally informed models of directed brain connectivity . Nature Reviews Neuroscience , pages 1 – 19 , 2024 . [32]. ↵ Guozhang Chen , Franz Scherr , and Wolfgang Maass . A data-based large-scale model for primary visual cortex enables brain-like robust and versatile visual processing . Science advances , 8 ( 44 ): eabq7592 , 2022 . OpenUrl CrossRef PubMed [33]. ↵ Bard Ermentrout . Neural networks as spatio-temporal pattern-forming systems . Reports on progress in physics , 61 ( 4 ): 353 , 1998 . OpenUrl CrossRef [34]. ↵ James M Shine . Neuromodulatory control of complex adaptive dynamics in the brain . Interface focus , 13 ( 3 ): 20220079 , 2023 . OpenUrl CrossRef PubMed [35]. ↵ Michael Breakspear . Dynamic models of large-scale brain activity . Nature neuroscience , 20 ( 3 ): 340 – 352 , 2017 . OpenUrl CrossRef PubMed [36]. ↵ Viktor K Jirsa , William C Stacey , Pascale P Quilichini , Anton I Ivanov , and Christophe Bernard . On the nature of seizure dynamics . Brain , 137 ( 8 ): 2210 – 2230 , 2014 . OpenUrl CrossRef PubMed Web of Science [37]. ↵ Efstathios D Kondylis , Michael J Randazzo , Ahmad Alhourani , Witold J Lipski , Thomas A Wozny , Yash Pandya , Avniel S Ghuman , Robert S Turner , Donald J Crammond , and R Mark Richardson . Movement-related dynamics of cortical oscillations in parkinsons disease and essential tremor . Brain , 139 ( 8 ): 2211 – 2223 , 2016 . OpenUrl CrossRef PubMed [38]. ↵ Dietmar Plenz , Tiago L Ribeiro , Stephanie R Miller , Patrick A Kells , Ali Vakili , and Elliott L Capek . Self-organized criticality in the brain . Frontiers in Physics , 9 : 639389 , 2021 . OpenUrl CrossRef [39]. Karl J Friston , Lee Harrison , and Will Penny . Dynamic causal modelling . Neuroimage , 19 ( 4 ): 1273 – 1302 , 2003 . OpenUrl CrossRef PubMed Web of Science [40]. ↵ David Sussillo and Larry F. Abbott . Generating coherent patterns of activity from chaotic neural networks . Neuron , 63 ( 4 ): 544 – 557 , 2009 . OpenUrl CrossRef PubMed Web of Science [41]. ↵ Isabelle D Harris , Hamish Meffin , Anthony N Burkitt , and Andre DH Peterson . Effect of sparsity on network stability in random neural networks obeying dale’s law . Physical Review Research , 5 ( 4 ): 043132 , 2023 . OpenUrl CrossRef [42]. ↵ Agostina Palmigiano , Francesco Fumarola , Daniel P Mossing , Nataliya Kraynyukova , Hillel Adesnik , and Kenneth D Miller . Common rules underlying optogenetic and behavioral modulation of responses in multi-cell-type v1 circuits . bioRxiv , pages 2020 – 11 , 2020 . [43]. ↵ Jesper R Ipsen and Andre DH Peterson . Consequences of dale’s law on the stability-complexity relationship of random neural networks . Physical Review E , 101 ( 5 ): 052412 , 2020 . OpenUrl [44]. ↵ David Dahmen , Moritz Layer , Lukas Deutz , Paulina Anna Dbrowska , Nicole Voges , Michael von Papen , Thomas Brochier , Alexa Riehle , Markus Diesmann , Sonja Grün , and Moritz Helias . Global organization of neuronal activity only requires unstructured local connectivity . eLife , 11 : e68422 , jan 2022 . ISSN 2050–084X . doi: 10.7554/eLife.68422 . URL https://doi.org/10.7554/eLife.68422 . OpenUrl CrossRef PubMed [45]. ↵ Jun Igarashi , Hatsuo Hayashi , and Katsumi Tateno . Theta phase coding in a network model of the entorhinal cortex layer II with entorhinal-hippocampal loop connections . Cognitive Neurodynamics , 1 ( 2 ): 169 – 184 , June 2007 . ISSN 1871-4099 . doi: 10.1007/s11571-006-9003-8 . URL https://doi.org/10.1007/s11571-006-9003-8 . OpenUrl CrossRef PubMed [46]. ↵ Gerhard Werner . Consciousness viewed in the framework of brain phase space dynamics, criticality, and the renormalization group . 55 : 3 – 12 . ISSN 0960-0779 . doi: 10.1016/j.chaos.2012.03.014 . OpenUrl CrossRef [47]. ↵ James A. Roberts , Leonardo L. Gollo , Romesh G. Abeysuriya , Gloria Roberts , Philip B. Mitchell , Mark W. Woolrich , and Michael Breakspear . Metastable brain waves . 10 ( 1 ): 1056 . ISSN 2041-1723 . doi: 10.1038/s41467-019-08999-0 . OpenUrl CrossRef PubMed [48]. ↵ Leonardo Dalla Porta and Mauro Copelli . Modeling neuronal avalanches and long-range temporal correlations at the emergence of collective oscillations: Continuously varying exponents mimic m/eeg results . PLoS computational biology , 15 ( 4 ): e1006924 , 2019 . OpenUrl CrossRef PubMed [49]. ↵ ML Mehta and Random Matrices . the statistical theory of energy levels , 1967 . [50]. ↵ Jonathan Kadmon and Haim Sompolinsky . Transition to chaos in random neuronal networks . Physical Review X , 5 ( 4 ): 041030 , 2015 . OpenUrl View the discussion thread. Back to top Previous Next Posted May 27, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Bridging Scales: Spectral Theory Reveals How Local Connectivity Rules Sculpt Global Neural Dynamics in Spatially Extended Networks Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Bridging Scales: Spectral Theory Reveals How Local Connectivity Rules Sculpt Global Neural Dynamics in Spatially Extended Networks Yuhan Huang , Keren Gao , Dongping Yang , Sen Song , Guozhang Chen bioRxiv 2025.05.23.655871; doi: https://doi.org/10.1101/2025.05.23.655871 Share This Article: Copy Citation Tools Bridging Scales: Spectral Theory Reveals How Local Connectivity Rules Sculpt Global Neural Dynamics in Spatially Extended Networks Yuhan Huang , Keren Gao , Dongping Yang , Sen Song , Guozhang Chen bioRxiv 2025.05.23.655871; doi: https://doi.org/10.1101/2025.05.23.655871 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Biophysics Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17697) Bioengineering (13895) Bioinformatics (41951) Biophysics (21456) Cancer Biology (18594) Cell Biology (25520) Clinical Trials (138) Developmental Biology (13381) Ecology (19903) Epidemiology (2067) Evolutionary Biology (24323) Genetics (15612) Genomics (22510) Immunology (17738) Microbiology (40401) Molecular Biology (17184) Neuroscience (88622) Paleontology (667) Pathology (2833) Pharmacology and Toxicology (4825) Physiology (7644) Plant Biology (15158) Scientific Communication and Education (2046) Synthetic Biology (4296) Systems Biology (9825) Zoology (2271)