## Abstract

One of the most relevant questions regarding the function of the nervous system is how sensory information is represented in populations of cortical neurons. Despite its importance, the manner in which sensory-evoked activity propagates across neocortical layers and columns has yet not been fully characterized. In this study, we took advantage of the distinct organization of the rodent barrel cortex and recorded with multielectrode arrays simultaneously from up to 74 neurons localized in several functionally identified layers and columns of anesthetized adult Wistar rats in vivo. The flow of activity within neuronal populations was characterized by temporally precise spike sequences, which were repeatedly evoked by single-whisker stimulation. The majority of the spike sequences representing instantaneous responses were led by a subgroup of putative inhibitory neurons in the principal column at thalamo-recipient layers, thus revealing the presence of feedforward inhibition. However, later spike sequences were mainly led by infragranular excitatory neurons in neighboring columns. Although the starting point of the sequences was anatomically confined, their ending point was rather scattered, suggesting that the population responses are structurally dispersed. Our data show for the first time the simultaneous intra- and intercolumnar processing of information at high temporal resolution.

## Introduction

The structure of the cerebral cortex, divided into 6 layers and composed of a large set of columnar modules, is one of the most salient characteristics of the mammalian brain (Mountcastle 1997). Motivated by the search of a ‘canonical’ microcircuit common to all areas of the neocortex, many studies have been aimed to describe the flow of excitation along the layers of a single cortical column. Based on anatomical data and single-cell recordings, the flow of excitation was initially proposed to be sequential, following the layer 4 (L4) → layer 2/3 (L2/3) → layer 5 (L5) pathway (Gilbert and Wiesel 1979; Armstrong-James et al. 1992; Feldmeyer 2012). However, new lines of evidence advocate for a rather simultaneous activation of the granular and infragranular layers, presenting the activity peak in the interval of 10–20 ms after sensory stimulus—thus suggesting that the thalamus activates 2 separate, independent “strata” of the cortex in parallel (de Kock et al. 2007; Sakata and Harris 2009; Constantinople and Bruno 2013).

Further experiments performed in the barrel cortex have provided insights in how single stimuli are propagated along a set of cortical columns. The rodent barrel cortex offers the advantage that the sensory periphery is represented in a well-known topographic map, where each single whisker has its corresponding barrel-related column (Feldmeyer et al. 2013). By using large-scale multielectrode arrays (Petersen and Diamond 2000) and voltage-sensitive dye (VSD) recordings (Petersen et al. 2003), it has been shown that single-whisker deflection leads to an initially confined activation of the respective barrel column after ∼7 ms, followed by a horizontal spread across the entire barrel field within ∼50 ms.

Despite the results presented in these previous reports, it is not yet known how neuronal activity propagates within the cortex across both layers and columns, and how the neuronal populations represent sensory stimuli in a collective manner. To address these questions, we recorded with multielectrode arrays the sensory-evoked responses of neuronal populations distributed across several layers of up to 4 columns in the barrel cortex of anesthetized adult rats in vivo (Csicsvari et al. 2003). In our datasets containing a total of 494 neurons, we segregated the recorded neurons into putative inhibitory (INH) and excitatory (EXC) cells, subdividing them further according to their laminar and columnar location (Yang et al. 2013). A well-established template-matching algorithm was subsequently applied to identify temporally precise repetitive spike sequences in these datasets, which have been proposed to be a substrate of both information transfer and transformation in cortical networks (Abeles and Gerstein 1988; Nadasdy 2000; Izhikevich et al. 2004; Sun et al. 2010).

With this approach, we found that although sensory-evoked activity was sparse and heterogeneous at the single neuron level, temporally precise repetitive spike sequences were consistently elicited by single-whisker stimulation. Those sequences which started shortly after stimulation were mainly led by a subgroup of putative INH neurons located at thalamo-recipient layers (4 and 5B) of the principal whisker column, thus constituting the first manifestation of feedforward inhibition at the level of large-scale cortical networks in vivo. The starting point of the sequences evolved subsequently toward infragranular layers of the neighboring columns (NCs), with later sequences mainly led by EXC neurons in infragranular layers. Although the starting point of the sequences was relatively well confined to specific cortical modules, the ending point was scattered across them, suggesting that the flow of cortical activity is anatomically dispersed, and therefore an intra- and intercolumnar processing of information is simultaneously performed. In this regard, we also show that specific L5 cell types might play a major role in synchronizing the intracolumnar network activity during the late phases of the evoked responses. In summary, our data provide for the first time a macroscopic view of the flow of activity within the barrel cortex at high temporal resolution, thus providing new insights into the possible mechanisms how sensory stimuli are represented in neuronal populations.

## Materials and Methods

### Surgical Preparation

All procedures were approved by the local ethics committee (#23177-07/G10-1-010) and followed the European and German national regulations (European Communities Council Directive, 86/609/ECC). The experiments were conducted in male Wistar rats, postnatal day 28–42 (80–200 g). Animals were anesthetized using urethane (1.5 g/kg) and head-fixed into a modified stereotaxic device using a metal chamber and dental cement. Lidocaine (4%, 0.5 mg) was administered at the site of incision. Depth of anesthesia was monitored by assessing eyelid and hindpaw reflexes and maintained at stage III/3–4 (Friedberg et al. 1999) by supplementary injections of urethane (10% of initial dose) (see below). A 3 × 3 mm2 craniotomy was performed over the barrel cortex of the left hemisphere (centered at 2.5 mm posterior and 5.5 lateral to Bregma). Extreme care was taken at all times not to damage the cortex, especially during the removal of the dura. One silver wire was inserted into the cerebellum, which served as a ground electrode. Body temperature was maintained at ∼37°C with a feedback temperature controller (World Precision Instruments, Sarasota, FL, USA). Ophtalmic ointment was placed on both eyes to prevent corneal drying. This experimental setup allowed VSD imaging and multielectrode electrophysiological recordings for up to 8 h.

The anesthesia level was post hoc confirmed by the most prevalent frequency of the recorded superficial local field potentials (LFPs). As previously shown (Friedberg et al. 1999), the stages of anesthesia are best subdivided using the dominant (most prevalent) frequency of the electrocorticogram (ECoG). As confirmed in 4 control experiments, the most superficial LFP was always almost identical to the ECoG recorded at the surface of the dura with a platinum ball electrode (data not shown) (Buzsaki et al. 2012). For this computation, the superficial LFP signal was filtered from 0.1 to 32 Hz so that our results could be directly compared with those previously described (Friedberg et al. 1999). The dominant frequency of the LFP was defined as the mean one showing the maximum power in the 4 s periods previous to whisker stimulation (see below). The majority of the animals (14 of 18, 77.8%) were generally at III-4 stage (dominant LFP frequency in the range of 0.5–2.5 Hz) and a minority (4 animals, 22.2%) at III-3 stage (in the range of 2.5–4.5 Hz). No animal presented a mean dominant frequency above or below these 2 ranges. In conclusion, our animals were commonly at rather medium and deep anesthesia stages (III-3 and III-4).

### Whisker Stimulation

Prior to vibrissae stimulation, all whiskers were trimmed to 8 mm length. Individual whiskers were briefly deflected in the vertical direction at 5 mm from their base by inserting each of them 3 mm into a capillary tube glued to a piezoelectric bimorph actuator (Physik Instrumente, Karlsruhe, Germany) that was controlled by a voltage pulse generator (Master-8, A.M.P.I., Jerusalem, Israel) (protocol similar to Petersen et al. 2003). The voltage pulse was an up-down square step function of 2 ms duration. The produced capillary movement had an amplitude of 150 µm and a rise time of 1.5 ms (measured at the capillary ending). In order to ensure isolated mechanical stimulation of the target whiskers, the remaining whiskers were completely trimmed before starting the recordings. Stimuli were applied to individual whiskers in blocks, each containing typically 200 trials (sweeps) with an intertrial interval of 5 s (occasionally 10, 15, or 30 s). This was repeated for each of the selected target whiskers (maximum 4). Afterward, a period of spontaneous activity was recorded for 500 to 2000 s.

### VSD Imaging

For VSD imaging, the VSD RH1691 (Optical Imaging, Rehovot, Israel) was dissolved at 1 mg/mL in a saline solution containing (in mM): 125 NaCl, 2.5 KCl, and 10 HEPES (pH 7.3 with NaOH). 150 µL of dye solution was topically applied to the exposed cortex and allowed to diffuse into the cortex over 1 h period. Afterward, the unbound dye was washed out for 15 min, the cortex was covered with 1% agar dissolved in Ringer's solution, and a glass coverslip was placed on top. VSD signals were collected using the same setup as previously described (Yang et al. 2013). Images from a field of view of 2.6 × 2.6 mm2 (100 × 100 pixels) were recorded with an unbinned whole frame sampling frequency of 500 Hz. Evoked activity following single-whisker stimulation was imaged in 2048 ms long imaging sweeps. Whisker deflection was produced 296 ms after the start of each sweep. Sequences of typically 6 sweeps were recorded for each whisker.

Images were subsequently displayed and analyzed for columnar localization using the MiCAM image-analysis software (BV-Analizer ver. 0909, BrainVision, Tokyo, Japan). First, the sequences related to a single whisker were averaged. Images were additionally smoothed by a 5 × 5 pixel spatial filter. To attenuate the slow signal components produced from bleaching of the VSD (Carlson and Coulter 2008), we afterward applied a weak high-pass filter of the MiCAM software (time constant τ = 66.67 ms) to the imaging data. From these signals, the fractional change in fluorescence (ΔF/F0) was represented for each pixel as the change of fluorescence intensity (ΔF) divided by the prestimulus fluorescence intensity (F0) in the same pixel. The time point for F0 was set to 4 ms before whisker stimulation. Only fluorescence changes with a maximal ΔF/F0 of at least 0.1% were considered as evoked responses. The frame corresponding to 10 ms after whisker deflection, which has been previously shown to be the initial source of neuronal activity in the barrel cortex (Petersen et al. 2003), was used to localize the barrel column associated to the stimulated whisker. The same time point was used to generate the contour plots representing the area of evoked responses in Figure 1C1. These contour plots included pixels with a value higher than 0.8 times the maximal ΔF/F0 amplitude contained in the same frame.

Figure 1.

Large-scale evoked responses acquired using VSD imaging and an 8 × 16 channels Michigan probe. (A) Schematic illustration of the experimental setup for selective mechanical stimulation of single whiskers and simultaneous multielectrode recordings in 3 barrel-related columns. (B) Illustration of the 8 × 16 channels probe showing shank and site spacing. Position of the barrels D1, C1, and B1, and the location of the different layers were estimated from the data presented in panels DF. (C1), Mean VSD responses averaged from 6 deflections of whiskers C1 and D1. These data were used to localize the barrels before probe insertion. White stars mark the insertion position of the 8 probe shanks. (C2), Photomicrograph of tangential section through L4 of the barrel field after processing the tissue with CO. Yellow points indicate the location of the 8 shanks, which were labeled with DiI before insertion. (D) Mean local field potentials (LFPs) and corresponding current source density (CSD) maps elicited by B1-whisker deflection averaged from 200 trials. Average LFP signals are plotted from 0 to 30 ms after stimulus at the location of the corresponding electrodes (see panel B). In the CSD maps, current sources and sinks are represented by red (positive) and blue (negative) colors, respectively (see color bar at the right side). Position of the D1, C1, and B1 barrel-related columns are indicated at the upper part. Location of the cortical layers and depth are represented at the left side. Gray lines represent the estimated separation of the cortical layers. (E and F) Similar as in panel D, but elicited by C1- and D1-whisker deflections, respectively.

Figure 1.

Large-scale evoked responses acquired using VSD imaging and an 8 × 16 channels Michigan probe. (A) Schematic illustration of the experimental setup for selective mechanical stimulation of single whiskers and simultaneous multielectrode recordings in 3 barrel-related columns. (B) Illustration of the 8 × 16 channels probe showing shank and site spacing. Position of the barrels D1, C1, and B1, and the location of the different layers were estimated from the data presented in panels DF. (C1), Mean VSD responses averaged from 6 deflections of whiskers C1 and D1. These data were used to localize the barrels before probe insertion. White stars mark the insertion position of the 8 probe shanks. (C2), Photomicrograph of tangential section through L4 of the barrel field after processing the tissue with CO. Yellow points indicate the location of the 8 shanks, which were labeled with DiI before insertion. (D) Mean local field potentials (LFPs) and corresponding current source density (CSD) maps elicited by B1-whisker deflection averaged from 200 trials. Average LFP signals are plotted from 0 to 30 ms after stimulus at the location of the corresponding electrodes (see panel B). In the CSD maps, current sources and sinks are represented by red (positive) and blue (negative) colors, respectively (see color bar at the right side). Position of the D1, C1, and B1 barrel-related columns are indicated at the upper part. Location of the cortical layers and depth are represented at the left side. Gray lines represent the estimated separation of the cortical layers. (E and F) Similar as in panel D, but elicited by C1- and D1-whisker deflections, respectively.

### Multielectrode Recordings

Neuronal activity was recorded with a 16- or 128-channel “silicon probe” inserted perpendicularly into the barrel cortex (NeuroNexus Technologies, Ann Arbor, MI, USA). The one-shank 16-channel probe presented a separation of 100 µm between recording sites (177 or 413 µm2 diameter) (Supplementary Fig. 1). The 128-channel probe contained 8 shanks separated by 200 µm, presenting a vertical spacing of 75 µm between the recording sites (413 µm2 diameter) (Fig. 1B).

The insertion points of the probe shanks were determined by the previous VSD imaging session and the position of the large macrovessels at the cortex surface. Probe insertion across or near the macrovessels was avoided; otherwise, their damage would result in bleeding, which subsequently led to the abortion of the experiment. In this regard, please note that the arteriolar and venular domains at the cortex surface do not co-localize with the underlying barrels, thus showing a large variability across animals (Woolsey et al. 1996; Blinder et al. 2013). For histological verification of tracks, the probes were labeled with DiI (1,1′-dioctadecyl-3,3,3′,3′-tetramethyl indocarbocyanine, Molecular Probes, Eugene, OR, USA) before insertion.

All data were continuously digitized at 20 kHz and stored for offline analysis using a multichannel extracellular recording system and the MC_RACK software (Multi Channel Systems, Reutlingen, Germany). The total duration of the recorded data varied between 45 min and 2 h.

### Histology

After the experiment, animals were deeply anesthetized with ketamine (120 mg/kg, ketamine, 50 mg/mL, Hameln Pharma, Hameln, Germany) and xylazine (5 mg/kg, Rompun 2%, Bayer, Leverkusen, Germany) and perfused through the aorta, first with 10 IU/mL heparin in phosphate-buffered saline (PBS), second with 4% paraformaldehyde. Subsequently, the brain was carefully dissected out of the skull. In order to visualize the entire barrel field representing the contralateral whiskers, the barrel cortex was tangentially cut and flattened, from ∼2 to 1 mm. After 24 h of ﬁxation in 4% paraformaldehyde, the cortex was rinsed in 0.1 M phosphate buffer (PB) and incubated in 30% sucrose (in PB) overnight. After washing with PB, the cortex was sectioned in 80-µm-thick sections and processed for cytochrome oxidase (CO) histochemistry. Sections were rinsed in PB, then incubated at 39°C in a solution of 0.6 mg cytochrome c, 0.5 mg DAB, and 44 mg saccharose per mL, with 0.3% catalase included (all from Sigma, Deissenhofen, Germany). Sufficient staining was achieved after 2–7 h of reaction. Finally, sections were intensified with 0.5% copper(II)sulfate (Sigma) for 2–3 min, air-dried, and mounted.

### Data Analysis

Unless otherwise stated, all analyses were performed offline using custom programs written in Matlab (Mathworks, Natick, MA, USA) and C. First, the raw files (digitized at 20 kHz) were accessed using functions included in the FIND toolbox (Meier et al. 2008). Signals were used at full sampling rate for spike detection and sorting (see below); further, they were downsampled to 1 kHz and re-stored for posterior analyses of the LFPs.

### Spike Detection and Sorting

The continuously recorded wide-band signals were high-pass filtered (0.8–5 kHz). Groups of 2 to 4 channels were selected and defined as “virtual tetrodes,” that is, groups of channels potentially recording from the same neurons (Harris et al. 2000; Einevoll et al. 2012). These groups of channels were nonoverlapping, that is, individual channels belonged only to a single group. Additionally and in order to ensure that neurons were present only in one group, at least one channel was left unused between neighboring groups. Spike detection was performed in each group of channels independently using amplitude-thresholding in the negative range, similarly to classical stereotrode- or tetrode-based techniques (Gray et al. 1995). Threshold level was set independently for each channel as −7.5 times the standard deviation (SD) of the signal. When a threshold crossing was detected in either of the channels within a group, all sampled amplitude values from all channels in the time range from −0.5 to +0.5 ms relative to the waveform negative peak were extracted. Feature vectors were then computed from the spike waveforms. These feature vectors contained 3 values for each channel, namely the negative peak amplitude plus the 2 first principal components derived from the waveforms. The (n × 3)-dimensional vectors (where n represents the number of channels within a group) were afterward sorted using KlustaKwik (http://klustakwik.sourceforge.net) and Klusters (http://klustakwik.sourceforge.net) (Harris et al. 2000; Hazan et al. 2006).

Several criteria were established in order to ensure the isolation quality of the sorted neurons (Hill et al. 2011). First, <1.5% of the spikes should occur within 2 ms of each other, thus accounting for a small number of refractory period violations (mean ± standard error of the mean [SEM] = 0.38 ± 0.02% averaged across the validated neurons). Second, the spontaneous firing rate (FR) of the isolated cells should be stable; therefore, the coefficient of variation (CV) of the FR computed over the whole recording period using 100 s bins had to be <1.5 (mean 0.7 ± 0.2 across validated neurons). Third, a minimum value of 10 was established for the “isolation distance” so that neurons were accepted for further analysis (mean 45.4 ± 2.4 across validated neurons) (Schmitzer-Torbert et al. 2005). In order to confirm our sorting validation criteria, spike auto-correlograms and multichannel spike waveforms were routinely inspected (Harris et al. 2000; Hill et al. 2011). A rather conservative criterion was used when accepting sorted units for further processing, thus only including clearly isolated neurons. This criterion was evident since (a) the number of single units was always less than a half than the number of channels for both probe configurations; and (b) the spike amplitudes ranged between ∼70 and 200 µV across all neuronal groups (see Results).

Cells were subsequently classified as putative INH and EXC neurons based on the mean spike waveform (Fig. 2B) (Sirota et al. 2008; Sakata and Harris 2009). For each neuron, 2 parameters were calculated, which have been shown to reliably separate between several identified neuronal populations (Royer et al. 2012): 1) the trough to right (late) peak latency (related to the repolarization of the intracellular action potential) and 2) the asymmetry index (possibly reflecting differences in the rate of fall of spike repolarization). Because of the clear separation observed between the 2 groups, we separated them by the line y = 2.8 × − 1.08 (see Fig. 2C).

