Current theories of cortical sensory coding are mostly rooted in the concepts that were established by Hubel and Wiesel (1) in the primary visual cortex, where single neurons were found to be selectively activated by specific stimuli such as oriented bars of light. Following these discoveries, understanding the feature lexicon and how it can be used to build elaborate internal representations has been one of the main aims of sensory neuroscience. This conceptual framework has been central in the whisker/barrel system, where early physiological studies in anesthetized rodents examined what whisker movement features were encoded by artificially varying one stimulus parameter at a time. Velocity (2), displacement amplitude (3), direction (4), deflection duration (3), whisker identity (5), and vibration frequency (6) were just some of the parameters found to be encoded by barrel cortex neurons. This level of feature diversity was in line with what had been found in the visual system and delineated a rich coding basis for building elaborate tactile representations.
More recently, improved stimulus control (7, 8) across the entire mystacial pad enabled reverse correlation techniques that were developed in the visual and auditory fields (9, 10) to be adapted to the whisker system. These approaches are geared toward anesthetized rodents, and rather than exploring the stimulus space one input parameter at a time, they are based on broadband stimulation of the whiskers and analysis of the spike-eliciting stimulus features in small prespike temporal windows. In barrel cortex, these techniques have identified a low-dimensional coding space containing linear combinations of only two basic stimulus features (11–13). These fundamental features were extracted by applying dimensionality reduction methods to the spike-triggering features observed across large populations of single neurons recorded in granular and infragranular layers of barrel cortex (12, 13). Despite the power of this approach, the derived features are oscillatory whisker movements that are difficult to relate to simple kinematic or mechanical parameters (i.e., direction, velocity, and curvature encoding), and the low dimensionality is at odds with the rich coding diversity observed when varying one stimulus parameter at a time.
Along with the conflicting results from different artificial stimulation approaches, a major obstacle in establishing a universal computational framework in the whisker system has been the difficulty in reconciling coding schemes observed under artificial conditions with those found in behaving animals. Under naturalistic conditions, high-speed videography and precise whisker tracking (14–16) have identified whisker curvature (17, 18), angular velocity (18, 19), and angular acceleration (20) as the whisker movement parameters most correlated with neural activity. High-velocity, high-acceleration events called stick-slips (21) have received notable attention because of their correlation with barrel cortex activity in rodents trained to whisk into surface textures (19) and discriminate them based on roughness (22). During other perceptual processes such as whisker-based object localization, barrel cortex neurons encode multifaceted features such as the phase of the whisking cycle in which object contact occurs, independent of the amplitude and set point of the ongoing whisker search behavior (23, 24). This suggests that barrel cortex cells can integrate inputs related not only to whisker-object contact forces but also to the state variables governing ongoing behavior (whisking set point, amplitude, and frequency). Subsequent “active sensing” studies examining natural whisker search (25) or electrically induced whisking (26) have found information about these state variables in various forms along the somatosensory neuraxis from the mechanoreceptors to the primary somatosensory and motor cortical areas (23, 25–27). Whether phase selectivity can be explained by a classical feed-forward feature-based coding framework remains unclear, but the presence of phase tuning in whisker follicle mechanoreceptors suggests that it can (26). While these naturalistic studies have obvious advantages, they are usually limited to one or a few highly stereotyped conditions, and natural behaviors can vary greatly in different sensory environments (28). These variations can make it difficult to generalize coding schemes outside of the limited conditions that are examined, and clarification of their implications is often necessary in precisely designed artificial contexts (29).
There is an immediate need for a unified model in the whisker system to tie together this constellation of schemes and form a general theory of whisker-based sensory processing that can help to carve out universal sensory coding principles. Such a model should draw on coding schemes uncovered by artificial stimulation approaches and natural behaviors requiring active or passive sensation and, at the same time, be able to explain a range of whisker-based perceptual processes from object localization to texture discrimination. To move in this direction, we begin by showing that the oscillatory features found in reverse correlation studies of barrel cortex (11–13) are compatible with a simple stick-slip model (based on naturalistic studies) in which neurons act as detectors of isolated high-velocity events. To investigate the diversity of stick-slip feature encoding, we designed two different stimulation conditions, one that has a flat power spectrum in the velocity domain and one that is spectrally flat in the acceleration domain. We find that these input spaces more than double the number of cells that show a preference for particular whisker motion features compared to what was found in prior studies with a stimulus that was spectrally flattened in the angular position domain (POSW). Beyond increasing the number of functionally selective granular and infragranular barrel cortex neurons, we disambiguate two major feature classes, sticks and sweeps. Sticks correspond to contact onsets and are the high-velocity, high-acceleration events that have been studied under natural conditions (19). Sweeps are a new functional class that is composed of high-velocity events that induce large angular displacements (>3°). These functional regimes mix in single deep-layer barrel cortex cells, where neurons can express one type of encoding or both, with or without direction selectivity. This functional integration in deep layers of barrel cortex is in line with what has been found in primary visual cortex, where parallel information streams that are segregated at the periphery come together in single deep-layer cortical neurons (30). These results establish a unified computational framework drawing on many different experimental approaches and delineate the two dominant feature classes encoded in deep layers of barrel cortex.
A simple stick-slip model produces the same feature spaces as reverse correlation techniques
One of the ways that feature encoding has been studied in the whisker system is using reverse correlation techniques (Fig. 1A, top). In general, the procedure for this approach is to (i) design a broadband stimulus that samples a large number of stimulus configurations, (ii) record the spiking activity of a neuron during presentation of this stimulus, and (iii) tabulate the spike-eliciting stimulus events to extract the features that are encoded using methods such as spike-triggered averaging (STA) or spike-triggered covariance (STC) (Fig. 1A, top). All studies using this approach in the barrel cortex have selected a stimulus that is POSW (11–13). These studies used the STA/STC methodology to show that the features encoded in barrel cortex are oscillatory whisker movements (see the two examples in Fig. 1A, right) (11–13), while studies under natural conditions suggest that barrel cortex neurons detect simple high angular velocity and acceleration events (19, 31). To test whether the oscillatory features that are obtained from reverse correlation approaches are in line with the results from naturalistic studies, we developed a simple stick-slip model. Our model hypothesizes that single barrel cortex neurons respond to isolated high-velocity or high-acceleration whisker movements at low latencies with some jitter (fig. S1A). The time window in which high-velocity, high-acceleration movements can affect firing is modeled as a skewed Gaussian (latency τ, jitter σ, and skewness η), and the respective contributions of absolute velocity and acceleration are weighted by the α parameter (Fig. 1A, bottom, and fig. S1A).
To determine biologically plausible parameter ranges for our model, we used data from a previous study (12) that extracellularly recorded spiking activity from a large population of granular and infragranular barrel cortex neurons during multiwhisker POSW stimulation in anesthetized rats. We explored a large range of values for each parameter (α, τ, σ, and η were varied systematically; see fig. S1C and Materials and Methods) and found the stick-slip model (combination of four parameters) that yielded the filter/filters most similar to the real filter/filters for each recorded cell (error and parameter distributions in Fig. 1B; red histograms, n = 560 cells). The error for each model filter set was computed by directly projecting the model filters into the space defined by the real filters. If this projection is close to 1, then the error is close to 0, and the filters are the same shape (fig. S1C). With the resulting biologically plausible parameter distributions, we simulated a population of 500 model cells respecting the parameter covariances observed with our optimization routine (Fig. 1B, blue histograms, and fig. S1B). If we compare the features that come from the real cells with those that come from the simulations (real cells and simulation; Fig. 1, C and D, respectively), then the resulting feature spaces are nearly identical. Principal components analysis of both feature spaces uncovers two dimensions that explain almost all the variance in the feature shapes. The two dominant principal components are oscillatory, and projecting the real feature principal components into the space defined by the simulations shows a high correspondence (red features in Fig. 1D are close to the unit circle). This shows that a population of simple stick-slip detectors, as defined by our model, is capable of generating the oscillatory feature spaces found in reverse correlation studies from infragranular and granular barrel cortex neurons.
Spectral equalization in the velocity domain improves sampling of stick-slip–like features
As studies under both naturalistic (19) and artificial broadband stimulation conditions (Fig. 1) suggest that barrel cortex neurons are mainly sensitive to isolated high angular velocity events, we next considered whether the POSW stimulus typically used for reverse correlation studies efficiently samples the input spaces preferred by barrel cortex neurons. An incomplete or inadequate set of input configurations could (i) preclude the identification of some of the features to which the neurons are tuned or (ii) mask key differences between encoded velocity features (make all neurons look like they encode the same feature). Looking at how high angular velocity events occur in a stimulus that is constrained to be POSW, it is clear that both of these concerns are warranted (Fig. 2A). First, the power spectrum of a POSW stimulus lacks low-frequency content in the velocity domain (Fig. 2A, bottom left inset) because, according to mathematical definitions, the first derivative of a signal with a flat power spectrum cannot also have a flat spectrum. To illustrate, if just the highest angular velocity whisker movements are selected (Fig. 2A, center and right), the high frequencies are prominent in the resulting short stimulus windows (Fig. 2A, right insets). The average high angular velocity events that occur under POSW conditions and their temporal correlations (Fig. 2A, right) suggest that these events are of limited variability and not well isolated from other high angular velocity movements (Fig. 2A, bottom right, many adjacent peaks). This implies that if neuronal spiking has some temporal imprecision, as expected at the level of the cortex, then it will not be possible to ascertain which high-velocity event triggered the spike.
One solution to this problem is to generate a stimulus that has a flat power spectrum in the angular velocity domain (VELW) (Fig. 2B). This involves adding low-frequency content in the angular position domain. (Fig. 2B, left insets, and fig. S2), which increases baseline drift and more closely resembles natural whisker movements under certain conditions (19). The main advantage of this VELW stimulus condition is that the high angular velocity events are well isolated from other high angular velocity events (Fig. 2B, bottom right, single peak with minimal side bands). This does come at a cost in the angular position domain, where high angular velocity events are now associated with long-range correlations in angular position (Fig. 2B, top right). This seems preferable to the short-range correlations in high angular velocity events present under the POSW condition if we consider that both naturalistic and artificial stimulation approaches indicate that the encoded features are in the velocity domain. Because some studies pointed toward angular acceleration as an important domain (20), we also designed a stimulus that has a flat power spectrum in the angular acceleration domain (ACCW), but this stimulation condition requires strong adjustments to keep the baseline drift within the operating range of our piezoelectric devices. These improved stimulus conditions should more efficiently sample the input spaces that are likely to be preferred by barrel cortex neurons.
Spectral equalization in the velocity or acceleration domain boosts detection of tactile feature encoding cells
To test the spectrally flattened velocity and acceleration stimuli, we performed extracellular recordings (244 well-isolated single units; Fig. 3A) with multielectrode silicon probes placed in granular and infragranular layers of barrel cortex (inserted to ~1400- to 1500-μm depth spanning 340 μm with four shanks or 1700-μm depth spanning 1.5 mm with a single shank) during multiwhisker stimulation. During the recordings, light isoflurane anesthesia was administered and carefully monitored (see Materials and Methods) to maintain functional responses in the cortex as has been described in prior studies (12). Identical stimuli (POSW, VELW, or ACCW) were presented to all 24 caudal macrovibrissae simultaneously (11–13).
Spikes were recorded from large populations of single units and sorted into well-isolated clusters using the Klusta suite (see Materials and Methods). After spike sorting, significant STA and STC filters were extracted for each unit under all three stimulation conditions (Fig. 3B, example cell). These significant filters represent the subspace of the stimulus distribution that is associated with the spiking activity of each cell. In line with prior barrel cortex studies (11–13), the STC filters recovered under the POSW stimulation condition had an oscillatory character (Fig. 3B, multiple lobes around rest position). In contrast, the significant STC filters from the VELW stimulation condition, now defined in the velocity domain, are simpler. They are either reduced to a single bump, indicating sensitivity to high velocity (fig. S4), or to a single oscillation, suggesting sensitivity to acceleration or position depending on the characteristic frequency of the oscillation (Fig. 3B and fig. S4B). For the ACCW stimulation condition, the first STC filter, shown in the acceleration domain, is a single oscillation indicating a sensitivity to high velocity, and the second STC filter is a uniphasic bump with small side bands, which implies sensitivity to acceleration (Fig. 3B and fig. S4C).
Looking at total filter yields and general spiking activity across the population, there are marked differences between the three stimulus conditions. As expected, the VELW and ACCW stimuli uncover much larger numbers of significant filters [79, 232, and 202 significant STC filters for POSW, VELW, and ACCW, respectively (Fig. 3C), and 56, 120, and 86 significant STA filters for POSW, VELW, and ACCW, respectively (fig. S4, D to F)] and many more functionally responsive neurons (43, 112, and 108 neurons for POSW, VELW, and ACCW, respectively; Fig. 3D) and elicit slightly but significantly higher firing rates (Fig. 3E) than the POSW stimulus. In terms of spiking reliability, the trial-by-trial cross-correlations for repeated white noise stimulations are also much higher for the VELW and ACCW stimulation conditions than for the POSW stimulation condition at the single neuron and population levels (Fig. 3F and fig. S3B). This suggests that the VELW and ACCW stimulation conditions more reliably activate the feature detectors that exist in barrel cortex.
To examine the structure of the STC filters coming from the entire population of neurons, we used principal components analysis (Fig. 3, G to I). For the POSW condition (Fig. 3G), we observed the same multiphasic feature space as past studies (12, 13), with two dominant dimensions (representing the basic features) that explain a similar amount of variance (Fig. 3G, see eigenvalue spectrum, more than 90%). When projected into this two-dimensional space, most of the individual features have similar weights (they lie near the unit circle) and are uniformly spread with no observable pattern. In contrast, for the VELW condition (Fig. 3H), the two dominant principal component dimensions explain very different amounts of variance (but still more than 80% of the variance between the two of them), and the nonuniform distribution of individual features in the space contains clusters and hotspots that could represent functional classes. For the ACCW condition (Fig. 3I), the results are more similar to the POSW condition with two dominant dimensions that have similar levels of explained variance.
Spectral equalization in the velocity domain unmasks feature selectivity
STC filters capture the subspace of a stimulus distribution in which a neuron is most responsive, but it is necessary to look where a cell’s individual spike-eliciting events fall within these filter dimensions to assess the selectivity of feature encoding. Under all three stimulus conditions, the feature spaces are well described in the two dimensions defined by the top two principal components of the population principal components analysis. Therefore, we used these principal component dimensions as common feature spaces and examined every cell’s spiking-eliciting stimulus events in these common spaces. This analysis reveals some cells with elevated firing rates in specific subregions (or a specific subregion) of the common feature spaces for all three stimulus conditions. For VELW and ACCW, example cell 1 (Fig. 4A, cell 1) has elevated firing behavior that is oriented along the most uniphasic parts of the feature spaces, which corresponds to high-velocity movements for VELW conditions and high-acceleration movements under ACCW conditions (fig. S4, B and C). The symmetry of example cell 1 in the VELW feature space indicates that it has elevated firing for high angular velocity movements in either direction (caudal or rostral), whereas in the ACCW feature space, it is more tuned to caudal acceleration. For the POSW condition, example cell 1 has slightly elevated firing rates for most of the perimeter of the space with one hotspot. The oscillatory feature associated with this hotspot is hard to interpret in terms of directional velocity and acceleration features. There are many cells that show no regional selectivity in the POSW feature space but that are selective in either the VELW or the ACCW feature space. Example cell 2 (Fig. 4A cell 2, same cell as Fig. 3B) is selective for high-velocity movements in the rostral direction in the VELW feature space and more symmetric but caudally biased in the ACCW space. Last, there are also some cells (Fig. 4A, cell 3) that are not selective for any region of the three feature spaces but do have elevated firing rates for all stimuli that fall near the unit circles.
If we look at the total subregional selectivity across the entire population of neurons in all three spaces, then we find markedly different amounts of feature selectivity under the three stimulus conditions. On the basis of the same stringent statistical tests used in past studies (12, 13), only 16% (26 of 162) of the cells that gave significant filters to any single stimulus type were selective to subregions of the POSW feature space, while 73% (118 of 162) of the same population of cells were selective to at least one subregion for the VELW feature space, and 41% were selective within the ACCW feature space (Fig. 4B). These differences, especially between the POSW and VELW conditions, certify that the VELW stimulation not only uncovers more feature encoding cells but also does better at disambiguating different encoded features in the velocity domain.
Definition of sweep- and stick-encoding functional classes
The nonuniform distribution of STC filters in the VELW feature space (Fig. 3H, bottom) and the large number of cells (n = 118 cells) that show regional selectivity (Fig. 4B) suggest that there might be functional classes present in the population. Because these effects were most pronounced during VELW stimulation, which could partially be due to the physical constraints that are imposed on the ACCW stimulus by our piezoelectric devices (low frequencies cut to keep drift within dynamic range of piezo; fig. S2D), we started by investigating the VELW feature space. We defined the stick axis in the VELW feature space (Fig. 5A, coral dashed line) as the dimension along which the features are uniphasic (one-sided bumps). The axis orthogonal to the stick axis (Fig. 5A, black dashed line), which we define as the sweep axis, is the dimension that corresponds to features that are a mix between sensitivity for velocity and position (fig. S4B, dashed filter). To quantify the different types of feature selectivity present in the population, we determined each cell’s orientation angle θ (−π/2 < θ < π/2) by projecting all spike-eliciting stimuli into the VELW feature space, computing the vector sum of these projections, converting the sum into polar coordinates (r, θ), and then multiplying the angle by 2 (see Materials and Methods). Because the feature space is symmetric (see coral and black feature shapes, which are inverted on opposing sides; Fig. 5A), we decided to examine all cells based on orientation angle alone first and treat the directionality separately. Across the population of cells, the resulting probability distribution of orientation angles is not uniform (Fig. 5B). The hotspots could correspond to functional classes that are not perfectly separated using this type of analysis.
To study how these putative classes might reveal their differences under other stimulation conditions, we considered each cell’s response to a sparse-noise stimulus. In sparse-noise stimulation, each whisker is moved independently (in a random order) and receives a ramp-hold-ramp movement trajectory where the velocity is constant during the movement periods and zero at all other times (Fig. 5C and fig. S2). This type of stimulation allows a fine temporal analysis of neural responses with respect to the onset of the stimulation of single whiskers (Fig. 5C, top right) and also a spatial assessment of which whiskers fall within the full receptive field (RF) of the neuron (Fig. 5C, bottom right). When examining the sparse-noise responses, we noticed that cells with elevated firing rates along the stick axis (coral axis) in the VELW feature space also respond with low-latency and low-jitter (sharply) to sparse-noise stimulation of at least one whisker. To quantify this, we detected sharp single-whisker sparse-noise responses by looking for high signal-to-noise ratio (SNR > 5 * basal firing rate) and low-latency (<15 ms after stimulus onset) responses across the population. We found 30 cells (25% of cells tested) with single-whisker sparse-noise responses that passed these criteria, where each cell had at least one sharp single-whisker response and as many as three. Example “stick-encoding” cells (Fig. 5D) are tuned to orientation angles that are near the stick axis in the VELW feature space (Fig. 5B, coral arrowheads), can be direction selective, and have sharp principal whisker (PW) sparse-noise responses [Fig. 5D, PW PSTHs (peristimulus time histograms)]. Notably, the cells do not respond as strongly (and sometimes not at all) to the return part of the sparse-noise ramp stimulus that moves the whisker back to its resting position. This is likely because the 10-ms hold time for our sparse-noise ramps is within the inhibitory period in barrel cortex that comes after stimulation (32). We quantified the sharpness (PWsh) of all PW stick-type responses as the peak of its PSTH divided by the half-width (Fig. 5C, top right, and fig. S5A).
Next, we observed that cells with an orientation angle in the VELW feature space that are pulled toward the sweep axis (black axis) have quite different responses to sparse-noise stimulation in terms of both temporal dynamics and whisker distribution. Example sweep-encoding cells (Fig. 5E) show selectivity for orientation angles that are biased toward the sweep axis in the VELW feature space (Fig. 5B, black arrowheads), can be direction selective, and have long-latency and high-jitter responses to sparse-noise stimulation (Fig. 5E, PW PSTHs). The sparse-noise responses are also more equally distributed across the whisker pad (Fig. 5E, full RFs). We quantified this RF spread (RFsp) (Fig. 5C, bottom right) in the population of cells by summing the Z-scored single-whisker responses (mean firing rate in 55-ms window after sparse-noise stimulation) and multiplying by a sharing parameter (maximal when response is shared equally across all whiskers). Note that if a sharp response was already detected for sparse-noise stimulation of a whisker, then this whisker was removed from the computation. Using this RFsp (Fig. 5G, bottom, and fig. S5B) and a shuffling approach (see Materials and Methods), we detected 84 sweep-encoding cells (71% of cells tested).
From the 30 stick-encoding and 84 sweep-encoding neurons, there are 20 cells (17% of cells tested) that passed both sets of criteria. These sweep-stick cells (Fig. 5F) generally had orientation angles in between the stick and sweep axes in the VELW feature space (Fig. 5B, light blue arrowheads), responded sharply to least one whisker during sparse-noise stimulation (Fig. 5F, PW PSTHs), and had extended multiwhisker RFs with less sharp temporal response properties (Fig. 5F, full RFs). One interesting observation is that the sweep-encoding component tended to span whiskers that are rostral on the whisker pad with respect to the stick-encoding component (Fig. 5F, full RFs, and fig. S5C).
To visualize the orientation angles across these three types of cells (stick, sweep, and sweep-stick), we normalized the PWsh and RFsp for each neuron and computed a scaled classification index (see Materials and Methods). When this index is plotted against the orientation angle in the VELW feature space, it reveals that exclusive stick-encoding neurons (Fig. 5G, 10 coral dots) are centered on the stick axis in the VELW feature space, exclusive sweep-encoding neurons (Fig. 5G, 64 black dots) are pulled toward the sweep axis, and mixed sweep-stick neurons (Fig. 5G, 20 light blue dots) mostly fall in the region in between. Thus, the spatiotemporal dynamics of the sparse-noise responses correlate with what feature a cell is selective for in the VELW feature space. These data together suggest that the distribution of orientation angles (Fig. 5B) reflects at least two overlapping functional classes, where single cells can belong to either class or both. It is also possible to see this functional separation using only the temporal dynamics of the sparse-noise response to the whisker that evokes the strongest response (most spikes) by computing a ratio of low-latency to high-latency spikes (Fig. 5H, w1:w2 ratio). Stick-encoding neurons have higher ratios than sweep-encoding neurons, and the mixed cells fall at intermediate values (Fig. 5I). With these classifications in place, we looked at the direction selectivity in a class-specific manner. The stick-encoding population (10 exclusive stick-encoding cells) contained cells with both rostral and caudal directional tuning when measured either in the VELW feature space or with sparse noise (Fig. 5J, left and middle), as did the sweep-encoding population (Fig. 5J, right). This directional tuning is in line with early barrel cortex studies (3) and was not recognizable in reverse correlation studies that only used POSW stimulation (11–13).
To understand how our newly defined cell classes respond to all of the different stimulus sets, we reexamined the spiking activity during POSW and ACCW types of white noise. If a cell encoded sticks (either stick cell or sweep-stick cell), then it always gave significant filters for all three stimulation types (POSW, VELW, and ACCW). However, from the cells that exclusively encode sweeps, less than half of them (30 of 64) had significant responses to POSW stimulation, while almost all of them yielded significant filters to either VELW (62 of 64) or ACCW (56 of 64) stimulation. This suggests that the high angular velocity movements under the VELW and ACCW conditions that are coupled with larger angular displacements (fig. S2) more strongly activate the sweep-encoding cells than the high-velocity movements under the POSW condition that are coupled with smaller angular displacements. The velocity ranges are matched for all three white noise stimulation types (See Materials and Methods). Last, we examined the cortical layer positions of every recorded neuron using current source density (CSD)–based layer identification (33). All neurons that encoded stick information were located in cortical layers known to receive lemniscal pathway inputs (layers 4 and 5b/6a; fig. S6). Sweep-encoding cells were spread across all of the recorded layers (layers 4, 5a, 5b, and 6a; fig. S6). This suggests that stick-related information propagates through the lemniscal pathway, while sweep information arrives through other sources.
Specialized sweep- and stick-encoding models
To generalize the conceptual differences between stick features and sweep features, we returned to our simple stick-slip model (Fig. 1). We examined what modifications are needed (if any) so that the model could generate the VELW feature shapes for both functional classes. These models are important because they clarify exactly what the observed filter shapes mean in physical terms outside of the stimulus distribution under our experimental conditions, which is not always intuitive.
Looking first at the stick-encoding functional class, we noticed that both the STA and the STC filters (Fig. 6A, right, and fig. S4E) from some of these cells contained two time scales on which they were sensitive to high-velocity events and these two time scales usually had the same directional preference (a double-bump one-sided filter). The simple stick-slip model (Fig. 1 and fig. S1) cannot generate filters similar to this because it only has one skewed Gaussian temporal window, which can never create a double peak in the same direction without also having a peak in the opposite direction in between (see fig. S4, A to C). To rectify this, we added a second skewed Gaussian window to the model and a directional velocity term, which introduces four new parameters (Fig. 6A). This stick-encoding model performed much better at generating the stick-filter shapes than the simple stick-slip model (coral versus gray histograms; Fig. 6C, top left) and also better than models that either only added a directional velocity term or a double window (fig. S7, A and C). The best fitting model parameters usually have both directional and absolute velocity components (Fig. 6C, top middle and top right). For the temporal aspects, the parameter distributions suggest that the inter–time scale intervals are ~10 ms (difference between τ1 and τ2; Fig. 6C, bottom left) and that the longer-latency window is less sharp (σ1 < σ2; Fig. 6C, bottom center). These double window filters could reflect input from adjacent whiskers that can facilitate the response to PW stimulation if the interwhisker stimulation delay is in the appropriate time window.
For the sweep-encoding population, the STC filters have a different character (Fig. 6B, right). Within the VELW stimulation conditions, the shapes of these filters suggest that the high-velocity movements of the whiskers need to induce large angular displacements (fig. S4B, dashed filter). To accommodate this, we added two terms to the simple stick-slip model (Fig. 1 and fig. S1), an absolute position term and a position-velocity interaction term (Fig. 6B). These two terms greatly improved the model’s ability to generate sweep filters (Fig. 6D, top left, black versus gray histograms), and a model with either term alone did not perform nearly as well (performance was also worse for models with acceleration terms) (fig. S7, B and D). As expected from the sparse-noise responses, sweep-encoding cells have long latencies (τ; Fig. 6D, bottom left) and large levels of spike jitter (σ; Fig. 6D, bottom center). Their windows are also always skewed left (η < 0; Fig. 6D, bottom right).
Note that whereas the simple stick-slip model (Fig. 1 and fig. S1) was sufficient to explain the feature shapes obtained during POSW stimulation, it provided very poor estimates for the VELW feature shapes, particularly for the sweep-encoding domain. The most plausible explanation for this is that from a feature space derived from POSW conditions; the feature encoding is blurred to a degree where a single model without functional domains can suffice. When the stimulation conditions change to allow the disambiguation of functional domains, new models are required that are domain specific (and perhaps sensory pathway specific). On the basis of the understanding derived from these two specialized models, we interpret the two functional classes to code for two different kinds of whisker movements. Stick encoding pertains to movement onsets (likely touch onsets; Fig. 6E). Stick detectors are activated at low latencies by high-velocity, high-acceleration events from a precise location on the whisker pad (~1 to 3 whiskers with usually one dominant PW). The sweep-encoding regime also encodes high-velocity movements, but rather than being activated at the onset of a movement from a single whisker, it is sensitive to the extent of the angular trajectory (Fig. 6F). That is, the high-velocity movements need to create large displacements (>3°) across the pad to activate the sweep detectors. These two functional regimes are integrated in deep layers of barrel cortex.
Our main goal in this study was to bridge the gap between the artificial stimulation tradition in barrel cortex research and the recent work in natural contexts. We began by clarifying the correspondence between the oscillatory feature shapes from past reverse correlation studies (11–13) and the stick-slip events (Fig. 1) that correlate with barrel cortex activity in animals trained to whisk into surface textures (19). Considering this correspondence, we noted that a POSW stimulus is not an ideal input space to sample the velocity-based features that correspond to stick-slip events (Fig. 2). To solve this problem, we designed stimuli that were spectrally equalized in the angular velocity and acceleration domains. At a glance, spectrally equalizing the velocity and acceleration domains of the stimulus results in whisker movement trajectories that are much more similar to natural whisker movement patterns (26) than the classically used POSW stimulus. These extended stimulus sets uncover at least twice as many functionally responsive cells and engender a more structured feature space than previously observed with reverse correlation methods in deep layers of barrel cortex (Figs. 3 and 4).
Our improved stimulus sets disambiguate two feature classes, which we define as sweeps and sticks (Fig. 5). Stick encoders emit precisely timed spikes at low latencies in response to high-velocity, high-acceleration whisker movements. They can be direction-selective (Fig. 5) and often integrate two time scales of high-velocity whisker movements with a latency difference of ~10 ms (Fig. 6). This time difference is in line with facilitation windows observed in multiwhisker integration studies in barrel cortex involving pairwise deflections of adjacent whiskers (34). All of these properties are consistent with stick encoding reflecting lemniscal pathway input. In support of this conclusion, we observed that all stick-encoding cells reside in cortical layers known to receive direct lemniscal inputs (fig. S6) (35, 36). The sweep encoders exhibit quite different characteristics. While sweep features also contain high-velocity movements, these movements must also result in large angular displacements (>3°). Sweep encoding is predominant in the barrel cortex population (71% of the classified cells encoded sweeps) and can also be direction-selective. The multiwhisker character of the RFs for sweep encoders and the long-latency, high-jitter responses to stimulation of single whiskers resemble the responses found in the posterior medial nucleus of the thalamus (37), which is the thalamic nucleus associated with the paralemniscal pathway. Together, these data suggest that cortical cells displaying mixed characteristics, the sweep-stick hybrid cells, might be integrating information from the different well-documented sensory pathways of the whisker system. Such an integrative process has also been documented in deep layers of primary visual cortex (30). While stick encoding is easy to reconcile with the “stick-slip” coding schemes that came from texture discrimination studies, sweep-stick coding can also be applied to other well-known perceptual processes, such as pole localization. From studies examining either electrically induced whisking (38) or active sensing in behaving rodents (23), whisking phase–selective encoding of whisker touches can be explained by combining directional sweep and stick information in single cells. For a given whisking set point, the protraction phase of the whisking cycle can be thought of as a rostral sweep, and the retraction phase can be thought of as a caudal sweep. If a sweep-stick cell is sensitive to rostral sweeps and caudal sticks, then it would selectively fire for contacts in the protraction phase. Therefore, the combination of these features can support multiple perceptual processes.
Our work resolves what has been more than a decade of uncertainty about how to interpret reverse correlation studies in the whisker system. The success of these methods in the functional assessment of the early visual and auditory systems led to their widespread use in other sensory systems (olfaction and touch), where they have had much less of an impact. This is likely due to nuances of each system that demand small alterations in the core methodology to sample the appropriate stimulus dimensions. It is then necessary to find adapted analysis methods to disentangle the functional implications of these alterations. In the whisker case, the broadband stimuli are generated by devices with different physical constraints than the ones used for vision or audition. These limitations force the results to depart from the underlying theory, which assumes ideal conditions (i.e., band-unlimited white noise), in misleading ways. We have now aligned the reverse correlation traditions in barrel cortex with early neurophysiological studies that found selectivity for simple stimulus parameters, as well as with coding schemes developed under natural conditions. Our results define the two dominant functional regimes found in the deep layers of barrel cortex, stick encoding, and sweep encoding. The identification of these functional regimes is an important step toward constructing a general computational theory for tactile information processing in barrel cortex.
MATERIALS AND METHODS
All experiments were performed in accordance with the French Ethical Committee (Direction générale de la recherche et de l’innovation) and European legislation (2010/63/EU). Procedures were approved by the French Ministry of Education in Research after consultation with the ethical committee #59 (authorization number 2015060516116339). The recordings were performed in three male Wistar rats (260 ± 16 g). The rats were first anesthetized with 3% isoflurane that was mixed with 80% N2O and 20% O2 and delivered at 1 liter/min. Once deeply anesthetized, animals were placed into a stereotactic device with ear bars and a nose clamp. The body temperature was monitored with a rectal probe and maintained at 37°C. The eyes were coated with an ophtalmic gel (Opthalon) to keep them from drying out. Lidocaine (1%) was injected below the skin 10 min before a straight incision was made along the sagittal suture of the skull. The incision was stretched and held with surgical clamps. The bone was cleaned, and a head fixation post was cemented to the right parietal bone. A small craniotomy (~0.5 mm) was then made over the left frontal lobe, and a small electrocorticography (ECoG) electrode was placed under the skull in contact with the surface of the brain. Through this electrode, the depth of the anesthesia was continually assessed throughout the experiment. A second craniotomy was performed over the barrel cortex (~5.7 mm lateral and 3.7 mm posterior to bregma) that was ~2 mm in diameter. After performing a durotomy, a four-shank (or single-shank) silicon probe was inserted at an angle of 58° with respect to the surface of the brain to a depth of ~1.4 to 1.5 mm. Once the electrode was in position, the anesthesia was lightened by gradually decreasing the isoflurane concentration until the ECoG was stable in stage III, plane 1 to 2. Throughout the experiment, the physiological state of the animal was controlled such that the respiration rate stayed in the range of 1 to 1.5 Hz, there was a lack of motor activity in both the paws and the eyes, and the ECoG showed fast oscillations (>5 Hz), as has been shown to maintain functional responses in the cortex (12).
In all experiments, 1,1′-dioctadecyl-3,3,3′,3′-tetramethylindocarbocyanine perchlorate (DiI) was deposited on a single-shank electrode, and the barrel for the delta whisker (identified electrophysiologically with single-whisker deflections) was marked with a single penetration before the recordings. Note that this is not always the same electrode on which the recordings were made. After the DiI marking, four-shank electrodes were lowered for the recordings targeting the D1 barrel and its neighbors. These shanks, although not DiI-labeled, are visible next to the DiI-marked barrel (fig. S8). When the recordings were finished, a lethal dose of pentobarbital was injected into the peritoneal cavity. The rat was perfused transcardially with saline and then 4% formaldehyde solution after clamping the descending aorta. The brain was dissected and then stored overnight in 4% formaldehyde solution at 4°C. The cortex was tangentially cut into slices with 80 μm in thickness to reveal the barrel map, which was stained with cytochrome oxidase.
Raw electrophysiological traces were recorded at 30 kHz with a Blackrock acquisition system. They were then processed using the Klusta suite (39). Briefly, the raw signals were filtered between 500 Hz and 3 kHz with a third-order Butterworth filter. The Spikedetekt2 and KlustaKwik2 extraction and sorting algorithms were then used to extract and cluster the identified spikes. Single units were then manually curated by looking at autocorrelograms (>1-ms refractory period) and cross-correlograms and making sure that there was no drift throughout the experiment. After the manual curation, signal-to-noise-ratio (SNR) was evaluated for each triggered spike with respect to root mean square noise of the recorded traces. Spikes with an SNR of <3 on all recording sites were discarded. To assess single-unit isolation, we computed the isolation distance (Id), of each cluster. Id quantifies the distance of the nth closest noise spike, if a cluster has n spikes, and measures type I (false-positive) error rates. Cells with an Id of <10 were discarded.
The 24 caudal vibrissae on the right mystacial pad were trimmed to 10 mm in length and inserted 3 mm into the tip of a piezoelectric deflector (one piezo for each whisker, 7 mm from the base). The rest angle of each piezoelectric device was carefully matched to the resting angle of each whisker. For a detailed description of the 24-whisker stimulator, see the original article (8). Four different types of stimulation were used in these experiments: POSW, VELW, ACCW, and sparse noise (one whisker at a time ramp-hold-ramp stimuli). Each stimulus epoch lasted 10.035 s, and the stimulus was continuously given for the full duration of the epoch. There was ~500 ms of rest time between each stimulus epoch where the whiskers were maintained at their rest positions. There were 653 epochs total, 200 for each of the three white noise types (150 nonrepeated and 2 × 25 repeated) and 53 epochs of sparse-noise stimulation. The POSW stimulus was generated by selecting points from a Gaussian distribution with a mean of 0 and a standard deviation (SD) of 0.58°. These points were connected by cubic splines and smoothed to obtain a whisker trajectory that fell within the technical specifications of our piezoelectric benders. The tails of the distributions were cut to avoid ringing artifacts. For VELW and ACCW stimulus types, the frequency spectrum was first fixed with a low-frequency cutoff (~4 Hz for VELW and ~ 16 Hz for ACCW) to make sure that the whisker movements did not exceed the ranges of movement of our piezoelectric benders. We then added variations to the power spectrum of individual epochs by multiplying each power in this ideal frequency spectrum by a factor that was randomly chosen from a standard normal distribution in both the real and imaginary domain. These real and imaginary parts were combined, conjugated, inverse Fourier-transformed, and then standardized by dividing by the SD of the resulting whisker trajectory. This value was then multiplied by the appropriate factor to generate whisker movements with a velocity range that was as closely matched as possible across the different classes of Gaussian stimuli (fig. S2). This was the choice that allowed the frequency spectra to be the most similar in the flattened domains of all three stimuli. Sparse noise stimulations were randomized deflections of one whisker every 50 ms in either the rostral or the caudal direction with a stimulus profile that was a 10-ms ramp –10-ms hold –10-ms ramp with an amplitude of 1.16° and a velocity of 116°/s.
Reverse correlation analysis
All reverse correlation analyses were performed offline and coded in the Python language. Code is available upon request.
First, small windows of stimulus (either 55 or 45 ms before until 5 ms after) around the spikes were tabulated to build the spike-triggered ensembles. For STA, these ensembles were averaged and normalized to make unit vectors. For STC, we downsampled the spike-triggered ensemble in time by a factor of 2 to reduce the number of dimensions. We then whitened the spike-triggered ensemble by multiplying each spike-eliciting stimulus with the inverse covariance matrix of the respective Gaussian white noise but with a ridge-regression smoothing factor applied along the diagonal (POSW, λ = 5 × 107; VELW, λ = 5 × 106; ACCW, λ = 5 × 105). Then, we carried out a singular value decomposition of the covariance matrix of the whitened spike-triggered ensemble. We tested for significance of the filters by applying the same procedure 200 times to shuffled spike trains that maintain the same interspike interval distribution. The singular vectors were considered statistically significant if their eigenvalues exceeded the mean of the shuffled distributions by more than 8 SDs. These methods are exactly the same as prior reports (12, 13).
Feature spaces and functional selectivity
The feature spaces for each GWN stimulus were generated by doing a principal components analysis on all significant filters for the respective stimulus. In all three stimulus types, this yielded a two-dimensional space with almost all variance explained by the top two principal components. Returning to the single cells, each spike-eliciting stimulus was then projected into the respective feature spaces. Each feature space was then binned in two dimensions with polar coordinate bins, and the spike-inducing 2D histogram was normalized by dividing by the average histogram obtained after 1000 shuffles of the spike times while maintaining interspike interval distributions. Feature space subregional selectivity was assessed by randomly shuffling the angle (in polar coordinate representation) of each spike-eliciting stimulus 1000 times. If the vector sum of all real spike-eliciting stimuli was outside of the null distribution (1000 shuffles, P < 0.001), then a cell was called regionally specific. To check for cells that were tuned for opposing regions in the feature space, we multiplied each angle by two and executed the same procedure. In general, the subregion orientation tuning for a cell was taken as the bidirectional value (multiply all angles by 2) because these were less noisy, except in special cases when the cells were extremely one-sided in their selectivity.
Simple stick-slip model
The simple stick-slip model is a linear-nonlinear Poisson (LNP) cascade with one alteration. First, the stimulus is temporally filtered by a skewed Gaussian windowing function that has three parameters (τ, σ, and η). These parameters correspond to the mean, SD, and skewness of the Gaussian window, respectively, with respect to the time of evaluation (τ0). Within this window, another parameter (α) determines the weighting between absolute velocity and absolute acceleration, which is the only departure of the simple stick-slip model from the LNP framework. This filtering part is technically not linear because we compute absolute velocities and accelerations. Because we use covariance approaches to extracting the model filters, this lack of directionality at this stage is not important. After filtering, the distribution of filter products is normalized by subtracting the mean and dividing by the SD of all filter products. Then, the nonlinear function is applied to these normalized values. The exponential nonlinear functions are generated to amplify only the tail of the filter product distribution. The basal rate (1 Hz) is the output unless the filter product is more than 1.5 SDs above the mean filter product value. If the filter product is 5 SDs above the mean, then the max rate (50 Hz) is the output. These values were chosen to be within physiological ranges for barrel cortex (12, 13). In between, the output is scaled to stay between the basal rate and max rate and grows quadratically (the square of the distance from the basal rate zone). It is basically a power 2 nonlinear function that is scaled to pass between a basal rate and a max rate, where the basal rate is increased only when the stimulus falls 1.5 SDs above the mean. After this step, the output is passed through a nonhomogeneous Poisson spike generator to give spikes.
The stick-encoding model is similar to the simple stick-slip model except that the filtering stage is altered. Instead, the stimulus is temporally filtered by a compound Gaussian windowing function that has six parameters (τ1, τ2, σ1, σ2, η1, and η2) defining two windows that are summed and normalized. Within this compound window, two parameters are defined to determine the weighting between directional velocity (α) and absolute velocity (β). The nonlinear stage is identical to the simple stick-slip model.
The sweep-encoding model is similar to the simple stick-slip model except that the filtering stage is altered. It uses the same windows as the stick-slip model, but within this window, three parameters are defined to determine the weighting between absolute velocity (α), absolute position (β), and the interaction between position and velocity (γ). The nonlinear stage is identical to the simple stick-slip model.
Model testing and error minimization
In all cases, the models were run repeatedly for many different combinations of parameters that span the biologically plausible ranges. For the simple stick-slip model, the ranges were α ∈ [0, 1], τ ∈ [−50, −1], σ ∈ [1, 18], and η ∈ [−5, 5]. For the stick-encoding model, the ranges were α ∈ [0, 1], β ∈ [0, 1], τ1 ∈ [−13, −1], σ1 ∈ [1, 9], η1 ∈ [−5, 5], τ2 ∈ [−21, −13], σ2 ∈ [5, 15], and η2 ∈ [−5, 5]. Although these are reported as the final ranges for the parameters, many other possibilities were tested before zooming in on the values that can best generate the observed filters. With eight parameters, there are a huge number of combinations, and to keep computation times low, the final ranges are just large enough to capture the variation seen in the cells. For the sweep-encoding model, the ranges were α ∈ [0, 1], β ∈ [0, 1], γ ∈ [0, 1], τ ∈ [−23, −1], σ ∈ [5, 21], and η ∈ [−5, 5]. The error between the actual filters of a neuron and the model filters resulting from a combination of parameters was computed in one of two ways. If there was a single significant filter, then only models that gave one significant filter were considered. The projection of the model filter on the real filter was computed for both the actual model filter and for its opposite (multiply model filter by −1), and the maximum of these two values was subtracted from one to give the error (fig. S1). This is necessary because STC analysis returns dimensions on which variance is large, which are equally described by a filter or its opposite (in terms of where the variance lives, the dimension is the same if it is multiplied by −1). If there were two significant filters, then this error computation was extended to two dimensions by computing the projection of each model filter on each real filter. Then, each model filter is expressed as an (x, y) coordinate in the two-dimensional space defined by the real filters. The radial distances to these two (x, y) points (one set of coordinates for each model filter) were then subtracted from one and added together to compute the error (fig. S1). The filter-matching procedure searched all model errors to find the combination of model parameters that resulted in the lowest filter errors.
Sparse noise analysis
For the ramp-hold-ramp type of stimulation (sparse noise), the whiskers were stimulated one at a time, and the spikes were then assigned into 1-ms bins with respect to the onset of the stimulation for every whisker separately (from −5 to 55 ms with respect to the stimulus). These PSTHs were then used to assess which whiskers can evoke firing in a cell. To detect stick-like responses, we first subtracted the basal firing rate of the cell from every time bin. The basal rate for each cell was taken as the outlier-removed maximum bin count from blanks that were included in the randomized whisker deflection sequence in the stimulus. Outliers were removed by considering the fifth largest value across all times when summing both stimulation directions. The total number of spikes that occurred above this baseline-subtracted value within 15 ms of the onset of the stimulus was then used as the test statistic to detect significant stick responses. The threshold was that this spike count had to be five times the baseline rate. For detected stick responses, we measure their strength using the PWsh, which is the peak firing rate of the PW PSTH divided by the half width of the peak. To quantify sweep responses, we computed the Z-scored mean firing rates in the 60 ms following stimulation for each whisker. The mean and SD used for the Z score were taken from blank stimulation trials. The Z scores that exceeded 0.5 were then summed and multiplied by a sharing parameter (HW). This sharing parameter was computed as
. This quantity, called the RFsp, was assessed for significance by carrying out 1000 shuffles of full multiwhisker RFs using Poisson firing statistics at the basal rate of each neuron that was determined from the blank stimulation. The RFsp was called significant if it exceeded all the shuffled RFsp values. The PWsh-to-RFsp index was computed by normalizing the indices (divide the entire population by the max) and then computing (PWsh – RFsp) / (PWsh + RFsp). This quantity is 1 if the cell only detects sticks and −1 if the cell only detects sweeps. Cells with intermediate values can be assessed on the basis of the relative strengths of their stick- and sweep-encoding compared to the population of stick and sweep encoders. The w1:w2 ratio was computed for each neuron by adding the total number of spikes that occurred from −5 to 15 ms from the stimulus onset and dividing that by the total number of spikes that occurs from 15 to 55 ms after stimulus onset only for the whisker that evoked the most spikes. This ratio was multiplied by 2 to account for the fact that the numerator is only 20 ms of time and the denominator is 40 ms.
Direction preferences for Gaussian stimuli and sparse noise
For the VELW stimulus, direction preferences were computed by dividing the feature space into two halves. For the stick-encoding cells, the midline was the line orthogonal to the stick axis (the uniphasic filter line defined in Fig. 5B). All values above and to the right of this line were considered caudal, and all values below and to the left of this line were considered rostral. The direction index was then computed as the sum of all caudal values minus the sum of all rostral values divided by the sum of all values. The procedure was similar for the sweep-encoding neurons, except that the midline was the stick axis, with everything below and to the right considered caudal and everything above and to the left rostral. For the sparse-noise stimulation, direction index was just the total spikes to the caudal ramp minus the total spikes to the rostral ramp divided by the total number of spikes to both.
Trial by trial cross-correlations for repeated white noise stimulations
Spikes were placed into 10-ms bins, and correlation coefficients were computed for every pair-wise combination of trials and averaged. For the population, spikes from all recorded single units were used.
CSD analysis and cell locations
Electrophysiological recordings were low pass–filtered with a cutoff of 150 Hz. The signals around the deflections of single whiskers (n = 400 trials for each whisker) were extracted, baseline-subtracted, and averaged. The second spatial derivative was then computed along the electrode sites to obtain the CSD. The cell locations were reconstructed using a barycentric coordinate system based on the energy of the spike waveform at every electrode site (33). Briefly, vectors pointing to each electrode site were scaled on the basis of the spike waveform energy at that site, and the summation of all these vectors gave the final position of the cell. Electrode sites with less than 5% of the total waveform energy were not included in the computation.
Acknowledgments: T. Deneux provided important guidance for the modeling approaches. L. Estebanez established the reverse correlation methods in the laboratory and provided valuable advice for the experiments and analysis. V. Ego-Stengel and I. Férézou gave insightful advice and feedback throughout the project and on the manuscript. We thank G. Hucher for performing the histology and A. Daret for providing the animal care for the rats. Funding: This work was supported by Equipe Fondation pour la Recherche Médicale (FRM) DEQ20170336761 (to D.E.S.). E.R.H. was supported by the Paris-Saclay University (Lidex NeuroSaclay) and by the International Human Frontier Science Program Organization (CDA-0064-2015 to B.B.). M.A.G. was supported by the Paris-Saclay University (Lidex NeuroSaclay) and the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement no. 702726. Author contributions: E.R.H. created the modeling framework, designed the stimuli, did the spike sorting and analysis, prepared the figures, and wrote the manuscript. M.A.G. did the extracellular recordings, tested the stimulus, did the analysis, contributed to the modeling framework, and edited the manuscript. B.B. oversaw the analysis, motivated the modeling framework, and wrote the manuscript. D.E.S. led the project and edited the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.