Transformation of Perception from Sensory to Motor Cortex.

To better understand how a stream of sensory data is transformed into a percept, we examined neuronal activity in vibrissal sensory cortex, vS1, together with vibrissal motor cortex, vM1 (a frontal cortex target of vS1), while rats compared the intensity of two vibrations separated by an interstimulus delay. Vibrations were "noisy," constructed by stringing together over time a sequence of velocity values sampled from a normal distribution; each vibration's mean speed was proportional to the width of the normal distribution. Durations of both stimulus 1 and stimulus 2 could vary from 100 to 600 ms. Psychometric curves reveal that rats overestimated the longer-duration stimulus-thus, perceived intensity of a vibration grew over the course of hundreds of milliseconds even while the sensory input remained, on average, stationary. Human subjects demonstrated the identical perceptual phenomenon, indicating that the underlying mechanisms of temporal integration generalize across species. The time dependence of the percept allowed us to ask to what extent neurons encoded the ongoing stimulus stream versus the animal's percept. We demonstrate that vS1 firing correlated with the local features of the vibration, whereas vM1 firing correlated with the percept: the final vM1 population state varied, as did the rat's behavior, according to both stimulus speed and stimulus duration. Moreover, vM1 populations appeared to participate in the trace of the percept of stimulus 1 as the rat awaited stimulus 2. In conclusion, the transformation of sensory data into the percept appears to involve the integration and storage of vS1 signals by vM1.

Introduction

In 1979, Whitfield [1] surveyed the first century of neuroscience and pointed to clear evidence that, after ablation of sensory and association regions of cerebral cortex, animals can still perform many forms of sensory discrimination. The conserved capacities must reflect the ability to access and act upon sensory information in the ascending, subcortical pathways. By contrast, cortical ablation causes a permanent deficit whenever the behavioral task requires the elemental sensory signals to be integrated, distributed, and stored, so as to connect with the traces of previous experience. In sum, intracortical processing serves to transform a representation of mere physical characteristics into the perception of things that are “out there” ([1], p. 146) in the world. In the intervening years, investigators have exploited ever-improving techniques to reinforce the notion that intracortical processing is critical to perceptual functions [2]. For instance, when rats [3] or monkeys [4] were required to report the presence of a tactile stimulus around the detection threshold, sensory cortical activity varied little across trials, whereas the triggering of all-or-none activity cascades in frontal cortex correlated with awareness of the stimulus [5].

In the present study, we examined the mechanisms and functions of intracortical processing by employing a perceptual task where the subject’s choice depended not on stimulus detection but on the statistical content of the stimulus. We trained rats to compare the intensities of pairs of “noisy” (stochastic) vibrations presented to the whiskers. Varying the durations of the two stimuli, we found that a longer-duration vibration was perceived as more intense. The divergence between the magnitude of the stream of incoming sensory data (which remained, on average, stationary) and the rat’s percept of stimulus intensity (which grew with stimulus duration) allowed us to explore the transformation of neuronal representations. In vibrissal primary somatosensory cortex (vS1)—also known as barrel cortex, the entry stage of whisker signals to cerebral cortex—neurons encoded the instantaneous speed of the whisker vibration; at this stage, we found a “local” representation of the stimulus. In primary vibrissal motor cortex (vM1)—the main frontal cortex target of vS1—neurons encoded the nonlinear temporal integral of stimulus speed. This second, “global” representation matched the rat’s percept.

Results

Perceived Intensity of a Vibration Depends on Duration

To gain insight into the neuronal processing that transforms a stream of sensory events into a percept, we measured rats’ judgements of tactile vibrations delivered to their whiskers (Figure 1A, left). To assess whether the perceptual phenomena generalize to humans, we carried out the same studies with stimuli delivered to the fingertip (Figure 1A, right). Subjects were required to compare the intensities of two vibrations—stimulus 1 of duration T1 and stimulus 2 of duration T2—separated by a delay of 2 s (Figure 1B) [6]. Vibrations were constructed by stringing together over time a sequence of velocity values, sp, sampled from a Gaussian distribution. We consider the stimuli as speed rather than velocity, so the distribution took the form of a folded half-Gaussian (right side of Figure 1B). A single vibration was thus defined by its nominal mean speed, denoted sp (equivalent to the SD of the Gaussian multiplied by ), and its actual, observed mean speed (sum of sp over entire stimulus, divided by T). Because the stimulus was stochastic, judgement of overall intensity required some form of integration over time.

Figure 1

Delayed Comparison Task: Stimuli and Performance

(A) The rat held the tip of its snout in the nose poke with whiskers resting on the motor-actuated plate. The human subject rested his/her arm on a cushion with left index fingertip in contact with the motor-actuated rod.

(B) Example easy (upper) and difficult (lower) stimulus pairs. Subjects compared stimulus 1 (mean speed sp1 and duration T1) to stimulus 2 (sp2 and T2), with intervening interstimulus delay. Vibrations were generated by stringing together velocity values sampled from a Gaussian distribution. On the right side of the panel, speed distributions (the folded, half-Gaussian of velocity) are illustrated for stimulus 1 (dashed) and stimulus 2 (solid). Horizontal lines indicate sp1 (dashed) and sp2 (solid). In these examples, the easy and difficult trials differ only in stimulus 2.

(C) Stimulus set for rats and humans. Dashed diagonal line is decision boundary, where sp1 = sp2. The sp1, sp2 pair for each trial was drawn from one box; arrows indicate the stimulus pairs from which the difficult and easy stimuli of (B) were drawn. Inside each box is the percent correct for that stimulus pair, averaged across all subjects. For pairs distant from the diagonal, accuracy was good except for (25, 47) in rats and (14, 21) in humans, where sp1 assumed its lowest value, most likely explained by “contraction bias” [6, 7, 8]. Green rectangles enclose the stimulus set used to generate psychometric curves.

(D) Average psychometric curves of rats and humans (black) and the ideal observer (gray). Error bars represent SEM across subjects.

The subject’s response was logged as correct according to the relative values of nominal mean speed: sp2 > sp1 or else sp1 > sp2. For a given value of sp1, a trial could be easy or difficult according to the value of sp2 (Figure 1B, upper and lower examples). Moreover, across trials, different combinations of sp1 and sp2 were presented such that neither stimulus 1 nor stimulus 2, taken alone, contained sufficient information to solve the task; the subject was required to execute a direct comparison on every trial [6, 9]. In Figure 1C, the coordinates of each box represent joint values of sp1 and sp2. Here, T1 and T2 both were 400 ms. The percent correct averaged across all sessions for 16 rats (left panel) and 12 humans (right panel) is given inside each box.

Trial difficulty—the difference between stimulus 1 and 2—was defined as the normalized speed difference (NSD), (sp2 − sp1) / (sp2 + sp1). NSD ranged from −0.3 to 0.3 for rats and from −0.2 to 0.2 for humans. As a measure of perceptual acuity, for the closely spaced stimulus pairs (inside the green rectangles, Figure 1C) we computed the percent of trials in which the subject judged stimulus 2 > stimulus 1 as a function of NSD and fit the data with a logistic function to generate psychometric curves. The average curves across all rats and humans are given in Figure 1D (black curves). These may be compared to the performance of an imaginary “ideal observer” (gray curve) who measures actual mean speed with no error. The ideal curve is not a perfect step function because vibrations were built stochastically; the nominal mean speed of a vibration, sp, and the actual, observed mean speed (sum of sp over entire stimulus, divided by T), could disagree, causing even the ideal observer to make errors.

