Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Neurosci. Author manuscript; available in PMC 2011 January 1.
Published in final edited form as:
PMCID: PMC2957657

Encoding and decoding cortical representations of tactile features in the vibrissa system


During behavior, rats and other rodents use their facial vibrissae to actively explore surfaces through whisking and head/body movement, resulting in complex sensory inputs that vary over a large range of angular velocities and temporal scales. How these complex sensory inputs manifest in the patterns of cortical firing events that ultimately form the perceptual experience is not well understood. Through single-unit cortical recordings of layer-4 neurons in S1 of the anesthetized rat, we systematically quantified the interactions between instantaneous velocity and timing of vibrissa motion, finding a strong interaction between angular velocity and timing of contacts on the 10s of milliseconds time scale. From the quantification of these joint tuning properties, a detailed nonlinear encoding model was formulated that was highly predictive of firing probability and timing characteristics of the sparse cortical representation of complex patterned tactile inputs. Within a Bayesian framework, the encoding model was then used to decode tactile patterns under simple transformations of the stimulus along dimensions of velocity and timing, as a demonstration of the lower-bound of the idealized perceptual capabilities of the animal.

Keywords: vibrissa, velocity, decoding, barrels


In the rodent vibrissa system, an array of facial vibrissae (whiskers) is endowed with acute tactile sensing mechanisms that serve as an important sensory modality in locomotion, balance, texture sensation, and wall-following (Vincent, 1912). The mechanical signals transduced by mechanoreceptors embedded in the soft tissue surrounding the vibrissa follicles are strongly shaped by the geometry of vibrissa-surface contact, the mechanics of the vibrissa and the follicle complex, and the kinematics of body motion and whisking. Recent behavioral studies have demonstrated that in ethologically relevant conditions in which the animals were actively moving their vibrissae across textures, the vibrissae undergo patterns of discrete high-velocity transients (Ritt et al., 2008; Wolfe et al., 2008), where the velocity and timing of the transients depend on the textural properties and the precise nature of the behavior. Only recently are the effects of variations in vibrissa input on cortical representations being explored in the context of the textural environment (von Heimendahl et al., 2007; Jadhav et al., 2009).

Studies of cortical activity in response to electrically-induced whisking across textured surfaces in anesthetized rats have revealed temporally precise, stimulus locked spiking activity, described as a “kinetic signature” of neural activity defined by the temporal pattern of vibrissa deflections (Arabzadeh et al., 2005, 2006). In acute studies, the mean firing rate of neurons in the vibrissa region of the primary somatosensory cortex has been shown to depend strongly upon the velocity of isolated, punctate vibrissa deflections (Ito, 1985; Simons, 1978; Pinto et al., 2000). Impinging suppressive influences on thalamocortical neurons further complicate this selectivity by inducing a highly nonlinear dependence of single cell responses on stimulus history for more complex stimuli, likely linked to the interplay between excitatory and inhibitory elements of the neural circuitry (Simons and Carvell, 1989; Kyriazi et al., 1994; Nicolelis and Fanselow, 2002; Higley and Contreras, 2006; Webber and Stanley, 2004; Boloori and Stanley, 2006; Higley and Contreras, 2007; Heiss et al., 2008). The functional consequence of the interaction between velocity sensitivity and nonlinear dynamics of the excitatory/inhibitory interplay has not been investigated, yet is a key to understanding the dynamics of encoding in complex, natural environments.

Here, in the vibrissa system we explore the cortical representation along two primary dimensions: the moment-by-moment angular velocity and the timing of transient events within the stimulus pattern. Through single-unit extracellular cortical recordings in the anesthetized rat, the effects of angular velocity on cortical firing probability and timing were systematically quantified using simplified stimuli, and the tuning properties were utilized to formulate a nonlinear encoding model that was highly predictive of cortical firing activity in response to the more complex tactile patterns, capturing the spike count, latency, and precision of firing in localized firing events. The encoding model was then utilized for decoding in a Bayesian framework, showing through systematic transformation of the stimulus patterns along dimensions of velocity and timing that the encoding model captures a range of dynamic properties of the pathway.

Materials and Methods

Surgical Preparation

For all acute experiments, all procedures were approved by the Animal Care and Use Committee at Harvard University, and are in agreement with guidelines established by the National Institutes of Health. See Webber and Stanley (2004) for a detailed description of experimental methods. Briefly, female albino rats (Sprague-Dawley; 250-330g) were sedated with 2% vaporized isoflurane, and anesthetized with sodium pentobarbital (50 mg/kg i.p., initial dose); supplemental doses of approximately 12 mg/kg/hr (Zhu and Connors, 1999) were given as needed to maintain a light level of anesthesia – confirmed by measurements of heart rate, respiration and eyelid/pedal reflexes to averse stimuli (toe or tail pinch). Following initial anesthesia, the animal was mounted on a stereotaxic device (Kopf Instruments, Tujunga, CA) in preparation for surgery and subsequent recordings. After the initial midline incision on the head, tissue and skin were resected and connective tissue was removed. A small craniotomy (approximately 2 mm in diameter) was made over the barrel field (stereotaxic coordinates: 1.0-4.0 mm caudal to the bregma, and 3.0-7.0 mm lateral to the midline; Paxinos and Watson, 1998). The dura layer was left intact for recordings with single tungsten micro-electrodes. For recordings made by tetrodes (see below), the dura layer was carefully reflected. In all cases a dam was created with bone wax around the craniotomy and the cortex was covered with mineral oil solution to stabilize the recording and to provide additional protection for the the cortical surface. Following the recording session, the animal was euthanized with an overdose of sodium pentobarbital solution (0.5 mL of Euthasol, Henry Schein, Inc., Melville, NY USA).


Single unit extracellular recordings were obtained using both single tungsten micro-electrodes (impedance range: 3-7 MΩ; FHC Inc., Bowdoinham, ME USA), and tetrodes (impedance range: 1-2 MΩ; Acute Probe Model: a 2×2-tet, NeuroNexus Technologies, Ann Arbor, MI USA). Data were collected using a 16-channel Plexon data-acquisition system (Plexon Inc., Austin, TX USA). Neuronal signals were filtered, amplified, and digitized at 40 kHz/channel. Both the tetrode and single-electrode recordings were analyzed using the OfflineSorter software suite (Plexon Inc.) to assign the recorded spike waveforms to single-units on the basis of standard template matching techniques and physiologically plausible refractory periods (Lewicki, 1998).

All recorded cells were located at stereotaxic depths of 800-1200 μm, commensurate with cortical layers IV-V. For all recorded neurons, spontaneous activity was collected prior to the beginning of the stimulus trials, and again after all stimuli were presented. Cortical receptive fields were determined manually by identifying all vibrissae that elicited spiking responses for each recorded cell. The primary vibrissa (PV) was identified as the whisker generating maximal responses.

Only cells classifiable as regular-spike units (RSUs) were reported in this study (n = 44 cells). Based on latency analyses (PV: 9.5 ms ± 0.6 ms), all included cells were determined to reside in monosynaptic thalamocortical recipient zones. The additional criteria for their classification as SI RSUs were: (1) stability of baseline responses to PV deflections over the duration of the experiment, used to ensure consistency over the long (approximately 3 hours) recording interval, (2) width of action potentials (typically 1.5 ms) (McCormick et al., 1985), used to exclude putative inhibitory units (i.e the fast spiking units, FSUs), (3) extracellular responses to at most 1-2 adjacent vibrissae (Simons and Carvell, 1989; Brumberg et al., 1996), in order to ensure that the neuron did not reside in the septal regions, and (4) shifts in response latency during periodic stimulation by up to a factor of three, again implying that they likely do not reside in the septa (Ahissar et al., 2000, 2001). For the large set of neurons (n = 44), a full velocity tuning protocol was completed. For a smaller subset (n = 19), in addition to the velocity tuning with punctate stimuli, the influence of angular velocity on pairwise interactions of vibrissa deflections was measured. Finally, for a smaller subset of neurons (n = 15) we were able to measure both the velocity tuning and velocity-dependent pairwise interactions described above, as well as the response to patterns of vibrissa deflections and their corresponding transformations (recording duration ~ 2.5 hours, see below). All data were analyzed using Matlab (Mathworks, Natick, MA, USA) and C programming languages.

Vibrissa Stimulation

A multi-layered piezo-electric bending actuator (range of motion: 1 mm, bandwidth: 200 Hz; Polytec PI, Auburn, MA) generated the vibrissa deflections used to drive SI neurons. The actuator was placed 10 mm from the face, and the vibrissa was inserted into a 4 cm section of a 20 μL glass pipette that was fixed to the end of the actuator (approximately 0.57° per 100 μm deflection). Input to the actuators was generated through a custom-made GUI program written in Visual C++ 6.0.

In designing the command input to the actuators, we accounted for the free-vibration dynamics of the piezo-electric bender to minimize the bender’s extraneous vibrations (i.e., “ringing”). The free-vibration dynamics of the piezo bender were estimated through white-noise analysis (see Marmarelis and Marmarelis (1978)): i.e., by presenting a white-noise input to the piezo, and recording the resulting vibrations with a photo-diode circuit (similar to those described by Andermann et al. (2004); Arabzadeh et al. (2005)). Briefly, as the tip of the actuator interrupts the path of the infrared beam generated by a slotted infrared switch (QVA11134, Fairchild Semiconductor), the output voltage of the optical switch varies in direct proportion to the tip displacement. This change in voltage was subsequently transformed to the corresponding change in actuator displacement using the inverse of the empirically determined relationship between small changes in beam occlusion (measured in terms of the displacement of the actuator tip), and the output voltage. The punctate pulse deflections consisted of exponential rising and falling phases, whose onset slopes (ω) corresponded to the desired velocity transient. These pulses were approximately 10 ms in duration in all cases. Patterns consisted of punctate deflections in the rostral or caudal directions, with variable angular speed.

