Functional interactions among neurons within single columns of macaque V1

  1. Ethan B Trepka
  2. Shude Zhu  Is a corresponding author
  3. Ruobing Xia
  4. Xiaomo Chen
  5. Tirin Moore
  1. Department of Neurobiology, Howard Hughes Medical Institute, Stanford University, United States
  2. Neurosciences Program, Stanford University, United States
  3. Center for Neuroscience, Department of Neurobiology, Physiology, and Behavior, University of California, Davis, United States

Abstract

Recent developments in high-density neurophysiological tools now make it possible to record from hundreds of single neurons within local, highly interconnected neural networks. Among the many advantages of such recordings is that they dramatically increase the quantity of identifiable, functional interactions between neurons thereby providing an unprecedented view of local circuits. Using high-density, Neuropixels recordings from single neocortical columns of primary visual cortex in nonhuman primates, we identified 1000s of functionally interacting neuronal pairs using established crosscorrelation approaches. Our results reveal clear and systematic variations in the synchrony and strength of functional interactions within single cortical columns. Despite neurons residing within the same column, both measures of interactions depended heavily on the vertical distance separating neuronal pairs, as well as on the similarity of stimulus tuning. In addition, we leveraged the statistical power afforded by the large numbers of functionally interacting pairs to categorize interactions between neurons based on their crosscorrelation functions. These analyses identified distinct, putative classes of functional interactions within the full population. These classes of functional interactions were corroborated by their unique distributions across defined laminar compartments and were consistent with known properties of V1 cortical circuitry, such as the lead-lag relationship between simple and complex cells. Our results provide a clear proof-of-principle for the use of high-density neurophysiological recordings to assess circuit-level interactions within local neuronal networks.

Editor's evaluation

This is an important paper which shows how high-density neurophysiological recordings in non-human-primates can be used to identify inter-neuronal interactions based on cross-correlations. This provides valuable insights such as the dependence of correlations on vertical distance and orientation tuning. Overall the techniques used here are compelling and set a standard for recordings in non-human-primates. The paper is of interest for a broad audience of neuroscientists that performs electrophysiological recordings or is interested in functional interactions among neuron pairs.

https://doi.org/10.7554/eLife.79322.sa0

Introduction

Understanding the functional logic of local neuronal microcircuits is among the more fundamental objectives in the study of neural systems, yet it is also among the most challenging. This seems particularly true for mammalian neocortical circuits involved in perceptual and cognitive functions, and most notably in nonhuman primate model systems for which the available tools to interrogate those circuits are the most limited. The columnar organization of the mammalian neocortex (Lorente de No, 1938; Mountcastle, 1957) and its distinctly layered structure within different cortical domains are both widely appreciated (Reviewed in DeFelipe et al., 2012; Horton and Adams, 2005; Mountcastle, 1997). In addition, several key principles of cortical circuitry, including constituent cell types (Harris and Mrsic-Flogel, 2013; Jiang et al., 2015; Kätzel et al., 2011; Markram et al., 2004; Network, 2021; Packer and Yuste, 2011; Yoshimura and Callaway, 2005), input-output organization (Callaway, 1998; Douglas and Martin, 2004; Lefort et al., 2009; Muñoz-Castañeda et al., 2021; Thomson and Bannister, 2003; Weiler et al., 2008) and local microcircuit motifs (Avermann et al., 2012; Frandolig et al., 2019; Karnani et al., 2016; Obermayer et al., 2018; Pfeffer et al., 2013; Pi et al., 2013) have emerged in recent years. Although it remains to be determined, such principles may turn out to generalize not only across neocortical areas, but also across species (Harris and Shepherd, 2015; Karten, 2015; Stacho et al., 2020) (see also Campagnola et al., 2022; Wildenberg et al., 2021). Yet, mapping complete cortical microcircuits within even a single cortical area remains a tremendous challenge (Adesnik and Naka, 2018).

Recent advances in recording technology have facilitated the development of large-scale, high-density micro-electrode arrays resulting in a substantial increment (>10x) in the number of neurons that can be studied simultaneously within a localized area of neural tissue. A prime example is the recent development of the Neuropixels probe (IMEC, Inc), which consists of a high-channel count Si shank with continuous, dense, programmable recording sites (~1000/cm). Numerous recent studies have demonstrated the advantages of such probes, such as their use in recording large neuronal populations within deep structures where optical approaches cannot be deployed (Jun et al., 2017; Steinmetz et al., 2019). In addition, the high-density capacity of such probes dramatically increases the quantity of single neurons that can be obtained within a localized area of neural tissue (Siegle et al., 2021), thus making them well-suited for investigations of local neuronal circuitry. Given that studies of local neuronal circuitry within the primate brain are notoriously difficult to achieve, high-density electrophysiological approaches may be particularly valuable. However, only a few electrophysiological studies of the primate brain using such probes have been carried out thus far (Hesse and Tsao, 2020; Paulk et al., 2022; Sun et al., 2022; Trautmann et al., 2019; Zhu et al., 2020).

To date, many studies have leveraged the covariation in spiking activity between simultaneously recorded neurons to elucidate underlying neural mechanisms in the primate brain with some success, particularly within the visual system (Briggs et al., 2013; Chu et al., 2014; Hansen et al., 2012; Hembrook-Short et al., 2019; Jia et al., 2013; Kohn and Smith, 2005; Koren et al., 2020; Krüger and Aiple, 1988; Maldonado et al., 2000; Smith and Kohn, 2008; Zandvakili and Kohn, 2015). In particular, temporally precise correlations in spiking activity have provided a unique means of assessing interactions among neurons in both local and distributed networks (Aertsen and Gerstein, 1985; Aertsen et al., 1989; Diba et al., 2014; Moore et al., 1970; Nelson et al., 1992; Nowak et al., 1999; Perkel et al., 1967; Siegle et al., 2021), and identification of such interactions has played an important part in understanding neural circuits in the mammalian visual system (Alonso and Martinez, 1998; Alonso et al., 1996; Alonso et al., 2001; Baker and Bair, 2012; Cohen and Kohn, 2011; Das and Gilbert, 1999; Denman and Contreras, 2014; Michalski et al., 1983; Nelson et al., 1992; Reid and Alonso, 1995; Schwarz and Bolz, 1991; Senzai et al., 2019; Siegle et al., 2021; Toyama et al., 1981a; Ts’o et al., 1986; Usrey et al., 1998; Usrey et al., 1999). However, the extent of circuit-level details addressable with crosscorrelation is greatly limited by the low incidences of simultaneous recordings from connected neurons when using conventional extracellular recording techniques (e.g. Alonso et al., 2001; Hembrook-Short et al., 2019; Nelson et al., 1992; Ts’o et al., 1986). The use of high channel-count probes should substantially mitigate that limitation by virtue of the large increment in recording yield. Moreover, the high-density of recordings should further increase the incidence of identifiable correlated neuronal activity by virtue of the proximity of recorded cells.

We assessed the capacity of high-density Neuropixels probes to identify functional interactions among pairs of neurons within cortical columns of primary visual cortex of macaque monkeys using established crosscorrelation approaches. Crosscorrelation assesses the statistical dependencies between spike trains of two or more neurons and has long played an important role in estimating how ensembles of neurons interact with one another (Casile et al., 2021; Okatan et al., 2005; Perkel et al., 1967). Although a number of different analyses have been employed (Casile et al., 2021; Keeley et al., 2020; Kobayashi et al., 2019), crosscorrelation is the most widely used, perhaps due to its simplicity. Significant crosscorrelations are broadly interpreted as identifying ‘functional connections’ (Aertsen and Gerstein, 1985; Dann et al., 2016; Denman and Contreras, 2014; Menz and Freeman, 2004; Schwarz and Bolz, 1991) or ‘functional interactions’ among neurons (Gochin et al., 1991; Ostojic et al., 2009), presumably to distinguish them from more direct measurements of synaptic effects. These functional interactions are interpreted as reflecting one of myriad putative circuit arrangements among neuronal pairs (Aertsen and Gerstein, 1985; Melssen and Epping, 1987; Moore et al., 1970; Ostojic et al., 2009), arrangements which include direct monosynaptic connections (Alonso and Martinez, 1998; English et al., 2017; Reid and Alonso, 1995) or neurons with common input (Constantinidis et al., 2001; Das and Gilbert, 1999; Denman and Contreras, 2014; Michalski et al., 1983; Ts’o et al., 1986; Türker and Powers, 2001), either of which can provide important insight into local network architecture. In this study, we identified 1000s of functionally interacting neuronal pairs during single recordings from neurons situated in different cortical layers. Our results demonstrate robust, systematic variations in the synchrony and strength of functional interactions within cortical columns. In addition, by leveraging the large numbers of interacting pairs, distinct classes of interactions could be identified within the full population.

Results

Identifying functional interactions within single columns of visual cortex

The activity of V1 neurons was recorded in two anesthetized macaque monkeys (M1, M2) using high-density, multi-contact Neuropixels probes (version 3A; IMEC Inc, Belgium; Figure 1a; Methods). Each probe consisted of 986 contacts (12 mm x 12 mm, 20 µm spacing) distributed across 10 mm, of which 384 contacts could be simultaneously selected for recording. Probes were inserted into the lateral operculum of V1 with the aid of a surgical microscope at angles nearly perpendicular to the cortical surface. The dense spacing between electrode contacts provided multiple measurements of the waveforms from individual neurons (mean = 4.52 measurements) (Figure 1a and b) and facilitated the isolation of large numbers of single neurons. In each of the 5 experimental sessions (3 in M1, 2 in M2), we measured the visual responses of 115–221 simultaneously recorded neurons to drifting gratings presented at varying orientations (total = 802 neurons). As expected, neurons were highly orientation selective, and exhibited both simple and complex cell properties (De Valois et al., 1982; Hubel and Wiesel, 1962; Hubel and Wiesel, 1968; Figure 1b). The ratio of simple to complex neurons, respectively, was 1:2.4; 236/802 neurons were simple, and 566/802 neurons were complex. As in previous studies (Briggs et al., 2013; Hembrook-Short et al., 2019; Jia et al., 2013; Kohn and Smith, 2005; Siegle et al., 2021; Smith and Kohn, 2008; Zandvakili and Kohn, 2015), we used the visually driven spike trains to measure crosscorrelations between simultaneously recorded neuronal pairs.

Identifying functional interactions within single columns of visual cortex.

(a) Upper left cartoon depicts the angle of Neuropixels probe penetrations made into the lateral surface and underlying calcarine sulcus of V1. Upper right, Neuropixels probe base and shank, and layout of electrode contacts for a section of the recording shank. Lower, raw voltage traces recorded from an exemplar section of channels and time period. (b) Example single-neuron recordings with Neuropixels probes, three simple cells (orange, blue, green) and one complex cell (purple). Top, spike waveforms recorded across multiple adjacent electrode contacts are shown for each neuron. Bottom, each neuron’s response to its preferred orientation (rasters and instantaneous spike rates) and their corresponding tuning curves. Red arrows (upper left) denote the drift direction of oriented gratings. (c) Example CCG between an example pair of V1 neurons. Corrected CCGs were generated from the difference between a jittered and an uncorrected CCG. Significance of each CCG was determined from comparisons between the peak and the noise distribution. (d) Distribution of ratios of CCG peaks to the noise (SD) for all recorded pairs. Shaded area denotes CCGs with peaks >7 SDs above the mean of the noise distribution. (e) Two example CCGs differing in both peak lag and peak efficacy.

To estimate the functional interactions between pairs of neurons recorded simultaneously within columns of V1, we computed cross-correlograms (CCGs) using the 802 visually responsive neurons recorded across sessions. CCGs were computed from the spike trains of 68,579 pairs of simultaneously recorded neurons (6,555–24,310 pairs/session, Methods). Each CCG was normalized by the firing rate (FR) and jitter-corrected to mitigate the influences of FR (Bair et al., 2001; Mastronarde, 1983) and correlated slow fluctuations (Harrison and Geman, 2009; Smith and Kohn, 2008), respectively, yielding a corrected CCG (Figure 1c). In addition, as in previous studies, we considered a CCG significant only if its peak occurred within 10ms of zero time lag, and if that peak was >7 standard deviations above the mean of the noise distribution (Siegle et al., 2021). Using these criteria, a total of 10,246 significant CCGs were obtained from all recording sessions (Figure 1d), with each session yielding 755–3,022 significant CCGs. The peak lag of each CCG, defined as the differences between zero and the time when the peak occurred, estimates the synchrony and/or direction of functional interactions between neuronal pairs; whereas the peak efficacy measures the strength of interactions (Figure 1e).

