Elementary motion sequence detectors in whisker somatosensory cortex.

How the somatosensory cortex (S1) encodes complex patterns of touch, such as those that occur during tactile exploration, is poorly understood. In the mouse whisker S1, temporally dense stimulation of local whisker pairs revealed that most neurons are not classical single-whisker feature detectors, but instead are strongly tuned to two-whisker sequences that involve the columnar whisker (CW) and one specific surround whisker (SW), usually in a SW-leading-CW order. Tuning was spatiotemporally precise and diverse across cells, generating a rate code for local motion vectors defined by SW-CW combinations. Spatially asymmetric, sublinear suppression for suboptimal combinations and near-linearity for preferred combinations sharpened combination tuning relative to linearly predicted tuning. This resembles computation of motion direction selectivity in vision. SW-tuned neurons, misplaced in the classical whisker map, had the strongest combination tuning. Thus, each S1 column contains a rate code for local motion sequences involving the CW, thus providing a basis for higher-order feature extraction.

Introduction

Tactile exploration generates complex spatiotemporal patterns of touch stimuli, whose representation in somatosensory cortex (S1) is poorly understood. Primary sensory cortex neurons classically integrate across sensory space and time to generate tuning for local spatiotemporal features[1-3]. But in whisker S1, the canonical touch system of rodents, a ‘one-whisker-one-column’ model dominates, in which each facial whisker maps to one cortical column and S1 neurons are considered single-whisker feature detectors for velocity, acceleration or bend of the columnar whisker (CW). Natural whisking generates complex, rapid patterns of multi-whisker contact[4-6], suggesting that S1 may compute tuning for multi-whisker patterns. How multi-whisker features are represented and organized in S1 remains unclear[7].

Several findings suggest that multi-whisker tuning may be prevalent in S1. S1 neurons spike weakly to single-whisker stimuli[8-10] and have broad subthreshold receptive fields suggesting multi-whisker integration[11,12]. In each column, most neurons are somatotopically tuned to the CW[13,14], but a surprising fraction prefer a surround whisker (SW) over the CW, and are thus ‘misplaced’ in the classic whisker map[15,16]. These SW-tuned cells may prefer multi-whisker features. Some neurons are tuned to highly structured multi-whisker features such as global motion wavefronts[17,18] and concentric center-surround motion contrast[19], but these engage a small fraction of neurons or elicit no more spikes than single-whisker stimuli. Tuning for extended whisker sequences exists, but its organization and effect on coding remain unclear[20]. We hypothesized that S1 neurons encode a more elementary spatiotemporal feature via tuning that is highly prevalent, organized systematically across S1 columns, and constructed through cross-whisker nonlinearities[21-26].

We tested for spatiotemporal feature tuning within an elementary subspace of multi-whisker stimuli, comprising local pair-wise whisker combinations and sequences, over a physiological range of deflection intervals (Δt). These represent local motion on the whisker pad, have tractable dimensionality, and are elementary building blocks of multi-whisker patterns. We mapped responses densely within this space, and found that most S1 neurons are 2-whisker sequence detectors. Within each S1 column, neurons preferred a diverse set of CW-containing sequences, usually in SW-leading-CW order, creating a robust population rate code for local motion. This reveals a novel organization for the S1 map, in which each column encodes a comprehensive set of local motion vectors involving the CW.

Results

To assess 2-whisker tuning, we measured spiking of single units in D1 or C1 columns of S1 in urethane-anesthetized mice (P28–45) using multi-site silicon probes. We deflected whiskers in a 3 × 3 array centered on the columnar whisker (CW) and including 8 immediate surround whiskers (SWs). Each whisker was deflected with a 40 ms triphasic waveform that optimally drives S1 neurons[19]. In Experiment 1, we surveyed for 2-whisker combination tuning by presenting single-whisker stimuli and two-whisker stimuli that included all CW-SW and all SW-SW combinations at 7 inter-whisker deflection intervals (Δt = 0, ±10, ±25, ±50 ms) (Fig. 1a). Tuning was analyzed in 142 whisker-responsive single units (out of 160 single units) in 8 recording sites in 5 mice. Units were distributed across S1 layers (L2/3 = 36, L4 = 48, L5 = 44, L6 = 14 units).

Figure 1.

S1 neurons preferentially encode CW-SW sequences. (a) Whisker stimuli applied via a 3×3 whisker array centered on the recorded S1 column (here, D1). Bottom: Whisker stimuli for Experiment 1. (b) Top, Mean single-whisker tuning across all units by layer; shading is SEM (L2/3 = 36, L4 = 48, L5 = 44 units). Bottom, Spiking response (z-scored within each unit) to CW deflection vs. best SW deflection for each unit. Histogram shows distribution of CW-tuned and SW-tuned units. (c) Tuning curves for 2-whisker sequences for 2 example units. Yellow and blue bars show the response to CW-SW or SW-SW combinations; left and right sectors show negative Δt (SW precedes CW) or positive Δt (CW precedes SW). Each bar shows the peak response across Δt’s for that whisker combination. Black bars show single-whisker responses. (d) Response of each unit to its best 2- whisker sequence vs. best single whisker. (e) Fraction of units whose best sequence contained the CW (n=89 CW-tuned and 53 SW-tuned units). Error bar is SE of the sample proportion. (f) Spiking response (z-scored within each unit) to best CW-SW vs. best SW-SW sequence for each unit. (g) Tuning sharpness, quantified as lifetime sparseness, among different sets of whisker sequences. Within each set, sequences were ranked from strongest to weakest response, and lifetime sparseness was calculated for increasing number of sequences ranked. Thus, tuning sharpness can be compared between the N best CW-SW and N best SW-SW sequences (at X=N on the x-axis). Shaded regions are SEM.

Single-whisker receptive fields in each S1 column

Classical whisker receptive fields, assessed from single-whisker deflections, were heterogeneous in each column. On average across units, the CW evoked the most spikes, followed by within-row SWs, matching standard S1 somatotopy[14] (Fig. 1b, top). But tuning was diverse at the single unit level, even within the same penetration (Supplementary Fig. 1). 61% of units were driven most strongly by the CW (CW-tuned units, n = 89), and 39% by an SW (SW-tuned units, n = 53). Some SW-tuned units showed nearly equivalent responses to the CW, but others strongly preferred the SW (Fig. 1b, bottom). CW- and SW-tuned neurons had comparable spontaneous activity (1.28±0.24 vs. 2.08±0.4 Hz) and maximal evoked firing (4.8±0.48 vs. 4.56±0.64 spikes/stimulus), but tuning was somewhat broader for SW-tuned neurons (Supplementary Fig. 1). Thus, a single column contains mostly CW-tuned neurons, but also a substantial fraction of SW-tuned neurons, as observed previously[15,16,27].

S1 units encode local motion involving multi-whisker sequences

