{"id":811317,"date":"2026-05-21T03:23:17","date_gmt":"2026-05-21T03:23:17","guid":{"rendered":"https:\/\/www.europesays.com\/us\/811317\/"},"modified":"2026-05-21T03:23:17","modified_gmt":"2026-05-21T03:23:17","slug":"a-critical-initialization-for-biological-neural-networks","status":"publish","type":"post","link":"https:\/\/www.europesays.com\/us\/811317\/","title":{"rendered":"A critical initialization for biological neural networks"},"content":{"rendered":"

All experimental procedures were conducted according to the Institutional Animal Care and Use Committee (IACUC) at Howard Hughes Medical Institute (HHMI) Janelia. Data analysis and simulations were performed in Python using pytorch and numpy, and figures were made using matplotlib and jupyter-notebooks61<\/a>,62<\/a>,63<\/a>,64<\/a>,65<\/a>.<\/p>\n

Data acquisitionAnimals<\/p>\n

All experimental procedures were conducted according to IACUC ethics approval received from the IACUC board at HHMI Janelia Research Campus. We performed 18 recordings in cortex in: (1) 12 mice bred to express jGCaMP8s66<\/a> in excitatory neurons: TetO-jGCaMP8s \u00d7 Camk2a-tTA mice (available as JAX 037717 and JAX 007004); (2) 3 mice bred to express jGCaMP8s66<\/a> in the somas of excitatory neurons: riboL1\u2013jGCaMP8s \u00d7 Slc17a7-Cre (similar to JAX 039267 without IRES; and JAX 037512) and (3) 3 mice bred to express tdTomato in inter-neurons VGAT-CRE \u00d7 Ai14 (JAX 016962; JAX 007914) with injections of a dual virus Thy1s:TTA (AAV9, 1.64\u00a0\u00d7\u00a01013 vector genomes\u2009ml\u22121) and TRE3G:RiboL1\u2013jGCaMP8s (AAV9, 2.45\u00a0\u00d7\u00a01013 vector genomes\u2009ml\u22121) as in ref. 67<\/a>, see also ref. 68<\/a>. We also performed eight recordings in hippocampal CA1 in six mice bred to express GCaMP6f in excitatory neurons: Thy1-GCaMP6f GP5.17 mice (JAX 025393)69<\/a>, as well as in two wild-type C57 mice (JAX 000664) with injections of the same RiboL1-jGCaMP8s virus combination described above. These mice were male and female, and ranged from 2 to 12\u2009months of age. Mice were housed in reverse light cycle, and were pair-housed with their siblings before and after surgery. Holding rooms are set to a temperature of 70\u2009\u00b0F\u2009\u00b1\u20092\u2009\u00b0F, and humidity of 50% relative humidity\u2009\u00b1\u200920%.<\/p>\n

Surgical procedures<\/p>\n

Surgeries were performed in adult mice (P35\u2013P125) following procedures outlined in refs. 70<\/a>,71<\/a>. In brief, mice were anaesthetized with isoflurane while a craniotomy was performed. Marcaine (no more than 8\u2009mg\u2009kg\u22121) was injected subcutaneously beneath the incision area, and warmed fluids plus 5% dextrose and buprenorphine 0.1\u2009mg\u2009kg\u22121 (systemic analgesic) were administered subcutaneously along with dexamethasone 7\u2009mg\u2009kg\u22121 by the intramuscular route. In the canula implants, the same total dexamethasone dose was administered tapered over 3\u2009days: 4\u2009mg\u2009kg\u22121 on the first day, 2\u2009mg\u2009kg\u22121 on the second day and 1\u2009mg\u2009kg\u22121 on the third day.<\/p>\n

For the visual cortical windows (which included the posterior parietal cortex), measurements were taken to determine the bregma\u2013lambda distance and location of a 4-mm circular window over the V1 cortex, as far lateral and caudal as possible without compromising the stability of the implant. For the sensorimotor windows, the craniotomy was centred at \u22120.75\u2009mm anteroposterior (AP) and 2.2\u2009mm mediolateral (ML) from bregma. A 4- and 5-mm double window was placed into the craniotomy so that the 4-mm window replaced the previously removed bone piece and the 5-mm window lay over the edge of the bone. For the hippocampal windows, the craniotomy was centred at 1.8\u2009mm AP and 2.0\u2009mm ML from bregma. Cortex was aspirated and a 3-mm glass coverslip attached to a stainless-steel was implanted over the dorsal CA1 region. CA1 surgeries were similar to those described in ref. 71<\/a>.<\/p>\n

After surgery, ketoprofen (5\u2009mg\u2009kg\u22121) was administered subcutaneously and the animal allowed to recover on heat. The mice were monitored for pain or distress and ketoprofen 5\u2009mg\u2009kg\u22121 was administered for 2\u2009days following surgery.<\/p>\n

Imaging acquisition<\/p>\n

We used a custom-built 2p mesoscope72<\/a> to record neural activity, and ScanImage73<\/a> for data acquisition. We used a custom online Z-correction module (now in ScanImage), to correct for Z and XY drift online during the recording. As described in ref. 70<\/a>, for the visual area and hippocampal recordings, we used an upgrade of the mesoscope that allowed us to approximately double the number of recorded neurons using temporal multiplexing74<\/a>.<\/p>\n

The mice were free to run on a styrofoam cylinder. Mice were acclimatized to running on the ball for several sessions before imaging, and one mouse was trained on a virtual reality task for 2\u2009weeks before the recording. The field of view was selected such that large numbers of neurons could be observed, with clear calcium transients. Recordings were performed for 100\u2013150\u2009min at a rate of 22\u2009Hz.\u00a0We performed one recording session in each of the 26 mice, and\u00a0did not perform a sample size analysis. Blinding and randomization were not used. Recordings from refs. 43<\/a>,75<\/a> were acquired at a rate of 3\u2009Hz. Recordings from refs. 2<\/a>,76<\/a> were acquired at a rate of 2.5\u20133\u2009Hz, on a Thorlabs Bergamo microscope.<\/p>\n

Processing of calcium imaging data<\/p>\n

Calcium imaging data were processed using the Suite2p toolbox (v.0.9.3)77<\/a>, available at www.github.com\/MouseLand\/suite2p<\/a>. Suite2p performs motion correction, ROI detection, neuropil correction and spike deconvolution as described elsewhere2<\/a>. We used a neuropil subtraction coefficient of 1.0. For the 22-Hz recordings, we used all ROIs output by Suite2p above a signal-to-noise ratio (SNR) threshold of 0.3, which included dendritic processes, to increase the number of units recorded. The SNR for the activity trace x for each ROI was defined as <\/p>\n

$$\\,{\\rm{S}}{\\rm{N}}{\\rm{R}}=1-\\frac{{\\rm{V}}{\\rm{a}}{\\rm{r}}[{x}_{t}-{x}_{t-1}]}{2\\,{\\rm{V}}{\\rm{a}}{\\rm{r}}[{x}_{t}]}$$<\/p>\n

(similar to ref. 78<\/a>). 61\u00a0\u00b1\u00a016% (mean\u00a0\u00b1\u00a0s.d.) ROIs had an SNR greater than 0.3, resulting in a range of 3,981\u201310,595 ROIs across recordings.<\/p>\n

We improved the spike deconvolution here by using the latest version of Suite2p, which will be described in an upcoming manuscript. Our approach was similar to that in ref. 78<\/a>, where a neural network is used to predict the ground-truth spikes from the noisy convolved traces. Unlike ref. 78<\/a> we trained the model on a large number of simulations single-neuron spiking activity convolved with GCaMP-like dynamics and we used a Unet predictive model79<\/a> with a style-vector to capture temporal context independently for each deconvolved trace80<\/a>. We verified that the deconvolved traces better capture neural activity on real data with less noise by evaluating the responses to visual stimuli presented at known times. To estimate the effect of binning in the GCaMP8s recordings (Extended Data Fig. 6e<\/a>), we binned the fluorescence traces by a factor of 7 and then performed deconvolution.<\/p>\n

Neuropixel recordings and processing<\/p>\n

As described in ref. 2<\/a>, eight-neuropixel electrode arrays were used to record simultaneously from up to 3,000 neurons across the brain in three mice60<\/a>. On the day of recording, mice were anesthetized briefly with isoflurane while eight small craniotomies were made with a dental drill. After several hours of recovery, mice were head-fixed in the International Brain Laboratory (IBL) task setup: seated in a plastic tube with their forepaws on a wheel, surrounded by three computer screens in a light-isolated enclosure55<\/a>,81<\/a>. The electrodes were advanced slowly (approximately 10\u2009\u03bcm\u2009s\u22121) to their final depth (4\u2009mm or 5\u2009mm deep), and allowed to settle for around 15\u2009min before recording. During the spontaneous part of the recording, the computer screens were black. Data were pre-processed by re-referencing to the common median across all channels82<\/a>. The probe locations were determined using https:\/\/github.com\/cortex-lab\/allenCCF<\/a>, and the brain mesh in Fig. 2c<\/a> was plotted using this tool, based on Allen Common Coordinate Framework data83<\/a>,84<\/a>.<\/p>\n

These recordings were re-processed with Kilosort4 (v.4.0.14), with default settings85<\/a>. The Kilosort4 spike sorter found 1,756, 2,837 and 2,962 neurons defined as \u2018good\u2019, with a refractory violation rate of less than 0.2 (default), from the three recordings. We excluded neurons with a firing rate of less than 0.01\u2009Hz during the spontaneous recording period, resulting in 1,716, 2,787 and 2,914 neurons in total in each recording. The detected neurons were located across cortex, the hippocampal formation, the striatum and other subcortical areas. For the grouping in Extended Data Fig. 5b,c<\/a>, the \u2018visual cortex\u2019 group consisted of primary and higher-order visual cortices; the \u2018sensorimotor cortex\u2019 group consisted of motor cortex and somatosensory cortex; the \u2018striatum\u2019 group consisted of the caudate putamen and lateral septal cortex and the \u2018subcortical areas\u2019 consisted of superior colliculus, thalamus and midbrain. Mice 1 and 2 had detected neurons in all five groups; mouse 3 had detected neurons in all groups except visual cortex. The spikes were binned at a rate of 22\u2009Hz\u2014the same acquisition rate as the calcium imaging. We also analysed bin sizes from 5\u2009ms to 100\u2009ms in Extended Data Fig. 5d,e<\/a>. The period of spontaneous activity in each recording was 22\u201342\u2009min long (only this period was used for analyses).<\/p>\n

Data analysis<\/p>\n

We normalized the neural activity to avoid fitting single-neuron statistics. We z-scored the activity of each neuron so that the mean activity of each neuron is 0 and its s.d. is 1. We ran Rastermap on the recordings with 100 clusters and 128 PCs, and visualized the sorted activity using 20 neurons per bin38<\/a>. The behavioural running state was estimated by interpolating the running trace to the recording frames, smoothing with a 25 frame Gaussian kernel, and setting a threshold of running\/not running at 1\/100 of the maximum smoothed running speed. We included sessions for analysis in which the mouse was running at least 10% of the time, and for the DMD analysis in which the mouse was running for at least 30\u2009min (Extended Data Fig. 7<\/a>).<\/p>\n

Eigenspectrum estimation<\/p>\n

We estimated the eigenvalues using the covariance between two halves of neurons from the recording. This is to avoid contaminating the eigenspectrum with single-neuron noise produced from the recording methods and from Poisson variability. We divided the recordings in half spatially, using a checkerboard of size 50\u2009\u03bcm for the calcium imaging recordings, and using sections of 40\u2009\u03bcm (eight contacts) on each Neuropixels probe in the electrophysiological recordings. We did not split the data into training and testing timepoints, as we found this inflated the power-law exponent estimates2<\/a> (Extended Data Fig. 4<\/a>).<\/p>\n

We estimated the power-law exponent of the eigenspectrum decay using a weighted linear regression in log space from rank 10 to 500, with weights as the inverse of the log of the rank86<\/a>. The eigenvalue spectra are normalized by the value of the linear regression fit at rank 1. If the length of the spectrum was less than 500 due to limited numbers of neurons, then we estimated the spectrum using 10 to half the length of the spectrum.<\/p>\n

We subsampled the number of neurons randomly in Extended Data Fig. 6f\u2013h<\/a>, and subsampled by brain region in the ephys data in Extended Data Fig. 5b,c<\/a>. We subsampled timepoints in the recordings using chunks of length 50 to 4,000, spaced in time by 4,000 time points.<\/p>\n

Estimating PC timescales from data<\/p>\n

To estimate the timescales of the PCs, we must take into account the noise. A naive estimation of the ACGs would find that later PCs have smaller time-lagged correlations, but that could be due simply to these PCs having lower SNR overall, and thus all their time-lagged correlations would be lower. Instead we take a similar approach to that from SVCA: we split the data into two random subsets of neurons using a checkerboard grid with size 50\u2009\u03bcm, and we calculate the components with a singular value decomposition (SVD) on the covariance between the two subsets of neurons. The resulting left and right singular vectors were used to project test data, and the correlograms were computed between the projections of one set of neurons and the other. The resulting estimates of the PC ACGs were then normalized to 1 at a time lag 0.<\/p>\n

Estimating rotational components from data<\/p>\n

When the dynamics matrix A is not symmetric, its eigenvalues are complex and the covariance of the multi-dimensional data is no longer related directly to these eigenvalues. Thus, to evaluate the complexity or rotational aspect of the dynamics directly, we cannot rely on the PCs alone. Instead, we fit linear predictive models that predict the neural population vector at time t\u00a0+\u00a0dt from the population vector at time t. Such models are referred to typically as DMD, with the small modification that we use a time lag dt\u00a0=\u00a00.23\u2009s instead of the more common dt\u00a0=\u00a01 time sample. The reason to use a time lag is to make estimation of the rotational modes more robust and, in particular, to avoid the potential influence of short timescale artefacts arising from the deconvolution of the calcium imaging data.<\/p>\n

To estimate DMD, we first reduced each dataset to 1,000 dimensions using PC analysis (PCA). We then used ridge regression with a penalty of 0.1 for the ephys and 0.01 for the 2p calcium imaging to predict Xt+dt from Xt in the reduced PCA space. Thus Xt+dt\u00a0\u2248\u00a0BXt, with B a square matrix mapping PCs to PCs. As usual, the neural activity was z-scored on a per neuron basis before applying PCA, but it was not re-normalized afterwards. Since PCA is an orthonormal projection, the eigenvalues of B are the same as would be expected in the full neuronal space, other than estimation errors. In the model, the matrix exponential describes the relation between Xt+dt and Xt, regardless of whether A is symmetric or not: <\/p>\n

$${X}_{t+{dt}}={\\exp }^{(A-I){dt}\/\\tau }{X}_{t}.$$<\/p>\n

Thus, the DMD matrix B we obtained from data is an estimate of \\({\\exp }^{(A-I)dt\/\\tau }\\), at least in the simulations. Looking at the complexity of the eigenvalues of B can thus indicate whether the dynamics are rotational or not. Note that the eigenvalues \\({\\lambda }^{{\\prime} }\\) of \\({\\exp }^{(A-I)dt\/\\tau }\\) are related to the eigenvalues \u03bb of A by \\({\\lambda }^{{\\prime} }={\\exp }^{(\\lambda -1)dt\/\\tau }\\). Thus the higher the complex part of \u03bb, the higher the complex part of \\({\\lambda }^{{\\prime} }\\). We can also now see why taking a larger dt is beneficial: when dt is very small relative to the timescale of the dynamics, the eigenvalues \\({\\lambda }^{{\\prime} }\\) approach 1, making it difficult to estimate their rotational component. The relation \\({\\lambda }^{{\\prime} }={\\exp }^{(\\lambda -1)dt\/\\tau }\\) could be inverted to obtain estimates of \u03bb. However, this multi-step process is likely to contain a lot of estimation error, and we preferred instead to directly compare the estimated distributions of \\({\\lambda }^{{\\prime} }\\) from data with those from appropriately matched simulations. In particular, we compute the number nrot10 of rotations per tenfold attenuation of the complex eigenvector \\({\\lambda }^{{\\prime} }\\): <\/p>\n

$$\\begin{array}{c}\\begin{array}{c}{n}_{\\mathrm{rot}10}=k\\cdot \\mathrm{angle}\\,({\\lambda }^{{\\prime} })\\\\ \\,|{\\lambda }^{{\\prime} }{|}^{k}=0.1\\end{array}\\end{array}.$$<\/p>\n

Thus: <\/p>\n

$${n}_{\\mathrm{rot}10}=\\frac{\\log 0.1}{\\log |{\\lambda }^{{\\prime} }|}\\cdot \\mathrm{angle}\\,({\\lambda }^{{\\prime} }).$$<\/p>\n

We used \u03c4\u00a0=\u00a020\u2009ms in the simulations with symmetric matrices to match approximately the timescales of the data. Longer or shorter \u03c4 in the simulations would simply contract all estimated eigenvalues towards 0, but otherwise leave the number of rotations per tenfold attenuation unchanged.<\/p>\n

We also estimated eigenvalues of dynamics using DMD on various recordings in which rodents and monkeys were performing tasks (Extended Data Fig. 8<\/a>). In terms of rodent experiments, we analysed one session from refs. 87<\/a>,88<\/a>, in which rats ran down a linear track, and 137 neurons from CA1 were recorded using silicon probes. We analysed one session from refs. 89<\/a>,90<\/a> in which mice were detecting a reward location in two different virtual reality corridors, and neural activity was recorded from CA1 using 2p calcium imaging at a rate of 15\u2009Hz. We analysed one session from refs. 91<\/a>,92<\/a> in which mice were performing a visual discrimination task in virtual reality, and neural activity was recorded from visual cortex using 2p calcium imaging at a rate of 3\u2009Hz. For DMD analysis in each of these recordings, we used a ridge penalty of 0.1, and a time delay of 2\u2009s. We also analysed ephys data collected from monkeys, compiled by the Neural Latents Benchmark challenge93<\/a>. We binned the spikes in each of these recordings at a rate of 50\u2009Hz, and ran DMD on the peri-stimulus time histograms (PSTHs) computed from each of the recordings, aligned to movement onset. In refs. 94<\/a>,95<\/a>, the monkey completed a maze task with one to many targets, with 108 different configurations\u2014in the figure we plotted example PSTHs from single target trials. In refs. 96<\/a>,97<\/a>, the monkey performed reaches between random elements of a grid\u2014we binned the reaches based on movement direction angle into 15 bins. In refs. 98<\/a>,99<\/a>, the monkey performed a centre-out reaching task in eight different directions (active trials), and on a subset of trials the joystick was moved (passive trials); recordings were performed in somatosensory Area 2. In each of these recordings, we performed DMD with a ridge penalty of 0.1 and a time delay of 200\u2009ms.<\/p>\n

Local correlation structure<\/p>\n

In Fig. 4<\/a>, we computed the correlations using the recording sampling rate (3\u2009Hz). For the simulations, we used the correlation matrix derived from the eigenvectors. For all the analyses we excluded neuron pairs within 20\u2009\u03bcm of each other. The pairwise correlations were binned in 200-\u03bcm bins (Fig. 4d<\/a>).<\/p>\n

We defined the top 1% of correlations per neuron as the \u2018strong pairs\u2019 for each neuron. We then computed the probability distribution of the strong pairs across spatial bins of 200\u2009\u03bcm. This distribution was normalized by the distribution of all other correlations across bins, producing the strong pair odds versus chance shown in Fig. 4f<\/a>. The strong pair odds, near versus far, was the ratio of the first bin of this curve (within 200\u2009\u03bcm) to the average of the last four bins of this curve (2,200\u20133,000\u2009\u03bcm).<\/p>\n

Dynamical systems analysis<\/p>\n

We assume that neural activity is governed by linear dynamics with independent stochastic inputs. In the case of normally distributed inputs, this model becomes the familiar Ornstein\u2013Uhlenbeck process, with connectivity A33<\/a>: <\/p>\n

$$\\tau \\frac{d{\\bf{x}}}{{dt}}=-{\\bf{x}}+A{\\bf{x}}+{{\\epsilon }}_{{\\bf{t}}},$$<\/p>\n

\n (1)\n <\/p>\n

where the noise is a Wiener process.<\/p>\n

The stationary distribution of the neural covariance matrix \u03a3 is defined by the Lyapunov equation37<\/a>: <\/p>\n

$$(A-I)\\Sigma +\\Sigma {(A-I)}^{{\\rm{\\top }}}=-I.$$<\/p>\n

When A is symmetric, the solution is given by <\/p>\n

$$\\Sigma =\\frac{1}{2}{(I-A)}^{-1}.$$<\/p>\n

\n (2)\n <\/p>\n

For the non-symmetric case, the solution must be calculated numerically100<\/a> (see also ref. 101<\/a> for an example of solving for A).<\/p>\n

For the symmetric case, we can derive the decay of the eigenvalue spectrum directly, under the assumption that A is a random symmetric matrix with an eigenspectrum distribution of a semicircle from \u22121 to 1. We assume here that the number of units is large enough to treat the eigenvalue distribution as an exact semicircle distribution and ignore finite size effects. The rank n of eigenvalue \\({\\lambda }_{A}^{n}\\) of A is defined by the integral of the semicircle distribution density p(\u03bbA) from \\({\\lambda }_{A}^{n}\\) to 1. We scale p(\u03bbA) to a maximum of 1 for convenience of the calculations. If we define \u03b8 as the angle subtended by \\(p({\\lambda }_{A}^{n})\\) on the semicircle, we can use geometric arguments to show that: <\/p>\n

$$\\begin{array}{c}\\begin{array}{c}n=\\pi \\left(\\frac{\\theta }{2\\pi }\\right)-\\frac{1}{2}\\cos \\theta \\sin \\theta \\\\ =\\frac{1}{2}(\\theta -\\cos \\theta \\sin \\theta )\\end{array}\\end{array}.$$<\/p>\n

The eigenvalues of the covariance \u03a3, denoted by \u03bb, are related to the eigenvalues \u03bbA of the connectivity matrix A: <\/p>\n

$$\\lambda =\\frac{1}{2(1-{\\lambda }_{A})}\\Rightarrow \\,{\\lambda }_{A}=1-\\frac{1}{2\\lambda }.$$<\/p>\n

Thus, we have <\/p>\n

$$\\theta ={\\cos }^{-1}({\\lambda }_{A})={\\cos }^{-1}\\left(1-\\frac{1}{2\\lambda }\\right).$$<\/p>\n

Plugging this into the equation for rank n and using \\(\\sin (\\theta )=\\sqrt{1-{\\cos }^{2}(\\theta )}\\) we have: <\/p>\n

$$\\begin{array}{l}2n={\\cos }^{-1}\\left(1-\\frac{1}{2\\lambda }\\right)-\\left(1-\\frac{1}{2\\lambda }\\right)\\sqrt{1-{\\left(1-\\frac{1}{2\\lambda }\\right)}^{2}}\\\\ \\,=\\,{\\cos }^{-1}\\left(1-\\frac{1}{2\\lambda }\\right)-\\left(1-\\frac{1}{2\\lambda }\\right)\\sqrt{\\frac{1}{2\\lambda }}\\sqrt{2-\\frac{1}{2\\lambda }}.\\end{array}$$<\/p>\n

We next used Wolfram Alpha102<\/a> to obtain the Taylor expansion for \\(\\sqrt{2-x}\\) and Puiseux expansion for \\({\\cos }^{-1}(1-x)\\) where x is small. Keeping only terms in which x is raised to a power less than 2, gives <\/p>\n

$$\\begin{array}{l}2n\\approx \\sqrt{\\frac{1}{\\lambda }}+\\frac{1}{24{\\lambda }^{3\/2}}-\\sqrt{\\frac{1}{\\lambda }}+\\frac{5}{8{\\lambda }^{3\/2}}+{\\mathcal{O}}\\left(\\frac{1}{{\\lambda }^{5\/2}}\\right)\\\\ \\,\\approx \\frac{2}{3{\\lambda }^{3\/2}}\\,\\Rightarrow \\lambda \\propto \\frac{1}{{n}^{2\/3}}.\\end{array}$$<\/p>\n

Thus, the eigenvalues decay approximately as a power-law with exponent 2\/3, which is very close to the value we found in simulations.<\/p>\n

As mentioned in the text, this is different from the value of 4\/3 found in ref. 10<\/a>; the discrepancy is due to estimating the \u2018long time window covariance\u2019, which assumes that the data are binned in infinitely long windows. The authors of ref. 10<\/a> argue that the formula converges for short windows (greater than\u00a050\u2009ms), which seems to use the single-neuron timescales as a reference. However, when the dynamical systems are close to critical, the emergent timescales are much longer, and thus a much longer window is needed to reach the stable state. Thus, the derivation in ref. 10<\/a>, although interesting, can capture only the infraslow timescales of neural activity and predicts a 4\/3 power-law on a rank plot in this case.<\/p>\n

We note that the 4\/3 exponent, while not explicitly calculated there, can be derived easily from Eq.\u2009(16) in ref. 10<\/a>: <\/p>\n

$$p(x)=\\frac{\\sqrt{2}}{\\pi }{x}^{-\\frac{7}{4}},$$<\/p>\n

which integrated gives: <\/p>\n

$${\\int }_{0}^{X}p(x)dx\\propto {x}^{-\\frac{3}{4}},$$<\/p>\n

which results in a power-law exponent of approximately 4\/3 in a rank plot10<\/a>.<\/p>\n

Relating timescales to eigenvalues<\/p>\n

In addition to estimating the decay of variances of the PCs, we also want to evaluate dynamical, temporal properties of the data and relate them to the model. In the model, the timescales of the system are related to the eigenvalues of A and therefore to the eigenvalues of \u03a3, following from the matrix exponential solution to the Lyapunov equation: <\/p>\n

$$X(t)={e}^{(A-I)t\/\\tau }X(0)+{\\int }_{0}^{t}{e}^{(I-A)(t-{t}^{{\\prime} })\/\\tau }\\sigma d{W}_{{t}^{{\\prime} }}.$$<\/p>\n

The second term on the right is a noise term that is independent of X(0). Using the SVD decomposition of A\u00a0=\u00a0U\u039bUT, with \u039b having the eigenvalues \\({\\lambda }_{A}^{i}=1-1\/(2{\\lambda }_{i})\\) on the diagonal, the multi-dimensional system can be divided into independent scalar equations: <\/p>\n

$${U}_{i}^{T}X(t)={e}^{({\\lambda }_{{A}^{i}}-1)t\/\\tau }{U}_{i}^{T}X(0)+{\\rm{noise}}.$$<\/p>\n

Thus, the PC component projections \\({U}_{i}^{T}X\\) have an auto-correlation that decays as<\/p>\n

$${e}^{({\\lambda }_{{A}^{i}}-1)t\/\\tau }={e}^{-\\frac{t}{2{\\lambda }_{i}\\tau }},$$<\/p>\n

and, thus, the timescales of the PCs are monotonic with the amplitude of the eigenvalue \u03bbi. For the purposes of the analysis in Fig. 3<\/a>, we need only observe the cross-correlation at time lag t, although its exact decay with \u03bbi cannot be predicted owing to the unknown single-unit timescale \u03c4.<\/p>\n

Simulations of dynamical systems<\/p>\n

We simulated 10,000 neurons governed by the dynamics in Eq.\u2009(1<\/a>). For the dynamics simulations, we performed integration using the Euler\u2013Maruyama method and used a step size of 2\u2009ms and a timescale \u03c4 for each neuron of 20\u2009ms using pytorch62<\/a>,103<\/a>. The random noise was drawn from a Gaussian with a mean of zero and s.d. of one for each neuron at each time step. We ran 80 simulations on a graphics processing unit in parallel, with random initial conditions drawn from a Gaussian with mean 0 and s.d.\u2009=\u20091, each consisting of 60,000 time steps, and discarded the first 4,000 time steps. To replicate the sampling rate in the data, we binned every 23 timepoints (approximately 22\u2009Hz). We also z-scored the unit activities, as in the data. When we visualized the rastermaps of simulations, we ran 40 simulations each consisting of 100,000 time steps, to obtain longer continuous segments of dynamics.<\/p>\n

For testing the eigenspectrum estimation methods (Extended Data Fig. 4<\/a>), we normalized the activity of each neuron by its s.d., applied a relu and then used the activity of the neuron as the mean of a Poisson process, scaled by different values to represent different levels of noise. We multiplied the activity by 0.7, 0.5 and 0.3, representing \u2018low\u2019, \u2018medium\u2019 and \u2018high\u2019 noise levels respectively. To simulate 2p recordings, we used the \u2018medium\u2019 noise Poisson activity traces and convolved each trace with an exponential filter with a decay time constant of 0.25\u2009s. All the activity traces were then scaled by a factor of eight and 400 was added to each trace, and then each trace was scaled individually by a random number drawn from an exponential distribution with 0.001 added to it, to approximate a distribution of SNR values across neurons, with a minimum SNR close to 0.3. We applied shot noise to each trace, by using the trace as the mean of a Poisson process, and then applied deconvolution as in the real neural data. The exponential distribution had mean values of 0.5, 0.2 and 0.08, representing \u2018low\u2019, \u2018medium\u2019 and \u2018high\u2019 shot-noise levels, respectively. The \u2018medium\u2019 shot-noise level had a mean SNR across recordings of approximately 0.56, similar to our neural recordings.<\/p>\n

Dense connectivity<\/p>\n

We drew the excitatory connections between neurons from a uniform random distribution from zero to two. We subtracted off the mean connectivity (one). We set the diagonal values of the matrix to zero. When this matrix is symmetric, its eigenspectrum distribution follows the semicircle law; for the non-symmetric case, it follows the circular law104<\/a>. We divided the matrix by a scalar so that the largest real value of the eigenvalues of the matrix was 0.998, setting A so that it is critically normalized.<\/p>\n

Sparse\/varied connectivity<\/p>\n

We created the sparse symmetric random matrices with random zero or one connections drawn from a Bernoulli with a mean varying from 2.4\u2009\u00d7\u200910\u22124 to 0.25 (Fig. 4a<\/a>). The mean of the Bernoulli was subtracted globally from the entire matrix, resulting in a matrix with zero mean connectivity. The diagonal of the matrix was set to zero. Random sparse symmetric matrices and graphs also follow the semicircle law, when the sparsity is not too high42<\/a>,105<\/a>,106<\/a>,107<\/a>.<\/p>\n

We created the clustered symmetric random matrices by setting the Bernoulli mean to 0.5 for within-cluster (local) connectivity, and the mean out-of-cluster (global) connectivity in a range from 2.4\u2009\u00d7\u200910\u22124 to 0.5 (Fig. 4b<\/a>). This results in a ratio of the probability of local connection to the probability of global connection (P(local)\/P(global)) ranging from 1 to 2,048. Each cluster consisted of 500 neurons. The mean of the Bernoulli distribution for each entry was subtracted (0.5 within-cluster, smaller outside), resulting in mean zero connectivity across the matrix. The diagonal of the matrix was set to zero.<\/p>\n

We created the locally connected symmetric random matrices using as the Bernoulli mean an exponential decay function of Euclidean distances dij in micrometres between neurons in the simulation (Fig. 4c<\/a>). The exponential function was defined as \\(0.5{e}^{-{d}_{ij}\/250}\\). The neurons were placed randomly on a torus of size 8,000\u2009\u00d7\u20098,000\u2009\u03bcm. The minimum value of the mean of the Bernoulli was set to a value from 2.4\u2009\u00d7\u200910\u22124 to 0.5, resulting in a range of P(local)\/P(global) from 1 to 2,048. Again, the mean of the Bernoulli was subtracted from each entry of the matrix, and the diagonal of the matrix was set to zero. To quantify the recovery of true connectivity (Fig. 4e<\/a>), we compute the fraction of true positive connections with P% of top pairwise correlations per neuron, and then normalize by the average probability of positive connection in the simulation (chance).<\/p>\n

As in the dense case, we divided each matrix by a scalar so that the largest eigenvalue of the matrix was 0.998. Because these connectivity matrices were symmetric, we computed the covariance matrix and the eigenspectrum directly from A using Eq.\u2009(2<\/a>).<\/p>\n

Simulations with inputs<\/p>\n

For Fig. 5<\/a> we added external inputs to the simulations described above. All simulation parameters were kept the same, and the input magnitude of m\u2009=\u20092.5 was chosen to be comparable with the amplitude of the ongoing, internally generated activity. Inputs were kept on for 50\u2009ms ending at simulation time 0\u2009s and the readout was averaged over 100\u2009ms of neural activity ending at the delay time on the x\u2009axis of Fig. 5d\u2013f<\/a>. For time-independent decoding, the 100-ms readout ended at a randomized time in the range between 120\u2009ms at a minimum and the time on the x\u2009axis of Fig. 5d\u2013f<\/a> at a maximum. Inputs were assumed to have 100 features, which were projected to the 10,000 neurons in the simulations either with completely random weights (Fig. 5d,e<\/a>), or with weights corresponding to the top 100 eigenvectors of spontaneous activity (Fig. 5f<\/a>). Inputs were presented in pseudo-random order at fixed intervals of 3\u2009s, except for the analyses in Fig. 5f<\/a> where the intervals were 9\u2009s. Mathematically: <\/p>\n

$$\\begin{array}{c}\\tau \\frac{d{\\bf{x}}}{dt}=-{\\bf{x}}+A{\\bf{x}}+d{\\bf{W}}+\\\\ \\,+\\,m\\cdot \\sum _{j}{{\\bf{v}}}_{{\\rm{\\sigma }}(j)}{\\delta }_{t\\in ({t}_{j}-0.05s,{t}_{j})}\\\\ \\,{{\\bf{v}}}_{i}=B{{\\bf{u}}}_{i},\\,{{\\bf{u}}}_{i}\\in {{\\bf{R}}}^{100}\\\\ {{\\bf{y}}}_{j}({t}^{{\\prime} })={\\int }_{{t}_{j}+{t}^{{\\prime} }-0.1s}^{{t}_{j}+{t}^{{\\prime} }}{\\bf{x}}(t)dt,\\end{array}$$<\/p>\n

where dW<\/b> is the Wiener process corresponding to the random noise simulation, v<\/b>\u03c3j is the random input added at time tj and corresponding to a random mapping \u03c3(j) to the range of inputs considered in each simulation. A is the critically normalized dynamics matrix as before and B is a 10,000\u2009\u00d7\u2009100 matrix of either random values or composed of the eigenvectors of the spontaneous activity generated in the absence of inputs. The readout \\({{\\bf{y}}}_{j}({t}^{{\\prime} })\\) for input j is calculated as the 100\u2009ms integral of activity x<\/b>(t) corresponding to a time offset of \\({t}^{{\\prime} }\\) from the input time tj. For the echo-state network version12<\/a>, we simply added a ReLu62<\/a> rectification on x at each Euler step of the same dynamical system. We scaled A to reach the edge of criticality, which in this case required an additional scaling factor of 1.4 to A.<\/p>\n

The Rastermap visualizations in Fig. 5b,c<\/a> were from 1,500\u2009s of simulations with three inputs. There were two randomly drawn 100-dimensional inputs in Fig. 5d<\/a>, with 1,000 repeats of each stimulus, half of which were used for training the decoder and half for testing. A one-nearest-neighbour decoder was used for this binary classification. There were 2,000 different, randomly drawn 100-dimensional inputs in Fig. 5e,f<\/a>, with 1,000 inputs used for training and 1,000 used for testing. To decode the input, we fit a regression problem from the output of the network to the inputs, using ridge regression with L2 regularizer of 10. In this case, a classification was considered correct if the predicted input features were closest to the ground-truth input features, as opposed to any of the other 999 inputs used in the testing set.<\/p>\n

Reporting summary<\/p>\n

Further information on research design is available in the\u00a0Nature Portfolio Reporting Summary<\/a> linked to this article.<\/p>\n","protected":false},"excerpt":{"rendered":"All experimental procedures were conducted according to the Institutional Animal Care and Use Committee (IACUC) at Howard Hughes…\n","protected":false},"author":3,"featured_media":811318,"comment_status":"","ping_status":"","sticky":false,"template":"","format":"standard","meta":{"footnotes":"","_share_on_mastodon":"0"},"categories":[8],"tags":[329685,10046,10047,329686,159,67,132,68],"class_list":["post-811317","post","type-post","status-publish","format-standard","has-post-thumbnail","category-science","tag-cellular-neuroscience","tag-humanities-and-social-sciences","tag-multidisciplinary","tag-network-models","tag-science","tag-united-states","tag-unitedstates","tag-us"],"share_on_mastodon":{"url":"https:\/\/pubeurope.com\/@us\/116610427129517973","error":""},"_links":{"self":[{"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/posts\/811317","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/comments?post=811317"}],"version-history":[{"count":0,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/posts\/811317\/revisions"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/media\/811318"}],"wp:attachment":[{"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/media?parent=811317"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/categories?post=811317"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.europesays.com\/us\/wp-json\/wp\/v2\/tags?post=811317"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}