Figure 2.

Multichannel spike sorting. (A) Left part, high-pass filtered signals (0.8–5 kHz) from 5 channels recorded with a 1 × 16 channels Michigan probe showing single-neuron spontaneous activity during 60 ms (data from same recording shown in Supplementary Fig. 1A). Stars mark the position where spikes from an excitatory (EXC, panel B2) and an inhibitory (INH, panel B3) neuron are detected and sorted. Right part, average spike waveform profiles of the exemplary EXC and INH neurons. (B1) Schematic drawing illustrating the 3 features extracted from the spike waveforms to perform neural classification: a and b, left (early) and right (late) baseline-to-peak amplitudes; c, trough to right peak latency. The first 2 features were used to calculate the peak amplitude asymmetry [(ba)/(b + a)]. (B2 and B3) Waveforms extracted from the exemplary EXC and INH neurons marked in A. (C) Separation of the 2 neuronal groups (putative EXC and INH neurons) based on the 2 parameters constructed as described in B1. Circles mark the exemplary EXC and INH neurons from panels B2 and B3.

Figure 2.

Multichannel spike sorting. (A) Left part, high-pass filtered signals (0.8–5 kHz) from 5 channels recorded with a 1 × 16 channels Michigan probe showing single-neuron spontaneous activity during 60 ms (data from same recording shown in Supplementary Fig. 1A). Stars mark the position where spikes from an excitatory (EXC, panel B2) and an inhibitory (INH, panel B3) neuron are detected and sorted. Right part, average spike waveform profiles of the exemplary EXC and INH neurons. (B1) Schematic drawing illustrating the 3 features extracted from the spike waveforms to perform neural classification: a and b, left (early) and right (late) baseline-to-peak amplitudes; c, trough to right peak latency. The first 2 features were used to calculate the peak amplitude asymmetry [(ba)/(b + a)]. (B2 and B3) Waveforms extracted from the exemplary EXC and INH neurons marked in A. (C) Separation of the 2 neuronal groups (putative EXC and INH neurons) based on the 2 parameters constructed as described in B1. Circles mark the exemplary EXC and INH neurons from panels B2 and B3.

### Current Source Density Maps

The cortical depth of the individual channels was assessed from the stereotaxically estimated depth of the probe tip and the vertical CSD profiles computed from the LFPs (see below). We used the early CSD sinks present at the thalamo-recipient L4 and L5B after sensory stimulation in order to assign the individual channels to specific cortical layers (see Fig. 1DF and Supplementary Fig. 1) (Mitzdorf 1985). The following procedure was used in each probe shank for layer identification. First, we established the upper border of the L4 sink as reference depth (500 µm). L4 thickness was determined by the upper and lower borders of its corresponding sink; similarly, the beginning of L5B was set as the upper border of the deep sink (see Supplementary Fig. 1). With this approach, the border between L2/3 and L4 was set as 500 µm depth; L5A started at an average depth of 798.6 ± 11 µm, and L5B at 1136.1 ± 12.7 µm (values averaged across experiments and given as mean ± SEM). Further, the somatic location of the spike-sorted neurons was estimated as the recording site containing the mean waveform with maximum negative peak amplitude.

CSD maps were computed using previously described methods using the averaged LFPs from up to 200 trials (Nicholson and Freeman 1975). In order to obtain the same number of CSD profiles and LFP signals, we first duplicated the uppermost and lowermost channels in each shank. Afterward, we smoothed the LFPs across spatially neighboring channels to reduce high spatial-frequency noise:

$ϕ¯(r)=14(ϕ(r+h)+2ϕ(r)+ϕ(r−h)),$
where $ϕ(r)$ is the LFP at depth r, and h is the sampling interval (100 or 75 µm depending on the probe) (Sakata and Harris 2009). Then, the second derivative was computed:
$D=1h2(ϕ¯(r+h)−2ϕ¯(r)+ϕ¯(r−h)).$

Data were finally interpolated and plotted as pseudocolor images, with current sources and sinks represented by red (positive) and blue (negative) colors, respectively.

In order to address the question whether septal-related areas could be characterized, we further computed CSD maps across the electrodes belonging to the same row within each shank (thus showing the horizontal CSD profiles rather than the vertical ones) (Supplementary Fig. 2). Even those few shanks (5 of 80), which were post hoc localized to septal-related areas in the histological sections, did not show any differences in terms of LFP or single-neuron related activity. Further, no differences were found between neurons recorded with shanks located near the center of the barrel, and those at the barrel periphery (Kerr et al. 2007)—suggesting that our methodology do not reach the spatial resolution needed to reliably characterize septal-related areas (see Discussion).

### Spike Train Analysis

Spontaneous activity was estimated for every neuron by computing their FR during the whole recording time, including the period of spontaneous activity (500–2000 s) recorded at the end of the experiment (Supplementary Fig. 3). The time windows from whisker stimulus to 1 s afterward were excluded from this analysis. Evoked activity was computed using the time window from 0 to 30 ms after stimulus. Further, the 100-ms window before whisker deflection was used to compare the post- to prestimulus activity in individual neurons (Figs 3B and 10B).

Figure 3.

Cellular response types to single-whisker stimulation in barrel cortex. (A) Typical responses of one inhibitory (INH) and 3 excitatory (EXC) neurons recorded in parallel. Each horizontal subplot represents the activity of an individual cell containing the raster plot aligned to whisker stimulation and the PSTH. Trials in the raster plot are presented in the same sequential order as recorded. Note that neurons upon whisker stimulation decreased their firing rate (FR) (ON cells), showed no change in FR (NR) or decreased their FR (OFF cells). (B) Comparison of mean FR of EXC (triangles) and INH neurons (circles) before versus after whisker stimulation, plotted in logarithmic scale. Prestimulus time is defined as the 100 ms period before whisker stimulation and poststimulus time as 0–30 ms afterward. Empty symbols represent neurons showing no significantly different activity in both periods, filled symbols neurons showing significantly different activity (P < 0.05). The neurons displayed in panel A are marked by their respective symbols.

Figure 3.

Cellular response types to single-whisker stimulation in barrel cortex. (A) Typical responses of one inhibitory (INH) and 3 excitatory (EXC) neurons recorded in parallel. Each horizontal subplot represents the activity of an individual cell containing the raster plot aligned to whisker stimulation and the PSTH. Trials in the raster plot are presented in the same sequential order as recorded. Note that neurons upon whisker stimulation decreased their firing rate (FR) (ON cells), showed no change in FR (NR) or decreased their FR (OFF cells). (B) Comparison of mean FR of EXC (triangles) and INH neurons (circles) before versus after whisker stimulation, plotted in logarithmic scale. Prestimulus time is defined as the 100 ms period before whisker stimulation and poststimulus time as 0–30 ms afterward. Empty symbols represent neurons showing no significantly different activity in both periods, filled symbols neurons showing significantly different activity (P < 0.05). The neurons displayed in panel A are marked by their respective symbols.

To display the overall changes in the activity of individual neurons, peristimulus time histograms (PSTHs) were computed using a time resolution of 5 ms (Fig. 3A) or 1 ms (Fig. 4A1–2). To display the mean ensemble activity of neuronal populations, mean spike density functions (SDFs) were computed using a Gaussian kernel of 1 ms SD (Fig. 10C). These parameters were chosen because they provided sufficient averaging of the neuronal discharges without losing important details.

Figure 4.

Responses of inhibitory (INH) and excitatory (EXC) neurons to principal whisker (PW) and neighboring whisker (NW) stimulation. (A1) Typical responses to PW deflection of 2 INH and 4 EXC neurons recorded simultaneously. Each of the neurons is represented individually in a horizontal subplot. Thick lines represent the PSTH of INH (EXC) cells. Activity is plotted from 10 ms before to 30 ms after whisker stimulus. Otherwise, same conventions are used as in Figure 3A. (A2) Activity of the same group of neurons evoked by NW deflection. (B1 and B2) Statistical analysis of first-spike latencies for the different putative neuronal groups. Each box plot represents median, interquartile, and range of latencies across all recorded neurons. Empty circles represent outliers. (C1 and C2) Precision (defined as the standard deviation of the first-spike latency in ms) for the different cell types after PW/NW deflection. Lower values represent neuronal groups with higher precision. Each bar represents mean and SEM values. (D1 and D2) Number of evoked spikes after PW/NW deflection (0–30 ms after stimulus) for the different cell types.

Figure 4.

Responses of inhibitory (INH) and excitatory (EXC) neurons to principal whisker (PW) and neighboring whisker (NW) stimulation. (A1) Typical responses to PW deflection of 2 INH and 4 EXC neurons recorded simultaneously. Each of the neurons is represented individually in a horizontal subplot. Thick lines represent the PSTH of INH (EXC) cells. Activity is plotted from 10 ms before to 30 ms after whisker stimulus. Otherwise, same conventions are used as in Figure 3A. (A2) Activity of the same group of neurons evoked by NW deflection. (B1 and B2) Statistical analysis of first-spike latencies for the different putative neuronal groups. Each box plot represents median, interquartile, and range of latencies across all recorded neurons. Empty circles represent outliers. (C1 and C2) Precision (defined as the standard deviation of the first-spike latency in ms) for the different cell types after PW/NW deflection. Lower values represent neuronal groups with higher precision. Each bar represents mean and SEM values. (D1 and D2) Number of evoked spikes after PW/NW deflection (0–30 ms after stimulus) for the different cell types.

### Temporally Precise Spike Sequences

Repetitive spike sequences were processed as described previously (Sun et al. 2010). In brief, repetitive spike sequences were detected in spike datasets by the use of a template-matching algorithm (Abeles and Gerstein 1988; Nadasdy et al. 1999; Rolston et al. 2007). Template spike sequences were generated using every spike in the dataset as starting point, each of them containing a sequence within an interval of 100 ms. These templates were compared with all subsequent spike sequences, which were generated using all following spikes as starting point and the same interval (100 ms). If 3 or more spikes were identical (with a precision of 2 ms) between the template and a subsequent spike sequence, then the matched elements were regarded as a candidate sequence. Thus, a candidate sequence should repeat at least 2 times in the spike dataset. Matches containing spikes from only one neuron, or with a delay of >10 ms between 2 sequential spikes, were excluded. Unless otherwise stated, a temporal precision of 2 ms was used in all searches because it gave enough resolution to account for the observed differences in neuronal latencies (see Fig. 4B1–2). This temporal precision is also reasonable in the light of recent findings, which prove by means of intracortical stimulation that differences in the timing of auditory cortex activity as small as 3 ms can guide decisions (Yang et al. 2008). Due to memory and computing time limitations, and because the present study was focused on sensory-evoked activity, we limited the searching algorithms to the time periods from 0 to 300 ms after whisker stimulation when the total number of neurons was >20.

Further, the statistical signiﬁcance of each sequence was evaluated by comparing its occurrence with that inside 2 surrogate sets generated by spike shufﬂing algorithms (Grun 2009). Only if the detected sequence occurred significantly more often in the original dataset than in any of the 2 surrogate sets (z-test, P < 0.001), it was accepted as a repetitive spike sequence. For the generation of each surrogate set a different procedure for spike shufﬂing was used, time jittering and neuron jittering—each of them produced from 100 shuffling trials (Sun et al. 2010). For time jittering, each spike was shifted by a random time interval following a Gaussian distribution with a mean of 0 and a SD of 5 ms. For neuron jittering, each spike was reassigned to a new neuron by a random step in the neuron index. The random step was given by a Gaussian distribution with a mean of 0 and a SD of 1. An additional control surrogate dataset generated by trial shuffling was taken into consideration (Supplementary Fig. 4).