We manipulated stimulus durations to determine whether perceived intensity was constant or time dependent. In the first experiment, stimulus 1 and stimulus 2 were of equal duration within a trial but duration varied randomly across trials (T1 = T2; balanced condition). Figure 2A shows the average psychometric curves of rats (left panel) with T1 and T2 set to 100, 200, 400, or 600 ms (progressively brighter shades of blue; same color code used in all panels) and humans (right panel) under the same conditions. Rats’ overall performance was inferior to that of humans, but both species yielded markedly steeper psychometric curves as durations of the stimulus pair increased from 100/100 up to 600/600 ms.

Figure 2

Experimental Evidence for Temporal Integration

(A) Psychometric curves averaged across 11 rats (left) and 6 humans (right) in the balanced condition (T1 = T2). Colors indicate stimulus duration: 100 ms (dark blue); 200 ms (blue); 400 ms (teal); and 600 ms (cyan). Error bars in this and all panels indicate SD across single subjects.

(B) Left: changes in performance (percent correct and maximum psychometric curve slope) for 200, 400, and 600 ms duration stimulus pairs in relation to the average value across all durations. Average among rats represented by circles, humans by squares. Colors indicate stimulus duration using the key of (A). Right: change in percent of trials judged stimulus 2 > stimulus 1 and change in the inflection point.

(C) Psychometric curves averaged across 11 rats (left) and 10 humans (right) in the unbalanced condition, where T1 = 400 ms and T2 ranged between 200 (blue), 400 (teal), and 600 (cyan) ms. Curves are shifted rightward for shorter T2 and leftward for longer T2.

(D) The same as (B) but for unbalanced condition where T1 was fixed to 400 ms.

(E) Psychometric curves averaged across ten rats (left) and three humans (right) in the unbalanced condition, where T2 = 400 ms and T1 ranged between 200 (blue), 400 (teal), and 600 (cyan) ms. Curves are shifted leftward for shorter T1 and rightward for longer T1.

(F) The same as (B) but for unbalanced condition where T2 was fixed to 400 ms.

To quantify the effect of duration on acuity, for each subject, we calculated the percent correct across the full psychometric curve as well as the psychometric curve maximum slope (a proxy for perceptual acuity). Then, we computed the change in both values for stimuli of 200, 400, and 600 ms duration in relation to the average values across all durations. The left panel in Figure 2B illustrates the average values across subjects. Both in rats and in humans, greater stimulus duration led to greater slope and accuracy. Besides undergoing a change in slope, psychometric curves may also shift laterally. To test for this, for each subject, we projected the psychometric curve inflection point down to its corresponding value of NSD. We also computed the percent of trials across the entire psychometric curve for which the subject judged stimulus 2 > stimulus 1, a second measure of the left/right shift, for if the percentage is decreased/increased then that curve is most likely shifted to the right/left, respectively. Then, we computed the change in both measures in the T1/T2 = 200/200 ms and 600/600 ms conditions, with respect to the 400/400 ms condition. The right panel in Figure 2B illustrates the average values across subjects. The plot reveals no systematic duration-dependent lateral shift.

We also conducted sessions where T1 was of 400 ms duration and T2 could be (randomly across trials) either 200, 400, or 600 ms. With “unbalanced” durations, rats (Figure 2C, left panel) exhibited clear shifts in the psychometric curves. The rightward shift of the T1/T2 = 400/200 ms curve means that truncating T2 decreased the likelihood of the choice stimulus 2 > stimulus 1; by the same token, the leftward shift of the 400/600 ms curve means that lengthening T2 increased the likelihood of the choice stimulus 2 > stimulus 1. The results for humans (Figure 2C, right panel) reveal similar effects. In sum, longer-duration stimuli were judged as more intense, shorter-duration stimuli as less intense. As before, we quantified the effect of duration on acuity for all subjects (Figure 2D, left panel). Neither in rats nor in humans did variation in T2 lead to systematic changes in slope or accuracy. By contrast, variation in T2 did lead to systematic shifts in psychometric curve position (Figure 2D, right panel).

We conducted sessions where T2 was of 400 ms duration and T1 could be (randomly across trials) of either 200, 400, or 600 ms duration. The effects of varying T1 (Figure 2E) were symmetric with those of varying T2 (Figure 2C), meaning that the effect of unbalanced durations on perception did not depend on stimulus order and was thus a perceptual bias rather than a choice bias. Again, we quantified the effect of duration on acuity for all subjects (Figure 2F, left panel) and found that variation in T1 did not lead to systematic changes in slope or accuracy. Variation in T1 led to systematic shifts in psychometric curve position (Figure 2F, right panel).

In summary, rats and humans can accumulate events from a stochastic tactile vibration to form an increasingly reliable estimate of magnitude as stimulus duration grows, provided the stimuli to be compared are of equal duration. In the case of unequal stimulus durations, the briefer stimulus is felt as less intense and the longer stimulus as more intense. The shift in perceived intensity in relation to duration will be the main focus of the study.

Perceived Intensity Is Accounted for by Nonlinear Summation of Instantaneous Speed

What form of integration leads to the perceptual judgements reported above? Because earlier work argued for a greater contribution of the early response in primary somatosensory cortex [10], we hypothesize that sp values are weighted by a function that decays exponentially after stimulus onset, a model we refer to as “summation by exponentially weighted primacy” (SEWP). SEWP is illustrated schematically in Figure 3A. The instantaneous speed of stimulus 1 at time t, sp (upper panel, black traces), is multiplied by the weighting function (middle panel), where is an individual subject’s intrinsic integration time constant. The final perceived intensity of stimulus 1 (lower panel) is proportional to the sum of the weighted values of speed from stimulus onset until vibration termination at time T. The perceived intensity of stimulus 2 is derived in the same manner and is compared to the memory of stimulus 1 (dashed gray line). In this example, sp2 is slightly greater than sp1; T2 is of 600 ms (cyan) or else 200 ms (blue). Because sp1 and sp2 are close, the subject’s choice would be affected by the relative durations, T1 and T2.

Figure 3

Model of Summation by Exponentially Weighted Primacy: Predictions and Tests

(A) Upper plot: two example trials to schematically illustrate the model. In both cases, T1 = 400 ms (black trace), but T2 = 600 (cyan) or 200 ms (blue). Middle plot: SEWP posits the existence of a weighting function that acts from t = 0 until T (solid black trace for stimulus 1 and cyan and blue traces for stimulus 2). Lower plot: summation from t = 0 until T of the product of instantaneous speed, sp, and the weighting function yields the final perceived intensity of the stimulus. The perceived intensity of stimulus 2 is derived in the same manner and is compared to the memory of stimulus 1 (dashed gray line).

(B) Dashed lines show the psychometric curve fit (as Figure 2) on experimental data of one rat and one human. For the same subjects, the psychometric curve fit on output of the SEWP model, averaged over 200 train/test cycles, is given by the solid lines. Stimulus durations are indicated by the color key.

(C) Investigation of the accuracy of the SEWP model in predicting experimental data. A set of rats and humans was tested using a larger set of stimulus durations. Predictions of the SEWP model closely fit the observed experimental data.

