The geometry and dimensionality of brain-wide activity

We recorded brain-wide Ca²⁺activity at a volume rate of 10 Hz in head-fixed larval zebrafish (Figure 2A) during hunting attempts (Methods) and spontaneous behavior using a Fourier light-field microscopy (Cong et al., 2017). Approximately 2000 ROIs (1977.3 ± 677.1, mean ± SD) with a diameter of 16.84 ± 8.51 µm were analyzed per fish based on voxel activity (Methods, Figure 2—figure supplement 1). These ROIs likely correspond to multiple nearby neurons with correlated activity. Henceforth, we refer to the ROIs as ‘neurons’ for simplicity.

We first investigate the dimensionality of neural activity ${\displaystyle D_{\text{PR}}}$ (Figure 1B) in a randomly chosen cell assembly in zebrafish, similar to multi-area Neuropixels recording in a mammalian brain. We focus on how ${\displaystyle D_{\text{PR}}}$ changes with a large sample size ${\displaystyle N}$. We find that if the mean squared covariance remains finite instead of vanishing with ${\displaystyle N}$, the dimensionality ${\displaystyle D_{\text{PR}}}$ (Figure 1B) becomes sample size independent and depends only on the variance ${\displaystyle {\sigma }_i^2}$ and the covariance ${\displaystyle C_{ij}}$ between neurons ${\displaystyle i}$ and ${\displaystyle j}$:

where ${\displaystyle \mathrm{E}(\dots )}$ denotes average across neurons (Methods and Dahmen et al., 2020). The finite mean squared covariance condition is supported by the observation that the neural activity covariance ${\displaystyle C_{ij}}$ is positively biased and widely distributed with a long tail (Figure 2—figure supplement 2A). As predicted, the data dimensionality grows with sample size and reaches the maximum value specified by Equation 1 (Figure 2D).

Next, we investigate the shape of the neural activity space described by the eigenspectrum of the covariance matrix derived from the activity of ${\displaystyle N}$ randomly selected neurons (Figure 2C). When the eigenvalues are arranged in descending order and plotted against the normalized rank ${\displaystyle r/N}$, where ${\displaystyle r=1,\dots ,N}$ (we refer to it as the rank plot), this curve shows an approximate power law that spans 10 folds. Interestingly, as the size of the covariance matrices decreases (${\displaystyle N}$ decreases), the eigenspectrum curves nearly collapse over a wide range of eigenvalues. This pattern holds across diverse datasets and experimental techniques (Figure 2F, Figure 2—figure supplement 2E–L). The similarity of the covariance matrices of randomly sampled neural populations can be intuitively visualized (Figure 2E), after properly sorting the neurons (Methods).

The scale invariance in the neural covariance matrix – the collapse of the covariance eigenspectrum under random sampling – is non-trivial. The spectrum is not scale invariant in a common covariance matrix model based on independent noise (Figure 2G). It is absent when replacing the neural covariance matrix eigenvectors with random ones, keeping the eigenvalues identical (Figure 2H). A recurrent neural network with random connectivity (Hu and Sompolinsky, 2022) does not yield a scale-invariant covariance spectrum (Figure 2I). A recently developed latent variable model (Morrell et al., 2024; Appendix 1—figure 6), which is able to reproduce avalanche criticality, also fails to generate the scale-invariant covariance spectrum. Thus, a new model is needed for the covariance matrix of neural activity.

Dimension reduction methods simplify and visualize complex neuron interactions by embedding them into a low-dimensional map, within which nearby neurons have similar activities. Inspired by these ideas, we use the ERM (Mézard et al., 1999) to model neural covariance. Imagine sprinkling neurons uniformly distributed on a ${\displaystyle d}$-dimensional functional space of size ${\displaystyle L}$ (Figure 3A), where the distance between neurons ${\displaystyle i}$ and ${\displaystyle j}$ affects their correlation. Let ${\displaystyle {\overrightarrow{x}}_i}$ represent the functional coordinate of the neuron ${\displaystyle i}$. The distance-correlation dependency is described by kernel function${\displaystyle f({\overrightarrow{x}}_i-{\overrightarrow{x}}_j)>0}$ with ${\displaystyle f(0)=1}$, indicating closer neurons have stronger correlations, and decreases as distance ${\displaystyle \Vert {\overrightarrow{x}}_i-{\overrightarrow{x}}_j\Vert }$ increases (Figure 3A and Methods). To model the covariance, we extend the ERM by incorporating heterogeneity of neuron activity levels (shown as the size of the neuron in the functional space in Figure 3A).

The variance of neural activity ${\displaystyle {\sigma }_i^2}$ is drawn i.i.d. from a given distribution and is independent of neurons’ position.

Figure 3 with 2 supplements see all

Download asset Open asset

This multidimensional functional space may represent attributes to which neurons are tuned, such as sensory features (e.g., visual orientation Hubel and Wiesel, 1959, auditory frequency) and movement characteristics (e.g., direction, speed Stefanini et al., 2020; Kropff et al., 2015). In sensory systems, it represents stimuli as neural activity patterns, with proximity indicating similarity in features. For motor control, it encodes movement parameters and trajectories. In the hippocampus, it represents the place field of a place cell, acting as a cognitive map of physical space (O’Keefe, 1976; Moser et al., 2008; Tingley and Buzsáki, 2018).

We first explore the ERM with various forms of ${\displaystyle f(\overrightarrow{x})}$ and find that fast-decaying functions like Gaussian and exponential functions do not produce eigenspectra similar to the data and no scale invariance over random sampling (Figure 3—figure supplement 1A–H and Appendix 2). Thus, we turn to slow-decaying functions including the power law, which produce spectra similar to the data (Figure 3C, D; see also Figure 3—figure supplement 2). We adopt a particular kernel function because of its closed-form and analytical properties: ${\displaystyle f(\overrightarrow{x})={ϵ}^{\mu }({ϵ}^2+\Vert \overrightarrow{x}{\Vert }^2{)}^{-\mu /2}}$ (Methods). For large distance ${\displaystyle \Vert \overrightarrow{x}\Vert \gg ϵ}$, it approximates a power law ${\displaystyle f(\overrightarrow{x})\approx {ϵ}^{\mu }\Vert \overrightarrow{x}{\Vert }^{-\mu }}$ and smoothly transitions at small distance to satisfy the correlation requirement ${\displaystyle f(0)=1}$ (Appendix 1—figure 3I, J).

To determine the conditions for scale invariance in ERM, we analytically calculate the eigenspectrum of covariance matrix ${\displaystyle C}$ (Equation 2) for large ${\displaystyle N,L}$ using the replica method (Mézard et al., 1999). A key order parameter emerging from this calculation is the neuron density ${\displaystyle \rho :=N/L^d}$. In the high-density regime ${\displaystyle \rho {ϵ}^d\approx 1}$, the covariance spectrum can be approximated in a closed form (Methods). For the slow-decaying kernel function ${\displaystyle f(\overrightarrow{x})}$ defined above, the spectrum for large eigenvalues follows a power law (Appendix 2):

