Fast cortical dynamics encode tactile grating orientation during active touch.

Touch-based object recognition relies on perception of compositional tactile features like roughness, shape, and surface orientation. However, besides roughness, it remains unclear how these different tactile features are encoded by neural activity that is linked with perception. Here, we establish a cortex-dependent perceptual task in which mice discriminate tactile gratings on the basis of orientation using only their whiskers. Multielectrode recordings in the barrel cortex reveal weak orientation tuning in average firing rates (500-ms time scale) during grating exploration despite high levels of cortical activity. Just before decision, orientation information extracted from fast cortical dynamics (100-ms time scale) more closely resembles concurrent psychophysical measurements than single neuron orientation tuning curves. This temporal code conveys both stimulus and choice/action-related information, suggesting that fast cortical dynamics during exploration of a tactile object both reflect the physical stimulus and affect the decision.

INTRODUCTION

Touch-based object recognition is essential for guiding behavior in a wide variety of environmental conditions. Reliable recognition generally depends on tactile search behavior executed with appendages like fingers for humans or the mystacial vibrissae for rodents (). The vibrissae, or whiskers, are rooted on the snout in densely innervated follicles, where mechanosensitive cells transduce whisker bending and contact forces into electrical signals (). The resulting sensory information has spatial (across whiskers) and temporal aspects that are integrated as it passes through several distinct somatosensory pathways before reaching the barrel cortex and other areas (). As the foremost recipient of primary somatosensory thalamic afferents (), the barrel cortex is seen as the major cortical hub for the processing of whisker-based tactile information (, ).

Numerous physiological studies on how barrel cortex neurons respond to simple, reliably targeted whisker stimuli have pointed toward a somato-topographical code based on high-velocity deflections of one or a few whiskers (). However, a simple velocity-based coding framework is not sufficient to support some of the reported perceptual functions of the barrel cortex. In head-fixed mice that use only a single row of spared whiskers (all others trimmed) to explore their proximal surroundings, barrel cortex neurons encode the precise location of a pole placed in reach of the whiskers and perturbation of barrel cortex activity causes mice to fail to recognize when the pole is in rewarded locations (, ). The encoding of pole location is thought to arise by combining information about contact events (high angular velocity and high angular acceleration whisker deflections) with the state variables of ongoing whisking behavior (i.e., phase and set point). These state variables fall in coding dimensions that do not exist when simple, reliably targeted stimuli are applied to otherwise passive whiskers. This highlights the need to study barrel cortex representations in conditions in which they serve a perceptual function. In this direction, it has been shown that the barrel cortex is not required to detect the simple presence or absence of objects in the proximal surroundings with a single spared whisker (, ), while the barrel cortex is essential in more demanding task conditions like the discrimination of sandpapers (–) or whisker-mediated gap crossing (). Together, perceptual studies suggest that the barrel cortex is critical for precisely localizing and recognizing tactile objects with the whiskers.

To turn these perceptual insights into an understanding of the barrel cortex coding principles underlying tactile object recognition, it is then necessary to finely dissect barrel cortex representations during the various perceptual feats that it is required to perform. So far, most studies have focused on how the coarseness of anisotropic surface textures (sandpapers) is encoded during exploration with one or a few whiskers (–). Even with such reduced whisker inputs, these studies found that object coarseness is encoded by temporal integration of whisker slip events, with higher rates of slip events causing higher firing rates in the barrel cortex (–). While this has provided insights into tactile coding, coarseness is just one of many features that can differ between objects. Along with variations in coarseness, natural objects also exhibit unique combinations of isotropic features, which means they can be decomposed into an arrangement of oriented surfaces (). While freely moving rats can discriminate oriented tactile gratings with their whiskers (), it is not known whether and on what time scale grating orientation is encoded in the barrel cortex during active sensation.

To address this, we developed a cortex-dependent Go/NoGo task in which head-fixed mice can use all of their whiskers to discriminate tactile gratings on the basis of orientation. While it is common practice in whisker-based perceptual studies to reduce the whisker input complexity by extensive trimming (usually one spared whisker and, in some cases, one row of spared whiskers), we decided to establish this task in animals with untrimmed whiskers. While this along with the three-dimensional nature of whisker search behavior complicates tracking of individual whisker movements (), it provides an opportunity to assess how information is encoded in the barrel cortex when it is not funneled through one or a few intact whisker-barrel channels. Multielectrode recordings during task performance revealed that while peak cortical firing rates occurred early on during grating exploration, stable orientation selectivity in the average firing rates of individual units in this period was weak. However, orientation information could be extracted from the time course of population activity (at a temporal resolution of 100 ms) during grating exploration by support vector machine (SVM) classifiers, which decoded the orientation category in line with concurrent psychophysical measurements. These results suggest that orientation information is first encoded in the temporal dynamics of population activity in the barrel cortex, and only later, if at all, do individual units develop orientation tuning in average firing rates (on 500-ms time scales). By analyzing perceptual errors, we found that cortical dynamics contain both sensory and choice/action-related information, confirming that, in perceptual tasks that require intact cortical processing, the barrel cortex encodes both sensory and decision-related information.

RESULTS

Mice categorize tactile gratings on the basis of orientation

To investigate whether mice are able to discriminate tactile gratings, we trained head-fixed, water-deprived mice (Fig. 1A and fig. S1) to report the orientation of a tactile grating by licking a tube to receive a water reward. The oriented gratings were presented in full dark conditions using a linear stage, and whisker interactions with the gratings were filmed with a high-speed infrared video camera (see Materials and Methods). For each trial, after no licking was detected on the reward port for at least 3 s, a 2-kHz sound was played to signify trial onset and a grating was translated into reach of the right whisker field (Fig. 1B). After a 1-s period of interaction with the grating, mice reported the orientation of the grating by either licking to obtain a water reward (Go trial) or refraining from licking to avoid punishment (Fig. 1B). In these trial conditions, mice were trained to perform a simple Go/NoGo discrimination between a vertically oriented grating (90°) and a horizontal grating (0°), with Go and NoGo stimulus types interchanged in different groups of animals (fig. S1 and Materials and Methods for all training details). After performance of simple Go/NoGo discrimination stabilized above 70% correct across 2 days, intermediate orientation angles spaced by 9° were gradually introduced and reinforced (Fig. 1B, see Materials and Methods). In this psychometric version of the task, the boundary between rewarded and nonrewarded orientations was 45°, and the fully ambiguous orientation was never presented.