(D) Same analysis as in (B) but for the unbalanced condition with T1 fixed at 400 ms and T2 of variable duration.

(E) Same analysis, with T2 fixed at 400 ms and T1 of variable duration.

See Figure S1 for SEWP model parameters, Figure S2 for accuracy of the SEWP model, and Figure S3 for the experiments in rats and humans to test the type of temporal integration.

The SEWP model can be implemented for individual rat and human subjects by obtaining the following five parameters: the exponential decay constant (τ), overall acuity , the bias (μ), and the upper and lower asymptotes. We fit the data by nonlinear least squares (see STAR Methods). Half of all trials served as the training dataset to estimate these parameters; then, for the other half, the model predicted the percent judged stimulus 2 > stimulus 1 for every stimulus difference (NSD) according to T1 and T2. The train/test cycle was repeated 200 times. Parameter values for all rats and humans are given in Figure S1.

Figure 3B shows psychometric curves fit to the model prediction (solid lines) and experimental data (dashed lines) for a single rat (left panel) and human (right panel) for different stimulus durations in the balanced condition. The close fit between the model and experiment attests to the model’s accuracy in predicting the improved acuity with increased stimulus duration. Comparison between model and experimental data for all subjects is given in Figure S2A.

In the balanced condition, we collected data from six rats and four humans with stimulus duration varying from 100 to 600 ms. Figure 3C shows the experimental results (points) and model prediction (curves) plotted simply as percent correct. Absolute value of NSD ranged from 0.1 (dark red) to 0.2 for humans and 0.3 for rats (both yellow). Model predictions closely fit the experimental data.

Next, we consider the unbalanced conditions. Following the same format as Figure 3B, Figure 3D shows the psychometric curves fit to experimental data and the model prediction for a rat (left panel) and a human (right panel) with T1 fixed at 400 ms and T2 of variable duration. Model and experimental data both indicate a duration-dependent bias—greater T2 led to an increase in the likelihood of the choice stimulus 2 > stimulus 1. The model also accurately predicted the unbalanced condition in which T2 was fixed at 400 ms and T1 was of variable duration (Figure 3E). Results from all subjects with variable T2 and variable T1 are shown in Figures S2B and S2C, respectively.

The SEWP model posits that the weighting function acting on sp is maximum at stimulus onset and then diminishes over time. However, a “recency” function, anchored to stimulus offset (as opposed to onset), would give on average the same quantity of weighted, summated speed. We carried out experiments in rats and humans to test which type of temporal integration is at play, and the observed psychometric curves (Figure S3) match the predictions of the primacy model, consistent with SEWP as the form of temporal integration in both species.

Cortical Encoding of the Whisker Vibration

In the remainder of this article, we examine neuronal activity in behaving rats from two regions (Figure S4)—vS1, the main entryway of sensory signals to the cerebral cortex [11, 12], and vM1, the target of projections from sensory cortical regions and a source of dense projections to major motor output systems [13, 14]. To elucidate the overall profile of firing, Figures 4A and 4B show the peristimulus time histograms (PSTHs) constructed from the pooled activity of single and multiunits from vS1 and vM1, respectively. The plots show trials with duration of stimulus 1 (green traces) set to 600 ms and stimulus 2 (red traces) set to 400 ms; trials with different values of sp were merged. Firing rate in vS1 was characterized by a sharp peak at stimulus onset, followed by a steady-state level 300–600 ms after stimulus onset, and a reduction during the interstimulus delay. vM1 also showed increased firing during stimulus presentation; unlike vS1, activity was maintained across the 600-ms interval between the end of stimulus 2 and the go cue.

Figure 4

Neuronal Coding of Vibration Speed

(A) Averaged peristimulus time histogram (PSTH) of 139 single and multiunits from vS1. The PSTH was generated with 1-ms bins and then convolved with a Gaussian kernel ( = 30 ms).

(B) PSTH generated from 319 vM1 single and multiunits.

(C) Left: PSTH depicting a vS1 multiunit response to stimulus 1 with trials separated according to sp1 (color bar). Only trials with T1 set to 400 ms are illustrated. Right: histogram of Pearson correlation coefficients between firing rate and vibration sp among all vS1 single and multiunits.

(D) Same as (C) but for vM1. In the left panel, an example response to stimulus 2 is shown.

See Figure S5 for comparison between stimulus 1 and stimulus 2 in vS1 and vM1

Both cortical regions encoded the vibration mean speed. The PSTH of a vS1 multiunit cluster in response to stimulus 1 (Figure 4C, left) reveals a higher firing rate for higher sp (dark red to yellow), with Pearson correlation coefficient ρ (computed between sp and firing rate from stimulus onset to offset) of 0.62 (p < 0.01; t test). We refer to a significant variation in overall firing rate in relation to sp as “global coding”. The distribution of correlation coefficients among vS1 single and multiunits (Figure 4C, right) demonstrates a high proportion of units characterized by “positive” global coding (ρ > 0). Responses to stimulus 1 are above the abscissa; significant correlations (p < 0.01; t test) are green (94 out of 208 units; 45%). Responses to stimulus 2 are below the abscissa; significant correlations are red (114 out of 208 units; 55%). Less than 3% of units showed significant “negative” coding (decreasing firing rate with increasing sp; ρ < 0; p < 0.01; t test).

The PSTH of a vM1 multiunit cluster in response to stimulus 2 (Figure 4D, left) also reveals a higher firing rate for higher sp (ρ = 0.75; p < 0.01; t test). The distribution of correlation coefficients among vM1 single and multiunits (Figure 4D, right) demonstrates a high proportion of units characterized by positive global coding (ρ > 0). As before, responses to stimulus 1 are above the abscissa; significant correlations are green (82 out of 348 units; 24%). Responses to stimulus 2 are below the abscissa; significant correlations are red (94 out of 348 units; 27%). Less than 6% of units showed significant negative coding. Individual neurons in vS1 and vM1 showed consistent strength of global coding for stimulus 1 and stimulus 2 (Figure S5).

Whereas neurons in both cortical regions encoded vibration speed, they differed in their integration timescales. The high degree of temporal precision characteristic of many vS1 neurons is evident in the upper raster plot of Figure 5A, which illustrates a single unit response to a vibration waveform presented in many different trials (“frozen noise”). The temporal alignment of spikes across trials reflects heightened firing probability in response to specific local stimulus features. The lower raster plot shows the same unit’s response for a different frozen noise waveform: with rearranged vibration waveform, the temporal configuration of the spike train was likewise rearranged. The lack of temporal precision characteristic of most vM1 single units is evident in Figure 5B. Presented with the same two instances of frozen noise, no temporal patterns emerged.

Figure 5

Comparison of Local Coding in vS1 and vM1

(A) A fixed stimulus waveform (“frozen noise,” upper trace) was presented in 40 trials, with different sp intermixed. The raster plot illustrates the spikes of a vS1 single unit. Below, a second example of frozen noise and the corresponding raster plot for the same vS1 single unit is shown.

(B) Same as (A) but for a vM1 single unit.

(C) Cross-correlogram (green) between the spike trains of the vS1 single unit (convolved with a Gaussian kernel of SD 2 ms) and the vibration waveform. Gray trace shows the cross-correlogram between the spike trains of the same vS1 single unit and the temporally reversed stimulus.

(D) Same as (C) for the vM1 single unit.

