Introduction

Neurons have characteristic preferences, or tuning, and their activity is transmitted via axonal projections to other brain areas. In the early stages of visual processing, visual information can be preserved. For example, the retina-to-lateral geniculate nucleus (LGN) network tends to preserve unmixed channels, by ensuring that axons from retinal ganglion cells with similar tuning tend to converge on individual LGN neurons (Liang et al., 2018). By contrast, the LGN-to-primary visual cortex (V1) network famously mixes channels to give neurons receptive fields with both dark-sensing and light-sensing subregions, and thus robust orientation tuning (Hubel and Wiesel, 1962). That said, discrete visual information can be transmitted from the retina to cortex. For example, when direction-selective neurons in the retina are genetically ablated, there is a decrease in direction-selective neurons (Rasmussen et al., 2020).

In the visual cortical system in mice, the primary visual cortex (V1) and its projections to multiple higher visual areas (HVAs) span millimeters (Wang and Burkhalter, 2007). Local networks within V1 can have precise local (< 50 microns) cellular-resolution functional connectivity (Ko et al., 2011). Studies of longer-range, millimeter-scale networks typically lack cellular resolution, but there are general biases observed. Neurons in V1 and HVAs respond to diverse visual stimuli (Yu et al., 2022; Vries et al., 2020) and are sensitive to a broad range of features including orientation and spatiotemporal frequencies (Andermann et al., 2011; Marshel et al., 2011). Although individual V1 neurons broadcast axonal projections to multiple HVAs (Han et al., 2018), the spatiotemporal frequency preferences of these feedforward projections generally match those of the target HVAs (Glickfeld et al., 2013; Han and Bonin, 2023). Feedback connections from HVAs carry frequency-tuned visual signals as well (Huh et al., 2018). Thus, there are HVA-specific spatiotemporal biases, but cellular resolution, millimeter-scale principles for cortical wiring remain to be elucidated.

In the current study, we investigated the degree of discreteness of the V1-HVA functional networks in direction and frequency information processing, by measuring the noise correlations (NC, also called spike count correlation Vinci et al. (2016)) between functional tuning groups of neurons. Functional tuning groups were defined using an unbiased clustering approach (Han et al., 2022; Yu et al., 2022; Baden et al., 2016). Large field-of-view (FOV) calcium imaging enables us to densely sample across millimeter cortical space, which allows us to simultaneously observe large and dense samples of neurons in these tuning groups to analyze NCs within- and-across cortical areas (Yu et al., 2021; Stirman et al., 2016). NCs are due to connectivity (direct or indirect connectivity between the neurons, and/or shared input), and thus provide a trace of connectivity (Cohen and Kohn, 2011; Vinci et al., 2016; Snyder et al., 2015). In particular, the connectivity that influences NCs is effective in vivo, during normal sensory processing, something that is not verified with purely anatomical measures. Activity-based estimates of neuronal networks can provide higher fidelity measures than anatomical-based studies (Randi et al., 2023). We find that NCs are a reliable measure at the population level. We also find that neurons can be categorized into six functional groups, and NCs are higher within these groups, both within and across cortical areas, indicating discrete channels in the network to preserve information. Moreover, we find that naturalistic videos draw upon the same functional networks, and modeling suggests that recurrent connectivity is critical for stabilizing these networks.

Results

Visual cortical neurons form six tuning groups

To measure neuronal activity, we used multi-region population calcium imaging of L2/3 neurons in V1 and four HVAs (lateromedial, LM; laterointermediate, LI; anterolateral, AL; and posteromedial, PM) using a multiplexing, large field-of-view two-photon microscope with sub-cellular resolution developed in-house (Stirman et al., 2016) (Figure 1A). Mice expressed the genetically encoded calcium indicator GCaMP6s (Madisen et al., 2015; Chen et al., 2013) in cortical neurons. We located the V1 and HVAs of each mouse using retinotopic maps obtained by intrinsic signal optical imaging (Marshel et al., 2011; Smith et al., 2017) (Figure S1A). We imaged neurons in two to four cortical areas simultaneously (Figure 1A), while mice viewed stimuli on a display system. We typically imaged neurons in V1 and one or more HVAs. Up to 500 neurons (V1: 129 ± 92; HVAs: 94 ± 72; mean ± SD.) were recorded per imaging region (500 × 500 μm²). The imaging regions were matched for retinotopy so that the neurons in the simultaneously imaged areas had overlapping receptive fields (RFs). Calcium signals were used to infer probable spike trains for each neuron, as our previous study (Yu et al., 2022). We mapped RFs for individual neurons and populations using small patches of drifting gratings (Figure S1B, C). Neurons in HVAs (LM, AL, PM and LI) had significantly larger RFs than V1 neurons (Figure S1D). Population RFs for a 500 × 500 μm² imaging region of HVAs covered significantly larger portions of visual space than that of V1 (Figure S1D), as expected given their differing magnification factors (Schulz et al., 2015). The overlap of population RFs confirmed that simultaneously imaged cortical areas (V1 and HVAs), each containing ∼102 neurons, responded to stimuli in the same region of the stimulus display system (Figure S1C). These experiments were repeated in 24 mice for a total of 17,990 neurons and noise correlations were measured for a total of 1,037,701 neuron pairs (Table 1).

Figure 1.