As shown in Figure 4 and utilized in Figures 5, 6, and 7, stimulus deflection patterns were created from a range of angular velocities observed using high-speed videography of rats whose vibrissae were in contact with textured surfaces (36 grit sandpaper, see below).

High-Speed Videography of Vibrissa Motion

In an isolated set of non-invasive behavioral experiments in a previously published study (Jenks et al., 2010), adult female rats (Long-Evans; 250g; Charles River Laboratories, Inc., Wilmington, MA) were trained to run back and forth on a linear track. Briefly, we trimmed rows A, B, and E of the vibrissa array to facilitate imaging. Tactile stimuli were either oriented spatial grating, with periodic spacing of 4.19 mm, or standard sandpaper (36 grit) surfaces, that were placed on either side of the track. Video was acquired using a high-speed CMOS camera (Phantom v5, Vision Research, Wayne, NJ), at 2000 frames per second and with a spatial resolution of between 700 × 400 and 950 × 512 pixels per frame, with a spatial resolution of between 5 and 8.35 pixels per mm. Tracking of head and vibrissa movements was implemented by modifying a software package provided by Knutsen, Derdikman, and Ahissar (Knutsen et al., 2005). The software uses code written in MatLab (v6.5, MathWorks, Natick, MA) and the C programming language. Vibrissae were tracked in the head frame of reference using a 3-point spline fitting algorithm (see (Knutsen et al., 2005) for details). The angular deflection of a vibrissa close to its base was estimated by calculating the angle of a straight line segment between the extreme points of a short (3mm long) spline fitted to a segment of the vibrissa close to its base (within 2 mm of the mystacial pad).

Response Measures

Peri-stimulus time histograms (PSTHs) of single cell responses to each stimulus type were generated using 60 repetitions of the stimulus, and a time resolution of 2 ms. Spike count was computed as the average number of spikes within the 3-30 ms post-stimulus interval, since RSUs have phasic responses to vibrissa deflections with latencies not exceeding 20-30 ms. In addition to spike count, the timing of single cell responses was studied with two measures: (1) latency, defined as the average time delay between the stimulus and evoked spikes, and (2) jitter (Butts et al., 2007), quantifying the temporal precision of the transient responses to punctate deflections. Latency was computed as the center-of-mass of the single cell PSTH response in the 3-30 ms post-stimulus time window. Jitter was defined as the standard deviation of the time to spike by the neuron over the 3-30 ms post-stimulus time window.

The response to a stimulus pair illustrates the simplest type of response interaction. In this case, the response to the second deflection (termed the Test) is attenuated as a function of the temporal separation between the Test deflection and the preceding Conditioning deflection. Earlier findings have described the suppressive effect of response integration in the SI cortex, where a single stimulus induces attenuation of the response to a second deflection over relatively long timescales (10-200ms) (Simons, 1985; Simons and Carvell, 1989; Brumberg et al., 1996; Kleinfeld and Delaney, 1996; Mirabella et al., 2001). The observed cortical dynamics have been attributed to several mechanisms, including intracortical inhibition, depression of the thalamocortical synapse (Chung et al., 2002), and inhibitory feedback between the thalamic reticular nucleus (nRT) and the ventro posterior medial nucleus (VPM) (Desîlets-Roy et al., 2002). The Conditioning-Test Ratio (CTR) curve (see Figure 2B) was used to quantify the attenuation of the Test spike count relative to its unattenuated (baseline) value. Previous studies (Simons and Carvell, 1989; Fanselow and Nicolelis, 1999; Webber and Stanley, 2004) have shown that the variation of the CTR curve with the inter-deflection interval u is approximately sigmoidal. We therefore used the following parametrization to fit all measured CTR curves: f(u)=12A[1+tanh(ut50τ)], where parameters {A, t50, τ} were estimated through nonlinear least-squares fits.

Velocity Tuning Curves

For each recorded cell, the dependence of spike count on instantaneous stimulus velocity was determined by measuring the mean number of spikes in the 3 to 30 ms window following vibrissa deflections of specified velocities. We averaged the resulting tuning curves for all recorded cells to obtain an averaged spike count tuning curve (see Figure 1D). We developed an analytical form to describe how spike count varies with deflection velocity. Several functional forms of the relationship were explored, including exponential, but a power-law was determined to best capture the behavior. The spike count is zero below the threshold of ±20 degs/sec. At higher velocities, the spike count increases with increasing velocities but the rate of increase gradually decreases:


where H(·) is the Heaviside step function, θ = 20°/sec is the velocity threshold, ωmax =850°/sec, and Nmax is the maximum spike count observed at the largest velocities. Note that the different values of m capture the observed asymmetry in the tuning curve. Analogous tuning curves for the dependence of latency μ (Figure 1E) and temporal jitter σ (Figure 1F) on the angular velocity of single vibrissa deflections were determined, and were fit analytically using nonlinear least-squares:



where A1 = 17 ms, A2 = 7.5 ms, λ1 = 0.004 sec/°, and A3 = 4 ms, A4 = 2 ms, λ2 = 0.004 sec/°.

Encoding Model of SI Neurons

We developed a functional encoding model of SI neurons (see block diagram in Figure 3), and fitted it to individual cells using experimental measurements. This dynamic model predicts single-neuron responses to sequences of deflections with varying instantaneous angular velocities. We should note that since the model is constructed from the separately measured tuning properties, velocity and history dependence are purposefully decoupled in this framework, but are likely intertwined in the context of continuous stimuli. The angular velocity of the vibrissae is obviously a continuously varying quantity, whose trajectory has a complex effect on the cortical response (see Discussion). We thus refer to the velocity dependence described here as instantaneous velocity dependence. The model here consists of three components, each of which is described in detail below: (1) Velocity tuning, (2) Suppression dynamics, capturing the effects of previous stimuli on the response to the current stimulus, and (3) Integrate-and-Fire, which captures the supra-threshold excitatory activity of the recorded cortical neuron.

Velocity Tuning Stage

Based on the functional form of the instantaneous velocity tuning for the overall magnitude of the response (Eq. 1), we defined a scaling of the response due to changes in instantaneous angular velocity:


where m = 0.4, now assuming symmetry with respect to positive and negative velocities, for simplicity. Note that h(tk) varies between 0 and 1. We then characterized the velocity stimulus sequence as an impulse train: S(t) = Σiωiδ(tti), where {ti} and {ωi} are the occurrence times and the angular velocities of the sequence of deflections. Since S(t) is a sequence of delta functions weighted by the angular velocities which range from −ωmax to ωmax, the resulting signal h(t) = γ(S(t)) is a sequence of delta functions with weightings scaled from 0 to 1, or h(t) = Σk γ (ωk)δ(ttk). In the absence of history dependence, the response of the neuron is determined by these scaling terms (see the front end of the neural encoding model in Figure 3).

Suppression Dynamics Stage

An SI neuron’s response to a given vibrissa deflection within a pattern is determined by a scaling relative to the response to a single vibrissa deflection, where the scaling captures the interactions between the suppressive effects of all preceeding stimuli. This scaling is denoted by a state variable x(t) that ranges from 0 (in the case of complete suppression of the suprathreshold response) to 1 (in the absence of any suppression), and is Markov in the sense that regardless of past history, the present “state” is the only quantity necessary to predict how the neuron will respond to the next deflection. Note that in some cases the state can also exceed 1, representing a situation where the response is facilitated above the isolated vibrissa deflection response. As we previously showed, the profoundly nonlinear suppression following the response to a punctate stimulus does not depend upon the firing of the neuron in question, but instead reflects network dynamics and the evolution of the excitatory and inhibitory balance (Webber and Stanley, 2004) - see Discussion.

As an example, suppose the vibrissa is stimulated with a pair of velocity pulses, S(t)=i=12ωiδ(tti). The state is assumed to be 1 prior to the first deflection, reflecting no suppression from prior stimuli (see the trajectory of the state prior to and following the first deflection, illustrated in the suppression stage of the block diagram in Figure 3). Following the first deflection, the state immediately drops to zero, and gradually recovers over a timescale described by the paired-pulse measurements. The state at the time of the second deflection, given the history of stimulation up to that time, is f(t2t1; ω1), where ω1 denotes the dependence of the CTR curve on the angular velocity of the conditioning stimulus. In the model, the velocity dependence for the deflection at time tk is captured by the scalar quantity h(tk), which varies from 0 to 1. We therefore modify the notation for the CTR function to incorporate this dependence, f(tk+1tk; h(tk)). A larger velocity ωk, and therefore a larger h(tk), results in a more profound suppression, and thus a lower value of f for the same time interval. This relationship is experimentally well-described as:


