Functional connectomics reveals general wiring rule in mouse visual cortex

Abstract

Understanding the relationship between circuit connectivity and function is crucial for uncovering how the brain computes. In mouse primary visual cortex, excitatory neurons with similar response properties are more likely to be synaptically connected^1-8; however, broader connectivity rules remain unknown. Here we leverage the millimetre-scale MICrONS dataset to analyse synaptic connectivity and functional properties of neurons across cortical layers and areas. Our results reveal that neurons with similar response properties are preferentially connected within and across layers and areas-including feedback connections-supporting the universality of 'like-to-like' connectivity across the visual hierarchy. Using a validated digital twin model, we separated neuronal tuning into feature (what neurons respond to) and spatial (receptive field location) components. We found that only the feature component predicts fine-scale synaptic connections beyond what could be explained by the proximity of axons and dendrites. We also discovered a higher-order rule whereby postsynaptic neuron cohorts downstream of presynaptic cells show greater functional similarity than predicted by a pairwise like-to-like rule. Recurrent neural networks trained on a simple classification task develop connectivity patterns that mirror both pairwise and higher-order rules, with magnitudes similar to those in MICrONS data. Ablation studies in these recurrent neural networks reveal that disrupting like-to-like connections impairs performance more than disrupting random connections. These findings suggest that these connectivity principles may have a functional role in sensory processing and learning, highlighting shared principles between biological and artificial systems.

Fig. 1. Overview of MICrONS dataset.

a, Depiction of functionally characterized volumes (left; GCaMP6s in green, vascular label in red) and EM data (right). Middle top, the overlap of the functional 2P (green) and structural EM (grey) volumes from which somas were recruited. Middle bottom, an example of matching structural features in the 2P and EM volumes, including a soma constellation (dashed white circles) and unique local vasculature (red arrows), used to build confidence in the manually assigned 2P–EM cell match (central white circle). All MICrONS data are from a single mouse. Scale bars, 5 μm. b, Deconvolved calcium traces from 100 imaged neurons. Alternating blue and white column overlay represents the duration of serial video trials, with sample frames of natural videos depicted below. Parametric stimuli (not pictured) were also shown for a shorter duration than natural videos. c, Schematic of the digital twin deep recurrent architecture. During training, video frames (left) are input into a shared convolutional deep recurrent core (orange and blue layers), resulting in a learned representation of local spatiotemporal stimulus features. Each neuron is associated with a location (spatial component) in the visual field (grey layer) to read out feature activations (shaded blue vectors), and the dot product with the neuron-specific learned feature weights (shaded lines, feature component) results in the predicted mean neural activation for that time point. CvT-LSTM, convolutional vision transformer, long short-term memory d, Top, depiction of 148 manually proofread mesh reconstructions (grey), including representative samples from layer 2/3 (red), layer 4 (blue), layer 5 (green) and layer 6 (gold). Bottom, presynaptic soma locations relative to visual area boundaries.

Fig. 2. Neurons with higher signal correlation are more likely to form synapses.