Additional control analyses were performed by characterizing the properties of “chance sequences” detected within surrogate data. Twenty surrogates were created for each original dataset by spike shuffling (5 ms SD), thus keeping the mean spike timings (i.e., the mean PSTH) of EXC and INH neurons similar [data not shown, but see Sun et al. (2010) for a demonstration that the shuffling led to slight shifting in the spike distributions]. Note that in contrast to the procedure used to detect real sequences (see above), chance sequences were detected within the surrogates without employing any method for statistical validation (Matsumoto et al. 2013).

### Statistics

The following approach was generally used to test the hypothesis that 2 samples represent similar distributions. First each distribution of data was tested for normality (D'Agostino-Pearson test) within each group. If the size of both samples was equal or higher than 20 and both samples were normally distributed, parametric tests (independent or paired t-test) were applied at the 5% significance level for between-group comparisons. Otherwise, nonparametric tests (Mann–Whitney U or paired Wilcoxon signed-rank test) were performed (again at the 5% significance level).

One-way or two-way ANOVA tests were used to test if there was any significant difference between 3 or more neuronal groups (5% significance level). In case there was a significant difference, post hoc Tukey's least significant difference (LSD) tests were employed in order to perform pairwise comparisons between the individual groups of data.

Values throughout this report were generally given as mean ± SEM. However, the median and interquartile range was used to present the data related to neural latencies—thus avoiding the bias introduced by these time-based values, which generally present skewed distributions.

## Results

We recorded with Michigan-type electrode arrays LFPs and multiple single neurons in identified barrel-related columns of anesthetized adult Wistar rats in vivo (Fig. 1A,B). Electrodes were inserted orthogonally into the barrel cortex after identification of the barrel-related columns of interest by VSD imaging (Fig. 1C1). The electrode position was verified post hoc by visualizing the DiI traces of the electrode shanks in cytochrome-oxidase-stained histological sections (Fig. 1C2). The LFP responses upon stimulation of up to 4 neighboring whiskers revealed a barrel- and layer-specific activation of the somatosensory cortex (Fig. 1D–F). We used the resulting current source density (CSD) depth profiles to identify the cortical layers (Fig. 1D–F and Supplementary Fig. 1), in addition to the stereotaxically estimated depth of the recording channels. L4 was identified by sinks occurring 6–8 ms after whisker stimulation in the CSD (Mitzdorf 1985). Localized sinks presenting similar latencies defined the border between L5A and L5B. In accordance with previous reports, the activity propagated subsequently toward L2/3 and infragranular layers of the principal barrel-related column, where maximal CSD amplitudes were observed after 10–15 ms (Mitzdorf 1985). In addition, the activity propagated from the activated column to first- and even second-order NCs between ∼12 and ∼25 ms after whisker stimulus (Fig. 1D–F). These results demonstrate the validity of our procedure, which allowed us to allocate the recorded neurons to specific cortical columns and layers.

Using a multichannel spike-sorting algorithm, we isolated a total of 494 neurons in the recordings from 18 animals included in the present study. In 8 animals, a 1 × 16-channels probe was used, which allowed simultaneous recordings from 5 to 9 neurons (median 6) in a single barrel-related column. An 8 × 16-channels probe was used in the remaining 10 animals, which allowed simultaneous recordings from 21 to 74 neurons (median 47) in 2 to 4 barrel-related columns (median 4). In general, we preferred to cover barrels in a row (6 animals) rather than in an arc (4 animals), unless the borders of the craniotomy or the position of the macrovessels at the cortex surface impeded a safe probe insertion (see Materials and Methods).

The isolated single neurons were further associated to the corresponding recording site in which they showed their maximum waveform amplitude (Fig. 2A). Putative INH and EXC neurons were identified according to waveform asymmetry and spike width (Fig. 2B,C) (Sakata and Harris 2009). From the 494 recorded neurons, 85.4% were classified as putative EXC (n = 422) neurons and 14.6% as putative INH neurons (n = 72).

The laminar distribution of the recorded neurons was nonuniform. Of the recorded 494 neurons, 7.1% (4 INH, 31 EXC) were located in L2/3, 15.6% (19 INH, 58 EXC) in L4, 42.5% (16 INH, 194 EXC) in L5A, and 34.8% (33 INH, 139 EXC) in L5B. In general, signals originating from L2/3 neurons showed small spike amplitudes, which hampered the isolation of single neurons from the multiunit activity. Both L2/3 INH (68.4 ± 5.02 µV, n = 4) and L2/3 EXC neurons (116.3 ± 10.9 µV, n = 31) showed significantly (P < 0.01 in all cases) lower spike amplitudes than L5A EXC (167.9 ± 7.8 µV, n = 194) and L5B EXC (205.3 ± 12.5 µV, n = 139) neurons. Further, L2/3 neurons fired very sparsely, a property that imposed difficulties to the sorting algorithms (Harris et al. 2000).

Across all cortical layers INH neurons showed higher spontaneous activity than EXC ones (Supplementary Fig. 3). The spontaneous FR was higher in L2/3 INH (0.8 ± 0.4 spikes per second [spks/s]) than in L2/3 EXC (0.3 ± 0.05 spks/s) neurons, in L4 INH (1.3 ± 0.2 spks/s) than in L4 EXC (0.6 ± 0.07 spks/s) neurons, in L5A INH (2.8 ± 1.02 spks/s) than in L5A EXC (1.7 ± 0.09 spks/s) neurons, and in L5B INH (2.07 ± 0.3 spks/s) than in L5B EXC (1.09 ± 0.1 spks/s) neurons. Accordingly, a two-way ANOVA indicated a significant effect of the layer (P < 0.001) and the type of neuron (INH/EXC) (P < 0.001) on the spontaneous FR, but not of their interaction (P = 0.8, n.s.). In conclusion, spontaneous FR was consistently higher in INH than EXC neurons, and it was highest in L5A and lowest in L2/3 neurons.

### Single Cells Responded Heterogeneously and Sparsely to Sensory Stimulation

In order to characterize the neuronal responses related to sensory stimulation, we first analyzed single-cell responses to repeated principal whisker deflections by computing raster plots and PSTHs in the time period between 100 ms before and 200 ms after stimulus (Fig. 3A). In agreement with previous reports (Simons and Carvell 1989; Wilent and Contreras 2004; de Kock et al. 2007), most neurons showed clear responses within 5–30 ms after whisker stimulation. While the majority of neurons increased the FR after the stimulus (classified as ON cells), in a considerable fraction of neurons, the FR did not change (nonresponding [NR] cells). In addition, we found some neurons that showed a clear reduction of their FR upon whisker stimulation (OFF cells).

To quantify these responses we compared the FR in 2 time windows, the 100 ms before to the 30 ms after whisker stimulus (Fig. 3B). Of the 263 EXC neurons that were located in the column of the stimulated principal whisker, 60.1% were classified as ON, 33.5% as NR, and 6.5% as OFF cells. Further, the proportion of INH cells was significantly different (χ2 test, P < 0.01) across the 3 groups (ON, NR, and OFF cells). The vast majority (83.6%) of the 55 INH neurons were ON cells, 14.5% NR, and only one (1.8%) was classified as OFF. These data indicate that neurons were heterogeneous in relation to their responses to whisker stimulation, and that INH cells formed a more homogeneous group, since they presented more often ON responses than EXC ones.

An interesting group was the one constituted by neurons showing OFF responses, which so far have been seldom characterized in previous studies (O'Connor et al. 2010; Gentet et al. 2012). Of the 18 OFF neurons, 14 were identified as EXC L5A cells, 3 as EXC L5B, and 1 as INH L4 cell. Thus, the majority of them (94.4%) were infragranular EXC neurons, which is in good agreement with a previous report (Oberlaender et al. 2012).

Next, we analyzed the latencies of cellular responses upon deflection of the principal or the neighboring whisker (Fig. 4A,B). For this computation only ON cells could be used, since the onset of the OFF responses could not be precisely determined from the limited number of trials recorded. Although a noticeable cell-to-cell variability was observed in the first-spike latencies (e.g., L5B EXC neurons presenting a total range of 18.1 ms), all recorded ON cells increased their activity within the first 30 ms after whisker stimulus (Supplementary Fig. 5). We observed significant differences in first-spike laten cies across the established neuronal groups (Fig. 4B1). L4 INH neurons showed significantly shorter latencies [median (interquartile range) = 7.8 (1.5) ms] when compared with all other groups (two-way ANOVA followed by LSD post hoc test, P < 0.05), indicating that they led the stimulus-evoked activity within the principal column (PC). In contrast, L2/3 EXC neurons presented the longest latencies within the PC [14.5 (3.8) ms]. A different sequence of activation was observed when the neighboring whisker was stimulated (Fig. 4B2). Under this condition, activity started in L4 and L5B INH neurons about 11 ms after stimulation, which presented significantly shorter latencies than all remaining EXC cell groups. Furthermore, the latency of L4 INH neurons, as well as L4 and infragranular EXC neurons, was significantly (P < 0.01) longer than after principal whisker stimulation. In summary, these data indicate a fast and nearly simultaneous intracolumnar processing of the sensory input led by L4 INH neurons ∼7.5 ms after whisker stimulus, which propagated rapidly to NCs ∼4 ms afterward.

To address the question, whether sensory-evoked activity may spread differently within a barrel arc or row, response latencies were analyzed considering those neural ensembles recorded from a barrel arc and those from a barrel row separately (data not shown). In both groups of experiments (arc and row), L4 INH neurons presented the shortest latencies (∼7.5 ms) of all neuronal groups when their principal whisker was stimulated. Further, when the neighboring whisker was stimulated the activity started in L4 and L5B INH neurons about 10–12 ms after stimulation in both groups of experiments, thus confirming the pooled data shown in Figure 4B.

Spike precision (defined as the SD of the first-spike in millisecond) was another parameter, which differed between the established neuronal groups. L4 INH neurons showed the most precise responses (2.2 ± 0.3 ms) upon principal whisker stimulation, being significantly (two-way ANOVA and LSD post hoc test, P < 0.05) more precise than all neuronal groups located in infragranular layers, and the L2/3 EXC neurons (Fig. 4C1). Further, L4 INH neurons showed also the highest precision (3.2 ± 0.4 ms) when the neighboring whisker was stimulated (Fig. 4C2); however, in this case, they were only significantly more precise than in L5A EXC neurons (4.9 ± 0.2 ms).

Next, we computed the number of evoked spikes per stimulus (spks/stim) in the range of 0–30 ms after whisker stimulation for the different neuronal groups, using all neurons located in the PC. Importantly, all INH neurons in granular and infragranular layers presented a significantly higher number of action potentials after principal whisker stimulation than the EXC neurons (Fig. 4D1). L4 INH neurons presented the highest number of spks/stim (0.67 ± 0.17) and L4 EXC neurons the lowest (0.15 ± 0.04). After neighboring whisker stimulation INH neurons in infragranular layers presented significantly higher activity than EXC ones (Fig. 4D2). In this case, L5B INH neurons presented the highest number of spks/stim (0.52 ± 0.09) and L4 EXC neurons the lowest one (0.15 ± 0.03). Thus, INH neurons presented not only higher spontaneous activity than EXC neurons in all layers, but also higher sensory-evoked activity in both the PC and the NC.

Similar to the response latencies, we also tested whether the number of spikes elicited per stimulus was different along the barrel arc or row (Supplementary Fig. 6). Whisker deflections elicited a higher number of spikes in L5B INH neurons of NCs located along the barrel row than along the arc. Accordingly, a significantly lower number of spikes were elicited in deep EXC neurons (L5B) in the “within-row” NCs than “across-row.” This result suggests that disynaptic lateral inhibition might be higher along the barrel row than along the arc within the deep layers. No other significant effect was found in any neuronal group when comparing responses along the arc or the row. Because all other features of the spread of activity were identical, and given that our purpose is to reveal basic principles of intra- and intercolumnar information transmission (independently of the precise barrels targeted), we performed our population-based analyses by pooling both sets of experiments together (see next section).

In summary, these results demonstrate that 1) stimulus-evoked responses are heterogeneous at the single-cell level, including neurons presenting ON responses, OFF responses, or no change in activity, 2) a fast processing of information occurs in the principal and NCs from ∼7.5 to ∼30 ms after whisker stimulus onset, being led by L4 and L5B INH neurons, 3) L4 INH neurons show the fastest and most precise responses, and 4) INH neurons reveal significantly less sparse responses than EXC ones in all layers, although they still respond sparsely to individual stimuli. In order to identify the information flow following whisker stimulation from these heterogeneous individual responses, we next analyzed whether we could detect temporally precise spike sequences, which might indicate distinct pathways of information propagation in a network (Abeles and Gerstein 1988; Nadasdy 2000; Izhikevich et al. 2004).

### Evoked Multineuronal Spike Sequences Were Variable and led by INH Cells

In order to characterize the flow of neuronal activity underlying sensory processing in the barrel cortex, a well-established template-matching algorithm was applied to search for stimulus-evoked repetitive spike sequences using a time precision of 2 ms (Abeles and Gerstein 1988; Sun et al. 2010) (see Materials and Methods). In brief, we considered only spike sequences that started within the first 30 ms poststimulus, containing at least 3 synchronous or consecutive spikes from at least 2 neurons and with a maximum delay of 10 ms between individual spikes. The statistical signiﬁcance of each sequence was evaluated by comparing its occurrence with that inside 2 surrogate datasets (generated by time and neuron jittering). A third surrogate method (trial shuffling) was used as a further control (Supplementary Fig. 4).

A total of 20 independent datasets from 10 animals were used for this analysis, each of them containing the responses to at least 60 (but typically 200) trials of individual whisker stimulation. Thus, 2 datasets were obtained from one animal if we recorded more than 60 trials upon stimulation of 2 different whiskers. Repetitive spike sequences appeared consistently in response to whisker stimulus and could be detected in all datasets (e.g., Fig. 5A,B). Typical spike sequences started in the PC at 7–10 ms after stimulus onset, propagating rapidly to adjacent ones within 10–20 ms after whisker stimulation (see below). In total, 16 174 different sequences were detected that started within the first 30 ms poststimulus.

Figure 5.

Occurrence of repetitive spike sequences after single-whisker defection. (A) Schematic illustration of the 8 × 16-channels probe covering the barrel-related columns C0 to C3. Barrels in the principal column (PC) and the neighboring columns (NCs) are marked with red and gray squares respectively. (B) Left, representative C1-whisker-evoked spike sequence #408 (shown in green) occurring in response to 3 individual stimuli. Black squares represent spikes from individual neurons occurring at 2-ms time bins after whisker stimulation. Right, same sequence presented on top of averaged activity over 100 trials. Color intensity represents the activity level in spikes per stimulus (see color bar at the right side). Note that the representative spike sequence occurs in 5 of the 100 trials. (C) Raster plot of 697 different spike sequences detected in 100 trials of C1 whisker stimulation. Spike sequences are ordered along the y-axis by their complexity (i.e., number of spikes constituting the sequence) and by the time of first appearance. Consecutive trials from 1 to 100 are represented in the x-axis. Upper subplots display for each trial (a) the percentage of new evoked spike sequences and (b) the number of evoked sequence occurrences. At the right side, 2 further subplots represent for each spike sequence (a) its complexity, that is, the number of spikes inside the sequence and (b) its response probability, that is, the total percentage of trials in which it appears. (D) Latency of spike sequence initiation after whisker stimulus. Central subplot represents for each detected spike sequence its response probability (coded color bars at the right side) along the time after whisker stimulation. Upper subplot depicts the total number of spike sequences initiated at a specific 2-ms time bin. Blue line represents sequences led by INH neurons, red by EXC ones, and black the summation of both. (E1) Number of detected spike sequences per dataset. Each bar represents the number of spike sequences starting at a specific time range after whisker stimulation. (E2) Contribution of INH (blue bars) and EXC (red bars) neurons to spike sequences starting at a specific time range after whisker stimulation. (E3) Contribution of spikes from INH and EXC neurons to spike sequences starting at a specific time range. (E4) Proportion of sequences led by INH and EXC neurons. Each bar represents the percentage of sequences in which the first neuron in the sequence is INH or EXC. Note that a particular sequence is considered to be led by both INH and EXC neurons if they fire simultaneously in the same 2-ms time bin used for spike sequence detection (see Materials and Methods). Therefore, the proportion of sequences led by INH neurons plus the proportion led by EXC neurons add up to more than 100%. Otherwise, same conventions as in previous panels.