Mouse V1 and HVA neurons exhibit diverse tuning preferences (i.e., biases) to drifting grating stimuli, in terms of spatiotemporal preferences and sharpness of orientation and direction tuning (Marshel et al., 2011; Andermann et al., 2011; Vries et al., 2020). Previous studies suggested that the projection from V1 to HVAs matches the spatiotemporal preferences of the target HVAs (Glickfeld et al., 2013). We sought to determine whether this was a general principle, that extended across V1 and HVAs. We recorded neuronal responses from V1 and multiple HVAs (LM, LI, AL, and PM) to sine-wave drifting grating stimuli with various spatiotemporal properties (8 directions x 3 spatial frequencies x 3 temporal frequencies for a total of 72 conditions; Figure 1B). HVAs exhibited similar responsiveness and reliability to the 72 different parameterized drifting gratings. V1 and LM were only marginally more reliable than other areas (Figure S1E). Neurons were partitioned into 65 tuning classes using an unbiased Gaussian Mixture Model (GMM) (Figure S1F, S2). This GMM classification was reliable, in that the center of the Gaussian profile of each class was consistent among GMMs of random subsets of neurons (Figure 1C). Neurons were consistently classified into the same class (Methods; Figure 1D).

To examine the spatiotemporal frequency selectivity of HVAs, we manually partitioned the 65 GMM classes into six spatial frequency (SF) - temporal frequency (TF) selective groups (Figure 1E). Groups 1, 2, and 3 all prefer low TF (1-2 Hz), and prefer low SF (0.02 cpd), medium SF (0.05 cpd), and high SF (0.19 cpd) respectively. Groups 4, 5, and 6 all prefer high TF (8 Hz) and prefer low SF, medium SF, and high SF respectively. Group 4 (low SF, high TF) was the only group that exhibited clearly increasing responses to the drift speed of the grating stimulus (drift speed = TF/SF, and is measured in deg/s). These groupings were robust and reliable (Figure S1G,H). While all visual areas had a similar fraction of neurons tuned to low TF and low SF (Group 1), fractions in other groups varied by area. AL had a larger fraction of neurons tuned to high TF, and low SF (Groups 4) (Figure 1E). PM and LI had a larger fraction of neurons tuned to high SF and low TF (Groups 3) (Figure 1E). Overall, there is a trend of increasing of preferred TF from the posterior-medial to the anterior-lateral visual cortex, and a trend of increasing of preferred SF from the anterior to the posterior visual cortex (Figure S3A, B).

Neurons in all six groups exhibited orientation and direction selectivity. The preferred directions of neurons were evenly distributed in V1 and HVAs, except high SF groups (Group 3 and 6) of AL, PM, and LI biased to cardinal directions (Figure S3C). The unbiased GMM approach revealed that the orientation selectivity index (OSI) and direction selectivity index (DSI) of visual neurons were jointly modulated by SF and TF. Neurons tuned to high SF and low TF (Group 3) exhibited lower OSI in all tested areas than all of the other groups (Group 3: mean OSI = 0.6; other groups ranged from 0.71 – 0.80; p < 0.0001, one-way ANOVA with Bonferroni correction; Figure 1E). Neurons tuned to high TF and medium-high SF (Groups 5 and 6) exhibited lower direction selectivity than other groups (Group 5, 6, mean DSI 0.38; other groups, mean DSI 0.45 - 0.54; p < 0.0001, one-way ANOVA with Bonferroni correction; Figure 1E).

NCs are robust measurements of functional networks

A unique aspect of this data set is the scale of the NC measurements, which allows us to measure NCs with individual neuron precision within dense local networks and over millimeter length scales, in awake mice. Pioneering work in this area focused on local populations, typically less than 1 mm across (Ko et al., 2011; Harris and Mrsic-Flogel, 2013; Lee et al., 2016; Wertz et al., 2015) or electrode studies over long distances with few neurons in each location (Clay Reid and Alonso, 1995; Siegle et al., 2021). To investigate the V1-HVA functional network, we computed the NCs of pairs of neurons within individual cortical areas (within-area NC), and NCs for pairs of neurons where the two neurons are in different cortical areas (inter-area NC) (Figure 2A). NCs are computed from the residual activity of individual neurons after subtracting the expected neuron firing on nominally identical trials. In this section, we evaluated the fidelity of our NC measurements. We considered potential measurement noise due to the imprecision of spike inference and the finite number of trials.

Figure 2.

We first evaluated the accuracy of NC calculations using inferred spikes from calcium imaging. We characterized the accuracy of spike inference using previously published data of simultaneous two-photon imaging and electrophysiological recording of GCaMP6s-positive neurons from mouse V1 (Chen et al., 2013). Consistent with a previous benchmark study on spike train inference accuracy (Theis et al., 2016), we found that the spike train inference methods used in the current study recovered 40-70% of the ground truth spikes (Figure S4A). We found that a similar fraction of spikes were missing regardless of the inter-spike interval. Nevertheless, the inferred spike train was highly correlated with the true spike train (Figure S4A; linear correlation, r = 0.80 ± 0.03 (n = 6)). Computing correlations between pairs of neurons using their inferred spike trains accurately reproduced the true correlation values (Figure S4B; linear correlation, r = 0.7). We further examined the fidelity of correlation calculations using modified spike trains that are missing spikes. We examined randomly deleting spikes, deleting isolated spikes, or deleting spikes within bursts (Figure S4D; Methods). We found that at the 1 s time scale, correlation calculations were tolerant to these spike train perturbations. The fidelity of correlation computations were >0.6 with up to 60% missing spikes (Figure S4E; Methods). Thus, with conventional spike inference accuracy, about 80% variance of the true correlation is recovered (Figure S4E). Thus, NCs are a robust measurement even with imperfectly inferred spike trains.

Next, we evaluated the robustness of NC measurements, given the finite number of trials that are feasible to obtain. We computed NCs for both within-area neuron pairs and inter-area neuron pairs (Figure 2A). NCs were computed using spike counts within 1 s bins, similar to previous work with electrophysiology (Cohen and Kohn, 2011; Smith and Sommer, 2013). Although both within-area and inter-area NCs had wide distributions (range: −0.2 – 0.6). The mean NCs across a population were positive and at least five times larger than control data, which are NCs computed after shuffling the trials (5 – 20-fold, 25 – 75% quantile; Figure 2B). The estimation of the population mean NC converges fast with increasing numbers of neurons, as suggested by both simulation and experimental data (e.g., the margin of error at 95% CI for mean NC is 0.008 for 100 neuron pairs, Figure 2C). While the population-level NC calculations are reliable, the NC estimation of individual neuron pairs is noisier due to the limited number of trials, albeit positively correlated (Figure 2D). A linear model explains about 53 ± 24 % of the variance between NCs for individual neuron pairs computed using different random subsets of trials. In summary, this evidence indicates that NCs can be accurately measured at the population level with our large FOV calcium imaging methods, despite imperfect spike train inference and a finite number of trials.