where r is the rank of the eigenvalues in descending order and ${\displaystyle p(\lambda )}$ is their probability density function. Equation 3 intuitively explains the scale invariance over random sampling. Sampling in the ERM reduces the neuron density ρ. The eigenspectrum is ρ-independent whenever ${\displaystyle \mu /d\approx 0}$. This indicates two factors contributing to the scale invariance of the eigenspectrum. First, a small exponent μ in the kernel function ${\displaystyle f(\overrightarrow{x})}$ means that pairwise correlations slowly decay with functional distance and can be significantly positive across various functional modules and throughout the brain. For a given μ, an increase in dimension ${\displaystyle d}$ improves the scale invariance. The dimension ${\displaystyle d}$ could represent the number of independent features or latent variables describing neural activity or cognitive states.

We verify our theoretical predictions by comparing sampled eigenspectra in finite-size simulated ERMs across different ${\displaystyle \mu }$ and ${\displaystyle d}$ (Figure 4A). We first consider the case of homogeneous neurons (${\displaystyle {\sigma }_i^2\equiv 1}$ in Equation 2, revisited later) in these simulations (Figures 3C, D, 4A), making ${\displaystyle C}$’s entries correlation coefficients. To quantitatively assess the level of scale invariance, we introduce a collapse index (CI, see Methods for a detailed definition). Motivated by Equation 3, the CI measures the shift of the eigenspectrum when the number of sampled neurons changes. The smaller CI values indicate higher scale invariance. Intuitively, it is defined as the area between spectrum curves from different sample sizes (Figure 4A, upper right). In the log–log scale rank plot, Equation 3 shows the spectrum shifts vertically with ${\displaystyle \rho }$.

Figure 4 with 2 supplements see all

Download asset Open asset

Thus, we define CI as this average displacement (Figure 4A, upper right, Methods), and a smaller CI means more scale invariant. Using CI, Figure 4A shows that scale invariance improves with slower correlation decay as ${\displaystyle \mu }$ decreases and the functional dimension ${\displaystyle d}$ increases. Conversely, with large ${\displaystyle \mu }$ and small ${\displaystyle d}$, the covariance eigenspectrum varies significantly with scale (Figure 4A).

Next, we consider the general case of unequal neural activity levels ${\displaystyle {\sigma }_i^2}$ and check for differences between the correlation (equivalent to ${\displaystyle {\sigma }_i^2\equiv 1}$) and covariance matrix spectra. Using the collapsed index (CI), we compare the scale invariance of the two spectra in the experimental data. Intriguingly, the CI of the covariance matrix is consistently smaller (more scale-invariant) across all datasets (Figure 4C, Figure 4—figure supplement 2C, open vs. closed squares), indicating that the heterogeneity of neuronal activity variances significantly affects the eigenspectrum and the geometry of neural activity space (Tian et al., 2024). By extending our spectrum calculation to the intermediate density regime ${\displaystyle \rho {ϵ}^d\ll 1}$ (Methods), we show that the ERM model can quantitatively explain the improved scale invariance in the covariance matrix compared to the correlation matrix (Figure 4—figure supplement 2B; Table 1).

Lastly, we examine factors that turn out to have minimal impact on the scale invariance of the covariance spectrum. First, the shape of the kernel function ${\displaystyle f(\overrightarrow{x})}$ over a small distance (small distance means f(x) near x = 0 in the functional space, Appendix 1—figure 3) does not affect the distribution of large eigenvalues (Appendix 1—figure 3, Table 3, Appendix 1—figure 2, Appendix 1—figure 1A).

This supports our use of a specific ${\displaystyle f(\overrightarrow{x})}$ to represent a class of slow-decaying kernels. Second, altering the spatial distribution of neurons in the functional space, whether using a Gaussian, uniform, or clustered distribution, does not affect large covariance eigenvalues, except possibly the leading ones (Appendix 1—figure 1B, Appendix 1). Third, different geometries of the functional space, such as a flat square, a sphere, or a hemisphere, result in eigenspectra similar to the original ERM model (Appendix 1—figure 1C). These findings indicate that our theory for the covariance spectrum’s scale invariance is robust to various modeling details.

So far, we have focused on random sampling of neurons, but how does the neural activity space change with different sampling methods? To this end, we consider three methods (Figure 5A1): random sampling (RSap), anatomical sampling (ASap) where neurons in a brain region are captured by optical imaging (Grewe and Helmchen, 2009; Gauthier and Tank, 2018; Stringer et al., 2019a), and functional sampling (FSap) where neurons are selected based on activity similarity (Meshulam et al., 2019). In ASap or FSap, sampling involves expanding regions of interest in anatomical space or functional space while measuring all neural activity within those regions (Appendix 1). The difference among sampling methods depends on the neuron organization throughout the brain. If anatomically localized neurons also cluster functionally (Figure 5A4), ASap ≈ FSap; if they are spread in the functional space (Figure 5A2), ASap ≈ RSap. Generally, the anatomical–functional relationship is in-between and can be quantified using the Canonical Correlation Analysis (CCA). This technique finds axes (CCA vectors ${\displaystyle {\overrightarrow{v}}_{\text{anat}}}$ and ${\displaystyle {\overrightarrow{v}}_{\text{func}}}$) in anatomical and functional spaces such that the neurons’ projection along these axes has the maximum correlation, ${\displaystyle R_{\text{CCA}}}$. The extreme scenarios described above correspond to ${\displaystyle R_{\text{CCA}}=1}$ and ${\displaystyle R_{\text{CCA}}=0}$.

Figure 5 with 7 supplements see all

Download asset Open asset

To determine the anatomical–functional relationship in neural data, we infer the functional coordinates ${\displaystyle {\overrightarrow{x}}_i}$ of each neuron by fitting the ERM using multidimensional scaling (MDS) (Cox and Cox, 2000) (Methods). For simplicity and better visualization, we use a low-dimensional functional space where ${\displaystyle d=2}$. The fitted functional coordinates confirm the slow decay kernel function in ERM except for a small distance (Figure 5—figure supplement 3). The ERM with inferred coordinates ${\displaystyle {\overrightarrow{x}}_i}$ also reproduces the experimental covariance matrix, including cluster structures (Figure 5—figure supplement 2) and its sampling eigenspectra (Figure 5—figure supplement 1).

Equipped with the functional and anatomical coordinates, we next use CCA to determine which scenarios illustrated in Figure 5A align better with the neural data. Figure 5B, C shows a representative fish with a significant ${\displaystyle R_{\text{CCA}}=0.327}$ (p-value = 0.0042, Anderson–Darling test). Notably, the CCA vector in the anatomical space, ${\displaystyle {\overrightarrow{v}}_{\text{anat}}}$, aligns with the rostrocaudal axis. Coloring each neuron in the functional space by its projection along ${\displaystyle {\overrightarrow{v}}_{\text{anat}}}$ shows a correspondence between clustering and anatomical coordinates (Figure 5B). Similarly, coloring neurons in the anatomical space (Figure 5C) by their projection along ${\displaystyle {\overrightarrow{v}}_{\text{func}}}$ reveals distinct localizations in regions like the forebrain and optic tectum. Across animals, functionally clustered neurons show anatomical segregation (Chen et al., 2018), with an average ${\displaystyle R_{\text{CCA}}}$ of 0.335 ± 0.054 (mean ± SD).

