Abstract
Somatosensory coding in rodents has been mostly studied in the whisker system and hairy skin, whereas the function of low-threshold mechanoreceptors (LTMRs) in the rodent glabrous skin has received scant attention, unlike in primates where the glabrous skin has been the focus. The relative activation of different LTMR subtypes carries information about vibrotactile stimuli, as does the rate and temporal patterning of LTMR spikes. Rate coding depends on the probability of a spike occurring on each stimulus cycle (reliability), whereas temporal coding depends on the timing of spikes relative to the stimulus cycle (precision). Using in vivo extracellular recordings in male rats and mice of either sex, we measured the reliability and precision of LTMR responses to tactile stimuli including sustained pressure and vibration. Similar to other species, rodent LTMRs were separated into rapid-adapting (RA) or slow-adapting based on their response to sustained pressure. However, unlike the dichotomous frequency preference characteristic of RA1 and RA2/Pacinian afferents in other species, rodent RAs fell along a continuum. Fitting generalized linear models to experimental data reproduced the reliability and precision of rodent RAs. The resulting model parameters highlight key mechanistic differences across the RA spectrum; specifically, the integration window of different RAs transitions from wide to narrow as tuning preferences across the population move from low to high frequencies. Our results show that rodent RAs can support both rate and temporal coding, but their heterogeneity suggests that coactivation patterns play a greater role in population coding than for dichotomously tuned primate RAs.
Significance Statement
Our sense of touch starts with activation of nerve fibers in the skin. Although response properties of various fiber types are well established in other species (e.g., primates), quantitative characterization in rats and mice is limited. To fill this gap, we performed a comprehensive electrophysiological investigation into the coding properties of tactile fibers in rodent nonhairy skin and then simulated these fibers to explain differences in their responses. We show that rodent tactile fibers resemble those from other species but that their heterogeneity at the population level may differ, with potentially important implications for encoding of touch. Simulations reveal intrinsic mechanisms that support this heterogeneity and provide a useful tool to explore somatosensation in rodents.
Introduction
Touch processing begins when low-threshold mechanoreceptors (LTMRs) in the skin convert mechanical stimuli into action potentials, or spikes, that are relayed centrally. Classic work classified LTMRs as slow-adapting (SA) or rapid-adapting (RA) based on responses to sustained pressure (Adrian, 1928). Subsequent work, mostly in primate glabrous skin, subclassified SAs and RAs based on preferred frequencies. SAs spike repetitively during sustained pressure and are most sensitive to low-frequency vibrations (Talbot et al., 1968; Freeman and Johnson, 1982; Muniak et al., 2007), whereas RAs spike transiently during sustained pressure and are subdivided into RA1s, which are most sensitive to <50 Hz vibrations perceived as flutter, and RA2/PCs (Pacinians), which are most sensitive to >100 Hz vibrations perceived as continuous vibration (Johnson, 2001). That said, the stimulus waveform is also important to consider when eliciting RA responses (Birznieks et al., 2019; Bhattacharjee et al., 2021).
Work in rats (Sanders and Zimmermann, 1986; Lewin and McMahon, 1991; Na et al., 1993; Leem et al., 1993a, b; Guzun et al., 2021) and mice (Cain et al., 2001; Heidenreich et al., 2011; Walcher et al., 2018; Prsa et al., 2019; Neubarth et al., 2020; Schwaller et al., 2021) has confirmed that rodent LTMRs resemble LTMRs in other species and has advanced our molecular understanding of somatosensation (for reviews, see Handler and Ginty, 2021; Meltzer et al., 2021). But the coding properties of rodent LTMRs have not been as thoroughly quantified as in primates, cats (Horch et al., 1977; Iggo and Ogawa, 1977), or raccoons (Pubols et al., 1971). The limited quantitative characterization of vibrotactile responses in rodents shows marked variability, especially in RAs (Leem et al., 1993a; Neubarth et al., 2020). Furthermore, despite emerging evidence of the importance of spike timing for tactile sensation (Mackevicius et al., 2012; Pruszynski and Johansson, 2014; Saal et al., 2016; Birznieks and Vickery, 2017; Sukumar et al., 2022), this has not yet been examined in the rodent glabrous skin. Technical advances have furthered the use of rodent models for studying touch (for reviews, see Basbaum et al., 2009; Handler and Ginty, 2021), in which the glabrous skin, in particular, is stimulated to assess pain (Chaplan et al., 1994; Deuis et al., 2017; Dedek et al., 2023). Therefore, a comprehensive investigation into the coding properties of LTMRs in the rodent glabrous skin is vital for our understanding of somatosensation and its pathological disruption.
An LTMR's response to vibration is influenced by stimulus intensity and frequency. Previous studies suggested that stimulus intensity is rate-coded and frequency is temporally coded (Talbot et al., 1968; Johnson, 1974; Jones et al., 2004; Butts et al., 2007; Muniak et al., 2007; Bensmaia, 2008; Mackevicius et al., 2012; Weber et al., 2013; Graczyk et al., 2016). Temporal coding depends on spike timing relative to the stimulus phase (i.e., spike precision), whereas rate coding reflects the probability of spiking on each stimulus cycle (i.e., spike reliability). Reliability and precision are both important for spike synchronization across neurons (Mainen and Sejnowski, 1995), which is important for information filtering (Grewe et al., 2017) and sensation (Bruno, 2011; Sagalajev et al., 2024). Examining the reliability and precision of LTMR responses is crucial for understanding somatosensory coding at a single-neuron level and how neurons work together at the population level.
Using rats and mice, we quantified the reliability and precision of LTMR responses to hind paw stimulation in vivo. Clustering based on reliability revealed that rodent RAs appear to form a continuum of frequency preferences. Conversely, increasing stimulus frequency universally increased spike precision. We confirmed that poorer precision at low frequencies was due to covariation of the stimulus waveform with frequency during sinusoidal stimulation. Generalized linear models (GLMs) fit to experimental data reproduced RA reliability and precision. Fitted parameters highlight key mechanistic differences supporting the spectrum of frequency preferences. Together, our results show that LTMRs in the rodent glabrous skin resemble those from other species, but their heterogeneity at the population level may differ, with potentially important implications for somatosensory coding.
Materials and Methods
In vivo experiments
All procedures were approved by the Animal Care Committee of The Hospital for Sick Children (protocol numbers 53451 and 47072) and were conducted in accordance with guidelines from the Canadian Council on Animal Care. Adult male Sprague Dawley rats (200–350 g) were obtained from Charles River. Male and female mice (10–12 weeks, 20–30 g) with a C57BL/6 background expressing channelrhodopsin-2 in sensory afferents were produced by crossing AdvillinCre/+ mice (kindly provided by Fan Wang, Duke University) or PVCre/+ mice (JAX:017320) with Ai32 mice (JAX:024109), but optogenetic stimuli were not tested in this study. Rats and mice were anesthetized with 1.2 g/kg urethane intraperitoneally and received 10% top-ups as needed to block the pedal reflex. Laminectomies were performed to expose the L4 and L5 dorsal root ganglion (DRG) for single-unit extracellular recordings as previously described in rats (Sagalajev et al., 2024) and mice (Al-Basha and Prescott, 2019). For rats, we used a multielectrode array with 16 contacts on four shanks (NeuroNexus, A4 type). For mice, we used a parylene-insulated tungsten microelectrode (A-M Systems, 573500). Signals were amplified, digitized at 40 kHz, and high-pass filtered at 300 Hz using an OmniPlex Data Acquisition System (Plexon). Gentle brushing was used to search for LTMRs innervating the glabrous skin of the hind paw. The hind paw, as opposed to the fore paw, was stimulated to remain consistent with rodent literature (Leem et al., 1993a; Chaplan et al., 1994; Cain et al., 2001; Deuis et al., 2017; Neubarth et al., 2020; Handler and Ginty, 2021). Receptive fields were then mapped using von Frey filaments. Once a unit was identified, precise and reproducible stimuli were applied with a flat-tipped indenter (1 mm diameter) using a computer-controlled mechanical stimulator (Aurora Scientific, model 300C-I) positioned with a micromanipulator. Stimuli included sustained pressure steps (25–225 mN), sinusoidal vibration (2–400 Hz, 50–225 mN), and trains of brief pulses (10–50 pps, 40% duty cycle, 50–225 mN) each applied for 1 s for each frequency–amplitude combination. Data in the 150–225 mN range were used for principal component analysis (PCA) and GLM fitting because spiking responses at these intensities were consistently robust (i.e., resulted in a large enough spike count to perform analysis on) across the 2–400 Hz frequency range. Stimulus intensity was defined here as the peak-to-peak amplitude of force variations. All stimuli were controlled by a Power1401 A-D board and Signal v5 software (Cambridge Electronic Design). Stimulus timing was sent as a trigger to Plexon to align neural responses to stimuli. Single units were isolated using the Offline Sorter v4 software (Plexon) and analyzed with MATLAB (MathWorks, R2021).
Spike train analysis
Reliability and precision of LTMR responses to periodic tactile stimuli were measured using entrainment and temporal dispersion, respectively. Entrainment was calculated as the number of spikes per cycle. A neuron spiking once on each stimulus cycle was considered 1:1 entrained with the stimulus. We plotted entrainment against stimulus frequency to identify the preferred frequencies (tuning) of each neuron. Temporal dispersion is inversely related to spike timing precision and, for periodic stimuli, is also related to phase-locking. Phase-locking is typically expressed in radians or degrees, but here we convert it to absolute time by normalizing by the stimulus period using the following formula:
Cluster analysis
For an unbiased characterization of RA frequency preferences, K-means clustering was performed on RA tuning curves. For the clustering and PCA, the features were individual stimulus frequency–intensity combinations (e.g., 10 Hz–150 mN), and the samples were the corresponding entrainment for each feature. The optimal number of clusters for each K-value (between 2 and 5) was determined using the mean silhouette value. Z-Scores were calculated for all features (based on the mean and standard deviation), and PCA was performed to visualize the high-dimensional clustering results. The first two PCs accounted for 55.3% of the total variance in the data (PC1, 34.7%; PC2, 20.6%). Clustering was first performed on rat data, and then mouse data were subsequently assigned to existing clusters.
Generalized linear models
The GLM structure uses code modified from Weber and Pillow (2017). GLMs have been used extensively to study sensory signaling (Truccolo et al., 2005; Pillow et al., 2008; Aljadeff et al., 2016). Briefly, model parameters include a stimulus filter (k), a postspike filter (h), and mean input (μ) that determines the baseline firing rate. GLMs were fit to the response of individual RAs (n = 20) to sinusoidal vibration across different combinations of frequency (5–400 Hz) and intensity (150, 175, 225 mN). To fit the GLM to experimental data, the stimulus and spike history (i.e., time-shifted spike trains) were each convolved with unique sets of basis functions. To allow for fine temporal structure near the spike, we used raised cosine basis vectors (bi) of the following form:
The output of the convolution of filters and basis vectors was summed and then put through a nonlinear exponential function (f = exp; corresponding to the inverse log link function) to produce an instantaneous firing rate (λ) as follows:
Filter analysis
To analyze the fitted stimulus (k) filter for the GLM, we calculated its width using a noise threshold (4× the STD of the filter at 40–50 ms before the spiking response) and then counted the total time that each filter's absolute value was above that threshold. We also measured each filter's excitation/suppression ratio by dividing the sum of the positive and negative components of each filter by its overall magnitude (sum of its absolute value over time). A power spectrum analysis was also performed for each stimulus filter using a discrete Fourier transform (fft function, MATLAB) which was converted to a double-sided power spectrum (P2) using, P2 = | fft(k) / L |, where L is filter length of 50 ms. The average (±SD) single-sided power spectrum for k-filters was plotted across different clusters. From each individual power spectrum (one per neuron), we measured the center frequency (CF), the bandwidth (BW), and the Q-factor (Q = CF / BW) and compared them with RAs. A high-Q indicates a narrow passband, and a low-Q indicates a wide passband. To analyze the postspike filter, we measured its width and area under the curve (AUC). The width of the postspike filter was measured from time 0 to the time the filter first passed zero on the y-axis. This postspike filter width is analogous to the neuron's refractory period following a spike. AUC was calculated using trapezoidal numerical integration (trapz function, MATLAB) for the entire filter, as well as its positive (+) and negative (–) components.
Mutual information analysis
Mutual information (MI) was calculated as follows:
The input to the decoder, ri ϵ R, was a feature-mapped representation of LTMR activity over 100 ms. A sliding window was not used; instead, 1 s was split into 10 intervals of 100 ms. For single neurons, the feature map involved discretizing a sequence of interspike intervals from single trials into 10 bins, with edges logarithmically spaced between 1 and 100 ms, and using the counts in each to form a 10-dimensional input vector. Feature mapping for population activity was accomplished by summing the feature vectors of each neuron in that population. This mapping both improved classification performance and standardized the dimensionality of the input. While this simple summation does not allow for combinatorial decoding, alternative representations of population activity which preserve neuronal identity had dimensionalities too high for training purposes given our low number of samples per stimulus. In total, there were 1,280 training samples per decoder and 128 stimuli and 10 intervals of 100 ms per stimulus. Eighty percent of these samples were randomly used to train the decoder, while the remaining 20% were used as a test set. Following training, the change in entropy between P(S) and P(S|R = ri) was computed for each test example (i), and these values were averaged together to produce the final measure of information between S and R. A new decoder was initialized and trained for each of the 3,800 {S, R} pairs: 2 groups (heterogeneous vs homogeneous; see below) × 100 permutations × 19 populations.
More homogeneous or more heterogeneous populations of LTMRs were created to analyze population-level coding. To control for the sum of individual neuron informativeness within each type of population (homogeneous vs heterogeneous), LTMRs were split into eight groups of two or three neurons with similar levels of MI per neuron. A neuron is selected from each group to create a population of N neurons. Each population type starts with the same first neuron, but all subsequently added neurons are chosen from randomly ordered groups. Neurons were chosen from each group based on their similarity in frequency tuning to the first neuron, as defined by their PC2 value. Specifically, homogeneous populations were created by choosing neurons with similar PC2 values and heterogeneous populations were created by choosing neurons with dissimilar PC2 values. This process was repeated 100 times per population size (i.e., 100 permutations; see Fig. 8B).
Statistical analysis
All statistical analysis was performed using GraphPad Prism 9.5.1 (GraphPad Software). Parametric tests were used when data were normally distributed (according to the Shapiro–Wilk test). The comparison of two groups was done using an unpaired t test with Welch's correction. When comparing more than two groups, a one-way ANOVA was performed, followed by post hoc Tukey's tests when appropriate. Pearson's correlation coefficient (R) and p values are stated for all correlations. χ2 and p values are reported when comparing proportions.
Code availability
Code for models is available on ModelDB (https://modeldb.science/2016218).
Results
Rodent LTMRs respond heterogeneously to tactile stimuli
To characterize the responses of rodent LTMRs to tactile stimulation, we applied various tactile stimuli to the hind paw of rats and mice while recording from somata in the DRG. Like other species, rat LTMRs were classified as SA, (n = 4) or RA (n = 31) if they spiked repetitively or transiently, respectively, to steps of sustained pressure (Fig. 1A). None of the RAs exhibited baseline (spontaneous) activity without stimulation, whereas three of four SAs did. In response to sinusoidal vibration, SAs consistently entrained best at low frequencies (<50 Hz), often firing >1 spike per stimulus cycle (Fig. 1B, green), whereas RAs were more heterogeneous (Fig. 1Bi–iv). Some RAs were tuned to low frequencies (<50 Hz; example i), others were broadly tuned (examples ii and iii), and some were tuned to high frequencies (≥50 Hz; example iv). Some RAs exhibited 2:1, 3:1, or even 4:1 entrainment to vibrotactile stimuli, especially those afferents with low-frequency preferences (e.g., 2:1 entrainment to 5 Hz by the RA in Fig. 1Bi). Whereas the frequency sensitivity of entrainment varied across RAs (Fig. 1A, top row), temporal dispersion consistently decreased (middle row) and vector strength stayed high (bottom row) as stimulus frequency was increased. Similar response properties were observed in mouse RAs (Fig. 2). A typical response to low (50 Hz) and high (200 Hz) frequency stimulation with spike times plotted as latency from the cycle onset (Fig. 1C, top and middle) or cycle phase (Fig. 1C, bottom) shows that responses with similar vector strength (i.e., distribution based on phase) have different temporal dispersion (i.e., distribution based on time) once phase data are normalized by stimulus frequency (Eq. 1). We proceeded to characterize RA responses in more detail.
Heterogeneity in rat LTMR responses. A, Sample responses to sustained pressure (left) and sinusoidal vibration at 10, 50, and 100 Hz from an SA and four RAs (i–iv) for four stimulus intensities (color density). Rasters depict each spike as a dot. B, Entrainment (top), temporal dispersion (TD, middle), and vector strength (r, bottom) as a function of stimulus frequency for units depicted in A. Dashed line corresponds to 1 spike per cycle (or 1:1 entrainment). The four RAs shown here are representative neurons from the four clusters in Figure 3. C, Example of spike precision at 50 Hz (r = 0.98; TD = 0.69 ms) and 200 Hz (r = 0.96; TD = 0.22 ms) for a single neuron across cycles of a 1-s-long vibrotactile stimulus. Cumulative probability functions show that although mean and standard deviation of spike phase across cycles are similar at 50 and 200 Hz (indicating similar vector strength; phase, 147.8 ± 12.8° for 50 Hz vs 134.8 ± 16.0° for 200 Hz), the mean and standard deviation of spike latency across cycles is smaller at 200 Hz (indicating smaller temporal dispersion; latency, 8.21 ± 0.71 ms for 50 Hz vs 1.87 ± 0.22 ms for 200 Hz). This demonstrates how responses with similar vector strength in phase space can have different dispersion in the temporal space.
Heterogeneity in mouse RA responses. Entrainment (top) and temporal dispersion (bottom) for four stimulus intensities (color density) as a function of stimulus frequency for four RAs (one from each cluster from Fig. 3).
RAs exhibit a continuum of frequency preferences
Tuning curves revealed that individual RAs have heterogeneous responses to vibration (Fig. 1B) but, as a population, had preferred frequencies (exhibited 1:1 entrainment) spanning the entire range of experimentally tested frequencies (2–400 Hz). To classify rat RAs based on their response to vibration and preferred frequencies, we applied K-means clustering to their tuning curves (Fig. 3). Clustering based on individual stimulus intensities (Fig. 3A) revealed three dominant clusters (filled circles) and a subset of neurons that changed clusters depending on stimulus intensity (open circles). Clustering based on all intensities (Fig. 3Bi) revealed four clusters, including three (Clusters 1, 2, 4) like those in panel A and a fourth (Cluster 3) representing neurons whose association to the other three clusters varied with stimulus intensity. Mouse RAs (n = 7; Fig. 3Bi, triangles) fell within clusters of rat RAs and were assigned accordingly. Sample tuning curves for an RA from each cluster are shown in Figure 1B (for rat) and Figure 2 (for mice). According to PC loading analysis, PC1 appears to capture tuning curve width (more positive PC1 values indicate broader tuning), and PC2 appears to capture frequency preference (more positive PC2 values indicate preference for lower frequencies; Fig. 3Biii). The correlation matrix of all RA tuning curves (Fig. 3Bii) showed some evidence of clustering, but correlations often extended to neurons in different clusters, with positive correlations smearing across the diagonal. Together, the clustering and correlation analyses suggest that rather than clearly separable clusters, RAs fall along a continuum, especially when considered in a single dimension (i.e., tuning breadth “or” preferred frequency). Comparing additional RA properties revealed that clustering was independent of the receptive field location (Fig. 3C) and there was little difference in CV (Fig. 3D), the presence of onset/offset spikes (Fig. 3E), or the threshold for onset/offset spikes (Fig. 3F) across the four clusters. These findings affirm that RAs cannot be subdivided into distinct groups defined by unique properties but, instead, represent a heterogeneous population with a continuum of frequency preferences. Henceforth, we refer to the four clusters when reporting response properties, but we also plot data against PC2, which represents frequency preference.
Rat RAs do not fall into discrete subtypes but, instead, exhibit a continuum of frequency preferences. A, K-Means clustering of RAs based on individual stimulus intensities (175, 200, 225 mN) suggested three clusters (blue, purple, orange), but certain neurons (open circles) did not fall consistently into the same cluster and instead varied depending on stimulus intensity. Optimal number of clusters was determined from the peak in silhouette values (red diamonds, inset). B, (i) Clustering based on all stimulus intensities suggested four clusters. Centroids are marked with ×. Cluster 3 comprises RAs depicted as open circles in panel A. Mouse RAs (triangles) were assigned to appropriate rat clusters after projecting data on PC1 and PC2. PC1 and PC2 reflect tuning curve width and preferred frequency, respectively (see iii). Collapsing points along PC1 (top box inset) or PC2 (right box inset) and plotting their cumulative distributions (black lines) demonstrate the spectrum of RAs portrayed by each PC. (ii) Correlation matrix of all RA tuning curves ordered by PC2 value. The matrix shows some clustering (evident as large blocks), but correlations tend to extend beyond clusters, smearing along the diagonal. N44 and N21 sit on opposite ends of the PC2 spectrum. Correlations (R) range from −1 to 1. (iii) PC loading for PC1 and PC2. Mean (±SEM) PC load is calculated from all stimulus intensities at each stimulus frequency (see Materials and Methods). C, Receptive field centers for individual RAs. D, CV differed significantly between some of the clusters (F3,27 = 6.062; **p = 0.0027; one-way ANOVA). Specifically, CV of Cluster 1 (17.08 ± 1.65 m/s) was higher than Cluster 2 (27.06 ± 1.44 m/s; *p = 0.0071; Tukey's post hoc test), but not Cluster 3 (22.01 ± 2.50 m/s; p = 0.37) or Cluster 4 (25.56 ± 1.35 m/s; p = 0.066). E, Proportion of neurons spiking transiently at the onset or offset of sustained pressure steps did not differ significantly across clusters (onset, χ2 = 2.0701; p = 0.56; offset, χ2 = 6.1742; p = 0.10). F, Threshold for onset and offset spikes often differed within a given RA, but not with any pattern that was consistent across clusters.
Models reproduce RA responses to vibration and pressure
To further study RA properties, we fit GLMs to responses evoked by a subset of vibrotactile stimuli (Fig. 4A). All successfully fitted models (n = 20) reproduced responses to vibrotactile stimuli not included in the training dataset (Fig. 4B). The models for RAs illustrated in Figure 1i–iv are shown in Figure 4Bi–iv, respectively. The entrainment and temporal dispersion of simulated responses did not on average deviate significantly from experimental data (see insets). Like experimental data, model RAs exhibited heterogeneous tuning curves and temporal dispersion decreased with increasing stimulus frequency (Fig. 4B).
Fitted GLMs quantitatively reproduced RA response properties. A, Schematic of GLM components modified from Weber and Pillow (2017). B, Sample spike trains (left) and corresponding entrainment (middle) and temporal dispersion (right) for GLMs fitted to units i–iv in Figure 1B. Experimentally observed rasters (colored background) are shown for comparison against three spike trains simulated from the corresponding GLM. Low and high frequencies used to simulate the rasters were as follows: i, 30, 100 Hz; ii, 30, 100 Hz; iii, 40, 100 Hz; and iv, 40, 100 Hz. Insets show entrainment and temporal dispersion of simulated (predicted, y-axis) responses plotted against experimental (observed, x-axis) responses. Data points, each representing a different frequency–intensity combination, tend to fall along the diagonal indicating equivalence. Correlations for reliability (i, R = 0.91; ii, R = 0.67; iii, R = 0.87; iv, R = 0.83; ***all p values < 1.0 × 10−8). Correlations for temporal dispersion (i, R = 0.80; ***p < 1.0 × 10−8; ii, R = −0.14; p = 0.29; iii, R = 0.26; *p = 0.044; iv, R = 0.69; ***p < 1.0 × 10−8). C, Threshold for onset and offset spikes in model RAs. Models divided into two groups, ON (circles) or OFF (triangles), based on spiking at the onset or offset to sustained pressure. There was no consistent alignment between ON/OFF model grouping and Clusters 1–4 (from Fig. 3) except for Cluster 1 RAs whose corresponding models were all ON. Deviation above or below diagonal (dashed line) indicates a preference (lower threshold) for onset or offset spiking, respectively, in experimental data from corresponding RAs. Models accurately predict this preference (i.e., ON models had lower experimental onset threshold and OFF models had lower offset threshold). D, Average stimulus filter for ON (black) and OFF (gray) models (shading = SEM). OFF models had a more negative filter (quantified as net AUC) compared with ON models (***p = 0.0004; Welch's t test), consistent with their preference for decreases in stimulus force (i.e., stimulus offset).
To further test the generalizability of our models, we also simulated their responses to sustained pressure to ensure they could capture the onset/offset spiking response that is quintessential to RAs. Models responded at the stimulus onset or offset (Fig. 4C) and are labeled ON or OFF, respectively; OFF models also responded at the onset of high-amplitude steps. Experimental data confirmed that neurons fitted with OFF models had a lower threshold for offset spikes than for onset spikes (i.e., prefer stimulus offset), whereas the opposite was true for neurons fitted with ON models. Since responses to sustained pressure were not included in the training dataset, reliable reproduction of these responses demonstrate that our models can accurately predict a neuron's preference for the stimulus onset or offset. Indeed, analysis of model stimulus filters revealed that OFF neurons have a more negative filter (quantified as AUC) compared with ON neurons (Fig. 4D; p = 0.00040) which supports their preference for downward deflections in stimulus position (i.e., the stimulus offset). Together, these simulations confirm that our models accurately reproduce RA responses to diverse tactile stimuli.
Stimulus filters capture the continuum of frequency preferences in RAs
Next, we analyzed RA filters to explore the mechanisms underlying heterogeneous frequency tuning. We hypothesized that the stimulus filter, which describes a neuron's preferred stimulus feature, would parallel the preferred frequency as determined from the entrainment tuning curves (Fig. 1). As predicted, ordering of RA stimulus filters according to experimental PC2 values revealed a continuum of filter widths (Fig. 5A). Visual inspection of the average stimulus filters demonstrated a transition from wide to narrow as tuning preference across the RA population moved from low to high frequencies (i.e., from positive to negative PC2 values; Fig. 5B). When formally measured, the width of individual stimulus filters is correlated with PC2 (R = 0.47; *p = 0.038; Fig. 5Ci) as is the excitation/suppression ratio of RA filters (R = 0.64; *p = 0.0024; Fig. 5Cii). The latter observation suggests that neurons with higher or lower-frequency preferences have more suppressive or excitatory components to their stimulus filters, respectively. Comparing our results in rodents with the RA/PC stimulus filters found in primates (Saal et al., 2015), we see strong similarities in the range of filter widths (from ∼10 to 60 ms) where neurons with higher-frequency preferences (i.e., RAs with low PC2 in rodents or PCs in primates) typically have narrower stimulus filters. Furthermore, in both datasets, neurons with low-frequency preferences (i.e., RAs with higher PC2 in rodents or RAs in primates) typically had filters with more excitatory (positive) components than those that had high-frequency preferences according to their excitation/suppression ratios. To demonstrate the role of the stimulus filter in determining a neuron's frequency preference, we swapped a narrow (N17, orange) and wide (N12, blue) stimulus filter which caused the model to transition from high- to low-frequency preference (Fig. 5D).
Stimulus filters of model RAs capture their spectrum of frequency preferences. A, Stimulus filters of RA models are shown aligned from smallest to largest PC2 value. Filters transition from wide to narrow as PC2 values move from positive to negative (i.e., as frequency preference shifts from low to high). B, Mean stimulus filters for each cluster (shading = SEM). C, The width (i; R = 0.47; *p = 0.038) and excitation/suppression (Exc./Supp.) ratio (ii; R = 0.64; *p = 0.0024) of stimulus filters in A are positively correlated with PC2. These data confirm that RAs with higher-frequency preferences have a narrower and more suppressive (negative) stimulus filter than neurons with a lower-frequency preference. Exc./Supp. ratio of 1, purely excitatory; 0, balanced; and −1, purely suppressive. D, Swapping a narrow stimulus filter (orange; N17) and wide stimulus filter (blue; N12) transitions the preference of the neuron from high (orange band) to low (blue band) frequencies. E, Average power spectrum of RA stimulus filters across each cluster. Peak power shifts from low frequency for RAs filters in Cluster 1 to high frequency for filters in Cluster 4. Shading = SEM. F, CF (i; R = −0.79; **p < 1.0 × 10−4) and Q-factor (ii; R = −0.69; **p < 1.0 × 10−3) of power spectra for individual neurons (from E) are negatively correlated with PC2. High-Q factor, narrow passband; low-Q factor, wide passband. This demonstrates that neurons with higher-frequency preferences (i.e., lower PC2 value) have a higher CF and narrower passband than neurons with a lower-frequency preference (i.e., a higher PC2 value). G, The BW of individual power spectra (from E) is positively correlated with PC1 (R = 0.51; *p = 0.02). This confirms that neurons with narrower tuning curves (i.e., lower PC1 values) have a narrower BW according to their stimulus filter.
To further quantify stimulus filter properties, we calculated the power spectrum for individual RA filters and averaged them across each cluster (Fig. 5E). Cluster 1 filters had peak power at lower frequencies (<50 Hz); Cluster 2 and 3 filters had broader, bimodal power spectra; and Cluster 4 filters had peak power at higher frequencies (>100 Hz). The CF of the power spectrum for each RA is correlated with the RA's PC2 value (R = −0.79; **p < 1.0 × 10−4; Fig. 5Fi), as is the Q-factor (R = −0.69; **p < 1.0 × 10−3; Fig. 5Fii), which represents the ratio of the spectrum's CF to its width. This finding demonstrates that RAs with a higher-frequency preference (based on their tuning curves) generally have stimulus filters with higher center frequencies and narrower passbands. Furthermore, plotting CF against PC2 demonstrated a gradual transition of peak filter power from high to low frequencies across PC2, which further confirms that the stimulus filter is capturing the spectrum of frequency preferences seen in RAs. As expected, we found that BW of the individual power spectra is correlated with PC1 (R = 0.51; *p = 0.020; Fig. 5G) but not with PC2 (R = −0.44; p = 0.053). In summary, analysis of model parameters revealed that cell-to-cell variability in the stimulus filter supports heterogeneity in RA tuning and further suggests that variations in stimulus filter width (i.e., length of the integration time window) underlie the continuum of RA frequency preferences.
Postspike filters capture the ability of RAs to spike repetitively
To separate the effects of waveform kinetics from the repetition rate on LTMR responses (since the maximum rate of change of force increases with sinusoid frequency), we also tested square pulses repeated at 10–50 pps. Unlike their response to low-frequency sinusoidal vibration (10–50 Hz), most RAs entrained 1:1 (or higher) to low-frequency pulse trains (10–50 pps) and did so with low temporal dispersion. Specifically, plotting entrainment and temporal dispersion for sinusoids versus pulse trains for the equivalent frequencies and intensities shows that entrainment was significantly higher (Fig. 6Ai; all ***p < 1.0 × 10−8) and temporal dispersion was significantly lower (Fig. 6Bi; C1, C2, C3, ***p < 1.0 × 10−8; C4, **p = 0.0082) for pulse trains. Comparing reliability with all neurons (except Cluster 4 neurons, which did not respond to sinusoids) revealed a positive correlation between PC2 and entrainment ratios (Fig. 6Aii; R = 0.39; **p < 1.0 × 10−10). This finding suggests that neurons with lower-frequency preferences (higher PC2 values) have similar responses (sine/pulse ratio ∼1) to low-frequency pulses and sinusoids, whereas neurons with higher-frequency preferences (lower PC2 values) exhibit higher entrainment (sine/pulse ratio <1) to pulses than to sinusoids. There was no correlation between PC2 and temporal dispersion ratios as precision was higher to pulses than sinusoids across all RAs (Fig. 6Bii; R = −0.025; p = 0.69). Whereas the previous trends (to sinusoids) mostly reflect differences in frequency tuning and stimulus filters across the RA population (Figs. 3 and 5), the trends in reliability and precision to pulse trains uniquely reflect the refractoriness of a neuron (i.e., how quickly it recovers from the last pulse in order to respond to the next pulse). In the GLM, refractoriness is reflected in the width of the postspike filter, which we analyzed next.
Response of RAs to pulse trains. A, Comparison of reliability to pulsed and sinusoidal vibration. i, Top, Each line depicts the mean entrainment for a unit averaged across four stimulus intensities. Dashed line represents 1:1 entrainment. i, Bottom, Entrainment is higher for pulse trains (yellow) than for sinusoidal stimulation (brown) at equivalent frequencies, as evident from data points (each representing a neuron's entrainment to pulse and sine stimuli at the same frequency and amplitude) sitting above the dashed line representing equivalence (all p values are ***p < 1.0 × 10−8; Wilcoxon signed-rank tests). ii, Top, Ratio of entrainment (sine/pulse) across all neurons (except Cluster 4 neurons, which did not respond to 10–50 Hz sinusoids). ii, Bottom, Positive correlation between entrainment ratio (sine/pulse) and PC2 (R = 0.39; **p < 1.0 × 10−10). B, Comparison of precision between pulses and sinusoids. i, Top, Each line depicts the mean temporal dispersion (TD) for a unit averaged across four stimulus intensities. i, Bottom, TD is lower for pulses than for sinusoids, as evident from data points (each representing a neuron's TD to pulse and sine stimuli at the same frequency and amplitude) sitting below the dashed line representing equivalence (C1, C2, and C3, ***p < 1.0 × 10−8; C4, **p = 0.0082; Wilcoxon signed-rank tests). ii, Top, Ratio of TD (sine/pulse) across all neurons (except Cluster 4 neurons, which did not respond to 10–50 Hz sinusoids). ii, Bottom, No correlation between TD ratio (sine/pulse) and PC2 (R = −0.025; p = 0.69).
Unlike the stimulus filter, the postspike filter of fitted GLMs did not change in a consistent manner across the frequency preference spectrum (Fig. 7). Ordering of postspike filters according to PC2 value revealed heterogeneity in this model parameter across the RA population (Fig. 7A). Furthermore, there was no significant correlation between postspike filter width and PC2 (R = −0.18; p = 0.47; Fig. 7B) nor between the postspike filter AUC and PC2 (−component, R = −0.01; p = 0.96; +component, R = 0.43; p = 0.07; Fig. 7C). To explore the role of the postspike filter on RA response properties, we simulated responses to square pulses repeated at different rates (Fig. 7D). Models entrained 1:1 with low-rate (≤50 pps) pulse trains. As the rate was increased, models could not maintain 1:1 entrainment as the interpulse interval dropped below their postspike filter width. For example, N12 (lighter blue) with a filter width of ∼20 ms could not entrain 1:1 with pulse trains >50 pps (i.e., cycle time <20 ms), whereas N39 (darker blue) had a shorter filter width of ∼10 ms and could entrain 1:1 at ≥100 pps. Both neurons had low entrainment at 200 pps because the cycle length (5 ms) was shorter than both their postspike filters.
Postspike filters of model RAs capture the ability to follow (i.e., entrain to) repeating stimuli. A, Postspike filters aligned by the neuron's PC2 value, from low- to high-frequency preference, show no consistent difference across clusters (i.e., heterogeneity is similar within and between clusters). Average postspike filters (±SEM, shading) are similar across Clusters 1–4. B, Postspike filter width (h-width) did not correlate with PC2 (R = −0.18; p = 0.47). C, Positive (+) and negative (−) components of postspike filter quantified as AUC did not correlate with PC2 (+AUC, R = 0.43; p = 0.07; −AUC, R = −0.01; p = 0.96). D, Width of postspike filter impacts model neuron's ability to entrain to pulse trains. N12 has a longer postspike filter width (∼20 ms, light blue) than N39 (∼10 ms, darker blue) and could not entrain 1:1 to pulse trains >50 pps (interpulse interval, 20 ms). N39 continues to entrain 1:1 at ≥100 pps (interpulse interval ≤10 ms). Neither neuron spikes reliability at 200 pps because the interpulse interval (5 ms) is shorter than their postspike filter widths.
Heterogeneous tuning improves the efficiency of population-level coding
Finally, to probe the functional consequences of heterogeneity in LTMR tuning, we compared the amount of MI on stimulus identity from either homogeneous or heterogeneous populations of LTMRs. This was done by estimating the lower bound of MI from ISI distributions and sums thereof for each population type. We hypothesized that, all else being equal, a heterogeneous population of RAs conveys more information about a stimulus than a homogeneous population because the individual neuron responses carrying that information are less redundant in the former case (Fig. 8A). In other words, if several neurons are coactivated by various stimuli, together those neurons will convey more information if their tuning is more different (less overlapping) than if their tuning is more similar (more overlapping). To test this, populations were created that either favored neurons with similar frequency preferences (homogeneous populations) or dissimilar ones (heterogeneous populations; Fig. 8B). As predicted, our analysis revealed that heterogeneous RA populations convey significantly more information about stimulus identity than homogeneous populations when populations contained between 6 and 17 neurons (Fig. 8C; *p < 0.001; #p < 0.01). Similarly, decoder accuracy was significantly higher for heterogeneous RA populations than homogeneous ones when the populations contained between 7 and 16 neurons (Fig. 8C, inset top). At most, the difference was 0.3 bits or 5% decoding accuracy. The amount of information converges because both populations were subsampled from the same set of 19 neurons. These results show that equivalent coding can be achieved using fewer neurons because of heterogeneous LTMR tuning.
Heterogeneous tuning improves the efficiency of population-level coding. A, Schematic outlining hypothesis that populations comprising an equal number of neurons with comparable informativeness per neuron (represented by circle area) will have greater MI on stimulus identity if constituent neurons have heterogeneous frequency preferences (yellow) than if they are more homogeneously tuned (green). Venn diagrams depict amount of overlapping information for each population type. The sum of neuron informativeness (circle area) must be controlled for because it is the upper bound on MI in a population code. B, Overview of procedure for sampling neurons that are more similar (homogeneous) or different (heterogeneous) based on frequency preference (PC2 value), when controlling for individual neuron informativeness. C, Heterogeneous populations (yellow) contain significantly more MI than homogeneous populations (green) comprising 6–17 neurons (left panel), which translates to significantly greater decoder accuracy (top-right panel; *p < 0.001; #p < 0.01; mean ± 2 SEM; n = 100 permutations; t tests with Bonferroni’s correction). The sum of neuron informativeness is equivalent for each population type (bottom-right panel).
Discussion
Our results show that rodent RA LTMRs innervating the glabrous skin respond heterogeneously to sustained and vibratory force (Fig. 1) and exhibit a continuum of frequency preferences (Fig. 3). Fitted GLMs reproduced experimental responses (Fig. 4), linking the continuum in frequency preference to a continuum in stimulus filter width (Fig. 5). The postspike filter did not correlate with frequency preference but helps explain (lack of) entrainment at high frequencies (Fig. 7). We also demonstrate that heterogeneity in tuning can improve the efficiency of population-level coding (Fig. 8).
Sensing of skin vibrations by LTMRs is vital for perception of object properties (e.g., texture; Hollins et al., 2002; Bensmaia and Hollins, 2003; Hollins and Bensmaia, 2007; O’Connor et al., 2021) and for sensorimotor control (e.g., slip detection; Talbot et al., 1968; Mountcastle et al., 1972; Knibestöl, 1973). In humans, monkeys, and cats, RA1s are most sensitive to ∼30 Hz, whereas RA2/PCs are most sensitive to ∼250 Hz (Talbot et al., 1968; Mountcastle et al., 1972; Iggo and Ogawa, 1977; Johansson et al., 1982). In contrast to this seemingly clear dichotomy, our data suggest that rodent RAs are more heterogeneous in their frequency preferences and sometimes have very broad tuning curves (Fig. 3). Previous studies in raccoons and cats also showed variability in RA tuning (Munger and Pubols, 1972; Gibson et al., 1975; Iggo and Ogawa, 1977), and this is also evident in recent mouse data (Neubarth et al., 2020). Conversely, RA tuning curves in monkey glabrous skin appear less heterogeneous (Mountcastle et al., 1967; Talbot et al., 1968) but are nonetheless quite wide, such that RA1s and RA2/PCs are coactivated by natural stimuli (Saal et al., 2015). Primate RAs are localized to skin areas necessary for gripping (Johansson and Flanagan, 2009), whereas RA2/PCs are rarely found in the rodent glabrous skin (Guzun et al., 2021) and are located instead in the joints and periosteum (Zelená, 1994; Prsa et al., 2019; Handler and Ginty, 2021). Notably, mouse fore and hind paws differ in LTMR innervation density and RA sensitivity (Walcher et al., 2018).
Recent work shows that mouse RA2/PCs are sensitive to kHz vibrations and that neurons in primary somatosensory cortex (S1) do not respond with precise spike timing at such high frequencies but, instead, that different frequencies are encoded by the activation of differently tuned neurons (Prsa et al., 2019), like in the auditory cortex (Wang et al., 2005; but see Downer et al., 2021). Precisely timed spikes are observed in rat (Ewert et al., 2008) and monkey (Harvey et al., 2013; Saal et al., 2015) S1 neurons for stimuli up to 700–800 Hz. Separate work suggests that rat LTMRs start to desynchronize at frequencies >350 Hz because of reduced entrainment (Sagalajev et al., 2024), and so encoding kHz vibrations may necessitate a different coding scheme than for lower vibration frequencies. Most textures evoke vibrations <400 Hz (Manfredi et al., 2014). We did not test >400 Hz for technical reasons, and none of our RAs had remarkably large receptive fields, suggesting we sampled RA1s exclusively. If so, then rodent RA1s can encode frequencies typically ascribed to RA2/PCs in other species (i.e., 100–400 Hz), leaving rodent RA2/PCs to encode even higher frequencies. This implies that rodents encode vibrations up to 400 Hz using heterogeneously tuned RA1s rather than dichotomously tuned RA1s and RA2/PCs. Independent control of the S1 spike rate and timing by RA1s and RA2/PCs, respectively, as demonstrated in monkeys (Saal et al., 2015), would not apply in rodents despite both rate and temporal codes being used in rodent S1 (Zuo et al., 2015; Lankarany et al., 2019; Kamaleddin et al., 2022).
Consistent with other species, rodent RAs responded transiently to sustained pressure (Fig. 1A) but had variable onset–offset responses. Heterogeneity in RA onset–offset responses was reported in humans and raccoons (Knibestöl, 1973; Pubols, 1980) but not in cats (Iggo and Ogawa, 1977). These studies agree that RAs favor the stimulus onset, and heterogeneity typically lies in the offset response (Knibestöl, 1973; Pubols, 1980). Most rodent RAs responded to both the stimulus onset and offset (Fig. 3E) but, based on differing thresholds, generally had a preference for one or the other (Fig. 4C). Recent work in mice identified two populations of RAs innervating Meissner corpuscles that respond at the stimulus onset and offset (TrkB+) or onset only (Ret+; Neubarth et al., 2020). Our models separated into two groups (Fig. 4C), but there was no consistent alignment between ON/OFF groups and the clusters in Figure 3. Cluster 1 neurons fire preferentially at the onset, but some can also fire at the offset (Fig. 4C), which is not consistent with the onset-only spiking exhibited by Ret+ neurons (Neubarth et al., 2020), and their slower CV is more consistent with TrkB+ neurons (Neubarth et al., 2020). Overall, our results suggest that physiologically defined RA “clusters” do not correspond with molecularly defined Ret+ or TrkB+ afferents and that functional heterogeneity is better conceptualized as a continuum.
To explore the basis for frequency preferences, we fit GLMs to RA responses (Fig. 4). We also simulated GLM responses to periodic data not used for fitting (Fig. 4B) and, although we did not use random input which would have provided a rich dataset for capturing neuronal responses (Chichilnisky, 2001; Pillow et al., 2008; Aljadeff et al., 2016), our models nonetheless generalized to nonperiodic input (i.e., sustained pressure; Fig. 4C), capturing quintessential characteristics of RA responses. Our GLMs revealed a continuum of stimulus filter widths (Fig. 5A) where narrower filters reflect a higher preferred frequency. Compared with RA/PC filters in primates (Saal et al., 2015), rodent RAs have strong similarities in filter width and excitation/suppression ratios (Fig. 5C), suggesting conserved intrinsic mechanisms for encoding across the frequency spectrum. A true coincidence detector has a filter with balanced positive and negative components, whereas an integrator has a broader filter with a positive bias (Ratté et al., 2013). This balance of excitation and suppression in RA stimulus filters was heterogeneous (Fig. 5Cii). Stimulus filter shape depends on the specialized end organs with which LTMRs associate (Handler and Ginty, 2021) plus the ion channels that convert receptor potentials into spike trains (Ratté et al., 2015); for example, Kv7.4 dysfunction increases RA sensitivity to low-frequency vibration (Heidenreich et al., 2011). In contrast, the postspike filter controls firing rate gain as a function of stimulus intensity (Gustafsson, 1984; Prescott and Sejnowski, 2008) and limits entrainment as stimulus frequency is increased (Fig. 7D). By controlling refractoriness, the postspike filter affects when a neuron can respond again to a repeating input, whereas the stimulus filter controls the preferred shape (kinetics) of that input. We did not observe any relationship between the stimulus and postspike filter shapes within a given neuron, but our results confirm that rodent RAs exhibit exquisite spike timing precision and can thus support temporal coding.
Broad overlapping tuning curves support combinatorial coding rather than labeled lines (Pouget et al., 2000; Butts and Goldman, 2006; Prescott and Ratté, 2012; Prescott et al., 2014; Saal and Bensmaia, 2014; Röth et al., 2021). Like trichromacy, where coactivation of broadly tuned cone photoreceptors supports color vision, touch appears to rely on the relative activation of broadly tuned LTMRs (which implies coactivation of differently tuned LTMRs) rather than relying on unique activation of particular, narrowly tuned LTMRs. Differently tuned LTMRs might fall into distinct groups (as seems to occur in primates) or they might fall along a continuum (as our data suggest for rodents). In either case, heterogeneously tuned populations can convey more information than homogeneously tuned ones (Padmanabhan and Urban, 2010) which appears to be true also for RAs (Fig. 8). One might expect sets of discretely tuned LTMRs based on association with distinct end organs, but heterogeneity in the end organs themselves, in other aspects of the mechanotransduction process (e.g., skin properties), or in LTMR excitability might all contribute to smearing the functional differences between transcriptionally defined afferent types (Usoskin et al., 2015). However, this does not mean molecular genetic differences are not important; for instance, receptive fields of Ret+ and TrkB+ RAs homotypically tile within subtypes, but heterotypically overlap across subtypes, which is beneficial for spatial acuity (Neubarth et al., 2020).
In summary, we found that rodent RAs respond heterogeneously to vibrotactile stimulation. Further analysis of their tuning curves suggests that RAs form a continuum of frequency preferences across the population, but more data are needed to strengthen this claim. Our modeling suggests that this frequency preference spectrum is supported by a continuum of stimulus filter widths. Rodent RAs also exhibit heterogeneity in their postspike filter, which is uncorrelated with variations in the stimulus filter. RAs can modulate their firing rate to support rate coding of stimulus intensity yet also produce precisely timed spikes as required for temporal coding of stimulus frequency. In that regard, individual rodent RAs are not so different from their primate counterparts. That said, rodent RAs are arguably more heterogeneous in their tuning, which might have implications for how information is decoded by downstream neurons.
Footnotes
This work was supported by Canadian Institutes of Health Research (CIHR) Foundation Grant 167276 (S.A.P.); CIHR Canada Graduate Scholar Doctoral (L.M., D.A.-B.) and Masters (A.H.) Awards; and Ontario Graduate Scholarship and SickKids Restracomp Awards (C.D.).
The authors declare no competing financial interests.
- Correspondence should be addressed to Steven A. Prescott at steve.prescott{at}sickkids.ca.