(E) Distribution of local coding index values—the difference between the peak of the spike train cross-correlogram with the actual and the reversed stimulus waveform. Open bars represent vS1 neurons, and filled bars represent vM1 neurons. Only neurons with significant index (p < 0.01; permutation test; bootstrapped 200 times) are illustrated.

(F) Distribution of Gaussian kernel widths that maximize the forward/reversed cross-correlation difference. Open bars represent vS1 neurons, and filled bars represent vM1 neurons.

See Figure S6 for detail explanation of cross-correlation analysis used.

We quantified the strength of local coding by measuring the cross-correlation between the spike train, convolved with a Gaussian kernel, and the sp vibration waveform (Figure S6; STAR Methods). A neuron that fired reliably after the occurrence of specific vibration features, like that of Figure 5A, will show a marked peak in the cross-correlogram (Figure 5C, green trace); as a control, we computed the cross-correlation with temporally reversed stimuli (Figure 5C, gray trace) where, by definition, local coding must be absent. Whereas the vS1 unit showed a large forward versus reversed cross-correlation difference, the vM1 unit showed no such difference (Figure 5D, green versus gray traces). We take cross-correlation difference as an index for the magnitude of local coding. In Figure 5E, the distribution of significant local coding index values for vS1 (above the abscissa) and vM1 (below the abscissa) is plotted. In vS1, 99 out of 208 (48%) showed significant local coding, and the index values were large. In vM1, 14 out of 336 (4%) showed significant local coding, and the index values were smaller.

With what temporal precision were vS1 and vM1 spike trains correlated with the vibration? Figure 5F shows, across the full dataset, the distribution of widths of the Gaussian kernel that maximized the forward versus reversed cross-correlation difference. Only instances where the forward versus reversed difference was significant (p < 0.01; permutation test) are illustrated. Whereas vS1 spike trains (represented by bars above the abscissa) tended to be correlated with the stimulus at a temporal precision of better than 10 ms, most vM1 spike trains (represented by bars below the abscissa) were correlated at a precision of 10–50 ms.

vM1, but Not vS1, Exhibits Temporal Integration Matching the Percept

The behavioral studies (Figures 1, 2, and 3) revealed that rats and humans have a robust capacity to judge stimulus intensity, but the intensity percept is modulated by stimulus duration. Therefore, the neuronal substrate for perception must entail not only a representation of sp but also of duration, T. In this section, we compare the vS1 and vM1 representations of both stimulus properties. First, we fit neuronal activity by a linear model in order to quantify the relations between firing rate and the following variables: sp1, sp2, and T1. The linear model (STAR Methods) reveals the extent to which variation in neuronal firing could be accounted for by variation in a selected variable. For this analysis, the dataset included trials with T1 of 200 or 600 ms; T2 was 400 ms. For sliding 300-ms time windows, we computed the percent of neurons (single and multiunits) whose firing encoded each variable of interest; Figures 6A and 6B show the results from vS1 and vM1. To allow depiction of 200- and 600-ms trials together, we split stimulus 1 into two sections: the left section is aligned to stimulus onset (denoted 0 ms) and the right section is aligned to stimulus offset (again denoted 0 ms). The percent of neurons expected to code stimulus features by chance (see STAR Methods) is subtracted so that even small values of percent are meaningful.

Figure 6

Stimulus Features Encoded by vS1 and vM1

(A) Percent of vS1 single and multiunits neurons with significant coding, as a function of time, of sp1 (green), sp2 (red), T1 (blue), and difference between sp2 and sp1 (black). Data points are plotted at the onset of a 300-ms sliding window. In this dataset, T2 was fixed (400 ms) and not treated as a task variable.

(B) The same as (A) for vM1 single units.

In vS1, 30%–40% of neurons significantly (p < 0.01; permutation test) encoded sp1 and sp2 during presentation of stimulus 1 and 2. A small number of vS1 neurons (under 10%) continued to encode sp1 up to 800 ms after its termination. During the interstimulus delay, neurons in vS1 did not encode the duration of the just-completed stimulus, T1. Coding of the difference between sp1 and sp2 is shown by the black trace. Stimulus-difference-coding neurons, which may have a role in comparing stimulus 1 and stimulus 2 in order to reach a decision [4], are identified by virtue of significantly encoding both sp1 and sp2 but with opposite sign. A very small set of vS1 neurons showed this property, perhaps reflecting frontal cortex feedback to vS1 [15].

In vM1, about 30% of neurons encoded sp1 during presentation of stimulus 1 and about 20% encoded sp2 during presentation of stimulus 2. In contrast to vS1, many neurons continued to show sp1-dependent firing throughout the interstimulus interval, signifying the participation of vM1 in working memory. Also, in contrast to vS1, T1 coding remained robust during the interstimulus delay. Finally, a much larger set of vM1 neurons encoded the difference between sp1 and sp2 (black trace) and thus may participate in decision making.

A correlated response variation between stimulus speed and stimulus duration—that is, increased firing rate both for increased sp and also for increased T—could contribute to a longer stimulus being perceived as stronger. Such correlated speed/duration effects were not common in vS1. An example neuron from vS1 is given in Figure 7A. Neuronal firing in a 300-ms window centered on the end of stimulus 1 is sorted according to T1. Higher values of sp1 evoked greater firing rate, but stimulus 1 of 600 ms duration (cyan) evoked a lower firing rate than did stimulus 1 of 200 ms (blue). This might be an outcome of adaptation, whereby neurons’ firing rates decrease as vibrations are extended in time [16]. Of the 41 neurons in vS1 that exhibited a statistically significant dependence of firing rate on sp as well as on T1, only 19% exhibited increased firing rate for the 600- versus the 200-ms stimulus (Figure 7A, pie chart inset).

Figure 7

Effect of Stimulus Duration on vS1 and vM1 Neurons

(A) vS1 single-unit response as a function of sp1 in the 300-ms window centered on stimulus offset. Responses are separated according to T1: 600 ms (cyan) and 200 ms (blue). The pie chart shows the proportion of neurons, among those with significant variation in firing as a function of sp, with significant positive (600 > 200; filled) together with negative (600 < 200; unfilled) firing rate correlation with T1. The measurement window for the pie chart includes the entire interstimulus interval.

(B) vM1 single-unit response 400–700 ms into the interstimulus delay, using the same conventions as (A). Contrary to the neuron in vS1, firing rate was higher when T1 was longer. The pie chart shows the proportion of neurons, among those with significant variation in firing as a function of sp, with significant positive (600 > 200; filled) and negative (600 < 200; unfilled) firing rate correlation with T1.

(C) Left: the projection of the vM1 population over time, optimized in each time window to represent differences in sp1 (see color bar above the panel). Right: the length of population projections for the time window 100–400 ms into the delay, with projections separated according to T1: 200 ms (blue) or else 600 ms (cyan). Greater T1 led to longer projection, mimicking the effect of higher sp1. Solid lines are the SEWP model fit on population response.

(D) Difference in population projection length for trials with T1 = 600 ms versus T1 = 200 ms. vS1 projection difference is gray; vM1 difference is green. Shaded area represents one SD over 500 training/test repetitions.