Figure 5.

Occurrence of repetitive spike sequences after single-whisker defection. (A) Schematic illustration of the 8 × 16-channels probe covering the barrel-related columns C0 to C3. Barrels in the principal column (PC) and the neighboring columns (NCs) are marked with red and gray squares respectively. (B) Left, representative C1-whisker-evoked spike sequence #408 (shown in green) occurring in response to 3 individual stimuli. Black squares represent spikes from individual neurons occurring at 2-ms time bins after whisker stimulation. Right, same sequence presented on top of averaged activity over 100 trials. Color intensity represents the activity level in spikes per stimulus (see color bar at the right side). Note that the representative spike sequence occurs in 5 of the 100 trials. (C) Raster plot of 697 different spike sequences detected in 100 trials of C1 whisker stimulation. Spike sequences are ordered along the y-axis by their complexity (i.e., number of spikes constituting the sequence) and by the time of first appearance. Consecutive trials from 1 to 100 are represented in the x-axis. Upper subplots display for each trial (a) the percentage of new evoked spike sequences and (b) the number of evoked sequence occurrences. At the right side, 2 further subplots represent for each spike sequence (a) its complexity, that is, the number of spikes inside the sequence and (b) its response probability, that is, the total percentage of trials in which it appears. (D) Latency of spike sequence initiation after whisker stimulus. Central subplot represents for each detected spike sequence its response probability (coded color bars at the right side) along the time after whisker stimulation. Upper subplot depicts the total number of spike sequences initiated at a specific 2-ms time bin. Blue line represents sequences led by INH neurons, red by EXC ones, and black the summation of both. (E1) Number of detected spike sequences per dataset. Each bar represents the number of spike sequences starting at a specific time range after whisker stimulation. (E2) Contribution of INH (blue bars) and EXC (red bars) neurons to spike sequences starting at a specific time range after whisker stimulation. (E3) Contribution of spikes from INH and EXC neurons to spike sequences starting at a specific time range. (E4) Proportion of sequences led by INH and EXC neurons. Each bar represents the percentage of sequences in which the first neuron in the sequence is INH or EXC. Note that a particular sequence is considered to be led by both INH and EXC neurons if they fire simultaneously in the same 2-ms time bin used for spike sequence detection (see Materials and Methods). Therefore, the proportion of sequences led by INH neurons plus the proportion led by EXC neurons add up to more than 100%. Otherwise, same conventions as in previous panels.

In all experiments, spike sequences occurred in response to whisker stimulation during the whole recording (Fig. 5C). The total number of different evoked spike sequences grew dynamically across consecutive and identical stimuli (trials), indicating that the specific subset of active cells in the network, and also the temporal order of the active neurons inside the sequences, changed from one trial to the next. However, the growth rate of new sequences decreased over trials, meaning that early trials presented a higher percentage of new sequences when compared with latter trials (Supplementary Fig. 7). This suggests that the recorded units can display only a limited number of activity patterns and, therefore, their variability can be better estimated when a higher number of trials are considered.

In this regard, and just taking into account the subset of active cells and their temporal order, the total (theoretical) number of possible 3-spike sequences in the recording presented in Figure 5C (38 neurons) would be 383 − 38 = 54 834 sequences (note that at least 2 different neurons should be integrated in the sequence in order to be detected by our searching algorithms, see Materials and Methods). As a corollary, the number of detected 3-spike sequences after 100 trials (532) represents only ∼1% of the total theoretical number. Further, and taking into consideration not only the 3-spike but also the 4-spike sequences, the total (theoretical) number of sequences would be (384 – 38) + (383 – 38) = 2 139 932 sequences. In this case, the number of detected sequences after 100 trials (682) represents only 0.03% of the total theoretical number. Obviously, when taking into account sequences with higher number of spikes, these proportions are even lower. Given the fact that the number of observed sequences was far below the combinatorial possible number of sequences, it is probable that the number of spike sequences involved in the representation of sensory information is limited.

On average, 809 ± 187 different spike sequences were detected in each of the datasets in the first 30 ms after whisker stimulus (n = 20 datasets). The majority of sequences (71.6 ± 3.1%) showed a complexity of 3 spikes, 24.4 ± 2.3% 4 spikes, and a minority (4 ± 1%) more than 4 spikes. One intriguing feature of these repetitive spike sequences was their sparse occurrence, with a response probability (i.e., the total percentage of trials in which they appear) of 2.8 ± 0.3% (n = 20 datasets). Accordingly, the number of sequence occurrences per stimulus showed a high trial-to-trial variability (mean 20.8 ± 5.2, n = 20 datasets), with 18.7 ± 4.6% of the stimuli evoking no sequences and 3.6 ± 1.8% evoking more than 100 sequences—reflecting once more the dynamical nature of the spike sequence occurrence.

Despite the variability in the contents of sensory-evoked spike sequences, a clear regularity could be observed in their latency and the type of neuron leading the sequences. The majority of the detected spike sequences (95 ± 2%, n = 20 datasets) started between 0 and 20 ms after whisker stimulus (Fig. 5D–E1). Interestingly, the fraction of early spike sequences led by INH neurons was considerably higher than the fraction of INH neurons (or their spikes) found in repetitive spike sequences. The proportion of INH neurons inside sequences was similar for sequences initiated before 10 ms (21 ± 3%), between 10 and 20 ms (18 ± 2%), and between 20 and 30 ms (17 ± 3%) after whisker stimulation (Fig. 5E2) (two-way ANOVA followed by LSD post hoc tests), reflecting the total proportion of INH neurons in our recordings (14.6%). Spikes from INH neurons also represented the minority of spikes inside sequences initiated before 10 ms (37 ± 4%), between 10 and 20 ms (27 ± 3%), and between 20 and 30 ms (24 ± 6%) after whisker deflection (Fig. 5E3). However, 63 ± 7% of the sequences initiated before 10 ms after whisker stimulus were led by INH neurons (Fig. 5E4). This proportion was significantly higher than in sequences initiated between 10 and 20 ms (47 ± 5%, P < 0.05, n = 20 datasets) and between 20 and 30 ms (34 ± 8%, P < 0.001, n = 17 datasets) after whisker stimulus. In contrast, a relatively low proportion of spike sequences initiated within 10 ms after whisker stimulus were led by EXC neurons (42 ± 7%), being significantly higher in sequences initiated between 10 and 20 ms (66 ± 4%, P < 0.01, n = 20 datasets) and between 20 and 30 ms (77 ± 7%, P < 0.001, n = 17 datasets) after whisker stimulus.

We next asked the question, from which layer did these early and late sequences start. The spike sequences starting within the first 10 ms after whisker stimulation and led by INH neurons were mainly initiated in L4 (35.6 ± 9.3%, n = 19 datasets) and L5B (43.9 ± 8.9%, Supplementary Fig. 8A). In the sequences starting within 20–30 ms after whisker stimulation and led by EXC neurons, the highest proportion of sequences was initiated in L5A (Supplementary Fig. 8B, 57.2 ± 6.7%, n = 16 datasets). In conclusion, early sequences were preferentially led by INH neurons in L4 and L5B, and late sequences by EXC neurons in infragranular layers—possibly shedding light into the routing of information processing across specific cortical pathways, in this case emphasizing its modular organization into layers (see also Fig. 9 and Discussion).

Next, we addressed the question whether INH neurons may serve as hub neurons, defined as neurons having significantly higher numbers of functional links with other neurons (Sun et al. 2010). In total, 85 of the 443 cells included in this analysis (19.19%) were identified as hub neurons (n = 20 datasets). A relatively high proportion of these 85 hub neurons were INH (30 cells, 35.3%). Taking into account that the total proportion of INH neurons in our recordings was 14.6%, these data indicate that INH cells act preferentially as hub neurons than EXC ones. Importantly, the INH hub cells led the majority of sequences (58 ± 6%) within 0–10 ms after stimulus, meaning that only ∼5% of the sequences were led by INH nonhub neurons. Thus, a subgroup of INH cells (30 of a total of 64 INH cells, 46.9%) led almost all early sequences led by INH neurons, being located majorly (74%) in L4 and L5B (Supplementary Fig. 8A). We conclude that a specific group of INH hub cells located at thalamo-recipient layers mediates the dominance of feedforward inhibition at the early phase of the evoked responses.

To further verify the specificity of the spike sequences, we performed control analyses from all 20 datasets by changing the temporal precision of the template-matching algorithm from 2 to 4 ms. Although, in these analyses, a higher number of sequences were detected in the first 30 ms after whisker stimulus (1709 ± 383), early sequences (from 0 to 10 ms after stimulation) were again preferentially initiated by INH neurons (61.9 ± 5.7% led by INH vs. 47.4 ± 6 led by EXC), and late sequences (from 20 to 30 ms after stimulation) by EXC neurons (39.9 ± 5.8 led by INH vs. 79.9 ± 5 led by EXC). The total proportion of INH neurons inside sequences was similar for sequences initiated before 10 ms (19.4 ± 2.7%), between 10 and 20 ms (17.2 ± 1.7%), and between 20 and 30 ms (21.8 ± 4.6%) after whisker stimulation. This control ensured that our results were not caused by a particular high temporal precision used for the template-matching algorithm.

Nevertheless, it might be that early sequences were led by INH neurons just because they are transiently and rapidly activated by whisker stimulation when compared with EXC cells. To address this possibility, a further control analysis was performed by characterizing the properties of chance sequences detected within surrogate data (see Materials and Methods). A significantly lower number of sequences per dataset were detected in the first 30 ms after whisker stimulus within the surrogate (108 ± 76) than in the original (809 ± 187) data (Paired Wilcoxon signed-rank test, P < 0.0001, n = 20 datasets). While, in the original data, all datasets contained at least 20 sequences [total range (20–2872)], 13 datasets contained no sequences in their surrogates [total range (0–1494)]. These results suggest that the level of synchrony in the activity (higher in the original when compared with the shuffled data) correlates with the number of sequences found per dataset. However and as in the original datasets, the majority of the detected spike sequences (86 ± 5.1%, n = 20 datasets) also started between 0 and 20 ms after whisker stimulus in the surrogates (Supplementary Fig. 9A). We next compared the proportion of sequences led by INH cells in the original datasets and surrogates within the different time periods. Importantly, the proportion of sequences led by INH neurons was significantly smaller in the surrogates (53.7 ± 5.2%) when compared with the original datasets (78 ± 7.6%) in the early time period (0–10 ms after whisker stimulus) (Paired Wilcoxon signed-rank test, P < 0.05, n = 7 datasets) (Supplementary Fig. 9B). However, no significant differences were found either in the 10–20 ms period (47.9 ± 5% in the surrogate vs. 55.5 ± 7.1% in the original data, Paired Wilcoxon signed-rank test, P = 0.3, n = 7 datasets) or in the late (20–30 ms) period (31.4 ± 7.1% vs. 47.1%, P = 0.87). In summary, although we still found an overrepresented proportion of early chance sequences led by INH neurons in the surrogate data, it was significantly smaller than that found in the original data. Therefore, we conclude that the overrepresentation of INH cells as leading early sequences represents real mechanisms underlying multineuronal sequences, rather than epiphenomena arising from spike patterns of individual neurons.

As mentioned above, a mean of 20.8 ± 5.2 sequences occurred simultaneously per stimulus (n = 20 datasets), showing a high trial-to-trial variability (Fig. 5C). We next examined how these sequences were combined within the different trials. In order to approach this question, we first investigated the similarity of whisker-evoked spike sequences by detecting core sequences (Sun et al. 2010; Matsumoto et al. 2013). To this end, we compared all possible pairs of spike sequences from all trials in the same dataset to identify elementary (core) spike sequences (Sun et al. 2010). On average, 9.4 ± 1.5% of spike sequences were identified as core sequences (n = 20 datasets). A mean of 6.7 ± 2.3 core sequences occurred per trial, thus showing a response probability of 7 ± 0.8%—which was significantly higher than that of all spike sequences, 2.8 ± 0.3% (P < 0.05, Mann–Whitney U test, n = 20 datasets). This result suggests that, among the dynamic sensory-evoked spike sequences, core sequences were particularly activated with higher reliability.

Core sequences could be jointly combined, so that the resulting number of sequence occurrences was higher than the number of evoked spikes—see Figure 6A, where 15 evoked spikes generated a total of 19 combined core sequence occurrences. Single core sequences could occur isolated (blue sequence in Fig. 6B and green sequence in Fig. 6C), or in combination (Fig. 6D). In order to test whether similar sets of sequences tended to repeatedly appear in fixed combinations or different sets were recruited in different trials at a stochastic level, we analyzed all sequence combinations using the Ward's method (Fig. 6E, obtained from the same dataset depicted in Fig. 5C) (Matsumoto et al. 2013). The resulting dendrogram indicated that spike sequences were clustered into 9 major subgroups (see right side), thus suggesting that groups of sequences appeared rather in fixed combinations. However, a large subgroup of sequences (light green group at the middle part) did not show a homogeneous structure—that is, sequences inside this group did not appear jointly within the same trials, as they do within the other 8 groups. The dendrograms from the remaining 19 datasets were similar to this illustrative one, containing up to 15 major groups, from which one contained rather unclustered data. Compatible results have been observed in vitro (Matsumoto et al. 2013). We conclude from these results that, although in general the sequences appeared in fixed combinations, a proportion of them did so at a stochastic level.

Figure 6.

Combinations of sequence occurrences within individual trials. (A–D) Representative C1-whisker-evoked core sequences occurring in response to 4 individual stimuli. Black dots represent spikes from individual neurons. Colored lines represent core sequences. Barrel-related columns to which the neurons belong are depicted at the right side. For comparison, trial #36 (panel A) is also depicted in Figure 5B. (E) Raster plot of 697 different spike sequences detected in 100 trials of C1-whisker stimulation (same sequences and trials as in Fig. 5C). Spike sequences are reordered along the y-axis, so that sequences that appear most often jointly within the same trials are grouped together. Consecutive trials from 1 to 100 are represented in the x-axis. Right side, Ward's dendrogram depicting the clustering of the spike sequences. Distance between the leaf nodes (sequences) is represented horizontally (i.e., along the x-axis). Sequences which belong to the same major group (colors in the dendrogram) are considered to appear more often jointly within the same trials.

