Functional groups of mouse visual neurons.

(A) Diagram of multi-region two-photon imaging of mouse V1 and HVAs, using a custom wide field-of-view microscope. Example imaging session of the simultaneous recording session of V1, LM, AL, and PM. Squares indicate 500 µm wide imaging regions. (B) Example responses from two neurons (mean calcium trace) to drifting gratings with eight directions at various SF-TF frequencies. (C) Neurons were distributed into 65 different classes using GMM (Fig. S1, S2). The mean correlation coefficients of the center of each class (in principal component space) between GMMs of 10 permutations of a random subset of neurons. (D) The confusion matrix shows that individual neurons are likely (>90%) to remain in the same class even when only a random subset of neurons is used to train the GMM (horizontal), compared to the full data set (vertical). (E) Center of individual neurons (left) overlay on an average visual cortex map. The average visual cortex map was generated by affine registration of visual area maps from all experiments. Neurons are colored by visual areas. Middle, average preferred TF exhibits spatial dependency over the visual cortex (TF: A→P, cor = −0.25, p =0.015, M→L, cor = 0.36, p = 0.0004). Right, the average preferred SF (right) exhibits spatial dependency over the visual cortex (SF: A→P, cor = 0.35, p = 0.0005, M→L, cor = −0.06, p = 0.54). Colored dots indicate the average TF and SF (computed with >30 neurons) within patches (180 µm x 180 µm local areas), overlaid on a map of V1 and HVAs. (F) These 65 classes were manually arranged into six tuning groups based on spatial frequency and temporal frequency (SF-TF) tuning preferences. Column 1, the fraction of neurons in different SF-TF groups. Dots represent individual sessions. Statistical significance was tested by the Ranksum test (*, p < 0.05, **, p < 0.01). Column 2, the characteristic SF-TF responses of each tuning group. Column 3, speed tuning of tuning groups. Column 4, distribution of cells’ orientation selectivity index (OSI) and direction selectivity index (DSI). The number of neurons belonged to the six tuning groups combined: V1, 5373; LM, 1316; AL, 656; PM, 491; LI, 334 (refer to Methods for neuron selection). These six groups provide a compact way of summarizing response diversity, but as shown later, the granularity of the 65 classes provides a superior match to the network properties (Fig. 4F).

Noise correlation measurements are reliable

(A) (left) An example covariance matrix of a simultaneous recording of neuronal activity in V1 and PM. (right) NC histogram of an imaging session (blue, note the large positive tail), compared to control, trial-shuffled data (red). (B) The observed population mean NC is always larger than control values (the NC after trial shuffling). Circles indicate the value of individual experiments. (C) Population NCs are precise measures. The margin of error at 95% confidence interval (CI) of the population mean NC reduces rapidly with an increasing number of neuron pairs (72-time bins x 10 trials). With the experimental size of the population (>100 neuron pairs), the estimation precision surpasses the 0.01 level. (D) (left) The NC of individual neuron pairs can be computed using different random subsets of trials, yet reliably converges on similar values (right) The variance of NC computed using a subset of trials is explained by the variance in the held-out subset of trials (53 ± 24 % variance explained; total 204 populations). Each subset contains half of all the trials. The variance explained is defined as the R2 of the linear model.

Factors that contribute to mesoscale NC.

(A) The variance of within-area and inter-area NCs (during grating stimuli) is explained by individual factors including firing rate (FR), signal correlation (signal corr.), neuron distance (neu. dist.), and receptive field (RF) overlap. (B) NCs of neurons with non-overlapping RF is modulated by orientation tuning similarity (within-area, PV 1 < 10−4 (N = 3401), PLM = 0.03 (N = 181), PAL = 0.019 (N = 284); inter-area: PV 1−LM = 0.019 (N = 650), PV 1−AL = 0.0004 (N = 998); t-test). (C) The variances of within-area NCs and inter-area NCs are explained by FR, signal corr., neu. dist., and RF overlap combined. The variance explained is a multi-linear regression model’s R2. (A, C) The error bar indicates the standard error of the mean of permutations. A subset of 100 neuron pairs was randomly selected for each permutation.