Tuning for whisker combinations was analyzed from CW-SW and SW-SW stimuli, initially independent of Δt. Fig. 1c shows two single units in the D1 column. Black bars show the response to each whisker deflected singly; colored bars show the response to each whisker combination at whichever Δt (0, 10, 25, 50 ms) drove the strongest response to that combination. Responses are spike counts in a broad, 125-ms window that includes all spikes evoked by both whiskers in the combination. Each unit spiked more to a specific 2-whisker combination than to other stimuli. Note that while one unit was CW-tuned and one was SW-tuned, both preferred a combination that contained the CW (D1). Of 142 units, nearly all spiked more to the best combination than the best single whisker (slope = 1.26, Fig. 1d). The best combination included the CW for 70% of all units, 85% of CW-tuned units, and 47% of SW-tuned units (chance = 22%) (Fig. 1e). Thus, the best CW-SW sequence evoked more spikes than the best SW-SW sequence (Fig. 1f) or the best single-whisker stimulus (Fig. 1d).

Whisker combinations are also sequences with a specific Δt. Neurons showed sharper tuning selectivity among CW-SW than SW-SW sequences (Fig. 1g). Selectivity was measured as lifetime sparseness[28] across sequences varying in both spatial whisker combination and Δt. Lifetime sparseness equals 0 when a neuron spikes equally to all sequences and 1 when a neuron spikes exclusively to one sequence, and thus represents spatiotemporal selectivity and discriminative capacity among sequences[29,30]. SW-tuned neurons had particularly sharp tuning for CW-SW sequences (Supplemental Fig. 2). Thus, neurons in each S1 column are preferentially tuned for, and discriminate, CW-SW sequences.

Dense spatiotemporal mapping of CW-SW sequence tuning

To characterize CW-SW sequence tuning with greater Δt resolution, we performed Experiment 2. Stimuli included all CW-SW combinations with Δt uniformly distributed between ±50 ms at 1 ms intervals. Because some S1 neurons are sensitive to CW-SW motion contrast[19], we presented both same- and opposite-direction CW-SW deflections (termed correlated and anti-correlated stimuli). Single whisker stimuli and a subset of SW-SW sequences were also applied (Fig. 2a). The average whisker deflection rate was ~3 Hz. 247 single-units were isolated in 12 mice, of which 224 were responsive to correlated and 227 to anti-correlated stimuli (L2/3: 54, 54 units [correlated, anti-correlated]; L4: 78, 79 units; L5: 71, 70 units; L6: 21, 24 units).

Figure 2.

Dense spatiotemporal mapping of CW-SW sequences reveals combination tuning in S1. (a) Stimulus set for Experiment 2. Whisker combinations were delivered at Δt distributed between ±50ms at 1 ms intervals. Top right, example CW-SW stimulus waveforms at +10ms Δt. (b) Example unit recorded in D1 column to correlated CW-SW stimuli. Each raster plot shows spiking across all Δt’s for each combination, with the D2-D1 combination enlarged below. The Δt tuning curve was calculated by averaging spikes in the spike count window for each Δt and smoothing (see Methods). (c) Columnar-whisker spatiotemporal receptive field (cwSTRF) for the unit in b. CSI (tuning sharpness across whisker combinations) was calculated for each Δt. (d) Combination tuning curves, rasters and PSTHs for 4 example combination-selective units, calculated within ± 5ms of each unit’s best Δt. Asterisk, best combination. The top-left unit is the same as in (b) and (c). (e) Distribution of CSI for significantly combination selective units (shuffle test, α ≤ .05, FDR control). (f) Average combination tuning for all combination-selective units (n = 187), and for the 25% most- and least-selective combination-selective units. Tuning curves were aligned by their peak response (asterisk). Gray, tuning after spike counts were shuffled across combinations. Green, spontaneous firing rate. Shaded regions are SEM. (g) Distribution of best CW-SW combination identity was non-uniform (Rayleigh test, p=4e-19; n = 451 units). Each bin is one CW-SW combination, organized by SW position on the face (D=dorsal, V=ventral, R=rostral, C=caudal).

For each single unit, we applied a model-based smoothing method on evoked spike counts to calculate a Δt tuning function for each whisker combination (see Methods) (Fig. 2b). The spatiotemporal response function (cwSTRF) summarizes the mean response to each CW-SW combination at each Δt, after Δt smoothing (Fig. 2c). The ‘best whisker combination’ and ‘best Δt’ were defined as those that jointly evoked the maximum response. We calculated a combination selectivity index (CSI) at each Δt as lifetime sparseness across the 8 CW-SW combinations at that Δt. Correlated and anti-correlated stimuli were analyzed separately for each single-unit, therefore we refer to each cwSTRF as a unit. A unit was considered combination selective if the CSI trace across Δt’s was significantly different than for spike count data shuffled across stimuli and trials (see Methods). All p-values are adjusted for multiple comparisons across units.

Combination tuning is strong and prevalent in S1

52% of S1 units (n=129) were significantly combination-selective, based on either correlated stimuli (35% of units) or anti-correlated stimuli (46% of units). Four example units are shown in Fig. 2d. The polar plots show the measured spiking response to different CW-SW combinations at the best Δt, relative to CW alone and to the linear sum of single CW and SW deflections. The latter is used as a null hypothesis for default tuning. All 4 units preferred a specific CW-SW sequence (asterisk) which evoked more spikes than the CW alone. Tuning for CW-SW sequences was consistently sharper than the linear prediction. Many S1 units showed sharp combination tuning (Supplemental Fig. 3).

We analyzed spatial combination tuning at best Δt for each unit. Combination-selective units varied in tuning sharpness (Fig. 2e) and maximal evoked firing rate (Supplemental Fig. 4). Across all combination-selective units, average tuning for CW-SW sequences was sharp, as shown by aligning combination-tuning curves to their peak (by rotating the polar plots). This tuning was real, because shuffling spike count data across stimuli and trials abolished tuning structure (Fig. 2f). Among combination-selective units, the most selective units had lower spontaneous and whisker-evoked firing than the least selective units (Fig. 2f). This relationship between firing rate and tuning selectivity is expected because spontaneous spiking is additive across the tuning curve and weakens tuning sharpness[30,31]. Substantial combination selectivity was also found among the most responsive combination-selective units and after subtraction of spontaneous activity (Supplemental Fig. 4). Thus, combination tuning is robust and present across a wide range of firing rates. At the population level, preferred stimuli included all CW-SW combinations, but with a strong bias toward the rostral within-row SW and SWs ventral to the CW (Fig. 2g). This also is evident in mean somatotopically aligned tuning (Supplemental Fig. 5), and may represent ethologically relevant axes of local motion.

Nonlinearities enhance rate coding of whisker combinations

CW-SW combination tuning was sharper than predicted from the linear sum of single-whisker responses for nearly all units (p=1e-50, paired t-test) (Fig. 3a), including both combination-selective and non-selective units (Fig. 3b). This was readily apparent in single example cells (Fig. 2d), and in average aligned combination tuning curves, which show strong sublinearity for non-optimal combinations, and near-linearity for the optimal combination (Fig. 3c–d). Combination selectivity existed in all layers, but was stronger in L2–4 than L5–6 (1-factor ANOVA, p=.0025) (Supplemental Fig. 6). Thus, tuning for CW-SW combinations is sculpted from broader tuning predicted by single-whisker input.

Figure 3.