Fig. 1.

Mice categorize tactile gratings on the basis of their orientation.

(A) A schematic showing the behavioral setup. (B) Task parameters for (1) simple Go/NoGo discrimination and (2) psychometric Go/NoGo grating orientation discrimination tasks. (C) Top: Lick raster from simple Go/NoGo discrimination between vertical (90°) and horizontal (0°) gratings. Middle: Lick probabilities from one session of simple Go/NoGo discrimination. Bottom: Mean learning curves (shaded areas are SEM, n = 13 animals) for All (black), Go (green), and NoGo trials (red). (D) Top: Lick probabilities from one session of psychometric Go/NoGo discrimination. Bottom left: Psychometric functions for two groups of mice where the Go/NoGo rules were interchanged. Bottom right: The psychometric functions controlled for motivation (error bars are SEM). Inset: Boxplot of d′ values.

From the beginning of simple Go/NoGo (0° versus 90°) training, mice quickly learned the appropriate time to lick, and after 10 to 15 days (~2000 trials), they discriminated between orthogonal grating orientations as measured by their licking behavior (Fig. 1C). Improved performance across time was mostly attributable to refraining from licking for the NoGo stimuli (Fig. 1C). After progressing to the psychometric version of the task, the ongoing motivational state of the animal, driven by thirst, determined whether false alarm or miss errors were more common. In most animals, we observed a more gradual change in licking behavior across orientation steps for NoGo than for Go orientations (Fig. 1D, lick histograms). Along with this, mice tended to make more false alarm errors than miss errors (Fig. 1D), indicative of a strategy aiming to minimize reward loss. This strategy results in asymmetric psychometric functions (Fig. 1D). To balance these curves (), we averaged across animals in which the Go and NoGo orientations had been interchanged, and this revealed that the discrimination performance controlled for motivation is almost perfectly symmetric (Fig. 1D). These results indicate that head-fixed mice can discriminate tactile gratings using only their whiskers, and they do so with high acuity.

Intact cortex is essential for discriminating oriented tactile gratings

We next asked whether cortex is essential to perform the simple Go/NoGo version of this task. Optogenetic manipulations can perturb performance even if a brain area is dispensable (), so we opted for a lesioning strategy. Mice were trained in the simple Go/NoGo version of the task until they reached stable performance above 70% correct across 2 days, after which thermocoagulation lesions () were applied that were centered on the contralateral postero-medial barrel field (Fig. 2A). As a control, another group of animals (sham group) underwent mock surgeries that involved the same duration of anesthesia, a large craniotomy centered over the barrel cortex, and the same process to reseal the exposed brain but with no lesion. The day after surgery, both lesion and sham groups performed the simple Go/NoGo task at chance levels (Fig. 2B), indicating that the general aftereffects of surgery and craniotomy have an impact on performance, as is the case even in whisker-based perceptual tasks that do not depend on cortical processing (). The same is not true for auditory discrimination tasks, where task performance is unperturbed the day after auditory cortex lesions if the task does not depend on intact cortical processing (). Over the ensuing days, the sham group steadily recovered performance, while the lesioned group continued to perform at chance levels (Fig. 2B and fig. S2). Lesions were examined post hoc in coronal sections to ensure that all postero-medial barrels (straddlers, A1-E4) in the whisker region of the primary somatosensory cortex had been removed (fig. S3). The lesions could also include other areas in the parietal regions of the cortex such as secondary somatosensory (S2) or posterior parietal cortex (PPC/PTLp), which also process whisker-related information, so the observed drops in performance could also reflect loss of these areas.

Fig. 2.

Cortex is required for discriminating tactile gratings.

(A) Top: Experimental approach and timeline. Bottom left: An example cortical lesion centered on the barrel cortex. Bottom right: the corresponding slice in the brain atlas. Abbreviations: third ventricle (3V), dentate gyrus (DG), and lateral ventricle (LV). (B) Simple Go/NoGo discrimination performance before and after surgery in lesion and sham groups (P = 0.0392, P = 0.0056, and P = 0.0076 correspond to the three asterisks from left to right where an asterisk indicates P < 0.05, bootstrap resample test, n = 5 for each group except on day +4, where there are four sham animals and two lesioned). (C) Whisker tracking during performance of the task. Top left: A binarized frame. Top right: A manually selected region of interest (ROI) containing the bases of the whiskers and the centroid (blue). Bottom left: The velocity of the centroid plotted across a trial. Bottom right: The resulting whisking envelope after rectification and smoothing. (D) Top: Average whisking envelopes across days for lesion and sham groups. Exploratory whisking was quantified between trial start and the time point indicated by the dashed line. (E) Baseline subtracted area under the curve (AUC) during the whisker search period across days for lesion and sham groups. (F) Performance broken down by trial type for lesion and sham groups. Error bars are standard errors.

Perturbation of barrel cortex activity is known to affect whisker movement control (, ). Therefore, we examined high-speed videos of whisker movements executed by the animals during task performance in sham and lesion groups. To quantify global whisker movements throughout a trial, we defined the whisking envelope as the rectified and smoothed centroid velocity of the binarized whisker image within a manually traced region of interest (ROI) around the whisker bases (Fig. 2C, see Materials and Methods). This envelope showed that whisking behavior is most pronounced between trial onset (the trial start sound cue) and the time when the grating is fixed and within reach of the whiskers (Fig. 2D, gray shaded rectangles). Surgery affected the average whisking envelope in both sham and lesion groups of animals, as quantified by integrating the whisking envelope between trial onset and grating halt (Fig. 2E, area under the whisking envelope curve). By day 3 after surgery, the total whisking of both groups returned to presurgical levels, but the behavioral performance only recovered in the sham group. Therefore, the drop in task performance after lesion cannot be explained by deficiencies in whisker control. Cortical lesions also did not affect performance by abolishing licking. On day 3 after surgery, hit rates and false alarm rates were equal in the lesioned animals (both at ~50%), indicating that mice randomly licked rather than never licking at all (Fig. 2F and fig. S2). These results suggest that intact cortical processing is required to discriminate oriented gratings with the whiskers, and this cannot be explained by changes in whisker search behavior or licking ability.

Exploratory whisking and peak cortical activity precede discriminative choice