Figure 6.

Combinations of sequence occurrences within individual trials. (A–D) Representative C1-whisker-evoked core sequences occurring in response to 4 individual stimuli. Black dots represent spikes from individual neurons. Colored lines represent core sequences. Barrel-related columns to which the neurons belong are depicted at the right side. For comparison, trial #36 (panel A) is also depicted in Figure 5B. (E) Raster plot of 697 different spike sequences detected in 100 trials of C1-whisker stimulation (same sequences and trials as in Fig. 5C). Spike sequences are reordered along the y-axis, so that sequences that appear most often jointly within the same trials are grouped together. Consecutive trials from 1 to 100 are represented in the x-axis. Right side, Ward's dendrogram depicting the clustering of the spike sequences. Distance between the leaf nodes (sequences) is represented horizontally (i.e., along the x-axis). Sequences which belong to the same major group (colors in the dendrogram) are considered to appear more often jointly within the same trials.

In summary, our data demonstrate that although responses to individual stimuli are heterogeneous and sparse at the single-cell level, temporally precise spike sequences occur repeatedly after single-whisker deflection. Further, a consistent regularity could be observed in 1) the latency and 2) the type of neuron leading the responses to sensory stimulation. Thus, a subgroup of INH hub neurons located at thalamo-recipient layers led the majority of early evoked sequences, thus revealing the presence of feedforward inhibition. In addition, similar sets of sequences tended to repeatedly appear in rather fixed combinations within individual trials.

### Evoked Spike Sequences Generally Show Dispersed Spatial Structure

Next, we characterized for the above-described 16 174 repetitive spike sequences their spatial trajectories, defined by the location of the implicated neurons (within cortical layers and columns) during the progression of the sequence. Therefore, we classified the spike sequences into 5 major classes (Fig. 7A): 1) sequences which contained only neurons belonging to a single column, 2) sequences restricted to 2 NCs (further subdivided into sequences starting in the principal or NC), 3) sequences containing only neurons located in one layer, 4) sequences containing only neurons located at a specific layer and column (designated as clustered with a diameter of 300–400 µm), and 5) sequences which could not be included in any of the aforementioned 4 groups (designated as dispersed).

Figure 7.

Spatiotemporal organization of spike sequences. (A) Representative maps of 5 types of spike sequences with typical spatial trajectories. Data are from the same experiment presented in Figure 5A–D. Maps are represented onto the schematic illustration of the 8 × 16-channels probe covering the barrel-related columns C0 to C3. Barrels in principal column (PC) and neighboring columns (NCs) are marked with red and gray squares respectively. Layers are marked at the right side. Sequence of neuronal activation is indicated by the red arrows. (B) Proportion of sequences presenting specific spatial trajectories (n = 16 174 spike sequences from 20 datasets, 10 animals). Sequences restricted to a single column or cluster were further subdivided by those occurring at the PC or the NC. Sequences included in the PC/NC group were separated into those starting in the PC and those starting in the NC. Layer-restricted sequences were subdivided into specific layers. (C1–4) Spatiotemporal evolution of sequences confined to a single column (C1), single cluster (C3), single layer (C4), or restricted to the PC and a NC (C2). Sequences were subdivided in each of these panels by those starting before 10 ms, those between 10 and 20 ms, and those between 20 and 30 ms after whisker stimulus.

Figure 7.

Spatiotemporal organization of spike sequences. (A) Representative maps of 5 types of spike sequences with typical spatial trajectories. Data are from the same experiment presented in Figure 5A–D. Maps are represented onto the schematic illustration of the 8 × 16-channels probe covering the barrel-related columns C0 to C3. Barrels in principal column (PC) and neighboring columns (NCs) are marked with red and gray squares respectively. Layers are marked at the right side. Sequence of neuronal activation is indicated by the red arrows. (B) Proportion of sequences presenting specific spatial trajectories (n = 16 174 spike sequences from 20 datasets, 10 animals). Sequences restricted to a single column or cluster were further subdivided by those occurring at the PC or the NC. Sequences included in the PC/NC group were separated into those starting in the PC and those starting in the NC. Layer-restricted sequences were subdivided into specific layers. (C1–4) Spatiotemporal evolution of sequences confined to a single column (C1), single cluster (C3), single layer (C4), or restricted to the PC and a NC (C2). Sequences were subdivided in each of these panels by those starting before 10 ms, those between 10 and 20 ms, and those between 20 and 30 ms after whisker stimulus.

A minority of the spike sequences was located entirely inside a single column (6.8 ± 1.2%), a single layer (11 ± 1.2%) or a single cluster (2.2 ± 0.4%). A higher proportion was restricted to the PC and a neighboring one (19 ± 4.2%). The majority of the sequences presented a dispersed trajectory rather than a local one (61.1 ± 4.7%), indicating that the nature of the spike sequences is spatially dispersed (Fig. 7B). Sequences restricted to a single cluster appeared at almost identical proportions in the principal (52.1 ± 10.9%) and the NC (47.9 ± 10.9%). Of the sequences restricted to a single column, 38.3 ± 7.8% appeared in the principal and 61.7 ± 7.8% in a NC (Fig. 7B). In contrast, the majority (72.9 ± 5.4%) of sequences spanning NCs started in the PC, while only 27.1 ± 5.4% started in the NC. Of the sequences occurring inside a single layer, the majority remained in L5A (51.5 ± 6.9%) or L5B (39.7 ± 7.9%). Only a minority was confined to L4 (7.1 ± 4.8%) or L2/3 (1.7 ± 1.1%). In accordance with this wide spatial distribution of spike sequences, control experiments performed with one-shank electrodes (i.e., covering only one column) revealed that only a low proportion of stimuli (11.6 ± 7%, n = 8 datasets) succeeded to evoke spike sequences. Nevertheless, even in these datasets, 64.4% of early spike sequences were led by INH neurons, but only 13.3% in later sequences—thus corroborating previously outlined conclusions (see above).

Next, we analyzed the temporal evolution of sequences belonging to the different trajectory classes. The majority of all the sequences starting within the first 10 ms after whisker stimulus were initiated in the PC (63 ± 8%) and only 37 ± 8% (P < 0.05) in the NC. This result was robust for those sequences confined to a single column (60.6 ± 9.9%, n = 19), 2 NCs (71.6 ± 7.4%, n = 19), or a single cluster (73.2 ± 10.8%; n = 15; Fig. 7C1–3). In contrast, a significantly lower proportion of sequences starting between 20 and 30 ms after whisker stimulation was initiated in the PC (34 ± 9%, one-way ANOVA test, followed by LSD post hoc tests, P < 0.01). Accordingly, also a lower proportion of late sequences confined to a single column (24.8 ± 13.7%, n = 9) or a local cluster (33.3 ± 21.1%, n = 6) were initiated in the PC. Only late sequences restricted to 2 NCs were still mainly initiated in the PC (68 ± 12.9%, n = 7) (Fig. 7C2). While a substantial fraction of early sequences confined to single layers was initiated in L4 (11.1 ± 7%, n = 18) or L5B (47.9 ± 9.8%, n = 18), the vast majority of sequences starting in the interval 20–30 ms after stimulations was initiated in L5A (76.4 ± 6.2%, n = 13, Fig. 7C4). Finally, the proportion of dispersed sequences was significantly lower in those initiated within the first 10 ms after stimulus (57.6 ± 5.6%, n = 20) than in the interval of 20–30 ms after whisker stimulation (67.7 ± 6.6%, n = 17) (P < 0.05).

In order to investigate the role of INH neurons in different spatial trajectory classes, we analyzed the proportion of spike sequences led by INH neurons in the 5 different trajectory classes. These analyses revealed no significant difference (ANOVA test, P > 0.05) in the proportion of spike sequences led by INH neurons between the different classes. For all trajectory classes, the proportion of spike sequences led by INH neurons was highest for sequences starting within 10 ms after stimulus (Fig. 8). We conclude from this result that the type of neuron leading the spike sequences does not determine their spatial trajectory (and vice versa).

Figure 8.

Inhibitory (INH) interneurons lead sequences independently of their spatial structure. Each symbol represents the proportion of spike sequences led by INH neurons for a different trajectory class. Sequences were subdivided by those starting before 10 ms, those between 10 and 20 ms, and those between 20 and 30 ms after whisker stimulus.

Figure 8.

Inhibitory (INH) interneurons lead sequences independently of their spatial structure. Each symbol represents the proportion of spike sequences led by INH neurons for a different trajectory class. Sequences were subdivided by those starting before 10 ms, those between 10 and 20 ms, and those between 20 and 30 ms after whisker stimulus.

In summary, these results demonstrate that although sensory-evoked temporally precise spike sequences are dispersed in their spatial structure, spike sequences initiated immediately after the stimulus are led by L4 and L5B INH neurons in the PC (Fig. 9A,C), whereas spike sequences that started at later poststimulus intervals are preferentially led by L5A EXC neurons in the NC (Fig. 9B,C).

Figure 9.

Subdivision of early and late sequences by the layer, column, and neuronal type of the leading neuron. (A) Spatial distribution of neurons leading temporally precise spike sequences initiated within the first 10 ms after whisker stimulation. Sequences were subdivided by 3 attributes of the leading neuron: its cortical layer, its location in the principal column (PC) or the neighboring column (NC), and its type (INH/EXC). Bars represent mean ± SEM. (B) Same for sequences initiated within 20–30 ms after whisker stimulus. (C) Schematic illustration depicting the major findings. For the early sequences (left side, panel A), a significant higher proportion of leading neurons were INH neurons located in L4 and L5B of the PC. In contrast, for late sequences (right side, panel B), a significantly dominant fraction was led by EXC neurons in L5A of the NC. Note that, in all cases, the majority of the sequences were spatially dispersed (see Fig. 7).

Figure 9.

Subdivision of early and late sequences by the layer, column, and neuronal type of the leading neuron. (A) Spatial distribution of neurons leading temporally precise spike sequences initiated within the first 10 ms after whisker stimulation. Sequences were subdivided by 3 attributes of the leading neuron: its cortical layer, its location in the principal column (PC) or the neighboring column (NC), and its type (INH/EXC). Bars represent mean ± SEM. (B) Same for sequences initiated within 20–30 ms after whisker stimulus. (C) Schematic illustration depicting the major findings. For the early sequences (left side, panel A), a significant higher proportion of leading neurons were INH neurons located in L4 and L5B of the PC. In contrast, for late sequences (right side, panel B), a significantly dominant fraction was led by EXC neurons in L5A of the NC. Note that, in all cases, the majority of the sequences were spatially dispersed (see Fig. 7).

### Role of Specific Layer 5 EXC Cell Types in Synchronizing Intracolumnar Activity

Despite the prominently dispersed nature of the temporally precise spike sequences described above, we further investigated whether specific L5 cell types could play a role in columnar synchronization—a hypothesis which is supported by their dendritic and axonal morphologies (Oberlaender et al. 2011, 2012).Therefore, we performed a classification of putative cell types based on their cortical depth and spontaneous activity. Neurons in infragranular layers showed higher spontaneous FRs than in supragranular and granular layers, reaching a maximum at a somatic depth of ∼1125 µm under the cortical surface, which coincides with the location of L5B previously reported (Fig. 10A) (Meyer et al. 2010; Oberlaender et al. 2012). We classified deep L5 neurons (1025–1200 µm) with relatively high spontaneous FRs (≥2 spks/s) as putative L5 thick-tufted (L5tt) neurons (n = 27) and more superficial (800–950 µm) neurons with lower spontaneous FRs (0.1–1.1 spks/s) as presumed L5 slender-tufted (L5st) neurons (n = 28).

Figure 10.

Properties of specific L5 cell types. (A) Spontaneous activity versus recording depth for the 263 EXC neurons stimulated by principal whisker deflection. Putative layers are marked at the upper part. Squares and triangles mark the selected populations of L5 thick-tufted (L5tt) and L5 slender-tufted (L5st) cells, respectively. (B) Comparison of mean firing rate (FR) of EXC neurons before versus after whisker stimulation, plotted in logarithmic scale. Pre- and poststimulus times are defined as in Figure 3B. Empty symbols represent neurons showing no significantly different FR in both periods, filled symbols neurons showing significantly different FR (P < 0.05). Rhomboids mark the selected L5tt cells shown in Figures 3A and 4A. Larger circle marks the selected L5st cell shown in Figure 4A. (C) Ensemble activity of L5tt and L5st cells, both evoked after principal (upper part) or neighboring whisker (lower part) stimulation. Continuous lines represent the mean spike density functions (SDFs) of both populations of cells (see Materials and Methods). Dashed lines represent 95% confidence intervals. Numbers of cells included in each of the ensembles are noted in the legend at the upper-right corner.

Figure 10.

Properties of specific L5 cell types. (A) Spontaneous activity versus recording depth for the 263 EXC neurons stimulated by principal whisker deflection. Putative layers are marked at the upper part. Squares and triangles mark the selected populations of L5 thick-tufted (L5tt) and L5 slender-tufted (L5st) cells, respectively. (B) Comparison of mean firing rate (FR) of EXC neurons before versus after whisker stimulation, plotted in logarithmic scale. Pre- and poststimulus times are defined as in Figure 3B. Empty symbols represent neurons showing no significantly different FR in both periods, filled symbols neurons showing significantly different FR (P < 0.05). Rhomboids mark the selected L5tt cells shown in Figures 3A and 4A. Larger circle marks the selected L5st cell shown in Figure 4A. (C) Ensemble activity of L5tt and L5st cells, both evoked after principal (upper part) or neighboring whisker (lower part) stimulation. Continuous lines represent the mean spike density functions (SDFs) of both populations of cells (see Materials and Methods). Dashed lines represent 95% confidence intervals. Numbers of cells included in each of the ensembles are noted in the legend at the upper-right corner.

Both classes of L5 neurons presented heterogeneous responses to principal whisker stimulation, thus showing ON and OFF responses, or a lack of them (NR cells) (Fig. 10B, compare to Fig. 3B) (distributions not significantly different, χ2 test, P = 0.86). L5tt neurons showed, however, a higher number of evoked spikes within 30 ms after principal whisker stimulation than L5st cells (0.44 ± 0.1 vs. 0.2 ± 0.06 spks/stim, Mann–Whitney U test, P < 0.05). This difference between both cell groups was even higher when the neighboring whisker was stimulated (0.5 ± 0.1 vs. 0.13 ± 0.03 spks/stim Mann–Whitney U test, P < 0.01). Despite these differences in the activity level of the evoked responses, neural latencies were not significantly different between both groups of L5 cells. Thus, after principal whisker stimulation L5tt cells showed similar first-spike latencies [median (interquartile range) = 11.5 (3.5) ms] as L5st cells [11.8 (2.4) ms] (Mann–Whitney U test, P = 0.95). After neighboring whisker stimulation both neuronal groups showed longer latencies, but not significantly different from each other [16.3 (9.8) vs. 12.7 (3.3) ms, Mann–Whitney U test, P = 0.16).

These functional features were noticeable when the ensemble activity was plotted (Fig. 10C). Both after principal and neighboring whisker stimulation, the activity peak at ∼12 ms after stimulus was higher in L5tt than in L5st cells. After neighboring whisker stimulation, the activity level within 20 and 40 ms after stimulus was ∼15 spks/s in L5tt cells, indicating a higher level of sustained activity in the neighboring than in the PC. These data suggest that, while both neuronal groups contributed significantly to the activity at the PC, only L5tt neurons extended their contribution to the NCs.