and f(tk+1tk; 1) represents the CTR curve obtained when h(tk) = 1 (i.e. ωkωmax, with ωmax being the velocity at which the velocity tuning curve saturates; see Figure 1). Note that the expression g[a, b] for a = h(tk) and b = f (tk+1tk; 1) can be intuitively interpreted. For a long time interval between deflections at tk and tk+1, the CTR curve f (tk+1tk; 1) goes to 1, as does b. In this case, g[a, b] goes to 1, yielding the maximal response, regardless of the strength of the scaled input h(tk). For a short time interval between deflections at time tk and tk+1, the CTR curve f (tk+1tk; 1) goes to 0, as does b. In this case, g[a, b] goes to 0, yielding no response, regardless of the strength of the scaled input h(tk). For intermediate values of the inter-deflection interval tk+1tk, the velocity of the deflection at time tk, captured by h(tk), will influence the amount of suppression induced. For example, for an inter-deflection interval such that the response is half of its maximum, and thus the CTR is 0.5, a strong preceding input of h(tk) = 1 will result in g[a, b] = 0.5. However, a weaker preceding input h(tk) = 0.5 will result in g[a, b] = 0.67, shifting the effective CTR to the left, or lessening the suppression on the response to the next deflection at tk+1, consistent with the experimental observation.

Generalizing to a sequence of pulses with varying velocities, the evoked neural responses can be predicted using only the knowledge of pair-wise interactions between stimuli – i.e., using only measurements of CTR curves. The following recursive relationship provides an accurate prediction of the state variable at the occurrence times of the stimulus {tn}:


which captures the nonlinear evolution of the suppression. Note that as we have previously described, the suppression does not simply “build-up”, but instead exhibits complex interactions (Webber and Stanley, 2004). Our previous studies have shown that in the context of sequences of transient, high-velocity vibrissa motion, the cortical response is strongly nonlinear, and is thus not well-described with the cascade of a linear filter and a simple, rectifying static nonlinearity (Webber and Stanley, 2004; Boloori and Stanley, 2006) - see Results and Discussion. The suppression of a response to a deflection at tk due to a prior deflection at tk−1 results in a weaker suppression of the response to a deflection at tk+1, which is captured by the above recursion.

Integrate-and-Fire Stage

The firing response of the cortical neuron is then determined by the “drive” d(t) = x(t) · h(t), as shown in Figure 3A. We employ a standard integrate-and-fire model of neuronal firing (Dayan and Abbott, 2001) to transform the sensory driven input to neuronal spikes:


where τ is the membrane time constant (set to 10 ms), d(t) is the sensory driven input scaled by a constant α, and n(t) is a noise term included to capture the timing variability in the model neuron response. When the potential reaches a threshold Vthresh, the model fires an action potential and resets (for implementation, Vrest was set to -70mV, and Vthresh was set to -40mV). The synaptic input to the integrator is the drive term αd(t), which is a pulse whose magnitude depends upon the state x(t), as described above. Depending on the past history of sensory input, a transient stimulus event produces a pulse input to the integrator, to which noise n(t) is added, resulting in a threshold crossing that is variable in its timing. The latency between the pulse input and the threshold crossing is determined by the magnitude of the drive input, resulting in varying the slope of the rising voltage. The noise was a Gaussian white noise process whose variance was set to yield the appropriate levels of timing jitter (the ratio of α to the standard deviation of n(t) was 2.6, and the signal-to-noise ratio (SNR) varied as a function of the magnitude of the drive d(t)). Importantly, stronger inputs produced shorter latencies, and decreased timing jitter, as seen in the experimental data, even though the variance of the noise was fixed, due to the changing relationship between the magnitude of the input and the fixed magnitude noise (Mainen and Sejnowski, 1995; Butts et al., 2007). As shown in Figure 3B-D, the integrate-and-fire mechanism produced spike counts, latencies, and jitters that were consistent with the experimental observations presented in Figure 1.

Decoding Neural Responses

The decoding of sensory inputs from the observed neural responses involves the detection of significant neural response events and Bayesian inference of the stimulus pattern from an assumed internal nonlinear model of neural encoding.

Detection of Response Events

The first step in decoding is the classification of single cell firing activity into significant response events or spurious noise. For the purpose of detecting the response events evoked by underlying deflections, neural firing rates were re-binned at 30 ms. From this point the measured neural firing rate (i.e. the PSTH) is referred to as the “firing rate”, whereas the 30 ms binned response is termed the “binned-rate”. An ideal observer of spike counts seeks to detect the presence of stimulus deflections on the basis of observing the evoked binned rate. Two types of error are associated with this discrimination task: “Misses” characterize occasions where the observer does not label a response event as significant, though an actual deflection had evoked its occurrence. A “False Alarm”, on the other hand, denotes mistaken identification of spurious neural activity as that due to a stimulus.

Following classical detection theory (Green and Swets, 1966), optimal detection of a “signal” (i.e., presence of a stimulus deflection S within a time-bin) that is corrupted by noise (i.e., due to spontaneous neural firing) uses two conditional probability distributions: (1) probability of the observed event size N, given the presence of a stimulus Pr[N[mid ]S], and (2) the PDF of the observed event size N given the absence of a stimulus Pr[N[mid ] ~ S]. For the observed spike count N′ in each time-bin, we constructed the likelihood-ratio L(N′) as:


The Neyman-Pearson criterion (Green and Swets, 1966) allows an observer to optimally detect the presence of stimuli subject to a maximum tolerated probability of False Alarm. In practice, we carried out the optimal detection through the likelihood-ratio test (Green and Swets, 1966): if L(N′) ≥ β, where β is a “threshold” value, then the observed response event of size N′ was attributed to a stimulus. For very small threshold values, almost all computed likelihood-ratio values were larger than the threshold, leading to identification of most response events as evoked by stimuli: both the False Alarm and correct-discrimination probabilities were therefore very high (close to 1). However, both the False Alarm and correct-discrimination probabilities decreased with increasing β. The optimal value of β depends on the specified probability of False Alarm: β was increased from zero until the likelihood-ratio test provided the maximum tolerated rate of False Alarm observations of 0.10.


The outcome of the “detection” step was a sequence of observed spike counts {Nk}, and their occurrence times {tk}. The decoding process consisted of the following steps that were carried out in sequence for each observed spike count Nk, starting with the first observation:

  1. For each observed significant spike count Nk, we estimated the drive d(tk) at the time of observation as d(tk) = Nk/Nmax, for all tk.
  2. Next, x(tk) and h(tk) were estimated recursively. For k = 1, x̂(t1) = 1, giving ĥ(t1) = d(t1). For k = 2, x̂(t2) = f (t2t1; ĥ(t1)), giving ĥ (t2) = d(t2)/x̂(t2), and so on, for all tk. Conceptually, this is based on an assumed internal model of the dynamics related to history dependence.
  3. Finally, the absolute value of the angular velocity was estimated as [mid ][omega with circumflex]k[mid ] = γ−1(h(tk)) for all tk, using equation 4. Note that the sign of the angular velocity, and therefore the direction of deflection, was ambiguous through the observation of a single neuron.


During behavior, rats and other rodents use their facial vibrissae to actively explore surfaces through whisking and head/body movement, resulting in complex sensory inputs that vary over a large range of angular velocities and temporal scales. Through single-unit cortical recordings of monosynaptic thalamorecipient neurons in S1 of the anesthetized rat, we first systematically quantified the interactions between velocity and timing of vibrissa motion to form a predictive encoding model which captures firing probability and timing characteristics of the cortical representation.

Velocity Tuning in SI Cortex

During behavior in which the animal is actively exploring a surface (Figure 1A), the vibrissae undergo patterns of discrete high-velocity transients (Ritt et al., 2008; Wolfe et al., 2008), where the velocity and timing of the transients depend on the textural properties, the intrinsic mechanics of the vibrissa complex, and the precise nature of the behavior. Figure 1B shows an example of the high-velocity transient events experienced by a single vibrissa in a behavioral setting, where the vibrissa motion was tracked using high-speed videography while the animal traversed a linear track lined with periodic textures (Jenks et al., 2010). In this ethologically relevant context, however, high-velocity transients (dotted boxes) are embedded within a highly complex sensory pattern which precludes a simple analysis. To investigate these high-velocity transients in controlled conditions, we first recorded extracellularly from S1 neurons while stimulating the primary vibrissa with isolated, short-duration pulse deflections of varying onset slopes (see the schematic at the top of Figure 1C). These simple stimuli allowed us to systematically parameterize the instantaneous velocity tuning effects on response timing and magnitude of topographically aligned neurons in cortical layer 4. Due to the short duration of the pulse (10 ms), the observed neural response was solely due to the onset slope, since the post-excitatory suppression induced by the rising phase completely suppressed the response to the falling phase (Webber and Stanley, 2004). Varying the onset slope allowed for the measurement of both the magnitude and temporal profile of single cell response as a function of velocity. Each trace in Figure 1C shows the average firing rate response (peri-stimulus time histogram, or PSTH) to a different pulse velocity presented at time zero (n = 44 cells; dashed black: 30°/sec, gray: 120°/sec, solid black: 920°/sec). As stimulus velocity of the isolated punctate stimulus increased, cortical neurons responded with higher magnitude, lower latency, and increasing temporal precision (decreasing jitter).

Figure 1
Responses of SI neurons systematically vary with vibrissa velocity