Next, we investigate the effects of different sampling methods (Figure 5A1) on the geometry of the neural activity space when there is a significant but moderate anatomical–functional correlation as in the experimental data. Interestingly, dimensionality ${\displaystyle D_{\text{PR}}^{\text{ASap}}}$ in data under anatomical sampling consistently falls between random and functional sampling values (Figure 5D). This phenomenon can be intuitively explained by the ERM theory. Recall that for large ${\displaystyle N}$, the key term in Equation 1 is ${\displaystyle {\mathrm{E}}_{i\ne j}(C_{ij}^2)}$. For a fixed number of sampled neurons, this average squared covariance is maximized when neurons are selected closely in the functional space (FSap) and minimized when distributed randomly (RSap). Thus, RSap and FSap ${\displaystyle D_{\text{PR}}}$ set the upper and lower bounds of dimensionality, with ASap expected to fall in between. This reasoning can be precisely formulated to obtain quantitative predictions of the bounds (Methods). We predict the ASap dimension at large ${\displaystyle N}$ as

Here, ${\displaystyle D_{\mathrm{P}\mathrm{R}}}$ is the dimensionality under RSap (Equation 1), ${\displaystyle k}$ represents the fraction of sampled neurons. ${\displaystyle R_{\mathrm{A}\mathrm{S}\mathrm{a}\mathrm{p}}}$ is the correlation between anatomical and functional coordinates along the direction where the anatomical subregions are divided (Methods), and it is bounded by the canonical correlation ${\displaystyle R_{\text{ASap}}\le R_{\text{CCA}}}$. When ${\displaystyle R_{\text{ASap}}=0}$, we get the upper bound ${\displaystyle D_{\text{PR}}^{\text{ASap}}=D_{\text{PR}}}$ (Figure 5D, dashed line). The lower bound is reached when ${\displaystyle R_{\text{ASap}}=R_{\text{CCA}}=1}$ (Figure 5A4), where Equation 4 shows a scaling relationship ${\displaystyle D_{\text{PR}}^{\text{ASap}}=D_{\text{PR}}^{\text{FSap}}\sim k^{2\mu /d}{D}_{\text{PR}}}$ that depends on the sampling fraction ${\displaystyle k}$ (Figure 5D, solid line). This contrasts with the ${\displaystyle k}$-independent dimensionality of RSap in Equation 1. Furthermore, if ${\displaystyle R_{\text{ASap}}}$ and its upper bound is not close to 1 (precisely ${\displaystyle R_{\text{ASap}}\le 0.84}$ for the ERM model in Figure 5D), ${\displaystyle D_{\text{PR}}^{\text{ASap}}}$ align closer to the upper bound of RSap. This prediction agrees well with our observations in data across animals (Figure 5D, Figure 5—figure supplement 6, and Figure 5—figure supplement 7).

Beyond dimensionality, our theory predicts the difference in the covariance spectrum between sampling methods based on the neuronal density ρ in the functional space (Equation 3). This density ρ remains constant during FSap (Figure 5A1) and decreases under RSap; the average density across anatomical regions ${\displaystyle ⟨\rho ⟩}$ in ASap lies between those of FSap and RSap. Analogous to Equation 4, the relationship in ρ orders the spectra: ASap’s spectrum lies between those of FSap and RSap (Methods). This further implies that the level of scale invariance under ASap should fall between that of RSap and FSap, which is confirmed by our experimental data (Figure 5E).

How does task-related neural activity shape the covariance spectrum and brain-wide functional organization? We examine the hunting behavior in larval zebrafish, marked by eye convergence (both eyes move inward to focus on the central visual field) (Bianco et al., 2011). We find that scale invariance of the eigenspectra persists and is enhanced even after removing the hunting frames from the Ca²⁺ imaging data (Figure 4C, Appendix 1—figure 7A, B, Appendix 1). This is consistent with the scale-invariant spectrum found in other datasets during spontaneous behaviors (Figure 5—figure supplement 1F, Figure 2—figure supplement 2G, H), suggesting scale invariance is a general phenomenon.

Interestingly, in the inferred functional space, we observe reorganizations of neurons after removing hunting behavior (Appendix 1—figure 7C, D). Neurons in one cluster disperse from their center of mass (Appendix 1—figure 7D) and decreases the local neuronal density ρ (Appendix 1 and Appendix 1—figure 7E). The neurons in this dispersed cluster have a consistent anatomical distribution from the midbrain to the hindbrain in four out of five fish (Appendix 1—figure 9). During hunting, the cluster has robust activations that are widespread in the anatomical space but localized in the functional space (Appendix 1, Appendix 1—Video 1).

Our findings suggest that the functional space could be defined by latent variables that represent cognitive factors such as decision-making, memory, and attention. These variables set the space’s dimensions, with neural activity patterns reflecting cognitive state dynamics. Functionally related neurons – through sensory tuning, movement parameters, internal conditions, or cognitive factors – become closer in this space, leading to stronger activity correlations.

What drives brain dynamics with a slow-decaying distance–correlation function ${\displaystyle f(\overrightarrow{x})}$ in functional space? Long-range connections and a slow decline in projection strength over distance (Kunst et al., 2019) may cause extensive correlations, enhancing global activity patterns. This behavior is also reminiscent of phase transitions in statistical mechanics (Kardar, 2007), where local interactions lead to expansive correlated behaviors. Studies suggest that critical brains optimize information processing (Beggs and Plenz, 2003; Dahmen et al., 2019). The link between neural correlation structures and neuronal connectivity topology is an exciting area for future exploration.

In the high-density regime of the ERM model, the rank plot (Equation 3) for large eigenvalues (${\displaystyle \lambda >1}$) follows a power law ${\displaystyle \lambda \sim r^{-\alpha }}$, with ${\displaystyle \alpha =1-\mu /d<1}$. The scale-invariant spectrum occurs when α is close to 1. Experimental data, however, align more closely with the model in the intermediate-density regime, where the power-law spectrum is an approximation and the decay is slower (for ERM model, Figure 4—figure supplement 1BC, and for data ${\displaystyle \alpha =0.47\pm 0.08}$, mean ± SD, n =6 fish). Stringer et al., 2019a found an ${\displaystyle \alpha \gtrsim 1}$ decay in the mouse visual cortex’s stimulus trial averaged covariance spectrum, and they argued that this decay optimizes visual code efficiency and smoothness. Our study differs in two fundamental ways. First, we recorded brain-wide activity during spontaneous or hunting behavior, calculating neural covariance from single-trial activity. Much of the neural activity was not driven by sensory stimulus and unrelated to specific tasks, requiring a different interpretation of the neural covariance spectrum. Second, without loss of generality, we normalized the mean variance of neural activity ${\displaystyle \mathrm{E}({\sigma }^2)}$ by scaling the covariance matrix so that its eigenvalues sum up to ${\displaystyle N}$. This normalization imposes a constraint on the spectrum. In particular, large and small eigenvalues may have different behaviors and do not need to obey a single power law ${\displaystyle \lambda \sim r^{-\alpha }}$ for all ${\displaystyle N}$ eigenvalues (Pospisil and Pillow, 2024) (Methods). Stringer et al., 2019a did not take this possibility into account, making their theory less applicable to our analysis.

We draw inspiration from the renormalization group (RG) approach to navigate neural covariance across scales, which has also been explored in the recent literature. Following Kadanoff’s block spin transformation (Kardar, 2007, Meshulam et al., 2019) formed size-dependent neuron clusters and their covariance matrices by iteratively pairing the most correlated neurons and placing them adjacent on a lattice. The groups expanded until the largest reached the system size. The RG process, akin to spatial sampling in functional space (FSap), maintains constant neuron density ρ. Thus, for any kernel function ${\displaystyle f(\overrightarrow{x})}$, including the power law and exponential, the covariance eigenspectrum remains invariant across scales (Appendix 1—figure 5A, B, D, E).

