Layer-Specific Physiological Features and Interlaminar Interactions in the Primary Visual Cortex of the Mouse

1A characteristic feature of the neocortex is its laminar organization. The cortical columnar microcircuitry is viewed as a stack of interconnected yet distinct neuronal networks in which each lamina possesses somewhat unique patterns with different specific inputs, projection targets, and feedback connections. How the laminar structure relates to mesoscopic physiological patterns, such as local field potential (LFP) oscillations and physiological interactions of single neurons across layers, is not well understood. Numerous oscillatory and transient LFP patterns of functional relevance have been described in various neocortical regions. However, their relationship to afferents, intracortical connectivity, and the firing patterns of individual neurons has been largely unexplored.

2Using multisite recording silicon probes that span all cortical layers, we sought to characterize layer-specific physiological patterns and neuronal cross-talk between layers in the primary visual cortex (V1) of freely behaving mice. The well-characterized anatomical connectivity and the diverse neuronal types of V1 make this cortical area ideal to investigate how the different cell types in different layers interact during physiological operations, such as sensory processing and offline states, such as sleep. Previous recordings in the hippocampus, using a similar approach, have led to a solid understanding of the relationship between extracellular signals and anatomical connectivity. LFP patterns can identify layers in vivo with <25–50 μm precision and have been used to identify the unique relationship between upstream activity levels and hippocampal spike outputs and relate such input-output transformation to behavior. Similar strategies have been followed in the V1 of waking monkeys and other cortical areas in rodents, mostly under anesthesia. In these previous studies, layer boundaries were estimated mainly by depth criteria, and the relationship among LFP depth profiles, neuronal activity, and interlayer interactions in different behavioral states was not addressed quantitatively.

3Our goal was to define layers functionally in V1 of the mouse and relate them to classical anatomical layers. Because spontaneous mesoscopic patterns in the primate V1 are characterized by gamma, alpha, and slow oscillations, we searched for their corresponding patterns in the mouse. Finally, we examined the brain state-dependent spike transmission probabilities among putative principal cells and inhibitory interneurons across cortical layers in waking and sleep to quantify their brain state dependence.

4Mice (n = 19) were implanted with a single-shank, 64-site linear silicon probe and recorded during free behavior in their home cage. The probe was placed parallel with the orientation of the apical dendrites of pyramidal neurons in the V1 (Figure 1). In our initial experiments, small electrolytic lesions were made to calibrate the positions of the recording sites with histological verification (n = 4; Figures 1A and S1). Subsequently, we used the characteristic depth distribution of LFPs and unit firing to identify five physiological landmarks that reliably anchored the recording sites in V1 across animals and recordings. The most prominent landmark was a large-amplitude peak of the depth profile of power between 500 Hz and 5 kHz corresponding to mid-layer 5 (Figure 1C, depth c). This peak likely reflected the aggregated power from the high firing rates of the large layer 5 pyramidal. The remaining four landmarks (Figures 1D and 1I, depths a, b, d, and e) were obtained from the depth distributions of the most prominent current source density (CSD) sinks and sources of spontaneous slow oscillations of non-rapid eye movement (REM) sleep.

5To investigate the origin of landmarks a–e, we compared the slow oscillation sinks and sources with those of light-evoked responses in the waking mouse (Figures 1E and 1J; 100 ms light pulses). Visual stimulation produced an early sink (b′) at the border between layers 3 and 4, which likely reflected afferent activation from the lateral geniculate body. The early sink was accompanied by a source in superficial layers (a′) and a sink-source pair corresponding to the border between layers 5 and 6 (sink, d′) and lower layer 6 (source, e′). Sink b, associated with the “UP” state of slow oscillation during non-REM sleep, was slightly but significantly deeper than the light-evoked sink b′ in the layer 3/4 border (39 ± 49 μm; p < 0.05, signed rank test; n = 13 mice; Figure 1F), whereas the other sinks and sources occurred at the same depth as those seen with visual stimulation (p > 0.05, signed rank test). Importantly, the distances between the five landmarks were consistent across animals and varied <100 μm from mouse to mouse (Figure 1G; n = 19). The consistent depth localization of spike power and CSD sinks and sources during non-REM sleep across animals allowed us to project all subsequent physiological data onto normalized depth coordinates and construct average maps across all animals (Figures 1I and 1J; STAR Methods).

6To test whether layers could be further disambiguated by the depth localization of LFPs, we examined coherence in the gamma band (30–100 Hz) across all recording sites (Figures 2A and S2). Each site served as a reference, and pairwise coherence with each of the remaining 63 sites was determined iteratively. Using a gradient-descent algorithm (; STAR Methods), the recording sites were clustered on the basis of their resulting cross-coherence matrix, resulting in six clusters whose five boundaries were consistent across animals when projected into normalized depth coordinates (Figure 2A; p > 0.05, signed rank test for cluster boundaries across animals). The depth distribution of the six gamma coherence clusters corresponded roughly to the six cortical layers in each brain state (Figures 2B and S2; p > 0.05, ANOVA for comparison across different states), suggesting that gamma oscillations in V1 have a highly specific laminar structure. To isolate layer-specific gamma oscillations by another method, we also used independent component analysis (ICA) decomposition of the band-passed (30–100 Hz) LFPs (; Figure S2). ICA yielded six main independent components (ICs) that were highly consistent across animals (p < 0.001, Spearman’s correlation for ICs spatial loadings), and their depth distribution was related to the anatomically defined layers (Figure 2C). This method has the additional advantage of removing volume-conducted currents that accounted for 0.34 ± 0.11 of V1 gamma LFP variance. The superficial gamma ICs had a significantly larger power compared with the deeper ones (top three versus bottom three; Figures 2D and S4; 0.14 ± 0.04 versus 0.08 ± 0.03 relative power for ICs 1–3 versus 4–6; p = 1.1 × 10−7, rank sum test).