(E) Four proposed models for weighting vS1 activity, from upper left to lower right: by firing rate with a uniform weighting function (“firing rate”); by summation of single spikes with a uniform weighting function (“spike count”); by summation of single spikes with an exponential weighting function anchored to stimulus offset (“recency spike count”); and by summation of single spikes with an exponential weighting function anchored to stimulus onset (“primacy spike count”). Solid lines are the neurometric curves that would result from decoding the vS1 population according to the corresponding model.

(F) SEWP prediction of vM1 population projection length (abscissa) versus the observed projection length of the population (ordinate) in different epochs of the trial.

A vM1 single unit expressing the correlated speed/duration effect is illustrated in Figure 7B. To uncover the effect of stimulus duration, we sorted neuronal firing in a window 400–700 ms after stimulus offset according to T1 (the same effect was seen in other windows). Higher values of sp1 evoked greater firing rate, but, in contrast to the vS1 neuron, stimulus 1 of 600 ms duration (cyan) evoked a greater firing rate than did stimulus 1 of 200 ms (blue). In general, it is this dependence of firing rate on duration that leads to a significant value of the T1 term in the linear model of Figure 6B. Following from this neuron’s temporal integration, a firing rate of 10 spikes/s could signify the occurrence of a sp1 of 80 mm/s for 600 ms or else 100 mm/s for 200 ms. Of the 38 neurons in vM1 that exhibited a statistically significant dependence of firing rate on sp1 as well as on T1, 79% exhibited increased firing rate for the 600- versus the 200-ms stimulus (Figure 7B, pie chart inset).

Having found vM1 single-neuron properties consistent with temporal integration, we extended the analysis to the population level. The data in Figure 7C were derived from 50 single units recorded in four rats. We represented population dynamics as high-dimensional trajectories in neuronal activity space [17] and identified the angle that maximizes variation in projection length as a function of sp1; projection length is then an estimate of population mean speed signal. Because the analysis method does not consider within-trial correlations among neurons, it was possible to pool data across 16 sessions. The left plot shows projection length, with trials grouped according to sp1 (dark red to yellow). T1 was either 200 or 600 ms; to combine all trials, the plot is aligned (set to 0 ms) at stimulus onset on the left and at stimulus offset on the right. The angle was optimized for each sliding 300-ms window. The main finding is that projection length increased as sp1 increased, and this population code was retained across the interstimulus delay.

To discern the role of temporal integration in the same cortical population, the right plot of Figure 7C illustrates projection length measured 100–400 ms after stimulus offset as a function of sp1 but now with 200- and 600-ms trials separated (blue and cyan, respectively). For both durations, projection length still increased with sp1 but was greater after the 600-ms stimulus than after the 200-ms stimulus. Thus, if the neuronal population were decoded according to projection length, a stimulus of a given sp would be perceived as more intense as duration increased. Figure 7D plots the difference in projection length between 200- and 600-ms stimuli for vM1 (green) and vS1 (gray; 52 units in five sessions) populations, over time. Stimulus offset is aligned to 0 ms. A projection difference of 0 is denoted by the red dashed line. In both vM1 and vS1, 600-ms projections were smaller than 200-ms projections in the final window of the stimulus presentation period, reflecting adaptation. During the interstimulus interval, the vS1 population showed equal projection length for 200- and 600-ms stimuli (projection difference = 0), indicating that the preceding T1 no longer affected the state of the population. By contrast, the vM1 population (even more so than vM1 single units; Figure 6B) showed greater projection length after 600-ms stimuli and maintained this difference across the entire interstimulus delay, providing a mechanistic account for the duration-dependent behavioral bias illustrated in Figures 2C–2F.

Figure 7E summarizes the relation of vS1 to the intensity percept. The points with error bars (identical across the four plots) indicate the choices made by two rats across five sessions. T2 was of either 200- or 600-ms (blue and cyan), whereas T1 was fixed at 400 ms. The activity of a population of 115 vS1 neurons is read out according to four posited integration schemes, each scheme providing predicted choices for the same set of trials (see STAR Methods). Predicted choices, in the form of neurometric curves, are shown by the solid lines.

The first candidate readout mechanism, mean firing rate, applies a uniform weighting function to vS1 firing. If rats’ choices are fit, temporal integration might already take place within vS1. Whereas the population neurometric curves replicate the rats’ capacity to judge the relative values of sp1 and sp2, the curves do not reflect their tendency to judge the longer stimulus 2 as more intense.

The next candidate readout mechanism, spike count, applies a uniform weighting function to vS1 firing but summates the spikes without division by duration. If rats’ choices are fit, then the primacy of the SEWP model might simply arise from the greater number of spikes at stimulus onset. But whereas the neurometric curves do predict the rats’ judgement of the longer-lasting stimulus 2 as more intense (because more spikes accumulate over 600 ms), they fail to replicate the comparison between sp1 and sp2.

The final two candidate readout mechanisms accumulate the spikes emitted by vS1 neurons after applying exponential weighting functions, either giving most weight to the response adjacent to the stimulus offset (labeled “recency spike count”) or adjacent to stimulus onset (“primacy spike count”). The neurometric curves based on recency capture both the rats’ tendency to judge the longer-lasting stimulus 2 as more intense as well as their capacity to compare the relative values of sp1 and sp2. The primacy-based neurometric functions even more closely fit the rats’ choices, capturing the duration dependence while conserving their sp dependence. These analyses, taken together, suggest that the behavioral phenomenon of SEWP does not originate in vS1, for the firing of this cortical region must be further weighted by primacy in order to match the observed behavior.

If vM1 contributes to the duration-dependent intensity percept, two properties would be expected: for some fixed duration, the neuronal population should encode sp in a linear manner, and that line should shift upward or downward in relation to stimulus duration. Whereas Figure 7B shows just such properties for a single neuron, we now ask whether the size of the shift in vM1 matches the prediction of the SEWP model. In Figure 7F, the response of vM1 is plotted as a population projection in sessions where stimulus 1 was of either 200 or 600 ms duration (same dataset as in Figures 7C and 7D). The color of the points represents sp, increasing from dark red to yellow. The vM1 firing is fit according to a function based on SEWP to give a predicted population projection length (abscissa). This is compared to the observed projection length of the population (ordinate). The model-experiment comparison is carried out across four intervals of the trial, denoted on each plot. The diagonal line corresponds to perfect alignment between the model and the observed data. In all plots, greater sp leads to a greater population projection, reflecting the predicted and observed stimulus coding in vM1 during the stimulus 1 and throughout the interval preceding stimulus 2. During the final 200 ms of stimulus 1, the observed vM1 population projection differs from the model prediction due to pronounced adaptation during the longer stimulus; no adaptation mechanism is incorporated in the model. At each of the other trial intervals, the observed vM1 data closely track the model prediction. These findings indicate that vM1 neuronal populations, but not vS1 populations, are aligned with the behavioral SEWP model and thus could be a component of the neuronal substrate for the vibration intensity percept.

Discussion

Stimulus Integration Causes a Non-stationary Intensity Percept

Psychophysical experiments revealed a fundamental mechanism for the integration of noisy sensory information over time, a mechanism that led to a distinction between the ongoing values of speed—which were, on average, stationary—and the final perceived intensity—which grew as a function of stimulus duration. In sum, longer stimuli were felt as stronger. The integration function was not a product of learning, for it did not depend on the training protocol. Time dependence has been reported in monkeys in loudness judgement [18] and tactile flutter frequency discrimination [10] and in humans for luminance and loudness judgment [19, 20]. To this previous work, we have added a simple model, SEWP, which accounted for psychophysical results across widely ranging stimulus conditions. A further novel finding is that humans and rats share this basic integrative mechanism: they differ in the overall acuity parameter (humans are more accurate) but are similar in , the time constant for integration.