a, Anatomical control selection schematic. For each presynaptic neuron (yellow), true postsynaptic partners (black) have controls drawn from unconnected neurons with non-zero axon–dendrite co-travel distance (ADP, red) or zero co-travel distance in the same cortical region (blue). b, Meshes showing presynaptic (yellow axon), postsynaptic (black dendrite) and ADP control (red dendrite) neurons. c, Presynaptic axons in EM space for all projection types (V1→V1, HVA→HVA, V1→HVA and HVA→V1), with soma centroids of connected partners (black), ADP controls (red), same-area controls (blue) and all other functionally matched neurons (grey). Orange triangles represent presynaptic soma. Dashed line is the V1–HVA boundary. Scale bars, 100 μm. Nucleus IDs: V1, 327859; HVA, 560530. d, Mean signal correlation differs between synaptic partners, ADP controls and same-region controls across projection types. Data are mean ± s.e.m.; two-sided paired t-test. Sample size in Supplementary Table 2. e, L_d increases with signal correlation (Δ co-travel distance and Δ signal correlation represent deviations from mean per presynatic neuron). V1→V1: mean L_d = 9.03 μm; HVA→HVA: mean L_d = 9.83 μm; V1→HVA: mean L_d = 4.17 μm; HVA→V1: mean L_d = 1.53 μm. Bands represent bootstrapped s.e.m. Sample sizes for GLMM statistics are shown in Supplementary Tables 3 and 4. f, Synapse density (N_syn/L_d) increases with signal correlation across projections. V1→V1: mean N_syn/L_d = 1.12 mm⁻¹; HVA→HVA: mean N_syn/L_d = 0.83 mm⁻¹; V1→HVA: mean N_syn/L_d = 1.55 mm⁻¹; HVA→V1: mean N_syn/L_d = 1.26 mm⁻¹. Bands represent bootstrapped s.e.m. Sample sizes for GLMM statistics are shown in Supplementary Tables 5 and 6. g, Meshes with small (896 voxels) and large (41,716 voxels) synapse cleft volumes. h, Synapse size (log₁₀ cleft volume in voxels) correlates with signal correlation (6,608 pairs, P = 3.997 × 10⁻²¹, linear regression on unbinned data). Residual signal correlation adjusted for L_d. Bars show bin-wise s.e.m. i, Multisynaptic meshes (yellow, presynaptic; black, postsynaptic). j, Signal correlations increase with synapse count (6,608 pairs, P = 0.009, linear regression on unbinned data). Residual signal correlation adjusted for L_d. Bars show bin-wise s.e.m. *P < 0.05, **P < 0.01, ***P < 0.001 for all figures; corrected for multiple comparisons by Benjamini–Hochberg procedure.

Fig. 3. Feature weight similarity predicts synaptic selectivity more accurately than RF centre distance.

a, Axon–dendrite co-travel distance increases with feature weight similarity and decreasing RF centre distance for within-area (V1→V1 and HVA→HVA), feedforward (V1→HVA) and feedback (HVA→V1) connectivity. Bands represent bootstrapped s.e.m. Sample sizes for GLMM statistics are shown in Supplementary Tables 7–10. b, Synapse density increases with feature weight similarity, but not with RF distance, except for HVA→V1 projections. Bands represent bootstrapped s.e.m. Sample sizes for GLMM statistics are shown in Supplementary Tables 11–14. c, Multiple synapses are associated with increasing feature similarity, but not RF centre distance for 6,608 pairs of connected neurons, after regressing out L_d. Error bars represent s.e.m. Feature weight similarity: P =  0.003; RF distance: P = 0.358; by linear regression. d, Only feature similarity (and not RF centre distance) is associated with an increase in cleft volume for 6,608 pairs of connected neurons after regressing out L_d. Error bars represent s.e.m. Feature weight similarity: P = 2.391 × 10⁻²¹; RF distance: P = 0.451; Benjamini–Hochberg corrected.

Fig. 4. Like-to-like effects are widespread but vary across brain areas, cortical layers and tuning similarity metrics.

a–f, Degree of like-to-like broken down by area and layer membership measured at axonal (a,c,e) and synaptic (b,d,f) scales. On the colour bar, red indicates more like-to-like and blue indicates less. Like-to-like coefficients are the coefficients of GLMMs fitted to predict axon–dendrite co-travel distance or synapse density with the corresponding functional similarity. Black borders indicate significance (P < 0.05); white border indicates P > 0.05; by GLMM fit after Benjamini–Hochberg correction for multiple comparisons. Details of statistical tests are presented in Supplementary Tables 15–20.

Fig. 5. Postsynaptic neurons with a common input are more functionally similar to each other than expected from a pairwise like-to-like rule.