Nonlinear sharpening of combination tuning enhances rate coding for whisker combinations. (a) Comparison of measured and linearly predicted CSI at best Δt. (b) Distributions of measured and linearly predicted CSI for combination-selective (n=187) and non-selective units (n=264). Inset show mean CSI±SEM, 2 sample t-test, p<0.005. (c-d) Average peak-aligned combination tuning across all combination-selective units (c), or the 25% most-selective and 25% least-selective combination-selective units (d). Asterisk indicates peak response. Shading is SEM. Sample sizes as in 2f. Black, mean response to the CW deflected individually. (e) Average performance of a neural population decoder that predicts CW-SW combination identity based on single-trial spiking activity of combination selective units. The decoding model for each ensemble size was fit by randomly sampling (with replacement) 1000 times from recorded units. (f) Mean confusion matrix of neural decoder with 187 combination-selective units. Entries along the diagonal are percent correct classifications for each CW-SW combination and rows sum to 1. Orange and gray bars are % correct classification for measured and linear responses respectively. (g) Average performance of a neural population decoder performing a binary discrimination between best and suboptimal CW-SW combinations. Stimuli for training and testing were either at each neuron’s best Δt, or in specific Δt ranges. Shaded regions and error bars in (e-g) are SD across 2500 bootstrap decoding trials.

The prevalence and diversity of combination tuning suggests a population rate code for whisker combinations in each column. We tested whether a neural population decoder could discriminate CW-SW combination identity (among the 8 possible CW-SW combinations) using spike count data from single trials. Each neuron in Experiment 2 was modeled as a one-vs-all classifier to predict the probability of each CW-SW combination from single-trial spike count data (see Methods). Training and testing data were from combinations with Δt ≤ ±5 ms from the unit’s best Δt. An ensemble of N neurons was sampled with replacement from all combination-selective units. The CW-SW combination with the highest summed probability across neurons was taken as the population prediction.

The average single-trial accuracy (computed across all 8 CW-SW combinations) increased with N, reaching 42% for N = 187. This is substantially higher than chance performance of 12.5% (1/8) (Fig. 3e). Accuracy was highest for the CW-SW combination involving the rostral within-row SW (82% for N=187), and lowest for the 3 dorsal-row SWs (23–28%) (Fig. 3f), matching the somatotopic bias of combination tuning. Training and testing the decoder on linearly predicted spike counts resulted in significantly weaker discrimination (Fig. 3e–f).

We also tested a separate, binary decoder trained to discriminate between each unit’s best CW-SW combination and suboptimal combinations (at the unit’s best Δt) using single-trial spike count data (Fig. 3g). This decoder achieved 70% accuracy for N=10 neurons and >90% for N=100 neurons, comparable to neural population decoders for texture discrimination[32,33] and stick/slip detection[8]. This suggests that higher-order neurons can achieve substantial decoding performance by pooling over 10 similarly tuned neurons. Performance was also strong for stimuli at specific Δt, rather than at each unit’s best Δt. Thus, the neural ensemble in a single S1 column accurately encodes CW-SW combinations on single trials and nonlinearities enhance the efficiency and accuracy of this code.

SW-tuned units enhance population coding of whisker combinations

SW-tuned units are ‘misplaced’ in the classical single-whisker columnar map, but may contribute importantly to multi-whisker feature coding. SW-tuned units were more sharply combination-tuned than CW-tuned units (p=3.4e-4, t-test) (Fig. 4a–b), consistent with Experiment 1 (Supplemental Fig. 2). This was achieved by strong sublinear sculpting of linearly predicted tuning (Fig. 4b). To determine how SW-tuned units contribute to population coding, we constructed 2 separate neural decoders, one based on CW-tuned units, and the other on SW-tuned units. Whisker combination identity could be decoded from SW-tuned units more accurately and efficiently than CW-tuned units, though the latter were still well above chance (Fig. 4c–d). Thus, SW-tuned units are sharply tuned for CW-SW combinations and enhance population coding of combinations in each column.

Figure 4.

SW-tuned units show strong rate coding for whisker combinations. (a) Probability distributions of combination-selectivity indices for CW (n = 169) and SW-tuned units (n = 282). Insets show mean CSI±SEM, *** 2-sample t-test, p=3.4e-4. (b) Mean normalized CW-SW combination tuning for these same CW and SW-tuned units. Tuning curves were aligned at the best combination. Asterisk indicates peak response. Shaded regions are SEM. (c) Average performance of neural decoders based on either CW- or SW-tuned units in predicting CW-SW combination identity. Decoders were built as in Fig. 3e. Shaded regions are SD across 2500 bootstrap trials. (d) Confusion matrix for the decoder trained on SW-tuned units only. Plotted as in Fig. 3f. Red and blue bars denote mean for SW and CW-tuned units respectively. Error bars are SD across 2500 bootstrap trials

Spatiotemporally specific linear integration enhances combination tuning in S1

CW-SW combination tuning is sculpted from broader, linearly predicted tuning which is, on average, weakly biased towards the best CW-SW combination (Fig. 3c). This linearly predicted tuning reflects asymmetry in the single-whisker receptive field. We analyzed how these two processes—non-uniform sublinear sharpening and the bias in predicted linear tuning—contribute to combination tuning.

To characterize the sharpening process, we quantified linearity of responses to CW-SW combinations at the best Δt. Across all stimuli and all units, combination responses were well predicted by a 64% scaling of the linear sum of single-whisker responses (R2=0.90) (Fig. 5a). Within this average relationship however, non-optimal combinations evoked strongly sublinear responses, while each unit’s best combination evoked a near-linear response (Fig. 5a). Non-optimal combinations had a linearity index of 0.63±0.003 (defined as measured response / linear prediction), and best combinations 0.95±0.02 (p=1.3e-147, 2 sample t-test) (Fig. 5b–c). This difference is not the trivial result of best combinations evoking stronger responses relative to the population linear fit in Fig. 5a. This was shown by a shuffling procedure in which pseudo-units were constructed by sampling 8 random points from Fig. 5a. The strongest responses of these pseudo-units (‘best combinations’) and the remaining, weaker responses (‘non-optimal combinations’) had very similar linearity (0.71 vs. 0.67) (Fig. 5c). Sublinearity for non-optimal combinations was not found just at the unit’s best Δt, but also at all Δt for these whisker combinations (Fig. 5c). Thus, near-linear integration is confined to the best combination and best Δt. Only 8% of units had significant supralinearity for the best combination. Thus, near-linear integration for the best combination and best Δt, plus sublinear suppression for other stimuli, enhances CW-SW combination tuning.

Figure 5.

Best combination responses at best Δt are enhanced relative to a global sublinear scaling. (a) Comparison of measured and linearly predicted combination responses at each unit’s best Δt. Black line is the line of best fit with the y-intercept set to 0. (n = 3,608 combinations, with 7 suboptimal and 1 optimal combination per unit). (b) Left, comparison between linearity (measured response divided by the linear prediction) and z-scored CW-SW combination response. Right, distribution of linearity indices for best and suboptimal combinations. Dashed lines are the mean linearity indices. (c) Average linearity indices (measured response divided by the linear prediction) for best and suboptimal combinations at best Δt for measured and shuffled data. The peak linearity index across all Δt’s is also shown for suboptimal stimuli. *** p < 0.001, 2-sample t-test. Panels (b) and (c) quantify the individual data points in (a).

Linear and nonlinear contributions to combination selectivity