Variation in the synchrony and strength of functional interactions within cortical columns

A number of previous studies using low-channel count probes or chronically implanted electrode arrays have shown that correlated activity in primate V1 declines with the horizontal distance separating pairs of neurons (Krüger and Aiple, 1988; Maldonado et al., 2000; Smith and Kohn, 2008) (see also Chu et al., 2014). Evidence from these studies suggest that correlations are greatest for pairs of neurons located within the same column, and diminish with greater columnar distance. Other evidence shows variation in the spike timing correlations between neuronal pairs located within different laminar compartments (Smith et al., 2013). However, considerably less is known about how the nature of correlations varies across the depth of individual columns where the degree of shared input and connectivity is at its highest. We therefore leveraged the large numbers of significantly correlated pairs obtained from high-density recordings to examine how the synchrony and strength of correlations depended on the vertical distance separating neurons within V1 columns. Figure 2a shows data from an example recording session in which 221 visually responsive neurons were recorded and 2,453 significantly correlated pairs were obtained. All neurons are shown along the ~2 mm depth of cortex. Shown also are two example correlated pairs whose CCGs are shown in Figure 1e. Of the two pairs, the vertical distance separating neurons in one pair was 138 µm greater than that of the other. In spite of this small difference, the CCG of the closer pair was both more synchronous and stronger than the more distant pair. This pattern of results was observed across all significantly correlated pairs and in all sessions (Figure 2b–c; Figure 2—figure supplements 12). The synchrony of correlated spiking diminished several fold across neuronal pair distance. This change could be fit with a linear function (r=0.42; p<10–5) in which the (absolute) peak lag increased at a rate of 1.3 ms / 500 µm of vertical distance. Peak efficacy of the significant CCGs also depended heavily on pair distance. This effect could be fit with an exponential decay function (r=–0.34; p<10–5) in which the peak efficacy decreased by half within 154 µm. Thus, both measures of functional interactions depended heavily on the vertical distance separating neuronal pairs. In addition, we confirmed that the effects of vertical distance on both the synchrony and strength of CCGs were independent of whether neuronal pairs were located within the same or different cortical layers (Figure 2—figure supplement 3).

Figure 2 with 3 supplements see all
Dependence of synchrony and strength of functional interactions on vertical distance within single cortical columns.

(a) Example session (M1, session 3) with 221 visually responsive neurons recorded simultaneously and their locations across cortical depth. (Horizontal axis is magnified for visualization). Cortical depth 0 denotes the boundary between Layer 4c and Layer 5. Laminar boundaries were determined using histological data and current-source-density (CSD) profile for each session (Methods). Two example correlated pairs from Figure 1e with varied CCGs are shown in color (blue-yellow pair and pink-red pair corresponding to Figure 1e top and bottom, respectively). (b) Linear dependence of synchrony on vertical pair distance. (c) Strength of CCGs decay with greater pair distance. In b and c, all significantly correlated pairs from all sessions are combined and each dot denotes the mean peak lag or median peak efficacy of significantly correlated CCGs within a (10% quantile) vertical distance bin. Error bars denote 95% confidence intervals. Black lines denote the linear and exponential fits in b and c, respectively; slope (b) and decay constant (λ) are shown. Red lines and bar plots show marginal distributions.

Dependence of synchrony and strength of functional interactions on tuning similarity

In addition to the dependence of correlated activity on the distance between neuronal pairs, many studies have shown that greater functional and synaptic connectivity typically occurs between neurons with similar stimulus preferences (Chu et al., 2014; Constantinidis et al., 2001; Cossell et al., 2015; DeAngelis et al., 1999; Denman and Contreras, 2014; Funahashi and Inoue, 2000; Lee et al., 2016; Ts’o et al., 1986) (but see Das and Gilbert, 1999; Maldonado et al., 2000). Within primate V1, stimulus selectivity is notably similar for neurons within the same column, particularly for orientation selectivity (Blasdel and Salama, 1986; Hubel and Wiesel, 1968; Hubel and Wiesel, 1974; Ts’o et al., 1990), and this was evident in our recording sessions, where the peak visual responses were largely aligned at the same stimulus orientation across cortical depth (Figure 3a). We considered that within orientation columns, functional interactions could be homogenous for populations of similarly tuned neurons. Alternatively, it could be that even small variations in tuning similarity could result in robust differences in the synchrony and strength of correlated activity. To address this, we examined the dependence of synchrony and strength on the similarity of visual properties of neurons within the same cortical column. As in previous studies (Shadlen and Newsome, 1998; Zohary et al., 1994), we quantified tuning similarity by computing signal correlations (rori) for each neuronal pair (Methods). Across the total number of neuronal pairs (N=68,579), the mean rori was 0.25. For the significantly correlated neuronal pairs, the mean rori was 0.33. Signal correlations for the two previous example neuronal pairs are shown in Figure 3b. The responses of both pairs are positively correlated, yet that correlation is much higher in the second, more proximal, pair (Figure 2a) and the one with a more synchronous and stronger CCG (Figure 1e). Overall, we found that both the peak lag and peak efficacy of CCGs for significantly correlated neuronal pairs varied monotonically with tuning similarity across the range of signal correlations (Figure 3c and d). Neuronal pairs with the highest signal correlations exhibited half the peak lags and twice the peak efficacies of uncorrelated pairs. This pattern was observed in each of the individual recording sessions (Figure 3—figure supplements 12).

Figure 3 with 2 supplements see all
Dependence of synchrony and strength of functional interactions on tuning similarity within single cortical columns.

(a) Heat map of visual responses across drift directions of oriented gratings and across cortical depth. The response tuning of each of 802 neurons was aligned to the overall preferred orientation shared by neurons recorded from the same session, and all sessions were combined. (b) Signal correlation between exemplar neurons. Left, Scatter plot of normalized responses to different stimulus orientations (n=36) for the two example pairs shown in Figure 1e and Figure 2a. Signal correlations (rori) are also shown. Right, each neuron’s orientation tuning curve. (c) Linear dependence of synchrony on the corresponding signal correlation. (d) Linear dependence of CCG strength on the corresponding signal correlation. (e) Difference in peak lag of distance-matched CCGs was negatively correlated with difference in signal correlation. (f) Difference in peak efficacy of distance-matched CCGs was positively correlated with difference in signal correlation. In c-f, all significantly correlated pairs from all sessions are combined and each dot denotes mean peak lag or median peak efficacy of significantly correlated CCGs within a (10% quantile) signal correlation bin. Error bars denote 95% confidence intervals. Black lines denote linear fits; slopes (b) are shown. Red lines and bar plots show marginal distributions.

We considered that the apparent relationship between the synchrony and strength of functional interactions and signal correlation might result indirectly from a collinear effect of vertical distance on CCGs (Figure 2). To address this, we examined differences in the peak lags and peak efficacies of CCGs between combinations of two neuronal pairs separated by comparable cortical distances. Specifically, we sorted all significantly correlated CCGs by their vertical distances, and then examined whether differences in signal correlations (rori) among adjacently sorted (distance-matched) pairs were still associated with differences in CCG peak lags and peak efficacies (Methods). Indeed, we found that the differences in peak lags of distance-matched CCGs were negatively correlated with signal correlation (Figure 3e) and the differences in peak efficacies of distance-matched CCGs were positively correlated with signal correlation (Figure 3f). These results indicate that signal correlations within the column predicted both the synchrony and strength of functional interactions independent of vertical pair distance. Nonetheless, the distance-matched correlations (Figure 3e–f) were smaller than their corresponding unmatched correlations (Figure 3c–d), suggesting that the vertical distance between neurons and their orientation signal correlations exhibit distinct, but overlapping, effects on the timing and strength of functional interactions within a single cortical column.

To quantify the distinct contributions of vertical pair distance and orientation signal correlation to the synchrony and strength of CCGs, we fit GLMs to predict CCG peak lag and peak efficacy using pair distance and signal correlation as predictors (Methods). Predictors were standardized (z-scored) so that their relative effects could be compared, and peak outliers (1.5*IQR criterion) were removed. The resulting regression equations were:

peaklag=2.48+0.94×pairdist.0.37×rori
peakefficacy=0.0180.0017×pairdist.+0.0029×rori

Regressions explained 19% of variance in peak lag (R2=0.191) and 20% of variance in peak efficacy (R2=0.195). Because predictors were standardized, the regression coefficients capture the change in peak lag/efficacy associated with a 1 standard deviation (SD) increase in pair distance/signal correlation. In the regression predicting CCG peak lag, a 1 SD increase in pair distance was associated with a 0.94ms increase in peak lag whereas a 1 SD increase in signal correlation was associated with a 0.37ms decrease in peak lag. Thus, for CCG peak lag, the coefficient of pair distance was nearly three times the coefficient of signal correlation. In contrast, for CCG peak efficacy, the coefficient of signal correlation was nearly twice that of pair distance. Thus, whereas signal correlation was less predictive of CCG peak lag, it was more predictive of CCG peak efficacy than pair distance.

Classification of functional interactions

CCG peak lags and peak efficacies are often the parameters of interest in cross correlations (Briggs et al., 2013; Hembrook-Short et al., 2019; Smith and Kohn, 2008), yet they are simplifications of the more complex, underlying crosscorrelation functions. The shape of these correlation functions may offer additional insights into the distinct types and properties of functional interactions present among neurons within a network. Several theoretical studies have suggested a correspondence between CCG shape and underlying pairwise connectivity (Aertsen and Gerstein, 1985; Melssen and Epping, 1987) that can be further influenced by overall network structure and background noise (Ostojic et al., 2009). For example, synchronous CCGs tend to correspond to pairs of neurons that receive input from a common source, while asynchronous CCGs tend to correspond to pairs that have direct synaptic connections (Ostojic et al., 2009). Moreover, synchronous CCGs with narrow peaks and synchronous CCGs with broad peaks may correspond to pairs of neurons that receive input from common sources with shorter and longer autocorrelation timescales, respectively (Ostojic et al., 2009). Experimental studies have corroborated these findings and identified similar CCG shapes in different cortical regions and species (Alonso and Martinez, 1998; Constantinidis et al., 2001; Hembrook-Short et al., 2019; Siegle et al., 2021). However, the distribution of these CCG shapes within a single cortical column remains unknown. Furthermore, whether that distribution within V1 corroborates other evidence about the functional and/or anatomical relationships among V1 laminae and cell types remains unclear.

To address these questions, we clustered the entire population of CCGs, taking advantage of the large number of significantly correlated pairs to identify robust CCG templates. To do this, we first normalized significant CCGs, and utilized t-distributed Stochastic Neighbor Embedding (t-SNE) to map CCGs to a lower dimensional space, and then clustered CCGs in the resulting space using k-means (Methods). To select a statistically reasonable number of clusters, we examined how the total variance explained by clustering and the silhouette score changed as a function of the number of clusters (Figure 4a). From this, we selected four as the optimal number of clusters given that silhouette score peaked ~3–4 clusters, and 4 clusters explained more variance than 3.

Identification of distinct classes of functional interactions within the full population.

(a) Explained variance (left) and silhouette score (right) as a function of number of clusters. The dashed vertical line indicates the selected number of clusters (n=4). (b) Scatter plot of dimensionality-reduced CCGs in the first two dimensions of t-SNE space. Randomly selected example CCGs are overlayed on the scatterplot in their corresponding location in t-SNE space. (c) CCG templates generated by averaging over all the CCGs in each cluster. The templates include a ‘sharply synchronous’ class (Ssync) with a narrow peak at τ=0, a ‘broadly synchronous’ class (Bsync) with a wide peak at τ=0, a ‘forward’ class (Fasync) (leading) with more probability density after τ=0, and a ‘reverse’ class (Rasync) (lagging) with more probability density before τ=0. Numbers denote the percentage of each class among all significantly correlated pairs.