Tuning similarity is a major factor in the V1-HVA functional network

Having established that NC measurements were reliable and robust, we examined potential NC-regulating factors, including firing rate, the physical distance between the neurons (laterally, across the cortex), signal correlation (SC, the similarity between two neurons’ average responses to stimuli), and RF overlap. We assessed the contributions of individual factors using a linear models. We found that both within- and inter-area NCs are similarly modulated by the aforementioned factors (Figure 3A; r² of linear regression model). SC is the most pronounced factor that explains about 10% of the variance of within-area NCs, and about 5% of the variance of inter-area NCs (Figure S5A). Fraction RF overlap contributes about 6% and 3% variance of within- and inter-area NCs, respectively (Figure S5B). The firing rate explained about 2% of the variance of within- and inter-area NCs (Figure S5C). Thus SC, RF overlap, and firing rate positively regulate both within and inter-area NC, with SC providing the strongest predictor.

Figure 3.

We then evaluated whether NCs are modulated by tuning similarity independent of RF overlap. In the subset of orientation-selective neurons, both within- and inter-area NCs were significantly modulated by orientation tuning selectivity. That is, neuron pairs that shared the same preferred orientation exhibited higher NCs (Figure S5D). NCs of a subset of neurons with non-overlapping RFs were significantly higher when the neurons shared the same preferred orientation (Figure 3B; t-test, p < 0.05 for V1, LM and AL neuron pairs, insufficient data for PM and LI). This result confirms that the connectivity between neurons is modulated by tuning similarity (SC) independent of RF overlap, over millimeter distance scales.

Overall, about 20% of the variance of within-area NCs, and 10% of the variance of inter-area NCs are explained by the aforementioned factors jointly (Figure 3C; r² of multi-linear regression model). Although inter-area NCs have a smaller mean and variance, it is less predictable by known factors (within-area NCs pooled over all tested area 0.012 ± 0.052, inter-area NCs between V1 and all tested HVAs 0.0063 ± 0.04; both t-test and F-test p < 10⁻⁴). In an expansion of prior work on local functional sub-networks (Lee et al., 2016; Wertz et al., 2015; Ko et al., 2011; Harris and Mrsic-Flogel, 2013), we find that signal correlation is the strongest factor regulating both within-area and inter-area NC networks, suggesting that neurons exhibiting similar tuning properties are more likely to form functional sub-networks across a broad spatial scale, spanning millimeters in the mouse V1-HVA network.

Neurons are connected through functionally discrete channels

Seeing as HVAs exhibited biased SF-TF selectivity (Figure 1E), and after finding that tuning selectivity (SC) is a major factor for functional connectivity even across the millimeter length scale (Figure 3A), we assessed the precision of this network. We performed additional analysis to determine whether the ST-TF biases in HVAs could be due to simple, weak biases in the NC network. Alternatively, there could be precise discrete channels of connectivity in the V1-HVA network to preserve information among similarly tuned neurons. We found evidence for this latter situation: discrete channels. Moreover, we found that the discrete channels consist of a greater number of neuron pairs with high NCs, rather than a small number of neuron pairs with very high NCs.

For this analysis, we focused on neuron pairs with high NCs, which we defined as NCs > 2.5 times the standard deviation of trial shuffled NCs for the population (Figure 4A). We focused on these high NC pairs, because they can represent high fidelity communication channels between neurons. Within each SF-TF group, for both V1 and HVAs, about 10-20% of neuron pairs exhibited high NCs, in contrast to 5% for inter SF-TF group connections (Figure 4A, B). The fraction of pairs that exhibit high NCs is relatively uniform across tuning groups and HVAs with a few exceptions (Figure 4B). For example, in HVA PM group 3 contains a higher fraction of high fidelity connections than the other HVAs. Overall, these results show that mixing between groups is limited, and instead group-specific high-NC sub-networks exist between neurons across millimeters cortical space.

Figure 4.

To determine whether it is the number of high NC pairs or the magnitude of the NCs, we performed further analysis. Prior findings from studies of axonal projections from V1 to HVAs indicated that the number of SF-TF-specific boutons —rather than the strength of boutons— contribute to the SF-TF biases among HVAs (Glickfeld et al., 2013). Though the functional connectivity is not completely defined by the feedforward axonal projections from V1 (Huh et al., 2018), this number vs. strength question is one that we can address with our data set. We found that the biased representation of SF-TF among HVAs is linearly related to the number of neuron pairs with high NCs (Figure 4C). That is, there are more high NC pairs that are in-group. By contrast, the fraction of pairs with high NCs does not account the SF-TF biases of HVAs (Figure 4D). That is, the in-group pairs do not tend to have higher NCs. Thus, the biases in SF-TF are likely related to the abundance of SF-TF-specific connections, but not the strength of the connections.

To this point, we have focused on the six SF-TF groups. The evidence supports group-specific channels among these neurons. However, these six groups originated with 65 classes from data-driven GMM clustering, which were then manually collected into the six SF-TF groups (Figure 1). The trends we see for groups may reflect general SF-TF biases. In that case, we would expect that the in-class NCs would exhibit similar distributions of NCs as the in-group NCs. However, there might be further precision in the specific channels not captured by the SF-TF groups. A hint towards that can be seen in the fact that orientation tuning can modulate NCs (Figure 3), because some of the 65 GMM classes differ by preferred orientation (Figure S2). Indeed, we found that when we plotted the NC distribution for in-class neuron pairs and compared it to the distribution for in-group neuron pairs, we found a pronounced positive tail for the in-class distribution (Figure 4E). Thus the GMM classes provide relevant, granular labels for neurons, which form functional sub-networks with discrete channels, which are more precise than predicted from simpler SF-TF biases or groups.