To define the computations underlying combination selectivity, we separated linear and nonlinear components of this process mathematically. We described combination-tuning as resulting from two sequential operations: i) a global scaling to 0.64 of linearly predicted responses, analogous to linear summation with divisive normalization; and ii) a combination-specific, nonlinear enhancement or suppression relative to this global sublinear 0.64x scaling (Fig. 6a). These components are useful to analyze separately but are not meant to represent distinct physiological processes.

Figure 6.

Linear and nonlinear computations underlying combination tuning. (a) Schematic showing global sublinear scaling (0.64x normalization), followed by combination-specific enhancement or suppression, that represent linear and nonlinear components for synthesis of combination tuning. See text for details. The polar tuning curve for an example combination-selective unit summarizes these steps (linear summation->scaling->combination-specific enhancement or suppression ->measured response). (b) Mean rank-ordered combination tuning curves for combination-selective (n = 187) and non-selective units (n = 264). Bottom plots show measured tuning, linear prediction, and the 0.64x-scaled linear prediction. Color legend is the same as (a). Top, combination-specific nonlinearity, defined as residual between measured response and 0.64x-scaled linear prediction. Unpaired t-test, * p=0.001 and p=0.002. (c) Mean combination-specific enhancement or suppression as a function of ranked strength of SW single-whisker response for combination-selective (n = 187) and non-selective units (n = 264). Error bars in (b-c) are SEM.

We compared linear and non-linear components in combination-selective vs. non-selective units. We ranked CW-SW stimuli by measured response strength in each unit, and then calculated the mean ranked combination tuning curve across units. For both combination-selective and non-selective units, the measured tuning curve was sharper than the 0.64-scaled linear prediction from single-whisker deflections (Fig. 6b, bottom). For both groups, the best combination was greater than the 0.64-scaled linear prediction (combination-selective: p=4.12e-25; non-selective: p=1.82e-39, paired t-test) and the worst combination was weaker (p=4.0e-21; p=1.7e-17). This sharpening cannot reflect additional, uniform subtractive inhibition, because this would simply shift the linear prediction along the y-axis. Thus, sharpening of combination tuning must involve combination-specific nonlinearities.

We estimated the magnitude of these combination-specific nonlinearities as the residual between measured tuning and the mean 0.64-scaled linear prediction (Fig. 6b, top). Both combination-selective and non-selective units showed significant enhancement for the best CW-SW combination, and modest suppression for the weakest CW-SW combinations, relative to the global nonlinearity. The magnitude of enhancement did not differ between groups (p=0.24, t-test), and suppression was only modestly greater in selective units (p=0.001 and p=0.002 for the two weakest CW-SW combinations, 2 sample t-test, 2 tailed Bonferroni correction). Instead, combination-selective and non-selective groups differed in two major ways. First, the linearly predicted receptive field was sharper in combination-selective cells (Fig. 3b; Fig. 6b). Second, when we examined the relationship between non-linear enhancement or suppression and the ranked strength of single-whisker SW responses, combination-selective units nonlinearly enhanced combinations that involved the strongest single SW (p=1.2e-8, 1-factor ANOVA), while non-selective units showed an opposite, much weaker trend (p=0.03) (Fig. 6c). Thus, combination-selective units achieve sharp combination tuning by having sharper linearly predicted tuning (i.e., spatially asymmetric single-whisker tuning) (Fig. 3c) and by enhancing responses relative to the global sublinearity for combinations involving the strongest SW (or, equivalently, shielding these responses from sublinear suppression). The individual measurements underlying these findings are in Supplemental Fig. 7.

Tuning for Δt reveals space-time inseparability

Combination-selective neurons were also tuned for Δt. For best CW-SW combinations, these units had sharp Δt tuning, with an average half-width of ~10 ms around the best Δt. Suboptimal combinations had weaker, broader Δt tuning around the units’ best Δt (Fig. 7a). Mean Δt tuning, calculated after standard Δt smoothing, showed a minimum response at −10 ms Δt, consistent with previous work on cross-whisker suppression (Fig. 7b). Best combinations evoked maximal responses at −30 to −50 ms Δt (SW leads CW, corresponding to inbound motion toward the CW), while suboptimal combinations evoked weaker responses overall, with maximal responses at +30 to +50 ms Δt (outbound motion). Thus, combination-selective units have inseparable spatial and temporal tuning and are joint spatiotemporal feature detectors.

Figure 7.

Δt tuning and preference for tactile sequences inbound to the CW. (a) Mean Δt tuning curves for combination-selective units (n = 187) aligned at best Δt. All Δt tuning curves were normalized to each unit’s peak response. Shaded regions are SEM. (b) Mean z-scored Δt tuning curves for best (n = 187) and suboptimal (n = 1309) CW-SW combinations. Z-scoring was done separately for each unit. (c) Cumulative variance explained by the principal components (PCs) of z-scored Δt tuning curves. (d) First three PCs. (e) Distribution of weights (or score) for PC 1 across best and suboptimal combinations, within combination-selective and non-selective units. Plots for PC 2 and 3 do not separate suboptimal stimuli by unit type. (f) CSI and alignment indices for combination-selective units that had a negative (n = 121) or positive best Δt (n = 56). Open circles are individual units, red line is mean and shaded region is SEM. The alignment index was the Pearson correlation coefficient between linearly predicted combination responses and combination-specific facilitation or suppression. 2-sample t-test, * p=0.03 (left) and p=0.01 (right). (g) Decoding of combination identity (as in Fig. 4) across different ranges of Δt. Shaded region is SD across 1000 bootstrap trials.

Order selectivity dominates Δt tuning

To understand the diversity of Δt tuning we performed principal components (PC) analysis on all significantly modulated Δt tuning curves (n= 357 tuning curves representing 79% of best combinations, and 1247 tuning curves representing 40% of suboptimal combinations). Tuning curves were z-scored before analysis. Three PCs captured 75% of variance in Δt tuning (Fig. 7c). PC1 was sigmoidal and centered around 0 ms, and thus captures differences in spike rate related to CW-SW order (Δt>0 vs. Δt<0). PC2 was U-shaped and nearly symmetrical around 0 ms, and thus captures differences in spike rate related to |Δt| (sequence speed). PC3 was a sharper, bimodal filter (Fig. 7d). In combination-selective cells, best and supobtimal combinations had very different PC1 weights (p=5.7e-16, t-test), with negative PC1 weights for best combinations and positive PC1 weights for suboptimal combinations. This is consistent with inbound order preference for best combinations. This pattern was nearly absent in non-selective cells (Fig. 7e). PC2 and PC3 weighting did not differ across cell or stimulus types (Fig. 7e). Thus, the Δt tuning feature that is most related to combination-selectivity is the inbound order preference for best combination stimuli in selective units. A neural decoder was able to report order from single-trial firing rates, indicating that sequence order is indeed represented in each S1 column (Supplemental Fig. 8).

Combination-selective units that preferred inbound order (best Δt < 0) were more sharply spatially tuned than those that preferred outbound order (p=0.03, t-test) (Fig. 7f). Inbound-preferring units also had better alignment between nonlinear and linear components of combination tuning (p=0.01, t-test), which may explain their sharper tuning (Fig. 7f). When the spatial combination decoder (Fig. 4) was retrained and retested using stimuli within specific Δt ranges, spatial combination decoding was better for inbound order stimuli than outbound order stimuli (Fig. 7g). Thus, the inseparability of temporal and spatial tuning yields improved spatial discrimination for inbound stimuli, suggesting that S1 preferentially encodes local spatial features of inbound motion to each column.