7Spectral analysis of the LFP traces revealed depth- and brain state-dependent variation of characteristic frequencies (Figures 3 and S3). Slow oscillation power (0.5–3 Hz) was dominant in non-REM sleep (see Supplementary Information; brain state comparison: p = 4.0 × 10−4, Kruskal-Wallis test; n = 19 mice), especially in deeper layers (layer comparison: p = 4.1 × 10−5, Kruskal-Wallis test; Figures 3C and S3). Slow oscillation was phase-coupled to gamma frequency power in superficial layers (69 ± 12 Hz; modulation index 2.7 ± 0.5 × 10−3) and to the high-gamma (>100 Hz) band in deep layers (Figure 3D; modulation index 2.3 ± 0.3 × 10−3). Slow gamma oscillations (30–60 Hz) dominated in superficial layers, especially in the waking state followed by REM and non-REM sleep (Figures 3 and S3; p = 0.028).

8Although LFP phase coherence rapidly decreased across layers, especially at higher frequencies (Figure S2A), phase-amplitude cross-frequency coupling revealed reliable interlaminar interactions. In the waking animal, the dominant 3–6 Hz band was phase-coupled to a 15–25 Hz band signal in multiple layers (22 ± 5 Hz; modulation index 4.8 ± 1.2 × 10−4), possibly corresponding to the spiky component of the pattern, to high gamma frequency (60–100 Hz) in the superficial layers (90 ± 13 Hz, modulation index 2.7 ± 0.7 × 10−4), and to the epsilon band (>100 Hz) in layer 5 (modulation index 6.4 ± 1.8 × 10−4), likely reflecting unit firing (Figure 3D).

9The most prominent LFP pattern in the waking mouse was a 3–6 Hz oscillation with largest power in layer 4 (Figure 3C, level b; p = 2.0 × 10−6, Kruskal-Wallis test; n = 19 mice). This strong power increase was due mainly to frequent epochs of a 3–6 Hz rhythm during waking immobility, in which the waveforms occasionally split into short “spike” and longer “wave” components with a depth profile distribution similar to the depth distribution of alpha waves (8–12 Hz) in monkey V1 (; Figures 3A and 3B). The 3–6 Hz band was cross-frequency phase-coupled to a 15–25 Hz band signal in multiple layers (22 ± 5 Hz, modulation index 4.8 ± 1.2 × 10−4), possibly corresponding to the spiky component of the pattern, to high gamma frequency (60–100 Hz) in the superficial layers (90 ± 13 Hz, modulation index 2.7 ± 0.7 × 10−4), and to the epsilon band (>100 Hz) in layer 5 (modulation index 6.4 ± 1.8 × 10−4), likely reflecting unit firing (Figure 3D).

10CSD analysis of the 3–6 Hz rhythm showed a depth distribution of sinks and sources similar to those of the slow oscillation (compare Figure 3E with Figure 1I). The duration of the non-spiking periods of this waking rhythm (50–150 ms) was also similar to the DOWN state of slow oscillation. However, whereas the spiking-associated phase of the waking 3–6 Hz rhythm was short, the duration of the UP state of slow oscillation varied extensively and showed a lognormal distribution. Both principal cells and interneurons were strongly phase-locked to the trough of the layer 4-recorded 3–6 Hz rhythm (Figure 3F).

11Volume-conducted theta power (6–9 Hz) from the hippocampus allowed us to distinguish between waking immobility and ambulation (STAR Methods). As soon as the animal started moving, the 3–6 Hz rhythm vanished, accompanied by increased slow gamma power. The antagonism between movement and the 3–6 Hz oscillation was quantified by the negative correlation between its power and the power of theta oscillation (p < 0.01, signed rank test; Figure S2C). Visual stimulation also blocked an ongoing 3–6 Hz rhythm. However, the offset of visual stimulus often induced a short rebound train 3–6 Hz oscillation (; Figure S2I).

12Because “layer” designation is traditionally based on the depth distribution of cell bodies, whereas our depth classification of the LFP patterns (Figure 2) reflects largely a combination of afferent-dendritic excitation and inhibition, we sought to establish a disciplined way to relate LFP-, CSD-, and ICA-based depth estimation to the recorded neurons in different anatomically defined layers.

13Unit clustering was performed as described previously using one session from each of the 19 mice, yielding a total of 1,472 units. The physiological identity of clustered units was determined using a multi-step approach (Figure 4). First, the unfiltered waveform was quantified by the trough-to-peak latency of the extracellular spike. This initial classification yielded a bimodal distribution of units with narrow-waveform spikes (putative interneurons, n = 251 I cells) and wide-waveform spikes (putative principal cells, n = 1,075 E cells; Figure 4A; for unclassified units, see Figure S4). We next took advantage of the simultaneously recorded units to physiologically identify them as E cells or I cells by their short-latency temporal interactions with other neurons. Putative monosynaptic connections are associated with precisely timed spike transmission at short-latency (<4 ms) offsets between two recorded neurons, as detected by narrow significant peaks (excitatory) or troughs (inhibitory) in the cross-correlogram (CCG; Figures 4B and 4C;, ). Using these criteria, 347 and 100 units were classified as CCG-based E and I cells, respectively (Figure 4A). These CCG-based classification results overlapped well with the waveform-based classification results (97.4% of CCG-based E cells were classified as waveform-based E cells, and 98.0% of CCG-based I cells were classified as waveform-based I cells), supporting the adequate separation on the basis of trough-to-peak latency feature.