The instantaneous velocity tuning curves in Figure 1D-F summarize the measurements and observations over a wide range of angular velocities. Figure 1D shows the average spike count per stimulus (computed by integrating the PSTH response over the 3-30 ms post-stimulus time window) as a function of velocity, with nearly all neurons displaying saturation for velocities greater than 500°/sec. The solid curve is a parametric power-law fit of the relationship between the angular velocity and the spike count, as described in the Methods. The average tuning curves for all recorded cells are approximately symmetric over rostral (positive) and caudal (negative) velocities. Additionally, almost all SI neurons studied possessed a velocity threshold below which they did not respond to pulse deflections – see inset high-lighting the low-velocity regime of the tuning curve. Given the average measured spontaneous firing rate of approximately 3 Hz, we estimate that the velocity threshold is approximately ±20°/sec.

As expected from the variation of PSTH responses with instantaneous velocity from Figure 1C, the response latency smoothly decreased with increasing instantaneous stimulus speeds, as shown Figure 1E, again with the solid curves representing a parametric exponential fit of the relationship. Response latency had a lower bound of approximately 10 ms at the highest measured velocities. The temporal precision of the response increased with instantaneous velocity, as evidenced by the decrease in temporal jitter in the bottom plot (see Materials & Methods). Importantly, the changes in spike count with angular velocity were accompanied by corresponding changes in latency and timing jitter, and the latency and jitter were well predicted by the spike count, as shown in Figure 1G-H. In each, the timing (latency or jitter) is plotted as a function of the spike count, and the corresponding parametric fits in D-E are combined to form explicit parametric relationships (solid curves).

History Dependence of SI Responses

In addition to the variations in the angular velocity of the vibrissa motion, the relative timing of the deflection transients is an important feature of the sensory input which is strongly influenced by behavior. The simplest way to systematically characterize the dependence of responses on the past history of stimulation – caused by post-stimulus suppression – is through the Conditioning-Test paradigm (Kyriazi et al., 1994; Webber and Stanley, 2004; Boloori and Stanley, 2006). Note that we have previously shown that this form of history dependence is stimulus driven, rather than an intrinsic, biophysical property of the neurons triggered by action potentials (such as prolonged refractoriness or after-hyperpolarization) (Webber and Stanley, 2004). Drawing inspiration from patterns of high-velocity transients in behavioral contexts (Figure 1B), a simple paired stimulus was used, consisting of a Conditioning deflection (C) followed by a Test deflection (T), as shown in Figure 2C (also conceptually highlighted by the high-velocity transients in Figure 1B).

Although the Conditioning-Test paradigm has been well studied, the relationship between velocity tuning and inter-deflection timing is not known. Figures 2A and B show the average PSTHs for a range of paired stimulus conditions. In Figure 2A, the PSTHs are shown for paired deflections in which the angular velocity of the second (Test) deflection was fixed at 690°/s, but the angular velocity of the first (Conditioning) deflection was varied, for a range of inter-deflection intervals (IDIs). The responses qualitatively show that the response to the second deflection is increasingly suppressed for shorter IDIs, and that a higher velocity Conditioning stimulus produces a more profound suppression of the response to the Test stimulus for a given IDI. Figure 2B shows the same for a fixed Conditioning stimulus at 690°/s, while varying the angular velocity of the Test stimulus over a range of IDIs. Qualitatively, for a given IDI, we observe that increasing angular velocity of the Test stimulus serves to allow the neuron to emerge from the suppression induced by the first deflection more quickly. It is important to note that the properties of the paired-deflection response shown in Figure 2 could arise from a very simple linear model/filter that involves an initial excitatory response, followed by a prolonged inhibitory tail. In this case, the response to a sequence of two deflections such as that in the figure would result in the response to the second deflection “riding” on the inhibitory tail, which when rectified, would result in an attenuated second response as a function of the inter-deflection interval, as we observe. However, we have previously shown through a simple three deflection sequence that the mechanism is not additive, but instead divisive in nature, and strongly nonlinear (Webber and Stanley, 2004, 2006; Boloori and Stanley, 2006), also demonstrated in intracellular cortical recordings (Higley and Contreras, 2003), and thus more powerful models are required for complex stimulus patterns.

Figure 2
Responses of SI neurons systematically vary with stimulus history

To better quantify these observations, the Conditioning-Test Ratio (CTR) was calculated for a given inter-deflection interval as the ratio of the spike count in response to the Test deflection preceded by the Conditioning deflection to the spike count in response to the Test deflection alone, where the spike count is taken in the 3 to 30 ms window following the deflection (Figure 2C). The CTR describes, as a function of the inter-deflection interval u, the degree to which the response to the test stimulus has been suppressed/excited by the conditioning stimulus, which is strongly nonlinear (Webber and Stanley, 2004; Boloori and Stanley, 2006). The timescale of post-stimulus suppression depends significantly on the velocity of the Conditioning stimulus. For a fixed Test stimulus, progressively higher Conditioning velocities evoked longer-lasting post-stimulus suppression in a majority of cells, as shown by the slower recovery time of the corresponding CTR curves (see Figure 2D). For a specific inter-deflection time-interval, therefore, the degree of suppression induced by a transient deflection is strongly dependent upon its angular velocity. Similarly, for a given velocity of Conditioning stimulus, increasing the velocity of the Test stimulus results in effectively shortening the time-course of the suppression induced by the first deflection, as shown in Figure 2E. Again, changes in response magnitude were well correlated with other properties of the response suppression; a decrease in response magnitude was accompanied by a decrease in precision and an increased latency.

From the perspective of an observer of single unit cortical activity, the above results point to ambiguities in the neural response when decoding the sensory input. On the one hand, the same isolated punctate stimulus deflection can result in dramatically different properties of the corresponding localized firing activity (spike count, precision), depending on past history, as illustrated by the CTR curves in Figure 2D-E. On the other hand, observation of a small response event (i.e. low spike count) could arise from a low-velocity vibrissa deflection, as described in Figure 1, or from a high-velocity stimulus whose response has been suppressed by a preceding deflection, as described in Figures 2D, or a range of dynamic trajectories which ultimately lead to that specific level of suppression/excitation when paired with a transient deflection. Thus, a transient deflection at a specific velocity does not result in a unique neural response, but is instead strongly dependent upon the past history. From both perspectives, for proper decoding (and perceptual interpretation), an encoding model which takes into account these properties is necessary.

An Encoding Model of SI Responses

To capture the encoding dynamics of the vibrissa pathway, we constructed a functional response model incorporating the dependence of SI response magnitude and timing on instantaneous deflection velocity and the nonlinear effects of stimulus history. Figure 3 schematizes the model, which transforms a sequence of vibrissa deflections into the firing activity of a single cell. The input to the vibrissa pathway can be envisioned as sequences of angular deflections, modeled here as punctate deflections with varying onset velocities (see Materials & Methods for a detailed description of the model). We should note that the biophysical processes that give rise to the observed cortical responses are complex, and that this model is not designed to capture this explicitly, but instead integrates the experimental observations described thus far within a functional framework that is analogous to functional models in other pathways (Paninski et al., 2004; Lesica et al., 2007; Mante et al., 2008; Pillow et al., 2008). Note, however, that at its core, the model presented here is highly nonlinear in form, differentiating it from sensory models based on linear filtering (see Discussion).

Figure 3
A dynamic encoding model for prediction of SI cortical responses to deflection patterns of varying velocity and timing

The velocity tuning stage of the model captures the scaling of the response as a function of the instantaneous angular velocity of the vibrissa deflection, quantified using a parametric fit to the experimentally-measured tuning curve (see Materials & Methods). The history dependence of the pathway captures the nonlinear dynamics of the suppression that follows a vibrissa deflection, modeled through an evolution of a “state” variable x(t) based on temporal interactions of the experimentally-measured CTR curves which were fit parametrically. As described in the Methods, this is distinctly different from suppression of activity following spiking, due to refractory or other spike-history dependent effects. This term instead reflects the relative level of the momentary balance between excitation and inhibition in the thalamocortical circuit, a highly nonlinear process. The resulting combination of velocity tuning and history dependence, which we denote the “drive” d(t), determines the firing activity through an integrate-and-fire stage of the model. In response to a transient deflection of the vibrissa, the input to the integrator is a pulse whose magnitude is determined by the drive d(t). The input, analogous to a synaptic current, is integrated, resulting in a depolarization that drives the mechanism to threshold, firing a spike. Noise is included in the input to the integrator, resulting in variability in timing of the integrate-and-fire model’s spiking.

Figure 3B shows the firing event parameters (spike count, latency, and jitter) in response to repeated presentations of a single transient drive input to the integrate-and-fire mechanism for varying levels of drive. The response characteristics of the integrate-and-fire model are qualitatively similar to those observed in the velocity tuning of actual neurons (see, e.g., Figure 1D-F). Increasing the magnitude of the input increases the spike count, decreases latency, and decreases timing jitter. Figure 3C shows the PSTH output of the integrate-and-fire mechanism for varying levels of drive, highlighted by the vertical lines in B. Note that latencies produced by the integrate-and-fire mechanism were superimposed on a latency of 10 ms to account for the synaptic delays from the periphery to cortical layer 4, since the integrate-and-fire model is conceptually meant to represent integration at the thalamocortical synapse. Finally, Figure 3D shows the relationships between the spike count and latency and spike count and timing jitter for the integrate-and-fire model (solid curves), along with the experimentally observed relationships from Figure 1G and H (symbols with error bars), illustrating that predicting the overall drive in the circuit captures not only the magnitude of the response, as measured through spike count, but timing aspects of the response as well.