When a judgement must be made about a single noisy stimulus, evidence can be accumulated over time [21, 22]. But with such single-stimulus judgements, the processes of sensory integration and decision making overlap in time as the brain moves toward a choice concurrent with stimulus presentation. In the present work, by virtue of uncoupling sensory integration from decision making and motor execution (i.e., no decision making operation can be initiated or executed during stimulus 1), it becomes clear that temporal integration occurs in the sensory-perceptual phase of cognition (Figure 7C).

The Neuronal Correlate of Perceived Intensity Resides in vM1 and Not in vS1

In tasks where rodents discriminate sequences of discrete whisker deflections, vS1 neurons fail to show substantial temporal integration [23, 24]. Similarly, in the present work, vS1 responses were characterized by a local timescale (Figures 5E and 5F). In vM1, a principal vS1 target in frontal cortex [13, 14], neuronal responses failed to show robust local coding; instead, they varied in relation to stimulus mean speed as well as in relation to stimulus duration. Firing rates in vM1 population did not “ramp up” along the course of stimulus 1 (Figure 4B), yet the population assumed a duration-dependent state immediately after termination of the vibration and maintained this state throughout the interstimulus delay (Figures 7C and 7D). During stimulus 2, vM1 neurons encoded the difference between sp1 and sp2, suggesting involvement in the operation of comparison (Figure 6B). Thus, vM1 participates in multiple steps of the task, similar to ventral premotor cortex and dorsolateral prefrontal cortex in monkeys [4, 25] and consistent with a prominent role in action selection [26, 27]. Brain regions we have not studied here, such as secondary vibrissal sensory cortex (vS2), posterior parietal cortex, and other frontal cortical areas may also play important roles, as in monkeys [4].

By what mechanism is the intensity judgment of humans and rats confounded by stimulus duration? The properties of vM1 single units and populations provide, for the first time, a possible answer: the effects of increasing stimulus speed and increasing stimulus duration were of the same sign. Due to the embedded, correlated coding of both properties, any given state of the population could be decoded as one of a (theoretically) infinite number of conjunctions of speed/duration. Thus, the representation of stimulus speed in vM1 (approximated by projection length in Figure 7) cannot exclude the influence of duration, even if the behavioral rule excluded duration as a factor.

The present work is consistent with the notion that the posterior-to-anterior transformation of information in the cerebral cortex is accompanied by a broadening of the timescale by which neurons encode sensory events [28, 29]. In behavioral tasks that require identification of instantaneous events, the critical stimulus features may be decoded from early sensory processing areas. In behavioral tasks that require extraction of longer-timescale stimulus properties, like the intensity of a stochastic vibration in the current work, the critical stimulus features may be decoded from higher-order, mixed-function processing areas, such as vM1.

STAR★Methods

Key Resources Table

Contact for Reagent and Resource Sharing

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Mathew E. Diamond (diamond@sissa.it).

Experimental Model and Subject Details

Twenty male Wistar rats (Harlan, San Pietro al Natisone, Italy) were housed individually or with one cage mate and maintained on a 14/10 light/dark cycle. Daily access to water was restricted to promote motivation in the behavioral task, yet weight gain followed a standard Wistar-specific curve, indicating that the quantity of water obtained during training and testing was comparable to the ad lib quantity. After each session, rats were placed for several hours in a large, multistory enriched environment with other rats. Twenty healthy human subjects (6 males and 14 females, ages 22–35) were tested after giving their informed consent. Protocols conformed to international norms and were approved by the Ethics Committee of SISSA and, for rat procedures, by the Italian Health Ministry (license numbers 569/2015-PR and 570/2015-PR).

Method Details

Behavioral task

Each trial began when the rat positioned its nose in the nose-poke (equipped with optic sensor) and placed its whiskers on the plate (Figure 1A). After a short delay (800 ms), stimulus 1 was presented, characterized by nominal mean speed, sp1, and duration, T1. After the interstimulus delay of 2 s, stimulus 2 (sp2 and T2) was presented (Figure 1B). The rat remained in the nose-poke throughout both stimuli and could withdraw only when the “go” cue sounded to terminate the poststimulus delay of 600 ms. Early withdrawal was considered an aborted trial and went unrewarded. At the go cue, the rat selected the left or right spout; reward location depended on the relative values of sp1 and sp2. Incorrect choices went unrewarded. Trials with sp1 = sp2 were rewarded randomly.

Rats learned the delayed comparison task by generalizing the comparison rule across the entire stimulus range. If sp1 were fixed across all trials and only sp2 shifted, rats might solve the task by ignoring stimulus 1and applying a constant threshold to stimulus 2. Likewise, if sp2 were fixed across all trials, rats might apply a constant threshold to the stimulus 1. To avoid such shortcut strategies, we used the stimulus generalization matrix (SGM; Figure 1C), whereby neither stimulus alone provided the information necessary for a correct choice [6].

Test sessions began after rats reached an overall performance of more than 75% correct on 10 successive training sessions of at least 200 trials each. During test sessions, the SGM included stimulus pairs with normalized speed difference of 0.3 and −0.3, slightly larger than the maximum NSD of 0.2 and −0.2 used in testing human subjects. This was done in order to provide easier trials to maintain the motivation of rats on the difficult psychometric stimulus pairs (green rectangles, Figure 1D).

During training, both T1 and T2 were 400 ms (“balanced” durations) in all except 3 rats; these received multiple combinations of T1 and T2 (“unbalanced” durations) from the earliest stage of training. In test sessions, the increase in steepness of the psychometric curves with longer duration in the balanced condition was quantitatively indistinguishable in 8 rats trained with T1 and T2 fixed at 400 ms versus 3 rats trained with variable T1 and T2 (comparison of change in maximum curve slope between the two groups, t test: from 100/100 to 600/600 ms, p = 0.59, permutation test). Likewise, the leftward and rightward curve shifts in the unbalanced duration were indistinguishable in the balanced-trained and unbalanced-trained rats (comparison of change in inflection point between the two groups, t test: from 200/400 to 400/400, p = 0.67; from 400/400 to 400/600, p = 0.91, permutation test). In sum, neither the improved performance with increasing duration nor the tendency to judge longer stimuli as stronger depended on training regime.

Both human and rat experiments were controlled using LabVIEW software (National Instruments, Austin, Texas). Vibrations were generated by stringing together sequential velocity values (v) at 10,000 samples/s, taken from a normal distribution. In all analyses we treat the stimuli as a sequence of discrete samples, although the motor moved through space continuously. The velocity time series for a given trial was taken randomly from among 50 unique sequences of pseudo-random values except in recording sessions where we wished to measure neuronal response to repeated, “frozen” noise (Figures 5A and 5B). Because stimuli were built by sampling a normal distribution, the statistical properties of an individual vibration did not perfectly replicate those of the distribution from which it was constructed. Converting v to its absolute value, sp, the distribution takes the form of a folded, half-Gaussian (see Figure 1B, to the right of vibration traces). Stimuli delivered to human subjects on the fingertip were the same as those used in rats except that the velocity values were halved. Human subjects received feedback (correct/incorrect) on each trial through a computer monitor.