To study the encoding of grating orientation in the barrel cortex during active discrimination, we made acute extracellular recordings (five recordings, 48 single units, and 152 multi-units) during task performance (Fig. 3A and fig. S4). Silicon probes with linearly spaced electrodes (spanning 775 μm) were lowered to 1-mm depth from the surface of the contralateral barrel cortex (targeted C2 whisker A/P: −1.5 mm, M/L: 0/3.3 mm). Electrode placement in the barrel cortex was histologically verified in tangential sections after the experiments (fig. S4), and most of the active cells that were recorded resided in deeper layers (fig. S4). All mice showed aggregated task performance above 70% correct on the day of the recording.

Fig. 3.

Exploratory whisking followed by an increase in barrel cortex activity precedes discriminative choice.

(A) Schematic of the recording setup. (B) Example hit trial showing licks, reward, whisking, and unit rasters. Whisking is normalized to baseline levels, which were computed as the average whisking envelope in the first 60 ms of the trial. (C) Same as (B) for a correct rejection trial. (D) Performance across trials and psychometric functions for five mice, example mouse in magenta, and average across animals in bold. (E) Licking, whisking, and total population spiking activity for an example mouse (left) or all mice (right). Shaded lick histograms are licks after reward/punishment. Shading around curves is SEM. tw is the time that the trial averages of whisking were at a peak, tfr is the peak firing rate, and td is the time that licking became discriminative. (F) Trial averages for whisker angle, absolute whisker velocity, and mean curvature across trial time. Tracking was done with WHISK software, significantly assessed with a permutation test that shuffles trial labels. For mean curvature (κ), asterisk indicates that P < 0.001 (difference in the minimums of the trial-averaged traces was never seen in 1000 shuffles).

In an example hit trial (Fig. 3B), the mouse initiated whisking before the grating came into reach and spiking activity increased once the grating was close enough to touch the whiskers. After ~500 ms of exploration, the mouse decided to lick and received a water reward, which triggered prolonged licking. In an example correct rejection trial (Fig. 3C), the mouse whisked into the grating, which produced peak cortical activity in the last ~250 ms before the grating stabilized at its fixed position in reach. Then, the mouse correctly withheld licking to avoid punishment. This same behavioral sequence was apparent when averaging across all trials in this animal or across all performing animals (Fig. 3, D and E). As the grating approached, mice executed whisker search behavior, which was followed by a burst of spiking activity in barrel cortex neurons that peaked just before the grating stopped near the snout. Licking behavior was initiated after the grating stopped and the average licking differences between orientation classes (> or <45°) could be discriminated ~590 ms after the peak of population activity (Fig. 3E, see Materials and Methods). After the decision to lick, low whisking levels were maintained and, in some mice, a rebound of whisking and barrel cortex activity was observed when the texture began to move away from the mouse (Fig. 3E).

To look closer at the whisker interactions with gratings during task performance, we applied several tracking algorithms (fig. S5) (, ). The WHISK software () provided the closest match to manual tracking (fig. S5) and allowed the extraction of whisker angles (averaged across whiskers in focus), absolute whisker velocities, and average curvatures across all trials (Fig. 3F). While we found no significant differences for whisker angles or velocities between grating orientations, there was a significant difference in mean curvature (the average curvature of all visible whiskers in a frame). Vertical orientations caused lower negative peaks in the average curvature at grating halt (Fig. 3F), indicating that whisker bending is more pronounced for grating angles >45°. Together, these data highlight the general pattern of behavior (search, find, and lick) and spiking activity in the barrel cortex during task performance and suggest that differential whisker bending could act as the cue to discriminate grating orientations.

Orientation tuning in average firing rates is weak during grating exploration

To study how barrel cortex activity encoded grating orientation, we first examined individual units. Some single units discharged many spikes at the onset of whisker interactions with the grating for all different orientations, and their firing rates decreased when the grating reached its fixed position (Fig. 4A and fig. S4, single unit 1). Other units had less pronounced responses during early exploration but still had elevated firing rates while the gratings were within reach (Fig. 4A and fig. S4, single unit 2). To quantify whether unit responses were selective for a particular grating orientation, we constructed tuning curves in 500-ms windows at different latencies with respect to trial onset (Fig. 4B and fig. S6, blue-magenta gradient). Tuning curves were computed for units by summing the spikes within a 500-ms window of interest for each trial of a given orientation and then dividing the total number of spikes for each trial by the size of the window (Fig. 4B, two examples; fig. S6, four examples). To assess tuning significance, we expressed the tuning curves in polar form and compared the magnitude of the vector sum to the vector sums obtained from shuffling the trial labels 200 times (Fig. 4B, bottom left). If the actual tuning vector was beyond the 95th percentile of the shuffled distribution, the unit was considered significantly tuned. For example, single unit 1 showed no selectivity during the first interactions with the grating (Fig. 4B, top left, blue) or after discriminative choice (Fig. 4B, top left, magenta). Single unit 2 was tuned to grating orientation (Fig. 4B, right). The tuning to 18° gratings originated in the second 500-ms time bin (before decisive licking in this animal measured across all trials, see Materials and Methods) after the grating came into reach of the whiskers, and it persisted through decision and feedback.

Fig. 4.

Orientation tuning is scarce during early grating exploration.

(A) Single-trial rasters and trial-averaged PSTHs from two example single units during performance of the psychometric Go-NoGo task. (B) Tuning curves from two example single units [same units from (A)] in various time windows color-coded by where the 500-ms time bin in which the curves were computed falls within the trial. Asterisks indicate that the unit is significantly tuned to orientation (P < 0.05 shuffle test) in the corresponding time window. (C) The percentage of responsive cells and orientation-tuned cells plotted versus trial time. Population firing rates and licking behavior are plotted below to serve as a reference. Number of single units per mouse was 10, 14, 14, 1, and 9, and the multi-unit counts were 40, 23, 43, 18, and 28, respectively.

Examining the orientation tuning across time for all units, we found that it was not above chance levels during the peak of firing that occurs during the first interactions with the grating (just before grating halt) and started to appear in the period between grating halt and discriminative licking (Fig. 4C). This suggests that in the exploratory period, information about grating orientation is not yet well encoded by average firing rates of single neurons, as might be the case in other sensory modalities ().