Model Prediction of Cortical Response to Stimulus Patterns

To evaluate the predictive capability of the model, we designed simplified artificial stimulus sequences based on the range of behaviorally observed distributions of angular velocities. Figure 4A shows the instantaneous velocity of a single vibrissa in contact with a coarse texture (36 grit sandpaper), as measured through high speed videography (Jenks et al., 2010) – see Methods. The full distribution of instantaneous angular velocity is shown by the black curve in Figure 4B. The horizontal dashed gray lines in Figure 4A represent a threshold imposed at ±160°/s, to identify transient stimulus events of relatively high angular velocity. The inter-deflection intervals were restricted to be greater than 40ms, resulting in the distribution of inter-deflection events shown in Figure 4C. The final set of stimulus events had the angular velocity distribution shown by the gray curve in Figure 4B. Note that these stimulus sequences were designed to provide a simple, controlled input to explore the dynamics of the pathway, while utilizing angular velocities and timings of transient vibrissa motions that are observed in more naturalistic settings.

Figure 4
Design of stimulus deflection sequence from behavioral observations

The top panel of Figure 5A shows a short segment of a rostral-caudal vibrissa deflection sequence generated by random draws from the distributions in Figure 4B and C. Each arrow indicates the time at which a punctate deflection was delivered, where the arrow height and direction indicate the velocity of the deflection (see Methods). The panels below the stimulus sequence show the experimentally measured PSTH response (PSTH averaged across a small subset of cells with similar dynamics, and therefore similar CTRs; n = 6) and the corresponding model prediction (inverted scale) for the Full Model, the model in the absence of the velocity dependence (No-Velocity, middle), and the model in the absence of the history dependence (No-History, bottom). The model was based on the average tuning curve and CTR measurements for the cells, and was utilized to predict cortical responses to stimulus sequences not used to fit the model. Localized Gaussian templates of the PSTH events were utilized for the plotting of the predicted response events, with the appropriate spike count, latency, and jitter as predicted by the encoding model. The Full Model prediction exhibits fairly good agreement with the experimentally measured data, predicting the presence of the firing events and the corresponding timing and magnitude, while also correctly predicting the absence of firing in response to some deflections (e.g. 5th stimulus event). In contrast, removing the velocity dependence results in a gross over-estimation of the firing activity, and an under-prediction of the latency and jitter (middle panel, No-Velocity). Most dramatically, however, when the strong nonlinearity of the paired-pulse suppression was excluded from the model, the model predicted significantly greater firing than measured experimentally, particularly predicting firing events in response to deflections during short IDIs that resulted in little or no actual firing (bottom panel, No-History).

Figure 5
Model prediction of cortical activity

For a quantitative assessment of the predictive capability of the model, event-by-event spike counts were evaluated for the measured responses and the model predictions, an example of which is shown for the Full Model in Figure 5B (measured spike count in black, model prediction in gray). The standard errors in the measurements of spike counts are not shown for clarity, but at almost all points they overlap with the values predicted by the model.

To quantify the goodness-of-fit over the longer stimulus trial, the correlation coefficient and mean-squared error (MSE) were calculated between the predicted and actual spike count for a larger sample of cells (n = 15). Note that the MSE reflects absolute deviations from the experimental observations and thus penalizes for improper scaling, serving as a complement to the correlation coefficient. We carried out two null-tests to determine the relative importance of the different constituent elements that comprise the model. In the first test, we predicted the cortical response by retaining history dependence but ignoring velocity tuning (by setting the tuning curve to a flat curve equal to the maximum spike count, as previously shown in Figure 5A middle). As shown in Figure 5C, top panel, the resulting correlations between the predicted and the measured values (denoted as No-Velocity) were the same as those observed for the full model (bars indicate 99% confidence intervals). This was due to the fact that accounting for history dependence allowed the model to track the time-course of suppression relatively well, even though the magnitudes of the predicted values were in error due to the neglect of velocity tuning. However, the resulting MSE for this test was significantly larger than that observed for the full model, showing that accounting for velocity tuning was important for the correct prediction of the magnitude of the responses (Figure 5C bottom, Full Model versus No-Velocity). In the second test, we predicted the cortical response by accounting only for velocity tuning, while ignoring history dependence, an example of which was previously shown in Figure 5A, bottom. The resulting correlation coefficient and MSE between these predictions and actual measurements are shown in Figure 5C (No-History). In comparison with the full model, ignoring history dependence drastically lowered the agreement of the predictions with the actual measurements in terms of both the linear correlation and the MSE. Similar results were obtained for other stimulus sequences generated (not shown). The above tests demonstrate the importance of both velocity tuning and history dependence to an accurate model. Comparison of the resulting correlation coefficients and prediction errors reveals that accounting for history dependence is necessary to appropriately predict the temporal pattern of the neural response, but that accounting for velocity tuning is critical in accurately predicting the magnitude of the firing events within the temporal pattern of the neural response. It should be noted that the predicted response from the full model achieved a correlation of 0.64 with the actual recorded response, which when compared to the correlations across subsets of the experimentally recorded responses (0.85) implies that the model captures approximately 75% of the explainable variance (see Discussion).

The predictive capability of the model here is consistent with findings for single unit model predictions at a number of sensory stages. Models of retinal ganglion cells predict 80-95% of the explainable variance (Pillow et al., 2008), under very well controlled, in-vitro conditions. LGN models typically predict less than 56% (and as low as 24%) of the explainable variance (Mante et al., 2008). In primary visual cortex, the best models to date predict approximately 40% of the explainable variance (personal communication, J. L. Gallant). The mismatch between the model prediction and the observed responses here could result from several factors. The encoding model imposed here is fairly parsimonious in nature, using only the power law fit to the velocity tuning (4 parameters), the timecourse of the sigmoidal fit to the CTR (3 parameters), and the integrate-and-fire model at the back end (3 parameters, excluding Vrest and Vthreshold). Further, this is fit to population-averaged data, thus the cell-to-cell variability in tuning properties could degrade the overall performance. Finally, the model does not take into account spike-history dependence, which may play a role in shaping the cortical response.

Bayesian Decoding of Tactile Patterns

Given the fidelity of the above encoding model, we can then ask the question, how much does the observation of the cortical response reduce the uncertainty about the sensory stimulus? The cortical response, as we have described, is temporally localized and event-like, linked to discrete features of the vibrissa input. As described in the Methods, the first step in this process is thus to identify response events, as shown in Figure 6A.

Figure 6
Bayesian decoding of tactile patterns

Shown are the average cortical PSTH (gray) and corresponding spike-count measure in 30ms windows (black) in response to a stimulus sequence of randomized punctate deflections of varying angular velocity and timing (shown in Figure 6C). As described in the Methods, the Neyman-Pearson criterion was utilized to detect response events from the localized spike-count measured, indicated with the gray bands in Figure 6B. The stimulus sequence could then be parsed into stimulus events that correspond to detected response events and those that do not (gray and black, respectively, in Figure 6C). The stimulus events that correspond to detected response events are thus those that are theoretically predictable from the observed cortical activity, while the remaining stimulus events cannot, and thus are not considered in the decoding of Figures 6D and E. Using the Bayesian decoding framework described in the Methods, the observed sequence of response events, the population encoding model, and the prior distribution of the stimulus (Figure 4B and D) were then utilized to infer the stimulus sequence that most likely induced the response. As previously noted, the angular velocities were both positive and negative (rostral and caudal), which creates an ambiguity in the decoding from a single neuron. For display purposes, here we assign the decoded stimulus velocity the correct sign (see Discussion), although the decoder yields a prediction that is sign-ambiguous. Figures 6D and E show the actual stimulus events (gray), and the decoded stimulus events (black), for the Full Model (E) and the No-History (D) cases. The No-History case, as before, is where the decoding is based on velocity tuning only, neglecting the nonlinear suppressive dynamics. As shown in comparing the two cases, neglecting the history dependence of the suppression causes the decoder to consistently underestimate the angular velocity of the stimulus events. This is summarized in a scatter plot of the actual stimulus angular velocity versus the corresponding prediction from the Bayesian decoder, as shown in Figure 6F. Decoding in the No-History case (gray) results in points that are significantly below the unity line for positive angular velocities, and points that are significantly above the unity line for negative angular velocities. The predictions based on the Full Model (black), however, are found to lie consistently closer to the unity line.

Encoding and Decoding Under Simple Transformations

The tuning properties shown in Figures 1 and and22 reveal strong dependencies of the cortical response on deflection velocity and timing, and further testing of the model revealed that accounting for both is necessary for accurate prediction of the cortical response. We have recently quantified the degree of variability in vibrissa deflection patterns when the vibrissa array is engaged with various surfaces during free locomotion (Jenks et al., 2010). Head and body movements create a significant amount of variability in the vibrissa deflection patterns for the same surfaces, primarily manifested in the timing and amplitude/velocity of vibrissa contacts. Given the importance of both of these factors, and the associated variability that we observed behaviorally, we designed a new set of stimuli by scaling along the two axes of velocity and timing, to test the encoding and decoding under these transformations.