14We also validated the physiological identification of inhibitory neurons using an optogenetic strategy. A subset of neurons in PV-Cre::Ai32 (n = 5) and VGAT-Cre::Ai32 (n = 9) mice were identified optogenetically as inhibitory neurons (Figures 4D and 4E; blue and green diamonds, respectively, in Figure 4A). Fourteen of 15 optogenetically identified PV cells (n = 11) and 40 of 52 VGAT cells were part of the inhibitory group defined by the spike waveform criterion. These results further support the adequate separation of I cells on the basis of waveform only and also show that a minor fraction of inhibitory neurons (n = 12 of 52) had wider waveforms and were initially falsely combined with the E cell group.

15To link neurons to physiological layers, we examined the relationship between the six gamma ICs and spiking of all neurons. Similar to the LFP-only analysis, spike-IC phase coupling (STAR Methods) also showed clear depth stratification, indicating that units were preferentially modulated by local gamma oscillations (Figures 5A and S5). Around and above depth level b (i.e., layer 3/4 border; first three components), unit-IC coupling of both putative E and I neurons had a peak in the slow gamma range. In deep layers, unit-IC coupling shifted toward the fast gamma band. This fast gamma (60–100 Hz) phase preference also segregated E and I neurons along the cortical depth (Figure S5). We therefore used the spike-gamma phase coupling profile of all units against the six ICs to perform K-means clustering (STAR Methods). This procedure yielded eight clusters (Figure 5B). By projecting back the K-means-clustered units to cortical depth, we found that the clusters largely corresponded to the anatomically defined cell body layers (Figure 5C). The top two clusters may correspond to neurons in layers 2 and 3 (blue and cyan), and the third cluster surrounded landmark b (green) and corresponded to the layer 3/4 border sink (Figure 1I). Three clusters surrounded landmark c (yellow, red, and pink), the presumed center of layer 5, with the lowest of the three clusters surrounding landmark d, the deep sink, likely corresponding to layer 5/6. The bottom two clusters (black and gray) were designated as layer 6 neurons. The remaining unit analyses were performed on these clustered groups of neurons.

16The various physiological features of the recorded principal neurons as a function of depth are summarized in Figure 6 (see Figure S6 for interneurons). Long-term mean firing rates of V1 principal cells showed a strongly skewed distribution in each layer. The fastest firing group corresponded to those in layers 5A and B (Figure 6, red and magenta groups), while the slowest ones were the superficial layers (Figure 6A; p = 3.7 × 10−83, Kruskal-Wallis test). Deep-layer neurons were significantly more active during waking compared with non-REM sleep (Figure 6B; see Figure S6 for other state comparisons).

17A characteristic physiological feature of cortical pyramidal neurons is their bursting property. Bursting was determined by calculating the number of spikes in the 3–10 ms bins of the spike autocorrelograms divided by the number of spikes in the 200–300 ms bin. In general, superficial layer neurons exhibited more bursting behavior, followed by layer 6 and layer 4 neurons, while layer 5 neurons had the lowest burst indices (Figure 6C; p = 5.2 × 10−37, Kruskal-Wallis test). Another measure of spike dynamic is the relationship of the neuron’s initial phasic response versus sustained firing to a given input. We calculated this transient/steady firing rate index by dividing the peak firing rates at 0–200 ms by the mean firing rates at 100–200 ms after the UP state onset of slow oscillation (Figures 6D and S2C). By this measure, deep layer 5 neurons showed the strongest transient responses (layer difference, p = 5.6 × 10−10). During the DOWN-UP state transition, layer 5B neurons fired first, followed by superficial and layer 6 neurons (,; Figure 6E; p = 1.3 × 10−38). The distribution of spike firing latency during the 3–6 Hz rhythm of waking displayed characteristic laminar differences (Figure 6F; p = 2.6 × 10−10), similar to spike latencies during slow oscillations. The mean vector length of phase locking to the 3–6 Hz rhythm was strongest in the superficial layers and weakest in layer 6 (Figure 6G; p = 7.6 × 10−50). Spike-gamma (30–100 Hz) coupling was stronger in the waking animal compared with non-REM in all layers but layer 5A (Figure 6H).

18The strongest spike-LFP coupling was observed in the 3–6 Hz band (Figures 3F and S7), followed by the slow gamma (30–60 Hz) band, each with characteristic cortical depth and brain state-dependent profiles (Figure S7).

19Of the 102,332 potential neuron pairs in all mice and sessions, we identified putative excitatory (E, n = 525) and inhibitory (I, n = 323) monosynaptic connections by the spike transmission probability method (Figures 7 and S7). This is a strongly under-sampled estimate, because the spike transmission probability-based method requires simultaneous recordings of spikes of both pre- and postsynaptic neurons, and large numbers of spikes are needed to detect weak synapses. Although many identified excitatory connections were detected within the same layer (layers in these calculations refer to IC-based clusters) for both E-I (Figures 7A and 7B; n = 173 of 385 connections) and E-E (n = 19 of 90) pairs as well as inhibitory I-E pairs (Figure S8; n = 98 of 323 connections), there were also many cross-layer excitatory and inhibitory pairs (I-I pairs were too few to draw meaningful conclusions). A considerable fraction of layer 2/3 principal cells had monosynaptic connections with layer 4 and 5 interneurons (54.7%) and principal cells (22.7%) (Figures 7E and 7G). Of the 172 putative excitatory connections of L2/3 pyramidal cells, 10 (5.8%) and 39 (22.7%) targeted L2/3 and L4/5/6 excitatory neurons, respectively, in addition to contacting 29 (16.9%) and 94 (54.7%) L2/3 and L4/5 interneurons, respectively. Most individual presynaptic principal neurons had only a single partner target in the recorded population but a small fraction targeted multiple partners (Figures 7B and 7C). A typical interneuron typically had one presynaptic partner, but a small fraction many, up to 16 (Figure 7C). Almost half of the E-I connections (168 of 392 [42.9%]) received reciprocal I-E connections when all layers were considered. Reciprocal E-I-E connections were almost twice more frequent in layer 5 (60%–70%) than in other layers (Figure S7C).