Combination tuning in S1 of awake mice

To test for combination tuning in wakefulness, we performed Experiment 3 in awake, head-fixed mice. Whisker stimuli were delivered via the piezo array. Stimuli were presented densely (every 0.42 ± 0.01 s) during 11.5-sec trials, and included all single whiskers plus all correlated CW-SW and SW-SW combinations at Δt = −40, 0 and +40 ms. To ensure alertness, mice were trained to lick for a water reward in response to a visual stimulus presented at the end of ~40% of trials (every 32.5 ± 1.38 s) (Fig. 8a). 4 mice were trained and performed the task during S1 recording, for 1–2 sessions of 1.28 ± 0.13 hours each. Mice licked to most visual stimuli (Go trials), but not during equivalent periods of non-visual (NoGo) trials (Fig. 8b). The overall lick rate during whisker stimulation periods was low (0.32 ± 0.10 Hz). We used a passive whisker task, because active whisking is incompatible with precise multi-whisker stimulus delivery. We recorded spiking from 81 well-isolated single units in S1.

Figure 8.

Combination-selectivity measured in awake S1. (a) Behavioral trial structure. 1- and 2-whisker stimuli are densely interleaved in each trial. On ~40% of trials, mice must detect a brief (2-ms) visual stimulus, to ensure alertness. (b) Behavioral performance. Mice lick in the response window after visual stimuli, but much less in no-visual stimulus trials. Error bars are SE of sample proportion. (c) Mean ranked whisker receptive field for CW-tuned units (n = 46) in awake mice (Experiment 3) vs. anesthetized (n = 89) (Experiment 1). 2-sample t-test, ** p < .001. Shaded regions are SEM across units. Inset, average single-whisker receptive field across all whisker-responsive single units (n=74) in awake mice. (d) Evoked response to best sequence vs. best single whisker in awake mice. Slope of line = 1.24. (e) CW-SW combination tuning for 3 example combination-selective units, compared to linear predicted tuning. Shaded regions are SEM across trials. (f) Mean peak-aligned CW-SW tuning curve across all combination-selective units (n=29), after normalizing to peak response in each unit. This is compared to linear predicted tuning or shuffled data representing no tuning. (g) Measured vs predicted linear response magnitude for each CW-SW combination at each unit’s best Δt. Shaded regions are SEM (h) Linearity index for best (n = 74) and suboptimal CW-SW combinations (n = 518) in (g). Shaded regions are SEM. 2-sample t-test, ***p=6.2e-17. (i) CW-SW tuning is sharper than linear predicted tuning for most units. (j) Peak CW-SW response at each Δt for each unit (n = 74). Responses were z-scored within each unit across all CW-SW combinations and Δt’s. Shading is SEM. (k) Decoding of best vs suboptimal CW-SW combinations using the binary decoder, for ensembles of combination-selective units. Chance performance is 50%. Line and shading are mean and SEM across 2500 bootstrap trials.

Recordings were made in L2/3 and L4 of C1, D1 and D2 columns. Average single-whisker tuning was sharp and peaked at the CW (Fig. 8c). Of whisker-responsive single units (n=74), 62% (46) were tuned to the CW (measured from single whisker deflections), similar to anesthetized mice. Among CW-tuned units, CW-evoked spiking was ~2-fold stronger, and sharper, than in anesthetized mice from Experiment 1 (Fig. 8c). Thus, signal-to-noise for sensory coding appears improved in awake mice.

Despite stronger single-whisker tuning, tuning for 2-whisker combinations was identical to anesthetized mice. The best 2-whisker sequence evoked more spikes than the best single-whisker stimulus in nearly all units (Fig. 8d). 39% of responsive units were significantly combination-selective (similar to 35% of units for correlated stimuli in Experiment 2). 76% of combination-selective units (and 80% of all responsive units) preferred a 2-whisker sequence that contained the CW. Combination tuning, measured at the best Δt for each unit, was narrower than predicted linearly, both for individual units (Fig. 8e) and in aligned tuning across all combination-selective units (Fig. 8f). Measured responses to the best CW-SW combination for each cell was roughly linear (linearity index 0.93 ± 0.03), while responses to suboptimal combinations were sublinear (linearity index 0.71 ± 0.009, p=6.2e-17 relative to best combinations, 2-sample t-test) (Fig. 8g–h). Combination tuning was sharper than linearly predicted tuning for virtually all responsive units (p = 9.6e-11 paired t-test) (Fig. 8i). Responses to inbound CW-SW sequences (Δt = −40 ms) were stronger than to outbound sequences (Δt = +40 ms) or simultaneous deflections (Δt = 0ms) (p = 2.2e-14, 1-factor ANOVA) (Fig. 8j). Decoding of best vs. suboptimal CW-SW combination identity from single-trial firing rate data was accurate and efficient (Fig. 8k). Thus, CW-SW combination tuning, and its computation by asymmetric sublinear suppression, were similar compared to anesthetized mice.

Discussion

We found highly prevalent tuning for local 2-whisker sequences in S1, sculpted from linearly predicted tuning by spatially asymmetric suppression. Prior studies reported prominent cross-whisker suppression and occasional facilitation in S1, and sensitivity to Δt[23,26,34]. By fully mapping these interactions across all local whiskers, we discovered that they generate sharp, spatiotemporally inseparable tuning for 2-whisker sequences, which represent local motion on the whisker pad. This tuning was somatotopically organized, with neurons in each column representing diverse inbound local motion vectors towards the CW of that column. This suggests a novel model of whisker coding and somatotopy in each S1 column.

Coding for space and time in multi-whisker stimuli

52% of S1 units were tuned for CW-SW combinations, with most having a single, sharp peak for one specific combination. 70% of combination-tuned units preferred a CW-containing sequence, with a strong preference for SW-leading-CW order. Individual units were also sharply tuned for Δt (10-ms resolution), but single-trial decoding of order (Δt < 0 vs. Δt > 0) was weaker than decoding of CW-SW combination identity (Fig. 3 and Supplemental Fig. 8), suggesting that the neural ensemble in each column represents time less well than space. Spatial and temporal tuning were inseparable. Thus, most S1 neurons act as spatiotemporal feature detectors, which provide a more compact representation of multi-whisker stimuli than classical single-whisker detectors. We recorded under adapted conditions similar to natural sensing, but Δt tuning could be sharper, and preferred Δt’s could be closer to 0 ms, in less adapted conditions, where stronger feedforward inhibition creates narrower integration time windows in pyramidal cells[35].

Within each column, neurons represented diverse CW-containing sequences, enabling accurate decoding of local motion sequences by firing rate in small neural populations (<50 units). The wide range of combination selectivity (Fig. 2) indicates that both single-whisker and local spatiotemporal features are represented in each column. The bias for inbound (SW leading CW) motion sequences may reduce redundancy between cortical columns in representing multi-whisker motion stimuli. The bias toward CW-SW combinations that include rostral and ventral SWs matches the prevalence of rostrocaudal contact sequences during natural whisking onto frontal objects[6], and the ventral tilt of forward whisking[36].