High-fidelity tuning-specific V1-HVA communication channels

(A) The distribution of NCs of a subset of LI neurons from tuning group 1 (blue) and NCs of a subset of LI neurons between group 1 and other tuning groups (red). In the next panel, we will focus on the positive tail: the portion of the distribution that is over 2.5x of the S.D. (B) Fraction of neuron pairs with high NCs (> 2.5*S.D. of trial-shuffled NCs) for within-group and inter-group pairs. Neurons within a group have a larger fraction of neuron pairs exhibiting high-fidelity connections (all comparisons, t-test, p < 0.0001). Distributions were generated with 100 permutations. (C) The fraction of high-fidelity connections is linearly related to the fraction of SF-TF groups in each HVA (r2 = 0.9; p < 0.0001). The X-axis indicates the number of high NC neuron pairs in each group in the simultaneously imaged V1-HVA populations divided by the total number of high NC pairs in the V1-HVA population imaged (sum over all groups). (right) This result is summarized in a diagram indicating that area-specific SF-TF biases correlate with the number of high-fidelity functional connections. (D) The average NC value (mean strength of connection) for each tuning group is not linearly related to the fraction of SF-TF groups in each HVA (r2 = 0.1; p = 0.1). (right) This result negates the hypothesis suggested by the diagram, where area-specific SF-TF biases correlate with the strength of functional connections. Instead, the number of connections (panel C) seems to account for the observed trends. (E) Density function plots of NCs for in-area (left) or inter-area (right) neuron pairs that shared the same GMM-based class (65 classes) or group (six SF-TF preference groups) indicate that the more granular, GMM-based class categorization accounts for the structure of the NC network with higher fidelity than the coarser SF-TF groups (full scale is inset, bin size is 0.00875). (F) The functional connectivity matrix for the V1-HVA network between GMM classes exhibits a modular structure. (right) Each module has a particular tuning selectivity and SF-TF bias. The value of this functional connectivity matrix was the fraction of high NC pairs.

Noise correlations (NCs) across different classes of stimuli are more stable than tuning, or signal correlations (SCs).

(A) (left) NCs measured during the naturalistic video are well correlated with NCs measured during drifting grating stimuli. The correlation between NCs across different stimuli is significantly higher than the correlation between corresponding SCs (corr(NCgrat, NCnat) = 0.22 ± 0.13; corr(SCgrat, SCnat) = 0.084 ± 0.065; t-test, p < 0.0001). Colored circles represent individual experiments. Gray dots represent trial-shuffled control (corr(NCshulfgrat, NCshulfnat) = 0.02 ± 0.06). The black/gray dot and error bars indicate the mean and SD for NC and SC. (right) The correlation between (top) NCs and (bottom) SCs during grating and naturalistic video stimuli in an example dataset (red arrowhead in the lefthand plot). (B) NCs to a naturalistic video are positively related to the SCs, as well as to the NCs to drifting gratings. The shaded area indicates SEM. (C) The percentage of NCnat variance is explained by a linear model of SCnat, NCgray, or both factors. NCgray is a better linear predictor compared to SCnat (NCgray, 5.3 ± 3%; SCnat, 4 ± 2%; t-test, p < 0.0001). Combining both factors predicts the NCnat even better (8 ± 3%; t-test, p < 0.0001). Variance explained is measured by R2 of the linear regression.

A network simulation shows that recurrent connectivity can contribute to the stability of the NC network.