After decision and feedback in the form of reward or punishment, differential licking and whisking behaviors can explain the apparent increase in orientation tuning (Fig. 4C) because if mice perform the task, these behaviors are strongly correlated with grating orientation. To investigate this, we fit the spiking data for each unit with a generalized linear model (GLM) containing two orientation terms (binary variables for horizontal or vertical gratings), terms for whisker angle, absolute whisker velocity, and mean whisker curvature (average values in 50-ms windows before each time point), licking terms, and interaction terms between these individual variables and grating position (fig. S7, Materials and Methods, all binned at 50-ms resolution). By comparing these full GLMs with reduced GLMs that did not contain orientation terms (just a single term for grating in reach), we found that 75/200 recorded units were modulated by orientation (likelihood ratio test, χ2 test with df = 1, P < 0.01), but that most, but not all, of these units were tuned for Go stimuli (fig. S7C). Together, these results suggest that orientation tuning in the average firing rates of individual units is scarce before decision and increases after decision; however, in this later stage, it is confounded by correlated changes in licking and whisking behavior.

Temporal decoders outperform rate decoders during grating exploration and provide the closest match to psychophysical measurements

Because tactile inputs occur in a sequence of multi-whisker contacts that evoke dynamic cortical responses, we hypothesized that during exploration before decision, grating orientation information could be carried by coordinated population activity sequences rather than by the firing rate of orientation-tuned units. To search for orientation-specific activity sequences, we looked at the co-firing patterns for all simultaneously recorded units in an example mouse. Co-firing matrices were computed separately for horizontal and vertical trials in the time period leading up to grating halt (Fig. 5A, top, 5× 100-ms bins), which is also the time when we detected curvature differences between grating orientations (Fig. 3F). Of the 37 simultaneously recorded units in this example mouse, 9 had stronger co-firing for horizontal gratings and 8 neurons had stronger co-firing for vertical gratings (Fig. 5A, top, red and blue regions in the co-firing matrix), while the rest of the population had equal co-firing for horizontal and vertical trials in the 500-ms period leading up to grating halt. We examined the time course of the trial-averaged response separately for these two subgroups of neurons, and we observed distinct temporal sequences of activity for vertical and horizontal trials in both subpopulations (Fig. 5A, bottom, strongest in the group coactivated by horizontal gratings). This observation suggests that coactivation sequences can carry some information about grating orientation.

Fig. 5.

Temporal decoders outperform rate decoders during grating exploration and provide the closest match to psychophysical measurements.

(A) Top: Co-firing matrices from an example population (n = 37 units) as the grating approached. After computing the difference between vertical and horizontal matrices, two groups of units were apparent. Bottom: Time course of group firing rate as grating approached. (B) Schematic showing how the SVM classifiers were trained and tested for one example mouse. Classifier performance is aligned with the licking behavior and the population firing rates. The end of the exploratory period is defined as the 500 ms before discriminative licking (across all trials) and the after-feedback period is the same size time window after reward or punishment. (C) Four classifier types defined by their bin arrangements were tested to assess the contribution of fast cortical dynamics in the sample populations of neurons. Performance, licking behavior, and population firing rates are all aligned. Total performance is the AUC that is above chance in a period. * means P < 0.05 in a bootstrap resample test. ns, not significant. (D) Psychometric and neurometric functions for the three main classifier types in the exploratory period or after feedback (inset).

To generally examine whether temporal patterns of population activity contain grating information, we trained SVM classifiers that take as input the number of spikes that each unit in the population produced in the five preceding 100-ms time bins (Fig. 5B, same example mouse). After leave-one-out cross-validation, the performance of this temporal decoder in classifying vertical versus horizontal gratings was above chance well before the onset of licking (Fig. 4). Classifier performance before decision was markedly reduced by removing the 17 units with orientation-specific co-firing patterns (Fig. 5B, brown), and training only on these 17 units gave performance that was indistinguishable from the whole population (fig. S8), confirming that the observed co-firing patterns provide orientation information during early grating exploration. To broadly evaluate the importance of temporal patterns of activity, we compared the ability of four different decoders to discriminate horizontal and vertical gratings in all task-performing mice: (i) the temporal decoder using the sequence of activity over the preceding 5× 100-ms time bins, (ii) an average firing rate decoder using the sum of activity over a single preceding 500-ms time bin, (iii) an instantaneous firing rate decoder taking activity over a single preceding 100-ms time bin, and (iv) a decoder using the same 5× 100-ms bin sequence as the temporal decoder but with the time bins shuffled across trials to control for dimensionality but destroy temporal patterning. Across all mice, the decoders with fine temporal resolution outperformed the other decoders during exploration, with the best performance coming from the temporal decoder with a five-bin history. Increasing the temporal resolution of the classifier binning did not affect these conclusions (fig. S8). This demonstrates that temporal sequences of population activity contain more information about grating orientation than the average firing rates before decision. This advantage vanished after decision and feedback, when differences in licking and whisking behavior drove discernible differences in cortical activity (Figs. 4 and 5C).

We next examined which decoder best followed concurrent psychophysical measurements by looking at classifier performance across different grating orientations (Fig. 5D). During exploration, the temporal decoder again outperformed rate decoders and showed the closest resemblance to the psychometric behavior (Fig. 5D, right). This indicates that precise temporal coactivation patterns provide a sensory coding space that can underlie this perceptual behavior. After decision, temporal and rate decoders were equally predictive of psychophysical measurements, which could reflect either the pronounced differences in licking behavior between grating classes in this time window or a reformatting of the sensory code.

Temporal sequences of cortical activity during grating exploration encode sensory and choice/action-related information

We next examined whether the activity preceding decision in the barrel cortex reflects the external sensory object or the upcoming choice of the animal by measuring to what extent it predicts single trial outcomes, which, in our case, were hits, misses, false alarms, and correct rejections. Because most animals performed very few misses, we concentrated on hit, false alarm, and correct rejection trials. In an example mouse, we observed that false alarm trials were associated with higher population firing rates than correct rejections and that the differences were not only attributable to licking (Fig. 6A, compare licking at the bottom to population firing rates on top, and Fig. 6B, left). However, these differences in population firing rate were no longer visible when looking across the cohort of five performing mice, with the only persevering differences across trial outcomes being more whisking and increased barrel cortex activity after punishment for false alarm trials (Fig. 6B, right). A more detailed look at the whisking parameters revealed that there were no significant differences in whisker angle, velocity, or curvature for different trial outcomes before decision, but there were differences after decision (Fig. 6C).

Fig. 6.

Temporal sequences of cortical activity during grating exploration encode sensory as well as choice-related information.