Morrell et al., 2024; Morrell et al., 2021 proposed a simple model in which a few time-varying latent factors impact the whole neural population. We evaluated if this model could account for the scale invariance seen in our data. Simulations showed that the resulting eigenspectra differed considerably from our findings (Appendix 1—figure 6). Although the Morrell model demonstrated a degree of scale invariance under functional sampling (or RG), it did not align with the scale-invariant features under random sampling, suggesting that this simple model might not capture all crucial features in our observations.

We emphasize that the covariance spectrum being a power law is distinct from the scale invariance we define in this study, namely the collapse of spectrum curves under random neuron sampling. The random RNN model in Figure 2I shows a power-law behavior, but lacks true scale invariance as spectrum curves for different sizes do not collapse. When connection strength ${\displaystyle g}$ approaches 1, the system exhibits a power-law spectrum of ${\displaystyle \lambda \propto {\left(\frac{r}{N}\right)}^{-\frac{3}{2}}}$. Subsampling causes the spectrum to shift by ${\displaystyle \lambda \propto k^{-\frac{1}{2}}{\left(\frac{r}{N}\right)}^{-\frac{3}{2}}}$, where ${\displaystyle k=N_s/N}$ is the sampling fraction (derived from Equation 24 in Hu and Sompolinsky, 2022).

The saturation of the dimensionality ${\displaystyle D_{\text{PR}}}$ at large sample sizes indicates a limit to neural assembly complexity, evidenced by the finite mean square covariance. This is in contrast with neural dynamics models such as the balanced excitatory–inhibitory (E–I) neural network (Renart et al., 2010), where ${\displaystyle {\mathrm{E}}_{i\ne j}(C_{ij}^2)\sim 1/N}$ resulting in an unbounded dimensionality (see Appendix 2). Our results suggest that the brain encodes experiences, sensations, and thoughts using a finite set of dimensions instead of an infinitely complex neural activity space.

We found that the relationship between dimensionality and the number of recorded neurons depends on the sampling method. For functional sampling, the dimensionality scales with the sampling fraction ${\displaystyle k:D_{PR}^{FSap}\sim k^{2\mu /d}{D}_{PR}}$. This suggests that if anatomically sampled neurons are functionally clustered, as with cortical neurons forming functional maps, the increase in dimensionality with neuron number may seem unbounded. This offers new insights for interpreting large-scale neural activity data recorded under various techniques.

Manley et al., 2024 found that, unlike in our study, neural activity dimensionality in head-fixed, spontaneously behaving mice did not saturate. They used shared variance component analysis (SVCA) and noted that PCA-based estimates often show dimensionality saturation, which is consistent with our findings. We intentionally chose PCA in our study for several reasons. First, PCA is a trusted and widely used method in neuroscience, proven to uncover meaningful patterns in neural data. Second, its mathematical properties are well understood, making it particularly suitable for our theoretical analysis. Although newer methods such as SVCA might offer valuable insights, we believe PCA remains the most appropriate method for our research questions.

It is important to note that the scale invariance of dimensionality and covariance spectrum are distinct phenomena with different underlying requirements. Dimensionality invariance relies on finite mean square covariance, causing saturation at large sample sizes. In contrast, spectral invariance requires a slow-decaying correlation kernel (small ${\displaystyle \mu }$) and/or a high-dimensional functional space (large ${\displaystyle d}$). Although both features appear in our data, they result from distinct mechanisms. A neural system could show saturating dimensionality without spectral invariance if it has finite mean square covariance but rapidly decaying correlations with functional distance. Understanding these requirements clarifies how neural organization affects different scale-invariant properties.

Our findings are validated across multiple datasets obtained through various recording techniques and animal models, ranging from single-neuron calcium imaging in larval zebrafish to single-neuron multi-electrode recordings in the mouse brain (see Figure 2—figure supplement 2). The conclusion remains robust when the multi-electrode recording data are reanalyzed under different sampling rates (6–24 Hz, Figure 2—figure supplement 4). We also confirm that substituting a few negative covariances with zero retains the spectrum of the data covariance matrix (Figure 2—figure supplement 3 and Methods).

The scale invariance of neural activity across different neuron assembly sizes could support efficient multiscale information encoding and processing. This indicates that the neural code is robust and requires minimal adjustments despite changes in population size. One recent study shows that randomly sampled and coarse-grained macrovoxels can predict population neural activity (Hoffmann et al., 2023), reinforcing that a random neuron subset may capture overall activity patterns. This enables downstream circuits to readout and process information through random projections (Gao et al., 2017). A recent study demonstrates that a scale-invariant noise covariance spectrum with a specific slope ${\displaystyle \alpha <1}$ enables neurons to convey unlimited stimulus information as the population size increases (Moosavi et al., 2024). The linear Fisher information, in this context, grows at least as ${\displaystyle N^{1-\alpha }}$.

Understanding how dimensionality and spectrum change with sample size also suggests the possibility of extrapolating from small samples to overcome experimental limitations. This is particularly feasible when ${\displaystyle \mu /d\to 0}$, where the dimensionality and spectrum under anatomical, random, and functional sampling coincide (Equations 3 and 4). Developing extrapolation methods and exploring the benefits of scale-invariant neural code are promising future research directions.

The handling and care of the zebrafish complied with the guidelines and regulations of the Animal Resources Center of the University of Science and Technology of China (USTC). All larval zebrafish (huc:h2b -GCaMP6f Cong et al., 2017) were raised in E2 embryo medium (comprising 7.5 mM NaCl, 0.25 mM KCl, 0.5 mM MgSO₄, 0.075 mM KH₂PO₄, 0.025 mM Na₂HPO₄, 0.5 mM CaCl₂, and 0.35 mM NaHCO₃; containing 0.5 mg/l methylene blue) at 28.5°C and with a 14-hr light and 10-hr dark cycle.

To induce hunting behavior (composed of motor sequences like eye convergence and J turn) in larval zebrafish, we fed them a large amount of paramecia over a period of 4–5 days post-fertilization (dpf). The animals were then subjected to a 24-hr starvation period, after which they were transferred to a specialized experimental chamber. The experimental chamber was 20 mm in diameter and 1 mm in depth, and the head of each zebrafish was immobilized by applying 2% low melting point agarose. The careful removal of the agarose from the eyes and tail of the fish ensured that these body regions remained free to move during hunting behavior. Thus, characteristic behavioral features such as J-turns and eye convergence could be observed and analyzed. Subsequently, the zebrafish were transferred to an incubator and stayed overnight. At 7 dpf, several paramecia were introduced in front of the previously immobilized animals, each of which was monitored by a stereomicroscope. Those displaying binocular convergence were selected for subsequent Ca²⁺ imaging experiments.