(A) A model with two leaky integrate and fire (LIF) neurons that are connected through excitatory synapses. The LIF neurons receive a fraction of shared input (red) and independent input (green) from a Poisson input layer. (B) The firing rate (left) and NC (right) of the two LIF neurons in a toy model (A) is regulated by the fraction of shared input and the strength of the recurrent connection. (C) Schematic of an LIF neuron network model with randomly connected LIF neurons and an input Poisson layer. The structure of the input connection and the strength of the recurrent connection are modulated in the simulation (D, E). (D) In networks with random input connection structures, increasing recurrent connection strength leads to higher cross-stimulus stability of the NC network. Among the values tested, recurrent connection = 0.2 (red) generated a network that was closest to the mouse L2/3 visual neurons (black). (E) In simulations with 0.2 recurrent connectivity strength, regulating the input structure does not change the cross-stimulus stability of the NC network but leads to higher cross-stimulus stability of the SC network. (F) Same LIF network as C with an additional source of top-down input. The dimensionality and strength of the top-down input are modulated in the simulation (G, H). (G) In a network that is mainly driven by sensory input and receives moderate top-down modulation, the sensory/top-down input ratio is 10, recurrent connection is required to reproduce the NC network’s cross-stimulus pattern. Changing the dimensionality of top-down input generated similar NC network patterns (overlapping dots in the figure). (H) In networks that receive strong top-down input, the sensory/top-down input ratio is no greater than 1.6, top-down modulation alone can generate cross-stimulus stable NC. Error bars in all panels indicate the SD of multiple randomly initialized simulations under the same condition.

Functional groups by multi-region two-photon calcium imaging.

(A) Example intrinsic signal imaging of mouse visual areas. (B) Moving square stimuli for quick RF mapping. (C) Example population RFs of simultaneously imaged populations. Blue and orange contours indicate the Gaussian profile of population RF of neurons from different visual areas, and blue shade indicates the overlap region of population RF of two simultaneous imaging regions. Values indicate the fraction of overlap. Upper right: example population RF of a quartic-region imaging. Lower right: summarize the fraction of population RF overlap of individual experiments (gray circle). Error bars indicate the mean and standard deviation. (D) Upper: short and long axes of the Gaussian profile of single neuron RF of all tested HVA neurons are longer than that of V1 (short, p < 0.0001; long, p < 0.0001; One-way ANOVA with Bonferroni correction). Bottom: population RFs of HVA are significantly larger than that of V1 (FWHM, p = 0.0003; Size: p < 0.0001. one-way ANOVA with Bonferroni correction). (E) The responsiveness of V1 and HVAs to the 72-condition sine-wave drifting grating stimuli. Left: the fraction of responsive neurons in HVAs is not significantly different (trial-to-trial Pearson correlation > 0.08; one-way ANOVA, p = 0.36). Right: distribution of neuron firing reliability (trial-to-trial Pearson correlation of inferred spike train at 500 ms bin). Only responsive neuron was considered. V1 and LM were slightly more reliable than AL, PM, and LI (one-way ANOVA with Bonferroni multiple comparisons, p = 1.7 × 10-7). (F) Number of neurons of each GMM class. Classes 26 and 34 were left unclassified, as they were largely insensitive to stimulus parameters, and classes 59 and 64 were small and thus excluded from further analysis. (G) The confusion matrix shows the joint probability of a neuron being identified as group A in GMM of randomly ordered data (shuffle the order of neurons, horizontal), and the neuron is classified as group B in GMM of the original data set (vertical). The diagonal indicates the probability of neurons being classified into the same group. The confusion matrix was generated by averaging a joint probability of 10 permutations. (H) The bar chart shows the probability of correctly allocating neurons into the same group in 10 permutations of GMMs of randomly ordered data.

t-SNE embedding of GMM classes.

(A) t-SNE embedding of all the neurons for the SF-TF analysis. Axes are arbitrary units. Neurons are colored by their GMM class ID (inserted number) and are organized into 6 groups by SF-TF preference. A polar plot shows the average tuning curve for the class is shown in the same color. (B) Four GMM classes are excluded due to being too small or having no characteristic feature selectivity. The response feature of each class is described by three panels: a polar plot shows the average tuning curve for cells in the class; the middle panel shows the normalized response to different joint combinations of TF (x-axis, Hz) and SF (colored line, blue 0.02 cpd, red 0.05 cpd, yellow 0.19 cpd); right panel shows the normalized response to different speed of gratings (x-axis, deg/s).