The GMM classes are widely distributed in all tested areas (Figure S7B). We constructed an inter V1-HVA connectivity matrix for the 65 classes (Figure 4F). The connection weight is defined by the fraction of pairs with high NCs. To investigate the modular structure of this network, we performed community detection analysis using the Louvain algorithm (Rubinov and Sporns, 2010). This analysis assigned densely connected nodes to the same module (Figure 4F). Overall, the connectivity matrix was split into four community modules (Figure 4F; S7C). Interestingly, the corresponding node in V1 and HVAs for each module exhibited similar direction and SF-TF preferences (Figure 4F). For example, the module 2 nodes exhibited narrow vertical direction tuning and preferred high SF and low TF. Module 1 exhibited high SF preference without direction bias. Area differences in the characteristic tuning selectivity of each module are small, suggesting that the GMM class channels are common across the V1HVA network. This is consistent with the overall broad-casting projection structure of V1 neurons (Han et al., 2018).

In summary, V1 and HVA neurons can be classified by their selectivity to oriented gratings, and they form precise, discrete channels or sub-networks. These subnetworks of neuron pairs with high NCs preserve selectivity by limiting inter-channel mixing. The organization of V1-HVA sub-networks exhibited properties that are consistent with that of V1-HVA feedforward projections in that the number of high-fidelity connections, rather than the strength of the connections, accounted for SF-TF biases among HVAs. Moreover, the precision of these networks extends beyond prior observations of general SF-TF biases, to include orientation and direction tuning.

Functional connectivity is stable across stimuli

Functional connectivity is dynamic and transient, which complicates its relationship with structural (i.e., anatomical) connectivity, yet can provide more accurate predictions for network dynamics than the latter (Randi et al., 2023). We performed additional analysis to determine whether the NC-based functional connectivity analysis we performed above provides fundamental insights into neuron circuits beyond a stimulus-specific transient. We compared NC measurements in response to drifting gratings (NC_grat) to NC measurements in response to naturalistic videos (NC_nat). This analysis was restricted to the subset of neurons that responded to both types of stimuli in a separate set of experiments.

So far, we have shown that SC (i.e., neuron tuning similarity) is the best predictor for NCs. However, a neuron pair that shares a high SC to drifting gratings does not guarantee a high SC to naturalistic videos (corr(SC_grat, SC_nat) = 0.084 ± 0.065). Thus, it is reasonable to expect that NCs in response to gratings do not predict the NCs in response to naturalistic videos. However, we were surprised to find that the correlation between NC to the two stimuli is significantly higher than that of SC (corr(NC_grat, NC_nat) = 0.22 ± 0.13; Figure 5A). Thus, NC across stimuli is more predictable than SC across stimulus types. To our knowledge, this is the first time this has been reported.

Figure 5.

We used SC to natural videos (SC_nat), and gratings NC_grat to predict NC_nat using linear regressions. Both predictors are positively related to the NC_nat (Figure 5B). We found that NC to gratings outperformed SC to naturalistic videos in predicting NC to naturalistic videos (t-test, p < 0.0001; Figure 5C). Meanwhile, combining both predictors almost linearly adds up the prediction power of the two factors in linearly predicting NC to natural videos (Figure 5C), suggesting that the cross-stimulus NC predictor adds an independent dimension to the SC predictor. These results are evidence that NC-assessed functional connectivity reflects a fundamental aspect of the architecture of neuronal circuitry that is independent of visual input.

Recurrent connection contributes to the stability of NC network

Surprised by the cross-stimulus stability of the NC-based functional connectivity of the visual cortical network among V1 and HVAs, we went on to investigate what could be an underlying mechanism. NCs can be due to both shared input and direct/indirection wiring. Indeed, using a simple model with two leaky integrate-and-fire (LIF) neurons, we found that the NC is positively regulated by a larger fraction of shared input as well as by the increasing recurrent connection strength (Figure 6A, B).

Figure 6.

We then asked how the two sources contribute to the cross-stimulus stability of the NC functional network using LIF neuronal network simulations (Figure 6C). The simulated neuronal network contained 80 excitatory neurons and 20 inhibitory neurons that are randomly connected. The input layer contains 1000 independent Poisson spiking neurons. The network parameters are determined based on previous work (Song et al., 2000) and all the simulations generated comparable LIF firing rates (4-6 Hz), as well as NCs (population mean: 0.05−0.25) and SC values (population mean: 0.01−0.15).

In the first set of simulations, the feedforward (FFD) connection from the input layer to the LIF network is random. Increasing recurrent connection strengths (ranging from 0.05 to 0.3) generated NC-based networks with higher cross-stimulus stability (Figure 6D). A recurrent connection strength of 0.2 best reproduced the mouse data. In the second batch of simulations, we fixed the recurrent connection strength to 0.2 but manipulated the input FFD connection structure ranging from random to increasingly wider bell shapes (Figure 6E). This means that the local neurons receive increasingly similar FFD input. We found that increasingly similar local FFD input does not lead to higher NC stability, but does increase SC similarity across stimuli (Figure 6E). Also, the random FFD input connection structure (0.18 FFD, red) reproduced the experimentally observed NC network the best (Figure 6E).

Thus, the LIF simulations showed that although both shared input and recurrent connections contributed to the NC, the recurrent connections are critical for generating the observed cross-stimulus stability of the NC functional network. Moreover, the simulations suggested that the network can receive relatively random FFD input.

Discussion