20To quantify the brain state-dependent changes in the spike transmission probabilities, we calculated the peak-to-baseline ratio of CCGs between E-I cells (Figures 7D and 7E) and between E-E cells (Figures 7F and 7G) as well as between I-E pairs (Figure S8) in different brain states. We observed that spike transmission strength between layer 2/3 E cells and deep layer I and E cells became weaker during non-REM sleep compared with waking (L2/3 E cells to L5 I cells, p = 0.0012; L2/3 E cells to L4 I cells, p = 0.0010; and L2/3 E cells to L5 E cells, p = 0.015; signed rank test; Figures 7D–7F, 7H, and 7I). In contrast, putative E-E spike transmission in the deep layers (layer 5/6 to layer 5/6 principal cells) significantly increased during non-REM sleep compared with waking (p = 0.0061, signed rank test; Figures 7G and 7I). These results indicate that communication from superficial to deep layers increases during waking, while excitatory interactions become stronger within deep layers during non-REM sleep. These findings were also reflected by unit-LFP coupling when all eight unit clusters were compared with the six gamma ICs (Figure 7I). When all I-E pairs were considered, spike suppression probability during waking was significantly stronger compared with non-REM (I-E dip wake, 24.1%; non-REM sleep, 22.8% relative to the CCG baseline; p = 0.006, signed rank test; Figure S8F).

21We found five neurons in four mice that were specifically and dominantly active during the DOWN state of non-REM slow oscillation (Figures 8A and 8B). All DOWN state-active neurons were found in deep layers 6. Four neurons had bursty autocorrelograms, while the fifth one was regular spiking. These five neurons were not only active during the DOWN state but mainly silent during the UP state (35.2 ± 23.1 Hz during DOWN, 4.3 ± 3.1 Hz during UP; p < 0.016, rank sum test). During light flash presentation, their spiking activity was further reduced at the time when the remaining neurons robustly responded (Figure 8C). During optogenetic activation of PV neurons or VGAT neurons, the DOWN state-active neurons were also suppressed, indicating that at least some classes of interneurons innervate them. However, ∼30 ms after the optogenetic stimulation, they rebounded to high level of activity, while the remaining population of principal neurons was still suppressed (Figure 8D). Cross-correlation between the DOWN state-active cells and other principal cells showed an inverse correlation not only during non-REM but also in the waking animal (Figures 8E–8H). In summary, the rare DOWN state-active neurons showed a robust inverse correlation with the activity of all other neurons during both sleep and waking.

22Depth profiles of unit power and sink-source distributions of slow oscillations of non-REM sleep provided consistent physiological landmarks in V1 across animals. We report the following findings: (1) Coherence and ICA of gamma oscillations (30–100 Hz) and spike-gamma LFP coupling identified six physiological layers and distinguished further sublayers. (2) Firing rates, burstiness, and other physiological features of neurons displayed distinct layer and brain state dependence. (3) Monosynaptic connections, quantified by a spike transmission probability method, revealed highly structured distributions within and across layers. Connection strengths were skewed, with a minority of highly connected hubs. Spike transmission between E-E pairs and E-I pairs from layer 2/3 to layer 5 was stronger during waking compared with non-REM sleep but stronger among deep-layer excitatory neuron pairs during non-REM sleep. (4) The most prominent LFP pattern during waking was a 3–6 Hz rhythm with characteristic phase preferences of spikes across layers. (5) Spiking of a small subset of neurons in deep layers was anticorrelated with all other neurons, and these neurons were most active in the DOWN state of slow oscillations. Our findings link mesoscopic LFPs and single-neuron interactions with multilayered anatomical organization in V1.

23The term “layer” has been used differently in different structures. For example, in the hippocampus, apical and basal dendritic layers are distinguished from somatic layers of the same neurons. In contrast, “layer” in the neocortex traditionally refers to histologically distinct somatic layers. However, a given somatic layer (e.g., layer 2) is also the apical dendritic layer of other neurons (e.g., layer 3 and layer 5 neurons). Afferents from different upstream regions can preferentially target neurons of a given layer or multiple layers. For example, although the highest density of thalamocortical afferents is present in layer 4, collaterals of these axons also innervate both superficial and deep layers, contacting several types of neurons. Nevertheless, afferents from different sources typically segregate on different segments of dendrites. Similarly, inhibitory interneurons of the same class converge on the same neuronal domains, whereas different types converge on distinct somadendritic domains of their targets. In addition, biophysical experiments demonstrate that different segments of pyramidal neurons are electrically isolated, and each functional class of inputs is initially processed in relative independence of the other. Our experiments identified these physiologically distinct layers.

24Previous work in the hippocampus has shown that coherence in the gamma frequency range, measured across electrodes in the same layer innervated by the same afferents, is high over long lateral distances, whereas coherence across electrodes placed just 100 μm apart but in different layers is low. Similar segregation by coherence has been described between superficial and deep layers in the V1 of the monkey. We hypothesize that segregation of excitatory and inhibitory inputs on the orderly arranged principal cell populations may be responsible for the observed layer-specific extracellular gamma currents.

25Both the gradient-descent algorithm performed on gamma oscillation and ICA identified six strata, which can be regarded as functionally distinct layers. Spike-gamma LFP phase-coupling, in turn, allowed us to relate these physiological layers precisely to the depth distribution of neuronal somata. This additional step effectively separated layer 6 neurons into deep and superficial groups and divided layer 5 neurons into three subgroups. By depth criteria, two of these three groups may be classified as layer 5A (cortico-cortical with thin apical dendrites) and 5B (cortical-subcortical with thick apical dendrites) neurons, respectively. Our third, most superficial group may represent a transitional form between layer 4 and layer 5 neurons. Our clustering method also identified two groups in layers 2/3, possibly corresponding to layer 2 and layer 3 neurons. Our physiology-based classification of principal neurons will require confirmation by future optogenetic experiments using available genetic markers of layer and sublayer-specific pyramidal neurons. Overall, our findings demonstrate that physiological properties of neurons, especially their relationship to gamma LFP, can be exploited to relate them to classical anatomy-based layer segregation. In turn, the spiking activity of the classified groups can be examined for their contribution to brain state-dependent collective network patterns.