(A) Top: Trial-averaged population firing rates separated by orientation and outcome for an example mouse (hits, misses, false alarms, and correct rejections are lime green, dark green, red, and purple, respectively). Bottom: Lick rasters for all trials for the same mouse. (B) Left: Trial-averaged whisking envelope, population firing rate, and licking behavior separated for hits (lime green), correct rejections (purple), and false alarms (red) for the example mouse. Right: Trial-averaged whisking envelope, population firing rate, and licking behavior separated for hits (lime green), correct rejections (purple), and false alarms (red) for all mice. (C) Average whisker angle, absolute velocity, and mean curvature across the trials for hits, false alarms, and correct rejections. (D) SVM classifier decoding for hits versus false alarms (black) and hits versus correct rejections (gold). Top: Total performance in the exploratory and after feedback periods for the different classifiers normalized to the baseline before trial start (P = 0.0214 and P = 0.035, asterisks indicate P < 0.05 using paired bootstrap resample test). Middle: Classifier performance across all time points relative to trial onset. Bottom: Licking behavior for hits and false alarms. (E) Same as (D) for false alarms versus correct rejections (P = 0.018 and P = 0.0134, asterisks indicate P < 0.05 using paired bootstrap resample test).

Next, we assessed whether choice/action-related activity was present by comparing classifier performance on hits and false alarm trials (Fig. 6D). Temporal decoders trained to discriminate hit versus false alarm trials performed worse in the exploratory period than decoders trained to discriminate hit versus correct rejection trials. This indicates that on trials where the mouse failed to withhold licking for the NoGo orientation, the classifiers were also more likely to err, suggesting that part of the cortical activity pertains to the choice to lick, lick preparation, or licking itself. The existence of residual performance for hit versus false alarm classification suggests that some of the activity was purely sensory, faithfully encoding the grating stimulus regardless of the decision. To confirm the existence of the sensory component, we trained temporal decoders to discriminate false alarm and correct rejection trials and found that this discrimination was also worse during grating exploration than hit versus correct rejection classifier performance (Fig. 6E). This further indicates that licking differences alone do not drive classifier performance. Together, these analyses establish that temporal sequences of population activity in the barrel cortex encode both sensory and choice/action-related information during grating exploration.

DISCUSSION

We have shown that head-fixed mice can discriminate tactile gratings on the basis of orientation using only their whiskers (Fig. 1) and that this perceptual task depends on cortical circuits (Fig. 2). Our findings indicate that leading up to decision, there is more information about grating orientation in the fine temporal dynamics (100-ms time scale) of barrel cortex population activity than in average firing rates (500-ms time scale) of individual orientation-selective neurons, and decoding this information at higher temporal resolution better reflects the psychophysics of an animal’s perceptual decisions (Figs. 4 to 6). This means that the higher temporal resolution afforded by electrophysiological recordings is likely to be crucial for understanding the recognition of compound tactile objects composed of sets of oriented edges, because as we have shown, fast temporal variations in barrel cortex activity provide important parts of the initial object representation in S1. This is in stark contrast with the well-studied single whisker detection tasks (, , ), which can easily be solved on the basis of average firing activity of single barrel cortex neurons. In more challenging task conditions, we propose that tactile information in S1 first appears in fast temporal fluctuations at a time scale of at most 100 ms, which are further integrated to produce a final decision. There are some precedents to this concept in both rodents (–) and primates. During whisker-based coarseness discrimination, precise spike timing and average firing rates in barrel cortex neurons convey stimulus information that is used by rats to guide decisions (). Rats are even highly sensitive to spike timing resulting from stimulation of single neurons (). In flutter discrimination tasks in macaques, most S1 neurons finely track the precise frequency of the vibration and thus encode with high temporal resolution (). However, a relatively smaller number of cells are less locked to the vibration and instead linearly increase their average firing rates with increasing vibration frequency (). Therefore, both temporal and average firing rate codes in S1 can provide the information required for the monkey to solve the task, and psychometric behavioral analysis suggests that the monkeys rely on a mix of the two, because they underperform the ideal observer of the temporal code and outperform the ideal observer of the rate code (). Along with these studies, our findings point toward a model where temporal information initially present in S1 is likely transformed either within the primary cortex itself or in concert with other areas to generate integrated firing rate representations, and both of these coding spaces are then linked with tactile perception.

To discriminate tactile gratings, mice need to execute a sequence of appropriate behaviors, and our head-fixed conditions ensure that these behaviors are precisely locked with trial cues. First, mice need to detect the incoming grating. In this pursuit, we found that mice whisked vigorously when the grating was approaching, and vertical gratings caused more pronounced whisker bending at grating halt (Fig. 3). This active search behavior generated temporally dynamic responses in barrel cortex neurons (Figs. 3 to 5). Freely behaving rodents can use head movements along with exploratory whisking to perform tactile search (), so the level of vigorous whisking observed here might be an adaptation to head-fixed task conditions that could emphasize temporal aspects of the coding. Once the approaching grating is fixed in reach, mice continue their search behavior and adaptively sample the grating. We focused our analysis on the time period between first contact and decision, which allowed us to show that orientation tuning in the firing rates of single neurons was scarce (Fig. 4), while temporal sequences of population activity could be used in decoding schemes that gave a closer match with concurrent psychophysical measurements. When we examined choice behavior in single trials, we found that the temporally dynamic code in the barrel cortex contains both sensory (what the object was) and choice/action-related (decision, lick preparation, or lick-related) information (Fig. 6); however, due to the nature of the Go/NoGo paradigm, we could not distinguish the contributions of licking or lick preparation and choice. To solve this problem, two-alternative forced choice (2-AFC) versions of the task will need to be developed. Despite this, we observed that false alarms were more difficult for classifiers to discriminate from hit trials in the exploratory period than correct rejections, and differences in licking behavior alone do not drive classifier performance (Fig. 6). This is consistent with detection, discrimination, and delayed match to sample tasks showing that barrel cortex neurons encode both stimulus- and choice-related information (, , –). Our results extend this finding to tasks involving the discrimination of oriented tactile gratings and suggest that the temporal structure of cortical dynamics can also contain choice-related information.