In order to further explore this issue, we analyzed the distribution and structure of the spike sequences originating in both L5 neuronal groups. For both L5tt and L5st cells, the majority of the sequences initiated within 10–20 ms after whisker stimulation originated in the PC (5.8 ± 2.7% and 2.8 ± 1%, respectively) rather than in the first- or second-order NCs (4.1 ± 1.7% and 0.97 ± 0.4%) (n = 16 174 sequences from 20 datasets, 10 animals) (Fig. 11A). Importantly, this distribution switched in the late period for L5tt cells only, during which more sequences originated in the NCs (5.4 ± 2.1%) than in the principal one (3.1 ± 1.4%).

Figure 11.

Differential involvement of specific L5 cell types in spike sequences and their possible role in intracolumnar synchronization. (A) Percentage of sequences led by L5 thick-tufted (L5tt) and L5 slender-tufted (L5st) in the principal column (PC) or the first-/second-order neighboring columns (NCs). Sequences are subdivided by those starting before 10 ms, those between 10 and 20 ms, and those between 20 and 30 ms after whisker stimulus. Percentages are computed by averaging the data from all datasets (20 datasets from 10 animals). Note that L5tt cells show a temporal evolution of the column from which the majority of the sequences are initiated (not the case for L5st cells). (B) Upper part, for each group of sequences led by L5tt neurons in the NCs and starting at a specific time period, percentage presenting a specific spatial trajectory. For this computation, all sequences led by L5tt cells in the NCs starting at the different time periods were pooled together (from early to late sequences n = 93; 743 and 169 sequences). Lower part, same for sequences led by L5st neurons in the PC (n = 76; 448 and 107 sequences). (C) For each group of sequences led by L5tt/L5st neurons in the PC/NCs and starting at a specific time period, percentage ending in neurons located in L2/3. (D) Schematic illustration depicting the proposed role of L5st and L5tt neurons in intracolumnar synchronization. For the sequences within 10–20 ms after stimulus (left side), a high proportion of the sequences led by L5st and L5tt neurons started at the PC, being mostly spatially dispersed. Later sequences (right side) led by L5st neurons showed a higher proportion of sequences restricted to the PC. In contrast, the majority of the sequences led by L5tt neurons started in the NCs, showing as well an overrepresentation of intracolumnar multineuronal sequences.

Figure 11.

Differential involvement of specific L5 cell types in spike sequences and their possible role in intracolumnar synchronization. (A) Percentage of sequences led by L5 thick-tufted (L5tt) and L5 slender-tufted (L5st) in the principal column (PC) or the first-/second-order neighboring columns (NCs). Sequences are subdivided by those starting before 10 ms, those between 10 and 20 ms, and those between 20 and 30 ms after whisker stimulus. Percentages are computed by averaging the data from all datasets (20 datasets from 10 animals). Note that L5tt cells show a temporal evolution of the column from which the majority of the sequences are initiated (not the case for L5st cells). (B) Upper part, for each group of sequences led by L5tt neurons in the NCs and starting at a specific time period, percentage presenting a specific spatial trajectory. For this computation, all sequences led by L5tt cells in the NCs starting at the different time periods were pooled together (from early to late sequences n = 93; 743 and 169 sequences). Lower part, same for sequences led by L5st neurons in the PC (n = 76; 448 and 107 sequences). (C) For each group of sequences led by L5tt/L5st neurons in the PC/NCs and starting at a specific time period, percentage ending in neurons located in L2/3. (D) Schematic illustration depicting the proposed role of L5st and L5tt neurons in intracolumnar synchronization. For the sequences within 10–20 ms after stimulus (left side), a high proportion of the sequences led by L5st and L5tt neurons started at the PC, being mostly spatially dispersed. Later sequences (right side) led by L5st neurons showed a higher proportion of sequences restricted to the PC. In contrast, the majority of the sequences led by L5tt neurons started in the NCs, showing as well an overrepresentation of intracolumnar multineuronal sequences.

Although the spatial structure of the evoked sequences was mostly dispersed (as shown in the previous section), we observed a temporal evolution in the proportion of sequences restricted to a single column (Fig. 11B). Thus, those sequences led by L5tt neurons in the NCs were more often restricted to a single column in the late (20–30 ms after stimulus) (34.9%) than in the early (0–10 ms) time period (11.8%, more similar to the generally observed ∼7%, see Fig. 7) (n = 169 vs. 93 sequences, respectively). The same effect was revealed in those sequences led by L5st neurons in the PC (24.3% vs. 7.9%, n = 107 vs. 76 sequences) (Fig. 11B). Also, the remaining sequences (those led by L5tt neurons at the PC and those led by L5st neurons at the NCs) showed a similar increase in the proportion of intracolumnar sequences (data not shown).

Concordantly, the proportion of sequences ending in L2/3 neurons increased for late sequences compared with the previous ones (30 ± 2% in the period 20–30 ms after stimulus vs. 12.2 ± 2.5 within 0–10 ms) (n = 384 vs. 559 sequences, respectively) (Fig. 11C). This effect was stable across both L5 cell populations in the PC (L5tt and L5st), and across the principal and the NCs for the L5tt cells. These data suggest that L5 neurons play a major role in synchronizing the intra-columnar network activity during the late phases of the evoked responses (Fig. 11D). While L5st neurons might play this important role in the PC, L5tt could be responsible for extending the intracolumnar synchronization towards the first- and second-order NCs, thereby mediating the propagation of evoked activity to more distant columns.

## Discussion

In the present study, we investigated at a large scale the flow of sensory-evoked activity in the adult rat barrel cortex in vivo using multielectrode arrays. We analyzed the spatiotemporal properties of functionally and topographically identified neurons and their integration into temporally precise spike sequences, which have been proposed as a signature for the flow of activity within a neural network (Abeles and Gerstein 1988; Nadasdy 2000; Izhikevich et al. 2004). Therefore, we characterized the distribution of these sequences within and between functional cortical columns. Several major conclusions can be drawn from our experiments: 1) temporally precise repetitive spike sequences are elicited by single-whisker stimulation, however with a low response probability; 2) repetitive spike sequences initiated immediately after the sensory stimulus are preferentially led by a subgroup of INH neurons in L4 and L5B of the PC; 3) later spike sequences are preferentially led by L5A EXC neurons in the NC; 4) similar sets of sequences tend to appear repeatedly mainly in stable combinations within individual trials, rather than in a stochastic manner; 5) although the starting point of the sequences is confined to individual cortical modules, their further spatiotemporal trajectories are distributed across the entire recorded area, indicating that the intra- and intercolumnar processing of information are performed in parallel; and 6) putative L5st and L5tt neurons play a major role in synchronizing the intra-columnar network activity during the late phases of the evoked responses. In summary, our data constitute the first large-scale characterization of neuronal activity propagation at high temporal resolution in the neocortex in vivo, describing the activity flow across different cortical layers, columns, and neuronal types.

### Spatial Resolution of the Multielectrode Recordings

In the present study, we used the early CSD sinks arising at the thalamo-recipient L4 and L5B after sensory stimulation in order to assign the individual channels to specific cortical layers (Supplementary Fig. 1) (Mitzdorf 1985; Csicsvari et al. 2003; Buzsaki et al. 2012). An obvious limitation of our procedure is the spatial separation between neighboring channels within a shank (75 µm in the 8 × 16 probe; 100 µm in the 1 × 16 probe). Moreover, the assignment of layer borders is nontrivial due to the complex 3D structure of the barrel cortex—particularly in the light of recent findings showing that cortical thickness and even the thicknesses and relative proportions of layers deviate between columns in different rows (Egger et al. 2012; Meyer et al. 2013). Thus, layer borders deviate between columns, but given the coarse spatial resolution of the present data, the layer borders determined here are within the wide range of previous reports (Meyer et al. 2010; Oberlaender et al. 2012; Feldmeyer 2012) and, more importantly, changing the layer borders within the given resolution (i.e., ∼100 µm) did not result in significant changes in the latencies or the level of the sensory-evoked responses computed across layers (data not shown).

The spontaneous activity of L5B EXC neurons in our dataset is smaller than that previously reported (de Kock et al. 2007; Oberlaender et al. 2012). Two possible explanations could account for this incongruence. First, it might well be that the depth referred to as L5B is more likely the area of L5 where L5st (L5A) and L5tt (L5B) neurons intermingle (Oberlaender et al. 2011). Alternatively, the area referred to as L5B could also include neurons actually located in layer 6A (L6A), which present lower spontaneous FRs (de Kock et al. 2007). In any case, the distribution of FRs across cortical depths was very similar to that reported by means of juxtacellular recordings [compare our Fig. 10A with Fig. 1D in de Kock and Sakmann (2009)], therefore supporting the validity of our procedure for layer separation.

The spatial resolution of the present method is not only determined by the separation of the probes, but also by the spatial spread of the LFP signal (∼200 µm) (Kajikawa and Schroeder 2011; Linden et al. 2011). Our current method might be thus unsuitable to pick up potentially important functional (Petersen et al. 2003) and structural (Meyer et al. 2013) differences across rows and arcs. In this regard, the reported preference of neural excitation spreading along the rows when compared with arcs in supragranular layers—shown by means of VSD imaging data (Petersen et al. 2003)—could not be observed in our evoked patterns of CSDs or LFP signals (data not shown). Furthermore, our recordings did not reach the spatial resolution necessary to reliably characterize septal-related areas (Supplementary Fig. 2). Future developments in silicon probe engineering could make possible to further address these important issues at high spatiotemporal resolution.

### Sparse and Heterogeneous Responses at the Single-Cell Level

A salient property of the neural responses was their sparse firing in response to sensory stimulation, an observation which is in good agreement with previous studies (Simons and Carvell 1989; de Kock et al. 2007; Kerr et al. 2007). However, INH cells showed a mean of up to ∼0.7 spks/stim in L4, supporting previous studies reporting that more reliable and precise action potentials can be elicited in fast-spiking (FS) interneurons (Bruno and Simons 2002; Sun et al. 2006). Moreover, the higher spike probability in INH neurons was maintained from the principal to the NC, being in this case higher in L5B than in L4. This is in agreement with previous reports indicating a stronger spread of activity in infragranular than in granular and supragranular layers (Sakata and Harris 2009; Wester and Contreras 2012). Our observation that INH neurons responded less sparsely than EXC neurons in the principal, but also in NCs, is compatible with their hypothesized role in shaping network activity within the cortex (Isaacson and Scanziani 2011). In addition, L4 INH neurons showed consistently the shortest spike latency, which is in good agreement with a recent study (Kimura et al. 2010). This result can be attributed to the specific innervation from the ventral posterior medial (VPM) thalamic nucleus to L4 (Swadlow 2003), intrinsic synaptic and membrane properties of INH neurons (Doischer et al. 2008), and a higher conduction velocity of thalamo-cortical axons innervating INH neurons (Kimura et al. 2010).

While the relative fraction of INH cells across all layers is around 14%, thus matching previous reports (Sirota et al. 2008; Sakata and Harris 2009), the relative proportions within each layer deviate from previous observations. Meyer et al. (2011) report the following proportions of GABAergic cells: L2/3, 17/9%; L4, 8%; L5A, 20%; L5B, 16%. In our study, the proportions are different (see Results): L2/3, 11%; L4, 25%; L5A, 8%; L5B, 19%. Of particular interest is the bias to record from more INH cells in L4 (25% vs. 6%). One explanation for this bias would be that VPM neurons preferentially innervate and/or activate INH cells in L4 after passive whisker deflection (an argument in line with the less sparse responses found in INH cells, see above) (Bruno and Simons 2002; Swadlow 2003). Alternatively, it might be possible that the current method fails to detect and/or discriminate EXC cells. However, we consider this possibility unlikely, since such a shortcoming should affect all recorded cortical layers, which was however not the case in supra- and infragranular layers.