Figure 7A shows a short segment of a rostral-caudal vibrissa deflection sequence (top) and the corresponding peri-stimulus time histogram (PSTH) of a typical cortical layer 4 single-unit response. As above, each arrow indicates the time at which a punctate deflection was delivered, where arrow length and direction indicate the velocity of the deflection (see Methods). Note that the 650 ms segment was taken from the middle of a longer 6 second segment. Figure 7B illustrates the same deflection pattern, preserving the temporal structure in A, but with all instantaneous velocities scaled up. Since this scaling was analogous to increasing the temporal “contrast”, we defined the deflection sequence in A as the Low Contrast (LC) sequence and that in B as the High Contrast (HC) sequence.

Figure 7
Cortical responses to simple transformations of tactile stimulus patterns

To independently explore the role of timing, we time-compressed the original LC sequence, while retaining the same instantaneous velocities of each LC stimulus event. The resulting Time Compressed (TC) sequence and the corresponding cortical response are shown in Figure 7C (TC, gray), along with the corresponding LC response for comparison (black, inverted scale). Note that the PSTH for the TC stimulus (gray) is shown at an expanded timescale to align both PSTHs to the same stimulus pattern.

The previously described encoding model was then used to predict the response stimulus patterns under this paradigm. Figure 8A shows the actual response (gray, positive axis) versus the model prediction (black, inverted axis) for the LC (left), HC (middle), and TC (right) cases (averaged actual response, n = 6 cortical neurons). The top row shows the performance of the Full Model, while the bottom row shows the same for the No-History case, where, as before, the predicted response is based solely on the velocity tuning properties.

Figure 8
Encoding and decoding under stimulus transformations

Removal of the nonlinear, suppressive history term resulted in a significant over-estimation of the firing activity. The scatter plots in 8B show the actual spike count versus the predicted spike count, for the Full Model (solid symbols) and the No-History case (open symbols), for LC, HC, and TC stimuli. The Full Model performed similarly well for the LC, HC, and TC cases, and the history term is critical for the model prediction (see figure caption for corresponding correlation coefficients). A single model thus captures the nature of the cortical response under different stimulus contexts.

Given knowledge of the encoding dynamics, we can then ask to what extent an ideal observer of the cortical activity can decode stimulus-related information from neural responses. As in Figure 6, the response events were first detected from the recorded responses, and then the encoding model was used in a Bayesian framework to infer the most likely stimulus sequence. Figure 8C shows the actual versus decoded event-by-event angular velocities. Note that again, the stimulus events shown are only those that correspond to a detectable response event, and the correct-sign of the angular velocity was assigned to the sign-ambiguous decoded angular velocity. Shown are the decoding results for LC, HC, and TC (circles, squares, triangles, respectively), for the Full Model (solid, black), and the No-History case (open, gray). Again the nonlinear, suppressive history term was critical in the decoding, with significant underestimates of the angular velocity when neglected. The decoding performed similarly well in each of the three stimulus classes (see figure caption for performance metrics). Theoretically, the three stimuli described here may give rise to different cortical activity, yet arise from transformations of the same surface/texture, as conceptually outlined in Figure 8D. The above results suggest that the same encoding model holds up under these hypothetical transformations, demonstrated through both predictive encoding and decoding in the stimulus space. However, to what extent the animal has access to internals models of self-motion, and how these might be combined with internal models of neural encoding, remains an open question.


It is the underlying tenet of this work that features of the external environment are represented in the spiking activity of neurons in early sensory pathways over a relatively fine temporal scale (on the order of tens of milliseconds), rather than firing activity averaged over large temporal windows. Neural activity at fine timing precision has been demonstrated in a wide range of sensory areas, including the retina (Berry and Meister, 1998; Liu et al., 2001; Passaglia and Troy, 2004), visual thalamus (Reinagel and Reid, 2000; Butts et al., 2007; Desbordes et al., 2008), auditory periphery (Cariani and Delgutte, 1996), and visual (Bair and Koch, 1996; Buracas et al., 1998) and auditory cortices (Wehr and Zador, 2003; DeWeese and Zador, 2006). Under acute conditions in the vibrissa pathway, in which precisely the same motion of the vibrissae can be repeated, cortical neurons exhibit temporally precise, stimulus-locked activity (Garabedian et al., 2003; Khatri et al., 2004; Boloori and Stanley, 2006) that suggests the possibility of such a neural code based on temporally precise neural activity reflective of textural properties, which has been described as a “kinetic signature” that is consistent with our findings here (Arabzadeh et al., 2005, 2006). However, recent chronic recordings from the awake animal engaged in a texture discrimination task have revealed patterns of cortical activity that are highly variable from trial-to-trial, which has led to the suggestion of a rate-based code in which the relevant information is encoded in the mean firing rate over hundreds of milliseconds (von Heimendahl et al., 2007). One central difference between these two experimental paradigms is that the acute studies permit the presentation of a precisely controlled tactile input that is repeatable across trials, whereas the behavioral studies are subject to trial-to-trial variations in sensory input even in the presence of the same tactile features of the external environment. It is important to note that although significant variability can exist from trial-to-trial for individual neurons, local populations on a single trial are all experiencing the same sensory input that has been shaped by the particular behavioral parameters during that trial; the degree of precision across populations remains an open question in this pathway. The lack of a temporally precise response across behavioral trials is not evidence against the importance of the precise features in the neural response, but instead makes more salient the need for internal models of self-motion for interpretation of the incoming, distributed sensory signal.

Velocity Sensitivity and History Dependence

The instantaneous velocity sensitivity of cortical neurons in the vibrissa system has been measured in a number of previous studies (Simons, 1978; Ito, 1985; Pinto et al., 2000; Shoykhet et al., 2000). The majority of work in this regard has focused on the spike count of the response to punctate or ramp-and-hold deflections, reporting a monotonic dependence of spike count on deflection velocity consistent with our present results. It has been shown that the timing of thalamic activity is important in the context of integration at the thalamocortical synapse (Usrey et al., 2000; Swadlow and Gusev, 2001), and asserted that the velocity dependence of the cortical spike count is established by thalamic timing (Pinto et al., 2000). Indeed, recent findings based on paired recordings of thalamic units and their cortical targets suggest that layer 4 cortical neurons in the vibrissa system are driven by weak input from a relatively large number of thalamic projections (Bruno and Sakmann, 2006). It is thus likely that the velocity sensitivity of the cortical spike count arises from the enhanced time-locking of thalamic neurons with increasingly strong stimuli. In addition, we showed that the velocity has a systematic effect on the latency and temporal precision of the cortical response, which together with the spike count shapes the cortical representation. Although cortical velocity sensitivity has been clearly established, the distribution of velocities experienced by the animal in more ethological settings has only recently been addressed in experimental and computational studies (Arabzadeh et al., 2005; Ritt et al., 2008; Wolfe et al., 2008). The simple patterns of velocity transients used here were designed to decouple the velocity and timing elements of the cortical representation of the vibrissa motion, rather than to explicitly represent naturalistic input to the pathway, which could potentially involve more continuous motion. However, there is recent evidence that vibrissae resonate transiently during discrete, high-velocity, and high-acceleration slip-stick events during whisking against surfaces of various types, and that the rate/pattern of these high-velocity vibrissa transients varies with the properties of the texture and self-motion (Ritt et al., 2008; Wolfe et al., 2008).

The history-dependent effects characterized by the conditioning test ratio measurements we obtained through paired-pulse stimuli are similar to those previously reported (Simons, 1985; Simons and Carvell, 1989; Brumberg et al., 1996; Kleinfeld and Delaney, 1996; Mirabella et al., 2001; Webber and Stanley, 2004, 2006; Boloori and Stanley, 2006). Barbiturate anesthesia enhances and mimics the action of GABA at the GABAA receptor complex, decreasing excitability at both cortical and subcortical stages, and thus likely affects the measurements presented here. There is a reported reduction in the amount of cortical adaptation with increased arousal, attributed to the overall depression of the thalamocortical synapse (Castro-Alamancos, 2004), although a more recent study showed that the thalamo-cortical synapse is not altered by the state of arousal, but instead the observed attenuation in response may result from a modulation in thalamic synchrony (Stoelzel et al., 2009). Paired-pulse electrical stimulation of the infraorbital nerve using a cuff electrode showed that the stimulus-induced suppression does persist in the aroused state, but over a shorter timecourse as compared to the awake but quiescent or anesthetized states (Fanselow and Nicolelis, 1999). The conceptual nature of the model proposed here thus is likely still relevant in behavioral contexts, although over shorter timescales (and thus higher frequencies), but this remains to be explicitly tested.

Encoding models of early sensory processing

General statistical techniques have long been used to model the functional properties of neurons in early sensory pathways (Jones and Palmer, 1987; Reid et al., 1997). Models of spiking activity of neurons in the retina, the lateral geniculate nucleus (LGN) of the thalamus, and cortical simple cells all rely on predominantly linear filtering of the sensory input in the prediction of the neuronal response passed through a static, rectifying nonlinearity (Chander and Chichilnisky, 2001; Mante et al., 2005; Lesica et al., 2007; Mante et al., 2008; Pillow et al., 2008). It has also previously been acknowledged that spike-history plays a significant role in the patterns of neural activity observed experimentally (Brillinger, 1992; Berry and Meister, 1998), and recent modeling efforts have extended the quasi-linear framework to capture absolute and relative refractoriness in the neuronal response of retinal ganglion cell activity (Keat et al., 2001; Pillow et al., 2008). The nonlinear history dependence we describe here is a network phenomenon, which we have previously shown to occur independently of the spiking of the measured neuron and is thus predicted directly from the stimulus input rather than the past history of spiking (Webber and Stanley, 2004). Further, the model presented here was constructed from systematically measured tuning properties, and uses filter properties that are fundamentally nonlinear in nature that are based on the propagation of the suppressive effects that we have previously reported (Webber and Stanley, 2004, 2006; Boloori and Stanley, 2006).