We developed a novel optomagnetic system that allows (1) precise control of the trajectory of the paramecium and (2) imaging brain-wide Ca²⁺ activity during the hunting behavior of zebrafish. To control the movement of the paramecium, we treated these microorganisms with a suspension of ferric tetroxide for 30 min and selected those that responded to its magnetic attraction. A magnetic paramecium was then placed in front of a selected larva, and its movement was controlled by changing the magnetic field generated by Helmholtz coils that were integrated into the imaging system. The real-time position of the paramecium, captured by an infrared camera, was identified by online image processing. The positional vector relative to a predetermined target position was calculated. The magnitude and direction of the current in the Helmholtz coils were adjusted accordingly, allowing for precise control of the magnetic field and hence the movement of the paramecium. Multiple target positions could be set to drive the paramecium back and forth between multiple locations.

The experimental setup consisted of head-fixed larval zebrafish undergoing two different types of behavior: induced hunting behavior by a moving paramecium in front of a fish (fish 1–5), and spontaneous behavior without any visual stimulus as a control (fish 6). Experiments were carried out at ambient temperature (ranging from 23 to 25°C). The behavior of the zebrafish was monitored by a high-speed infrared camera (Basler acA2000-165umNIR, 0.66×) behind a 4F optical system and recorded at 50 Hz. Brain-wide Ca²⁺ imaging was achieved using XLFM. Light-field images were acquired at 10 Hz, using customized LabVIEW software (National Instruments, USA) or Solis software (Oxford Instruments, UK), with the help of a high-speed data acquisition card (PCIe-6321, National Instruments, USA) to synchronize the fluorescence with behavioral imaging.

To validate our findings across different recording methods and animal models, we also analyzed three additional datasets (Table 2). We include a brief description below for completeness. Further details can be found in the respective reference. The first dataset includes whole-brain light-sheet Ca²⁺ imaging of immobilized larval zebrafish in the presence of visual stimuli as well as in a spontaneous state (Chen et al., 2018). Each volume of the brain was scanned through 2.11 ± 0.21 planes per second, providing a near-simultaneous readout of neuronal Ca²⁺ signals. We analyzed fish 8 (69,207 neurons × 7890 frames), 9 (79,704 neurons × 7720 frames), and 11 (101,729 neurons × 8528 frames), which are the first three fish data with more than 7200 frames. For simplicity, we labeled them l2, l3, and l1(fl). The second dataset consists of Neuropixels recordings from approximately ten different brain areas in mice during spontaneous behavior (Stringer et al., 2019b). Data from the three mice, Kerbs, Robbins, and Waksman, include the firing rate matrices of 1462 neurons × 39,053 frames, 2296 neurons × 66,409 frames, and 2688 neurons × 74,368 frames, respectively. The last dataset comprises two-photon Ca²⁺ imaging data (2–3 Hz) obtained from the visual cortex of mice during spontaneous behavior. While this dataset includes numerous animals, we focused on the first three animals that exhibited spontaneous behavior. While this dataset includes numerous animals, we focused on the first three animals that exhibited spontaneous behavior:spont_M150824_MP019_2016-04-05 (11,983 neurons × 21,055 frames), spont_M160825_MP027_2016-12-12 (11,624 neurons × 23,259 frames), and spont_M160907_MP028_2016-09-26 (9392 neurons × 10,301 frames) (Stringer et al., 2019b).

To begin, we multiplied the inferred firing rate of each neuron (see Methods) by a constant such that in the resulting activity trace ${\displaystyle a_i}$, the mean of ${\displaystyle a_i(t)}$ over the nonzero time frames equaled one (Meshulam et al., 2019). Consistent with the literature (Meshulam et al., 2019), this step aimed to eliminate possible confounding factors in the raw activity traces, such as the heterogeneous expression level of the fluorescence protein within neurons and the nonlinear conversion of the electrical signal to Ca²⁺ concentration. Note that after this scaling, neurons could still have different activity levels characterized by the variance of ${\displaystyle a_i(t)}$ over time, due to differences in the sparsity of activity (proportion of nonzero frames) and the distribution of nonzero ${\displaystyle a_i(t)}$ values. Without normalization, the covariance matrix becomes nearly diagonal, causing significant underestimation of the covariance structures.

The three models of covariance in Figure 2G–I were constructed as follows. For model in Figure 2G, the entries of matrix ${\displaystyle G}$ (with dimensions ${\displaystyle N\times T}$) were sampled from an i.i.d. Gaussian distribution with zero mean and standard deviation ${\displaystyle \sigma =1}$. In Figure 2H, we constructed the composite covariance matrix for fish 1 achieved by maintaining the eigenvalues from the fish 1 data covariance matrix and replacing the eigenvectors ${\displaystyle U}$ with a set of random orthonormal basis. Lastly, the covariance matrix in Figure 2I was generated from a randomly connected recurrent network of linear rate neurons. The entries in the synaptic weight matrix are normally distributed with ${\displaystyle J_{ij}\sim \mathcal{N}(0,g^2/N)}$, with a coupling strength ${\displaystyle g=0.95}$ (Hu and Sompolinsky, 2022; Morales et al., 2023). For consistency, we used the same number of time frames ${\displaystyle T=7200}$ when comparing CI across all the datasets (Figures 4B, C, 5D, E, Figure 4—figure supplement 2C). For other cases, we analyzed the full length of the data (number of time frames: fish 1 – 7495, fish 2 – 9774, fish 3 – 13,904, fish 4 – 7318, fish 5 – 7200, and fish 6 – 9388). Next, the covariance matrix was calculated as ${\displaystyle C_{ij}=\frac{1}{T-1}\sum _{t=1}^T\left(a_i(t)-{\overline{a}}_i\right)\left(a_j(t)-{\overline{a}}_j\right)}$, where ${\displaystyle {\overline{a}}_i}$ is the mean of ${\displaystyle a_i(t)}$ over time. Finally, to visualize covariance matrices on a common scale, we multiplied matrix C by a constant such that the average of its diagonal entries equaled one, that is, ${\displaystyle \mathrm{Tr}(C)/N=1}$. This scaling did not alter the shape of covariance eigenvalue distribution, but set the mean at 1 (see also Equation 8).

To maintain consistency across datasets, we fixed the same initial number of neurons at ${\displaystyle N_0=1,024}$. These ${\displaystyle N_0}$ neurons were randomly chosen once for each zebrafish dataset and then used throughout the subsequent analyses. We adopted this setting for all analyses except in two particular instances: (1) for comparisons among the three sampling methods (RSap, ASap, and FSap), we specifically chose 1024 neurons centered along the anterior–posterior axis, mainly from the midbrain to the anterior hindbrain regions (Figure 5DE, Figure 5—figure supplement 6). (2) When investigating the impact of hunting behavior on scale invariance, we included the entire neuronal population (Appendix 1).

We used an iterative procedure to sample the covariance matrix ${\displaystyle C}$ (calculated from data or as simulated ERMs). For RSap, in the first iteration, we randomly selected half of the neurons. The covariance matrix for these selected neurons was a ${\displaystyle N/2\times N/2}$ diagonal block of ${\displaystyle C}$. Similarly, the covariance matrix of the unselected neurons was another diagonal block of the same size. In the next iteration, we similarly created two new sampled blocks with half the number of neurons for each of the blocks we had. Repeating this process for ${\displaystyle n}$ iterations resulted in 2ⁿ blocks, each containing ${\displaystyle N:=N_0/2^n}$ neurons. At each iteration, the eigenvalues of each block were calculated and averaged across the blocks after being sorted in descending order. Finally, the averaged eigenvalues were plotted against rank/${\displaystyle N}$ on a log–log scale.