CCG shape was relatively heterogenous within each of the four clusters (Figure 4b). Nonetheless, by averaging over all CCGs in each cluster, we could construct CCG templates that summarized key characteristics of the clusters (Figure 4c). Within the full population, we identified two synchronous classes of functional interactions, a ‘sharply synchronous’ class (Ssync) with a narrow peak at τ=0 and a ‘broadly synchronous’ class (Bsync) with a wide peak at τ=0. In addition, two asynchronous classes were identified, a ‘reverse’ class (Rasync) (lagging) and a ‘forward’ class (Fasync) (leading) with more probability density before and after τ=0 (median τ=3ms), respectively (Figure 4c). Aside from clear differences in peak lags between subsets of the putative classes (e.g. synchronous vs. asynchronous), CCGs of different classes also differed in their peak efficacies; synchronous classes exhibited higher average peak efficacies than asynchronous classes (median peak efficacy: Ssync 0.021, Bsync 0.020, Fasync 0.015, Rasync 0.014). Importantly, our objective was not to find the exact number of distinct classes of functional interactions in V1 or to perfectly categorize every interaction into a homogenous cluster. Instead, we sought to identify at least one set of clusters that is consistent with that expected in local microcircuits.

Corroboration of putative CCG classes with V1 microcircuitry

We next examined the extent to which the putative CCG classes were also distinguishable from one another along anatomical and functional lines given other known properties of V1 microcircuits. First, we considered that the identified classes might differ in their vertical pair distances and signal correlations. Indeed, we found that vertical pair distances were larger and orientation signal correlations were smaller in asynchronous (Fasync and Rasync) than in synchronous (Ssync and Bsync) classes (Figure 5a–b) (significant pairwise comparisons: p<10–5). Given that both the peak lag and peak efficacy components of CCGs were clearly predicted by distance and signal correlation (Figures 23), the observed difference between the synchronous and asynchronous classes is expected. However, additional differences emerged between the identified synchronous classes. For example, we found that in spite of exhibiting similar CCG peak lags, Bsync pairs were separated by greater vertical distances than Ssync ones (Figure 5a; p<10–5). Furthermore, in spite of being separated by a greater cortical distance, Bsync pairs exhibited higher signal correlations than Ssync pairs (Figure 5b; p<10–5). These findings thus provide some validation of the apparent subtypes of CCGs.

Figure 5 with 1 supplement see all
Corroboration of putative classes of functional interactions with V1 microcircuitry.

(a, b) Boxplots of vertical pair distance (a) and orientation signal correlation (b) across the 4 identified CCG classes. Boxplots illustrate the medians, first and third quartiles, and non-outlier (1.5*IQR method) minima and maxima. Asterisks denote significant differences in medians between pairs of classes (Wilcoxon rank-sum test; P<0.05, Bonferroni corrected). (c) Percentage of CCGs in each class composed of two neurons from the same (‘within layer’) or different layers (‘between layer’). In a-c, the reference neuron in a neuronal pair was selected randomly. (d) Distribution of putative CCG classes among neuronal pairs composed of a simple cell as the reference neuron and complex cell as the target neuron. Only the pairwise comparison between Fasync and Rasync is illustrated to show the direction of functional interactions between simple and complex cells.

Next, we examined whether the identified classes of functional interactions differed in their laminar distributions. Indeed, we found that different classes were differentially distributed across V1 layers such that one or more of the identified classes was often overrepresented among functional interactions within particular layers (Figure 5—figure supplement 1; p<10–5). To simplify this result, we compared the proportion of CCGs in each class composed of two neurons in the same layer or different layers. We found that most of the asynchronous pairs were composed of neurons from different layers, while most of the synchronous pairs, particularly the Ssync ones, were composed of neurons from the same layer (Figure 5c) (within vs between proportion: Ssync [0.65 vs 0.35], Bsync [0.52 vs 0.48], F [0.29 vs 0.71], R [0.27 vs 0.73], one proportion z-test: Ssync, F, R: P<10–5; Bsync: p=0.022). This observation dovetails the relationship between pair distance and CCG class described above. Nonetheless, we found that cortical layer had an independent effect of distance on CCG class assignment among nearby (~200 μm) pairs of neurons. The location of neuronal pairs within the same or different, nearby layers predicted whether pairs belonged to the Bsync class and the asynchronous classes (Fasync or Rasync), but not the Ssync class, when controlling for the effects of vertical distance (Table 1; logistic regression, p<10–2). More specifically, CCGs composed of two neurons within the same layer had a higher probability of falling in the Bsync class and a lower probability of falling into the asynchronous classes than CCGs with comparable vertical distances composed of two neurons in different cortical layers.

Table 1
Dependence of putative classes on laminar pairing and vertical distance for pairs of neurons separated by 86–310 μm.
Dependent variable (in/out of cluster)PredictorCoefficientStandard errorp-Value
SsyncDistance–0.0043 /μm0.0014.00*10–5
Layer0.104/layer0.1260.41
BsyncDistance–0.0017 /μm0.00079.24*10–3
Layer0.244/layer0.0854.11*10–3
Fasync or RasyncDistance0.0032 /μm0.00054.23*10–7
Layer–0.262/layer0.08021.13*10–3
  1. Coefficients, standard errors, and p-values from logistic regressions predicting class membership using the distance between pairs of neurons and whether pairs were located in the same or different cortical layer(s). Only pairs of neurons with pair distances greater than the 5% of pairs located in different cortical layers (>86 μm) and less than 5% of pairs located in the same cortical layer (<310 μm) were included. Significant predictors are highlighted.

In addition to the laminar organization, V1 neurons exhibit clear differences in their receptive field properties. In particular, V1 neurons classically fall into two broad functional types: simple (S) and complex (C) cells (De Valois et al., 1982; Hubel and Wiesel, 1962; Hubel and Wiesel, 1968; Movshon et al., 1978; Skottun et al., 1991) (see also Chance et al., 1999; Mechler and Ringach, 2002; Priebe et al., 2004). Among the significant CCGs, a majority were comprised of pairs of complex cells (S/S=5.2%, C/C=63%; S/C=31.8%; one-proportion z-test: p<10–5). Complex cells appear to receive converging input from groups of simple cells, and thus simple cells should lead rather than lag complex cells in their CCGs. To test this in our data, we compared the distribution of putative CCG classes among significantly correlated neuronal pairs composed of a simple cell as the reference neuron and a complex cell as the target neuron (Figure 5d). We found that the proportion of forward (Fasync) CCGs was larger than the reverse (Rasync) class (p<10–5). Notably, although the dominant lead-lag relationship between simple and complex cells is consistent with established models of V1 (Alonso and Martinez, 1998; Martinez and Alonso, 2001; Yu and Ferster, 2013), there were also many CCG pairs in which complex cells led simple cells or where the pair fired synchronously. This heterogeneity in functional interactions between simple and complex cells is consistent with studies suggesting that simple and complex cells might arise from variations in a continuous process as opposed to being two clearly distinct populations (Chance et al., 1999; Kim et al., 2021; Mechler and Ringach, 2002; Priebe et al., 2004).

Corroboration of different classes with V1 input and local circuitry

Previous studies have characterized the anatomical organization of dorsal lateral geniculate nucleus (dLGN) input to V1 in extensive detail (Blasdel and Lund, 1983; Hendrickson et al., 1978; Hubel and Wiesel, 1972). In the macaque brain, dLGN magnocellular and parvocellular axons primarily project to V1 layers 4cα and 4cβ, respectively, along with inputs that terminate in layer 6 (reviewed in Briggs and Usrey, 2011; Lund, 1988; Merigan and Maunsell, 1993; Nassi and Callaway, 2009; Figure 6a). However, the extent to which functional interactions within layers of V1 reflect these anatomical projections remains unclear. Thus, we examined the distribution of CCG classes across pairs of V1 input layers 4cα, 4cβ, and 6. We found that for the 4cα-4cα, 4cβ-4cβ, and 6–6 pairings, Ssync CCGs were observed much more frequently than other CCG classes (Figure 6b) (chi-squared test; 4cα-4cα, 4cβ-4cβ: P<10–5, 6–6: P<10–2). This overrepresentation of Ssync CCGs may reflect the fact that neurons in 4cα, 4cβ and 6 receive common and converging input from the dLGN. Furthermore, it is noteworthy that the Ssync class was overrepresented in V1 input layers, but the Bsync class was not.

Distribution of different putative classes of functional interactions within V1 input layers.

(a) Diagram of dLGN input to V1 layers 4cα, 4cβ, and 6; dLGN axons terminate in layers 4cα and 4cβ and layer 6, and layer 6 projects to layers 4cα and 4cβ. (b) Percentage of CCGs in each class composed of two neurons in layer pairings of 4cα-4cα, 4cβ-4cβ, or 6–6, out of all the pairwise layer pairing combinations. (c) Percentage of CCGs in each class composed of one neuron in layer 4cα and one in 4cβ (4cα-4cβ or 4cβ-4cα). In (b-c), the reference neuron in a neuronal pair was selected randomly. (d) Percentage of CCGs in each class composed of a reference neuron in layer 6 and a target neuron in 4cβ (left) or 4cα (right). For (b-d), error bars denote 95% confidence intervals. Large asterisks denote significant chi-squared test across all classes (p<0.05, Bonferonni corrected). Small asterisks denote significant chi-squared tests between a particular class and the class with the maximum percentage (p<0.05, Bonferroni corrected).

In contrast to the overrepresentation of Ssync CCGs within the input layers, this class of CCGs was not overrepresented in functional interactions between input layers. Of the four CCG classes, the proportions of each found among pairs composed of one neuron in layer 4cα and one neuron in layer 4cβ were statistically indistinguishable (Figure 6c) (chi-squared test; p=0.20). The lack of an overrepresentation of Ssync CCGs among 4cα/4cβ pairs could reflect the lack of synchrony between magnocellular (fast) and parvocellular (slow) inputs to V1. This result is noteworthy given that the average distance between neuronal pairs across 4cα and 4cβ was comparable to the distances between neuronal pairs within V1 input layer 6 (mean distance: 4cα/4cβ = 106 µm; 6/6=80 µm). In examining functional interactions between layers 4c and 6, we considered that a temporal offset between layer 6 and 4c neurons might exist given extensive projections from layer 6 pyramidal neurons to layer 4c (Wiser and Callaway, 1996). To test this, we examined the 6–4cα and 6-4cβ pairs in which the layer 6 neuron was the reference neuron in the crosscorrelation function. Indeed, in addition to observing that Ssync CCGs were poorly represented, we found that the Fasync class was significantly overrepresented in these pairs (Figure 6d) (6–4cα: p<10–5, 6-4cβ: p<10–3).

Lastly, we examined the distribution of CCG classes across pairs of neurons involving layer 2/3 neurons. A wealth of evidence indicates that layer 2/3 neurons provide a major source of output to other neocortical areas (reviewed in Callaway, 1998; Douglas and Martin, 2004; Felleman and Van Essen, 1991; Harris and Shepherd, 2015; Thomson and Lamy, 2007). In macaque V1, layer 2/3 neurons send projections to higher visual areas such as V2 (Livingstone and Hubel, 1984; Rockland, 1992; Sincich and Horton, 2005) and V4 (Yukie and Iwai, 1985), and receive inputs from all the deeper cortical layers, including layers 4cα, 4cβ, 4A, 4B, 5, and 6 (Blasdel et al., 1985; Callaway, 1998; Callaway and Wiser, 1996; Fitzpatrick et al., 1985; Kisvarday et al., 1989; Lachica et al., 1992; Lund and Boothe, 1975; Sawatari and Callaway, 2000; Vanni et al., 2020; Wiser and Callaway, 1996; Yarch et al., 2017; Yoshioka et al., 1994; Figure 7a). Consequently, one might predict that a predominant proportion of projections to 2/3 neurons from other layers might be forward ones (Callaway, 1998; Mejias et al., 2016; Schmidt et al., 2018). Consistent with this prediction, we found that the forward (Fasync) class was overrepresented among functional interactions from layers 6, 5, 4cβ, 4cα, and 4A/B to layer 2/3 (Figure 7b) (chi-squared test; 6: p<10–5; 5: p<10–5; 4cβ: p<10–3; 4cα: p<10–2; 4A/B: p<10–5). In contrast, functional interactions within layer 2/3 exhibited a very different pattern. Within the same layer, the classes of 2/3-2/3 CCGs were more evenly represented, in stark contrast to the pattern of within-layer CCGs observed in the input layers (Figure 6b). Within layer 2/3, the Ssync and Bsync CCGs were overrepresented among functional interactions (Figure 7c) (chi-squared test: p<10–4), and there was an equal representation of Ssync and Bsync CCGs among 2/3-2/3 pairings.