26The most prominent LFP pattern in the waking V1 was a 3–6 Hz oscillation. We hypothesize that this rhythm is an evolutionary precursor of the primate alpha activity in the visual cortex. First, the 3–6 Hz oscillation is the dominant rhythm in the waking mouse as is the case for the 8–12 Hz alpha oscillation in the primate. Its power is highest in layer 4, and its current sink-source depth distribution is similar to the alpha waves in the monkey. Another similarity to alpha oscillations is the strong phase-power and power-power coupling between 3–6 Hz and gamma oscillations. Second, the behavioral correlates of the 3–6 Hz oscillations are similar to alpha. Both are most prominent during quiet waking, occur in several seconds-long bouts, and are reduced or eliminated by movement and arousal and reduced in non-REM sleep (,; present observations). Yet their probability can be increased by the offset of visual stimulation in the mouse (; present findings) and rat and by milk reward in the cat. When visual cues are presented during the 3–6 Hz oscillations, the responses of V1 pyramidal neurons are reduced, similar to attenuated visual evoked responses in humans. Intracellular recording in layer 2/3 pyramidal cells, parvalbumin and somatostatin-expressing neurons in V1 of mice revealed large amplitude, rhythmic hyperpolarizations at 3–6 Hz after the offset of visual stimulation, and spiking between the hyperpolarization, as observed for interneurons, including parvalbumin-expressing cells in our study. Similar rhythmic hyperpolarizations were observed in layer 2/3 pyramidal neurons at 2 Hz under anesthesia and blocked by visual stimulation. Third, we observed the strongest spike-LFP coupling in all cortical layers in the 3–6 Hz band. The earliest and strongest firing was present in layers 4 and 5B, coinciding with the LFP sink in layer 4 followed by neurons in the superficial layers, similar to neuronal patterns during alpha oscillations in different cortical layers in the monkey.

27We found similar sink-source distributions during the 3–6 Hz rhythms and slow oscillations of non-REM sleep. It is generally believed that the neuronal substrate of several slow-frequency oscillations in the neocortex is the thalamocortical resonant network and that the distinct patterns, such as slow oscillations, alpha rhythm, mu rhythm, and sleep spindles, and their timing circuits are set largely by the combination of subcortical neuromodulatory inputs affecting both thalamic and neocortical targets. The increased tone of cholinergic and monoaminergic afferents during movement, arousal, and visual stimulation are known to suppress thalamocortical oscillations, including alpha oscillations in humans, 6–9 Hz thalamocortical oscillations in rats, and the 3–6 Hz rhythm in mice. In sum, although the frequency of the primate alpha is more than twice higher than the 3–6 Hz rhythm in the mouse, many of their physiological features are remarkably similar. In support of the evolutionary conservation hypothesis, the frequency of the V1 rhythm increases with brain size from mouse (4 Hz; ) to cat (6–12 Hz;, ) to human (8–12 Hz; ).

28Both firing rates and bursting properties of principal cells varied across different layers. In agreement with previous observations in different cortical areas and species, layer 2/3 pyramidal neurons had the lowest median firing rates, followed by layer 4 and layer 6 neurons, whereas the highest firing rate cells resided in layer 5. The lower firing rates of layer 2/3 excitatory neurons may, at least in part, result from their resting membrane potentials being ∼10 mV hyperpolarized relative to superficial layer pyramidal neurons, according to in vitro measurements. Layer 2/3 pyramidal neurons may therefore require substantially more excitatory synaptic input to drive them to action potential threshold compared with L5 pyramidal neurons. Deep-layer neurons increased their excitability during waking, as quantified by their firing rates, whereas no such state dependence was present in superficial neurons.

29Layer 5 neurons are often segregated into intrinsically bursting and regular firing neurons on the basis of in response to intracellular current injection. In our experiments, we found that layer 5 neurons, in general, had the lowest propensity to induce spike clusters at short interspike intervals. On the other hand, layer 5 neurons had the highest index of transient/steady ratio during the UP state of slow oscillation, reminiscent of the neuron’s response to sustained intracellular depolarization.

30We identified monosynaptic connections between pairs of neurons both within and across layers using short-term spike transmission probability. The magnitude of both spike transmission probability and single neuron-induced excitatory postsynaptic potentials (EPSPs) (in E-I and E-E connections) and IPSPs (I-E connections) showed a lognormal distribution, implying that a very small number of strong connections are mixed with a large number of weak connections. Detection of monosynaptic connections by the spike transmission probability method depends on the combination of the strength of the synapse and the number of spikes available in a given recording session, therefore it is biased toward neuronal pairs with stronger connections, such as E-I connections and I-E pairs. This relationship may explain why only a small fraction of all possible E-E connections was detected and why more E-I than E-E pairs were identified in our study. Yet our under-sampled population can still be regarded as representative for the relative distribution of the connectivity features.

31Surprisingly, most E-I and I-E connections were across layers, including deep-to-superficial and superficial-to-deep links. Layer 2/3 neurons innervated >3 times as many layer 4/5 principal cells and interneurons as their peers in the home layer. Connection strength tended to cluster around a few hub neurons. We also found that reciprocal connections between E and I neurons varied across layers and were almost twice as frequent in layer 5 compared with other layers. These findings indicate that cortical connectivity is highly structured, and interlaminar interactions are fundamental to multi-layer cortical operations. The use of higher density, multishank probes in future experiments will undoubtedly increase the yield of synaptically connected neurons and provide more reliable quantification than we found in the present experiments. This expected progress is facilitated by the fact that the number of neuron pairs increases as the square of the number of units recorded.