In the case of ASap and FSap, the process of selecting neurons was different, although the remaining procedures followed the RSap protocol. In ASap, the selection of neurons was based on a spatial criterion: neurons close to the anterior end on the anterior–posterior axis were grouped to create a diagonal block of size ${\displaystyle \frac{N}{2}\times \frac{N}{2}}$, with the remaining neurons forming a separate block. FSap, on the other hand, used the RG framework (Meshulam et al., 2019) to define the blocks (details in Appendix 1). In each iteration, the cluster of neurons within a block that showed the highest average correlation (${\displaystyle {\mathrm{E}}_{i\ne j}(C_{ij}^2)}$) was identified and labeled as the most correlated cluster (refer to Figure 5D, Figure 5—figure supplement 6, and Figure 5—figure supplement 7).

In the ERM model, as part of implementing ASap, we generated anatomical and functional coordinates for neurons with a specified CCA properties as described in Methods. Mirroring the approach taken with our data, ASap segmented neurons into groups based on the first dimension of their anatomical coordinates, akin to the anterier–posterior axis. FSap employed the same RG procedures outlined earlier (Appendix 1).

To determine the overall power-law coefficient of the eigenspectra, α, throughout sampling, we fitted a straight line in the log–log rank plot to the large eigenvalues that combined the original and three iterations of sampled covariance matrices (selecting the top 10% eigenvalues for each matrix and excluding the first four largest ones for each matrix). We averaged the estimated α over 10 repetitions of the entire sampling procedure. ${\displaystyle R^2}$ of the power-law fit was computed in a similar way. To visualize the statistical structures of the original and sampled covariance matrices, the orders of the neurons (i.e., columns and rows) are determined by the following algorithm. We first construct a symmetric Toeplitz matrix ${\displaystyle \mathcal{T}}$, with entries ${\displaystyle {\mathcal{T}}_{i,j}=t_{i-j}}$ and ${\displaystyle t_{i-j}\equiv t_{j-i}}$. The vector ${\displaystyle \overrightarrow{t}=[t_0,t_1,\dots ,t_{N-1}]}$ is equal to the mean covariance vector of each neuron calculated below. Let ${\displaystyle \overrightarrow{c_i}}$ be a row vector of the data covariance matrix; we identify ${\displaystyle \overrightarrow{t}=\frac{1}{N}\sum _{i=1}^{N}D(\overrightarrow{c_i})}$, where ${\displaystyle D(\cdot )}$ denotes a numerical ordering operator, namely rearranging the elements in a vector ${\displaystyle \overrightarrow{c}}$ such that ${\displaystyle c_0\ge c_1\ge \dots \ge c_{N-1}}$. The second step is to find a permutation matrix P such that ${\displaystyle \Vert \mathcal{T}-PC{P}^T{\Vert }_F}$ is minimized, where ${\displaystyle \Vert {\Vert }_F}$ denotes the Frobenius norm. This quadratic assignment problem is solved by simulated annealing. Note that after sampling, the smaller matrix will appear different from the larger one. We need to perform the above reordering algorithm for every sampled matrix so that matrices of different sizes become similar in Figure 2E.

The composite covariance matrix with substituted eigenvectors in Figure 2H was created as described in the following steps. First, we generated a random orthogonal matrix ${\displaystyle U_r}$ (based on the Haar measure) for the new eigenvectors. This was achieved by QR decomposition ${\displaystyle A=U_{r}R}$ of a random matrix ${\displaystyle A}$ with i.i.d. entries ${\displaystyle A_{ij}\sim \mathcal{N}(0,1/N)}$. The composite covariance matrix ${\displaystyle C_r}$ was then defined as ${\displaystyle C_r:=U_r\mathrm{\Lambda }{U}_r^T}$, where ${\displaystyle \mathrm{\Lambda }}$ is a diagonal matrix that contains the eigenvalues of ${\displaystyle C}$. Note that since all the eigenvalues are real and ${\displaystyle U_r}$ is orthogonal, the resulting ${\displaystyle C_r}$ is a real and symmetric matrix. By construction, ${\displaystyle C_r}$ and ${\displaystyle C}$ have the same eigenvalues, but their sampled eigenspectra can differ.

In this section, we introduce the participation ratio (${\displaystyle D_{\text{PR}}}$) as a metric for effective dimensionality of a system, based on Recanatesi et al., 2019; Litwin-Kumar et al., 2017; Gao and Ganguli, 2015; Gao et al., 2017; Clark et al., 2023; Dahmen et al., 2020. ${\displaystyle D_{\text{PR}}}$ is defined as:

Here, ${\displaystyle {\lambda }_i}$ are the eigenvalues of the covariance matrix ${\displaystyle C}$, representing variances of neural activities. ${\displaystyle \mathrm{T}\mathrm{r}(\cdot )}$ denotes the trace of the matrix. The term ${\displaystyle {\mathrm{E}}_{i\ne j}(C_{ij}^2)}$ denotes the expected value of the squared elements that lie off the main diagonal of ${\displaystyle C}$. This represents the average squared covariance between the activities of distinct pairs of neurons.

With these definitions, we explore the asymptotic behavior of ${\displaystyle D_{\text{PR}}}$ as the number of neurons ${\displaystyle N}$ approaches infinity:

This limit highlights the relationship between the PR dimension and the average squared covariance among different pairs of neurons. To predict how ${\displaystyle D_{\text{PR}}}$ scales with the number of neurons (Figure 2D), we first estimated these statistical quantities (${\displaystyle {\mathrm{E}}_{i\ne j}(C_{ij}^2)}$, ${\displaystyle \mathrm{E}({\sigma }^2)}$, and ${\displaystyle \mathrm{E}({\sigma }^4)}$) using all available neurons, then applied Equation 5 for different values of ${\displaystyle N}$. It is worth mentioning that a similar theoretical finding is established by Dahmen et al., 2020. The transition from increasing ${\displaystyle D_{\text{PR}}}$ with ${\displaystyle N}$ to approaching the saturation point occurs when ${\displaystyle N}$ is significantly larger than ${\displaystyle D_{\text{PR}}}$.

We consider the eigenvalue distribution or spectrum of the matrix ${\displaystyle C}$ at the limit of ${\displaystyle N\gg 1}$ and ${\displaystyle L\gg 1}$. This spectrum can be analytically calculated in both high- and intermediate-density scenarios using the replica method (Mézard et al., 1999). The following sketch shows our approach, and detailed derivations can be found in Appendix 2. To calculate the probability density function of the eigenvalues (or eigendensity), we first compute the resolvent or Stieltjes transform ${\displaystyle g(z)=-\frac{2}{N}{\mathrm{\partial }}_z⟨\ln det(zI-C{)}^{-1/2}⟩}$, ${\displaystyle z\in \mathbb{C}}$. Here, ${\displaystyle ⟨\dots ⟩}$ is the average across the realizations of ${\displaystyle C}$ (i.e., random ${\displaystyle {\overrightarrow{x}}_i}$’ s and ${\displaystyle {\sigma }_i^2}$’ s). The relationship between the resolvent and the eigendensity is given by the Sokhotski–Plemelj formula:

where ${\displaystyle \mathbf{I}\mathbf{m}}$ means imaginary part.