Discriminating oriented gratings with multiple whiskers is a much more challenging task than detecting the vibration or contact of a single whisker, which is reflected by the essential role of the cortex in this task and the temporally dynamic encoding that is present in the barrel cortex during task performance. Stimulus complexity aside, when all whiskers can be used, sensory information is spread across the whisker pad. While this lack of control is not ideal for precisely quantifying single whisker contributions to perception, we are still able to recover both sensory and choice/action-related information from sampled barrel cortex populations. This spread of information across sensors more closely resembles the computational problems that the barrel cortex faces in natural conditions than the vibration of a single whisker. The fact that perceptual decisions can still be linked with the discriminability of sampled barrel cortex activity suggests that both sensory and choice/action-related information are widespread in the cortex even if the same whiskers might not be involved on every trial. In summary, our results establish a cortex-dependent tactile discrimination task in which the fast cortical dynamics are informative and lay out how temporally structured co-firing events in subpopulations of neurons can support grating orientation perception. There is much to learn about the circuits that are responsible for transforming temporal population codes into stable firing rate codes () that can be found in downstream areas () and how the interaction between cortical areas facilitates tactile object recognition.

MATERIALS AND METHODS

Animal care

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 9714-2018011108392486). Mice were housed in cages in groups of two to four individuals with food available ad libitum on a 12-hour light/12-hour dark cycle with temperature kept at 23°C. During behavioral training, animals received 800 to 1000 μl of water through the lick port per day. To make sure this was the case, reward size was regularly calibrated (at least once per week), and if animals did not receive rewards that add up to above 800 μl total during training, then they were given supplemental water through the lick port at the end of the session. Bottles were placed in the home cages after behavioral training on Fridays to give the mice water ad libitum until Saturday afternoon, when bottles were removed.

Behavioral setup

Mice were trained in a custom-built behavioral setup that was interfaced using a National Instruments (NI) card (USB-6343) to control a linear stage (Newmark eTrack series) that brought the gratings within reach of the whiskers and an Arduino Uno to control stepper motors (Makeblock) for adjusting the orientation angle of the grating and a solenoid valve (LVM10R1-6B-1-Q, SMC) for delivering water rewards (5 to 8 μl). Sound cues were played with loudspeakers (Labtec spin 85 speakers). Licking signals were acquired and digitized using a capacitive sensor (Sentronic AG, SK-3-18/2,5-B-VA/PTFE) before being fed into the NI card. Software to carry out the training protocols and log the licking data was coded in Matlab using the data acquisition toolbox.

Headpost implantation

To stabilize the animals in the behavioral apparatus, a head-fixation post was implanted along the midline of the skull. Mice (C57BL/6) that were 6 to 8 weeks old (20 to 26 g) were anesthetized by intraperitoneal injection of a mix of ketamine (Ketasol; 80 mg/kg) and medetomidine (Domitor; 1 mg/kg). Once the mice were insensitive to hindpaw pinch, they were placed on a nose clamp and their eyes were kept moist with Ocry-gel (TVM Lab). Body temperature was maintained at 36° using a thermal blanket. Xylocaine was injected under the skin in the center of the skull near bregma. Fur in the surgical location was removed using Veet, and a long incision was made in the skin along the midline of the skull 10 min after Xylocaine injection. After being fully exposed, the dorsal surface of the skull was scratched with a scalpel to create striations. The scratched skull was then cleaned with hydrogen peroxide. A head-fixation post was glued in place along the midline using cyanoacrylate and then the exposed skull and base of the post were covered with Super-Bond (C&B, Sun Medical Co. Ltd.). The implant and all exposed surfaces were then embedded in dental cement. After everything had solidified, the mice were injected in one of the hindlimbs with 15 μl of atipamezole (Antisedan, Orion Pharma) and transferred to a recovery cage that was placed on a heating blanket. Mice recovered for at least 1 week before any further manipulation.

Orientation discrimination training protocol

Mice were weighed every day during water deprivation periods to make sure they did not fall below 80% of pre-deprivation body mass. For 2 days before training, mice were fully water-deprived. On the first day of training, the mice were placed on the head-fixation post for 10 min in the dark with the lick port in reach. They were then given single water rewards (5 to 8 μl) randomly until they started to lick regularly at the lick port. Once they were comfortable licking the lick port for water reward, a protocol was launched that made one reward possible every 10 s if the animal licked to initiate the delivery, for up to a maximum of 100 rewards. After this habituation (1 to 2 days, 1 hour per day), the animals were given trials only with the Go grating until they licked regularly at the correct time within single trials. The trial timeline is shown in Fig. 1. For the first 40 trials, rewards were given automatically 1 s after the grating came into reach of the whiskers. After these free rewards, the mice had to lick in a 2-s window that started 1 s after the grating came into reach to receive the reward. The starting threshold to trigger reward was a single lick, which was then increased to as high as four (two to four across all mice) licks to trigger a reward. If animals performed three misses in a row, the next Go trial automatically was rewarded, and this “miss” counter was reset while the trial was still scored as a miss. Once the rewards were action-contingent within the trial framework, performance was tracked. When the animals were able to perform 70% correct across an entire training session, a NoGo stimulus was introduced the next day interleaved pseudo-randomly with the Go grating at a ratio of three Go trials for every one NoGo trial. If the addition of NoGo trials and their associated punishments (white noise at 60 to 80 dB and time out of 5 to 30 s) did not cease reward-seeking behavior, the ratios were equilibrated (50% Go 50% NoGo) on the next day of training. The first NoGo stimulus was a flat surface (a small circle of printer paper glued on a disk the same dimensions as the gratings) with no grating (fig. S1). Once the animals discriminated this flat surface from the Go grating (fig. S1, performed 70% correct across 200 trials in a single day), the NoGo stimulus was changed to a grating orthogonal to the Go grating. Punishments (loudness of the white noise and length of the time out) and lick thresholds were increased if animals could not refrain from licking for NoGo gratings. Punishments began with only 60-dB white noise, which was increased until 80 dB if a mouse did not reduce false alarm rates. If sounds were not enough, then time outs were incorporated starting at 5 s and increased to 30 s if needed. After 2 days of 70% performance in discriminating orthogonal gratings, intermediate grating orientations were introduced. At first, only four intermediate orientations (9°, 18°, 72°, and 81°) were given, but then another 4 (27°, 36°, 54°, and 63°) were added after performance stabilized above 70% correct. For the full psychometry, a single training session contained 40 trials for each extremity (0° and 90°) and 20 trials for each intermediate grating, for a total of 240 trials.

Task performance and psychophysics analysis