32Spike transmission probability between superficial and deep layers was stronger in the waking animal for both E-I and E-E pairs, compared with non-REM sleep. This finding, combined with the suggestion that thalamocortical synapses are not affected by brain state, supports the hypothesis of top-down control of visual processing. In contrast, spike transmission strength increased between E-E pairs and I-E pairs during non-REM sleep, preserving the E-I balance even during the UP state of slow oscillation, associated with strong recurrent excitation of deep layer pyramidal neurons. The present study compared only distinct brain states, such as waking and sleep. However, the translaminar recordings combined with quantified spike transmission probability measures can be also used to monitor within and interlaminar interactions in visual or other cortical areas for studying cortical processing in perception, learning, and motor planning.

33A small subset of layer 6 neurons showed a robust inverse correlation with the activity of other principal cells and interneurons in all brain states, which was particularly striking during non-REM sleep. We named these cells DOWN state-active neurons because they are specifically and selectively active during the DOWN state of slow oscillations. The identity of the DOWN-state active neuron is not known at present. On the basis of their location and brain state correlates, they might correspond to a subset of GABAergic interneurons that express neuronal NOS (nNOS). These nNOS-expressing neurons increase their c-Fos expression during non-REM sleep and are distributed mostly in deep cortical layers. Optogenetic identification and manipulation of these neurons in future experiments may reveal their exact physiological contribution to slow oscillations and other network functions.

34Further information and requests for reagents and resource may be directed to, and will be fulfilled by the Lead Contact, Dr. György Buzsáki (gyorgy.buzsaki@nyumc.org).

35All animal handling procedures were approved by the Institutional Animal Care and Use Committee of New York University Medical Center. We used several lines of transgenic mice Mice (n = 19; n = 3 Ai35; n = 5 PV-Cre::Ai32; n = 9 VGAT-Cre::Ai32; n = 1 CaMKII-Cre::Ai35; n = 1 CaMKII-Cre::Ai32) for optogenetic tagging of the recorded units.

36A total of 19 male mice (28-35 gr, 3-8 months old) were implanted with recording electrodes under isoflurane anesthesia. These procedures were performed in two steps.

37In the first step, a ground electrode (100-μm diameter tungsten wire) was implanted in the contralateral cerebellum and head plate base was placed around the implantation target area (the primary visual cortex, V1). After 13 day of recovery from the first step, mice underwent the silicon probe implantation procedure. A single shank silicon probe (Cambridge NeuroTech H3 64x1 probe) was mounted on a movable microdrive for recording the activity of multiple single-units and local field potentials (LFPs) in V1. The high-density silicon probe had 64 recording sites (100k–1MΩ impedance each site), aligned on the linear edge of the probe (20-μm vertical separation). After the craniotomy above the target implantation site, the probe was implanted at anteroposterior: +1.0mm, mediolateral: +2.5 mm, with a 21° angle from dorso-ventral axis and a 10° angle from the mediolateral axis, so that the probe was perpendicular to the brain surface at the target site. The probe was lowered to 1.0 mm below the brain surface. For optogenetic tagging of specific cell types, an optic fiber was placed right above the skull over the implantation site. The back end of the fiber was coupled to a laser diode (450 nm blue, Osram).

38After the second operation (< 30 min), mice were allowed to recover overnight in their home cage before the recording session.

39We recorded from the mice while they slept or walked around freely moving in the home cage for of 68 hr. Electrophysiological data were acquired using an Intan RHD2000 system (Intan Technologies LLC) digitized with 20 kHz rate. The wide-band signal was down-sampled to 1.25 kHz and used as the LFP signal.

40For optogenetic tagging of specific neuron types, blue laser light pulses were delivered above the V1. The maximum light power at the tip of the optic fiber was 1 to 3 mW (450 nm, Osram Inc). 10 ms light pulses with 40%, 70% and 100% of the maximum power were delivered (n = 500 times at each intensity, 1 Hz). At the beginning or the end of the recording session, light from laser diode (520 nm, Osram) was delivered in the home cage to induce light-evoked responses of V1 neurons (n = 500 times, 400 ms duration pulses, 0.5 Hz).

41In four mice, a small current was delivered through two of the recording sites to produce an electrolytic lesion for subsequent histological verification of the recording depths. The mouse was anesthetized with isoflurane and placed in a stereotaxic apparatus. A thin metal bar (anal electrode) lubricated with vaseline was inserted in the anus of the animal and 4 to 10 μA current was applied between the anal electrode and target electrode channel on the silicon probe for 5-10 s per channel. After electrolytic lesioning procedure, mice were returned to their home cage. After two days, the animal was sacrificed.

42At the end of the recording session or two days after the electrolytic lesion, mice were overdosed with pentobarbital injection (100 mg/kg body weight), perfused with saline and 4% paraformaldehyde before their brains were rapidly removed. After overnight post fixation in 4% paraformaldehyde solution, the brain was washed in PBS three times. Coronal sections (50 μm) were cut on a vibratome (Leica, VT1000S) and brain slices were collected in PBS. To identify the borders between neocortical layers, fluorescent Nissl staining was performed. The procedure consisted of one time 10 min wash with PBS-0.1% Triton solution (PBS-T) for permeabilization, 20 min incubation in NeuroTrace 500/525 Green Fluorescent Nissl Stain diluted by 300-fold, two times of 10 min wash in PBS-T, three times of 5 min wash in PBS, and 2 hr wash in PBS. After the washing procedures, brain sections were mounted in Fluoromount with DAPI (Sigma) and imaged with a confocal laser-scanning microscope (Zeiss, LSM 700).