Spatial modulation on SF-TF and orientation tuning.

(A) Polar plots of averaged preferred directions of six tuning groups of V1 and HVAs. Polar plots were generated with >30 neurons. Black and gray lines indicate the mean and SEM of normalized preferred directions. (B) The heatmap shows the fraction of orientation (OSI > 0.85, left) and direction-selective neurons (DSI > 0.85, right) of each tuning group in V1, LM, AL, and PM. (C) Mean preferred SF (left) and TF (right) of tested visual areas. AL prefers lower SF and higher TF than the other tested visual areas (p < 0.0001, one-way ANOVA with Bonferroni multi-comparison).

Tolerance of noise correlation to missing spikes.

(A) Left, comparing the inferred spike train and ground truth spike train (cell-attached recording) of one example neuron. Spike inference recovered 50% of the spikes of this neuron, the linear correlation between the inferred spike train and the true spike train is 0.79 (bin 1 s). Right, the correlation between the inferred spike train and the true spike train at various time bins. (B) The inter-neuron cross-correlation computed by true spike train and inferred spike are linearly correlated (r = 0.7). (C) The ground truth spike trains (top) and spike train after different types of perturbations of example neurons from the spikefinder dataset (Methods). (D) Compute correlation of residual spike count at 1 s time bin after spike perturbations from left to right: random missing spikes; missing isolated spikes with inter-spike-interval (ISI) > 0.03 s; missing all spikes within a burst ISI < 0.01 s; missing 60% spikes within a burst with ISI <0.02 s. (E) Fidelity (left) and variance explained (right) of correlation calculation with spike train perturbation. The fidelity was defined as the linear correlation between spike count correlation before and after perturbation. Variance explained was measured as r2 of a linear regression between true correlation and perturbed correlations. The colored text in the figure indicates the ISI thresholds.

Factors contribute to the variance of NCs.

(A) Within- and inter-area noise correlations are positively related to signal correlation. (B) Within-area (left) and inter-area (right) NC is significantly higher in neuron pairs with shared RF (within-area, PV 1 < 0.0001, PLM = 0.98, PAL = 0.006, PPM = 0.03, PLI < 0.0001; inter-area:PV 1−LM = 0.3, PAL = 0.0007, PV 1−PM = 0.15, PV 1−LI = 0.82; t-test). Overlapping groups and non-overlapping groups are defined as neuron pairs share > 60% RF, and <20% RF, respectively. (C) Noise correlations of V1 and HVAs are positively related to joint spike count (For all within- and inter-area correlation, r = 0.09-0.18, p < 0.0001). Mean joint spike count is the geometric mean of the spike count to all stimuli. (D) Plot within- and inter-area noise correlation as a function of difference in preferred orientation. Only orientation-selective neurons (OSI > 0.5) were included.

Distance-dependence of inter-area NC explained by retinotopic map.