Analysis of behavioral data

To characterize the behavior, we computed the proportion of trials in which subjects judged stimulus 2 > stimulus 1 on closely spaced pairs (Figure 1C, pairs within green rectangles). We fit the data with a four-parameter logistic function using the nonlinear least-squares fit in MATLAB (MathWorks, Natick, MA):where NSD is normalized speed difference, (sp2 − sp1) / (sp2 + sp1), is the lower asymptote, is the upper asymptote, 1/ is the maximum slope of the curve and is the NSD at the curve’s inflection point.

Behavioral model for temporal integration of speed

To characterize working memory performance, we computed the proportion of trials in which subjects judged stimulus 2 > stimulus 1 as a function of sp2 and sp1 [6]. We fit the data with a generalized linear model in MATLAB, as follows:Where w1 is the sp1 regressor, w2 is the sp2 regressor, and w is the baseline regressor that captures the overall (stimulus-independent) bias of the subject in judging stimulus 2 > stimulus 1. γ and λ are the lower and upper asymptotes, respectively. These regressors were then used to normalize the relative impact of each stimulus on individual subject’s choice.

Starting with the simple idea that the subject integrates both stimuli and then compares them, choice can be modeled as follows. For the general case when the degrees of noise (variability) of two signals are not necessarily equal, discriminability can be quantified as:

At this stage, we substitute parameters of the standard by our model of the subject’s percept (Figure 3A). Given a half-Gaussian distribution of speed values in each trial, we hypothesize that subjects accumulated instantaneous speed sp over time, after applying a weighting function:andwhere is the weight at time t (in ms after stimulus onset) and T is the total stimulus duration. Next, we calculate the expected value and variance of the perceived stimulus, as follows:

The expected value of a sum is equal to the sum of the expected values. Therefore,

Since  =  = sp, then:

The solution of the sum is:

Multiplying both sides by we get:and

Therefore, we can rewrite the expected value of perceived intensity as:

Next, we calculated the variance of perceived intensity:

One property of variance is that

As covariance between speed samples is zero (due to the fact that samples are randomly taken from a half-normal distribution), we can express the variance as:

We then rewrite the variance as follows:

As variance of speed is equal to we can again rewrite it as:

Finally we can rewrite as: is converted to choice by a sigmoid function similar to that used in behavioral analysis,where is taken as a mathematical expression of the stimulus sensitivity carried by brain signals, γ is the lower asymptote, λ is the upper asymptote, 1/ is the maximum slope of the curve and μ is overall bias. We then split the data for each subject into training and testing set (50/50 percent).

The training trials were used to estimate (nonlinear least-squares, MATLAB) the values of the parameters of the model for each individual subject. The percent of test trials judged stimulus 2 > stimulus 1 was then predicted by the model and compared to the observed choice. This procedure was repeated 200 times.

Electrode implantation and data acquisition

Rats (n = 6) were anesthetized with 2%–2.5% Isoflurane delivered with oxygen under controlled pressure through a plastic snout mask. They received an implant either in vS1 (n = 1), in vM1 (n = 4), or in both (n = 1). The two target regions (Figure S4) were accessed by craniotomy, using standard stereotaxic technique (vS1: centered 2.8 mm posterior to bregma and 5.8 mm lateral to the midline, vM1: centered 2.5 mm anterior to bregma and 1 mm lateral to the midline). Dura mater was removed over the entire craniotomy of both areas with a small syringe needle. The remaining pia mater, even if usually not considered to be resistant to penetration, nevertheless presents a barrier to the entry of the microelectrode arrays. This resistance leads to dimpling of the brain at the moment of penetration. To minimize such dimpling, after applying the petroleum based antibiotic (Bimixin, Sanofi-Aventis) to the center of craniotomy, we applied cyanoacrylate adhesive directly to the pial surface bordering the edge of the cranial opening. This procedure fastens the pia mater to the overlying bone and the resulting surface tension prevents the brain from compressing under the advancing electrodes. The electrode arrays (Tucker-Davis Technologies (TDT), Alachua FL) were configured as 16 electrodes (2 rows of 8, with 250 μm within-row spacing and 375 μm between-row spacing) or 32 electrodes (4 rows of 8, equivalent spacing). Electrodes were sharply tapered and shaft diameter was 50 μm. The electrode array was inserted by slowly advancing a Narashige micromanipulator. After inserting the array(s), the remaining exposed cortex was covered with biocompatible silicon (KwikSil; World Precision Instruments). Two rats were implanted with a novel movable electrode array (SISSA microdrive, TDT) only in vM1. In all rats, 5 small screws were fixed in the skull as support for dental cement. One of the screws served as a ground electrode.

One hour after the beginning of the anesthesia, atropine (2 mg/kg) was injected (s.c.) to avoid secretions in the respiratory tract and maintain a stable heart rate. They were given the antibiotic (Baytril; 5 mg/kg; i.p.) and analgesic (Rimadyl; 2.5 mg/kg; i.m.) one hour before conclusion of the operation. After surgery, a local antibiotic (Isaderm) was applied around the wound to help the healing. In addition, both the antibiotic and the analgesic were delivered through the water bottle for 24 hr after completion of surgery. During this recovery time, rats had unlimited access to water and food. Recording sessions in the apparatus began thereafter.

Extracellular activity of both vS1 and vM1 was sorted into single units and multiunits, as verified through the spike waveform and the refractory period observed in interspike interval histogram using a MATLAB-based software, UltraMegaSort 2000 [30]. In total 90 single- and 300 multiunits were recorded in vS1 of 2 rats and 217 single- and 398 multiunits were recorded in vM1 of 5 rats. Two criteria had to be met for a neuron to be included in the analyses. First, overall firing rate within the session was at least 2 Hz. Second, overall firing rate per trial for the entire session, with many different stimulus conditions intermixed across trials, did not show a significant non-zero linear correlation (p < 0.05) over time. The second criterion was aimed at excluding unstable recordings. In the analyses, we were able to pool the data from vS1 single- and multiunits whenever appropriate, due to the similar coding properties indentified among single units. We did not pool vM1 single- and multiunits because of coding heterogeneity.

Measure of local coding

To measure to what extent the neuronal response varied according to the local speed events in the stimulus waveform, we computed the cross-correlation between the vibration waveform, sp, and the spike train (Figure S6). The complete speed waveform for one vibration is denoted SPW. When the spike train temporal resolution is sufficiently fine (in this analysis 0.1 ms), each bin can contain a single spike (assigned the value 1) or else no spikes (assigned the value 0). The spike train during the first stimulus of the first trial can be depicted as a vector of 0’s and 1’s, RV1, with size s of T1 (ms) × 10. The spike train accompanying stimulus 1 of trials 1 to n was concatenated to produce a long response vector, RVstimulus 1 = (RV1,RV2,RV3,…,RVn). The size of the vector RVstimulus 1 is the sum of all T1s. After carrying out the same procedure for all presentations of stimulus 2, RVstimulus 1 and RVstimulus 2 were combined to produce a long vector of spikes trains, RV RV was convolved with a Gaussian Kernel, K(t,σ) to generate a spike density function SDF. The vibration waveform for the first trial and stimulus, SPW1 was concatenated with the vibration waveforms of other trials (replicating the spike train procedure) to produce a long vector SPWstimulus 1 = (SPW1, SPW2, SPW3, …, SPWn). SPWstimulus 1 and SPWstimulus 2 were then combined, preserving the stimulus sequence corresponding to the concatenated spike train.