43Spike sorting was performed semi-automatically, using Kilosort. This was followed by manual adjustment of the waveform clusters using the software Phy. Following parameters were used for the Kilosort clustering.ops.Nfilt 6 ∗ numberChannelsops.nt0 64ops.whitening ‘full’ops.nSkipCov 1ops.whiteningRange 64ops.criterionNoiseChannels 0.00001ops.Nrank 3ops.nfullpasses 6ops.maxFR 20000ops.fshigh 300ops.ntbuff 64ops.scaleproc 200ops.Th [4 10 10]ops.lam [5 20 20]ops.nannealpasses 4ops.momentum 1./[20 800]ops.shuffle_clusters 1ops.mergeT 0.1ops.splitT 0.1ops.initialize ‘no’ops.spkTh 31

44ops.Nfilt 6 ∗ numberChannels

45ops.nt0 64

46ops.whitening ‘full’

47ops.nSkipCov 1

48ops.whiteningRange 64

49ops.criterionNoiseChannels 0.00001

50ops.Nrank 3

51ops.nfullpasses 6

52ops.maxFR 20000

53ops.fshigh 300

54ops.ntbuff 64

55ops.scaleproc 200

56ops.Th [4 10 10]

57ops.lam [5 20 20]

58ops.nannealpasses 4

59ops.momentum 1./[20 800]

60ops.shuffle_clusters 1

61ops.mergeT 0.1

62ops.splitT 0.1

63ops.initialize ‘no’

64ops.spkTh 31

65Cross-correlograms of spike trains of neuron pairs can reveal putative synaptic connections between them. This takes the form in the cross-correlogram of short time-lag (1–4 ms) positive or negative deviations from baseline indicating putative excitatory or inhibitory connections, respectively. Such detection is based on testing the null hypothesis of a homogeneous baseline at short time-scale. To this end, cross-correlograms binned in 0.5-ms windows were convolved with a 7-ms standard deviation Gaussian window resulting in a predictor of the baseline rate. At each time bin, the 99.9999 percentile of the cumulative Poisson distribution (at the predicted rate) was used as the statistical threshold for significant detection of outliers from baseline. A putative connection was considered significant when at least two consecutive bins in the cross-correlograms within +1.5 to +4 ms passed the statistical threshold. Corrected spike transmission probability was calculated as described in.

66Units with maximum waveform amplitude with positive sign were classified as ‘positive waveform units’, and those with negative sign as ‘negative waveform units’. Putative fiber potential was identified by the kurtosis of the normalized waveform (maximum amplitude = 1), calculated as the second derivative of negative waveform units. If the kurtosis was larger than ± 5SD of the mean, they were classified as putative fiber units. For the remaining negative waveform units, units with the trough-to-peak latency (TP latency) > 0.55 ms were tentatively classified as putative E cells and those with TP latency ≤ 0.55 ms were classified as putative I cells.

67To optogenetically tag PV cells and VGAT cells in V1, PV-Cre::Ai32 mice and VGAT-Cre::Ai32 mice were used, respectively. 10 ms light pulses were delivered every 1 s for 500 times. Light-triggered post-event histogram was binned at 1 ms width. If the peak firing rate in a 1-6 ms window after the light delivery (average of the first and second maximum value) was larger than the ± 8 SD of the mean of the baseline (–100-0 ms before the light delivery), the unit was defined to be optogenetically activated.