Here we follow the field-theoretic approach (Mézard et al., 1999), which turns the problem of calculating the resolvent to a calculation of the partition function in statistical physics by using the replica method. In the limit ${\displaystyle N\to \mathrm{\infty }}$, ${\displaystyle L^d\to \mathrm{\infty }}$, ρ being finite, by performing a leading order expansion of the canonical partition function at large ${\displaystyle z}$ (Appendix 2), we find the resolvent is given by

In the high-density regime, the probability density function (pdf) of the covariance eigenvalues can be approximated and expressed from Equations 6 and 7 using the Fourier transform of the kernel function ${\displaystyle \tilde{f}(\overrightarrow{k})}$:

where ${\displaystyle \delta (x)}$ is the Dirac delta function and ${\displaystyle \mathrm{E}({\sigma }^2)}$ is the expected value of the variances of neural activity. Intuitively, Equation 8 means that ${\displaystyle \lambda /\rho }$ are distributed with a density proportional to the area of ${\displaystyle \tilde{f}(\overrightarrow{k})}$’ level sets (i.e., isosurfaces).

In Results, we found that the covariance matrix consistently shows greater scale invariance compared to the correlation matrix across all datasets. This suggests that the variability in neuronal activity significantly influences the eigenspectrum. This finding, however, cannot be explained by the high-density theory, which predicts that the eigenspectrum of the covariance matrix is simply a rescaling of the correlation eigenspectrum by ${\displaystyle \mathrm{E}({\sigma }_i^2)}$, the expected value of the variances of neural activity. Without loss of generality, we can always standardize the fluctuation level of neural activity by setting ${\displaystyle \mathrm{E}({\sigma }^2)=1}$. This is equivalent to multiplying the covariance matrix ${\displaystyle C}$ by a constant such that ${\displaystyle \mathrm{Tr}(C)/N=1}$, which in turn scales all the eigenvalues of ${\displaystyle C}$ by the same factor. Consequently, the heterogeneity of ${\displaystyle {\sigma }_i^2}$ has no effect on the scale invariance of the eigenspectrum (see Equation 8). This theoretical prediction is indeed correct and is confirmed by direct numerical simulations and quantifying the scale invariance using the CI (Figure 4—figure supplement 2A).

Fortunately, the inconsistency between theory and experimental results can be resolved by focusing the ERM within the intermediate density regime ${\displaystyle \rho {ϵ}^d\ll 1}$, where neurons are positioned at a moderate distance from each other. As mentioned above, we set ${\displaystyle \mathrm{E}({\sigma }^2)=1}$ in our model and vary the diversity of activity fluctuations among neurons represented by ${\displaystyle \mathrm{E}({\sigma }^4)}$. Consistent with the experimental observations, we find that the CI decreases with ${\displaystyle \mathrm{E}({\sigma }^4)}$ (see Figure 4—figure supplement 2B). This agreement indicates that the neural data are better explained by the ERM in the intermediate density regime.

To gain a deeper understanding of this behavior, we use the Gaussian variational method (Mézard et al., 1999) to calculate the eigenspectrum. Unlike the high-density theory where the eigendensity has an explicit expression, in the intermediate density the resolvent ${\displaystyle g(z)}$ no longer has an explicit expression and is given by the following equation:

where ${\displaystyle ⟨\dots {⟩}_{\sigma }}$ computes the expectation value of the term within the bracket with respect to σ, namely ${\displaystyle ⟨\dots {⟩}_{\sigma }\equiv \int \dots p(\sigma )\mathrm{d}\sigma }$. Here and in the following, we denote ${\displaystyle \int \mathrm{D}\overrightarrow{k}\equiv \int \frac{{\mathrm{d}}^d\overrightarrow{k}}{(2\pi {)}^d}}$. The function ${\displaystyle G(\overrightarrow{k},z)}$ is determined by a self-consistent equation,

We can solve ${\displaystyle \int \mathrm{D}\overrightarrow{k}G(\overrightarrow{k},z)}$ from Equation 10 numerically and below is an outline, and the details are explained in Appendix 2. Let us define the integral ${\displaystyle \mathcal{G}\equiv \int \mathrm{D}\overrightarrow{k}\tilde{G}(\overrightarrow{k},z)}$. First, we substitute ${\displaystyle z\equiv \lambda +i\eta }$ into Equation 10 and write ${\displaystyle \mathcal{G}=\mathbf{R}\mathbf{e}\mathcal{G}+i\mathbf{I}\mathbf{m}\mathcal{G}}$. Equation 10 can thus be decomposed into its real part and imaginary part, and a set of nonlinear and integral equations, each of which involves both ${\displaystyle \mathbf{R}\mathbf{e}\mathcal{G}}$ and ${\displaystyle \mathbf{I}\mathbf{m}\mathcal{G}}$. We solve these equations at the limit ${\displaystyle \eta \to 0}$ using a fixed-point iteration that alternates between updating ${\displaystyle \mathbf{R}\mathbf{e}\mathcal{G}}$ and ${\displaystyle \mathbf{I}\mathbf{m}\mathcal{G}}$ until convergence.

We find that the variational approximations exhibit excellent agreement with the numerical simulation for both large and intermediate ρ where the high-density theory starts to deviate significantly (for ${\displaystyle \rho =256}$ and ${\displaystyle \rho =10.24}$, ${\displaystyle ϵ=0.03125}$, Figure 4—figure supplement 1). Note that the departure of the leading eigenvalues in these plots is expected, since the power-law kernel function we use is not integrable (see Methods).

To elucidate the connection between the two different methods, we estimate the condition when the result of the high-density theory (Equation 8) matches that of the variational method (Equations 9 and 10; Appendix 2). The transition between these two density regimes can also be understood (see Equation 22 and Appendix 2).

Importantly, the scale invariance of the spectrum at ${\displaystyle \mu /d\to 0}$ previously derived using the high-density result (Equation 3) can be extended to the intermediate-density regime by proving the ρ-independence using the variational method (Appendix 2).

Finally, using the variational method and the integration limit estimated by simulation (see Methods), we show that the heterogeneity of the variance of neural activity, quantified by ${\displaystyle \mathrm{E}({\sigma }^4)}$, indeed improves the collapse of the eigenspectra for intermediate ρ (Appendix 2). Our theoretical results agree excellently with the ERM simulation (Figure 4—figure supplement 2A, B).

Throughout the paper, we have mainly considered a particular approximate power-law kernel function inspired by the Student’s t distribution

To understand how to choose ${\displaystyle ϵ}$ and ${\displaystyle \mu }$, see Methods. Variations of Equation 11 near ${\displaystyle x=0}$ have also been explored; see a summary in Table 3.

It is worth mentioning that a power law is not the only slow-decaying function that can produce a scale-invariant covariance spectrum (Figure 3—figure supplement 2). We choose it for its analytical tractability in calculating the eigenspectrum. Importantly, we find numerically that the two contributing factors to scale invariance – namely, slow spatial decay and higher functional space – can be generalized to other nonpower-law functions. An example is the stretched exponential function ${\displaystyle f(\overrightarrow{x})=e^{-\Vert \overrightarrow{x}{\Vert }^{\eta }}}$ with ${\displaystyle 0<\eta <1}$. When ${\displaystyle \eta }$ is small and ${\displaystyle d}$ is large, the covariance eigenspectra also display a similar collapse upon random sampling (Figure 3—figure supplement 2).

This approximate power-law ${\displaystyle f(\overrightarrow{x})}$ has the advantage of having an analytical expression for its Fourier transform, which is crucial for the high-density theory (Equation 8),