It has recently been shown that, unlike neurons in the early visual pathway (including simple cells in V1), when driving the system with white noise vibrissa stimulation, spike-triggered averaging (STA) does not reveal a first order relationship between the vibrissa motion and the cortical response (Maravall et al., 2007). However, consideration of the second-order structure of the spike-triggered distribution of the stimulus does reveal stimulus features to which the system is sensitive. Estimation of the input-output relationship relative to these stimulus features reveals a symmetry of the corresponding nonlinearities, capturing the responsiveness of a barrel cortex neuron to vibrissa motion in opposite directions. Similar work in the early visual pathway suggests that the covariance methods provides complex feature selectivity (Fairhall et al., 2006), so it is possible that this approach would be predictive of the nonlinear interactions we demonstrate here on an event-by-event basis. It is also possible that STC-based models cannot accurately predict cortical responses to relatively strong stimuli. There is good reason to believe that the stimuli used in studies of other sensory pathways (e.g., spatiotem-poral white noise in visual studies) are not strong enough to push the system into a strongly nonlinear regime, and that the punctate deflections of the vibrissae are strong in comparison. This would also explain why the linear models of the early visual pathway perform poorly in the context of natural visual scenes, where strong spatial correlations, and sudden transients, drive sparse, episodic firing at various stages of processing (Butts et al., 2007), suggesting the possible existence of similarly strong, suppressive mechanisms as those described here.

Cortical Representations of the Tactile Scene

Representations of the tactile input in ethological conditions would be distributed across vibrissae and across the cortical population, ultimately placing limits on what we might expect to uncover through activity of a single neuron. For example, it is impossible for one to recover information on stimulus direction from observation of single cell responses. However, we might envision downstream neurons receiving input from a population of neurons within a cortical column with a diversity of directional tuning properties, likely to collectively provide unambiguous information concerning stimulus direction. Such a population code would be more easily envisioned if the directional tuning map were arranged in a manner that would facilitate such a readout, as recently suggested in several studies (Bruno et al., 2003; Lee and Simons, 2004; Andermann and Moore, 2008). Secondly, consideration of responses from multiple barrels (perhaps as are available to supra-and infra-granular SI regions and higher cortical areas) would also yield information on stimulus direction: the order in which two adjacent vibrissae are deffected contains information about the direction of contact. Population coding in the vibrissa pathway has not been well studied until recently (Petersen and Diamond, 2000; Panzeri et al., 2003; Petersen et al., 2003), but likely plays an important role in the representation of more complex tactile sensations. The experiments and analyses here also provide an opportunity to speculate upon the possible perceptual implications of the described phenomena. In a passive setting, it would trivially be the case, for example, that a texture with a high spatial frequency presented at a low velocity would be indistinguishable from a texture with a low spatial frequency presented at a high velocity. Beyond this simple example, however, it is important to note that any past history of sensory input that places the system in the same “state” would have the same neural response to the same subsequent sensory input, despite vast differences in the stimulus history. It is possible that manipulating the stimulus history prior to a cued stimulus used in a discrimination task would potentially provide some clues about the time window over which the past stimulus history influences the percept (Gerdjikov et al., 2010), but this remains to be fully explored experimentally.