Another property of the recorded neurons was the heterogeneity in their responses. The majority of the EXC neurons showed ON or no significant responses to principal whisker stimulation, both response types which have been previously reported (Simons and Carvell 1989; Armstrong-James et al. 1992; Petersen and Diamond 2000). Interestingly, a small but considerable proportion of cells showed significantly lower FRs after whisker stimulus (OFF cells), an effect which has been already observed (see Fig. 8A in Oberlaender et al. 2012). Under urethane anesthesia, there is almost a lack of such OFF responses to passive whisker touch across all layers and cell types of the barrel cortex (Simons and Carvell 1989; de Kock et al. 2007). Interestingly, these rare OFF responses are not present during free-whisking behavior (de Kock and Sakmann 2009; Oberlaender et al. 2012). However, OFF responses have been found by means of cell-attached recordings in mice performing vibrissa-based object localization tasks (O'Connor et al. 2010), and also in identified groups of interneurons during quite wakefulness in response to whisker stimulation (Gentet et al. 2012). This would be in accordance with the strong reduction in the amplitude of the passively evoked responses in L2/3 during active whisking recorded by means of VSD imaging (Ferezou et al. 2006). In conclusion, and in agreement with recent reports obtained in the visual cortex (Haider et al. 2013) and auditory cortex (Sakata and Harris 2012), it is very likely that the presence of inhibition (and thus OFF responses) in granular and supragranular layers of the barrel cortex is enhanced during wakefulness and lighter stages of anesthesia. Further, responses of individual neurons in the rodent barrel cortex are likely as heterogeneous as in the visual cortex (Hubel and Wiesel 1959). In line with this, individual cells have been shown to present different responses depending on stimulus direction (Simons and Carvell 1989) or multiwhisker interactions (Le Cam et al. 2011).

### Flow of Multineuronal Activity in the Barrel Cortex

In the present study, temporally precise spike sequences were identified as a manifestation of the activity flow across cortical layers and columns. Recurring sequences of neural activity, which have been proposed as a potential substrate of both information transfer and transformation in cortical networks, have been observed in the intact brain (Nadasdy et al. 1999), in brain slices (Ikegaya et al. 2004), in dissociated cortical cultures (Rolston et al. 2007; Sun et al. 2010), and in silico (Izhikevich et al. 2004). This is in agreement with the concept that a certain degree of synchrony must be achieved by the activity in order to propagate through the cortical network (Aertsen et al. 1996; Diesmann et al. 1999; Bruno and Sakmann 2006; Kumar et al. 2010; Stanley 2013). In contrast to optical methods, spike sequences identified with multielectrode arrays will represent only a small part of the entire stimulus-evoked synfire chain (Aertsen et al. 1996), due to the limited number of recorded neurons (Rolston et al. 2007; Sun et al. 2010). On the other hand, electrophysiological recordings offer the advantage of a higher temporal resolution and simultaneous recordings from supragranular, granular, and infragranular layers in vivo.

The issue of multineuronal sequences has been controversial in the literature, with several studies asserting that electrophysiological evidence for such sequences can be a stochastic artifact (Mokeichev et al. 2007). In the present study, the statistical signiﬁcance of each detected sequence was evaluated by comparing its occurrence with that inside 3 surrogate sets generated by spike shufﬂing algorithms (see Materials and Methods and Supplementary Fig. 3). These algorithms (time and neuron jittering and trial shuffling) conserved specific features of the multineuronal activity while destroying others. For instance, time jittering destroyed the precise temporal resolution of the activity, keeping the mean spike timing (i.e., coarse PSTH) and trial order unchanged; in contrast, trial shuffling destroyed pairwise (noise) correlations among neurons, while maintaining intact the single-cell PSTHs and their precise temporal resolution (for review see Grun 2009). The fact that a specific sequence appeared statistically more often in the experimental than in the surrogate data indicated its validity—thus deriving from true mechanisms underlying multineuronal sequences, rather from epiphenomena (or side effects) that naturally arise from spike patterns of individual neurons.

One of the most important properties of the detected spike sequences was their dynamic nature, both in their composition and spatial extent. Individual spike sequences did not repeat very consistently, indicating that the identity of the neurons constituting the sequences and their order changed from trial to trial. This property, which has been already reported for L2/3 neurons by means of calcium imaging (Kerr et al. 2007), is however expected from the low number of spikes elicited per stimulus by single neurons and by their low spike precision. The dynamics of ongoing activity can also contribute to the variability in sensory-evoked neuronal responses (Arieli et al. 1996). We conclude from this observation that the defined intra- and intercolumnar flow of neuronal activity does not rely on the activity of specific neuronal connections. In contrast, we assume that neuronal activity is distributed via a dispersed network of individual connections, with interneurons shaping the immediate response properties (see below). However, it is important to note that the number of observed sequences was far below the total number of possible combinations, suggesting that the number of spike sequences involved in the representation of sensory information within the subset of cortical neurons is rather limited. As a corollary, the population of recorded cells in the barrel cortex most likely constitutes a circuit with a subset of possible dynamic states, despite the trial-to-trial variability found in their responses.

While the majority of sequences showed a dispersed distribution (thus spanning at least 2 functional columns), only ∼7% of the spike sequences were confined to a single column. In this respect, it must be considered that these spike sequences can last for up to 50–70 ms poststimulus. The dispersed nature of the spike sequences may thus reflect the distributed nature of cortical information processing during later stages of sensory processing, which could be mediated by neurons spanning different cortical columns (Petersen and Diamond 2000; Petersen et al. 2003). Intracortical horizontal connections can mediate these large-scale interactions (Luhmann et al. 1990a, 1990b), a concept which has been proven by means of dynamic photo stimulation in vitro (Nawrot et al. 2009; Boucsein et al. 2011). However, if only spike sequences that start immediately after the stimulus were considered, the majority of these sequences were initiated within the PC. Thus, our results are in line with the hypothesis that cortical columns are a fundamental functional element present in a diversity of cortical areas (Mountcastle 1997) including the barrel cortex (Feldmeyer et al. 2013), which implies that the initial cortical responses were mainly restricted to L4 and L5B of the whisker-related column (Petersen and Diamond 2000; Petersen et al. 2003). In accordance with the postsynaptic targets of lemniscal thalamo-cortical fibers, we found that the majority of early spike sequences were indeed initiated in L4 and L5B (Petersen and Diamond 2000; Petersen et al. 2003; Meyer et al. 2010; Oberlaender et al. 2012).

The columnar organization is further supported by the finding that putative L5 cell types led a relatively high proportion of late sequences restricted to just one column. The paralemniscar pathway could play a major role in the PC, where L5st cells showed an overrepresentation of intracolumnar activity; at the same time, the lemniscar pathway could be responsible for propagating the intracolumnar synchronization horizontally—thus attributing different but complementary roles to both pathways in columnar organization (Chagnac-Amitai et al. 1990; Ahissar et al. 2000; Feldmeyer et al. 2013). In addition, these assumptions are very much in line with previous functional and anatomical data on L5st and L5tt cells. L5tt cells have broad passive deflection receptive fields, that is, high surround whisker-evoked spiking (de Kock et al. 2007), while L5st cells have narrow receptive fields. Furthermore, L5tt axons innervate surrounding columns only within L5, but remain confined to the PC in L2–4. The opposite is true for L5st axons that remain confined to the PC in L5 but densely innervate L2/3 in the PC and surround columns (Oberlaender et al. 2011). Thus, it has been previously concluded that strong coupling between L5tt cells may result in synchronous activation of L2/3 neurons in multiple columns after principal whisker deflection (Varga et al. 2011), while L5st neurons may serve to lock L2/3 dendrites across columns to the phase of the whisking cycle during periods of free whisking in awake animals (Oberlaender et al. 2011).

L4 and L5B INH neurons led the majority of spike sequences that were initiated immediately after whisker stimulation, despite the fact that they neither represent the majority of neurons nor fire the majority of spikes involved in spike sequences—thereby revealing the presence of feedforward inhibition. The higher spike time precision of INH neurons, in particular in L4, may contribute to their dominance as leading neurons in spike sequences. This dominance as leading neuron, as well as their spike precision (Doischer et al. 2008), suggests a possible role of INH interneurons in determining information flow and processing in the brain (Bonifazi et al. 2009; Isaacson and Scanziani 2011). Two possible mechanisms (not mutually exclusive) could account for the possible role of INH neurons in the formation of temporally correlated assemblies: 1) L4 and L5B INH neurons could narrow the subsequent time window for neural excitation (also called window of opportunity) (Kimura et al. 2010) or 2) the firing of L4 SOM INH interneurons could inhibit local FS interneurons, subsequently increasing the firing probability of EXC cells (Xu et al. 2013).

We consider the presence of an additional EXC cortical neuron driving these INH neurons implausible, since they have first-spike latencies of ∼7.5 ms. Comparable short latencies have been also found in previous studies using intra- or juxtacellular recordings upon whisker stimulation (Wilent and Contreras 2004; de Kock et al. 2007), preferentially in INH cells (Kimura et al. 2010). Under consideration of the trisynaptic whisker to barrel cortex pathway (Feldmeyer et al. 2013), this short latency does not justify the integration of an additional presynaptic neuron.

On the other hand, the observation that the vast majority of spike sequences during later poststimulus intervals were led by infragranular EXC neurons is in line with the dominance of spikes from L5tt neurons during this interval (de Kock et al. 2007)—which send cortico-cortical projections not only monosynaptically, but also through higher order feedforward pathways via the thalamus (Guillery and Sherman 2002). Thus, early sequences led by L4/5B INH neurons are relayed to the cortex, while later sequences starting in L5 EXC neurons most likely propagate not only intracortically, but also to subcortical structures—thereby exhibiting the routing of information processing across specific thalamo-cortical pathways.

### Cortical Multineuronal Spike Sequences in Different Brain States

As discussed above, it is likely that the inhibition in granular and supragranular layers is stronger during wakefulness and lighter stages of anesthesia (Ferezou et al. 2006; Sakata and Harris 2012; Haider et al. 2013), which ultimately produces lower FRs (i.e., sparser activity) in pyramidal cells. This might be a mechanism by which neurons in superficial layers enhance the extraction of behaviorally relevant signals from noisy brain activity (Sakata and Harris 2012). Importantly, in the same study (Sakata and Harris 2012), both INH and EXC neurons in deeper layers showed heterogeneous effects in their sensory-evoked activity during lighter states of anesthesia, thus requiring a different (and probably more complex) interpretation.

The majority of studies characterizing precise spatiotemporal patterns of activity from larger (>25) populations of neurons have been done either in vitro (Mao et al. 2001; Ikegaya et al. 2004; MacLean et al. 2005) or in silico (Abeles and Gerstein 1988; Aertsen et al. 1996; Diesmann et al. 1999; Izhikevich et al. 2004; Schrader et al. 2008). Those studies relating the synchronized activity to the brain state of the animal in vivo have covered smaller populations of cells (≤15); for instance, a key study has shown that patterns of coincident action potentials occur in the monkey motor cortex in relation to specific stimuli or movements during a behavioral task (Riehle et al. 1997). This has been extended to specific time delays between pairs of neurons (Shmiel et al. 2006), which are also present in close neighboring neurons in the premotor cotex (Sakurai and Takahashi 2006). Concordant results have been obtained in the auditory cortex of freely moving rats by showing that particularly precise multineuronal spike sequences predict future behavioral responses in a go/no-go task (Villa et al. 1999). The probability of detecting repetitive patterns increase however with the number of neurons recorded in parallel (Tetko and Villa 1997). Therefore, the higher occurrence of precise spatiotemporal activity in the cortex of awake, behaving animals might have been obscured by the fact that a lower number of cells have been recorded under these conditions.

Consistent with this argument, the synchronized activity of large neuronal groups in the gamma-frequency band (30–100 Hz) is enhanced by attentional processes (Fries et al. 2001), thereby shortening the behavioral response times to attended stimuli (Womelsdorf et al. 2006). The properties of the multineuronal firing rhythm are determined by the fine timescales of EXC and INH currents, and by the balance between the mean recurrent excitation and inhibition (Brunel and Wang 2003). Therefore, the presence of synchronized gamma activity reveals a temporally precise coordination in the population activity—which is in agreement with the observed increase in the reliability of firing sequences during beta/gamma oscillations in the cat visual cortex (Havenith et al. 2011). Future studies should address systematically to which extent the trial-to-trial variability of the temporally structured multineuronal activity is decreased during more awake and/or aroused brain states.

## Conclusions

In summary, our data provide a macroscopic characterization of the flow of cortical activity at high temporal resolution, revealing the presence of column-transcending feedforward inhibition, and the propagation of activity from thalamo-recipient layers in the PC to infragranular layers in the NC. We conclude from our experiments that the intra- and intercolumnar processing of information is performed in parallel, and that neurons belonging to the same functional network can be anatomically dispersed across nonadjacent cortical modules.

## Funding

This work was supported by Deutsche Forschungsgemeinschaft (German Research Foundation) grants to H.J.L.

## Notes

We thank Beate Krumm for excellent technical assistance; Jenq-Wei Yang, Shuming An, Alexandros Kyriakatos and Evgeny Bobrov for their help on experimental procedures; and Barbara Imbrosci, Jérôme Mordel and the anonymous reviewers for their helpful comments and suggestions. Conflict of Interest: None declared.

## References