Distribution of different putative classes of functional interactions with layer 2/3.

(a) Diagram of input to layer 2/3 from V1 laminar compartments. (b) Percentage of CCGs in each class composed of a reference neuron in layer 6, 5, 4cβ, 4cα, or 4A/B and a target neuron in layer 2/3. (c), Percentage of CCGs in each class composed of pairs of neurons in layer 2/3. The reference neuron in a neuronal pair was selected randomly. For b-c, conventions are the same as in Figure 6.

Discussion

Using high-density recordings from single neocortical columns of macaque V1, we identified 1000s of functionally interacting neuronal pairs using established crosscorrelation approaches. The results demonstrate clear and systematic variations in the synchrony and strength of functional interactions within single V1 columns. Notably, we observed that in spite of residing within the same column, the functional interactions between pairs of V1 neurons depended heavily on their vertical distance within the column; both the peak lag and peak efficacy of CCGs between neuronal pairs changed dramatically within only a few hundred micrometers of vertical distance within the column. In addition, we found that the synchrony and strength of CCGs also depended on laminar location and the similarity of orientation tuning between neuronal pairs. We leveraged the statistical power provided by the large numbers of functionally interacting pairs to categorize functional interactions between neurons based on their crosscorrelation functions. These analyses identified distinct classes of functional interactions within the full population. Those distinct classes exhibited different distributions across defined laminar compartments, and those differences were consistent with known and/or expected properties of V1 cortical circuitry. The results demonstrate a novel utility of high-density neurophysiological recordings in assessing circuit-level interactions within local neuronal networks. Below, we discuss both the implications and the limitations of this approach.

Effect of cortical distance on functional interactions

A wealth of previous evidence has established a clear effect of cortical distance on functional interactions, yet a majority of past studies have focused on the effect of horizontal distance across which large changes in shared input between neurons are expected. Evidence that spiking correlations and synchrony decline with horizontal cortical distance within V1 has been shown in cats (Das and Gilbert, 1999; Gray et al., 1989; Hata et al., 1991) (but see Samonds et al., 2006; Schwarz and Bolz, 1991), monkeys (Chu et al., 2014; Krüger and Aiple, 1988; Maldonado et al., 2000; Smith and Kohn, 2008), and in mice (Denman and Contreras, 2014). Very few studies have examined crosscorrelations among pairs of neurons within a single column, where the feedforward input is largely shared (e.g. DeAngelis et al., 1999). Longer timescale, spike count (‘noise’), correlations, which have been widely assessed in studies of primate visual cortex (Averbeck et al., 2006) have been shown to be layer dependent within macaque V1 where weaker correlations occur in layer 4 (Hansen et al., 2012). However, no evidence that such correlations depend on distance independent of layer was observed. In contrast, measurement of crosscorrelations in earlier studies of V1 columns in cat indeed suggest that functional interactions are restricted to local regions across cortical depths (Toyama et al., 1981b). Within rat auditory cortex, functional interactions diminish dramatically within ~300 µm of vertical columnar distance (Atencio and Schreiner, 2013), similar to what is observed in rat somatosensory cortex (Khateb et al., 2021). A dependence of interactions on vertical distance is further supported by evidence from multiple whole-cell recordings in mouse visual cortex which demonstrates that connection probability decreases sharply within a distance of 250 µm (Jiang et al., 2015). Our observation that the peak efficacy of CCGs was greatly diminished within <200 µm within macaque V1 is thus consistent with estimates from other sensory cortices and species.

Effect of orientation tuning similarity on functional interactions

Many previous studies have reported an effect of tuning similarity on functional interactions, yet most of these studies focused on interactions between pairs of neurons in different cortical columns. Evidence from horizontal array recordings in macaque V1 suggests that pairs of neurons selective for similar orientations (Chu et al., 2014; Kohn and Smith, 2005; Smith and Kohn, 2008), or color and luminance (Chu et al., 2014) tend to exhibit stronger functional interactions compared to pairs with dissimilar tuning, perhaps reflecting the functional specificity of horizontal connections (Gilbert and Wiesel, 1989). For example, one study found that the strength of functional interactions between neuronal pairs with the highest orientation signal correlations was nearly twice that of uncorrelated pairs (Smith and Kohn, 2008). Consistent with these previous studies, we also found that variations in tuning similarity resulted in robust differences in the synchrony and strength of correlated activity. However, our results demonstrate that an interdependence of tuning similarity and functional interaction exists even within single orientation columns. In addition, in contrast to comparisons across columns where tuning similarity appears to be the primary factor (Chu et al., 2014), we found that within single columns, vertical distance and orientation tuning similarity exerted distinct effects on the synchrony and strength of functional interactions.

Distinct classes of functional interactions

Upon clustering the full population of significant CCGs, we identified four putative classes of functional interactions. Notably, these classes of CCGs depicted the set of pairwise interactions that one might logically expect, namely two asynchronous classes (forward and reverse) and two synchronous ones (sharply and broadly synchronous). More importantly, the clustering-based evidence of distinct classes of CCGs was corroborated by the observation of highly differential distributions of those putative classes across cortical layers. For example, asynchronous CCG classes were more often observed among neuronal pairs within different layers, whereas synchronous pairs more often resided within the same layer. This corroboration of distinct classes extended to functional properties of V1 neurons as well in that simple cells were more often paired in a forward manner with complex cells, whereas the reverse was true for complex cells. Nevertheless, the existence of exactly four distinct classes among V1 pairs is by no means certain. Indeed, the choice of three classes was almost as valid as that of four in the clustering procedure (Figure 4). Yet, given the clear evidence of two asynchronous classes, and the differential distributions of broadly and sharply synchronous pairs across cortical layers, the choice of only three classes of CCGs seems less parsimonious than four. Although there appeared to be less evidence for the existence of five or more classes of functional interactions, that possibility cannot be ruled out either. For example, additional distinct classes of CCGs might be present, but significantly less frequent or weaker than the other four. Indeed, given their low incidence, our selection criteria already excluded CCGs with significant inhibitory peaks. As in previous studies, the frequency of excitatory CCGs in our dataset was considerably higher than that of inhibitory CCGs (Aertsen and Gerstein, 1985; Hembrook-Short et al., 2019; Table 2). Consequently, these additional classes were eliminated by the statistical threshold employed to identify significant CCGs. It is likely that additional distinct classes, excitatory or inhibitory, were also eliminated and/or simply fell within a mixture of the more dominant four classes that exceeded the statistical criterion. Future work will therefore be needed to more extensively characterize the distribution of distinct classes of spiking crosscorrelations among neuronal pairs in cortical columns.

Table 2
Number of CCGs with peaks or troughs significantly above or below noise.
Stdev. above/below noiseNumber of pairs with peaks above noiseNumber of pairs with troughs below noise
133,50215,582
233,50211,347
331,6604,537
425,1801,555
518,739603
613,757265
710,246136
  1. The number of recorded pairs with a peak or trough at least 1–7 standard deviations (SD) above or below noise is shown. Only CCGs with peaks or troughs within 10ms of zero time lag were considered. 136 CCGs had troughs that were more than 7 SD below noise whereas 10,246 CCGs had peaks that were more than 7 SD above noise.

Classes of functional interactions and V1 microcircuitry

We found that the four putative classes of functional interactions were differentially distributed across the cortical column and across functional pairs of neurons. Most notably, the different CCG classes were observed in different proportions across V1 layers. In spite of those differences, it need not follow that the relative distribution of any specific putative class (e.g. sharply synchronous) fits with the known (or predicted) connectivity between different V1 neurons. For example, our observation that asynchronous CCGs (Fasync and Rasync) were considerably more frequent among neuronal pairs situated in different laminae, and that synchronous CCGs (Ssync and Bsync) were found among neurons in the same laminae, would not be expected if the differences in synchrony resulted primarily from measurement noise. Likewise, the observed disproportionality of Fasync and Rasync CCGs among functional interactions between simple and complex cells would not be expected if the two asynchronous classes were indistinguishable in our measurements. Instead, not only were the two classes disproportionate among simple and complex cells, but the overall direction of disproportionality was consistent with the known connectivity between the two functional classes of cells (Alonso and Martinez, 1998; Martinez and Alonso, 2001; Yu and Ferster, 2013). Thus, overall, we found that the pattern of differential distributions of the putative classes of CCGs across the column and across functional pairs of neurons was largely consistent with known properties of V1 microcircuitry.

Functional interactions and synaptic connectivity

Analyses of the statistical dependencies between spike trains of two or more neurons have long played an important role in estimating how ensembles of neurons interact with one another (Casile et al., 2021; Okatan et al., 2005; Perkel et al., 1967). We interpret the functional interactions identified here in the same manner. Using statistical criteria employed in a recent study of mouse visual cortex (Siegle et al., 2021), we found that ~15% (10,246/68,579) of neuronal pairs within columns of macaque V1 exhibited significant functional interactions. The proportion of significant pairs is within the range observed in previous studies in macaque V1 (Chu et al., 2014; Hembrook-Short et al., 2019; Kohn and Smith, 2005; Smith and Kohn, 2008). This similarity with previous studies exists in spite of notable differences in the electrophysiological approach and statistical criteria. In addition, the proportion of significant pairs is also similar to previous CCG measurements made in the cat (Alonso and Martinez, 1998) and mouse V1 (Denman and Contreras, 2014; Siegle et al., 2021).

In addition, it has long been understood that the statistical interdependence of spike trains is not necessarily an indicator of synaptic connectivity among neurons (Ostojic et al., 2009; Perkel et al., 1967). Indeed, spike crosscorrelations are in no way a substitute for more direct measurements of synaptic connectivity, for example multi-patch recordings (Cadwell et al., 2020; Jiang et al., 2015; Song et al., 2005) or electron microscopy (Bock et al., 2011; Lee et al., 2016). Notably, only a subset of the identified CCG classes here are consistent with, although not necessarily indicative of, monosynaptically connected neuronal pairs. Specifically, pairs with non-zero time lags, such as the asynchronous (Fasync and Rasync) classes we describe are those most likely to be monosynaptically connected. In our data, the asynchronous class amounted to ~9% of the total pairs tested. In contrast, the synchronous classes, particularly the Ssync class, which are more likely to reflect pairs with common input (Ostojic et al., 2009; Perkel et al., 1967), amounted to ~6% of the total pairs tested. However, for the former group, it is uncertain of course how many of the total neuronal pairs represent genuine monosynaptic connections. Past studies have employed widely varying criteria to identify monosynaptic connections in a number of different brain structures and species (Alonso and Martinez, 1998; English et al., 2017; Hembrook-Short et al., 2019; Liew et al., 2021; Reid and Alonso, 1995; Senzai et al., 2019). Those criteria varied in the stringency with which CCGs are labeled as monosynaptic and include criteria that necessarily yield false negatives (Senzai et al., 2019).