a, Left, schematic illustrating the null hypothesis that postsynaptic neurons (postsyns, dark grey circles) of a common presynaptic neuron (presyn, yellow circle) have no additional feature similarity with each other beyond their like-to-like similarity with their common presynaptic neuron. Neurons are shown embedded in functional space, where units with similar functional properties are closer together. In this scenario, postsynaptic neurons are distributed uniformly around the presyn in the like-to-like region of functional space (large grey circle). Right, schematic illustrating the alternative hypothesis that the postsynaptic neurons are closer in functional space than predicted from a pairwise like-to-like rule, equivalent to being clustered non-uniformly within the like-to-like region. b, Schematic illustrating the functional connectivity model used to simulate the null hypothesis in a. Left, pairwise functional measurements—including signal correlations, feature weight similarity and RF location distance—were passed through a function relating functional similarity to connection probability. Right, then, within this modelled network, we computed the pairwise similarity of all postsynaptic neurons downstream of a common presynaptic neuron. c, Observed mean pairwise signal correlations between postsynaptic neurons compared with those expected from model predictions across projection types (n = 52 presynaptic neurons for HVA→HVA, n = 38 for HVA→V1, n = 17 for V1→HVA, n = 35 for V1→V1). Box plots show median, interquartile range (box) and 1.5× interquartile range (whiskers); points indicate outliers. Three of four projection types showed significantly higher similarities than predicted. Two-sided Wilcoxon signed-rank test; exact P values are presented in Supplementary Table 21.

Fig. 6. Like-to-like connectivity in an RNN.

a, A vanilla RNN with 1,000 hidden units was provided images as inputs and weights were trained so that a readout of the final state identified the input’s label. b, Mean signal correlations among all and connected neuron pairs for the same RNN before and after training. Neurons were classified as connected when their weights exceeded a fixed threshold. c, Connection probability as a function of signal correlation for the same network before and after training. d, Test accuracy of the network as a function of the number of connections that were ablated when ablating random or like-to-like connections. Connections were classified as like-to-like whenever the weight and signal correlation both exceeded a fixed threshold. e, Mean postsynaptic–postsynaptic signal correlations and the expected postsynaptic–postsynaptic signal correlation given a pairwise model similar to Fig. 5c before and after training (two-sided sign test: P = 9.85 × 10⁻⁵). The postsynaptic–postsynaptic signal correlations are computed for all 1,000 hidden units in the RNN. Box plots show median, interquartile range (box) and 1.5× interquartile range (whiskers); points indicate outliers.

Extended Data Fig. 1. Example proofread presynaptic axons in EM cortical space and their connected, ADP, and same region controls.

The axon for every presynaptic (presyn) neuron is shown twice, once as a “local” projection type and again as a “long-range” type (even if the neuron has no local or long-range projections). The six digit ID from Table “nucleus_detection_v0” is displayed above both plots. For each plot, the soma centroids of connected neurons, ADP controls, and same region controls are plotted in black, red, and blue, respectively. Gray dots are soma centroids of all other functionally matched neurons not used as controls for that presyn. The dashed gray line represents the V1-HVA boundary. Scale bar = 100 μm. a, Example fully proofread presynaptic axons with somas in V1. “Fully proofread” neurons are those where a proofreader attempted to extend every axonal branch to completion. b, Example fully proofread presynaptic axons with somas in HVA c, Example partially proofread presynaptic axons with somas in HVA. “Partially proofread” neurons are those where a proofreader only extended axonal branches that were pre-screened for whether they projected inter-areally (specifically to enrich for feedback connections).

Extended Data Fig. 2. The digital twin signal correlations align better with the in vivo benchmark than in vivo signal correlations generated with less data.

a, Correlation of in vivo signal correlations generated with 6 video clips and varying numbers of repeats to in vivo signal correlations generated with 6 clips and 30 repeats, for two animals. 10 repeats (red marker) reasonably approximates the saturation point and is the number used for all other analyses. b, Signal correlation matrices of 1000 neurons generated from in vivo responses to 6 video clips (left), in vivo responses to 30 video clips (benchmark, middle) and digital twin responses to 250 video clips (in silico, right). The benchmark matrix is ordered by Ward’s hierarchical clustering. The in vivo and in silico signal correlation matrices are ordered in the same order as the benchmark matrix. The fine structure of the in silico matrix is qualitatively more similar to the benchmark than the in vivo matrix generated with 6 video clips is to the benchmark. c, 2D heatmaps of signal correlations from the benchmark (same benchmark as in b) vs in vivo responses to 6 video clips (left) and in silico responses to 250 clips (right). The correlation of in silico signal correlations to the benchmark is higher than the correlation of in vivo signal correlations generated with 6 video clips to the benchmark (0.69 vs 0.40). Colorbar: 2D bin counts in log scale. d, The correlation of in silico signal correlations to the benchmark vs the correlation of in vivo signal correlations generated with 6 video clips to the benchmark for three animals. Error bars are standard deviations estimated through resampling. All data points are in the upper left corner indicating that in silico signal correlations outperform in vivo signal correlations generated with 6 video clips. (p-value < 0.001 for all three animals).