We used large scale two-photon calcium imaging across cortical areas to show that NC-based assessments of functional connectivity exhibited tuning-specific organization with single neuron precision, across millimeter length scales. This connectivity is consistent with the axonal projection patterns obeserved in the mouse visual cortex (Glickfeld et al., 2013; Han et al., 2018), but goes further, to include the precise GMM classes we found– 65 different classes in response to grating stimuli. Moreover, we found that NC-based networks are consistent across stimulus classes. Thus, V1 broadcasts high-fidelity channels of information to HVAs. The projections preserve fidelity by minimizing mixing among channels, i.e., tuning groups or classes.

Multi-scale functional connectivity and structural connectivity

To better understand this aspect of network organization, we examined the connectivity of the nervous system at multiple scales (Sporns, 2016), and considered the coefficient of variance (CV, = standard deviation / mean) at each level and process (Figure 7). At the small scale, neurons interact through the release of neurotransmitters, which can affect ion channel activity. At the mesoscale, e.g. the current study, neurons are organized into networks or circuits with specific structures and hierarchies to serve a function in information processing, processing stimuli and/or guiding behavior. At the macroscale, brain regions can have diverse functions and interact to generate adaptive behavior. Assessing neuronal dynamics at these different scales requires a variety of technologies (Dorkenwald et al., 2022; Harris et al., 2019; Silver et al., 2003; Stringer et al., 2019).

Figure 7.

We compiled data on multiple scales and modalities of mechanistic connectivity of the mouse nervous system. To compare across modalities and scales we focused on the CV, called dispersion hereafter, of various parameters (Figure 7). Dispersion (CV) is unitless and allows us to bring disparate data into the same axis. We compared the dispersion of multi-scale connectivity to gain insights into the computations bridging multi-scale interactions. Notice that both functional and structural connectivity are available at the microscopic scale, but only one of the two is available for mesoscale and macroscale networks (Figure 7). Functional connectivity has been measured at the macroscale in mice using one-photon techniques (Ren and Komiyama, 2021) and fMRI (Kim et al., 2023), but these approaches lack cellular resolution, and the analyses to date do not provide a quantitative measurement of dispersion that can map onto the same axis, thus they are omitted from this discussion. Still, we can provide a large range of dispersion measurements (Figure 7).

The dispersion of components of multi-scale connectivity of the nervous system ranges from 0.1 to 1.5 (Figure 7). Synaptic interactions are characterized by the number and size of synapses, and the post-synaptic potential. The size of synapses and the post-synaptic potential exhibited larger dispersion compared to the number of synapses, suggesting that neuronal interactions are regulated by the size and strength of synapses rather than the number of synapses (Figure 7). The combination of high-resolution electron microscopy (EM) and slice electrophysiology has un-covered a linear relationship between synaptic size and strength in the mouse L2/3 somatosensory cortex (Holler et al., 2021).

Overall, dispersion reduces at larger spatial scales (Figure 7). The measurement at these larger scales could mask some key sources of variance in connectivity. This becomes clear when comparing the dispersion of the tuning channel-specific connectivity with that of the mixing-channel connectivity measured by the current study. The dispersion of the latter reduces for both inter- and within-area connectivity, suggesting that tuning channel-specific communication is a key contributor to defining the functional connectivity between brain regions (Figure 7).

Similarly, axonal projections measured without tuning channel-specific information provide the connectivity between all cortical regions and shed light on laminar differences, but do not discriminate brain areas thoroughly (Figure 7). In particular, the dispersion of V1-HVA L2/3 projections is small and similar to that of the L2/3 channel-mixing NC functional network, emphasizing the importance of tuning specificity in organizing the network in these cortical layers. L1 being the major destination for cortical feedback projections (Harris and Shepherd, 2015), exhibited the largest dispersion in the anterograde projection connectivity.

Additional efforts are needed to bridge analysis of multi-scale connectivity. For example, models have been developed to understand how small-scale processes, such as synaptic interactions, collectively give rise to population neuron-level behaviors (Gerstner and Kistler, 2002). Synaptic interactions are described at various levels of detail and biological realism in these models to serve specific research purposes. Population-level neuronal network models can be scaled up to represent brain regions or systems, which involves adding more neurons, specifying cell types, and connecting multiple networks (Markram et al., 2015; Sporns, 2014). Integrating functional connectivity with anatomical connectivity would help develop a system model with functional causality (Honey et al., 2010; Sporns, 2016). These models can simulate complex tasks or behaviors by connecting various brain regions and networks, accounting for sensory inputs, decision-making, and motor outputs.

In summary, mesoscale two-photon imaging techniques open up the window of cellular-resolution functional connectivity at the system level. How to make use of the knowledge of functional connectivity remains unclear, given that functional connectivity provides important constraints on population neuron behavior. One approach to address the gap is to link functional connectivity with structural connectivity to learn how the two constrain each other and integrate the functional connectivity into network models.

Acknowledgements

Funding was provided by grants from the NIH (R01EY024294, R01NS091335), the NSF (1707287, 1450824), the Simons Foundation (SCGB325407), and the McKnight Foundation to SLS; a Helen Lyng White Fellowship to YY; a career award from Burroughs Wellcome to JNS; and training grant support for CRD (T32NS007431).

Author Contributions

All experiments and analyses were performed by YY. The imaging system was built by JNS. Analysis was assisted by CRD. Study design and supervision by SLS. YY and SLS wrote the paper.

Disclosures

SLS serves as a consultant for optics and neuroscience companies and is a founder of Pacific Optica.

Methods

Animals and surgery