Tuning for 2-whisker combinations may be a fundamental aspect of S1 sensory coding that underlies tuning for more complex, higher-order combinations. Consistent with this idea, 2-whisker interactions are predictive of tuning for global motion across the whisker pad[17], and strongly predict higher order interactions for vibrotactile sequences of pairs of whiskers[34]. Responses to optimal higher-dimensional (3+ whisker) sequences are no stronger than to 2-whisker sequences[20]. Tuning for motion sequences is also highly prevalent in primate S1[37].

A revised functional role for the S1 column

These findings suggest a novel functional role for the S1 column: rather than solely representing single-whisker features, each column also represents the set of all elementary, inbound local motion sequences that involve the CW. This CW-centered, multi-whisker view of the column explains the somatotopic breadth of subthreshold input received by each cell, and is consistent with tuning to concentric motion coherence or motion contrast relative to the CW[19]. It also provides an explanation for the presence of SW-tuned neurons, which do not fit classic somatotopy defined by single-whisker tuning. Half of SW-tuned units preferred a CW-containing sequence, and thus contribute to representing local motion involving the CW. SW-tuned neurons had particularly strong combination tuning, and were essential for accurate population coding of motion sequences (Fig. 4). Since SW responses are strongest at column edges[15], we predict that combination selectivity will be greatest at column edges, where a ring of multi-whisker selectivity has been found[18]. The dominance of SW-leading-CW order selectivity predicts a functional boundary between neighboring columns based on opposite order tuning.

Sublinear sculpting of local motion tuning

Prior studies of 2-whisker integration mostly studied a single arbitrarily chosen CW-SW pair, and reported largely sublinear suppression[21,23,25,26,34,38]. Linear or supralinear integration was sometimes observed[22,24], particularly with temporally dense stimulation[20,26]. By mapping all CW-SW combinations we discovered that individual neurons suppress most CW-SW combinations, but fail to suppress a single best CW-SW combination, thus constructing sharp combination tuning. A small fraction of units (8%) showed supralinear facilitation to the best combination stimulus.

Combination tuning was separated mathematically into a global sublinear integration process, analogous to divisive normalization[39], plus combination-specific nonlinearities. This analysis revealed that combination tuning partially reflects spatial asymmetry in single-whisker receptive fields, but also involves nonlinear enhancement relative to the global sublinearity (or protection from sublinear suppression) by the CW-SW combination that contains the strongest single SW (Fig. 6). While this is not a mechanistic model, it suggests one simple neuronal implementation for sharpening of combination tuning. Sublinear integration likely reflects a combination of cross-whisker inhibition and common, sublinear summation of EPSPs in pyramidal cell dendrites[31]. Inputs representing the preferred CW-SW sequence may synapse nearby each other on dendrites, and thus preferentially recruit nonlinear boosting, including by NMDA receptors[40]. This NMDA receptor-dependent boosting mechanism also contributes to angular tuning in S1[41]. Combination tuning could be partially synthesized subcortically, but global motion tuning[42] and inter-whisker nonlinearities are generally weaker in thalamus[21,43], suggesting that cortical circuits are a likely site of computation.

A shared computation for local visual and tactile motion

The sculpting of whisker combination tuning by spatially asymmetric suppression strongly resembles the computation of visual motion direction selectivity (DS) in mammalian retina. Barlow-Levick motion detectors synthesize DS by inhibitory suppression of responses to the null direction of visual motion while the preferred direction lacks inhibition[44-46]. DS can also involve selective facilitation for the preferred direction[47-49]. Sublinear integration for non-optimal CW-SW combinations and near-linear integration for the preferred combination (with modest facilitation in a small fraction of units) is analogous to the Barlow-Levick model. Thus, there is a common computational basis for spatiotemporal feature extraction in the visual and whisker systems, although the underlying circuit mechanisms are likely to be distinct.

Methods

Surgical preparation and in vivo electrode placement

All procedures were approved by the UC Berkeley Animal Care and Use Committee and meet NIH guidelines. Male C57BL/6 mice (age: P28–45) used. Experiments 1 (Fig. 1) and 2 (Fig. 2–7) followed the same surgical and data acquisition protocols. Mice were anaesthetized with urethane and chlorproxithene (1.3 g/kg and 0.02 mg in 10 mL saline). Body temperature was maintained at 36.5C using a feedback-controlled heating pad (FHC, 40–90-8D). Anesthetic depth was assessed via toe pinch and supplemental urethane (10% of initial dose) was provided as needed. The skull was exposed, cleaned and a stainless steel head-post was implanted. A 2 mm craniotomy was made over S1 (coordinates: 1.5 mm rostral, 3.3 mm lateral of bregma). The target column (C1 or D1) was localized using receptive field mapping of multi-unit activity in L4, recorded with a tungsten microelectrode.