Abeles
M
Gerstein
GL
.
1988
.
Detecting spatiotemporal firing patterns among simultaneously recorded single neurons
.
J Neurophysiol
.
60
:
909
924
.
Aertsen
A
Diesmann
M
Gewaltig
MO
.
1996
.
Propagation of synchronous spiking activity in feedforward neural networks
.
J Physiol Paris
.
90
:
243
247
.
Ahissar
E
Sosnik
R
Haidarliu
S
.
2000
.
Transformation from temporal to rate coding in a somatosensory thalamocortical pathway
.
Nature
.
406
:
302
306
.
Arieli
A
Sterkin
A
Grinvald
A
Aertsen
A
.
1996
.
Dynamics of ongoing activity: explanation of the large variability in evoked cortical responses
.
Science
.
273
:
1868
1871
.
Armstrong-James
M
Fox
K
Das-Gupta
A
.
1992
.
Flow of excitation within rat barrel cortex on striking a single vibrissa
.
J Neurophysiol
.
68
:
1345
1358
.
Blinder
P
Tsai
PS
Kaufhold
JP
Knutsen
PM
Suhl
H
Kleinfeld
D
.
2013
.
The cortical angiome: an interconnected vascular network with noncolumnar patterns of blood flow
.
Nat Neurosci
.
16
:
889
897
.
Bonifazi
P
Goldin
M
Picardo
MA
Jorquera
I
Cattani
A
Bianconi
G
Represa
A
Ben-Ari
Y
Cossart
R
.
2009
.
GABAergic hub neurons orchestrate synchrony in developing hippocampal networks
.
Science
.
326
:
1419
1424
.
Boucsein
C
Nawrot
M
Schnepel
P
Aertsen
A
.
2011
.
Beyond the cortical column: abundance and physiology of horizontal connections imply a strong role for inputs from the surround
.
Front Neurosci
.
5
:
32
.
Brunel
N
Wang
XJ
.
2003
.
What determines the frequency of fast network oscillations with irregular neural discharges? I. Synaptic dynamics and excitation-inhibition balance
.
J Neurophysiol
.
90
:
415
430
.
Bruno
RM
Sakmann
B
.
2006
.
Cortex is driven by weak but synchronously active thalamocortical synapses
.
Science
.
312
:
1622
1627
.
Bruno
RM
Simons
DJ
.
2002
.
Feedforward mechanisms of excitatory and inhibitory cortical receptive fields
.
J Neurosci
.
22
:
10966
10975
.
Buzsaki
G
Anastassiou
CA
Koch
C
.
2012
.
The origin of extracellular fields and currents - EEG, ECoG, LFP and spikes
.
Nat Rev Neurosci
.
13
:
407
420
.
Carlson
GC
Coulter
DA
.
2008
.
In vitro functional imaging in brain slices using fast voltage-sensitive dye imaging combined with whole-cell patch recording
.
Nat Protoc
.
3
:
249
255
.
Chagnac-Amitai
Y
Luhmann
HJ
Prince
DA
.
1990
.
Burst generating and regular spiking layer 5 pyramidal neurons of rat neocortex have different morphological features
.
J Comp Neurol
.
296
:
598
613
.
Constantinople
CM
Bruno
RM
.
2013
.
Deep cortical layers are activated directly by thalamus
.
Science
.
340
:
1591
1594
.
Csicsvari
J
Henze
DA
Jamieson
B
Harris
KD
Sirota
A
Bartho
P
Wise
KD
Buzsáki
G
.
2003
.
Massively parallel recording of unit and local field potentials with silicon-based electrodes
.
J Neurophysiol
.
90
:
1314
1323
.
de Kock
CP
Bruno
RM
Spors
H
Sakmann
B
.
2007
.
Layer- and cell-type-specific suprathreshold stimulus representation in rat primary somatosensory cortex
.
J Physiol
.
581
:
139
154
.
de Kock
CP
Sakmann
B
.
2009
.
Spiking in primary somatosensory cortex during natural whisking in awake head-restrained rats is cell-type specific
.
.
106
:
16446
16450
.
Diesmann
M
Gewaltig
MO
Aertsen
A
.
1999
.
Stable propagation of synchronous spiking in cortical neural networks
.
Nature
.
402
:
529
533
.
Doischer
D
Aurel Hosp
J
Yanagawa
Y
Obata
K
Jonas
P
Vida
I
Bartos
M
.
2008
.
Postnatal differentiation of basket cells from slow to fast signaling devices
.
J Neurosci
.
28
:
12956
12968
.
Egger
R
Narayanan
RT
Helmstaedter
M
de Kock
CPJ
Oberlaender
M
.
2012
.
3D reconstruction and standardization of the rat vibrissal cortex for precise registration of single neuron morphology
.
PLOS Comput Biol
.
8
:
12
.
Einevoll
GT
Franke
F
Hagen
E
Pouzat
C
Harris
KD
.
2012
.
Towards reliable spike-train recordings from thousands of neurons with multielectrodes
.
Curr Opin Neurobiol
.
22
:
11
17
.
Feldmeyer
D
.
2012
.
Excitatory neuronal connectivity in the barrel cortex
.
Front Neuroanat
.
6
:
24
.
Feldmeyer
D
Brecht
M
Helmchen
F
Petersen
CCH
Poulet
JFA
Staiger
JF
Luhmann
HJ
Schwarz
C
.
2013
.
Barrel cortex function
.
Prog Neurobiol
.
103
:
3
27
.
Ferezou
I
Bolea
S
Petersen
CC
.
2006
.
Visualizing the cortical representation of whisker touch: voltage-sensitive dye imaging in freely moving mice
.
Neuron
.
50
:
617
629
.
Friedberg
MH
Lee
SM
Ebner
FF
.
1999
.
Modulation of receptive field properties of thalamic somatosensory neurons by the depth of anesthesia
.
J Neurophysiol
.
81
:
2243
2252
.
Fries
P
Reynolds
JH
Rorie
AE
Desimone
R
.
2001
.
Modulation of oscillatory neuronal synchronization by selective visual attention
.
Science
.
291
:
1560
1563
.
Gentet
LJ
Kremer
Y
Taniguchi
H
Huang
ZJ
Staiger
JF
Petersen
CCH
.
2012
.
Unique functional properties of somatostatin-expressing GABAergic neurons in mouse barrel cortex
.
Nat Neurosci
.
15
:
607
612
.
Gilbert
CD
Wiesel
TN
.
1979
.
Morphology and intracortical projections of functionally characterised neurones in the cat visual cortex
.
Nature
.
280
:
120
125
.
Gray
CM
PE
Wilson
M
McNaughton
B
.
1995
.
Tetrodes markedly improve the reliability and yield of multiple single-unit isolation from multi-unit recordings in cat striate cortex
.
J Neurosci Met
.
63
:
43
54
.
Grun
S
.
2009
.
Data-driven significance estimation for precise spike correlation
.
J Neurophysiol
.
101
:
1126
1140
.
Guillery
RW
Sherman
SM
.
2002
.
Thalamic relay functions and their role in corticocortical communication: generalizations from the visual system
.
Neuron
.
33
:
163
175
.
Haider
B
Hausser
M
Carandini
M
.
2013
.
Inhibition dominates sensory responses in the awake cortex
.
Nature
.
493
:
97
100
.
Harris
KD
Henze
DA
Csicsvari
J
Hirase
H
Buzsáki
G
.
2000
.
Accuracy of tetrode spike separation as determined by simultaneous intracellular and extracellular measurements
.
J Neurophysiol
.
84
:
401
414
.
Havenith
MN
Yu
S
Biederlack
J
Chen
NH
Singer
W
Nikolic
D
.
2011
.
Synchrony makes neurons fire in sequence, and stimulus properties determine who is ahead
.
J Neurosci
.
31
:
8570
8584
.
Hazan
L
Zugaro
M
Buzsáki
G
.
2006
.
Klusters, NeuroScope, NDManager: a free software suite for neurophysiological data processing and visualization
.
J Neurosci Methods
.
155
:
207
216
.
Hill
DN
Mehta
SB
Kleinfeld
D
.
2011
.
Quality metrics to accompany spike sorting of extracellular signals
.
J Neurosci
.
31
:
8699
8705
.
Hubel
DH
Wiesel
TN
.
1959
.
Receptive fields of single neurones in the cat's striate cortex
.
J Physiol
.
148
:
574
591
.
Ikegaya
Y
Aaron
G
Cossart
R
Aronov
D
Lampl
I
Ferster
D
Yuste
R
.
2004
.
Synfire chains and cortical songs: temporal modules of cortical activity
.
Science
.
304
:
559
564
.
Isaacson
JS
Scanziani
M
.
2011
.
How inhibition shapes cortical activity
.
Neuron
.
72
:
231
243
.
Izhikevich
EM
Gally
JA
Edelman
GM
.
2004
.
Spike-timing dynamics of neuronal groups
.
Cereb Cortex
.
14
:
933
944
.
Kajikawa
Y
Schroeder
CE
.
2011
.
How local is the local field potential?
Neuron
.
72
:
847
858
.
Kerr
JN
de Kock
CP
Greenberg
DS
Bruno
RM
Sakmann
B
Helmchen
F
.
2007
.
Spatial organization of neuronal population responses in layer 2/3 of rat barrel cortex
.
J Neurosci
.
27
:
13316
13328
.
Kimura
F
Itami
C
Ikezoe
K
Tamura
H
Fujita
I
Yanagawa
Y
Obata
K
Ohshima
M
.
2010
.
Fast activation of feedforward inhibitory neurons from thalamic input and its relevance to the regulation of spike sequences in the barrel cortex
.
J Physiol
.
588
:
2769
2787
.
Kumar
A
Rotter
S
Aertsen
A
.
2010
.
Spiking activity propagation in neuronal networks: reconciling different perspectives on neural coding
.
Nat Rev Neurosci
.
11
:
615
627
.
Le Cam
J
Estebanez
L
Jacob
V
Shulz
DE
.
2011
.
Spatial structure of multiwhisker receptive fields in the barrel cortex is stimulus dependent
.
J Neurophysiol
.
106
:
986
998
.
Linden
H
Tetzlaff
T
Potjans
TC
Pettersen
KH
Grun
S
Diesmann
M
Einevoll
GT
.
2011
.
Modeling the spatial reach of the LFP
.
Neuron
.
72
:
859
872
.
Luhmann
HJ
Greuel
JM
Singer
W
.
1990a
.
Horizontal interactions in cat striate cortex: II. a current source-density analysis
.
Eur J Neurosci
.
2
:
358
368
.
Luhmann
HJ
Singer
W
Martínez-Millán
L
.
1990b
.
Horizontal interactions in cat striate cortex: I. anatomical substrate and postnatal development
.
Eur J Neurosci
.
2
:
344
367
.
MacLean
JN
Watson
BO
Aaron
GB
Yuste
R
.
2005
.
Internal dynamics determine the cortical response to thalamic stimulation
.
Neuron
.
48
:
811
823
.
Mao
BQ
Hamzei-Sichani
F
Aronov
D
Froemke
RC
Yuste
R
.
2001
.
Dynamics of spontaneous activity in neocortical slices
.
Neuron
.
32
:
883
898
.
Matsumoto
K
Ishikawa
T
Matsuki
N
Ikegaya
Y
.
2013
.
Multineuronal spike sequences repeat with millisecond precision
.
Front Neural Circuits
.
7
:
112
.
Meier
R
Egert
U
Aertsen
A
Nawrot
MP
.
2008
.
FIND—a unified framework for neural data analysis
.
Neural Netw
.
21
:
1085
1093
.
Meyer
HS
Egger
R
Guest
JM
Foerster
R
Reissl
S
Oberlaender
M
.
2013
.
Cellular organization of cortical barrel columns is whisker-specific
.
.
110
:
19113
19118
.
Meyer
HS
Schwarz
D
Wimmer
VC
Schmitt
AC
Kerr
JN
Sakmann
B
Helmstaedter
M
.
2011
.
Inhibitory interneurons in a cortical column form hot zones of inhibition in layers 2 and 5A
.
.
108
:
16807
16812
.
Meyer
HS
Wimmer
VC
Hemberger
M
Bruno
RM
de Kock
CPJ
Frick
A
Sakmann
B
Helmstaedter
M
.
2010
.
Cell type-specific thalamic innervation in a column of rat vibrissal cortex
.
Cereb Cortex
.
20
:
2287
2303
.
Mitzdorf
U
.
1985
.
Current source-density method and application in cat cerebral cortex: investigation of evoked potentials and EEG phenomena
.
Physiol Rev
.
65
:
37
100
.
Mokeichev
A
Okun
M
Barak
O
Katz
Y
Ben Shahar
O
Lampl
I
.
2007
.
Stochastic emergence of repeating cortical motifs in spontaneous membrane potential fluctuations in vivo
.
Neuron
.
53
:
413
425
.
Mountcastle
VB
.
1997
.
The columnar organization of the neocortex
.
Brain
.
120
:
701
722
.
Z
.
2000
.
Spike sequences and their consequences
.
J Physiol Paris
.
94
:
505
524
.
Z
Hirase
H
Czurko
A
Csicsvari
J
Buzsaki
G
.
1999
.
Replay and time compression of recurring spike sequences in the hippocampus
.
J Neurosci
.
19
:
9497
9507
.
Nawrot
MP
Schnepel
P
Aertsen
A
Boucsein
C
.
2009
.
Precisely timed signal transmission in neocortical networks with reliable intermediate-range projections
.
Front Neural Circuits
.
3
:
1
.
Nicholson
C
Freeman
JA
.
1975
.
Theory of current source-density analysis and determination of conductivity tensor for anuran cerebellum
.
J Neurophysiol
.
38
:
356
368
.
Oberlaender
M
Boudewijns
ZS
Kleele
T
Mansvelder
HD
Sakmann
B
de Kock
CP
.
2011
.
Three-dimensional axon morphologies of individual layer 5 neurons indicate cell type-specific intracortical pathways for whisker motion and touch
.
.
108
:
4188
4193
.
Oberlaender
M
de Kock
CPJ
Bruno
RM
Ramirez
A
Meyer
HS
Dercksen
VJ
Helmstaedter
M
Sakmann
B
.
2012
.
Cell type-specific three-dimensional structure of thalamocortical circuits in a column of rat vibrissal cortex
.
Cereb Cortex
.
22
:
2375
2391
.
O'Connor
DH
Peron
SP
Huber
D
Svoboda
K
.
2010
.
Neural activity in barrel cortex underlying vibrissa-based object localization in mice
.
Neuron
.
67
:
1048
1061
.
Petersen
CCH
Grinvald
A
Sakmann
B
.
2003
.
Spatiotemporal dynamics of sensory responses in layer 2/3 of rat barrel cortex measured in vivo by voltage-sensitive dye imaging combined with whole-cell voltage recordings and neuron reconstructions
.
J Neurosci
.
23
:
1298
1309
.
Petersen
RS
Diamond
ME
.
2000
.
Spatial-temporal distribution of whisker-evoked activity in rat somatosensory cortex and the coding of stimulus location
.
J Neurosci
.
20
:
6135
6143
.
Riehle
A
Grun
S
Diesmann
M
Aertsen
A
.
1997
.
Spike synchronization and rate modulation differentially involved in motor cortical function
.
Science
.
278
:
1950
1953
.
Rolston
JD
Wagenaar
DA
Potter
SM
.
2007
.
Precisely timed spatiotemporal patterns of neural activity in dissociated cortical cultures
.
Neuroscience
.
148
:
294
303
.
Royer
S
Zemelman
BV
Losonczy
A
Kim
J
Chance
F
Magee
JC
Buzsaki
G
.
2012
.
Control of timing, rate and bursts of hippocampal place cells by dendritic and somatic inhibition
.
Nat Neurosci
.
15
:
769
775
.
Sakata
S
Harris
KD
.
2009
.
Laminar structure of spontaneous and sensory-evoked population activity in auditory cortex
.
Neuron
.
64
:
404
418
.
Sakata
S
Harris
KD
.
2012
.
Laminar-dependent effects of cortical state on auditory cortical spontaneous activity
.
Front Neural Circuits
.
6
:
109
.
Sakurai
Y
Takahashi
S
.
2006
.
Dynamic synchrony of firing in the monkey prefrontal cortex during working-memory tasks
.
J Neurosci
.
26
:
10141
10153
.
Schmitzer-Torbert
N
Jackson
J
Henze
D
Harris
K
Redish
.
2005
.
Quantitative measures of cluster quality for use in extracellular recordings
.
Neuroscience
.
131
:
1
11
.
S
Grün
S
Diesmann
M
Gerstein
GL
.
2008
.
Detecting synfire chain activity using massively parallel spike train recording
.
J Neurophysiol
.
100
:
2165
2176
.
Shmiel
T
Drori
R
Shmiel
O
Ben Shaul
Y
Z
Shemesh
M
Teicher
M
Abeles
M
.
2006
.
Temporally precise cortical firing patterns are associated with distinct action segments
.
J Neurophysiol
.
96
:
2645
2652
.
Simons
DJ
Carvell
GE
.
1989
.
Thalamocortical response transformation in the rat vibrissa/barrel system
.
J Neurophysiol
.
61
:
311
330
.
Sirota
A
Montgomery
S
Fujisawa
S
Isomura
Y
Zugaro
M
Buzsáki
G
.
2008
.
Entrainment of neocortical neurons and gamma oscillations by the hippocampal theta rhythm
.
Neuron
.
60
:
683
697
.
Stanley
GB
.
2013
.
Reading and writing the neural code
.
Nat Neurosci
.
16
:
259
263
.
Sun
JJ
Kilb
W
Luhmann
HJ
.
2010
.
Self-organization of repetitive spike patterns in developing neuronal networks in vitro
.
Eur J Neurosci
.
32
:
1289
1299
.
Sun
QQ
Huguenard
JR
Prince
DA
.
2006
.
Barrel cortex microcircuits: thalamocortical feedforward inhibition in spiny stellate cells is mediated by a small number of fast-spiking interneurons
.
J Neurosci
.
26
:
1219
1230
.
HA
.
2003
.
Fast-spike interneurons and feedforward inhibition in awake sensory neocortex
.
Cereb Cortex
.
13
:
25
32
.
Tetko
IV
Villa
AEP
.
1997
.
Fast combinatorial methods to estimate the probability of complex temporal patterns of spikes
.
Biol Cybern
.
76
:
397
407
.
Varga
Z
Jia
H
Sakmann
B
Konnerth
A
.
2011
.
Dendritic coding of multiple sensory inputs in single cortical neurons in vivo
.
.
108
:
15420
15425
.
Villa
AEP
Tetko
IV
Hyland
B
Najem
A
.
1999
.
Spatiotemporal activity patterns of rat cortical neurons predict responses in a conditioned task
.
.
96
:
1106
1111
.
Wester
JC
Contreras
D
.
2012
.
Columnar interactions determine horizontal propagation of recurrent network activity in neocortex
.
J Neurosci
.
32
:
5454
5471
.
Wilent
WB
Contreras
D
.
2004
.
Synaptic responses to whisker deflections in rat barrel cortex as a function of cortical layer and stimulus intensity
.
J Neurosci
.
24
:
3985
3998
.
Womelsdorf
T
Fries
P
Mitra
PP
Desimone
R
.
2006
.
Gamma-band synchronization in visual cortex predicts speed of change detection
.
Nature
.
439
:
733
736
.
Woolsey
TA
Rovainen
CM
Cox
SB
Henegar
MH
Liang
GE
Liu
DQ
Moskalenko
YE
Sui
J
Wei
L
.
1996
.
Neuronal units linked to microvascular modules in cerebral cortex: response elements for imaging the brain
.
Cereb Cortex
.
6
:
647
660
.
Xu
H
Jeong
HY
Tremblay
R
Rudy
B
.
2013
.
Neocortical somatostatin-expressing GABAergic interneurons disinhibit the thalamorecipient layer 4
.
Neuron
.
77
:
155
167
.
Yang
JW
An
SM
Sun
JJ
Reyes-Puerta
V
Kindler
J
Berger
T
Kilb
W
Luhmann
HJ
.
2013
.
Thalamic network oscillations synchronize ontogenetic columns in the newborn rat barrel cortex
.
Cereb Cortex
.
23
:
1299
1316
.
Yang
Y
DeWeese
MR
Otazu
GH
AM
.
2008
.
Millisecond-scale differences in neural activity in auditory cortex can drive decisions
.
Nat Neurosci
.
11
:
1262
1263
.

## Author notes

Vicente Reyes-Puerta and Jyh-Jang Sun contributed equally to this work.