Learning curves across trials were calculated by dividing the number of correct responses (hits + correct rejections) in the preceding 25 trials by 25. Across days, the curves were stitched together and smoothed with a Gaussian kernel. If the animals ceased licking for more than 15 trials, the trials were removed from the learning curves, as blocks of inactivity of this size indicate that the mouse is distracted or satiated. Discriminative licking was detected by Wilcoxon rank-sum tests on the licking histograms (100-ms bins) generated for each trial (significance for P < 0.01) comparing horizontal trials (<45°) with vertical trials (>45°) at each time bin. The first bin with a significant difference was taken as the “discrimination time.” Psychometric functions in Fig. 1D were taken from 2 days of task performance (480 trials). The criteria for selecting these days were that total performance was above 70% correct across the entire day and the miss/false alarm ratio was between 0.5 and 1.5, indicating a balance between thirst and punishment avoidance. d′ was computed as Z(Hit/(Hit + Miss)) − Z(FA/(FA + CR)), where Z(p), p ∈ [0,1] is the inverse of the cumulative Gaussian distribution ().

Cortical lesions and histology

After all mice learned the full psychometric version of the task except for one that was only trained on the simple version of the task, they were anesthetized (1.5% isoflurane delivered with SomnoSuite, Kent Scientific) and placed in a nose clamp. A thermal blanket kept body temperature above 36°C. Ocry-gel (TVM Lab) was applied to the eyes to keep them from drying out. The location of the C2 barrel had been marked on the skull (A/P: −1.5 mm, M/L: 0/3.3 mm) from the headpost implantation surgery in these mice, and this mark was used as the center of a 3- to 4-mm-diameter craniotomy. Thermocoagulation lesions were carried out with a fine-tipped cauterizer, making sure not to touch the surface of the brain, but to bring the cauterizer just close enough to blacken the exposed cortical tissue containing the barrel field. The craniotomy was then covered with Kwik-Cast (World Precision Instruments) and then sealed with dental cement. Sham animals underwent the same surgical procedure except they did not receive thermocoagulation lesions. After surgery and recuperation (~1 hour in a recovery cage), mice were given 250 μl of water and returned to their home cages. The behavioral testing began again the day after surgery. When behavioral testing was complete, lesioned mice were transcardially perfused with saline followed by a 4% formaldehyde solution in 0.1 M phosphate buffer. Brains were dissected and then postfixed overnight at 4°C. After washing with phosphate-buffered saline, brains were cut into 80-μm coronal slices. Slices were mounted and then imaged using a Nikon Eclipse 90i microscope (Intensilight, Nikon) and Nikon Plan UW objectives (1×/0.04 W.D 3.2 or 2×/0.06 W.D. 7.5). Slices were then manually aligned with the Paxinos mouse brain atlas and the lesioned areas were tracked along the anterior-posterior axis to make sure they covered the posterior-medial barrel field (fig. S2). Sham mice were used later for electrophysiological recordings during task performance, after which their brains were treated in the same way, except they were sliced tangentially to reveal electrode locations with respect to the barrels (fig. S3). Electrode tracks in these preparations were visible because 1,1’-dioctadecyl-3,3,3′,3′-tetramethylindocarbocyanine perchlorate (DiI) was placed on the shanks before they were inserted into the brain.

Whisker movement tracking

During some sessions, high-speed videos of the whisker interactions with the gratings were filmed with an infrared video camera (Baumer; 500 fps). The frames were grabbed on the same clock as the stimulus presentation to assure synchronization. Four different methods of whisking tracking were then carried out. In manual tracking, a whisker was selected that remained within focus for the duration of a trial, and every five frames, the base and the half length of the whisker were manually demarcated. For the centroid ROI method, an ROI was manually selected around the bases of the whiskers that were in focus. In this ROI, the centroid of the binarized whiskers was computed, and this centroid was then projected onto a line that was perpendicular to the rostral whiskers to give a single coordinate. The velocity of the centroid coordinate across frames was rectified and smoothed to give the whisking envelope. This quantifies the global rostral-caudal movement of all the whiskers. This procedure is graphically displayed in Fig. 2C. For the pivot ROI method, an ROI containing the bases of all the whiskers and a pivot point along the snout of the animal were demarcated, and the ROI was rotated about this pivot point to find the optimal match to the next frame (). Another approach was using the WHISK software (). For each frame, angles, velocities, and curvatures of detected whiskers (~12 per frame) were averaged to obtain a global measurement for the whisker pad. Normalization to whisking levels in the first 30 frames (first 60 ms of a trial) was sometimes applied to compare across mice with different levels of baseline whisking activity.

Electrophysiological recordings

On the day of the recording, mice were briefly anesthetized (30 min, 1% isoflurane delivered with SomnoSuite, Kent Scientific), and the dental cement that was covering the craniotomies from the sham surgery (n = 5 sham animals) was removed. In four other experiments, fresh craniotomies were drilled following the same protocols described in the lesion section above (except no lesions). After durectomy, the exposed cortical surface was moistened with fresh Ringer’s solution and then covered with Kwik-Cast (World Precision Instruments), which was secured in place with cyanoacrylate. The mice were then allowed to recover for 2 to 3 hours in a cage that was placed on a heating blanket. Mice were then placed in the behavioral setup and the Kwik-Cast was carefully removed, making sure not to damage the brain in the process. Multielectrode silicon probes (A2x32 5mm-25-200-177, NeuroNexus) that had been coated with DiI were then slowly lowered into the left hemisphere barrel cortex at about 2 μm/s. Once they reached a depth of 800 to 1000 μm and sufficient spiking activity was seen across all channels, the preparation was stabilized for 20 min before the behavioral protocol was launched, with periodic water rewards given to keep the mice awake and unstressed. In three mice, intermediate orientations were rewarded or punished and the number of trials for each orientation followed the protocol detailed in the orientation discrimination training section. In two mice, intermediate orientations were given as catch trials, and in these experiments, fewer intermediate orientation trials were given (90 horizontal trials, 90 vertical trials, and 5 catch trials for each of four intermediate orientations). Psychometric data were pooled across these five mice for the electrophysiological dataset. For the behavior alone (Fig. 1), all animals followed the same protocols that are described in the orientation discrimination training section.

Data processing and analysis for electrophysiological recordings