The cross-correlation between SPW and RV indicates by how many ms RV must be shifted to make it identical to SPW. When the two signals match, then the quantity of cross-correlation is maximum. The corresponding shift can be considered the time lag between the two signals. As SPW contains different trials (therefore different global sp values), the peak of the correlation might be related to trial to trial changes in sp and not to local sp. As a control, we computed the same analysis for temporally reversed waveforms (each individual stimulus waveform reversed before concatenation), where the local correlation must be negligible, by definition. The width of the Gaussian kernel, K(t,σ) was varied and, for a given Gaussian kernel width, forward and reversed cross-correlograms were compared and standardized by subtracting the mean and dividing by the standard deviation of the whole vector. When the Gaussian kernel width produced maximal cross-correlogram difference, then the size of the difference as well as the kernel width were defined as measures of local coding.

Linear model and neuronal correlates of task events

All analyses and statistical tests were done with custom codes written in MATLAB. To specify how the neuronal populations encoded relevant aspects of the behavioral task, we quantified the relations between neurons’ (single- and multi-units) firing rate and the variables of interest – mean speed of both stimuli (sp1 and sp2) and durations of both stimuli (T1 and T2) by using a multi-variable, linear regression,where rate(t) is the firing rate of each neuron measured in a 300-ms time bin, sliding in 5-ms steps, expressed as a z-score with respect to whole-session mean firing rate at each time bin, and are the sp1 and sp2 regressors, and are the T1 and T2; all regressors are fit as a function of time t. is a constant set according to the neuron’s overall excitability at each time time t, independently of task variables. Values of the regressors indicate how strongly, at any given moment, the neuron’s response varied in relation to the associated variable. In some sessions, T2 was fixed (400 ms) and not treated as a variable; we can therefore take to be 0 while in other sesssions was set to 0 as T1 was kept fixed. The test for considering a neuron to express significant sp1, sp2, T1 or T2 coding was to have significant linear coefficients (t test, p < 0.05) in 6 out of 10 sequential windows. The number of neurons with significant coding was then normalized by subtracting the number of neurons expected by chance to show significant linear coefficient (criterion of chance level: 4 standard deviations from the mean of a shuffled distribution for each bin).

Neurons that encoded both sp1 and sp2 with opposite sign but similar absolute magnitude, would pro-duce a differential signal, a signal that could reflect a comparison process [9]. Similar criteria to those described above were used to define the significance of each neuron’s involvement in the stimulus com-parison process (Figure 6).

Analysis of cortical population response

The analysis of population response focused on stimulus 1 so that coding of sp1 and T1 could be examined without mixed signals related to the stimulus comparison, choice, and action. Following the methods for demixed dimension reduction [17], the normalized firing rate per time bin was averaged over each condition (i.e., sp1 and T1) and expressed as a z-score with respect to the average rate in that time bin. The time bin was 300 ms and was shifted in steps of 5 ms.

The size of the population averaged matrix, PopM, is the numberofneurons × conditions × timebins. We performed a principal components analysis (PCA) to the PopM to find the dimensions that best explained the variability of neuronal responses. We used the first n principal components that yielded at least 75% of the explained variability as a “de-noising” matrix, D. We then used linear regression to define the stimulus- and task-relevant axes of sp1 and T1 and obtained the de-mixed population trajectory by projecting the population activity onto these task-related dimensions. It is important to take into account that the linear regression matrix used to obtain the task-related dimension is different for each time bin, so that the optimal population projection is likely to be different for each time bin. The dataset was divided into the training set and test set. The population trajectories obtained from the test data -set are depicted in Figure 7. Similar linear regression as Equation 16 was used to define stimulus- and task-relevant axes of sp1 and T1. The de-noising matrix was obtained by combining the first n principal components [17],

We then obtained the de-mixed population trajectory by projecting the population activity onto the task related dimensions. The de-noised population activity was obtained by multiplying the population averaged response to the de-noised linear coefficients, Wdenoised.where Wdenoised (t) is the denoised regression weights matrix (size of N neurons × timebins) and D is the “de-noising” matrix and Wsp1 (t).

The final step was to apply the de-noised linear coefficient matrix (Wdenoised) to the population averaged matrix PopM to obtain the task-relevant dimension:

Prediction of population activity by SEWP model

Analysis of vS1 readout

The analysis of vS1 population coding was restricted to neurons defined as responsive to stimulus presentation compared to the control period (200 ms before the stimulus onset, t test, p < 0.05). On each trial, the response of each neuron during the presentation of stimulus 1 and stimulus 2 was computed based on 4 different readout mechanisms (Figure 7E). According to the relative magnitudes of these responses, the neuron was considered to have reported stimulus 2 > stimulus 1 or else stimulus 1 > stimulus 2. For the full population, the decoded choice on that trial corresponded to the choice of the majority of neurons weighted by their Pearson correlating coefficient. To generate neurometric curves, trials were combined according to the two vibrations’ relative mean speeds and the percent of trials classified as stimulus 2 > stimulus 1 was fit with a sigmoid function to generate neurometric curves (Figure 7E):where NSD is normalized speed difference, (sp2 − sp1) / (sp2 + sp1), γ is the lower asymptote, λ is the upper asymptote, 1/ν is the maximum slope of the curve and μ is the NSD at the curve’s inflection point.

Author Contributions

Conceptualization, A.F., A.A., V.S., and M.E.D.; Methodology – Formal Analysis, A.F., A.A., V.S., and M.E.D.; Investigation – Animal Subject Training, Testing, Neurophysiology, Histology, A.F. and F.P.; Investigation – Human Subject Training & Testing, A.F.; Writing – Original Draft, A.F. and M.E.D.; Writing – Review & Editing, A.F., A.A., F.P., V.S., and M.E.D.; Funding Acquisition, M.E.D.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Chemicals, Peptides, and Recombinant Proteins
Isoflurane	Merial	AP/DRUGS/220/96
Epigel	Ceva	N/A
Atropine	ATI	AIC no. 101948014
Antibiotic (Baytril)	Bayer	AIC no. 100155062
Analgesic (Rimadyl)	Zoetis	AIC no. 102191119
Topic antibiotic	Dechra	AIC no. 102881012
Local anesthetic (lidocaine)	Molteni Farmaceutici	AIC no. 005638010
Paraformaldehyde	Sigma-Aldrich	158127
Cresyl Violet	Sigma-Aldrich	C5042
Experimental Models: Organisms/Strains
Wistar rats	Harlan - Envigo	16808M
Software and Algorithms
LabView 2014	National Instruments	http://www.ni.com/download/labview-development-system-2014/4735/en/
MATLAB v 2015b	MathWorks	https://www.mathworks.com/products/matlab/
UltraMegaSort 2000 for spike sorting	By DN Hill, David Kleinfeld	http://neurophysics.ucsd.edu/software.php
Other
Digital TDT recording system	Tucker-Davis Technologies (TDT)	model: RZ2 BioAmp Processor
Hyperdrive	Custom built; Designed by SISSA Mechatronics Lab and Tucker-Davis Technologies	http://www.tdt.com/tdt-microdrive.html