68Spectrograms for brain state scoring were constructed with a 1 s sliding 10 s window FFT of 1250 Hz LFP data at log-spaced frequencies between 1 and 100 Hz. Three types of signals were used to score state in our recordings: broadband LFP, narrowband theta frequency LFP, and electromyogram (EMG). For broadband LFP signal, principal components analysis (PCA) was applied to the z-transformed (1-100 Hz) spectrogram. The first PC in all cases was based on power in the low (< 20 Hz) frequency range and had oppositely weighted power in the (> 32 Hz) higher frequencies. Theta dominance was taken to be the ratio 5-10 Hz and 2-16 Hz power from the spectrogram. EMG was extracted from the intracranially recorded signals by detecting the zero time-lag correlation coefficients (r) between 300-600 Hz filtered signals (using a Butterworth filter at 300 – 600 Hz with filter shoulders spanning to 275 – 625 Hz) recorded at all sites. The state scoring algorithm was performed by a series of divisions with thresholds set at the trough between the peaks of distributions in these three metrics (https://github.com/buzsakilab/buzcode/tree/master/detectors/detectStates/SleepScoreMaster). First, non-REM time points were extracted via a threshold at the trough between the peaks in the broadband PC1 histogram. Among the remaining time points, REM was recognized by finding time points with both EMG values below the bimodal dip in that metric and theta in the upper mode of the distribution. Waking was defined as 7 min or longer arousals. After automated brain state scoring, all states were manually reviewed by the experimenter and minor corrections were made when discrepancies between automated scoring and user assessment occurred.

69Slow waves were detected using the coincidence of a two-stage threshold crossing in two signals: a drop in high gamma power (100-400 Hz, representative of spiking) and a peak in the delta-band filtered signal (0.5-8 Hz). The gamma power signal was smoothed using a sliding 80-ms window, and locally normalized to account for non-stationaries in the data. Two thresholds were used for event detection in each LFP-derived signal: a “peak threshold” and a “window threshold.” Time epochs in which the delta-filtered signal crossed the peak threshold were taken as putative slow wave events, with start and end times at the nearest crossing of the window threshold. To determine the delta threshold, all peaks in the delta-filtered signal greater than 0.25 standard deviations were detected as candidate delta peaks and binned by peak magnitude. The peri-event time histogram (PETH) for spikes from all cells was calculated around delta peaks in each magnitude bin and normalized by the mean rate in all bins. The smallest magnitude bin at which spiking (i.e., the PETH at time = 0) was lower than a set rate threshold was taken to be the peak threshold. The window threshold was set to the average delta value at which the rate crosses this threshold in all peak magnitude bins. The gamma power threshold was calculated similarly but using a drop below a gamma power magnitude instead of peaks above a delta magnitude. Once the thresholds were calculated, candidate events were then detected in the delta and gamma power signals, and further limited to a minimum duration of 40 ms. Slow wave events were then taken to be overlapping intervals of both the gamma and delta events.

70To calculate standardized depth coordinates across animals, we aligned depth profiles using electrophysiological landmarks. These included the largest amplitude peak of the depth profile of high-frequency LFP power (500 Hz - 5 kHz) corresponding to mid-layer 5 and four prominent sinks and sources from the averaged DOWN-UP CSD maps for each animal (see Figure 1). Depth profiles from the different animals were aligned according to these landmarks and the distance between each of them was “warped” (either extended or stretched) accordingly. Linear interpolation of inter-electrode distance was performed to obtain distance in μm. Estimated pia surface was 80 μm above landmark a.

71To separate the different sources that contribute to the LFP mixed signal, we employed independent component analysis (ICA) analysis as has been described and validated previously for hippocampal recordings. ICA is a blind source-separation technique that can isolate spatially segregated stable patterns of activity in a mixed signal recorded with an array of sensors. Applied to linear profiles of LFP it can separate physiologically meaningful sources. Here, we applied ICA to spatially contiguous LFP channels after filtering in the gamma band (30 – 100 Hz). Prior to application of the ICA algorithm, we performed a principal component analysis (PCA) reduction and maintained only the first ten PCs for subsequent ICA decomposition. The ICA algorithm (runinca; ) takes a time series of data with dimension equal to the number of recording sites, and returns a time series of the same dimensionality, but rotated so that each dimension represents a different IC. The inverse of the mixing matrix that transforms the LFP data into the ICs gives the channel weight of each component that is captured for each electrode. When projected to the anatomical location of the electrodes, this corresponds to the spatial voltage loadings of each IC. We ranked the components by the amount of variance they explain in the original data (relative power).

72Using coherence as similarity measure, an interaction-energy based clustering was implemented as previously done for hippocampal LFP7. Every site served as a reference against all the other referred sites. The resulting values were clustered using a gradient-descent algorithm, so that each site was merged with that cluster for which the resulting coherence gain after merging was the largest. Starting from random initial assignments, the clustering algorithm formed stable but fewer clusters corresponding to a local energy minimum. Energy of cluster A is defined as:Where, Cij is the coherence between ith and jth sites and NA is the number of recording sites in cluster A. The energy gap between two different assignments to cluster A and B of site i is:If the energy gap is positive, site i is moved into cluster B, otherwise it remains in cluster A. Since the method results only local minima of energy and stochastic components, such as the random initial condition and update order affects, clustering consistency was verified by repeating the process several times.

73To perform spectral analysis at a high resolution in time and frequency, the complex wavelet transform (CWT) of the LFP (or ICs) was calculated using complex Morlet wavelets. Wavelets were calculated using a logarithmically spaced frequency vector in the band of interest. Phase-amplitude cross-frequency coupling for a given LFP recording was assessed using the modulation index measure (MI; ). Phase time-series were binned into phase intervals and the mean wavelet amplitude was calculated for each of them. The MI was obtained by measuring the divergence of the observed amplitude distribution from the uniform distribution. The statistical significance of the MI values (p value) was assessed by a surrogate analysis (n = 1000 surrogates) with random shifts between the phase and amplitude time series. Mean-subtracted spectral analysis was obtained by calculating the mean power or coherence across all V1 channels for each session and subtracting it from each channel value. For the presented plots, grand averages were calculated as the mean across all animals.

74The phase-locking of spikes to LFP features at each frequency was measured for individual units using the wavelet phase from 30-100 Hz (20 logarithmically spaced wavelet scales) at the time of each spike9. Modulation indices were calculated using the mean resultant length of the phases, and significance was estimated using the Rayleigh test for non-uniformity (p < 0.05) using circular statistics. Preferred frequency of modulation was determined as the largest mean vector length of each significantly modulated neuron.

75To classify V1 units, we took advantage of the layer-specificity of gamma oscillations (Figure 2). For each unit, we first calculated their spike-phase modulation by the 6 gamma ICs. Next, we used the matrix of spike-IC phase coupling of all V1 units (excitatory and inhibitory) against all six ICs to perform k-Means clustering. The ‘cityblock’ metric was used as distance measure. We reduced the dimensionality of the data with t-Distributed Stochastic Neighbor Embedding (t-SNE) method by optimizing the Kullback-Leibler divergence of distributions between the original space and the embedded space, computing the Mahalanobis distance. We found the optimal number of clusters in the data to be eight by employing the silhouette criterion.

76Burstiness was determined by calculating the average number of spikes in the 1.5–13.5 ms bins of the spike autocorrelogram divided by the average number of spikes in the 200–300 ms bins.

77All statistical analyses were performed with standard MATLAB functions. No specific analysis to estimate minimal population sample were used, but the number of animals, trials, and recorded cells were larger or similar to those employed in previous works. Unless otherwise noted, for all tests, non-parametric two-tailed Wilcoxon rank-sum (equivalent to Mann-Whitney U-test), Wilcoxon signed-rank or Kruskal-Wallis one-way analysis of variance were used. For multiple comparisons, Tukey’s honesty post hoc test was employed. Boxplots represent median and 25th/75th percentiles and their whiskers the data range. In some of the plots outlier values were not shown but they were always included in the statistical analysis.

78LFP and spike data have been deposited in buzsakilab.com and are freely available for further analyses.