Unfortunately, there exists little or no ground-truth measurements of the rate of synaptic connectivity, or common input, between neurons within single V1 columns in any species, though such measurements may be imminent (https://www.microns-explorer.org/). Nonetheless, the observed rate of asynchronous classes identified in our data is consistent with estimates of the connection probability among V1 neurons obtained in multi-patch, slice recordings. For example, recent studies identified synaptic connectivity in 5–12% of excitatory neurons in mouse/rat V1 (Cadwell et al., 2020; Song et al., 2005). Although estimates from such studies clearly underestimate the rate of connections given that many connections are cut in the slice preparation, they nonetheless provide a plausible lower bound. Thus, in spite of the inherent arbitrariness of the statistical criteria employed in identifying significant CCGs, the observed proportion of asynchronous interactions among pairs in our data appears largely consistent with previous estimates. As a means of comparing the criteria used here to those of prior studies in monkey V1, we applied one criterion used in Hembrook-Short et al., 2019 to our data, which required that CCG peaks be sharp and narrow (<5ms) to be considered monosynaptic. Using this criterion, we found that 86% of our asynchronous CCGs were categorized as consistent with monosynaptic connections. Thus, it seems likely that our asynchronous (Fasync and Rasync) classes captured a reasonable proportion of monosynaptically connected neurons within the V1 column.

Future studies

We found that the relative instances of different types of crosscorrelations observed among large populations of neuronal pairs may provide a means of constraining models of cortical microcircuits. This approach could prove particularly valuable in less experimentally tractable model systems such as nonhuman primates, or perhaps even in the human brain, where direct interrogation of microcircuits is difficult or not yet possible. In such cases, the ability of high-channel count, high-density, probes to dramatically increase the number of identifiable functional interactions within a local network of neurons is among their greater benefits. Our results thus far suggest that this approach works well and could be extended to examine higher-order interactions among larger sets of neurons, and to identify neuronal ensembles with distinct functional properties (Fujisawa et al., 2008; Miller et al., 2014; See et al., 2018). In addition, future studies should be able to compare local interactions across different putative cell types estimated from their spike waveforms (Johnston et al., 2009; Lee et al., 2021; Mitchell et al., 2007; Wilson et al., 1994) and/or spiking patterns (Onorato et al., 2020). Combined with measurements of functional interactions, such an approach could be used to constrain models of microcircuit architecture from neurophysiological data obtained from any number of uniquely evolved primate brain structures.

Methods

Experimental model and subject details

Anesthetized recordings were conducted in 2 adult male rhesus macaques (Macaca Mulatta, M1, 13 kg; M2, 8 kg). All experimental procedures were in accordance with National Institutes of Health Guide for the Care and Use of Laboratory Animals, the Society for Neuroscience Guidelines and Policies, and with approved Institutional Animal Care and Use Committee (IACUC) protocol (#APLAC-9900) of Stanford University.

Electrophysiological recordings

Prior to each recording session, treatment with dexamethasone phosphate (2 mg per 24 hr) was instituted 24 hr to reduce cerebral edema. After administration of ketamine HCl (10 mg/kg body weight, intramuscularly), monkeys were ventilated with 0.5% isoflurane in a 1:1 mixture of N2O and O2 to maintain general anesthesia. Electrocardiogram, respiratory rate, body temperature, blood oxygenation, end-tidal CO2, urine output and inspired/expired concentrations of anesthetic gases were monitored continuously. Normal saline was given intravenously at a variable rate to maintain adequate urine output. After a cycloplegic agent was administered, the eyes were focused with contact lenses on a CRT monitor. Vecuronium bromide (60 μg/kg/hr) was infused to prevent eye movements.

With the anesthetized monkey in the stereotaxic frame, an occipital craniotomy was performed over the opercular surface of V1. The dura was reflected to expose a small (~3 mm2) patch of cortex. Next, a region relatively devoid of large surface vessels was selected for implantation, and the Neuropixels probe was inserted with the aid of a surgical microscope. Given the width of the probe (70 μm x 20 μm), insertion of it into the cortex sometimes required multiple attempts if it flexed upon contacting the pia. The junction of the probe tip and the pia could be visualized via the (Zeiss) surgical scope and the relaxation of pia dimpling was used to indicate penetration, after which the probe was lowered at least 3–4 mm. Prior to probe insertion, it was dipped in a solution of the DiI derivative FM1-43FX (Molecular Probes, Inc) for subsequent histological visualization of the electrode track.

Given the length of the probe (1 cm), and the complete distribution of electrode contacts throughout its length, recordings could be made either in the opercular surface cortex (M1) or within the underlying calcarine sulcus (M2), by selecting a subset of contiguous set of active contacts (n=384) from the total number (n=986). Receptive fields (RFs) from online multi-unit activity were localized on the display using at least one eye. RF eccentricities were ~4–6° (M1) and ~6–10° (M2). Recordings were made at 1–3 sites in one hemisphere of each monkey. At the end of the experiment, monkeys were euthanized with pentobarbital (150 mg kg−1) and perfused with normal saline followed by 1 liter of 1% (wt/vol) paraformaldehyde in 0.1 M phosphate buffer, pH 7.4.

Visual stimulation

Visual stimuli were presented on a LCD monitor NEC-4010 (Dimensions = 88.5 (H)* 49.7 (V) cm, pixels = 1360 * 768, frame rate = 60 Hz) positioned 114 cm from the monkey. Stimuli consisted of circular drifting Gabor gratings (2 deg./sec., 100% Michelson contrast) positioned within the joint RFs of recorded neurons monocularly. Gratings drifted in 36 different directions between 0–360° in 10° steps in a pseudorandom order. Four spatial frequencies (0.5, 1, 2, 4 cycle/deg.) were tested and optimal SFs were determined offline to categorize V1 neurons into simple or complex cell. The stimulus in each condition was presented for 1 second and repeated 5 or 10 times. A blank screen with equal luminance to the Gabor patch was presented for 0.25 s during the stimulus interval.

Data acquisition and spike sorting

Raw spike-band data was sampled and recorded at 30 kHz. It was then median-subtracted and high-pass filtered at 300 Hz during the pre-processing stage. Spike-sorting was carried out with Kilosort2 (https://github.com/MouseLand/Kilosort; Pachitariu et al., 2022) to find spike times and assign each spike to different units. The raw sorted data was then manually curated in Phy (https://github.com/cortex-lab/phy; Rossant et al., 2022) to remove spikes with atypical waveforms and perform minimal merging and splitting. One potential issue with the template-matching approach used by Kilosort2 is that the algorithm will occasionally fit a new template from the residual after subtracting the first template from the original data. This artificial template/neuron will share an abnormal number of double-counted spikes with the real neuron, resulting in a high zero-time lag synchrony between those two neurons. To examine whether this issue may affect our results, we used a criteria suggested by a previous study (Siegle et al., 2021) to identify double-counted spikes by counting spikes with peak times within 5 samples (0.167ms) and from pairs of neurons within 50 μm (~5 channels). Out of the total 68,579 pairs of neurons included in this current study, we found that only 7 pairs shared more than 20% overlapping spikes. Considering that we identified 1408 pairs of neurons in the sharply synchronous class, those potentially artificial pairs contributed to less than 0.5% of this class. Although we did not remove those pairs from the data, we believe their contributions are negligible. Here, we list key parameters in Kilosort2 that may affect the ‘double-counting’: Ops.th=[10,4]; Ops.lam=20; Ops.AUCsplit=0.9; Ops.ThPre=8; Ops.spkTh=-6. Moreover, only neurons with a minimum firing rate of 3 spikes/s were included in the study.

Layer assignment

The laminar location of our recording sites was estimated based on a combination of functional analysis and histology results. For each recording, we first performed the current source density (CSD) analysis on the stimulus-triggered average of local field potentials (LFP). LFP were low-pass filtered at 200 Hz and recorded at 2500 Hz. LFP signals recorded from each four neighboring channels were averaged and realigned to the onset of visual stimulus. CSD was estimated as the second-order derivatives of signals along the probe axis using the common five-point formula (Nicholson and Freeman, 1975). The result was then smoothed across space (σ=120 μm) to reduce the artifact caused by varied electrode impedance. We located the lower boundary of the major sink (the reversal point of sink and source) as the border between layer 4c and layer 5/6. Based on this anchor point, we assign other laminar compartment borders using the histological estimates.

Single neuron properties

To characterize the visual properties of each neuron, the stimulus evoked activity was assessed using mean firing rate (spikes/s) over the entire stimulus presentation period, offset by a response latency of 30ms. Only responses to the preferred spatial frequency were used. Modulation ratio was defined as F1/F0, where F1 and F0 are the amplitude of the first harmonic at the temporal frequency of drifting grating and constant component of the Fourier spectrum to the neuron’s response to preferred orientation. Simple cells were defined as cells with modulation ratio larger than 1, and complex cells have modulation ratios smaller than 1 (De Valois et al., 1982; Skottun et al., 1991).

Signal correlations

To measure the similarity of orientation tuning between neuronal pairs, we computed an orientation signal correlation (rori). The orientation signal correlation was defined as the Pearson’s correlation coefficient between the mean responses of two neurons to each of the 36 stimulus orientations (Smith and Kohn, 2008). For each neuron and orientation, a single mean response was computed by averaging spiking activity over the entire duration of stimulus presentation (1 s) across all trials with a particular orientation.

Cross-correlograms (CCGs)

To measure correlated firing, we computed the crosscorrelation between spike trains of all pairs of simultaneously recorded neurons (Jia et al., 2013; Siegle et al., 2021; Smith and Kohn, 2008; Zandvakili and Kohn, 2015). We focused on the spiking activity within the 0.4–1 s window of each visual stimulus presentation, which ensured that the analysis was not affected by the transient response to stimulus onset. To mitigate firing rate effects, we normalized the cross-correlation for each pair of neurons by the geometric mean of their firing rates. Thus, the cross-correlogram CCG for a pair of neurons (j,k)was defined as follows:

CCG(τ)jk=1Mi=1Mt=1Nτxji(t)×xki(t+τ)θ(τ)λjλk

where M is the number of trials, N is the number of time bins within a trial, τ is the time lag, xjit is one if neuron j fired in time bin t of trial i and zero otherwise, and λj is the mean firing rate of neuron j computed over the same bins used to compute the CCG at each time lag. θτ is a triangular function, θ(τ)=N|τ|, that corrects for the difference in the number of overlapping bins at different time lags. We denote the CCG computed with neuron j as the first (reference) neuron and k as the second (target) neuron in the correlation function as j-k.

To correct for correlation due to stimulus-locking or slow fluctuations in population response (e.g. gamma-band activity), we computed a jitter-corrected cross-correlogram by subtracting a jittered cross-correlogram from the original cross-correlogram:

CCGcorrected=CCGoriginal-CCGjittered

The jittered cross-correlogram (CCGjittered) reflects the expected value of cross-correlograms computed from all possible jitters of each spike train within a given jitter window (Harrison and Geman, 2009; Smith and Kohn, 2008). The jittered spike train preserves both the PSTH of the original spike train across trials and the spike count in the jitter window within each trial. As a result, jitter correction removes the correlation between PSTHs (stimulus-locking) and correlation on timescales longer than the jitter window (slow population correlations). Here, a 25 ms jitter window was chosen based on previous studies (Jia et al., 2013; Siegle et al., 2021; Zandvakili and Kohn, 2015).

We classified a CCG as significant if the peak of the jitter-corrected CCG occurred within 10ms of zero and was more than seven standard deviations above the mean of the noise distribution. The noise distribution for a CCG was defined as the flanks of the jittered-corrected CCG ({CCG(τ)100|τ|50}). This significance criterion was chosen based on that of Siegle et al., 2021. All analyses presented here involve only significant, jitter-corrected cross-correlograms. Note that the criterion identifies only positive peaks in the CCG and excludes significant inhibitory correlations. However, consistent with earlier studies (Aertsen and Gerstein, 1985; Hembrook-Short et al., 2019), we found that the frequency of CCGs with significant troughs was approximately 40x lower than those with significant peaks (Table 2).

Classification of Cross-correlogram

To identify distinct classes of cross-correlation functions, we clustered significant crosscorrelations. We only analyzed crosscorrelation functions between τ=-10 and τ=10 such that our input CCGs had 21 features, corresponding to the 21 crosscorrelation values between τ=-10 and τ=10. For clustering, we included two crosscorrelation functions for each pair of neurons j,k , one computed using the above CCG function with neuron j as the reference neuron j-k and the other with neuron k as the reference neuron k-j . This was done in order to avoid introducing biases in the direction of the CCG templates. We z-scored each CCG prior to clustering to encourage clustering based on the shape of the correlation function rather than its magnitude. For subsequent statistical analyses, only a single neuron in the pair was used as the reference.

To simplify the clustering problem, we used t-distributed stochastic neighbor embedding (t-SNE) to reduce our input data with 21 features to 3 features (tsne, MATLAB R2019a). t-SNE was used instead of principal component analysis (PCA) because it is more robust to outliers since it captures neighbor relationships in the input space. We clustered the dimensionality-reduced data using k-means with k=1 to 10 (50 replicates, 100 max iterations, kmeans MATLAB R2019a). To determine the optimal number of clusters, we used two complementary approaches, the elbow method and silhouette method. The elbow method selects k based on the magnitude of the change in the variance explained by clustering as k increases. For a set of points S=s1,s2,,sn divided into k clusters S1,S2,,Sk , the percent of variance explained by clustering ηk is:

ηk=TSS-WCSSkTSS
TSS=xS||xmean(S)||
WCSSk=i=1kxSi||xmean(Si)||

where TSS denotes the total sums of squares and WCSSk denotes the sum of within cluster sums of squares over all clusters. The optimal number of clusters occurs at the point where the percent of explained variance plateaus (or ‘elbows’) as the number of clusters increases. The silhouette criterion captures how similar a point is to its own cluster versus how different it is from the nearest cluster that it is not a member of. We computed the silhouette criterion using MATLAB’s ‘evalclusters’ function with default parameters (MATLAB R2019a).

Statistical analyses

The effects of vertical pair distance and orientation signal correlation on CCG peak lag and peak efficacy were fit using linear and exponential functions. In linear regressions predicting CCG peak lag, all significant CCGs were included, and mean squared error was used as the cost function for regressions. In linear and exponential regressions predicting CCG peak efficacy, only significant CCGs with non-outlier peaks (1.5*IQR criterion) were included, and mean absolute error was used as the cost function for regressions to encourage fit of the plotted median peak efficacies.

The relationships between classes of functional interactions and signal correlation/pair distance were evaluated using Wilcoxon rank-sum tests, and the relationship between functional class and layer/cell type pairings was assessed using chi-squared tests or one proportion z-tests. Finally, the dependence of functional class on whether a CCG was composed of two neurons within the same or different cortical layer(s) with comparable vertical distance was assessed using logistic regression.

Distance matching

Distance matching was used to compare the effects of orientation signal correlation on CCG peak lag and peak efficacy among neuronal pairs with comparable cortical distances. To match pairs with comparable distances (Figure 3C), we sorted significant CCGs by cortical distance, then paired the CCGs with the smallest and second smallest distances and paired the CCGs with the third and fourth smallest distances and so forth. Thus, every significant CCG was paired with exactly one other significant CCG, resulting in 5122 pairs. To verify that this procedure effectively matched pairs of CCGs with comparable cortical distance, we examined the difference in cortical distance for distance-matched pairs. More than 99% (5067/5122) of the distance-matched pairs had a difference in cortical distance of less than 2 µm. Finally, we examined the correlation between the difference in CCG peak lag or peak efficacy and difference in signal correlation for matched pairs to determine whether signal correlation predicts peak lag or peak efficacy when controlling for distance.

Data availability

All the raw data generated as part of this study are publicly accessible. All the raw code generated for analyzing the data has already been deposited to GitHub and is currently freely accessible (https://github.com/et22/functional_connections_macaque_v1, copy archived at swh:1:rev:39dda9dffab8f6e54069fb6514d2230200412472).

The following data sets were generated
    1. Trepka E
    2. Zhu S
    3. Xia R
    4. Chen X
    5. Moore T
    (2022) Dryad Digital Repository
    Functional Interactions Among Neurons within Single Columns of Macaque V1.
    https://doi.org/10.5061/dryad.x3ffbg7p2

References

    1. Horton JC
    2. Adams DL
    (2005) The cortical column: a structure without a function
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 360:837–862.
    https://doi.org/10.1098/rstb.2005.1623
    1. Lorente de No R
    (1938)
    Physiology of the Nervous System
    291–339, Cerebral cortex: architecture, intracortical connections, motor projections, Physiology of the Nervous System, London, Oxford University Press.

Decision letter

  1. Martin Vinck
    Reviewing Editor; Ernst Strüngmann Institute (ESI) for Neuroscience in Cooperation with Max Planck Society, Germany
  2. Timothy E Behrens
    Senior Editor; University of Oxford, United Kingdom
  3. Jens Kremkow
    Reviewer; Charite-Universitätzmedizin Berlin, Germany

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Functional Connections Among Neurons within Single Columns of Macaque V1" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Timothy Behrens as the Senior Editor. The following individual involved in the review of your submission has agreed to reveal their identity: Jens Kremkow (Reviewer #1).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1. Although the paper does provide rich information on interactions within local cortical circuits, the main weakness of the paper is using the term "functional connection" in an imprecise manner. Cross-correlograms (CCG) of spike trains of pairs of neurons show different shapes depending on the underlying connectivity and not all significant peaks in CCGs reflect functionally connected neuron pairs. For example, CCGs of synaptically connected neuron pairs show a transient peak that is offset from the 0-ms lag due to the synaptic delay. CCGs with this shape thus reflect "functionally connected neuron pairs". In contrast, common inputs to pairs of neurons can induce significant peaks in CCGs, despite the fact that these neurons are only correlated but not functionally connected (e.g. Ostojic et al. 2009). Therefore, taking the shape of significant CCGs into account is important when discussing "functionally connected neuron pairs". While the authors mention this point in the paper, the term "functional connection" is nonetheless used irrespective of the CCG shapes which can be confusing to the reader. Moreover, the authors claim that the method allows identifying "1000s of functionally connected neuronal pairs". This statement is likely not fully supported by the data, evident by the fact that CCGs with the shape of mono-synaptic connections (transient and non-zero lag peak) are not among the distinct classes of CCGs shown in Figure 4.

The term "functionally connection" implies a synaptic connection between a pair of neurons. However, the majority of CCGs shown in the paper likely reflect correlated activity at fast timescales rather than true functionally connections. Because this can be confusing to readers I recommend defining the terminology more precisely and also using the term "functionally connected neuron pairs" only in cases where this is justified. This is particularly important because one claim of the paper is that Neuropixels probes allow identifying "1000s of functionally connected neuronal pairs". Please either show that indeed this large number of synaptically connected neurons is in the dataset, or change the title/abstract to match the results.

2. Likewise, it is surprising that CCGs that reflect mono-synaptic connections are not among the distinct classes shown in Figure 4? Why is that? Given the large number of tested interactions, and the claims of the paper, we would have expected that mono-synaptic connections form one of the distinct classes. It could be interesting and useful to specifically identify and study mono-synaptic connections, e.g. by employing methods reported in Liew et al. 2021 JNP, or by other approaches.

3. It is well established that fast-spiking neurons in the cortex receive stronger inputs from neurons in the local circuitry as compared to regular spiking neurons. Because fast-spiking and regular spiking neurons can be distinguished based on their spike waveform in extracellular recordings it could be interesting, but not required, to see whether this cell-type dependent connection strength is also evident in this dataset. This could add to the significance of the work and provide another angle to investigate circuit interactions. Again, this is only a suggestion.

4. It would be useful to show the CCGs in Figure 1e in a higher temporal resolution around the peak such that the lag of the peak is visible.

5. It seems odd to treat the forward and reverse correlation functions as two distinct types since the labeling is arbitrary. That is, swapping the label for the reference neuron and its partner would flip the correlation function from one category to the other. Since the choice of the reference neuron is arbitrary, it is not clear why there are two categories. Please provide justification.

6. Separate from this issue, the labeling looks erroneous to the best of our knowledge. That is, 'forward' correlation functions have more mass at negative lags. This is opposite to the conventional definition. Mass at negative lags means the partner neuron tends to fire before the reference neuron, which is not a forward connection (reference driving partner) but the opposite. Please clarify.

7. Please provide the species used.

8. Page 10: It was not obvious that the regression coefficients could be directly compared since the covariates (pair distance and r_ori) have very different magnitudes. Please clarify.

9. Please provide an explanation of how the spike sorting was done. One concern about the sharp correlation peaks is that they are artifactual, resulting from issues in Kilosort (https://github.com/MouseLand/Kilosort/issues/29). In brief, Kilosort can match the same spike waveform to two different templates, one template capturing most of the waveform shape and the second capturing the residual. This issue is discussed on the Github Kilosort page. Are the authors sure that this does not contribute to/drive the sharp peak phenomenology? One way to provide reassurance would be to provide more information on the shape and width of the sharp peaks. If they are ~1 ms wide, it seems more likely to be artifactual than biological.

10. Please provide the proportion of simple and complex neurons in the data set.

11. Discussion: The authors do a good job of discussing the limitations of correlation analysis for identifying synaptic connections. However, they refer to their R and F shapes (9% of cases) as potentially capturing such connections. Traditionally (e.g. Reid and Alonso), correlation function consistent with synaptic connections must pass a more stringent description than a simple asymmetry with respect to zero time lag. Namely, they need to have a sharp peak clearly offset from zero. This difference in stringency should be incorporated into the Discussion if the authors wish to propose that R and F types are related to connectivity.

12. Line 635: The definition of the correlation function seems different from previous work. The correlation function appears normalized by a quantity that depends on the time lag of one neuron and not the other (the denominator). Please clarify the notation.

13. Figure 2B: There is an oscillation in the marginal histogram of peak delays. This is also evident in Figure 3C. Please clarify.

https://doi.org/10.7554/eLife.79322.sa1

Author response

Essential revisions:

1. Although the paper does provide rich information on interactions within local cortical circuits, the main weakness of the paper is using the term "functional connection" in an imprecise manner. Cross-correlograms (CCG) of spike trains of pairs of neurons show different shapes depending on the underlying connectivity and not all significant peaks in CCGs reflect functionally connected neuron pairs. For example, CCGs of synaptically connected neuron pairs show a transient peak that is offset from the 0-ms lag due to the synaptic delay. CCGs with this shape thus reflect "functionally connected neuron pairs". In contrast, common inputs to pairs of neurons can induce significant peaks in CCGs, despite the fact that these neurons are only correlated but not functionally connected (e.g. Ostojic et al. 2009). Therefore, taking the shape of significant CCGs into account is important when discussing "functionally connected neuron pairs". While the authors mention this point in the paper, the term "functional connection" is nonetheless used irrespective of the CCG shapes which can be confusing to the reader. Moreover, the authors claim that the method allows identifying "1000s of functionally connected neuronal pairs". This statement is likely not fully supported by the data, evident by the fact that CCGs with the shape of mono-synaptic connections (transient and non-zero lag peak) are not among the distinct classes of CCGs shown in Figure 4.

The term "functionally connection" implies a synaptic connection between a pair of neurons. However, the majority of CCGs shown in the paper likely reflect correlated activity at fast timescales rather than true functionally connections. Because this can be confusing to readers I recommend defining the terminology more precisely and also using the term "functionally connected neuron pairs" only in cases where this is justified. This is particularly important because one claim of the paper is that Neuropixels probes allow identifying "1000s of functionally connected neuronal pairs". Please either show that indeed this large number of synaptically connected neurons is in the dataset, or change the title/abstract to match the results.

We thank the reviewer for raising this point as we agree that it is important to use the correct terminology. In particular, as with any terminology, we think it is crucial to use terms as they are used in the literature. As we did not invent the term “functional connectivity”, it is thus important to consider how it has been used previously. The reviewer states that not all significant CCGs reflect “functionally connected” pairs thus implying that the phrase “functionally connected” must be synonymous with “synaptically connected”, or “mono-synaptically connected”. However, this is not how the term has generally been used in the literature. For example, Aertsen and Gerstein (1985) noted that: “The commonly used framework for investigating connectivity in a neural network is the cross-correlation function of simultaneously recorded spike trains. The presence of significant departures from a flat background is interpreted as indicative of a functional connection, either direct (through one or more synapses) or indirect (shared input, either neural or stimulus-induced)”. (Gerstein of course was one of the pioneers of the use and interpretation of CCGs.) Similarly, Friston and colleagues define “functional connectivity” as “statistical dependencies among remote neurophysiological events” (Friston, 1994, 2011).

In addition, many other studies have used the term “functional connectivity” similarly to the way we have and have treated CCGs consistent with common input as a type of “functional connection”. These studies include: Toyama, 1988; Schwarz and Bolz, 1991; Menz and Freeman, 2004; Denman and Contreras, 2014; Dann et al., 2016. In Menz and Freeman (2004), the authors define three types of functional connections:

“monosynaptic connection”, “polysynaptic connection” and “common input connection”. Denman and Contreras (2014), used the term “functional connectivity” throughout their paper and differentiated two types of functional connectivity: zeros-spanning and offset synchrony, similar to our “synchronous” and “asynchronous” categories. Aertsen and colleagues (1989) use the term “effective connectivity” and explicitly argue that “…The notion of ‘effective connectivity’ usually can be split into two types of circuitry: direct interactions and shared input.” A similar argument is made by Constantinidis et al., (2001).

In the above sense, “functional connectivity” is statistical in nature and refers to the temporal correlations between time series from different neuronal spike trains. This sense of the term avoids conclusions about the presence or absence of an underlying causal effect, which strictly speaking cannot be inferred from any CCG as it is not a ground truth measure. In the logical and grammatical sense, the word “functional” in “functional connectivity” simply qualifies the word “connection”, since by itself, the word “connection” must mean a synaptic one. In other words, in the literature, the term was used by investigators to distinguish statistically determined “connectivity” from synaptic connectivity, the latter of which must be directly determined by other means. Thus, our use of the term “functional connectivity” neither contradicts prior uses and definitions in the literature nor is it illogical or imprecise, and we are by no means the first to include CCGs consistent with common input as a type of “functional connection”. Indeed, this latter type of CCG is particularly important in our study as it appears to validate the expectation of common LGN input to neurons within particular thalamorecipient compartments (Figure. 6a-b).

There are indeed some prior studies that use the term “functional connectivity” to imply (mono)synaptic connections, but they are relatively few in number. For example, the term has been used this way to describe CCGs between simple and complex cells (Alonso and Martinez, 1998), thalamocortical connections (Bereshpolova, Hei, Alonso and Swadlow, 2020) and retinocollicular connections (Sibille et al., 2021). However, notably, the term itself is not clearly defined anywhere in these papers, and so the precise meaning, and its correspondence with uses in other studies, is left ambiguous. The reviewer cites Ostojic et al., 2009 to argue that neuronal pairs with significant CCGs induced by common input are only correlated but cannot be called “functionally connected”. However, Ostojic et al. do not use or define the term “functional connection” or “functional connectivity” at all in their paper. Instead, they use the term “functional interaction” and state in their introduction that “…statistically significant cross-correlations arise from the presence of a direct synaptic connection and/or from common or correlated inputs to the two neurons”.

To summarize, we agree that accurate use of terminology is important. Indeed, our use of “functional connectivity” is consistent with its use in the prior literature. Moreover, in the original manuscript, we had explicitly defined our use of the term in the discussion, and clarified which subset of identified CCG classes might be consistent with mono-synaptically connected pairs. Nevertheless, we also agree that it may be more important to provide the clearest possible use of terminology, given that many readers will be less familiar with prior terms used in the literature. Thus, we have expanded the introduction of the revised paper to both summarize past interpretations of CCGs and the terminology used in the literature and to clarify upfront our own interpretations and terminology to minimize confusion. In addition, to further avoid confusion regarding our interpretation of CCGs generally, we have adopted the term “functional interaction” as an alternative, since it avoids the word “connectivity” altogether. In the revised paper, we address the extent to which particular subclasses of significant CCGs are, or are not, consistent with monosynaptic connections in their respective instances.

2. Likewise, it is surprising that CCGs that reflect mono-synaptic connections are not among the distinct classes shown in Figure 4? Why is that? Given the large number of tested interactions, and the claims of the paper, we would have expected that mono-synaptic connections form one of the distinct classes. It could be interesting and useful to specifically identify and study mono-synaptic connections, e.g. by employing methods reported in Liew et al. 2021 JNP, or by other approaches.

First, it is important to emphasize that the distinct classes shown in Figure 4 are templates of CCG clusters (i.e. not individual CCGs) and that it is not necessarily the case that we should expect any of those templates to be monosynaptic-looking. Correlations between pairs of neurons within cortical columns likely reflect more than a single form of interaction, e.g. monosynaptic connection plus common inputs, or mono and disynaptic interactions, etc. Thus, when clustering CCGs, the dominant shapes of the correlations captured in the CCG templates (from large numbers of individual CCGs) likely reflect those mixtures. Unlike simple circuits, such as thalamocortical ones where one expects a single dominant, monosynaptic interaction between neuronal pairs (i.e., thalamic principal cell to layer IV cell), intracortical circuits are considerably more complex, and thus it may be overly simplistic to assume that one dominant, and singularly monosynaptic-looking CCG shape should emerge from the clustering.

Second, the presumption that mono-synaptic-looking connections are not part of our asynchronous cluster is wrong. When we apply the Liew et al. criteria (peak lag of 1-4 ms, peak > 3.5 s.d. trial-shuffled noise) to the asynchronous CCGs, those considered to be monosynaptic connections comprise 63% of the total group. Moreover, when we apply an alternative criterion to capture CCGs with sharp peaks (Hembrook-Short et al.) to the asynchronous CCGs, those considered to have sharp peaks consistent with monosynaptic connectivity comprise 86% of the total group. Thus, if we assume that the applied criteria are correct, a majority of the asynchronous class of CCGs may consist of a single, monosynaptic interaction between neuronal pairs. Note also that if we average across a large set of monosynaptic CCGs with varying peak lags (i.e. 1-4ms), the resultant template will necessarily have a broader peak than individual CCGs. Moreover, it is important to note that monosynaptic CCGs likely look quite different in different cortical regions, species, and in thalamocortical versus cortico-cortical connections. The study cited by the reviewer (Liew et al. 2021), for example, addressed monosynaptic connections between VPM and barrel cortex in rodents. The time lag and complexity of cortico-cortical interactions in macaque V1 likely differ substantially from that of thalamocortical connections in rodent barrel cortex (Stratford et al. 1996; Gil, Connors, Amitai, 1999; Gilman, Medalla, Luebke, 2017).

3. It is well established that fast-spiking neurons in the cortex receive stronger inputs from neurons in the local circuitry as compared to regular spiking neurons. Because fast-spiking and regular spiking neurons can be distinguished based on their spike waveform in extracellular recordings it could be interesting, but not required, to see whether this cell-type dependent connection strength is also evident in this dataset. This could add to the significance of the work and provide another angle to investigate circuit interactions. Again, this is only a suggestion.

We appreciate the reviewer’s suggestions. Indeed, we have previously classified the neurons described here into different categories based on their waveforms (Zhu et al., 2021, BioRxiv). We found that in addition to fast-spiking (FS) neurons, at least two-subclasses of regular-spiking (RS) neurons, namely regular-spiking medium (RSM) and regular-spiking long (RSL) neurons, exist within V1. This is generally consistent with evidence of the presence of more than two distinct classes of spike waveforms within mammalian cortex (e.g. Munoz, Tremblay and Rudy, 2014; Trainito, von Nicolai, Miller and Siegel, 2019) in addition to classes defined by burst propensity (Onorato et al., 2020). We found that the density of neurons from within these 3 putative waveform classes varied significantly across layers, and these laminar distributions were different among the 3 classes. Thus, although it would be interesting to examine whether fast-spiking and regular-spiking neurons exhibit differences in connection strength, we are cautious about overinterpreting the results from this analysis because (1) the 2 class FS vs RS categorization is probably an oversimplification of V1 population, and (2) the interactions between FS and RS neurons are likely confounded by their distinct laminar distributions and interactions with the local circuitry.

Moreover, although a number of previous studies have shown that fast-spiking (FS) neurons receive stronger and more dense input as compared to regular-spiking (RS) neurons, these studies were of thalamocortical connections and not cortico-cortical connections (Swadlow 2003; Callaway 2004; Cruikshank 2007; Bereshpolova, Hei, Alonso, and Swadlow 2020). As a result, in our data, the differences in the strength of input to RS and FS neurons may only be reflected indirectly in common-input pairs. It is upon future studies to more directly measure the strength of thalamic input to FS and RS neurons in monkey V1.

Nonetheless, as suggested by the reviewer, we examined differences in cell-type dependent connection strength between fast-spiking and regular-spiking neurons, separately for each cluster. We found that sharply synchronous pairs composed of two fast-spiking neurons (FS-FS) had significantly larger peak efficacies than sharply synchronous pairs composed of two regular-spiking neurons (RS-RS) (median peak efficacy: RS-RS = 0.030, FS-FS = 0.020; Wilcoxon rank sum test; p<105). In contrast, peak efficacy was not significantly different between FS-FS and RS-RS pairs in either the asynchronous or broadly synchronous classes (median peak efficacy: Async: RS-RS = 0.014, FS-FS = 0.016; Bsync: RS-RS = 0.022, FS-FS = 0.019; Wilcoxon rank sum test; Async, p=0.93, Bsync, p=0.01). Although the difference in peak efficacy for sharply synchronous RS-RS and FS-FS pairs observed here may reflect a difference in the strength of thalamic input to FS and RS neurons, we did not add these results to the manuscript due to the aforementioned concerns.

4. It would be useful to show the CCGs in Figure 1e in a higher temporal resolution around the peak such that the lag of the peak is visible.

We agree and have edited the new Figure 1e to show the CCGs in a higher temporal resolution around the peak.

5. It seems odd to treat the forward and reverse correlation functions as two distinct types since the labeling is arbitrary. That is, swapping the label for the reference neuron and its partner would flip the correlation function from one category to the other. Since the choice of the reference neuron is arbitrary, it is not clear why there are two categories. Please provide justification.

We agree with the reviewer that the forward and reverse distinction is arbitrary, and in the revised paper, we have modified and clarified our treatment of forward and reverse CCGs. As now described, at the CCG clustering stage, there is no arbitrary choice regarding which neuron is the reference. As we clarify (Section of “Classification of Cross-correlogram” under Methods), the forward and reverse clusters emerge as a result of including both directions of the CCG. In other words, either neuron was chosen as the reference neuron. Thus, at the clustering stage, the number of ‘pairs’ was twice the number of significant CCGs. At this stage, we merely treat the forward and reverse distinction as expected from the clustering.

By contrast, in subsequent analyses where we compare the proportions of the different classes across layers or functional cell types, the choice of the reference neuron is done either randomly (Figures. 5a-c, 6b-c, 7c) or deliberately (Figures. 5d, 6d, 7b, Figure 5—figure supplement 1), e.g., in comparing the proportion of forward and reverse CCGs among simple-complex cell pairs, simple cells were chosen as the reference. These analyses, along with figures 5-7 have also been revised for emphasis and clarity. Crucially, all of the results remain unchanged.

6. Separate from this issue, the labeling looks erroneous to the best of our knowledge. That is, 'forward' correlation functions have more mass at negative lags. This is opposite to the conventional definition. Mass at negative lags means the partner neuron tends to fire before the reference neuron, which is not a forward connection (reference driving partner) but the opposite. Please clarify.

We thank the reviewer for noticing the mislabeling. We have corrected this in the revised manuscript.

7. Please provide the species used.

We have added the species details (Macaca mulatta) to the section of “Experimental Model and Subject Details” under Methods in the revised manuscript.

8. Page 10: It was not obvious that the regression coefficients could be directly compared since the covariates (pair distance and r_ori) have very different magnitudes. Please clarify.

We thank the reviewer for noting this. Indeed, the covariates of pair distance and rori have very different magnitudes. As a result, we have standardized (z-scored) these two independent variables, pair distance and rori, prior to fitting the regression. Thus, the resulting regression coefficients can be interpreted as the change in the dependent variable associated with a one standard deviation (SD) change in the independent variable (see Kutner et al., Applied Linear Statistical Models, chpt. 7.5, pg. 272 for reference). Because we did not standardize the dependent variable, the units of the regression coefficients are the units of the dependent variable. For example, in the regression predicting peak lag, the coefficient of 0.94 for pair distance indicates that a one SD increase in the distance between two neurons is associated with a 0.94 ms increase in lag. We have amended the Results section “Dependence of synchrony and strength of functional interactions on tuning similarity” to clarify this.

9. Please provide an explanation of how the spike sorting was done. One concern about the sharp correlation peaks is that they are artifactual, resulting from issues in Kilosort (https://github.com/MouseLand/Kilosort/issues/29). In brief, Kilosort can match the same spike waveform to two different templates, one template capturing most of the waveform shape and the second capturing the residual. This issue is discussed on the Github Kilosort page. Are the authors sure that this does not contribute to/drive the sharp peak phenomenology? One way to provide reassurance would be to provide more information on the shape and width of the sharp peaks. If they are ~1 ms wide, it seems more likely to be artifactual than biological.

We thank the reviewer for pointing out the potential issue with spike sorting using Kilosort. Indeed, as the developer Marius Pachitariu acknowledged, this is an inherent problem using a template matching approach for spike sorting. Currently, there are no ideal solutions in Kilosort that fix this issue, so the users need to fix this themselves during the post-processing. We have added a section “Data acquisition and spike sorting” in the Methods to clarify this.

We adopted some common criteria in determining the occurrence of “double-counted spikes” (Siegle et al., 2021). For each pair of neurons with peak waveforms that are located within 50 μm (around 5 channels), we counted the number of spikes from one neuron that “overlap” with the other neuron, and then calculated the proportion of the “overlapping” spikes out of all the spikes in the given neuronal pair. Overlapping is defined as two spikes from the neuronal pair occurring within 5 samples (0.16 msec). We place no restriction on the amplitude relationship between the two neurons, although usually the double counted neuron from a residual exhibits much smaller amplitude.

Author response table 1
Percentage of overlapping spikes>5%>10%>20%
N_pair_overlapped26167
N_pair_overlapped
/N_pair_sharp synchronous
26/1408
=1.85%
16/1408
=1.14%
7/1408
=0.5%

As Joshua Siegle mentioned (https://github.com/MouseLand/Kilosort/issues/29), we should expect to see some perfectly overlapping spikes, just not very many. However, the exact boundary between true overlapping spikes and artificial overlapping spikes is not clear. In our case, we counted neuronal pairs within 50 μm that had at least 5%, 10%, and 20% overlapping spikes, and found that the number of pairs were 26, 16 and 7, respectively. Considering that we identified 1408 pairs (without considering directionality) that are sharply synchronous, those potentially artificial pairs contributed less than 2% of this class. Thus, although we didn’t remove these potentially contaminated pairs, we believe their contributions are negligible.

10. Please provide the proportion of simple and complex neurons in the data set.

The ratio of simple to complex neurons, respectively, was 1:2.4 in the dataset (236/802 neurons were simple, 566/802 neurons were complex). We have added this detail to the Results section (“Identifying functional interactions within single columns of visual cortex”) in the revised manuscript.

11. Discussion: The authors do a good job of discussing the limitations of correlation analysis for identifying synaptic connections. However, they refer to their R and F shapes (9% of cases) as potentially capturing such connections. Traditionally (e.g. Reid and Alonso), correlation function consistent with synaptic connections must pass a more stringent description than a simple asymmetry with respect to zero time lag. Namely, they need to have a sharp peak clearly offset from zero. This difference in stringency should be incorporated into the Discussion if the authors wish to propose that R and F types are related to connectivity.

We appreciate the reviewer’s suggestion and have revised the Discussion section to detail the variation in stringency of criteria used in past studies to identify monosynaptic connections (Section “Functional interactions and synaptic connectivity” under Discussions). The reviewer cites two such examples, namely the Liew et al. study (comment #2) and the Reid and Alonso studies (e.g., Reid and Alonso, 1995, 1999) (comment #11), which use different criteria. As we mention in the revised discussion, other, even more stringent, criteria have also been described, e.g., English et al., 2017; Senzai et al., 2019.

In the revised discussion, we also point out the obvious pitfalls of simply applying a more stringent criterion in order to identify monosynaptic connections. Although a key objective is to avoid false positives, the likelihood of false negatives is also a concern, particularly given that the variability in the parameters of such connections within a cortical circuit (e.g., synaptic latencies, synaptic strength) is not known. For example, consider the highly stringent criteria used by Buzsaki et al. (English et al., 2017; Senzai et al., 2019). Senzai et al., 2019 classified only ~0.5% of recorded pairs (525/102,332) as monosynaptic connections in mouse V1, which is considerably lower than the estimated rate of 5-12% “ground-truth” synaptic connections found using multi-patch slice recordings (e.g., Song et al., 2005; Cadwell et al., 2020). As the authors themselves explain in their discussion: “…detection of monosynaptic connections by the spike transmission probability method depends on the combination of the strength of the synapses and the number of spikes available in a given recording session, therefore it is biased toward neuronal pairs with stronger connections… This relationship may explain why only a small fraction of all possible E-E connections was detected” (Senzai., et al., 2019).

In addition, we note that the criteria used in the other example study cited by the reviewer (Liew et al., comment #2), is actually less stringent than ours in terms of the SD threshold. Using their criteria, which only requires the peak to be >3.5 SD above the shuffled background, our data yield 22,103 significant asynchronous connections. This is in stark contrast to the 6,124 significant asynchronous connections (asynchronous types) obtained with our 7 SD threshold, which is based on the threshold used by Siegle et al., 2021.

12. Line 635: The definition of the correlation function seems different from previous work. The correlation function appears normalized by a quantity that depends on the time lag of one neuron and not the other (the denominator). Please clarify the notation.

The difference is primarily notational, as we clarify here. The cross-correlation function used in previous work (e.g., Smith and Kohn 2005, 2008; Siegel et al. 2021) is,

CCG(τ) = 1Mi=1Mt=1Tτxi1(t) xi2(t+τ)(T|τ|)λ1λ2 (Equation 1)

where M is the number of trials, T is the trial duration in milliseconds, xi1 is the spike train of neuron 1 in trial i, and λ1 is the firing rate of neuron 1. In this formulation, the cross-correlation is normalized by the product of the geometric mean of the firing rates of two neurons (λ1λ2) and a triangle function of lag (T|τ|). This triangle function is a quantity that depends on the time lag of one neuron. Our definition of the normalization term is equivalent to that of previous work if the firing rates λ1 and λ2 are computed using only the bins used to compute the CCG at each time lag, i.e.,

λ1(τ)=i=1Mt=1Tτxi1(t)(Tτ)×M and λ2(τ)=i=1Mt=1Tτxi2(t+τ)(Tτ)×M.

This can be verified via substitution of λ1(τ) and λ2(τ) into Equation 1.

The difference between computing the firing rate normalization across all time bins or only time bins used to compute the CCG is typically negligible. However, the latter is more desirable when computing the CCG between neurons with non-stationary PSTHs that fluctuate at the beginning or end of a trial. To clarify the connection between our definition and that of previous work, we have modified the notation of the cross-correlation function in the revised manuscript (Section “Cross-correlograms (CCGs)” under Methods).

13. Figure 2B: There is an oscillation in the marginal histogram of peak delays. This is also evident in Figure 3C. Please clarify.

The peak delay is an integer between 0 and 10 ms for significant CCGs because CCGs were computed using spikes binned with 1 ms resolution. The oscillations in the marginal distribution were an artifact of applying a KDE function to this discrete distribution.

In the revised manuscript, we have fixed this in Figure. 2b, 3c, and 3e by illustrating raw histograms in the marginals instead of KDE functions of the histograms.

https://doi.org/10.7554/eLife.79322.sa2

Article and author information

Author details

  1. Ethan B Trepka

    1. Department of Neurobiology, Howard Hughes Medical Institute, Stanford University, Stanford, United States
    2. Neurosciences Program, Stanford University, Stanford, United States
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Supervision, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    Contributed equally with
    Shude Zhu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3002-800X
  2. Shude Zhu

    Department of Neurobiology, Howard Hughes Medical Institute, Stanford University, Stanford, United States
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    Contributed equally with
    Ethan B Trepka
    For correspondence
    shude@stanford.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8674-9607
  3. Ruobing Xia

    Department of Neurobiology, Howard Hughes Medical Institute, Stanford University, Stanford, United States
    Contribution
    Conceptualization, Software, Formal analysis, Investigation, Visualization, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
  4. Xiaomo Chen

    1. Department of Neurobiology, Howard Hughes Medical Institute, Stanford University, Stanford, United States
    2. Center for Neuroscience, Department of Neurobiology, Physiology, and Behavior, University of California, Davis, Davis, United States
    Contribution
    Resources, Data curation, Funding acquisition, Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
  5. Tirin Moore

    Department of Neurobiology, Howard Hughes Medical Institute, Stanford University, Stanford, United States
    Contribution
    Conceptualization, Supervision, Funding acquisition, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    Competing interests
    Senior editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3345-2930

Funding

National Institute of Neurological Disorders and Stroke (NS116623)

  • Tirin Moore

National Eye Institute (EY014924)

  • Tirin Moore

National Eye Institute (EY029759)

  • Xiaomo Chen

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Jonathan C Horton for extensive help with the recordings and histology. We thank Tim Harris and Karel Svoboda for providing the Neuropixels probes, Shellie Hyde and Sam Baker for technical assistance.

Ethics

All experimental procedures were in accordance with National Institutes of Health Guide for the Care and Use of Laboratory Animals, the Society for Neuroscience Guidelines and Policies, and with approved Institutional Animal Care and Use Committee (IACUC) protocol (#APLAC-9900) of Stanford University.

Senior Editor

  1. Timothy E Behrens, University of Oxford, United Kingdom

Reviewing Editor

  1. Martin Vinck, Ernst Strüngmann Institute (ESI) for Neuroscience in Cooperation with Max Planck Society, Germany

Reviewer

  1. Jens Kremkow, Charite-Universitätzmedizin Berlin, Germany

Version history

  1. Preprint posted: February 21, 2022 (view preprint)
  2. Received: April 7, 2022
  3. Accepted: October 30, 2022
  4. Accepted Manuscript published: November 2, 2022 (version 1)
  5. Version of Record published: November 14, 2022 (version 2)

Copyright

© 2022, Trepka, Zhu et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,290
    Page views
  • 211
    Downloads
  • 3
    Citations

Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Ethan B Trepka
  2. Shude Zhu
  3. Ruobing Xia
  4. Xiaomo Chen
  5. Tirin Moore
(2022)
Functional interactions among neurons within single columns of macaque V1
eLife 11:e79322.
https://doi.org/10.7554/eLife.79322

Further reading

    1. Computational and Systems Biology
    2. Neuroscience
    Huu Hoang, Shinichiro Tsutsumi ... Keisuke Toyama
    Research Article Updated

    Cerebellar climbing fibers convey diverse signals, but how they are organized in the compartmental structure of the cerebellar cortex during learning remains largely unclear. We analyzed a large amount of coordinate-localized two-photon imaging data from cerebellar Crus II in mice undergoing ‘Go/No-go’ reinforcement learning. Tensor component analysis revealed that a majority of climbing fiber inputs to Purkinje cells were reduced to only four functional components, corresponding to accurate timing control of motor initiation related to a Go cue, cognitive error-based learning, reward processing, and inhibition of erroneous behaviors after a No-go cue. Changes in neural activities during learning of the first two components were correlated with corresponding changes in timing control and error learning across animals, indirectly suggesting causal relationships. Spatial distribution of these components coincided well with boundaries of Aldolase-C/zebrin II expression in Purkinje cells, whereas several components are mixed in single neurons. Synchronization within individual components was bidirectionally regulated according to specific task contexts and learning stages. These findings suggest that, in close collaborations with other brain regions including the inferior olive nucleus, the cerebellum, based on anatomical compartments, reduces dimensions of the learning space by dynamically organizing multiple functional components, a feature that may inspire new-generation AI designs.

    1. Neuroscience
    Max L Sterling, Ruben Teunisse, Bernhard Englitz
    Tools and Resources Updated

    Ultrasonic vocalizations (USVs) fulfill an important role in communication and navigation in many species. Because of their social and affective significance, rodent USVs are increasingly used as a behavioral measure in neurodevelopmental and neurolinguistic research. Reliably attributing USVs to their emitter during close interactions has emerged as a difficult, key challenge. If addressed, all subsequent analyses gain substantial confidence. We present a hybrid ultrasonic tracking system, Hybrid Vocalization Localizer (HyVL), that synergistically integrates a high-resolution acoustic camera with high-quality ultrasonic microphones. HyVL is the first to achieve millimeter precision (~3.4–4.8 mm, 91% assigned) in localizing USVs, ~3× better than other systems, approaching the physical limits (mouse snout ~10 mm). We analyze mouse courtship interactions and demonstrate that males and females vocalize in starkly different relative spatial positions, and that the fraction of female vocalizations has likely been overestimated previously due to imprecise localization. Further, we find that when two male mice interact with one female, one of the males takes a dominant role in the interaction both in terms of the vocalization rate and the location relative to the female. HyVL substantially improves the precision with which social communication between rodents can be studied. It is also affordable, open-source, easy to set up, can be integrated with existing setups, and reduces the required number of experiments and animals.