Extended Data Fig. 3. Synaptic connectivity increases with empirical signal correlations measured directly in vivo rather than via the digital twin.

a, Mean in vivo signal correlation is different (mean ± sem, paired t-test) for connected pairs, ADP controls, and same area controls for all projection types, as in Fig. 2d. b, Axon-dendrite co-travel distance (μmL_d) increases in a graded fashion with in vivo signal correlation for all projection types, as in Fig. 2e. c Synapse density (N_syn/mmL_d) increases in a graded fashion with signal correlation, for all projection types, as in Fig. 2f. The shaded regions in b and c are bootstrap-based standard deviation. d, Synapse size (log₁₀ cleft volume in voxels) is positively correlated with in vivo signal correlation after regressing out L_d (p-value by linear regression), as in Fig. 2h. e, In vivo signal correlations increases with number of synapses after regressing out L_d (p-values by linear regression), as in Fig. 2j. f, Area/layer joint membership breakout as in Fig. 4 for in vivo signal correlations at axonal scale. g, Area/layer joint membership breakout as in Fig. 4 for in vivo signal correlations at synaptic scale. h, Comparison of the observed and expected postsynaptic functional similarity as in Fig. 5 for in vivo signal correlations. (For all panels, * = p-value < 0.05, ** = p-value < 0.01, *** = p-value < 0.001, multiple comparison correction by BH procedure. For statistics and sample sizes, see Supplementary Tables 23–28).

Extended Data Fig. 4. Comparison of in silico and in vivo receptive field centers.

a, Visual comparison of Spike-Triggered Average receptive fields (STAs) generated from in vivo responses to a sparse noise stimulus (top row) vs STAs generated from in silico responses to the same stimulus (bottom row) for three animals (blue, orange, and green). The black cross represents the model readout location. Examples are randomly chosen from the top 44% of neurons remaining after a threshold on in vivo STA quality is applied. b, Model readout location vs in vivo STA center for azimuth coordinate (left) and elevation coordinate (right). c, in silico STA center vs in vivo STA center for azimuth coordinate (left) and elevation coordinate (right). d-i, Retinotopic maps for animal id: 29755. d, Retinotopic maps generated from in vivo STA centers with top 44% of neurons after an in vivo STA quality threshold is applied. Left: Azimuth retinotopic map, each dot represents one neuron in the cortical space, the color represents the azimuthal visual angle of its receptive field center. Middle: Elevation retinotopic map. Right: The coherence of the retinotopic map is visualized as a scatterplot of the pairwise cortical distance vs. the pairwise retinotopic distance for 10,000 randomly selected neuron pairs within the dotted circular region in the retinotopic maps. The coherence is quantified as the Spearman’s rank correlation coefficient between the distances. e, Retinotopic maps generated from in vivo STA centers with the bottom 25% neurons based on the quality of the in vivo STA. f, g, Retinotopic maps generated from in silico STA centers for the same neurons in d and e. h, i, Retinotopic maps generated from the digital twin model readout location for the same neurons in d and e. Colorbar: degree of visual angle for both azimuth and elevation coordinates. Anatomical axes: A = anterior, P = posterior, M = medial, L = lateral. Scale bar: 100 μm. j, k, Analysis in Fig. 3a,b repeated with in silico STA centers instead of model readout location.

Extended Data Fig. 5. Example of connected, ADP, and same area controls for the area/layer analysis in Fig. 4.