(A) The distance-dependent decrease of within-area NCs (blue) was not different among V1, LM, AL, and PM with a linear decreasing rate of NC = −0.02 ± 0.009 mm−1 (mean ± S.D.). In these areas, SCs (purple) also exhibited significant distance-dependent decreases. (B) Inter-area NC (blue) and SC (purple) exhibited region-specific distance-dependence patterns that were often non-monotonic. (A-B) Solid lines and shaded areas indicate the mean and standard error of the mean. Each distance bin contains >50 data points. Pearson correlation r is shown in the figure and stars indicate significance (*, p < 0.05; **, p < 0.01; ***, p < 0.0001). (C) Example affine transformation of ISOI maps. The left shows the original V1 map, the middle is the V1 map after affine transformation and the right is the original LM map. (D) Left, is a cartoon of two recurrent layers with an aligned retinotopic map. Right: neuron location on the visual cortex before and after warping. To examine whether the distance-dependent increase of V1-LM NCs to drifting gratings could be explained by retinotopy, we aligned V1 and LM retinotopy by affine transformations of intrinsic signal retinotopic maps of V1 to match those of LM. After retinotopic warping, V1 and LM were treated as two layers with aligned retinotopies (E) Distance-dependent increasing of V1-LM NC to sine-wave drifting gratings before (left) and after retinotopic warping (right). Individual experiments with significant distance dependence are shown in colored curves. The black curve shows the population mean and standard error (Pearson correlation, before warping, r = 0.066, p < 0.0001; after warping, r = −0.026, p < 0.0001). (F) After warping, the inter-area NCs between V1 and LM exhibited a distance-dependent decrease (linear correlation, r = −0.26, p = 3 × 10-6), similar to the trends found for within-area NCs in V1 and LM. (G) Example affine transformation of ISOI maps of V1 and AL. The left shows the original V1 map, the middle is the V1 map after affine transformation and the right is the original AL map. (H) Distance-dependent decreasing of inter V1-AL NC after retinotopic alignment (linear correlation, r = −0.05, p < 0.0001).

Connectivity between GMM classes

(A) The density function of within-area (left) and inter-area (right) NC for neuron pairs from the same tuning group, or the same GMM class, shared high SC, or from the local neighborhood. The SC or the distance threshold is defined such that the same number of neuron pairs is included as the in-class neuron pairs. (B) The spread of neurons from GMM class 1 on a registered map for visual areas. (C) The modularity of the V1-HVA connectivity between GMM classes is regulated by a spatial smooth parameter γ. We reported the modular structure of the connectivity matrix when γ = 0.85, as it generated the largest deviation from a random connectivity matrix. The left side shows the modular structure of V1-HVA GMM classes (upper), and that of a random matrix preserving the degree distribution (lower).

Gain modulation does not explain the NC connectivity

(A) Example time series of pupil dynamics: the pupil area (pixel), x,y position of the pupil centroid, and the first PC of the pupil video. (B) Decomposition of the variance of population neural activity to stimulus, pupil dynamics, and network interaction. Pupil dynamics explained less variance than stimulus and network interactions (***, p < 0.0001, paired t-test). A multi-dimensional pupil dynamics model explains more variance than a 1D model (Multi-D 3.5 ± 2%, 1D 1.5 ± 1%, p < 0.0001, paired t-test). The rank of the multi-dimensional model is the best-fit model of reduced rank regression, the rank dimension is 13.4 ± 3. The 1D model is generated by the most pronounced eigenvector of the reduced rank regression model. (C) To examine the relative contributions of stimulus and pupil dynamics, we analyzed each experiment and plotted a histogram of the ratio of neural activity variance explained by stimulus input and pupil dynamics (arousal modulation). Note that the x-axis is in log10 scale. The red line indicates the mean of the distribution, which indicates that the stimulus modulates neuronal activity approximately 10-fold more than pupil dynamics / arousal. The analysis was only performed on the subset of experiments with pupil recordings (16 mice, 34 recording sessions, total of 91 cortical areas). (D) We reproduced the qualitative results in Figure 5A using only the residual neural activity after subtracting the arousal-modulated portion. The analysis was only performed on the subset of experiments with pupil recordings.

Entire data list.

The Animal ID is a simple identifier number. Note that some animals were used for multiple imaging configurations. On the left is the information for mice that were imaged during viewing of the drifting grating visual stimuli only. On the right is the information for the mice that were imaged both during viewing of gratings and during viewing of the naturalistic video stimuli. On the bottom right are summary figures for the total numbers of animals, neurons, and unique neuron pairs (imaged simultaneously to permit the computation of noise correlations).