We would like to thank Daniel Butts for comments at various points of this work. This work was supported by the National Institutes of Health (NIH R01NS48285).


  • Ahissar E, Sosnik R, Haidarliu S. Transformation from temporal to rate coding in a somatosensory thalamocortical pathway. Nature. 2000;406:302–306. [PubMed]
  • Ahissar E, Sosnik R, Bagdasarian K, Haidarliu S. Temporal frequency of whisker movement. II. laminar organization of cortical representations. Journal of Neurophysiology. 2001;86:354–367. [PubMed]
  • Andermann ML, Moore CI. A somatotopic map of vibrissa motion direction within a barrel column. Nat Neurosci. 2008;9:543–551. [PubMed]
  • Andermann ML, Ritt J, Neimark MA, Moore CI. Neural correlates of vibrissa resonance: Band-pass and somatotopic representation of high-frequency stimuli. Neuron. 2004;42:451–463. [PubMed]
  • Arabzadeh E, Zorzin E, Diamond ME. Neuronal encoding of texture in the whisker sensory pathway. PLoS Biology. 2005;3:155–165. [PMC free article] [PubMed]
  • Arabzadeh E, Panzeri S, Diamond ME. Deciphering the spike train of a sensory neuron: counts and temporal patterns in the rat whisker pathway. Journal of Neuroscience. 2006;26:9216–9226. [PubMed]
  • Bair W, Koch C. Temporal precision of spike trains in extrastriate cortex of the behaving macaque monkey. Neural Comput. 1996;8:1185–1202. [PubMed]
  • Berry MJ, Meister M. Refractoriness and neural precision. J Neurosci. 1998;18:2200–2211. [PubMed]
  • Boloori AR, Stanley GB. The dynamics of spatiotemporal response integration in the somatosensory cortex of the vibrissa system. J Neurosci. 2006;26:3767–3782. [PubMed]
  • Brillinger DR. Nerve cell spike train data analysis: A progression of technique. Journal of the American Statistical Association. 1992;87:260–271.
  • Brumberg JC, Pinto DJ, Simons DJ. Spatial gradients and inhibitory summation in the rat whisker barrel system. Journal of Neurophysiology. 1996;76:130–140. [PubMed]
  • Bruno RM, Sakmann B. Cortex is eriven by weak but synchronously active thalamocortical synapses. Science. 2006;312:1622–1627. [PubMed]
  • Bruno RM, Khatri V, Land PW, Simons DJ. Thalamocortical angular tuning domains within individual barrels of rat somatosensory cortex. J Neurosci. 2003;23:9565–9574. [PubMed]
  • Buracas GT, Zador AM, DeWeese MR, Albright TD. Effcient discrimination of temporal patterns by motion-sensitive neurons in primate visual cortex. Neuron. 1998;20:959–969. [PubMed]
  • Butts DA, Weng C, Jin J, Yeh CI, Lesica NA, Alonso JM, Stanley GB. Temporal precision in the neural code and the time scales of natural vision. Nature. 2007;449:92–95. [PubMed]
  • Cariani PA, Delgutte B. Neuronal correlates of the pitch of complex tones. J Neurophysiol. 1996;76:1698–1716. [PubMed]
  • Castro-Alamancos MA. Absence of rapid sensory adaptation in neocortex during information processing states. Neuron. 2004;41:455–464. [PubMed]
  • Chander D, Chichilnisky EJ. Adaptation to temporal contrast in primate and salamander retina. J Neurosci. 2001;21:9904–9916. [PubMed]
  • Chung S, Li X, Nelson SB. Short-term depression at thalamocortical synapses contributes to rapid adaptation of cortical sensory responses in vivo. Neuron. 2002;34:437–446. [PubMed]
  • Dayan P, Abbott LF. Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems. MIT Press; 2001.
  • Desbordes G, Weng C, Jin J, Lesica NA, Stanley GB, Alonso JM. Timing precision in population coding of natural scenes in the early visual system. PLOS Biology. 2008;6:2672–2682. [PMC free article] [PubMed]
  • Desîlets-Roy B, Varga C, Lavallée P, Deschênes M. Substrate for cross-talk inhibition between thalamic barreloids. The Journal of Neuroscience. 2002;22:RC218. [PubMed]
  • DeWeese MR, Zador AM. Non-gaussian membrane potential dynamics imply sparse, synchronous activity in auditory cortex. J Neurosci. 2006;22:12206–12218. [PubMed]
  • Fairhall AL, Burlingame CA, Narasimhan R, Harris RA, Puchalla JL, Berry MJ. Selectivity for multiple stimulus features in retinal ganglion cells. J Neurophysiol. 2006;96:2724–2738. [PubMed]
  • Fanselow EE, Nicolelis MAL. Behavioral modulation of tactile responses in the rat somatosensory system. The Journal of Neuroscience. 1999;19:7603–7616. [PubMed]
  • Garabedian CE, Jones SR, Merzenich MM, Dale A, Moore CI. Band-pass response properties of rat SI neurons. Journal of Neurophysiology. 2003;90:1379–1391. [PubMed]
  • Gerdjikov TV, Bergner CG, Stüuttgen MC, Waiblinger C, Schwarz C. Discrimination of vibrotactile stimuli in the rat whisker system: behavior and neurometrics. Neuron. 2010;65:530–540. [PubMed]
  • Green DM, Swets JA. Signal Detection and Psychophysics. New York, NY: John Wiley & Sons, Inc; 1966.
  • Heiss JE, Katz Y, Ganmor E, Lampl I. Shift in the balance between excitation and inhibition during sensory adaptation of S1 neurons. J Neurosci. 2008;28:13320–13330. [PubMed]
  • Higley MJ, Contreras D. Nonlinear integration of sensory responses in the rat barrel cortex: An intracellular study in vivo. The Journal of Neuroscience. 2003;23:10190–10200. [PubMed]
  • Higley MJ, Contreras D. Balanced excitation and inhibition determine spike timing during frequency adaptation. The Journal of Neuroscience. 2006;26:448–457. [PubMed]
  • Higley MJ, Contreras D. Frequency adaptation modulates spatial integration of sensory responses in the rat whisker system. J Neurophys. 2007;97:3819–3824. [PubMed]
  • Ito M. Processing of vibrissa sensory information within the rat neocortex. J Neurophys. 1985;54:479–488. [PubMed]
  • Jadhav S, Wolfe J, Feldman D. Sparse temporal coding of elementary tactile features during active whisker sensation. Nat Neurosci. 2009;12:792–800. [PubMed]
  • Jenks RA, Vaziri A, Boloori AR, Stanley GB. Self-motion and the shaping of sensory signals. J Neurophysiol. 2010 in press. [PubMed]
  • Jones JP, Palmer LA. The two-dimensional spatial structure of simple receptive ?elds in cat striate cortex. J Neurophysiol. 1987;58:1187–1211. [PubMed]
  • Keat J, Reinagel P, Reid RC, Meister M. Predicting every spike: a model for the responses of visual neurons. Neuron. 2001;30:803–817. [PubMed]
  • Khatri V, Hartings J, Simons D. Adaptation in thalamic barreloid and cortical barrel neurons to periodic whisker deflections varying in frequency and velocity. Journal of Neurophysiology. 2004;92:3244–3254. [PubMed]
  • Kleinfeld D, Delaney KR. Distributed representation of vibrissa movement in the upper layers of somatosensory cortex revealed with voltage-sensitive dyes. Journal of Comparative Neurology. 1996;108:375–389. [PubMed]
  • Knutsen P, Derdikman D, Ahissar E. Tracking whisker and head movements in unrestrained behaving rodents. Journal of Neuroscience. 2005;93:2294–2301. [PubMed]
  • Kyriazi HT, Carvell GE, Simons DJ. Off response transformations in the whisker/barrel system. Journal of Neurophysiology. 1994;72:392–401. [PubMed]
  • Lee SH, Simons DJ. Angular tuning and velocity sensitivity in different neuron classes within layer 4 of rat barrel cortex. Journal of Neurophysiology. 2004;91:223–229. [PubMed]
  • Lesica NA, Jin JZ, Weng C, Yeh CI, Butts DA, Stanley GB, Alonso JM. Adaptation to stimulus contrast and correlations during natural visual stimulation. Neuron. 2007;55:479–491. [PMC free article] [PubMed]
  • Lewicki MS. A review of methods for spike sorting: the detection and classification of neural action potentials. Network: Computation in Neural Systems. 1998;9:R53–R78. [PubMed]
  • Liu RC, Tzonev S, Rebrik S, Miller KD. Variability and information in a neural code of the cat lateral geniculate nucleus. J Neurophysiol. 2001;86:2789–2806. [PubMed]
  • Mainen ZF, Sejnowski TJ. Reliability of spike timing in neocortical neurons. Science. 1995;268:1503–1506. [PubMed]
  • Mante V, Frazor RA, Bonin V, Geisler WS, Carandini M. Independence of luminance and contrast in natural scenes and in the early visual system. Nat Neurosci. 2005;8:1690–1697. [PubMed]
  • Mante V, Bonin V, Carandini M. Functional mechanisms shaping lateral geniculate responses to artificial and natural stimuli. Neuron. 2008;58:625–638. [PubMed]
  • Maravall M, Petersen RS, Fairhall AL, Arabzadeh E, Diamond ME. Shifts in coding properties and maintenance of information transmission during adaptation in barrel cortex. PLoS Biology. 2007;5:e19. [PMC free article] [PubMed]
  • Marmarelis PZ, Marmarelis VZ. Analysis of Physiological Systems. New York, New York, USA: Plenium Press; 1978.
  • McCormick DA, Connors BW, Lighthall JW, Prince DA. Comparative electrophysiology of pyramidal and sparsely spiny stellate neurons of the neocortex. Journal of Neurophysiology. 1985;54:782–806. [PubMed]
  • Mirabella G, Battison S, Diamond ME. Integration of multiple-whisker inputs in rat somatosensory cortex. Cerebral Cortex. 2001;11:164–170. [PubMed]
  • Nicolelis MAL, Fanselow EE. Thalamocortical optimization of cortical processing according to behavioral state. Nat Neurosci. 2002;5:517–523. [PubMed]
  • Paninski L, Pillow JW, Simoncelli EP. Maximum likelihood estimation of a stochastic integrate-and-fire neural encoding model. Neural Computation. 2004;16:2533–2561. [PubMed]
  • Panzeri S, Petroni F, Petersen RS, Diamond ME. Decoding neuronal population activity in rat somatosensory cortex: Role of columnar organization. Cerebral Cortex. 2003;13:45–52. [PubMed]
  • Passaglia CL, Troy JB. Information transmission rates of cat retinal ganglion cells. J Neurophysiol. 2004;91:1217–1229. [PubMed]
  • Paxinos G, Watson C. The Rat Brain in Stereotaxic Coordinates. 4 New York: Academic Press; 1998.
  • Petersen CCH, Grinvald A, Sakmann B. Spatiotemporal dynamics of sensory responses in layer 2/3 of rat barrel cortex measured in vivo by voltage-sensitive dye imaging combined with whole-cell voltage recordings and neuron reconstructions. The Journal of Neuroscience. 2003;23:1298–1309. [PubMed]
  • Petersen RS, Diamond ME. Spatial-temporal distribution of whisker-evoked activity in rat somatonsensory cortex and the coding of stimulus location. The Journal of Neuroscience. 2000;20:6135–6143. [PubMed]
  • Pillow JW, Shlens J, Paninski L, Sher A, Litke AM, Chichilnisky EJ, Simoncelli EP. Spatiotemporal correlations and visual signaling in a complete neuronal population. Nature. 2008;454:995–999. [PMC free article] [PubMed]
  • Pinto DJ, Brumberg JC, Simons DJ. Circuit dynamics and coding strategies in rodent somatosensory cortex. Journal of Neurophysiology. 2000;83:1158–1166. [PubMed]
  • Reid RC, Victor JD, Shapley RM. The use of m-sequences in the analysis of visual neurons: Linear receptive ?field properties. Vis Neurosci. 1997;14:1015–1027. [PubMed]
  • Reinagel P, Reid RC. Temporal coding of visual information in the thalamus. J Neurosci. 2000;20:5392–5400. [PubMed]
  • Ritt JT, Andermann ML, Moore CI. Embodied information processing: vibrissa mechanics and texture features shape micromotions in actively sensing rats. Neuron. 2008;57:599–613. [PubMed]
  • Shoykhet M, Doherty D, Simons DJ. Coding of deflection velocity and amplitude by whisker primary afferent neurons: implications for higher level processing. Somatosens Mot Res. 2000;17:171–180. [PubMed]
  • Simons DJ. Response properties of vibrissa units in rat SI somatosensory neocortex. Journal of Neurophysiology. 1978;41:798–820. [PubMed]
  • Simons DJ. Temporal and spatial integration in the rat SI vibrissa cortex. Journal of Neurophysiology. 1985;54:615–635. [PubMed]
  • Simons DJ, Carvell GE. Thalamocortical response transformation in the rat vibrissa/barrel system. Journal of Neurophysiology. 1989;61:311–330. [PubMed]
  • Stoelzel CR, Bereshpolova Y, Swadlow HA. Stability of thalamocortical synaptic transmission across awake brain states. J Neurosci. 2009;29:6851–6859. [PMC free article] [PubMed]
  • Swadlow HA, Gusev AG. The impact of ‘bursting’ thalamic impulses at a neocortical synapse. Nat Neurosci. 2001;4:402–408. [PubMed]
  • Usrey WM, Alonso JM, Reid RC. Synaptic interactions between thalamic inputs to simple cells in cat visual cortex. J Neurosci. 2000;20:5461–5467. [PubMed]
  • Vincent SB. The function of the vibrissae in the behavior of the white rat. Behav Monographs. 1912;1:7–85.
  • von Heimendahl M, Itskov PM, Arabzadeh E, Diamond ME. Neuronal activity in rat barrel cortex underlying texture discrimination. PLOS Biology. 2007;5:2696–2708. [PMC free article] [PubMed]
  • Webber RM, Stanley GB. Nonlinear encoding of tactile patterns in the barrel cortex. J Neurophysiol. 2004;91:2010–2022. [PubMed]
  • Webber RM, Stanley GB. Transient and steady-state dynamics of cortical adaptation. Journal of Neurophysiology. 2006;95:2923–2932. [PubMed]
  • Wehr M, Zador AM. Balanced inhibition underlies tuning and sharpens spike timing in auditory cortex. Nature. 2003;27:442–446. [PubMed]
  • Wolfe J, Hill DN, Pahlavan S, Drew PJ, Kleinfeld D, Feldman DE. Texture coding in the rat whisker system: slip-stick versus differential resonance. PLOS Biol. 2008;6:1661–1677. [PMC free article] [PubMed]
  • Zhu JJ, Connors BW. Intrinsic firing patterns and whisker-evoked synaptic responses of neurons in the rat barrel cortex. Journal of Neurophysiology. 1999;81:1171–1183. [PubMed]