A silicon laminar probe (NeuroNexus, 32 channel, 1 shank, poly2 or poly3 channel geometries (A1×32-Poly2–5mm-50s-177-A16 and A1×32-Poly3–6mm-50–177-A32) was then inserted radially into the target column via a small durotomy. The probe was slowly advanced until the deepest recording pad was in L4. Simultaneous L2/3 and L4 recordings were made at this depth. Subsequently, in a subset of animals, the probe tip was further advanced to L6, and simultaneous L6 and L5a/b recordings were made at this position. L2/3, L4, L5 and L6 were defined by microdrive depths as 100–417, 418–587, 588–889 and 890–1154 μm below the pia[50].

For all recording penetrations, we confirmed post-hoc that silicon probe multi-unit activity in L4 was tuned to C1 or D1, based on single-whisker deflections. In a subset of cases, recording location was confirmed by coating the recording electrode with DiI, perfusing the mouse, and recovering DiI staining in cytochrome oxidase stained flattened tangential sections, which show the L4 barrels.

Whisker stimulation

Calibrated deflections were applied independently to a 3 × 3 grid of whiskers, centered on the columnar whisker for the recorded column (either C1 or D1). Stimuli were controlled using custom software in Igor Pro (Wavemetrics). Each whisker was trimmed to 8 mm length, and inserted into a glass tube carried on a piezoelectric bender actuator. The piezo was positioned to deflect the whisker at 5 mm distance from the face. Each whisker was deflected rostrocaudally with triphasic waveform that was shown previously to optimally drive S1 neurons and captured most of the evoked response variance in S1 (first common filter)[19] (Fig. 1a). The waveform was 40 ms duration, 300 μm peak amplitude (3.8° deflection), and had a mean frequency content of 53Hz. The waveform resembles a linear combination of position and velocity filters, which are two kinematic features that are strongly represented across the whisker-to-barrel system[19,51,52].

Experiments 1 and 2 differed in the stimulus set that was applied. Experiment 1 had a total of 261 unique stimuli; 9 single-whisker deflections and 252 2-whisker sequences (that is, 9 choose 2 spatial combinations, each of which was delivered at 7 Δt’s: 0, ±10, ±25 and ± 50ms) (Fig. 1a). Each single-whisker stimulus was delivered 125–175 times (136.67 ± 14) and each 2-whisker sequence 40–150 times (82.76 ± 35). Thus, Experiment 1 sampled all possible 2-whisker combinations within the 3×3 grid at a subset of Δt’s. All stimuli were randomly interleaved at 0.6 s inter-stimulus interval, yielding an overall average deflection rate for any whisker of 3.24 Hz. Each recording was 3.5–4.5 hours.

Experiment 2 had a total of 3000 unique stimuli; 18 single-whisker deflections (9 whiskers with peak deflection amplitude in either the rostral or caudal direction), 1600 CW-SW sequences (8 correlated and 8 anti-correlated CW-SW combinations at 100 Δt’s within ± 50 ms range) and 1400 SW-SW sequences (7 correlated and 7 anti-correlated SW-SW combinations at 100 Δt’s within ± 50 ms range). Only SW-SW combinations that involved the SW located rostral and within the same row as the CW were presented. Thus, for D1 column recordings, all SW-SW combinations involving D2 were presented, but no other SW-SW combinations were presented. This was necessary to achieve a tractable number of unique stimuli for this experiment. Each single-whisker stimulus was delivered 140–300 times (222.52 ± 44.72) and each unique 2-whisker sequence 0–50 times (8.42 ± 4.2). Thus, Experiment 2 sampled all possible CW-SW combinations, plus a subset of SW-SW combinations focused on one chosen SW, with denser and uniform sampling compared to Experiment 1. All stimuli were randomly interleaved at 0.6 s inter-stimulus interval, yielding an overall average deflection rate for any whisker of 3.12 Hz. Each recording was 4.5–6 hours in duration. In both experiments, Sham stimuli (blank trials in which no whisker was deflected) were also interleaved to quantify spontaneous spiking.

Data acquisition, preprocessing and spike sorting

Recordings were amplified and bandpass filtered (Plexon Instruments PBX2/16sp-G50, × 1,000 amplification, 0.3–8 kHz bandpass) and digitized at 31.25 kHz. Noise was reduced by common average referencing[53]. Poly 2 electrode sites were divided into groups of 4 spatially adjacent channels (tetrodes; adjacent sites were 50μm apart). Poly 3 electrode sites were divided into groups of 4 channels, selecting channels with the maximal signal-to-noise ratios located within a 50μm depth range. Negative-going spikes were detected using an amplitude threshold (2.8–3.2s.d. of noise floor), with a shadow period of 0.66ms after each threshold-crossing. Detected spikes were clipped (1.5 ms) for spike sorting.

Spike sorting used UltraMegaSort2000[54], implemented in Matlab. Clusters were excluded if they had < 1000 spikes, >0.8% refractory period violations (defined as inter-spike interval < 1.5ms), or > 30% estimated missed spikes (based on Gaussian fit of detected spike amplitudes relative to the detection threshold). For each recording site, this initial spike sorting process was used to identify the mean spike waveforms of all detectable, sortable units. To increase efficiency in spike detection, we then re-performed spike detection using a template-matching method in which spikes were identified in the raw voltage recording based on similarity in shape to these mean spike waveforms, irrespective of amplitude[55]. This method successfully detected more spikes than the initial threshold-based method, particularly spikes that had been suppressed by the shadow period. On average, 10% more spikes evoked by weak stimuli, and 40–70% more spikes evoked by strong stimuli, were found by the template matching method[55]. This method was validated using computationally modeled extracellular voltage data (data not shown). All spikes detected by the template method were then sorted using the spike-sorting process described above.

Analysis of neuronal data

Spiking responses for single-whisker stimuli were quantified as the average number of spikes that occurred within 125ms following whisker deflection onset. Responses to 2-whisker sequences used the same window, following deflection of the first whisker in the sequence. This broad window captures essentially all spikes evoked by both whiskers. Units that showed statistically significant responses relative to baseline to any of the 9 stimulated whiskers were classified as whisker responsive. Statistical significance was assessed via likelihood ratio test that assumed Poisson statistics across the 9 whisker and sham stimulation. This test is analogous to performing a one-way analysis of variance (ANOVA) across the 10 stimulus classes but assumes that spiking responses follow a Poisson rather than Gaussian-distribution. Thus, the test outputs the probability that spiking is equal across the 10 classes assuming Poisson statistics. Neurons were classified as responsive if the p-value < .05.

Δt smoothing

2-whisker sequence responses were quantified as the number of spikes that occurred within a 125ms time window after the deflection of the first whisker involved in the sequence. The large number of stimulus repetitions in Experiment 1 allowed us to compute mean sequence responses by simply averaging across stimulus repetitions (mean: 82 per Δt). Experiment 2 had a significantly lower number of stimulus repetitions (mean: 8.4 per Δt) preventing a simple averaging procedure to accurately capture neural responses. To address this we adapted a model based smoothing method originally used to smooth post-stimulus-time-histograms (PSTHs)[56]. This method assumes that 2-whisker combination responses follow an inhomogeneous Poisson process whose firing rate varied as a function of Δt.

Tuning curves were first calculated at a 1ms resolution by averaging single-trial spiking responses at each Δt. This was done in 90% of the total observations (training set) and the log Poisson likelihood of the rest of the data (test set) was calculated based on this initial, un-smoothed tuning curve. Afterwards the tuning curve was smoothed by convolving it with a Gaussian kernel with a standard deviation that ranged from 0.5 to 30ms at 0.2ms resolution. The training and test sets were selected via 10-fold cross-validation, repeated 50 times. Finally, the Gaussian kernel with a standard deviation that maximized the cross-validated log-likelihood was used to smooth the 1ms resolution tuning curve calculated on all the data. This process effectively low-pass filters the Δt tuning curve with a cutoff frequency chosen via maximum likelihood and cross-validation.

Tuning curves whose likelihood monotonically increased across all smoothing (likelihood function had no maximum) were classified as not significantly Δt modulated, and the Δt tuning curve was assumed to be constant and equal to the mean response across all Δt’s. 33% of Δt tuning curves had monotonically increasing likelihoods. Tuning curves whose non-smoothed likelihood exceeded the likelihood of all smoothed tuning curves were classified non-modulated as well. These tuning curves would imply that S1 units are sensitive to Δt at a 1ms resolution, which is not physiologically likely in S1, and a likelihood ratio test revealed that none of these tuning curves had a statistically significant fit to the data. This happened in 23% of tuning curves, hence 56% (a total of 2004 out of 3608 tuning curves) were assumed to be constant across Δt.

Combination-selectivity

For each single-unit, the 8 CW-SW smoothed Δt tuning curves in either correlated or anti-correlated directions were grouped to build two separate columnar whisker spatiotemporal receptive field (cwSTRF) that described spiking responses as a function of Δt and CW-SW combination identity. The cwSTRF’s in correlated and anti-correlated directions were analyzed separately and labeled as a unit. A combination selectivity index (CSI) for each Δt was calculated via the lifetime sparseness equation[28]: where r is the spiking response to each combination, E denotes the mean across combinations and n the number of combinations which was 8 CW-SW combinations for all units. See Ref [57] for interpretation of lifetime sparseness and its advantages over other methods of calculating stimulus selectivity. To establish statistically significant selectivity we shuffled single-trial spiking data across combinations 5000 times for each Δt, computed a non-smoothed Δt tuning curve at a 1 ms resolution and then convolved each Δt tuning curve with its corresponding Gaussian kernel for smoothing (the standard deviations were derived from non-shuffled, measured data as described in the previous section). We then built an cwSTRF from the smoothed shuffled data and computed the CSI across combination responses for each Δt. We then computed a distribution of expected CSI traces expected from shuffled data by computing the Mahalanobis distance of the 5000 CSI traces to the mean CSI trace of the shuffled data. With this distribution we computed the p-value for the Mahalanobis distance of the measured CSI trace (example shown in Fig. 2c) from the shuffled data. The threshold for statistical significance was 0.05 and controlling for multiple comparisons across units via false discovery rate control[58].

Neural decoders

Neural decoders were constructed to predict either CW-SW combination identity or order (-Δt vs. +Δt) from single-trial spike counts of recorded units. For combination identity decoders, each unit was represented by a one-versus-all (ova) classifier that was trained by logistic regression to report the probability of each combination based on the number of spikes within the evoked time window on one trial (0–125 ms after the deflection onset of the first whisker in the sequence). Training trials were selected randomly from recorded trials for that unit. For decoders based on linearly predicted spike counts, these spike counts were generated for each CW-SW combination by sampling single-trial, single-whisker CW- or SW-evoked spike counts and adding them together. Each classifier comprised eight logistic functions, one for each CW-SW whisker combination. For CW-SW order decoders, an identical process was used, except each unit was represented by a single classifier comprising a single logistic function, and spike counts were sampled from trials with −50 ≤ Δt ≤−10 ms and 10 ≤ Δt ≤ 50 ms ranges, which represent the two possible orders.

Logistic functions were fit over 1000 iterations using logistic regression and 10-fold cross-validation. Model fitting was performed using a randomly chosen subset of the recorded trials (90%), and decoder performance was assessed on the remaining trials. All sampled units in each model were treated as a single, virtual population. The population stimulus prediction was calculated by summing the probabilities of each stimulus over all units and selecting the stimulus with the maximal summed probability. This approach models the population as if all single units were independent. In Fig. 3–4, population decoding was assessed 2500 times for each ensemble size, and mean performance is reported. Different ensemble sizes were built by randomly sampling the population of model units (or ova classifiers) with replacement. For Fig. 7, a different model was fit for each Δt range, and population decoding was assessed 1000 times for each model.

Testing for statistically significant supralinearity

If a neuron’s mean spiking response to a CW-SW combination is due to linear integration then μCW-SW = μCW + μSW, where μCW-SW, μCW and μSW are the neuron’s mean spiking responses to CW-SW combination and single-whisker deflections respectively. It will also be true that , where σCW-SW, σCW and σSW are the standard errors for the respective means. Thus, we tested whether a response was statistically supralinear by testing whether μCW-SW > μCW + μSW by z-scoring the combination response relative to the linear prediction and using a one-tailed t-test to determine whether the z-score was significantly different from 0. The z-score is equal to: .

Recordings in awake behaving mice

Surgery, behavioral training, and neural recording methods were as in Ref [59]. Briefly, 4 C57BL/6 mice were implanted with a lightweight headpost at P28-P35. A week later, mice were accommodated to handling and trained operantly to lick for water reward while head-fixed on the behavioral rig. Behavioral training took place daily, and general behavior, weight, and water consumption were monitored carefully to ensure mouse health. At the start of each session, mice were transiently anesthetized (1.5% isoflurane), head-fixed via the implanted headpost, and whiskers were placed in the 3×3 piezoelectric actuator array, positioned 2.5 mm from the whisker pad. Isoflurane was halted, and behavioral trials began. Each trial lasted 11.5 sec. Mice were required to suppress licking for 1 sec to initiate the trial. In the first 10 sec of the trial, 1- and 2-whisker stimuli were delivered, randomly interleaved, every ~0.3 sec. Stimuli were identical to Experiment 1, except that amplitude was ~100 μm (2.5° deflection), and only Δt=0 and ±40 ms were presented. Only same-direction (correlated) stimuli were presented. In ~40% of trials (Go trials), a blue LED flashed (2-ms duration) at 10 sec. Mice were trained to lick in response to this visual cue. Licks during the subsequent 1.5-sec response window were detected by infrared beam break, and triggered water delivery (4–5 μL). In ~60% of trials, the visual cue was not presented, and water reward was not available (NoGo trials). After a 1-sec intertrial interval, the next trial could start. The fraction of Go trials was adjusted daily to maximize the number of trials performed. Mice were alert, did not visibly whisk, and licked only occasionally outside of Go trial responses (quantified in Fig. 8b, during recording sessions). Mice were trained until they achieved >85% lick rate on Go trials, <10% lick rate on NoGo trials, and could perform at least 250 trials per day, corresponding to >90 min of behavioral task performance. This typically took 2 weeks.

Once the mouse achieved behavioral criteria, recording sessions began. Prior to the first recording session, the C1, D1, or D2 column was localized by intrinsic signal imaging, and a craniotomy was performed. On each recording day, a laminar polytrode (NeuroNexus A1×16-Poly2–5mm-50 s-177-A16) was inserted into the targeted column in S1 while the mouse was under anesthesia in the behavior rig. Recordings were made in L2/3 and L4 during behavioral trials. Voltage signals were digitized (24.4144 kHz), and stored using a TDT RZ5D system. Spikes were isolated by off-line bandpass filtering (0.3–6 kHz) and common average referencing, using custom MATLAB code. Spike detection and sorting was performed as for anesthetized recordings. Epochs within ±200 ms of a lick were excluded, to avoid contamination by lick-related movements and lick-related spiking. Data were analyzed from recordings whose L4 multi-unit tuning clearly matched the target whisker column, defined as mean CW-evoked spiking > 1 SEM more than any other whisker. 6 recording sessions in 4 mice passed this criterion and were analyzed. After recording, the craniotomy was resealed to allow another recording session 1–2 days later.

Statistical Methods and Experimental Design

Information on experimental design and statistics, including sample size, replication, randomization, blinding, and exclusions, is in the Life Sciences Reporting Summary. Statistical tests were implemented in Matlab. Standard error (SE) of sample proportion, used in some figure panels, was calculated as SE = sqrt( (p*(1-p)/N ), where p = proportion of units with a specific binary characteristic and N = number of units. No statistical methods were used to pre-determine sample sizes, but our sample sizes are similar to those reported in previous publications[19,24,27,33]. In cases where parametric statistics are reported, data distribution was assumed to be normal, but this was not formally tested. Randomization is not relevant for this study because there are no animal treatment groups. Spike sorting was performed blind to any stimulus or location-related information, which ensures that unit selection is not biased to specific tuning properties.

Data were excluded from recording sessions in which the recording electrode was not positioned in the target cortical columns for the study. In Experiments 1 and 2, one session (i.e., one mouse) out of 18 was excluded for this reason. During spike sorting, single units with less than 1000 total spikes were excluded from analysis, which is necessary to accurately estimate tuning properties. All other units are included in the study. Experiment 3 used the same analysis and exclusion criteria, however because targeting was more difficult in the awake experiments, six sessions were excluded from analysis due to the electrode not being in the target cortical column.

Data availability

Data are available in the Collaborative Research in Computational Neuroscience repository, CRCNS.org.

Code availability

The Matlab code for performing the statistical and data analysis are available upon request to the