Here, ${\displaystyle K_{\alpha }(x)}$ is the modified Bessel function of the second kind, and ${\displaystyle \mathrm{\Gamma }(x)}$ is the Gamma function. We calculated the above formulas analytically for ${\displaystyle d=1,2,3}$ with the assistance of Mathematica and conjectured the case for general dimension ${\displaystyle d}$, which we confirmed numerically for ${\displaystyle d\le 10}$.

We want to explain two technical points relevant to the interpretation of our numerical results and the choice of ${\displaystyle f(\overrightarrow{x})}$. Unlike the case in the usual ERM, here we allow ${\displaystyle f(\overrightarrow{x})}$ to be non-integrable (over ${\displaystyle {\mathbb{R}}^d}$), which is crucial to allow power law ${\displaystyle f(\overrightarrow{x})}$. The nonintegrability violates a condition in the classical convergence results of the ERM spectrum (Bordenave, 2008) as ${\displaystyle N\to \mathrm{\infty }}$. We believe that this is exactly the reason for the departure of the first few eigenvalues from our theoretical spectrum (e.g., in Figure 3). Our hypothesis is also supported by ERM simulations with integrable ${\displaystyle f(\overrightarrow{x})}$ (Figure 3—figure supplement 1), where the numerical eigenspectrum matches closely with our theoretical one, including the leading eigenvalues. For ERM to be a legitimate model for covariance matrices, we need to ensure that the resulting matrix ${\displaystyle C}$ is positive semidefinite. According to the Bochner theorem (Rudin, 1990), this is equivalent to the Fourier transform (FT) of the kernel function ${\displaystyle \tilde{f}(\overrightarrow{k})}$ being nonnegative for all frequencies. For example, in 1D, a rectangle function ${\displaystyle \mathrm{rect}(x)=\{\begin{array}{ll}1,& \text{if}|x|\le \frac{1}{2}\\ 0,& \text{otherwise}\end{array}}$ does not meet the condition (its FT is ${\displaystyle \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}(x)=\frac{\sin (x)}{x}}$), but a tent function ${\displaystyle \mathrm{tent}(x)=\{\begin{array}{ll}1-|x|,& \text{ if }|x|\le 1\\ 0,& \text{otherwise}\end{array}}$ does (its FT is ${\displaystyle {\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}}^2(x)}$). For the particular kernel function ${\displaystyle f(\overrightarrow{x})}$ in Equation 11, this condition can be easily verified using the analytical expressions of its Fourier transform (Equation 12). The integral expression for ${\displaystyle K_{\alpha }(x)}$, given as ${\displaystyle K_{\alpha }(x)={\int }_0^{\mathrm{\infty }}{e}^{-x\cosh t}\cosh (\alpha t)dt}$, shows that ${\displaystyle K_{\alpha }(x)}$ is positive for all ${\displaystyle x>0}$. Likewise, the Gamma function ${\displaystyle \mathrm{\Gamma }(x)>0}$. Therefore, the Fourier transform of Equation 11 is positive and the resulting matrix ${\displaystyle C}$ (of any size and values of ${\displaystyle {\overrightarrow{x}}_i}$) is guaranteed to be positive definite.

Building upon the theory outlined above, numerical simulations further validated the empirical robustness of our ERM model, as showcased in Figures 3B–D–4A. In Figure 3B–D, the ERM was characterized by the parameters ${\displaystyle N=1024}$, ${\displaystyle d=2}$, ${\displaystyle L=10}$, ${\displaystyle \rho =10.24}$, and ${\displaystyle \mu =0.5}$ and ${\displaystyle ϵ=0.03125}$ for ${\displaystyle f(\overrightarrow{x})}$. To numerically compute the eigenvalue probability density function, we generated the ERM 100 times, each sampled using the method described in Methods. The pdf was computed by calculating the pdf of each ERM realization and averaging these across the instances. The curves in Figure 3D showed the average of over 100 ERM simulations. The shaded area (most of which is smaller than the marker size) represented the SEM. For Figure 4A, the columns from left to right were corresponded to ${\displaystyle \mu =0.5,0.9,1.3}$, and the rows from top to bottom were corresponded to ${\displaystyle d=1,2,3}$. Other ERM simulation parameters: ${\displaystyle N=4096}$, ${\displaystyle \rho =256}$, ${\displaystyle L=(N/\rho {)}^{1/d}}$, ${\displaystyle ϵ=0.03125}$, and ${\displaystyle {\sigma }_i^2=1}$. It should be noted that for Figure 4A, the presented data pertain to a single ERM realization.

We quantify the extent of scale invariance using CI defined as the area between two spectrum curves (Figure 4A, upper right), providing an intuitive measure of the shift of the eigenspectrum when varying the number of sampled neurons. We chose the CI over other measures of distance between distributions for several reasons. First, it directly quantifies the shift of the eigenspectrum, providing a clear and interpretable measure of scale invariance. Second, unlike methods that rely on estimating the full distribution, the CI avoids potential inaccuracies in estimating the probability of the top leading eigenvalues. Finally, the use of CI is motivated by theoretical considerations, namely the ERM in the high-density regime, which provides an analytical expression for the covariance spectrum (Equation 3) valid for large eigenvalues.

we set ${\displaystyle q_1}$ such that ${\displaystyle \lambda (q_1)=1}$, which is the mean of the eigenvalues of a normalized covariance matrix. The other integration limit ${\displaystyle q_0}$ is set to 0.01 such that ${\displaystyle \lambda (q_0)}$ is the 1% largest eigenvalue.

Here, we provide numerical details on calculating CI for the ERM simulations and experimental data.

Here we briefly explain the CCA method (Knapp, 1978) for completeness. The basis vectors ${\displaystyle {\overrightarrow{v}}_{\text{func}}}$ and ${\displaystyle {\overrightarrow{v}}_{\text{anat}}}$, in functional and anatomical space, respectively, were found by maximizing the correlation ${\displaystyle R_{\text{CCA}}=corr(\{{\overrightarrow{v}}_{\text{func}}\cdot {\overrightarrow{x}}_i\},\{{\overrightarrow{v}}_{\text{anat}}\cdot {\overrightarrow{y}}_i\})}$. These basis vectors satisfy the condition that the projections of the neuron coordinates along them, ${\displaystyle \{{\overrightarrow{x}}_i\cdot {\overrightarrow{v}}_{\text{func}}\}}$ and ${\displaystyle \{{\overrightarrow{y}}_i\cdot {\overrightarrow{v}}_{\text{anat}}\}}$, are maximally correlated among all possible choices of ${\displaystyle {\overrightarrow{v}}_{\text{func}}}$ and ${\displaystyle {\overrightarrow{v}}_{\text{anat}}}$. Here, ${\displaystyle \{{\overrightarrow{x}}_i\}}$, ${\displaystyle \{{\overrightarrow{y}}_i\}}$ represent the coordinates in functional and anatomical spaces, respectively. The resulting maximum correlation is ${\displaystyle R_{\text{CCA}}}$. To check the significance of the canonical correlation, we shuffled the functional space coordinates ${\displaystyle \{{\overrightarrow{x}}_i\}}$ across neurons’ identity and re-calculated the canonical correlation with the anatomical coordinates, as shown in Figure 5—figure supplement 4.