For all panels, a single presynaptic neuron skeleton is shown as in Fig. 2c with its postsyns (black dots), ADP controls (red dots), same area controls (blue dots), and all other functionally matched neurons (gray dots). The location of the presynaptic soma is indicated by the orange triangle. The vertical gray dotted line is the V1 and HVA boundary, with V1 on the left, and HVA on the right. The scale bars are 100 μm a, The V1 L2/3 neuron reproduced from Fig. 2c with its V1 → V1 postsyns and controls (left), and shown again with only its V1 L2/3 (left middle), V1 L4 (right middle) and V1 L5 (right) postsyns and controls. b, The V1 L2/3 neuron reproduced from Fig. 2c with its V1 → HVA postsyns and controls (left), and shown again with only its HVA L2/3 (left middle), HVA L4 (right middle) and HVA L5 (right) postsyns and controls. c, The HVA L5 neuron reproduced from Fig. 2c with its HVA → HVA postsyns and controls (left), and shown again with only its HVA L2/3 (left middle), HVA L4 (right middle) and HVA L5 (right) postsyns and controls. d, The HVA L5 neuron reproduced from Fig. 2c with its HVA → V1 postsyns and controls (left), and shown again with only its V1 L2/3 (left middle), V1 L4 (right middle) and V1 L5 (right) postsyns and controls.

Extended Data Fig. 6. Distribution of pairwise functional measurements.

Density distribution of connected pairs (black), ADP control pairs (red) and same region control pairs (blue) for in vivo signal correlations (a), in silico signal correlations (b), feature weight similarity (c), and RF center distance (d) for all projection types.

Extended Data Fig. 7. Pairwise functional measurements across varying levels of model predictive performance.

Mean of in vivo signal correlations (a), in silico signal correlations (b), feature weight similarity (c), and RF center distance (d) for all projection types across 4 quantiles of model predictive performance (CC_abs). All panels share a base filtering for visual responsiveness (CC_max > 0.4, 90% of neurons pass this threshold). Presynaptic neurons are filtered to CC_abs > 0.2 (4 did not pass this threshold).

Extended Data Fig. 8. In silico orientation tuning is consistent with in vivo orientation tuning.

a, Sample frame from global directional parametric stimulus (“Monet”) used to characterize orientation and direction selectivity. Directional motion was orthogonal to orientation, and was tested at 22.5° intervals. b, Schematic of domain validation experimental design. In a single scan in a new animal, neuronal responses are collected in response to sufficient stimuli to both train the digital twin model (natural stimuli) and characterize orientation tuning (“Monet”) from in vivo responses. Later, in silico orientation tuning is extracted from model responses to parametric stimuli, and compared against in vivo orientation tuning for the same neurons. c, Comparison of in silico (red dotted line) and in vivo (black solid line) mean responses per stimulus direction, fitted tuning curves, and extracted preferred orientation (vertical lines) for three neurons of various gOSI levels. d, 95th percentile difference in preferred orientation between in silico and in vivo fitted responses as a function of gOSI threshold. Dotted lines correspond to gOSI > 0.25 threshold applied for all analyses and resulting 95th percentile difference in preferred orientation ≈19.7° across all three animals imaged. Lines correspond to individual animals (gray) or cumulative across all animals (black). e, f, Two-dimensional histogram of in silico versus in vivo preferred orientation for all neurons across three animals (e) and only neurons with gOSI > 0.25 (f).

Extended Data Fig. 9. Analysis in Fig. 3 repeated with in silico orientation preference.

a, Difference in preferred orientation (Δ Ori) derived from in silico responses to parametric stimuli for tuned (gOSI > 0.25) neurons along with both feature weight similarity and receptive field center distance (reproduced from Fig. 3) at axonal scale. b, same as in a, at synaptic scale. For analysis with synapse size, see Extended Data Fig. 13) c, Area/layer joint membership breakout as in Fig. 4 for in silico Δ ori at axonal scale. d, As in c but at synaptic scale. All analyses are centered per presyn by accounting for the presyn mean (e.g. Δ feature weight similarity). For statistics and sample sizes, see Supplementary Tables 7–14, 29–34.