Extracellular signals were acquired at 20 kHz with an Intan RHD2000 recording system. The raw data were median-filtered to remove common mode noise from all channels and then passed into Kilosort2 for spike detection and clustering. Clusters were manually curated to pick out waveforms with physiological shapes that decay with distance from a primary electrode (electrode with the largest magnitude waveform). The units that passed visual inspection and entered the analysis pipeline were both single units and multi-units depending on the refractory periods found in their autocorrelograms. Data from single and multi-units were pooled for all analyses, with single units representing the activity of single neurons and multi-units representing the activity of local populations that were close enough to the electrode to yield waveforms that resembled single units. Trial-averaged spiking histograms were created by binning spikes in 50-ms bins (Fig. 3). Normalized firing rates were computed by dividing by the baseline firing rate, which was taken as the mean firing rate across 500 ms beginning 1 s before trial onset.

Orientation tuning and response detection

Orientation tuning curves were constructed by breaking trials up into 500-ms blocks. For each unit and each 500-ms block, the trial-averaged firing rates for a stimulus of a given orientation determined the magnitude of the vector pulling in that direction in a polar coordinate system where all the orientation angles were multiplied by 2. The vector sum of these 10 (or 6) oriented vectors (0°, 18°, 36°, 54°, 72°, 108°, 126°, 144°, 162°, 180° or 0°, 36°, 72°, 108°, 144°, 180°) was compared to the distribution of vector sums obtained by shuffling the trial labels 200 times. If the actual vector sum was outside of the sphere defined by 95% of the 200 shuffles (P < 0.05), then the cell was called orientation tuned in that 500-ms block. False-positive rates were thus kept at 5%. Cells were deemed significantly responsive if evoked firing rates were 5 SDs above the baseline firing rate. For detecting modulation of neural activity by orientation using GLMs (fig. S7), we used a modeling approach that was established to separate head direction and place tuning in entorhinal cortical cells () and adapted to the whisker system (). The following terms were explicitly fit in the models: (1) a binary term for vertical orientation in reach, (2) a binary term for horizontal orientation in reach, (3) average whisker angle in the last 50 ms, (4) average absolute whisker velocity in the last 50 ms, (5) average mean whisker curvature in the last 50 ms, (6) lick counts two bins before (50-ms bin size), (7) lick counts one bin before, (8) lick counts in the same bin as a spike, (9) lick counts in the bin after the spikes, (10) lick counts two bins after the spikes, and (11 to 18) interaction terms for grating in reach and terms 3 to 10. In reduced GLMs, all terms were the same except there was only a single binary term for grating in reach rather than two orientation specific terms (1 and 2 above). All data were binned at 50-ms resolution (spikes, licks), and whisker movement data (angle, absolute velocity, and curvature) were averaged across 50-ms windows. With these two models, likelihood ratio tests were performed using the χ2 distribution (df = 1) and any cell that performed better with orientation terms (P < 0.01) was deemed modulated by orientation (). We built the design matrices on the basis of all trials (fig. S7), and they were fit using the scikit-learn glm.fit() method with a Poisson link function.

Defining the exploratory period and the after-feedback period

Significant differences in licking behavior were assessed by binning the digital lick signal counts into 100-ms bins. Then, the distributions of Go trial licks and NoGo trial licks were compared at each time point relative to trial onset using Wilcoxon rank-sum tests, and the first time point in the trials that gave a significant result with P < 0.01 is where the mouse was said to have licked discriminatively. For each mouse, the 500 ms before this time point was counted as the early period. The late period was a fixed period after reward or punishment that was chosen to maximize differential licking behavior, after licking had ceased for false alarm trials. This was purposefully chosen to show how a code driven by marked differences in licking behavior would present itself.

Co-firing analysis

Spikes were binned into 100-ms bins, and the 5× 100-ms bins leading up to grating halt were used to calculate co-firing. To do this, for each pair of units, the dot product of their activity profiles (five bins of activity on a given trial) was computed and averaged across trials of a stimulus type (horizontal or vertical). Then, the difference was computed between trial types (vertical − horizontal), and this matrix was manually clustered into neurons with different co-firing levels for the respective stimulus classes.

SVM classifiers

Spikes were placed into 100-ms bins to generate population vectors of various types for each trial (Fig. 5A). The trials were divided either by grating orientation (Fig. 5) or by trial outcome (Fig. 6). Binary nonlinear SVMs were then trained using the scikit-learn module in python along with the leave-one-out protocol in the model selection subdirectory of this module. The nonlinear classifiers used a gamma function with an input parameter of 1/n_features (the “auto” option from the sklearn documentation). The population vectors were moved one step at a time (always the smallest step present, 100 ms for Figs. 5 and 6), and for each time step in the trials, the classifiers were retrained on the basis of the corresponding subspaces of the population vectors that ended at that time step. The performance was the percentage of all trials correctly classified. Each trial was left out only once. Performance curves were smoothed with a three-bin kernel that took the average across the three bins and assigned that value to the central bin. For the shuffling procedures, a population vector was generated for each trial at a particular time point in the trial as an nUnits × nTimeBins matrix (rows = units × columns = time bins in Fig. 5C). The columns of this matrix were shuffled differently for each trial, so if 1, 2, 3, 4, 5 is the correct column order, random permutations for each trial would look like: {trial 1: 2, 1, 5, 3, 4}; {trial 2: 4, 3, 1, 5, 2} etc., while the rows were left in their normal order. This shuffling procedure changes the temporal order of spiking activity while maintaining total spikes across time for each unit (row) in the matrix, which should not have an effect on performance if grating orientation is encoded by average firing rates at the scale of ~500 ms. To compute classifier performances when trial counts were imbalanced (Fig. 6), we computed the performance separately for each trial type and then averaged this to obtain total performance.

Paired bootstrap resample test

For small sample sizes (n = 5) that are common in challenging experimental conditions such as these, the most accurate statistical test is nonparametric bootstrap resampling. Wilcoxon and Mann-Whitney tests rely on approximations for small sample sizes (typically when n < 20). To carry out this test, we resampled 1000 times with replacement from the pool of N (usually five) mice and permuted the labels of what was being tested (lesion versus sham, temporal decoders versus average firing rate decoders, etc.). When appropriate, the permutations were done while keeping the measurements paired. If the difference of the mean values obtained was > or <95% of the shuffled resampled mean differences, then the measurement was deemed significant with P < 0.05. Exact P values are provided as averages of five different resamples composed of 1000 shuffles each.