All animal procedures and experiments were approved by the Institutional Animal Care and Use Committee of the University of North Carolina at Chapel Hill or the University of California Santa Barbara and carried out in accordance with the regulations of the US Department of Health and Human Services. GCaMP6s expressing transgenic adult mice of both sexes were used in this study. Mice were 110 - 300 days old for data collection. GCaMP6s expressing were induced by the triple crossing of TITL-GCaMP6s line (Allen Institute Ai94), Emx1-Cre line (Jackson Labs #005628), and ROSA:LNL:tTA line (Jackson Labs #011008)(Madisen et al., 2015). Mice were housed under a 12-h/12-h light-dark cycle, and experiments were performed during the dark cycle of mice. Mice were anesthetized with isoflurane (1.5 - 1.8%) and acepromazine (1.5 - 1.8 mg/kg body weight) when performing visual cortex craniotomy. Carpofen (5 mg/kg body weight) was administered prior to surgery. Mice’s body temperature was maintained using physically activated heat packs during surgery. Mouse eyes were kept moist with ointment during surgery. The scalp overlaying the right visual cortex was removed, and a custom head-fixing imaging chamber with a 5-mm diameter opening was mounted to the skull with cyanoacrylate-based glue (Oasis Medical) and dental acrylic (Lang Dental). A 4-mm diameter craniotomy was performed over the visual cortex and covered with #1 thickness coverslip.

Locating visual areas with intrinsic signal optical imaging (ISOI)

ISOI experiments were carried out similarly as previously (Stirman et al., 2016; Smith et al., 2017; Smith and Trachtenberg, 2007). Briefly, the pial vasculature images and intrinsic signal images were collected using CCD camera (Teledyne DALSA 1M30) at the craniotomy window. A 4.7 × 4.7 mm² cortical area was imaged at 9.2 μm/pixel spatial resolution and at 30 Hz frame rate. The pial vasculature was illuminated and captured through green filters (550 ± 50 nm and 560 ± 5 nm, Edmund Optics). The ISO image was collected by focusing 600 μm down from the pial surface. The intrinsic signals were illuminated and captured through red filters (700 ± 38 nm, Chroma and 700 ± 5 nm, Edmund Optics). Custom ISOI instruments were adapted from Kalatsky and Stryker (2003). Custom acquisition software for ISOI imaging collection was adapted from David Ferster (Stirman et al., 2016). During ISOI, mice were 20 cm from a flat monitor (60 × 34 cm²), which covers the visual field (110°x 75°) of the left eye. Mice were lightly anesthetized with isoflurane (0.5%) and acepromazine (1.5–3 mg/kg). The body temperature was maintained at 37 °C using a custom electric heat pad (Stirman et al., 2016). Intrinsic signal response to the vertical and horizontal drifting bar was used to generate azimuth and elevation retinotopic maps (Figure S1A). The retinotopic maps were then used to locate V1 and HVAs. Borders between these areas were drawn at the meridian of elevation and azimuth retinotopy manually (Marshel et al., 2011; Smith et al., 2017). The vasculature map then provided landmarks to identify visual areas in two-photon imaging.

In vivo two-photon calcium imaging

Two-photon imaging was carried out using a custom Trepan2p microscope controlled by custom LabView software (Stirman et al., 2016). Simultaneous dual-region imaging was achieved by splitting the excitation beam and temporally multiplexing laser pulses (Stirman et al., 2016). Two-photon excitation light from an 80 MHz Ti:Sapph laser (Newport Spectra-Physics Mai Tai DeepSee) was split into two beams through polarization optics, and one path was delayed 6.25 ns relative to the other. The two beams were independently directed with custom voice-coil actuated steering mirrors and tunable lenses, such that the X, Y, Z planes of the two paths are independently positioned within the full field (4.4 mm diameter). Both beams were scanned by the resonant scanner (4 kHz, Cambridge Technologies), and a single photon signal was collected by a photo-multiplier tube (PMT) (H7422P-40, Hamamatsu), and demultiplexed using outboard electronics prior to digitization. In the current study, two-photon imaging regions of 500 × 500 μm² were collected at 13.3 Hz for two-region imaging or 6.67 Hz for quad-region imaging. Imaging was performed with <80 mW of excitation (910 nm) laser power, as measured out of the front of the objective. Mice recovered in their home cage for at least 2 days after surgery, before acquiring two-photon imaging. Mice were head-fixed ∼ 11 cm from a flat monitor, with their left eye facing the monitor, during imaging. Approximately 70° x 45° of the left visual field was covered. If not otherwise stated, two-photon images were recorded from quiet awake mice. For anesthetized experiments, mice were lightly anesthetized under 1% isoflurane. During two-photon imaging, we monitored mouse pupil position and diameter using a custom-controlled CCD camera (GigE, National Instruments) at 20 - 25 fps. No additional light stimulation was used for pupil imaging.

Visual stimuli

Visual stimulation was displayed on a 60 Hz LCD monitor (9.2 × 15 cm²). All stimuli were displayed in full contrast. For course population RF and single neuron RF mapping (Figure S1B-D), a rectangular (7.5°x 8.8°) bright moving patch containing vertical drifting grating (2 Hz, 0.05 cpd) on a dark background was displayed. The moving patch appeared and disappeared on a random position of the full monitor in pseudo-random order without interruption by a gray screen, and presented on each position for 5 sec.

To characterize the value and structure of the correlation of V1 and HVAs, we showed mice full-screen sine-wave drifting grating stimuli in 8 directions (0 – 315°, in 45° step), with an of 0.02, 0.05 or 0.19 cpd, and a TF of 1, 2 or 8 Hz (72 conditions in total). Each of the sine-wave drifting grating stimuli was presented for 2 s in pseudo-random order. Stimuli with the same SF and TF were presented successively without interruption. A gray screen was presented for 3 seconds when changing the SF and TF of stimuli.

In a subset of experiments, we also characterized the cross-stimulus stability of functional networks using combo stimuli with naturalistic videos and full contrast drifting gratings (at 2 Hz, 0.05 cpd). Two naturalistic videos, each lasting for 32 s were generated by navigating a mouse home cage using a GoPro camera.

Calcium imaging processing

Calcium imaging processing was carried out using custom MATLAB codes (Yu et al., 2022). Two-photon calcium imaging was motion corrected using Suit2p sub-pixel registration module (Pachitariu et al., 2016). Neuron ROIs and cellular calcium traces were extracted from imaging stacks using custom code adapted from Suit2p modules (Pachitariu et al., 2016). Neuropil contamination was corrected by subtracting the common time series (1st principal component) of a spherical surrounding mask of each neuron from the cellular calcium traces (Harris et al., 2016). Neuropil contamination corrected calcium traces were then deconvolved using a Markov chain Monte Carlo (MCMC) method (Pnevmatikakis et al., 2013). For each calcium trace, we repeated the MCMC simulation for 400 times and measured the signal-to-noise of MCMC spike train inference for each cell. For all subsequent analyses, only cells that reliable spike train inference results were included. Neurons with low responsiveness were excluded for subsequent analysis (trial averaged spike count to preferred spatiotemporal frequency summed over all orientations < 1; or trial averaged spike count to a 32 s naturalistic video < 1).

Receptive field

We mapped RFs by reverse correlation of neuronal responses with the locations of the moving patch of drifting grating stimulus. For population RF mapping, population neuronal responses of simultaneously recorded neurons from a 500 × 500 μm² imaging window were reverse correlated with the stimulus locations.

Gaussian mixture model

To characterize the tuning properties unbiasedly, neurons were clustered using a Gaussian mixture model (GMM) based on the trial-averaged responses to the drifting gratings. Only reliable responsive neurons were included for GMM analysis (trial-to-trial Pearson correlation of the inferred spike trains > 0.08, spike trains were binned at 500 ms). Neuronal responses of the whole population pooled over all texted areas, were first denoised and reduced dimension by minimizing the prediction error of the trial-averaged response using the principle components (PC). 45 PCs were kept for population responses to the drifting gratings. We also tested a wide range number of PCs (20 − 70), and we found the tuning group clustering was not affected by the number of PCs. Neurons collected from different visual areas and different animals were pooled together in training GMM. GMMs were trained using MATLAB build function fitgmdist with a wide range number of clusters. A model of 65 classes was selected based on the Bayesian information criterion. We also examined models with smaller (20, 30, and 45) or even larger numbers of classes (75), the overall results held regardless of the number of GMM classes. Figure S2 show the response pattern of GMM classes of drifting gratings. The size of each GMM class is shown in figure S1F.

To summarize the spatial-temporal tuning properties of neurons (Figure 1E), we manually organized the 65 GMM classes into 6 groups based on their preferred SF, TF (Figure S2). Group 1 prefers low SF and low TF (LSLT, 0.02 cpd, 1-2 Hz), group 2 prefers medium SF and low TF (MSLT, 0.05 cpd, 1-2 Hz), group 3 prefers high SF low TF (HSLT, 0.19 cpd, 1-2 Hz), group 4 prefers low SF high TF (LSHT, 0.02 cpd, 8 Hz), group 5 prefers medium SF and high TF (MSHT, 0.05 cpd, 8 Hz), group 6 prefers high SF and high TF (HSHT, 0.19 cpd, 8 Hz) and group 7 not specific. Group 7 included 4 classes that did not exhibit specific response features, among them two classes are extremely small (each contains <5 neurons), and the other two contain neurons with small response strength (mean spike count < 0.5 spikes/s to preferred stimulus). As we have been inclusive in data selection for the GMM training and included low-firing neurons, the latter two classes contain about 1500 neurons in total. It is justifiable to exclude low-response neurons from further analysis. Thus, the whole group 7 was excluded for further analysis.

GMM classification accuracy

We examined the accuracy of GMM classification for neuron responses to drifting gratings. We performed GMM clustering on 10 random subsets of neurons (90% of all neurons). We found the center of the Gaussian profile of each class was consistent (Figure 1C). The same class of different GMMs was identified by matching the center of the class. Then we asked whether a neuron was classified in the same class in each GMM model. We found neurons were consistently classified into the same class in GMMs of a random subset of data (Figure 1D). We also performed GMM on population data after randomly shuffling neuron identity (10 permutations). Classes were identified by matching the center of the class and then grouped following the previous definition. We found that neurons are allocated into the same SF-TF group in GMMs of randomly ordered data (Figure S1G, H). These analyses suggested that GMM provided a reliable classification of neurons.

Orientation and direction selectivity

The direction and orientation selectivity of each neuron were computed using neuron response to its optimal spatial-temporal frequency of drifting grating stimuli. The direction selectivity index and orientation selectivity index were computed using the following equations.

The polar plots of tuning groups were generated by averaging responses to the preferred direction of each neuron within a tuning group, and normalized to one (Figure S3C). For neurons with high direction selectivity, neuron responses to preferred direction was considered, while for neurons with low direction selectivity (DSI < 0.5), neuron responses to both preferred and null directions were included.

ISOI warping

We spatially registered ISOI map of V1 to align with that of LM or AL. We first segmented the ISOI map by color segmentation using K-means clustering and then determined the center of each color segment. Then we performed the affine transformation of color band centers of V1 to match that of LM or AL. The transformation matrix M was determined by minimizing the distance between transformed V1 centers and LM or AL centers using Matlab function fminsearch.

Correlation calculation

Noise correlation was defined as the trial-to-trial correlation of residual spike count (1 s time window, if not otherwise stated) after subtracting the mean response to each stimulus of the 72-condition sine-wave drifting gratings. Residual spike count to all stimuli (eg. gratings with different directions and SFs and TFs), and all trials were concatenated into one column vector per neuron (u_i, i = 1,2… N, neuron identity). The noise correlation r_sc was computed as the Pearson correlation of u_i and u_j.

i, j indicate neuron identity. Signal correlation was defined as the neuron-to-neuron Pearson correlation of mean responses. Mean response was a 72-element column vector, computed by trial averaging responses to sine-wave gratings with 72 conditions. To examine the relation between noise correlation and joint firing rate between a pair of neurons. We computed the mean joint spike count (geometric mean spike count average over all stimuli.

We computed inter-area NCs with simultaneously recorded regions that shared greater than 40% of population RF. We kept this criterion even though we did not detect a relationship between the inter-area NC and the fraction of population RF overlap within the tested range (p = 0.37).

Fidelity of noise correlation measurement

Tolerance of correlation calculation to inaccuracy in spike train inference

We quantify the spike train inference accuracy using a previously published data set with simultaneous cell-attached recording and two-photon imaging of GCamp6s from mouse V1 (Chen et al. (2013); http://crcns.org/data-sets/methods/cai-1). We performed spike train inference on the recordings with stable baseline and good correspondence between calcium trace and electrophysiology recording (linear correlation, r > 0.1; bin 0.1 s; Figure S4A, B). The signal-to-noise (SNR) of the calcium trace of the calibration data is 12.3 ± 5. It is comparable with the SNR of the calcium signal of the current study (8.7 ± 1.8).

We further evaluated how the correlation calculation was affected by inaccurate spikes train recovery. We took publicly available electrophysiology recordings of mouse V1 neurons (Theis et al. (2016); http://spikefinder.codeneuro.org/), and computed residual spike count correlation at 1 s time bin after perturbations on the ground truth spikes train. We did four types of perturbations, (1) randomly missing spikes; (2) missing isolated spikes as the signal-to-noise of the calcium signal of isolated spikes may be low; (3) missing all spikes within a burst; (4) missing 60% spikes within a burst (Figure S4C). We identified isolated spikes or burst spikes by thresholding the inter-spike-interval of each spike. A spike that was >t s distance away from spikes franking itself was a t isolated spike. A spike that was <t s distance away from another spike was a t burst spike. The residual spike count correlation computed with perturbed spikes trains was linearly correlated with ground truth (Figure S4D) and exhibited good tolerance to up to 60% missing spikes by all types of perturbation (fidelity > 0.6; Figure S4E).

Significance of noise correlation

Since the value of noise correlations was small, we tested whether these values were significantly above zero. We compared the noise correlation with trial-shuffled noise correlation, the latter was computed using trial-shuffled data (the order of trials were randomized for each neuron independently). The population-mean noise correlation computed with trial-aligned data was significantly higher than that of the trial-shuffled data with the size of the experimental population (Figure 2B).

Accuracy of noise correlation

We investigated the accuracy of noise correlation estimation with both data and model. The individual noise correlations of the same set of neurons varied when computing using a different random subset of trials (Figure 2D). We computed the population mean value of the noise correlation of a random subset of neuron pairs and calculated the confidence interval for estimating the population mean noise correlation. The accuracy of population-mean estimation increases with the number of neurons, even with a limited number of trials (Figure 2C). We further characterized the estimation accuracy by simulating correlated neuron population (Macke et al., 2009), which allows an arbitrary number of trials. The expected firing rate and expected population mean correlation match our experimental data. To achieve an accurate estimation (1/10 standard error/mean value) of the population mean correlation converges with >100 neurons even using experimental level trial numbers (Figure 2C).

Community module analysis

We constructed a V1-HVA connectivity matrix using the fraction of high NC (NC > mean + 2.5*SD of trial-shuffled NC) pairs between each GMM class. We performed community detection analysis using the Louvain algorithm (Rubinov and Sporns, 2010), which assigned densely connected nodes to the same module. The spatial smooth parameter γ that generated the largest deviation from a random connectivity matrix is picked. The analysis was performed using the Brain Connectivity Toolbox (brain-connectivity-toolbox.net).

Leaky integrate-and-fire neuron network simulation

LIF simulations were carried out using the Brian2 simulation engine in Python (Stimberg et al., 2019). The LIF neuron network model was defined similarly as Song et al. (2000). In brief, the membrane potential of LIF neurons was given by the equation below:

Where τ_m corresponds to the membrane time constant (20 ms). ge, gi and Ee, Ei are the excitatory and inhibitory synaptic conductance and their respective reversal potential (Ee = 0 mV, Ei = −80 mV). the membrane potential was simulated with a time resolution of dt = 0.1 ms. El (−70 mV) corresponds to the resting potential. Connections between LIF neurons occurred with probability p = 0.02, and the strength of the connections is defined as W_ij, (i, j indicate source and target neuron ID, i = j). W_ij = J_max or 0 if not connection. J_max defined specifically in each simulation setting (Figure 6D, E). The dynamics of synaptic conductance were given by exponential decay functions ge/dt = −ge/τ_e and gi/dt = −gi/τ_i. τ_e (5 ms) and τ_i (10 ms) are the decay time constants for excitatory and inhibitory synapses.

The LIF network received feedforward input from Poisson neurons (N = 1000 in the network simulation, and N = 80 in the toy model), whose firing follows time-varying Poisson processes (0-30 Hz). The Poisson input neurons are connected to LIF neurons with probability p = 0.2, and the strength of the connections is defined as S_ij. S_ij = 0.2 or 0 if not connection. In the toy model, the connectivity S_ij is defined by the fraction of shared input. In the network simulation, the connectivity S_ij is defined by the equation below:

a, b, and c are parameters that manipulate the structure of the FFD connection, ranging from fully random to fully bell-shaped (Figure 6E).

In the toy model, the Poisson neuron firing constant at 5 Hz. In the network simulations, the instant firing rate of Poisson neurons is defined by a combination of five Gaussian profiles (Figure 6C).

Supplementary Information

Figure S1.

Figure S2.

Figure S3.

Figure S4.

Figure S5.

Figure S6.

Figure S7.

Supplementary Table

Table 1.