Extended Data Fig. 10. Distribution of in silico orientation preference and comparison to previous literature.

a, Distribution of orientation preference of tuned neurons (gOSI > 0.25) derived from in silico responses to parametric stimuli (see Methods). Note the cardinal bias in orientation preference distribution, in which orientation preference for 0 and 90 degree angles is overrepresented. Gold: presynaptic neurons, Gray: all other neurons. b, As in a but for tuned neurons in V1 L2/3. Difference in preferred orientation (Δ Orientation) for neurons in V1 L2/3 for connected pairs (c, f), unconnected pairs (d, g), and the ratio of connected / unconnected (“connection probability”, e, h) for our study vs Lee et al. (c-e) and vs Ko et al. (f-h). The connected V1 L2/3 neurons in our study show a strong like-to-like effect, consistent with both Lee et al. and Ko et al. (c, f), however unlike Lee et al. and Ko et al., the unconnected neurons in our study also show a strong like-to-like effect (d, g) indicating that the like-to-like effect seen in connected pairs results from an orientation preference bias. This bias likely explains why we do not observe significant a like-to-like effect between V1 L2/3 neurons at axonal scale or synaptic scale in Extended Data Fig. 9, (i.e. when pairs are tested against region-matched controls).

Extended Data Fig. 11. Performance of various functional metrics in predicting axon-dendrite co-travel distance (L_d, Axonal scale) or synapse density (N_syn/mmL_d, Synaptic scale).

Model performance of GLMMs (Nakagawa’s conditional R²) for predicting axon-dendrite co-travel distance (L_d): a, b, c and synapse density (N_syn/mmL_d): d, e, f, for all coregistered neurons: a, d, all visually responsive, well predicted neurons: b, e, and neurons tuned to oriented stimuli: c, f. The GLMMs are fit to predict axon-dendrite co-travel distance or synapse density independently with each functional metric, the projection type, and the interaction between the two while considering the interaction term of projection type and presynaptic neuron identity as random effects. The baseline models were not fitted with information about functional metrics. They predict axon-dendrite co-travel distance or synapse density with the projection type alone while considering the interaction term of projection type and presynaptic neuron identity as random effects.

Extended Data Fig. 12. Pairwise functional measurements for connected and unconnected pairs vs soma-soma distance for V1 → V1 connections.

Mean ± SEM of in vivo signal correlations (a), in silico signal correlations (b), feature weight similarity (c), and RF center distance (d) vs soma-soma distance for connected and unconnected pairs for V1 → V1 connections.

Extended Data Fig. 13. In silico Δ Ori, RF location similarity, and feature weight similarity vs synapse size density.

a, Analysis in Extended Data Fig. 9a,b repeated with synapse size density, rather than synapse density. Synapse size density is computed similarly to synapse density except that the numerator N_syn is replaced with the summed synaptic weight (sum of all synapse cleft volumes for a pair of neurons in 4 × 4 × 40 nm³ voxels).

Extended Data Fig. 14. Postsynaptic neurons with a common input are more similar to each other than expected by a pairwise like-to-like rule at both axonal and synaptic scale.

a, Mean pre-post signal correlations in the data (dark gray, “observed”) and the model (blue, “expected”) are not significantly different, indicating that the model reproduces the expected pairwise like-to-like rule b, Mean pairwise in silico signal correlation of postsyns, reproduced from Fig. 5c. The observed data shows significantly higher postsyn to postsyn similarity than predicted by the model fit with only a pairwise rule, for three out of four projection types. c, As in a, but at “Axonal” scale. d, As in b, but at “Axonal” scale. e, As in c, but at “Synaptic” scale. f, As in d, but at “Synaptic” scale. For statistics and sample sizes, see Supplementary Tables 21, 22.

Extended Data Fig. 15. Signal correlation distributions for connected neurons vs all neurons in the RNN before and after training.

a, Signal correlation distribution for connected neurons vs all neurons in the RNN before training. A neuron pair was classified as connected if the associated weight was in the top 35^th percentile of all weights. b, Same as a except after training.