pressure is uniform (the wave length of the vibrations is more than 30 cm at 1 kHz), the skin average displacement over the surface of the belt aperture was of the order of 0.3µm (see **Table 5**). For periodic stimuli, a two ways ANOVA (no interaction term) revealed a significant effect of NF (*p <* 0*.*001) but not of velocity on the vibrations magnitude. The larger the force, the more intense was the vibration. A linear regression revealed a significant slope of 0.13 of the logarithm of intensity as a function of the NF. For sandpapers, we found no effect of velocity.

The logarithm of the intensity is plotted on **Figure 7** as a function of the five stimuli's spatial period **(A–C)** and the three sandpaper's grit size **(B–D)**. The range of RMS values varied a lot across subjects. Before pooling the data of all subjects, the data of each individual were normalized with respect to a specific condition: *λ* = 0*.*4, NF = 2 N and *v* = 100 mm/s for periodic stimuli and P240, NF = 1 N and *v* = 100 mm/s for sandpapers. Normalization factors are shown in **Table 5** for each subject.

**Figure 7** presents data for all speeds and forces pooled together in black lines. The effect of NF (**Figures 7A,B**) and finger velocity (**Figures 7C,D**) is shown in colored lines. For periodic stimuli, the logarithm of the RMS magnitude increased linearly with the logarithm of the spatial period (slope = 0.24, *p <* 0*.*001, all forces pooled). The slope was significant for each level of force separately (*p* = 0.029, *p <* 0*.*001, and *p* = 0.003 for NF equal to 1, 2, and 4 N, respectively). There was a significant pair wise difference

corresponds to the theoretical value. The x-axis has a logarithmic scale. Logarithm bin width equals 0.015.

between each of P1-P2-P3 and P5 (*p <* 0*.*05, Tukey's test) and between P1 and P4 (*p <* 0*.*05, Tukey's test).

For sandpapers (**Figure 7B**), the grit size induced a significant increase of the vibration RMS (*p <* 0*.*05).

#### **COEFFICIENT OF FRICTION**

The coefficient of dynamic friction is plotted on **Figure 8** as a function of the five stimuli's spatial period **(A)** and the three sandpaper's grit size **(B)**. For the periodic stimuli the coefficient of friction was small (around 0.5) and constant across the different stimuli (*p* = 0.86).

The dynamic friction was much higher for the sandpapers, ranging from 0.9 for P800 to 1.3 for P80. We found a significant pair wise difference between P80 and P800 (*p* = 0.04, Tukey's test).

#### **DISCUSSION**

Libouton et al. (2012) suggested that the somatosensory system was able to collect roughness from other sources of information than fingertip afferences. The authors highlighted that the tactile roughness discrimination performance was unaffected by entrapment or traumatic section of the median nerve at the wrist. As a consequence, they concluded that, if the innervation of the finger pad was compromised, information about textures could be captured and encoded by remote mechanoreceptors located in more proximal tissues where the innervation was intact. We attempted to explain the mechanism through which humans can detect textured surfaces and discriminate between them when the innervation of the finger pad is radically compromised.

To this end, we recorded vibrations in the range from 30 Hz upto 1.2 kHz taking place in subjects' wrist tissues while they explored textured surfaces. Our results demonstrate the propagation of vibratory waves produced by the finger pad interaction

#### **Table 4 | Principal modes for the estimate of spatial period** *λ***ˆ for periodic and non-periodic stimuli.**

*The percentage of observations corresponds to the fraction of the observations that have a spatial period within* ±*20% of the spatial period of the stimuli.*

#### **Table 5 | Normalization factors for pressure variation and skin displacement for each subject (corresponding to 0 in the log scale of Figure 7).**

*The skin displacement was computed according to ∂V* = −*βV ∂p, with β being the air compressibility, V the belt volume and p the pressure. We considered the stethoscope belt as a cylinder, and the skin acting like a piston.*

with a textured surface to at least the wrist regions. We found that the vibration spectral energy was located on the fundamental frequency plus harmonics for periodic stimuli. This frequency, together with the fingertip velocity, corresponded to the spatial period of the stimulus. We did not find any obvious features that could reflect the grit size in the spectra for sandpapers. We also found that the coefficient of dynamic friction varied according to the sandpapers' grit sizes but remained constant across the different polyurethane gratings. These results corroborate the findings of Wiertlewski et al. (2011b) who measured the vibrations taking place at the finger-surface interface and found that the transformation of the surface geometry into vibrations was strongly non-linear, causing both frequency spread and background noise.

The present results show that the biomechanics of hand tissues can transmit surface interaction signals far away from the contact, up to the forearm, and probably much beyond. The

**FIGURE 8 | Comparison of coefficient of dynamic friction on the different explored surfaces. (A)** Periodic gratings: *λ* is the spatial period of the periodic stimulus. The x-axis has a logarithmic scale. **(B)** Sandpapers: from smooth (P800) on the left to rough (P80) on the right. Error bars indicate the standard error of the mean.

likely mechanical paths are bone conduction (Corso, 1963), tendon conduction (Pourcelot et al., 2005), surface Rayleigh waves (Liang and Boppart, 2010), shear waves in the bulk (von Gierke et al., 1951); or a combination thereof. Our results establish the possibility for the somatosensory system to combine information from multiple modalities and from distributed locations in order to assess properties of the characteristics of an unknown texture, using the fingertips as a distant probe and through a variety of different mechanisms yet to be described.

When rubbing a finger against sand paper, the spatial periodicity of the vibration signal in the wrist was in the range of the fingerprint periodicity, i.e., near 0.3 mm (see **Table 4**), suggesting that the fingerprints could make a spectral selection due to their intrinsic periodicity (Martinot, 2006; Prevost et al., 2009; Scheibert et al., 2009). For the periodic stimuli, some estimates of the spatial period were around 0.4 mm, especially for 0.8 mm gratings (see **Figure 4D**), which could also be linked to the fingerprint periodicity although a more likely explanation is simply a frequency doubling effect in a preferred range due to the presence of two edges per period of the stimulus.

It was hypothesized that the Pacinian channel could encode vibratory intensity (Hollins and Bensmaïa, 2007). This population of receptors is indeed a good candidate to encode scanned textures through remotely transmitted signals. Moreover, the exquisite sensitivity of the Pacinian corpuscule makes them good candidates for picking up signals having travelled over considerable distances (Hunt, 1961), see Hamann (1995) for a review. Significant vibrations in the wrist imply that, in healthy subjects, very large populations of receptors can be stimulated throughout the hand and the arm.

In the present study we have demonstrated that vibrations generated during the scanning of textured surfaces propagate through the finger and hand and stimulate receptor populations

error of the mean.

in regions far away from the contact region at least up to the wrist. The spectrum and the magnitude of the measured vibrations have to be considered as the best predictors of roughness discrimination. These findings have important implications in providing realistic sensory feedback for prosthetic-hand users. Indeed, it suggests that advanced prosthetic arms, equipped with sensory feedback, could partially restore the tactile sensation of amputees, through the activation of the remote mechanoreceptive system. Therefore, the transmission and reception of vibratory stimuli related to texture to the remaining stump becomes a major priority in the conception of future prosthetic hands.

#### **CONCLUSION**

We showed the presence of high frequency vibrations in the wrist of subjects exploring rough surfaces. The vibratory power is the best predictor of the texture roughness sensed remotely. It accounts for both periodic and non-periodic surfaces, with very different frictional properties. This result is consistent with perceptual studies carried on either with fingertip touch or with

#### **REFERENCES**

*Teleoperator Systems. World Haptics 2007*, (Tsukuba, Japan: IEEE), 66–71.

probe scanning. Moreover, we showed that the frequency content of periodic stimuli contained also potential information to encode roughness. Thus, many questions remain unclear about the contribution of the spectral content of the vibrations in the perceived roughness. Finally, friction may contribute to refine perception of roughness.

#### **ACKNOWLEDGMENTS**

The authors would like to thank the subjects for their kind participation, Simon Johnson and Bill Ranken from Unilever R&D Port Sunlight for providing the gratings, Vincent Théate for his technical support and Anne Klöcker for revision of the manuscript. This research was supported by a grant from Prodex, ESA (European Space Agency), PAI, FNRS, FRSM, ARC (BELGIUM), and a grant from the Commission of the European Union under Contract N◦ EU-FP7-NMP4-SL-2009-228844, NanoBioTouch. Additional funding was from the European Research Council, Advanced Grant PATCH, agreement No. 247300 awarded to Vincent Hayward.

coherence elastography. *IEEE Trans. Biomed. Eng.* 57, 953–959.

Yoshioka, T., Bensmaïa, S. J., Craig, J. C., and Hsiao, S. S. (2007). Texture perception through direct and indirect touch: an analysis of perceptual space for tactile textures in two modes of exploration. *Somatosens. Mot. Res.* 24, 53–70.

Yoshioka, T., Craig, J. C., Beck, G. C., and Hsiao, S. S. (2011). Perceptual constancy of texture roughness in the tactile system. *J. Neurosci.* 31, 17603–17611.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 28 March 2012; paper pending published: 16 April 2012; accepted: 19 June 2012; published online: 06 July 2012.*

*Citation: Delhaye B, Hayward V, Lefèvre P and Thonnard J-L (2012)* *Texture-induced vibrations in the forearm during tactile exploration. Front. Behav. Neurosci. 6:37. doi: 10.3389/ fnbeh.2012.00037*

*Copyright © 2012 Delhaye, Hayward, Lefèvre and Thonnard. This is an openaccess article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

## Rasch-built measure of pleasant touch through active fingertip explorations

#### **Anne Klöcker 1 , Carlyne Arnould2 , Massimo Penta1 and Jean-Louis Thonnard1\***

1 Institute of Neuroscience, Université Catholique de Louvain, Brussels, Belgium

2 Department of Occupational Therapy, Paramedical Category, Haute Ecole Louvain en Hainaut, Montignies-sur-Sambre, Belgium

#### **Edited by:**

Pavel M. Itskov, Champalimaud Foundation, Portugal

#### **Reviewed by:**

Helena Backlund Wasling, Gothenburg University, Sweden Ekaterina Vinnik, Champalimaud Foundation, Portugal

#### **\*Correspondence:**

Jean-Louis Thonnard, Institute of Neuroscience, Université Catholique de Louvain, Tour Pasteur, 53 Avenue Mounier (B1.53.04), 1200 Brussels, Belgium. e-mail: jean-louis.thonnard@ uclouvain.be

**Background:** Evidence suggests that somatic sensation has a modality for pleasant touch. **Objective:** To investigate pleasant touch at the fingertip level (i.e., glabrous skin site) through the elaboration of a linear unidimensional scale that measures (i) various materials according to the level of pleasantness they elicit through active fingertip explorations and (ii) subjects according to their pleasantness leniency levels. **Subjects:**We enrolled 198 healthy subjects without any neurological disease. **Methods:** Blindfolded subjects actively explored 48 materials with their index fingertips and reported the perceived pleasantness of each on a 4-level scale. The fingertip moisture levels on each subject were measured before the experimental session. Data were analyzed using the Rasch model. **Results:** We elaborated unidimensional linear scale that included 37 materials according to their pleasantness of touch. The pleasantness level of 21 materials was perceived differently, depending on the fingertip moisture levels of the subjects. **Conclusion:** Based on our findings, we formulated a Pleasant Touch Scale. Fingertip moisture levels appeared to be a major factor for (un)pleasant feelings during active exploration.

**Keywords: active touch, pleasantness, Rasch model, latent variable measurement, friction, moisture**

#### **INTRODUCTION**

The hand is the key organ of the sense of touch (Lundborg, 2005), making sensation fundamental to hand function. Four major modalities of somatic sensation have been described: (i) discriminative touch, (ii) proprioception, (iii) nociception, and (iv) temperature sense (Lundborg, 2005), although studies over the last decade have suggested the existence of a supplementary modality, namely *pleasant touch* (Johansson et al., 1988; Nordin, 1990; McGlone et al., 2007; McGlone and Reilly, 2010).

Studies investigating pleasant touch at hairy skin sites have suggested that C-tactile nerve fibers (CT-fibers) play a fundamental role in the detection and transmission of pleasant stimuli applied to the skin (Olausson et al., 2002, 2010;McGlone et al., 2007; Essick et al., 2010; Morrison et al., 2011). These CT-fibers were identified during microneurography recordings from infra- and supraorbital nerves (Johansson et al., 1988; Nordin, 1990). CT-fibers were later observed in other hairy skin areas, such as the arm (Vallbo et al., 1999). However, they seem to be missing from glabrous skin sites, such as the fingertips (McGlone and Reilly, 2010; Olausson et al., 2010). Studies investigating pleasant touch at glabrous skin sites frequently used either magnitude estimation (ME) or categorical rating procedures. The results of such studies suggested that soft and smooth materials were perceived as more pleasant than rough materials (Major, 1895; Ripin and Lazarsfeld, 1937; Ekman et al., 1965; Verrillo et al., 1999).

Pleasant touch may be described as a positive sensation induced by cutaneous stimulation. Hence, pleasant touch cannot be directly observed or measured and may be regarded as a *latent variable* rather than an observable variable. Observable variables can be directly quantified and typically generate linear measures (e.g., the grain size of sandpapers can be measured in micrometers). Latent variables can only be measured indirectly (e.g., pain, intelligence, or pleasant touch), generally by using a questionnaire or a set of stimuli (Thurstone, 1928; Rasch, 1960; Tesio et al., 2007). Questionnaires usually provide ordinal scores, rather than linear scores, with no constant unit. For example, pleasant touch could be assessed by presenting various stimuli to a subject and asking the subject to rate these stimuli on the basis of a 4-level response scale, e.g., "very unpleasant" (scored 0), "unpleasant" (scored 1), "pleasant" (scored 2), and "very pleasant" (scored 3). Two major issues would, however, prevent the widespread uniform use of numbers resulting from such a response scale. First, the scores would be ordinal, meaning that equal differences (e.g., from scores of 0–1 and 2–3) would not necessarily represent the same change/quantity of pleasant touch. Second, the scores attributed to the stimuli would not necessarily be unidimensional; thus, adding these individual scores may result in a total score that conceals numerous unrelated dimensions (e.g., color, hardness, roughness, temperature, mental representation). As a consequence, quantitative comparisons of these ordinal scores are not possible (Merbitz et al., 1989; Wright and Linacre, 1989). Nevertheless, probabilistic measurement models can be used to determine linear and unidimensional measures from ordinal scores, the most promising being the Rasch model (Rasch, 1960). The Rasch model is a prescriptive model, rather than a descriptive one, requiring that solely the location of the item (e.g., the amount of pleasantness elicited by a stimulus) and the ability of a subject (e.g., the "satisfaction level" induced by perception of a stimulus) determine the probabilities of category choices on the response scale. When applied to pleasant touch, this measurement framework can be

used to construct a linear and unidimensional *Pleasant Touch Scale* assessing (i) each stimulus according to its level of pleasantness and (ii) each subject according to his/her leniency level, defined as a subject's "satisfaction level" induced by perception of the stimulus.

The aim of this study was to investigate pleasant touch at a glabrous skin site through the elaboration of a gradual hierarchy of stimuli eliciting pleasant touch (i.e., the *Pleasant Touch Scale*). This unidimensional linear scale was designed (i) to quantify the pleasantness level elicited by various materials during direct active fingertip explorations and to grade these materials by pleasantness level and (ii) to quantify the pleasantness leniency level of subjects and to grade them by pleasantness leniency level. The elaboration of such a *Pleasant Touch Scale* may contribute to knowledge of pleasant touch at glabrous skin sites. Indeed, the *quantitative* information obtained by this scale may provide objective insight of the characteristics (of subjects and/or materials) that influence pleasant touch. Consequently, future research on pleasant touch may refer to such scale and specifically examine the characteristics that influence pleasantness or unpleasantness.

#### **MATERIALS AND METHODS**

The present study was approved by the Biomedical Ethical Commission of the Faculty of Medicine of the Université catholique de Louvain.

#### **SUBJECTS**

We enrolled 198 healthy subjects, 74 males and 124 females, aged between 20 and 70 (mean age, 39.5 ± 13.4 years).

#### **INSTRUMENT**

To elaborate the *Pleasant Touch Scale*, we selected 48 materials, encountered in everyday life, based on three perceptual tactile dimensions (i.e., hard vs. soft, rough vs. smooth, and sticky vs. slippery; Hollins et al., 1993, 2000; Yoshioka et al., 2007). Of these 48 stimuli, 26 had been used in studies dealing with the discriminative aspect of touch (Major, 1895; Hollins et al., 1993, 2000; Bergmann Tiest and Kappers, 2006;Yoshioka et al., 2007;Yoshioka and Zhou, 2009). The remaining 22 materials included everyday life materials such as baking paper, cast, cling film, and tights. Each selected material was glued onto an aluminum plate, 77 mm long and 32 mm wide.

#### **PROCEDURE**

Before each experiment, the index fingertip moisture level of each subject was measured using the Corneometer®CM 825. The subjects were blindfolded, and each of the 48 materials was placed in random order in front of the subjects. The subjects were instructed to place the pulp of their dominant index finger on the selected material and to explore its surface through slow lateral sliding movements. No specific instructions on the normal force to apply to each material or on the speed of sliding were provided. Exploration of each material could be repeated as many times as required, after which each subject was asked to rate the material's pleasantness level on a 4-level scale: (0) very pleasant, (1) pleasant, (2) unpleasant, or (3) very unpleasant (**Figure 1**). Subjects hesitating between categories (1) pleasant and (2) unpleasant were

instructed to choose "unpleasant" as soon as the material was not perceived as being pleasant or "pleasant" as soon as the material was not perceived as unpleasant. Each subject's total score could range from 0 (i.e., all materials rated as very pleasant, 48 × 0) to 144 (i.e., all materials rated as very unpleasant, 48 × 3). Similarly, the total score of each material could range from 0 (i.e., all subjects rated it as very pleasant, 198 × 0) to 594 (i.e., all subjects rated it as very unpleasant, 198 × 3).

#### **DATA ANALYSIS**

Data were analyzed using the Rasch Unidimensional Measurement Models software (RUMM2020). The rating scale model was preferred to the partial credit model, since it yields better fit of the data to the model (Wright, 1999).

The Rasch model is a *probabilistic* model elaborated by the Danish mathematician Rasch (1960). The Rasch model allows the elaboration of unidimensional and linear scales for latent variables (e.g., pain, intelligence, or pleasant touch) by locating subjects and items along *a single* underlying linear scale (Rasch, 1960; Wright and Stone, 1979; Tennant et al., 2004; Tennant and Conaghan, 2007). This model was originally elaborated to analyze dichotomous data (i.e., two response categories per item such as "yes/no" or "pass/fail"). The model formulates the probability of success of a subject with a given ability on an item with a given difficulty. No other parameter besides the subject's ability (i.e., his/her location) and the item's difficulty (i.e., its location) is needed to determine the expected response. In this sense the model is unidimensional as it involves only one latent dimension of the subjects.

Later on, the model was expanded to accommodate polytomous response formats (i.e., more than two response categories per item such as"very low/low/moderate/high/very high";Wright and Masters, 1982; Andrich, 1988). For polytomous response formats, the model prescribes that the probability of endorsing any response to an item depends solely on the subject's ability, the item's difficulty, and the threshold difficulties, where thresholds (i.e., "boundaries" between successive categories) are the locations on the underlying scale at which successive responses become the most probable if subjects of increasing ability answer a given item (Andrich, 1978a,b; Wright and Masters, 1982; Tennant et al., 2004; Tennant and Conaghan, 2007).

When applied to touch related pleasantness perception, the Rasch model prescribes that the probability that a subject selects a given response category (i.e., pleasantness category) for a given item (i.e., material) depends solely on (i) the subject's pleasantness leniency level, (ii) the material's pleasantness level, and (iii) the threshold locations. Based on the estimated subject leniency and the material pleasantness, the expected response of a subject to an item can be computed by the model. The similarity between the observed and expected responses to any item is reported by the software, through a χ 2 fit statistic (Andrich et al., 2004). Materials with a significant χ 2 value indicate that the difference between the observed and expected responses was too high to be random; they do not fit the unidimensional scale defined by the other materials and can be eliminated. In addition, the hypothesized order of the response categories (i.e., "very pleasant," "pleasant," "unpleasant," and "very unpleasant") must be empirically verified. Items presenting disordered thresholds between successive response categories indicate that the subjects were unable to discriminate the four response categories, allowing these items to be eliminated (Andrich, 1996). The resulting scale is defined in logits (log-odds-units), a unit defined as the natural logarithm of the odds of success (i.e., the pass/fail probability ratio) of a subject to an item and which is constant throughout the measurement scale. Consequently, at any level of the measurement scale, a 1-logit difference in an item's pleasantness level indicates a constant ratio of its odds of being pleasant (e1 = 2.71) for any subject, whereas a 2-logit difference indicates the odds of being pleasant at a ratio of e2 = 7.39, etc. The origin of the scale is conventionally set at the average material pleasantness level.

#### **MATERIAL SELECTION**

Based on the 48 original materials, successive analyses were performed to select materials (i) presenting ordered response categories and (ii) fitting a unidimensional scale.

#### **Ordered response scale**

The subjects were asked to report their pleasantness perception on a 4-level scale: (0) very pleasant, (1) pleasant, (2) unpleasant, or (3) very unpleasant. To control the order of response categories, whether successive categories represent decreasing pleasantness levels was verified for each material. If the response categories were ordered, the thresholds were located in the anticipated order. Provided that the anticipated order of response categories was verified, less lenient subjects should have selected a higher response

for any given material,whereas subjects selecting a higher response to a given material should be less lenient. However, a reversed order of the thresholds between successive response categories indicated that the response scale was not used as expected (Tennant, 2004). If the categories were not distinguished by the sample, the concerned categories were collapsed into a single category (Tennant and Conaghan, 2007).

#### **Unidimensional scale**

Different statistic methods are available to test the fit of the data to the model requirements. This study used item χ 2 fit statistics. In this statistical method, the squared standardized residuals (i.e., the difference between the response expected by the model and the observed response by the subject) of all subjects are summed, leading to a χ 2 value for each material. Subsequently, a significance test was used to evaluate whether the χ 2 was too high to be attributed to random variation. A *p*-value of the item fit statistics of <0.05 indicated a misfitting item that may threaten the concept of unidimensional (Tennant and Conaghan, 2007). There are two types of misfits, *underfit* and *overfit*. Underfitted items violate the unidimensionality concept as the subjects' answers to these items are influenced by variables other than the measured one (Penta et al., 2005). An overfitted item does not necessarily represent a threat to the unidmensionality concept (Penta et al., 2005). Indeed, the response patterns of overfitted items are more deterministic than those predicted by the Rasch model. However, in contrast to underfitted items, the pleasantness levels of overfitted items can be predicted from each subject's pleasantness leniency level (Penta et al., 2005). As a consequence, only underfitted items (i.e., materials) were eliminated from this study.

#### **DIFFERENTIAL ITEM FUNCTIONING**

Once satisfactory metric properties were achieved for the whole sample, the invariance of the pleasantness hierarchy of items within the sample was tested through a differential item functioning (DIF) analysis. DIF occurs if subjects of distinct subgroups (e.g., males vs. females) with the same leniency level perceive any given material differently (Tennant et al., 2004). In this study, DIF was investigated according to (i) gender, (ii) age (≤37 vs. >37 years, the median age of the subjects), and (iii) fingertip moisture level (≤70 vs. >70 arbitrary units, or "a.u.," the median moisture level). To investigate DIF, each subgroup was divided into five class intervals (CI) of decreasing pleasantness leniency levels, and two-way analysis of variance (ANOVA) was computed on the standardized residuals of the different CIs (Andrich et al., 2004). Factors analyzed in two-way ANOVA included (i) subject subgroups (e.g., males vs. females) and (ii) CIs of decreasing pleasantness leniency. A significant subgroup's main effect indicated a uniform DIF, occurring if the standardized residuals change significantly from one subgroup to another and evolve in parallel. In such a case, the relative difficulty of the item would differ between these subgroups (Smith, 1992). An item displaying a uniform DIF may be deleted or split into as many specific items as there are subject subgroups; each specific item would therefore have a difficulty peculiar to the corresponding subject's subgroup (Tennant et al., 2004).

#### **RELIABILITY OF THE PLEASANT TOUCH SCALE**

To assess the reliability of the *Pleasant Touch Scale*, a separation reliability index was computed for the items and subjects (Wright and Masters, 1982). The Item Separation Index enables the number of material pleasantness levels that may statistically be distinguished within the sample to be calculated. The Person Separation Index enables the number of statistically distinguished subject pleasantness leniency levels using the selected materials to be calculated. A higher item/person separation index indicates better separation of the measures (Wright and Masters, 1982).

#### **RESULTS**

Successive Rasch analyses were performed to construct the final *Pleasant Touch Scale*. Of the 48 original materials (*hereafter* items), 11 showed disordered thresholds, indicating that our 198 subjects were unable to discriminate among the four levels of pleasantness for these items. As the subjects were unable to discriminate between "unpleasant" and "very unpleasant," these categories were merged into one category, "unpleasant," and the entire data set was reanalyzed using a 3-level scale: (0) very pleasant, (1) pleasant, and (2) unpleasant. Two items, however, adhesive UHU® Patafix and microfiber dust cloth, still presented disordered thresholds and were deleted from further consideration. Rasch analysis showed that six items (i.e., tissue carpet, thermal isolation, sponge, antislide carpet, velcro, and plastic door mat) did not fit a unidimensional scale and were therefore also eliminated, which resulted in a 40-item scale. No significant DIF based on age and gender was observed for any item, but invariance investigation of the items according to the subjects' levels of fingertip moisture showed that 17 items shared a common location for low and high moisture level subgroups, whereas the other 23 items presented a DIF. These 23 items elicited different levels of pleasantness when touched by subjects with dry (i.e., ≤70 a.u.) and wet (i.e., >70 a.u.) skin. Each of these 23 items was therefore split into two different items with locations specific to fingertip moisture level, one each for subjects with dry and wet skin. After the splitting procedure, four split items and one common item showed underfitting (i.e., teflon high moisture, teflon low moisture, crepe paper high moisture, crepe paper low moisture, and corduroy), and were therefore removed. Consequently, the final *Pleasant Touch Scale* includes 58 items, of which (i) 16 items share a common location in low and high moisture level subgroups and (ii) 42 items have locations specific to the fingertip moisture level.

#### **METRIC PROPERTIES OF THE PLEASANT TOUCH SCALE**

The calibration of the 58 items of the *Pleasant Touch Scale* is presented in **Table 1**, in which the items were ordered from the most unpleasant at the top ("sandpaper\_LM") to the most pleasant at the bottom ("paper\_160 g/m2 "). Items followed by "\_LM" and "\_HM" are split items with pleasantness locations specific to subjects with low and high fingertip moisture levels, respectively. The pleasantness levels of the 58 materials covered a range of 6.91 logits (range: −4.47 to 2.44 logits), indicating that the odds of pleasing any particular subject was in a ratio higher than 1000:1 (i.e., e6.91:1) between the most and least pleasant items. This pleasantness range was arbitrarily centered at 0 logits. **Table 1** also shows the pleasantness level of each material expressed in percent, with 0 and 100% indicating the least and most pleasant materials, respectively. Furthermore, the table presents the standard errors (SE) associated with the estimations of the pleasant levels of the different items (mean: 0.18 logits; range: 0.36 logits). The item fit statistics (χ 2 with associated *p*-values) indicate that all 58 items contribute to the definition of a unidimensional scale. Finally,**Table 1** also shows the perceptual tactile dimension of each material.

#### **DESCRIPTION OF THE PLEASANT TOUCH SCALE**

**Figure 2** illustrates the structure of the final *Pleasant Touch Scale*. The lower part shows the relationship between the total ordinal raw score and the subjects' linear pleasantness leniency levels. As the data were analyzed using a 3-level scale, the subjects' ordinal raw scores ranged from 0 (i.e., subjects finding all 58 materials very pleasant, 0 × 58) to 116 (i.e., subjects finding all 58 materials unpleasant, 2 × 58). The relationship between total score and pleasantness leniency levels was roughly linear for total scores between 25 and 90. In this area, an increase of 1 point in total raw score represents always the same progression in pleasantness measures. For example, increases from 40 to 41 and from 60 to 61 correspond to increases in pleasantness leniency of 0.05 logits. Outside this central range, a 1 point increase in total ordinal raw score did not correspond to the same progression in linear pleasantness measures. Indeed, an increase in total raw score from 114 to 115 corresponded to 0.55 logits, whereas an increase from 5 to 6 corresponded to 0.28 logits. This large difference in the progression of the pleasantness leniency levels for a same increase in total ordinal raw score illustrates the non-linearity of the raw scores (Wright and Stone, 1979; Arnould et al., 2004).

The upper part of **Figure 2** shows that the distribution of pleasantness leniency levels in subjects ranged from −2.62 to 3.42 logits, indicating that their odds of leniency levels are in a ratio over 400:1 (i.e., e6.04:1) between the most and the least lenient subject. Subjects on the left represent the most lenient (i.e., "the most easy to please") and those on the right represent the least lenient (i.e.,"the least easy to please"). The average subject location was 0.68 logits (SD: 0.88 logits), indicating that overall, the subjects found the materials more unpleasant than pleasant. In addition, the distribution shows that the pleasantness leniency of subjects with wet fingers (blank bars in **Figure 2**) and dry fingers (black bars in **Figure 2**) alternate on the *Pleasant Touch Scale*. This indicates a lack of relationship between fingertip moisture and pleasantness leniency [*t*-test, *t* = −0.849; *p* (two-tailed) = 0.397].

The middle part of **Figure 2** shows items of the *Pleasant Touch Scale*, representing the most unpleasant at the top to the most pleasant at the bottom. Each line presents the most probable response of the subjects to the different items of the *Pleasant Touch Scale* as a function of their pleasantness leniency levels. For example, to find the material "tights" very pleasant, a subject would have to have a leniency level of ≤−0.378 logits. A subject with a leniency level ≥2.516 logits would be expected to find this same material unpleasant. Subjects with leniency levels between −0.378 and 2.516 logits would be expected to find "tights" pleasant.

The expected responses of the subjects to each item are obtained by comparing the pleasantness leniency levels of the subjects to the pleasantness level of each item. For example, a subject with a total raw score of 58, corresponding to a pleasantness

#### **Table 1 | Description of materials of the pleasant touch scale.**

(Continued)

#### **Table 1 | Continued**

Items are ordered from least to most pleasant.

SE, standard error; df, degree of freedom; LM, low fingertip moisture; HM, high fingertip moisture.

leniency level of 0 logit, has the highest probability of finding (i) items 1 ("sandpaper\_LM") through 9 ("rough sponge\_HM") and 12 ("leather chamois") unpleasant, (ii) items 10 ("carbon paper\_HM") through 11 ("linen"), 13 ("tesa tape\_HM") through 46 ("paper 250 g/m2\_LM"), and 48 ("cotton tissue\_HM") through 53 ("synthetic tissue") pleasant, and (iii)items 47 ("marble\_LM") and 54 ("velvet") through 58 ("paper\_160 g/m2 ") very pleasant (see middle part of **Figure 2**).

#### **RELIABILITY OF THE PLEASANT TOUCH SCALE**

The material separation coefficient was equal to 0.98, indicating that 10.14 material pleasantness levels could be statistically distinguished by the sample. The reliability subject separation coefficient was equal to 0.88, indicating that 3.97 subject pleasantness leniency levels could be statistically distinguished using the selected materials (Wright and Masters, 1982).

#### **DISCUSSION**

We have described here the construction of a linear, unidimensional *Pleasant Touch Scale*, which quantifies (i) the level of pleasantness elicited in subjects by various materials during direct active fingertip explorations and (ii) the pleasantness leniency level of subjects. To construct this scale, 198 healthy subjects successively explored 48 different materials and reported their perception of pleasantness using a 4-level pleasantness scale. The collected data were analyzed using Rasch model. The initial analysis indicated that subjects did not discriminate between "unpleasant" and "very unpleasant."Consequently, these two categories were merged, and the data were analyzed using three response levels. In addition, several materials had to be eliminated because they did not fit the unidimensional scale defined by the other items. This indicates that the pleasantness levels of these materials were influenced by dimensions other than pleasantness, suggesting that subjects associated the sensation of these materials with, for example, positive or negative events in their lives. The investigation of the invariance of the scale showed that the levels of pleasantness of more than half of the items differed according to the moisture levels on the subjects' fingertips. Indeed, each of these DIF items was perceived differently, depending on whether subjects had dry or wet skin, and could therefore be split into two items with different locations specific to fingertip moisture level. The final *Pleasant Touch*

*Scale* orders 58 items – 16 with a common location in subjects with low and high moisture levels and 42 with locations specific to the fingertip moisture level – arranged according to their levels of pleasantness and 198 subjects arranged according to their levels of pleasantness leniency.

In this study, the data were collected using a category rating (CR) scale. Indeed, each stimulus was initially rated on a 3-level scale, yielding ordinal scores. An alternative data collection procedure would be ME. In contrast to CR, ME is an unlimited rating scale method, allowing subjects to freely choose a number reflecting their perception of a (set of) stimulus (stimuli), without presetting rating categories (Stevens, 1975). ME was initially designed, i.a. to yield ratio level scales, although this remains unclear (Wills and Moore, 1994). Indeed, both ME and CR have been found to yield ordinal scores (Wills and Moore, 1994), with both methods generating data/scores lacking essential psychometric properties (i.e., linearity and unidimensionality), precluding the objective and quantitative comparison of the measured variable. Therefore, in order to allow objective comparisons of the measured variable (e.g., pleasantness leniency level), the data/scores obtained through both scaling methods must be transformed into linear and unidimensional measures. Furthermore, compared with CR scale procedures, it has been suggested that the ME methods increase response variability, lowering the power of a test (Mellers, 1983). We therefore collected data using a CR scale, followed by transformation into linear unidimensional measures through Rasch model analysis.

The hierarchy of items sharing the same pleasantness level regardless of the fingertip moisture levels of the subjects indicates that rough and sticky items are more unpleasant to explore than smooth, slippery, or soft items. Indeed, the most unpleasant materials were "latex," "wax," "cling film," and "linen," while the most pleasant ones are "silk," "viscose tissue," "foam," "tights," "paper (70 g/m2 )," "baking paper," "synthetic tissue," "velvet," "cellular rubber," "paper (80 g/m2 )," and "paper (1690 g/m2 )." Similarly, other studies have found that soft and smooth materials were rated as most pleasant to touch, whereas rough materials were perceived as unpleasant (Major, 1895; Ekman et al., 1965; Essick et al., 1999, 2010; Verrillo et al., 1999).

Subjects with high fingertip moisture levels perceived *rough and sticky materials* as *more pleasant* than subjects with low fingertip

**FIGURE 2 | (A)** Subject distribution according to pleasantness leniency levels (expressed in logits). The blank bars correspond to subjects with high fingertip moisture and the black bars to subjects with "low fingertip moisture. **(B)** Classification of materials by pleasantness, increasing from top to bottom.

Items followed by "\_LM" and "\_HM" are split items with pleasantness locations specific to subjects with low and high fingertip moisture levels, respectively. **(C)** Relationship between the total ordinal raw score and linear pleasantness measures.

moisture levels. Increased fingertip moisture was found to increase the dynamic coefficient of friction (i.e., the ratio between friction and normal force during slippage; André et al., 2011) until a certain moisture threshold (Nacht et al., 1981; Tomlinson et al., 2007; Gerhardt et al., 2008). Similarly, skin hydration can increase finger's adhesion to a contacting surface, decreasing the probability of occurrence of slippage (André et al., 2011). Above the before mentioned moisture threshold, however, a further increase of skin hydration decreases the coefficient of friction (Tomlinson et al., 2007), due to the accumulation of a thin layer of moisture between the explored surface and the skin. Consequently, it may be hypothesized that rough and sticky materials, perceived as more pleasant if fingertip moisture was high, were materials that allowed fingertip moisture to accumulate to a very high level on their surfaces, resulting in the formation of a thin layer of moisture between the explored material and the skin of the fingertip. This may decrease the dynamic coefficient of friction of sticky materials and the roughness perception of rough materials, increasing the perceived pleasantness. This last assumption is strengthened by studies highlighting that the rough nature of rough materials are less perceived with higher skin hydration (Verrillo et al., 1998, 1999).

Conversely, subjects with high fingertip moisture levels perceived *smooth materials* as *less pleasant* than subjects with low fingertip moisture levels. It may be hypothesized that these smooth materials assimilated lower levels of moisture on their surfaces and failed to reach the moisture threshold above which the coefficient of friction decreases. The accumulated moisture may increase the

#### **REFERENCES**

dynamic coefficient of friction,leading to a less pleasant perception of the explored materials, which are perceived as stickier.

Finally, hardness was found to be unrelated to the perception of pleasantness, with some (e.g., tile\_LM, cast\_HM, marble\_LM) perceived as pleasant and others (e.g., wood\_LM, argil\_LM, plexiglass\_HM) as unpleasant. Rather, pleasantness was more dependent on roughness and stickiness than on hardness.

The aim of this study was to investigate pleasant touch of glabrous skin. Results suggest that a subject's pleasantness leniency level, elicited by active touch, may be defined as a latent trait common among subjects. The *Pleasant Touch Scale* is a unidimensional and linear scale, ordering materials according to their level of pleasantness and healthy subjects according to their level of pleasantness leniency. Consequently, our findings provide *objective* insights into the characteristics of materials and subjects that influence pleasant touch. Indeed, we found that (i) smooth and soft materials are perceived as more pleasant than rough and sticky materials and (ii) the level of fingertip moisture influences the perception of pleasantness. Furthermore, this scale may be useful in future studies, aimed, for example, at determining correlations between the pleasantness levels of materials and the dynamics of active touch.

#### **ACKNOWLEDGMENTS**

We are grateful for the kind participation of all subjects. This work was supported by a grant from the Commission of the European Union under Contract No. EU-FP7-NMP4-SL-2009-228844, NanoBioTouch.

multidimensional scaling. *Percept. Psychophys.* 62, 1534–1544.

foundations of misinference. *Arch. Phys. Med. Rehabil.* 70, 308–312.

tactile sensation. *Somatosens. Mot. Res.* 15, 93–108.

modes of exploration. *Somatosens. Mot. Res.* 24, 53–70.

Yoshioka, T., and Zhou, J. (2009). Factors involved in tactile texture perception through probes. *Adv. Robot.* 23, 747–766.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 13 March 2012; paper pending published: 30 April 2012; accepted: 30 May 2012; published online: 21 June 2012.*

*Citation: Klöcker A, Arnould C, Penta M and Thonnard J-L (2012) Rasch-built measure of pleasant touch through active fingertip explorations. Front. Neurorobot. 6:5. doi: 10.3389/fnbot.2012.00005*

*Copyright © 2012 Klöcker, Arnould, Penta and Thonnard. This is an openaccess article distributed under the terms of the Creative Commons Attribution Non Commercial License, which permits non-commercial use, distribution, and reproduction in other forums, provided the original authors and source are credited.*

## Whisker encoding of mechanical events during active tactile exploration

#### *Yves Boubenec 1, Daniel E. Shulz 1 and Georges Debrégeas 2 \**

*1 Unité de Neurosciences Information et Complexité, UPR 3293, Centre National de la Recherche Scientifique, Gif-sur-Yvette, France 2 FRE3231 Laboratoire Jean Perrin, Centre National de la Recherche Scientifique, UPMC Université Paris 6, Paris, France*

#### *Edited by:*

*Blythe Towal, California Institute of Technology, USA*

#### *Reviewed by:*

*Mitra J. Hartmann, Northwestern University, USA Garrett B. Stanley, Georgia Institute of Technology, USA Jason Ritt, Boston University, USA*

#### *\*Correspondence:*

*Georges Debrégeas, FRE3231 Laboratoire Jean Perrin, Centre National de la Recherche Scientifique, UPMC Université Paris 6, 24 rue Lhomond, 75005 Paris, France. e-mail<EMAIL_ADDRESS>Rats use their whiskers to extract a wealth of information about their immediate environment, such as the shape, position or texture of an object. The information is conveyed to mechanoreceptors located within the whisker follicle in the form of a sequence of whisker deflections induced by the whisker/object contact interaction. How the whiskers filter and shape the mechanical information and effectively participate in the coding of tactile features remains an open question to date. In the present article, a biomechanical model was developed that provides predictions of the whisker dynamics during active tactile exploration, amenable to quantitative experimental comparison. This model is based on a decomposition of the whisker profile into a slow, quasi-static sequence and rapid resonant small-scale vibrations. It was applied to the typical situation of a rat actively whisking across a solid object. Having derived the quasi-static sequence of whisker deformation, the resonant properties of the whisker were analyzed, taking into account the boundary conditions imposed by the whisker/surface contact. We then focused on two elementary mechanical events that are expected to trigger significant neural responses, namely (1) the whisker/object first contact and (2) the whisker detachment from the object. Both events were found to trigger a deflection wave propagating upward to the mystacial pad at constant velocity of ≈3–5 m/s. This yielded a characteristic mechanical signature at the whisker base, in the form of a large peak of negative curvature occurring ≈4 ms after the event has been triggered. The dependence in amplitude and lag of this mechanical signal with the main contextual parameters (such as radial or angular distance) was investigated. The model was validated experimentally by comparing its predictions to high-speed video recordings of shock-induced whisker deflections performed on anesthetized rats. The consequences of these results on possible tactile encoding schemes are briefly discussed.

**Keywords: exploration, rat, resonance, tactile, vibration, vibrissae, whiskers, whisking**

#### **INTRODUCTION**

The vibrissal system of the rat is one of the prominent model systems for investigating the mechanisms of sensory information processing in the tactile modality. Rats use their whiskers to sense their close environment and gather information about object features such as location (Krupa and Nicolelis, 2001; O'Connor et al., 2010a), shape (Brecht and Merzenich, 1997; Polley et al., 2005), texture (Carvell and Simons, 1990; Morita et al., 2011), and size (Anjum et al., 2006). Active movements of body, head (Milani et al., 1989; Carvell and Simons, 1990; Towal and Hartmann, 2006; Mitchinson et al., 2007) or whiskers themselves during whisking (Welker, 1964; Carvell and Simons, 1990; Berg and Kleinfeld, 2003) induce contact between the whisker and the probed environment. This mechanical interaction elicits sequences of whisker deflection. Each whisker is embedded in a follicle in the skin (Ebara et al., 2002), where mechanoreceptors transduce whisker base deflections into neural signals (Lichtenstein et al., 1990; Szwed et al., 2006). Neurons along the

trigeminal pathway respond to various aspects of the whisker base movements, such as high acceleration events (Jadhav et al., 2009; Lottem and Azouz, 2009; Jadhav and Feldman, 2010), whisker speed (Arabzadeh et al., 2004), average noise level (Arabzadeh et al., 2005), or characteristic features of the whisker motion spectra (Hipp et al., 2006). These different properties of the whisker base dynamics are then processed by the central nervous systems to extract relevant features of the environment. In order to decipher the underlying neural code, one needs to relate the dynamics of the whisker to the physical and geometrical characteristics of the contacting surface. This amounts to understanding the way the whisker carries and shapes information from the contact point to the follicle.

Consider a typical exploration task during which a rat whisks across an object (**Figure 1**). This sequence can be broken up into three consecutive phases. First, the whisker rotates freely in air; second, it slides over the object and is submitted to a frictional contact at the whisker tip; the whisker then detaches and pursues

its motion in air (no contact). Recently, several groups have studied how the rat could extract the contour of an object from the time evolution of the torque at the base of each of its whiskers (Kaneko et al., 1998; Scholz and Rahn, 2004; Clements and Rahn, 2006; Kim and Möller, 2007; Solomon and Hartmann, 2010). The encoding mechanism requires proprioceptive information (the angular position of the whisker with respect to the snout) and the knowledge of the relationship between the radial distance and the resulting torque in the follicle. This relationship was derived by computing the successive equilibrium profiles of the whiskers for a non-frictional contact. Due to the slender geometry of the whisker, any change in the contact configuration (during the first contact or following the whisker detachment from the object) is expected to trigger a burst of whisker oscillations ignored in these quasi-static descriptions. During the sliding phase itself, stickslip instabilities are bound to occur, which should also result in brief vibrating episodes. Several studies have suggested that these mechanical events are encoded by specific mechanoreceptors. In particular, Szwed et al. (2003) established that distinct populations of trigeminal ganglion neurons specifically respond to first contact and/or detachment. Other works (Ritt et al., 2008; Wolfe et al., 2008; Jadhav et al., 2009; Lottem and Azouz, 2009) have further suggested that surface roughness may in large part be encoded by the rate of discrete high acceleration events elicited by stick-slip instabilities at the whisker/object contact. Notice that analogous mechanical events, separating different action phases, have been shown to play a central role in the planning and control of manipulation tasks in the context of active human touch (Johansson and Flanagan, 2009).

In the present article, we focused on the encoding of the first contact and detachment events by the whisker. A biomechanical model is developed to quantitatively predict the sequence of whisker deformation that these events elicit in realistic conditions of tactile exploration. Specific predictions of the model were validated through high-speed video recordings of shock-induced whisker deflections. By investigating how various contextual parameters (such as radial and angular distance or friction coefficients) control the mechanical signature at the whisker base, we aimed to understand what information can be extracted by the neural system. The consequence of these results on possible tactile encoding schemes are discussed, which we believe may be of interest to the fields of whisker tactile perception as well as neurorobotics.

### **RESULTS**

Our biomechanical approach is based on a decomposition of the whisker dynamics into rapid small amplitude resonant oscillations superimposed onto a slow (quasi-static) sequence of deformation. The oscillating term is further decomposed along a series of resonant modes whose spatial and temporal properties are computed numerically.

In a first section, we describe the quasi-static evolution of the whisker as it is swept across a rectangular obstacle. In the second section, the mode decomposition scheme used to describe the rapid dynamics of the whisker is presented. Two distinct mechanical events are then successively examined: (1) the initial contact (shock) between the whisker and the object; (2) the detachment of the whisker from the object. In both cases, the precise timesequence of the whisker dynamics can be accurately predicted. A particular focus is put on the time-evolution of the moment at the whisker base (in the follicle) as it constitutes the relevant peripheral input for the mechanoreceptors.

#### **QUASI-STATIC EVOLUTION OF A WHISKER SCANNED ACROSS A RECTANGULAR OBJECT**

In this first section, the quasi-static evolution of the whisker, i.e., the series of equilibrium configurations, is calculated as it is swept across a rectangular object. Several robotics studies analyzed the shape of a whisker submitted to a contact force. Most of those works used numerical solutions to determine whisker profiles (Scholz and Rahn, 2004; Clements and Rahn, 2006; Solomon and Hartmann, 2010). In the limit of small deflections, a small angle approximation can be used which yields an analytical solution to this mechanical problem (Birdwell et al., 2007). A large majority of those works, however, ignored any frictional interaction between the whisker and the substrate in their theoretical derivations. Nevertheless a few studies performed with artificial whiskers investigated the influence of frictional interactions for predicting radial distance (Solomon and Hartmann, 2008) and local object shape (Schroeder and Hartmann, 2012) during whisker touch. Here, we numerically derive the quasi-static sequence of whisker deflection as it is swept across an obstacle, taking into account the frictional force. In line with physiological observations made in rats (Voges et al., 2012) and other whiskerbearing animals (Williams and Kramer, 2010), the whisker is modeled as a truncated tapered rod. We denote *L* the length of the non-truncated cone, *b* the maximum (base) radius and α = *b/L* the cone angle (**Figure 2A**). In order to simplify the equations governing the whisker mechanics, the whisker profiles are described using a curvilinear coordinate *s* defined as the normalized arc length such that the cone tip (the end of the non-truncated whisker) position defines the origin *s* = 0 and the whisker base is located at *s* = 1. The radius of the whisker at position *s* thus reads *r(s)* = α*sL*. In this coordinate, the whisker

physical tip is located at *s*tip = 1 − *L*whisker*/L*. Note that, owing to the truncation of the whiskers, *L* can be significantly larger than the actual whisker length *L*whisker.

The whisker profile is described in curvilinear angular coordinates as θ = θ*(s)* (**Figure 2B**). The whisker is embedded in the pad tissue down to 4–5 mm (Ebara et al., 2002), which strongly constrains the angular position of the whisker base. Although the pad does exhibit some level of elastic compliance, as established by Jenks et al. (2010) in awake rats (Jenks et al., 2010), here we assume the anchorage to be strictly rigid and thus impose θ*(*1*)* = θ*b*. The whisking process is modeled by imposing a rotation of the whisker around its base at constant rate γ. The base angle θ*b* = θ*(s* = 1*)* thus grows linearly with time such that θ*b(t)* = γ*t*. The whisker is considered linearly elastic with a uniform Young's modulus *E* and density ρ. Any intrinsic (spontaneous) curvature and out-of-plane deformations are also ignored. When in contact with the object, the whisker is submitted to a frictional force *F* which is assumed to apply at a single point located at *s* = ε along the whisker. The present analysis is restricted to configurations where ε *> s*tip, which constrains the whisker to be locally tangent to the object surface. The orientation of the force *F* with respect to the direction normal to the whisker at *s* = ε is set by the friction angle φ = tan*(*μ*)* where μ is the friction coefficient (Persson, 2000). Within these hypothesis, the equation governing the whisker bending moment equilibrium reads (see "Methods" for the detailed derivation):

$$(s^4\theta')' - \widetilde{F}\cos(\phi - \theta + \theta(\varepsilon)) = 0\tag{1}$$

where -*F* = 4*F/* πα4*EL*2 is an adimensional force. The single contact point hypothesis further imposes that the moment is null at the contact point such that θ *(*ε*)* = 0. For given values of θ*b* (imposed by the rotation of the whisker base), the friction angle φ and the contact point ε, one can numerically compute a series of equilibrium whisker profiles by imposing different values of the contact angle θ*(*ε*) <* θ*b*. This method is first used to derive the quasi-static evolution of the whisker as it rolls over the edge of the rectangular obstacle. In this regime, the position of the contact point on the object is fixed. This condition yields, for each whisking angle θ*b*, a unique solution associated with a contact location ε along the whisker. As θ*b* increases, the contact angle θ*(*ε*)* decreases and eventually vanishes. This time marks the onset of a second phase during which the whisker slides along the surface of the object. In this second regime, the contact angle θ*(*ε*)* is null while the radial distance *D* (the distance from the whisker base to the free surface) remains constant. Again, a unique solution, associated with a contact location ε, is obtained for each whisking angle θ*b*. When the contact point reaches the edge of the object, the whisker snaps off and then continues to rotate in air at a constant rotation rate.

**Figure 3A** shows the quasi-static evolution of the whisker profile for a friction coefficient μ = 0*.*4, a radial distance *D* = 0*.*83 L, and a rotation rate γ = 400◦/s. During this sequence, the contact point position along the whisker varies within a small range 0*.*08 *<* ε *<* 0*.*15 (**Figure 3A2**). The graphs **A3** and **A4** display the evolution of the base moment κ*(t)* = θ *(s* = 1*)* and its time-derivative κ˙*(t)*. In all graphs, and throughout the article, double scales are used in order to show the data both in reduced (all lengths being normalized by *L*) and physical units. For the latter, a typical whisker length *L* = 3 cm is used. **Figure 3B** displays similar traces for various friction coefficients μ = {0*.*2*,* 0*.*4*,* 0*.*6*,* 0*.*8}. This range should encompass most physical situations (Persson, 2000). Although the associated profiles appear quite similar, increasing the friction coefficient yields a significant amplification of the base moment signal κ*(t)* as shown in **Figure 3B4**. One may notice that the effect becomes significant when the whisker is sufficiently deformed while the different graphs collapse in the early moments following the initial contact, i.e., when the whisker is essentially straight. The friction coefficient also controls the time (or base angle) at which the whisker detaches from the object. Although these frictional effects are significant, the base moment κ*(t)* appears to be mostly controlled by the radial distance *D* as shown in **Figure 3C**. Reducing this distance by 10% yields a five-fold increase of the maximum base moment experienced during the exploratory sequence. It thus seems unlikely that this quantity provides significant cue for the discrimination of surfaces exhibiting different frictional properties.

#### **WHISKER RESONANT DYNAMICS**

The preceding section addressed the quasi-static evolution of the whisker as it is swept across a rectangular object. This sequence is expected to be valid for a massless whisker or at infinitely slow scanning speed (Quist and Hartmann, 2012). For a real whisker, however, inertia effects will induce significant deviations. In particular, the contact and detachment processes, which mark the

transition between distinct mechanical conditions at the whisker tip, will trigger brief episodes of oscillations.

These dynamic modulations are treated perturbatively in the form of a small displacement *u(s,t)* normal to the quasi-static profile sequence (see **Figure 4A**). We assume that both the quasistatic deformation and dynamic oscillations amplitude remain sufficiently moderate such that a small angle approximation (θ*(s)* 1) can be implemented. The validity of this hypothesis, for both mechanical events, will be discussed *a posteriori*. In this limit, the classical Euler–Bernoulli equation that governs the force equilibrium normal to the whisker can be expanded around the quasi-static profile, yielding for *u(s,t)* (see Weaver et al., 1990):

$$\frac{\partial^2}{\partial s^2} \left( EI \frac{\partial^2 u}{\partial s^2} \right) + \rho A \frac{\partial^2 u}{\partial t^2} = 0 \tag{2}$$

where *A* = π*r*2 is the whisker section area. This equation is rewritten in reduced coordinates (all distances being expressed in unit of the ideal non-truncated whisker length *L*) in the form:

$$\frac{\partial^2}{\partial s^2} \left( s^4 \frac{\partial^2 u}{\partial s^2} \right) + k^2 s^2 \frac{\partial^2 u}{\partial t^2} = 0 \tag{3}$$

where *k* = 2 √ρ*/E L/*α is a time-scale characterizing the mechanical resonance of the isolated whisker: the fundamental resonance frequency of the freely vibrating whisker reads *f*FRF = 1*.*39/k. In this expression, all lengths are expressed in units of *L*. Equation (3) is classically solved by separation of time and space variables: *u(s,t)* = *V(s)q(t)*. The spatial term *V(s)* obeys the following equation:

$$(s^4V'')'' - k^2 \omega^2 s^2 V = 0\tag{4}$$

amplitude deformation *u(s,t)* normal to the quasi-static profile. **(B)** First spatial modes for a whisker in contact at *s* = ε = 0*.*1 with boundary conditions *V (*ε*)* = *V (*ε*)* = 0. **(C)** Resonant reduced angular frequency and resonant frequency (double scale) for the first two modes as a function of the contact location ε. **(D)** First spatial modes for a freely oscillating whisker (*V (*ε*)* = *(s*4*V )(*ε*)* = 0).

The choice of boundary conditions (discussed below) sets the series of admissible angular frequencies ω*i(i* = 1*,* 2*,...)* = κ*i/k* and corresponding resonant spatial modes *Vi(s)(i* =1*,* 2*,...)*. Each mode is associated with a harmonic equation of motion that reads:

$$
\ddot{q} + \alpha\_i^2 q = 0 \tag{5}
$$

In the absence of any knowledge on the underlying mechanism, dissipative processes are accounted for by introducing a linear damping term in the time-dependent Equation (5) with a mode-independent damping ratio ζ. The dynamic equation is thus rewritten as:

$$
\ddot{q} + 2\zeta\omega\_i \dot{q} + \omega\_i^2 q = 0 \tag{6}
$$

Within this hypothesis, the general solution for the freely oscillating whisker finally reads:

$$u(s,t) = \sum\_{i} V\_i(s) \left(\alpha\_i \cos(\sqrt{1-\xi^2}\omega\_i t) + \beta\_i \sin(\sqrt{1-\xi^2}\omega\_i t)\right) e^{-\xi\omega\_i t} \tag{7}$$

#### **DEPENDENCE OF THE WHISKER RESONANCE FREQUENCY ON THE CONTACT POINT LOCATION**

The resonant spatial modes *Vi(s)* and associated angular frequencies ω*i* depend on the boundary conditions. As already mentioned, the whisker is assumed to be rigidly anchored at its base, which imposes *V(*1*)* = *V (*1*)* = 0. The different phases of the exploration correspond to distinct boundary conditions at the whisker tip.

For each value of the contact location ε in the range 0*.*01 *<* ε *<* 0*.*2, the first five modes *V*1−5*(s)* and associated adimensional angular frequencies κ1−5 are numerically computed using Mathematica v8.0 (Wolfram Research). **Figures 4B**,**D** display the spatial resonant modes obtained for ε = 0*.*1 and boundary conditions (a) and (c), respectively. Their strong asymmetry results from the tapered geometry of the whiskers. The resonant normalized angular frequencies κ1 and κ2 are shown in **Figure 4C** as a function of ε for boundary conditions (a). Both angular frequencies are found to increase as the contact point moves toward the base. The arrows indicate the corresponding angular frequencies for an isolated (freely vibrating) whisker. In order to express the resonant frequencies *fi* = κ*i/(*2π*k)* in physical units (right axis scale), a typical time-scale *k* is computed using data from the literature (Hartmann et al., 2003). This value *k* = 30*.*4 ms is conserved throughout the article.

The slender geometry of the whisker confers it the property of a resonant oscillator, which thus acts as a mechanical band-pass filter. Neimark et al. observed resonance when shaking the whisker near its tip with a piezoelectric actuator (fixed–fixed boundary conditions) (Neimark et al., 2003). Although the boundary conditions were different between those data (fixed–fixed) and our model (fixed-pinned), the frequencies they reported for β and C1 fell within our predicted range (**Figure 4C**). Owing to their various lengths, the resonant frequencies of the freely vibrating whiskers span a wide range across the pad. This observation led Neimark et al. to propose a tonotopic scheme for texture encoding in which each whisker would transduce one particular spatial wavelength (Neimark et al., 2003). The present model may in part explain why this encoding scheme hypothesis failed to receive experimental validation so far. The effective resonant frequencies varies with the location of the contact point along the whisker and the way that the whisker is pinned and/or fixed at the contact point (see **Figure 4C**). The optimal transduction frequencies in real sensing situations, rather than being whisker specific, are thus expected to vary by a factor of up to 3 over the course of a single sensing task. It is tempting to suggest in reverse that the instantaneous resonant frequency may used by the rat to extract information about the distance from the pad to the touched object. However, the dependence being relatively weak (owing to the tapered geometry of the whisker), the spectral characteristic of the whisker dynamic is unlikely to play a significant role in the precise determination of radial distance.

#### **SHOCK AGAINST THE OBJECT'S EDGE**

The resonant modes are now used to investigate the dynamics induced by the shock of the whisker against the object's edge. The whisker initially rotates in air at constant angular velocity γ around the whisker base. At time *t* = 0, the whisker makes contact with the object's edge at a position *s* = ε along the whisker. In line with previous observations (Hartmann et al., 2003), the collision is assumed to be inelastic such that the whisker tip remains in constant contact with the object at *t >* 0 (no rebound). The duration of the shock-induced oscillation is expected to damp out over a time period of order τ*/*ζ, where τ is the period of the fundamental mode. During this time, the contact point location ε, deduced from the quasi-static sequence, varies by less than 0*.*02 (0*.*002 after one resonant period). We ignore this minute change and assume ε to be constant, which allows one to describe the oscillating dynamics on a well-defined series of resonant modes.

The whisker profile *U(s,t)* at time *t >* 0 is decomposed as *U(s,t)* = *U*qs*(s,t)* + *u(s,t)* where *U*qs*(s, t)* is the quasi-static profile evolution and *u(s,t)* characterizes the shock-induced dynamics. By introducing this decomposition in Equation (3), *u(s,t)* is found to obey the dynamic equation (see "Methods"):

$$\frac{\partial}{\partial s^2} \left( s^4 \frac{\partial^2 u}{\partial s^2} \right) + k^2 s^2 \frac{\partial^2 u}{\partial t^2} = -k^2 s^2 \ddot{U}\_{\text{qs}}(s, t) \tag{8}$$

One thus needs to compute the second time-derivative of the quasi-static profile *U*¨ qs*(s,t)*. Prior to the shock, the whisker experiences a solid rotation at constant rotation rate γ. Immediately after the shock, within the small deflection approximation, the quasi-static profile also evolves linearly in time (see "Methods"). The expression of *U*¨ qs*(s,t)* can thus be written as:

$$
\ddot{U}\_{\rm qs}(s, t) = 8(t)\chi \ddot{U}(s) \tag{9}
$$

where δ*(t)* is the Dirac function and *U*¯ *(s)* is a normalized profile that depends on ε and γ (see "Methods" for analytical derivation). The dynamic term *u(s,t)* is decomposed onto the resonant modes *Vi(s)* corresponding to boundary conditions *V(*ε*)* = *V(*ε*)* = 0, in the form *u(s,t)* = *i qi(t)Vi(s)*. Projecting expression (Equation 8) onto each mode *Vi(s)*, and using the orthogonality of the resonant modes, the mode amplitudes *qi(t)* are found to obey the dynamic equation [see Equations (36–41) in "Methods"]:

$$\ddot{q}\_i + 2\zeta\omega\_i\dot{q}\_i + \omega\_i^2 q\_i = -\int\_{\varepsilon}^{1} s^2 \dot{U}\_{\text{qs}}(s, t) V\_i(s) ds$$

$$= -\delta(t) \chi \int\_{\varepsilon}^{1} s^2 \bar{U}(s) V\_i(s) ds \qquad (10)$$

Notice that, as indicated before, a damping term is added to account for dissipative processes. After integration, the complete shock-induced dynamics reads:

$$U(s, t > 0) = \gamma t \,\bar{U}(s) - \gamma \sum\_{i} V\_i(s) \left( \int\_{\varepsilon}^{1} s^2 \,\bar{U}(s) V\_i(s) ds \right)$$

$$\times \frac{e^{-\xi \alpha t} \sin \left( \sqrt{1 - \xi^2} \alpha\_i t \right)}{\sqrt{1 - \xi^2} \alpha\_i} \tag{11}$$

**Figure 5A** shows the successive whisker profiles for a contact located at ε = 0*.*2. The damping factor is set at ζ = 0*.*1, consistent with values reported in the literature for the fundamental mode (Hartmann et al., 2003; Neimark et al., 2003). **Figure 5B** displays the same sequence in the reference frame of the whisker base (the imposed rotation of the whisker base has been subtracted). This graph illustrates how the shock at the tip of the whisker triggers a wave of deflection that travels up to the whisker base. Notice that the first maximum (indicated by an arrow in **Figure 5B**) is negative, i.e., opposite to the long time scale whisker deflection induced by the object. The position of the first minimum displays a linear dependence with the time elapsed since the shock (**Figure 5C**). This constant wave velocity results from the tapered geometry, since a √*t* dependence is expected in the case of a cylindrical rod (Audoly and Neukirch, 2005). From dimensional analysis of Equation (3), the wave velocity is expected to be of order *c*wave = α*c* where α is the cone angle and *c* = √*E/*ρ is the sound velocity in the bulk material. The apparent wave velocity obtained by linear fitting on graph C is consistently found to be 2*.*54*c*wave (5*.*02 m/s in physical units). **Figure 5D** shows the evolution of the moment κ*(t)* = *∂*2*U/∂s* 2*(s* = 1*, t)* at the base of the whisker (black solid line). The graph also displays both the dynamic (κdyn*(t)* = *∂*2*u/∂s* 2*(s* = 1*,t)*, gray solid line) and quasi-static (κqs*(t)* = *∂U*2 qs*/∂s* 2*(s* = 1*, t)*, dotted line) components for comparison. Although the maximum amplitude of the deflection wave is of the order of 15μm, i.e., a small fraction

shock-induced oscillations are visible through the varying density of the profiles. **(B)** Whisker profiles in the reference frame of the whisker base (δ*t* = 0*.*2 ms). The distance indicated on each graph corresponds to the

2*.*54*c*wave . **(D)** Time-evolution of the quasi-static κqs*(t)* (dotted line) and dynamic κdyn*(t)* (solid line) base moment. **(E)** Evolution of the time-derivative of the base moment κ˙*(t)*.

of the whisker base diameter, it yields a significant negative dip in the whisker base curvature dynamics. The time-derivative κ˙*(t)* signal (**Figure 5E**) in turn exhibits a clear signature of the whisker/object contact in the form of large amplitude oscillation with a peaked maximum occurring at a time τpeak ≈ 4 ms after the shock.

Equation (11) indicates that the peak amplitude κ˙max should be linearly proportional to the rotation rate γ. It also increases as the contact point moves toward the whisker base i.e., for shorter radial distances (**Figures 6A**–**C**). Notice that the peak amplitude remains larger than the quasi-static value κ˙qs*(*τpeak*)* (Inset). The last graph of **Figure 6** shows the delay τpeak between the shock event and the arrival of the mechanical signal at the whisker base as a function of the contact location. The dependence is quasi-linear within the range of ε explored, which allows us to extract an approximate velocity 0*.*69*c*wave (1*.*36 m/s in physical units) of the same order as the wave velocity determined in **Figure 5**.

In the range of parameters explored, the maximum angular deflection of the whisker during the process is 0*.*06. This value validates the small angle hypothesis underlying the present analysis (Birdwell et al., 2007).

#### **EXPERIMENTAL MEASUREMENTS OF SHOCK-INDUCED OSCILLATION**

These predictions are tested experimentally on an anesthetized rat using high-speed videography (see "Methods"). The rat, and thus the whisker base, are maintained fixed while a thin bar is moved at constant speed *V*bar = 60 mm/s against the whisker tip (**Figure 7A**). We used a constant speed and not a tap since in that case a sudden acceleration would happen, which is a different condition and ethologically more unlike to happen. The bar is vertical and fixed by its edge to a rectilinear motor. During the early instants following the shock, which are analyzed here, the whisker did not slip onto the bar such that the contact point along the whisker can be considered invariant. Whisker movements are captured using a bird's-eye view high-speed camera operating at 2.5 kHz (**Figure 7A**). Whisker centerline profiles are tracked within each frame between the fur and the bar with a custom-designed semi-automatic script. The shock-induced whisker deflection profile *U(s,t)* is obtained by subtracting, for each frame, the whisker intrinsic profile as determined from images recorded before the shock.

The resulting sequence of shock-induced whisker deflections is shown in **Figure 7B** for a β whisker (length 46 mm) during the first 3 ms from the instant of collision. Consistent with the

model, a deflection wave is observed characterized by a negative maximum deflection propagating upward. As expected, the wave appears to propagate at constant speed as indicated by the linear dependence of the location of the maximum deflection with time (**Figure 7C**). This allows us to extract an effective wave velocity equal to 5.6 m/s, in close agreement with the typical value predicted before.

In order to more quantitatively compare these results with the biomechanical model, we decomposed the whisker dynamics into rapid resonant oscillations superimposed onto a slow quasistatic sequence of deformation imposed by the moving bar. The quasi-static evolution *U*qs*(s,t)* is first evaluated using the profiles measured at long time scale, i.e., when the relative contribution of the dynamic oscillation is expected to be negligible (*t >* 8 ms). As expected, in this regime, *U*qs*(s,t)* can be linearized as *U*qs*(s,t)* = *V*bar*tU*¯ *(s)*. The normalized profile *U*¯ *(s)* is well fitted by the static theoretical profile in the limit of small deflection (see "Methods"). The best fit yields a value of ε = 0*.*25. Using the long time-scale normalized profile *U*¯ *(s)*, the complete sequence of whisker deflection *U(s,t)* is computed using the same scheme as described earlier [adapted to the linear displacement configuration (see "Methods")]. The comparison, for each parameter, is shown in **Figures 7B**–**D**. The model quantitatively captures not only the wave propagation dynamics (*r*2 = 0*.*96, *p <* 2*.*10−9) (**Figure 7C**), but also the amplitude of the maximum deflection (*r*2 = 0*.*94, *p <* 2*.*10−8) (inset in **Figure 7C**) and, most importantly, the base curvature signal (*r*2 = 0*.*94, *p <* 2*.*10−8) (**Figure 7D**).

Identical measurements were performed for a C1 whisker (length 36 mm; ε = 0*.*15). The wave dynamics was found consistent with the model prediction (*r*2 = 0*.*95, *p <* 10−5) with a wave velocity of 5.8 m/s, close to the value obtained for the β whisker. The measured wave amplitude appeared significantly lower than predicted during the first ms following the shock, resulting in low values of correlation for this parameter (*r*2 = 0*.*41, *p* = 0*.*06). After this initial period, however, a very consistent match was recovered (*r*2 = 0*.*97, *p <* 5*.*10−5). The measured curvature signal at the whisker base being rather unaffected by the early deflections of the whisker tip, it was found to agree with the prediction for all the duration of the process (*r*2 = 0*.*65, *p* = 0*.*0027).

of the whisker/object initial contact (sampling rate 2*.*5 kHz) between the vertical bar and the β whisker. The result of the whisker tracking is superimposed (white line) on the frame. The direction of movement of the bar is indicated by an arrow. **(B)** Whisker profiles in the reference frame of the whisker base (δ*t* = 0*.*8 ms). The bottom panel is a magnified view of the

### maximal deformation as a function of the time elapsed since the shock (see arrow in **B**). The gray line is the best linear fit and corresponds to a wave velocity 5*.*6 m/s. **(D)** Time-evolution of the experimental (solid line) and predicted (gray line) κ*(t)* base curvature signal.

#### **CONSEQUENCE FOR EVENT-BASED OBJECT POSITION ENCODING**

As the rat repetitively whisks onto an object, it produces a series of shocks. The present work demonstrates that each of them triggers, a few ms after contact, a characteristic signature in the base curvature signal which can be quantitatively predicted using a first order mechanical model of the whisker. The resulting mechanical stimulation at the whisker pad should be sufficiently intense to trigger a clear neural response. As displayed in **Figure 7**, the whisker peak base curvature is of the order of 10−3 mm−1 and its maximum time-derivative varies up to 1 mm−1*.*s−1. Stimuli of comparable intensity have been shown to elicit reliable cortical discharges (O'Connor et al., 2010b; Huber et al., 2012).

Notice that the peak in the base curvature time derivative has an opposite sign and a larger amplitude compared to the longtime component, which may explain how so-called touch cells may specifically respond to first contact (Szwed et al., 2003). It has been recently proposed that these cells may mediate the coding of object angular position with respect to the pad through the precise timing of the shock event within the whisking cycle (Knutsen and Ahissar, 2009). The present biomechanical analysis allows us to estimate how the pre-neural whisker transduction contributes to the horizontal resolution of such an encoding scheme.

If one assumes that the neural response is triggered by the maximum of κ˙*(t)*, the jitter *T* in the mechanoreceptor response (the dispersion in spike timing) should be a fraction of the fundamental resonant period, i.e., a few ms. With a characteristic rotation rate during whisking of order γ = 400◦, this yields an angular resolution of γ*T* ≈ 1◦ consistent with available behavioral data (Knutsen et al., 2006). This crude evaluation implicitly ignores the lag between the shock and the arrival of the deflection wave at the whisker base that elicits the mechanoreceptors' response. This delay was shown to be of the order of *c*wave*L*contact-base where *L*contact-base is the arc length between the contact point and the whisker base. This lag effect induces an additional error of up to a few degrees on the angular position of the object if its radial distance is unknown.

One may therefore suggest that the radial distance is evaluated using a parallel coding channel. One possible scenario, as originally proposed by Szwed et al. (2003), relies on the intensity of the mechanical signal which appears to decay rapidly with the radial distance, as shown in **Figure 6C**. By combining both information (timing and intensity of the shock-induced whisker base mechanical signal), a precise localization of the object can be recovered (Knutsen and Ahissar, 2009).

#### **DETACHMENT**

We now turn to the detachment process that occurs when the whisker tip reaches the second edge of the object. We note *t* = 0 the time at which the whisker snaps off. At *t >* 0, the whisker detaches from the object then oscillates freely in air. The profile is decomposed as *U(s,t)* = *U*qs*(s,t)* + *u(s,t)* where *U*qs*(s, t)* is the quasi-static evolution. The continuity of the position and velocity profiles at time *t* = 0 imposes:

$$
\mu(s, t = 0) = \Delta U\_{\rm qs}(s) \tag{12}
$$

$$
\dot{u}(s, t=0) = \Delta \dot{U}\_{\text{qs}}(s) \tag{13}
$$

where *U*qs*(s)* = *U*qs*(s,* 0−*)* − *U*qs*(s,* 0+*)* and *U*˙ qs*(s)* = *U*˙ qs*(s,* 0−*)* − *U*˙ qs*(s,* 0+*)*. The term *u(s,t)* is decomposed along the resonant modes with boundary conditions *V(*ε*)* = *(s* 4*V) (*ε*)* = 0 as in Equation (7):

$$u(s,t) = \sum\_{i} V\_i(s)(\alpha\_i \cos(\sqrt{1-\xi^2}\omega\_i t))$$

$$+ \beta\_i \sin(\sqrt{1-\xi^2}\omega\_i t))e^{-\xi\alpha\_i t} \tag{14}$$

Projecting Equations (12) and (13) on each spatial mode, the coefficients α*i* and β*i* can be written, in the limit ζ 1:

$$\alpha\_{i} = \int\_{\varepsilon}^{1} s^{2} \Delta U\_{\text{qs}}(s) V\_{i}(s) ds \tag{15}$$

$$\beta\_i = \frac{1}{\alpha\_i} \int\_{\varepsilon}^{1} s^2 \Delta \dot{U}\_{\text{qs}}(s) V\_i(s) ds \tag{16}$$

In order to obtain an expression of *U*qs and *U*˙ qs, the quasistatic profiles are computed at time {−δ*t,* 0−*,* 0+*,* δ*t*} with δ*t <<* 1*/*γ. This allows us to obtain linearized expressions of the quasistatic profile sequences around *t* = 0, before and after the shock, in the form:

$$U\_{\rm qs}(\text{s, t} < 0) = U\_{\rm qs}(\text{s, 0}^{-}) + \frac{U\_{\rm qs}(\text{s, 0}^{-}) - U\_{\rm qs}(\text{s, } - \text{\&t})}{8t}t \quad (17)$$

$$U\_{\rm qs}(s, t > 0) = U\_{\rm qs}(s, 0^+) + \frac{U\_{\rm qs}(s, \delta t) - U\_{\rm qs}(s, 0^+)}{\delta t}t \tag{18}$$

This yields:

$$
\Delta U\_{\rm qs}(s) = U\_{\rm qs}(s, 0^{-}) - U\_{\rm qs}(s, 0^{+}) \tag{19}
$$

$$
\Delta \dot{U}\_{\rm qs}(s) = \frac{1}{8t} \left( U\_{\rm qs}(s, 0^{-}) - U\_{\rm qs}(s, -\delta t) - U\_{\rm qs}(s, \delta t) \right)
$$

$$
+ U\_{\rm qs}(s, 0^{+}) \Big) \tag{20}
$$

**Figure 8A** shows the evolution of the whisker profile just after the detachment, for different base angle θ*b* at detachment. When the quasi-static evolution is subtracted (**graph B**), one recovers a wave propagation mechanism qualitatively similar to that observed after the shock (wave speed of order 3.5 m/s). This event produces a characteristic signature at the whisker base shown in **Figures 8B**,**C**. Notice that the mechanical signal is rather insensitive to the whisker base angle prior to the detachment.

As indicated earlier, the present analysis is based on the assumption that the whisker remains weakly deformed during the whole process. As shown in **Figure 8A**, this hypothesis becomes valid only a few hundreds of microseconds after the event is triggered. At very early time, and for large angular and small radial distances to the object, the tip of the whisker exhibits large deflection angle. However, owing to the conical shape of the whisker, the large deformation appears to be confined to the very end of the whisker: in the range of configurations explored, 80% of the whisker, i.e., 99% of the mass, displays a deflection angle smaller than 0.5. The small angle approximation should therefore provides a reasonable first order approximation of the detachment dynamics for all situations explored.

The whisker detachment from the object also triggers a deflection wave that propagates upward at a speed comparable to that observed for the shock event. The resulting signal at the whisker base appears to be ≈10 times more intense than for the shock. It should thus trigger significant neural response, as confirmed by the existence of so-called detachment cells whose response is specifically triggered by such events (Szwed et al., 2003). Our study shows that an increase in the friction coefficient delays the detachment process (**Figure 3B**). The precise timing of the detachment-induced signal may thus provide an indication on the surface frictional properties. Such a scenario would however, require the independent knowledge of the angular position of the object edge.

#### **VIBRISSAE CAN BE ASSUMED TO BE UNIFORMLY DAMPED OVER ALL SPATIAL MODES**

A key assumption in the development of the present model is that the damping coefficient of the vibrissa is constant over all spatial modes. Previous experimental studies, however, have shown that magnification ratio increases with mode number (Hartmann et al., 2003; Neimark et al., 2003). Specifically, the bottom plot of Figure 4 of Hartmann et al. (2003), and Figure 4 of Neimark et al. (2003) (e.g., the β whisker) both show an increase in peak magnitude with mode number. Note that Neimark et al. chose to plot the curve assuming constant amplitude (instead of constant acceleration), so each peak must be scaled by the frequency squared to demonstrate the increase with mode number. This experimental finding has until now been unexplained, although Hartmann et al. (2003) suggested that it was likely in part due "to nonlinear effects of viscous and/or hysteretic damping."

Here we show that, at odds with the authors' proposed interpretation, a linearly damped model of whisker with a modeindependent damping ratio, as assumed in the present work, allows one to correctly capture this measurement. It also illustrates how the mode-decomposition approach introduced in the present study, may allow to quantitatively re-analyze earlier data on whisker dynamics.

Throughout the calculation, we use the same notation as in the present article. The text indicates that the whisker has a length *L*whisker = 53 mm, a base radius *b* = 105μm and a tip radius *r*tip = 5*.*5μm. Assuming a conical shape, this yields = *r*tip*/b* = 0*.*05. Using the Young modulus *E* = 3*.*02 GPa and ρ = 1*.*14 mg/mm3 indicated in the article, the values for the (undamped) resonant frequencies can be computed and read: *Fi* = {37*.*9*,* 93*.*4*,* 175*.*4 Hz}. The whisker base is forced with a

sinusoidal motion [**Figure 9A** reproduced from Hartmann et al. (2003)], starting at *t* = 0:

$$A\_{\text{base}}(t) = A\_{\text{base}} \sin(\omega t) \tag{21}$$

The whisker response is derived using the same scheme as detailed in the present article. The displacement sequence is broken into *U(s,t)* = *U*qs*(s,t)* + *u(s,t)* where the first term represents the massless whisker response: *U*qs*(s,t)* = *A*base sin*(*ω*t)*. The dynamic term can thus be written as:

$$u(s,t) = \sum\_{i} V\_i(s) \int\_{\epsilon}^{1} s^2 V\_i(s) ds \int\_0^t G(t-t') \dot{U}\_{\text{qs}}(s,t') dt' \quad \text{(22)}$$

$$= -A\_{\text{base}} \sum\_i V\_i(s) \int\_{\epsilon}^{1} s^2 V\_i(s) ds \int\_0^t G(t-t') \alpha^2$$

$$\times \sin(\alpha t') dt' \tag{23}$$

where *G(t)* is the Green's function of the resonant system (Equation 49).

The time-dependent term in Equation (23) can be analytically integrated. It oscillates at the driving frequency with an absolute peak-to-peak amplitude that can be written as 2*F(*ω*/*ω*i,*ζ*)*. Due to the resonant nature of the system, the spectrum is dominated by the closest resonant mode. We thus approximate the maximal amplitude of the response signal with the sum of the maximal amplitude of each mode signal. The ratio between the whisker tip oscillation amplitude (the peak-to-peak amplitude of signal *U(s,t)* measured at *s* = ε) and the driving base amplitude 2*A*base thus reads:

$$R(\boldsymbol{\omega}) \approx 1 + \sum\_{i} \left[ V\_i(\boldsymbol{\epsilon}) \int\_{\boldsymbol{\epsilon}}^{1} s^2 V\_i(\boldsymbol{s}) ds \right] F(\boldsymbol{\alpha}/\boldsymbol{\omega}\_i, \boldsymbol{\xi}) \tag{24}$$

The dependence of the maxima of the response signal with mode number thus depends primarily on the series of prefactors *Vi()* 1 *s* 2*Vi(s)ds*, such that no conclusion can be drawn without prior knowledge of the resonant spatial modes. In the present case, this prefactor reads, for the first three modes: {3*.*7*,* 6*.*5*,* 7*.*6} and thus does increase with the mode number. Notice that the evolution of the peak amplitude results from the conical shape of the whisker—it would decrease for a cylindrical rod. **Figure 9B** compares the prediction obtained for a damping ratio of 0*.*15 (a characteristic value given by the authors) identical for all three modes. Although no adjustable parameter has been used to produce this plot, it does captures the evolution of the whisker tip

response. The slight over-estimation of the two first modes can be due to the fact that the whisker displays significant spontaneous curvature, such that the small angle limit is not strictly valid.

It thus appears that the data reported in previous experimental studies (Hartmann et al., 2003; Neimark et al., 2003) are, to first order, consistent with a linear damping model with a unique damping ratio for all modes.

#### **DISCUSSION**

The biomechanical model presented here provides a generic method to describe the mechanical transduction operated by the whisker in realistic conditions of exploration. The theoretical approach consists of successively computing the slow (quasistatic) evolution of the whisker and the rapid resonant dynamics. It provides prediction of the complete whisker deflection, and in particular of the mechanical signal elicited in the follicle which conveys the relevant peripheral input for the mechanoreceptors. It is therefore directly amenable to experimental comparison as illustrated in the present article in the case of first contact. In this configuration, we were able to characterize the dynamic of the deflection wave with micron-scale resolution by subtracting the static configuration of the whisker prior to the event. For the detachment case, the deflection wave proved to be difficult to observe with such precision owing to the continuous evolution of the whisker quasi-static profile prior to the event.

Beyond the specific context in which it is implemented in the present study, this biomechanical model may constitute a useful tool for understanding the neural encoding of tactile information. It may in particular guide the type of stimuli that needs to be played at the whisker base during electrophysiological recordings in the somatosensory cortex. It should also find direct applications in robotic rats by easing the development of decoding algorithms allowing one to extract relevant physical information from the whisker inputs (Kim and Möller, 2007). This approach, unlike finite-element models, allows one to clearly identify the consequences of whiskers' resonant properties.

The model is based on several simplifying assumptions and should therefore be considered as a first order description of the whisker dynamics. First, the whisker spontaneous curvature and out-of-plane deflections are ignored. For moderately curved whisker, the spontaneous curvature can be accounted for when deriving the quasi-static sequence by modifying Equation (1) as explained in the "Methods" section. The prediction of the dynamic model itself does not depend on the resting state. However, for largely deformed whisker, the paradigmatic scenario envisioned here might be changed. In particular, the hypothesis of a shock normal to the whisker axis might be unrealistic for strongly concave-forward whiskers (Quist and Hartmann, 2012). For such vibrissae, the out-of-plane deflections are bound to be also significant.

Second, the whisker is assumed to exhibit a conical shape and to be linearly elastic, with uniform density and elastic modulus. It has been recently reported that the Young's modulus *E* may vary significantly along the whisker, owing to its internal structure (Quist et al., 2011), although conflicting results exist (Carl et al., 2012). However, one should notice that the modulus appears in the equation through the bending stiffness which scales as *r*4*E* (*r* whisker radius), so that the most critical hypothesis should actually lie in the conical shape assumption. Indeed, a 10% deviation to the conical shape is expected to induce a larger deviation than a 40% variation of *E* along the whisker.

Third, the model assumes linear damping, characterized by a unique parameter ζ in the dynamic equation. This *ad-hoc* introduction of a damping term is not based on any sort of physical modeling and reflects the absence of prior knowledge on the nature of dissipative processes. At least five different mechanisms may contribute to damping: visco-elastic dissipation in the whisker itself, mechanical radiation in the tissue, air friction, interfacial dissipation at the whisker/substrate contact, viscoelastic dissipation in the follicle. In the absence of well-established mechanical descriptions of the dissipative mechanism, it is common for lightly damped systems (i.e., effective damping ratio *<*0*.*2) to assume proportional damping (Weaver et al., 1990), as done in the present work. One should notice that typical damping ratios reported in the literature for the fundamental mode lies in the range 0*.*1 *<* ζ *<* 0*.*2 such that the resonant system is well within the underdamped limit. As a result, the resonant frequency is not strongly affected by the precise choice of the damping ratio. During the review process, one of the referees pointed to us that the linear damping hypothesis with a unique damping ratio appears at odds with resonance curves experimentally obtained by Hartmann et al. on isolated whiskers (Hartmann et al., 2003). When re-analyzed using the resonant modes description, these data are in fact consistent with this hypothesis as demonstrated in the "Methods" section.

Within this set of hypothesis, it appears that the whisker transduction mechanism depends on the whisker intrinsic properties through only three independent parameters: (1) the length of the equivalent non-truncated whisker *L* (2) the wave velocity *c*wave = α √*E/*ρ that characterizes the speed at which information is conveyed from the whisker tip to the follicle, and (3) the damping factor ζ that reflects the dissipative processes (the intrinsic timescale *k* is the product 2*c*wave*L*). Notice that the Young's modulus *E*, the density ρ, and the cone angle α vary in a small range across the whisker pad (Hartmann et al., 2003; Carl et al., 2012; Voges et al., 2012) so that *c*wave should weakly depend on the whisker identity. Similarly, the values reported for the damping factor appear to lie within a rather small range (Hartmann et al., 2003; Neimark et al., 2003). When displayed in reduced coordinates (all lengths being expressed in units of *L*), the present results should therefore provide a reasonable description of the rapid whisker dynamics, regardless of the whisker identity.

The last and most limiting hypothesis of the proposed dynamic model is that its validity requires the whisker deflection to remain small. For large deflections, the dynamic equation should include additional terms that reflect the coupling between the quasistatic (deformed) state and the oscillating dynamic. This correction might become significant when one attempts to describe situations in which the whisker vibrates around strongly predeformed configurations, such as encountered during exploration of textured surfaces for instance. At the expense of significant mathematical developments, it should be possible to derive the resonant spatial modes in this regime and then implement the same theoretical scheme. However, as illustrated with the detachment process, the whisker tapered shape tends to confine the large deflection to the very end of the whisker even for large angular distances. One may thus hope that this first order description provides a good approximation in most realistic exploratory conditions. We postpone the discussion of this regime of whisker oscillations during sliding to a forthcoming publication.

#### **METHODS**

#### **EQUILIBRIUM PROFILE OF A FRICTIONAL WHISKER**

The whisker is modeled as a slender tapered rod of length *L* and maximum base radius *b*. The cone angle *b/L* is noted α. The whisker base is rigidly clamped with a fixed angle θ*(s* = 1*)* = θ*b* and submitted to a frictional force *F* assumed to be applied at a single point located at *s* = ε along the whisker. The present analysis is restricted to configurations where ε *> s*tip which constrains the whisker to be locally tangent to the object surface. The orientation of the friction force *F* with respect to the direction normal to the whisker at *s* = ε is set by the friction angle φ = tan*(*μ*)* where μ is the friction coefficient (see **Figure 2B**). Forces balance along the whisker yields an expression for the normal *N(s)* and tangential force *T(s)* at location *s >* ε:

$$N(\mathbf{s}) = F \cos \left(\phi - \theta + \theta(\mathbf{s})\right) \tag{25}$$

$$T(\mathbf{s}) = F\sin\left(\phi - \theta + \theta(\mathbf{s})\right) \tag{26}$$

The momentum equilibrium further imposes:

$$\frac{\partial M}{\partial s} - N(s) = 0\tag{27}$$

where the bending moment can be expressed as *M(s)* = *EI*θ *,I* = *π*α4*s* 4 4 is the area moment of inertia and *E* is the Young's modulus. From Equations (25) and (27), one obtains the dimensionless static equilibrium equation:

$$(s^4\theta')' - \widetilde{F}\cos\left(\phi - \theta + \theta(\varepsilon)\right) = 0\tag{28}$$

where all distances are scaled by the length of the non-truncated cone *L* and the reduced force -*F* = 4*F/ π*α4*EL*2 . Notice that the intrinsic in-plane whisker curvature θ int*(s)* could be taken into account by substituting the first (derivative) term by *(s* 4*(*θ − θ int*))* . In the present study, the intrinsic whisker curvature is supposed to be null. Equation (28) is solved numerically using Mathematica v8.0 (Wolfram Research) through integration of the following equation:

$$\frac{\partial}{\partial s} \left( \frac{(s^4 \theta')'}{\cos \left( \phi - \theta + \theta(\varepsilon) \right)} \right) = 0 \tag{29}$$

We assume the whisker to be rigidly clamped in the pad such that θ*(*1*)* = θ*b*. The friction force is assumed to apply at a single location so that the bending momentum is null at the contact point: θ *(*ε*)* = 0. For any values of the friction coefficient μ, maximum angle θ*b*, and contact point location ε, one finds a unique solution when further imposing the contact angle θ*(*ε*) <* θ*b*.

For weakly deformed whiskers, a small angle approximation of Equation (29) can be used which reads, in Cartesian coordinates:

$$(s^4 U'')'' = 0\tag{30}$$

#### **ORTHOGONALITY OF THE RESONANT MODES**

Multiplying Equation (4) for mode *i* by *sVj* and integrating between ε and 1 gives the equality:

$$\int\_{\underline{\varepsilon}}^{1} \left(s^{4}V\_{i}^{\prime\prime}\right)^{\prime\prime}V\_{j}ds = k^{2}\alpha\_{i}\int\_{\underline{\varepsilon}}^{1}s^{2}V\_{i}V\_{j}ds\tag{31}$$

Integrating the left member by parts (twice) yields:

$$[(s^4V\_i'')'V\_j]\_\varepsilon^1 - [s^4V\_i''V\_j']\_\varepsilon^1 + \int\_\varepsilon^1 s^4V\_i'''V\_j'' ds = k^2 \alpha\_i \int\_\varepsilon^1 s^2 V\_i V\_j ds \tag{32}$$

The boundary conditions involve that the bracket terms are null. One can rewrite the same expression by swapping *i* and *j*. If we now subtract both equations, we obtain:

$$(\boldsymbol{\alpha}\_i^2 - \boldsymbol{\alpha}\_j^2) \int\_{\boldsymbol{\varepsilon}}^1 s^2 V\_i V\_j ds = 0 \tag{33}$$

The spatial mode *Vi(s)* are normalized such that 1 ε *s* 2*V*2 *i* = 1. This yields the orthogonality property of the resonant modes:

$$\int\_{\ell}^{1} s^2 V\_i V\_j ds = 0 \text{ for } i \neq j$$

$$\int\_{\ell}^{1} s^2 V\_i V\_j ds = 1 \text{ for } i = j \tag{34}$$

This relationship has practical consequences. First, let's con- sider a dynamic evolution of the whisker in the form: *u(s,t)* = *i qi(t)Vi(s)ds*. The amplitude of mode *i* can be directly computed as:

$$q\_i(t) = \int\_{\ell}^{1} s^2 V\_i(s) u(s, t) ds\tag{35}$$

Second, if we now consider a situation where the whisker oscillates through the application of a distributed normal force gradient *f(s,t)*, Equation (3) reads:

$$\frac{\partial}{\partial s^2} \left( s^4 \frac{\partial^2 u}{\partial s^2} \right) + k^2 s^2 \frac{\partial^2 u}{\partial t^2} = f(s, t) \tag{36}$$

Introducing the mode decomposition in this equation, we write:

$$\sum\_{i} \left[ (s^4 V\_i^{\prime\prime})^{\prime\prime} q\_i + k^2 s^2 V\_i \dddot{q}\_i \right] = f(s, t) \tag{37}$$

The modes *Vi* being solutions of Equation (4), the equation can be written as:

$$\sum\_{i} s^{2} V\_{i} \left[ \ddot{q}\_{i} + \omega\_{i}^{2} q\_{i} \right] = \frac{1}{k^{2}} f(s, t) \tag{38}$$

Projecting this equation on the spatial mode *Vi(s)* and using the orthogonality property, one obtains the dynamic equation for the mode amplitude *qi(t)*:

$$
\ddot{q}\_i + \omega\_i^2 q\_i = \frac{1}{k^2} \int\_{\ell}^1 f(s, t) V\_i(s) ds \tag{39}
$$

#### **COUPLING THE WHISKER'S RAPID DYNAMICS TO ITS QUASI-STATIC EVOLUTION**

The whisker displacement *U(s,t)* is decomposed in the form *U(s,t)* = *U*qs*(s,t)* + *u(s,t)* where *U*qs*(s,t)* describes the quasi-static evolution. By introducing this decomposition in Equation (3), one obtains the dynamic equation for *u(s,t)*:

$$\frac{\partial^2}{\partial s^2} \left( s^4 \frac{\partial^2 u}{\partial s^2} + s^4 \frac{\partial^2 U\_{\text{qs}}}{\partial s^2} \right) + k^2 s^2 \left( \frac{\partial^2 u}{\partial t^2} + \frac{\partial^2 U\_{\text{qs}}}{\partial t^2} \right) = 0 \tag{40}$$

We assume *U*qs*(s,t)* to obey the small deflection equilibrium Equation (30): *(s* 4*U* qs*)* = 0. The resonant component of the whisker deflection thus obeys the dynamic equation:

$$\frac{\partial^2}{\partial s^2} \left( s^4 \frac{\partial^2 u}{\partial s^2} \right) + k^2 s^2 \frac{\partial^2 u}{\partial t^2} = -k^2 s^2 \dot{U}\_{\text{qs}}(s, t) \tag{41}$$

Using the mode-decomposition *u(s,t)* = *i qi(t)Vi(s)* and the orthogonality of the resonant modes [see Equations (36) and (39)], the amplitude of each mode *qi(t)* is found to obey the dynamic equation:

$$
\ddot{q}\_i + 2\xi\omega\_i \dot{q}\_i + \omega\_i^2 q = -\int\_{\varepsilon}^{1} s^2 \dot{U}\_{\text{q}^\natural}(s, t) V\_i(s) ds \tag{42}
$$

This approach is implemented in the case of a shock. Prior to the shock, the whisker experiences a solid rotation at constant rotation rate γ which simply reads: *U*qs*(s,t <* 0*)* = γ*(*1 − *s)t*. Integrating Equation (30) with boundary conditions *U(*1*)* = 0*, U (*1*)* = γ*t, U(*ε*)* = *U(*ε*)* = 0, one obtains the evolution of the quasi-static profile immediately after the shock:

$$U\_{\rm qs}(s, t > 0) = \gamma t \frac{(1 - s) \left(2s^2 + \varepsilon^2 (1 + s) - \varepsilon s(3 + s)\right)}{2(1 - \varepsilon)^2 s^2} \tag{43}$$

Based on these two expressions, one can compute the second time-derivative of the quasi-static profile *U*¨ qs*(s,t)*:

$$
\ddot{U}\_{\rm qs}(s, t) = \delta(t)\gamma \ddot{U}(s) \tag{44}
$$

where δ*(t)* is the Dirac function and *U*¯ *(s)* is a normalized profile that reads:

$$\ddot{U}(\mathbf{s}) = \dot{U}\_{\rm qs}(\mathbf{s}, \mathbf{0}^{+}) - \dot{U}\_{\rm qs}(\mathbf{s}, \mathbf{0}^{-}) = \frac{\varepsilon (1 - \mathfrak{s})^{2} (\mathfrak{e} - \mathfrak{z}\mathfrak{s} + 2\mathfrak{s}\mathfrak{e})}{2(1 - \mathfrak{e})^{2}\mathfrak{s}^{2}} \tag{45}$$

In order to integrate Equation (42), one needs to compute the Green's function *G(t)*, i.e., the response of the dynamic system, initially at rest, to a Dirac of unit force:

$$\left(\frac{d^2}{dt^2} + 2\zeta\alpha\_i\frac{d}{dt} + \alpha\_i^2\right)G(t) = \delta(t) \tag{46}$$

with boundary conditions *G(t <* 0*)* = *G (t <* 0*)* = 0. For *t >* 0, the general solution reads:

$$G(t>0) = \left(\alpha\cos(\sqrt{1-\xi^2}\alpha\_i t) + \beta\sin(\sqrt{1-\xi^2}\alpha\_i t)\right)e^{-\xi\alpha\_i t} \tag{47}$$

In order to determine α and β, we integrate Equation (46) between 0− and 0+, once and twice:

$$\left(G'(0^+) - G'(0^-)\right) + 2\xi\omega\_i \left(G(0^+) - G(0^-)\right) = 1$$

$$\left(G(0^+) - G(0^-)\right) = 0\qquad(48)$$

These two equations and the boundary conditions impose *G(*0+*)* = 0 and *G (*0+*)* = 1. The solution thus finally reads:

$$\text{For } t > 0 \colon G(t) = \frac{1}{\omega\_i \sqrt{1 - \xi^2}} e^{-\xi \omega\_i t} \sin \left(\omega\_i \sqrt{1 - \xi^2} t\right)$$

$$\text{For } t < 0 \colon G(t) = 0 \tag{49}$$

The mode amplitude *qi(t)* of the dynamic component associated with the event occurring at time *t* = 0 simply follows, for *t >* 0:

$$q\_i(t) = -\left(\int\_{\varepsilon}^{1} s^2 \bar{U}(s) V\_i(s) ds\right) G(t) \tag{50}$$

#### **SURGICAL PREPARATION AND WHISKER STIMULATION FOR** *In vivo* **EXPERIMENTS**

Experiments were conducted in conformity with French (JO 2001-464) and European legislation (86/609/CEE) on animal experimentation. An adult male Wistar albino rat (318 g) was anesthetized with urethane (1.5 g/kg, i.p.). Body temperature was maintained at 37°C. The scalp was retracted after a subcutaneous injection of a local anesthetic (Lidocaine 1%). The skull was then cemented to a metal bar. A 1 mm thick metal bar was moved by a linear stepper motor (Linear Motor MLL302, Systro Gmbh) at a speed *V*bar = 60 mm/s. The bar contacted the whisker near its tip. The whiskers were trimmed at the end of the experiment, and their diameters were measured at different positions, from base to tip. This allowed us to determine the length *L* of the ideal (non-truncated) cone.

#### **HIGH-SPEED VIDEOGRAPHY OF WHISKER DEFLECTIONS**

Whisker movements were recorded at a frame rate of 2.5 kHz with a high-speed camera (Photron Fastcam SA3/105 mm f-2.8 DG Macro Sigma; pixel resolution 55 μm) mounted vertically above the animal. Whiskers were illuminated from below using a backlight (SSLUB, Phlox and PP520, Gardasoft). The camera produces bird's-eye view movies of the whiskers. We choose to study whiskers β and C1 whose curvature plane are essentially horizontal and can thus be entirely imaged with the present optical configuration (**Figure 7A**).

#### **WHISKER TRACKING**

Detecting the shock-induced deflection wave, whose typical amplitude is of the order of a few tens of microns, required subpixel determination of the whisker centerline profile in each frame. The whisker deflected by the object was tracked using a semi-automated algorithm written in Python. We subtracted the average background image from each frame in the movie. Whisker centerline profiles were then extracted: we scanned all columns of each movie frame in order to extract the pixel of maximal intensity defining the approximate position of the whisker. The precise position of the whisker centerline in each column of the frame was then defined as the barycenter (each pixel is

#### **REFERENCES**

encoding of texture in the whisker sensory pathway. *PLoS. Biol.* 3:e17. doi: 10.1371/journal.pbio.0030017

weighted by its intensity) of the pixels surrounding the pixel of maximal intensity (center position ±2 pixel) . In order to extract the shock-induced deflection sequence (**Figure 7B**), the whisker profile determined before the shock was subtracted to each postshock profile. The deflection sequence was then smoothed using a five pixel wide sliding window. This procedure yields a spatial resolution of ≈1 μ on the shock-induced whisker deflection.

#### **ANALYSIS OF EXPERIMENTAL DATA**

Our experimental configuration slightly differed from that considered in the model. In the latter, the whisker rotated at constant speed across a fixed obstacle. In the experiment, the whisker was fixed, and the obstacle was moved linearly at constant speed *V*bar. For this shock configuration, it is easy to show that Equation (11) becomes:

$$U(s, t > 0) = V\_{\text{bar}} t \bar{U}(s) - V\_{\text{bar}} \sum\_{i} V\_{i}(s) \left( \int\_{e}^{1} s^{2} \bar{U}(s) V\_{i}(s) ds \right)$$

$$\times \frac{e^{-\xi \omega\_{l} t} \sin \left( \sqrt{1 - \xi^{2}} \omega\_{l} t \right)}{\sqrt{1 - \xi^{2}} \omega\_{l}} \tag{51}$$

The back-illumination geometry combined with the bird's-eye optics results in a shadow region in the vicinity of the object that prevented the tracking of the whisker down to the whisker/object contact point. Similarly, the fur hampered a proper observation of the whisker near the skin. As a consequence, the positions *s* = ε and *s* = 1 could not be determined from the images, and were thus left as free parameters. The wave speed was determined by linear fitting the position of the first deflection minimum as a function of the time elapsed since the first contact (**Figure 7C**). The four first resonant modes were used to fit the data. Fitting was performed with Matlab's fminsearch function. Correlation between experimental and theoretical data was assessed with Pearson's correlation test on the first 3 and 5 ms for C1 and β whiskers, respectively.

#### **ACKNOWLEDGMENTS**

We thank D. Golomb and S. Neukirch for stimulating discussions, V. Ego-Stengel and L. Bourdieu for helpful comments on the manuscript and E. Coudert for drawing. This work was supported by CNRS, ANR (project TRANSTACT 0294) and EU 7th Framework Program under grant N◦ 269921 (BrainScaleS).

Characterization of statical properties of rat's whisker system. *IEEE Sens. J.* 12, 340–349.

F. L. (2002). Similarities and differences in the innervation of mystacial vibrissal follicle-sinus complexes in the rat and cat: a confocal microscopic study. *J. Comp. Neurol.* 449, 103–119.

and performance. *J. Neurosci.* 26, 8451–8464.

*Nanoscience and Technology*. New York, NY: Springer.

vibrissal touch. *J. Neurophysiol.* 95, 791–802.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 28 March 2012; accepted: 18 October 2012; published online: 06 November 2012.*

*Citation: Boubenec Y, Shulz DE and Debrégeas G (2012) Whisker encoding of mechanical events during active tactile exploration. Front. Behav. Neurosci. 6:74. doi: 10.3389/fnbeh.2012.00074*

*Copyright © 2012 Boubenec, Shulz and Debrégeas. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

### Sensory prediction on a whiskered robot: a tactile analogy to "optical flow"

#### *Christopher L. Schroeder 1 and Mitra J. Z. Hartmann1,2\**

*1 Department of Biomedical Engineering, Northwestern University, Evanston, IL, USA 2 Department of Mechanical Engineering, Northwestern University, Evanston, IL, USA*

#### *Edited by:*

*Robyn Grant, University of Sheffield, UK*

#### *Reviewed by:*

*Quan Zou, University of Nevada, USA Gerald Loeb, University of Southern California, USA*

#### *\*Correspondence:*

*Mitra J. Z. Hartmann, Mechanical Engineering Department, Biomedical Engineering Department, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA. e-mail: m-hartmann@ northwestern.edu*

When an animal moves an array of sensors (e.g., the hand, the eye) through the environment, spatial and temporal gradients of sensory data are related by the velocity of the moving sensory array. In vision, the relationship between spatial and temporal brightness gradients is quantified in the "optical flow" equation. In the present work, we suggest an analog to optical flow for the rodent vibrissal (whisker) array, in which the perceptual intensity that "flows" over the array is bending moment. Changes in bending moment are directly related to radial object distance, defined as the distance between the base of a whisker and the point of contact with the object. Using both simulations and a 1×5 array (row) of artificial whiskers, we demonstrate that local object curvature can be estimated based on differences in radial distance across the array. We then develop two algorithms, both based on tactile flow, to predict the future contact points that will be obtained as the whisker array translates along the object. The translation of the robotic whisker array represents the rat's head velocity. The first algorithm uses a calculation of the local object slope, while the second uses a calculation of the local object curvature. Both algorithms successfully predict future contact points for simple surfaces. The algorithm based on curvature was found to more accurately predict future contact points as surfaces became more irregular. We quantify the inter-related effects of whisker spacing and the object's spatial frequencies, and examine the issues that arise in the presence of real-world noise, friction, and slip.

**Keywords: rat, vibrissae, prediction, tactile flow**

#### **INTRODUCTION**

As an animal moves through the environment, the spatial and temporal gradients of sensory data it acquires are related through the velocity of its moving sensory surfaces. This relationship is represented by the "complete derivative" (Munson et al., 2009), and represents a mathematically inviolate description of information flow over moving sensory surfaces.

In the field of visual neuroscience, the complete derivative has been termed the "optical flow" equation. The optical flow equation relates spatial and temporal intensity (brightness) gradients to the velocity of the animal (Barron et al., 1994; Beauchemin and Barron, 1995; Horn and Schunck, 2003). Recent papers have formalized the idea of "tactile flow" to describe the flow of strain energy density across the hand (Bicchi et al., 2008; Scilingo et al., 2008), and bending moment (torque) at the whisker base across the rodent vibrissal (whisker) array (Gopal and Hartmann, 2007; Hartmann, 2009).

Notably, the optic and tactile flow equations are typically used by assuming that the animal makes use of spatial and temporal gradients of brightness (vision) or strain energy (tactile) to compute the velocity of its sensor array.

We recently proposed a complementary scheme: namely, if the animal already knows its own velocity, then it can use the complete derivative to predict future sensory data (Gopal and Hartmann, 2007; Hartmann, 2009). Computing the complete derivative at multiple spatial scales would allow the animal to predict the stimulus that it will measure in the next sensory instant and provide a mechanism to distinguish between externallygenerated and self-generated motion.

In the present work, we used the rat vibrissal system as a model to examine the plausibility of using the complete derivative to predict upcoming sensory data. Rats actively brush and tap their whiskers (vibrissae) against objects to tactually explore the environment. During detailed exploration of objects, rats move their vibrissae rhythmically, between 5 and 25 Hz (Welker, 1964; Berg and Kleinfeld, 2003). During navigation behaviors, however, rats often hold their whiskers out in a relatively static position as they follow along a wall or tunnel.

The present work was designed to investigate tactile flow across the vibrissal array in this type of translational navigation behavior. We used both simulations and a 1 by 5 array of artificial vibrissae (hardware) to investigate the plausibility of using the complete derivative to predict upcoming whisker-object contact points.

#### **METHODS**

Our initial investigations were performed in simulation, and we then investigated the real-world issues that arise with implementation on a 1 × 5 array of robotic whiskers.

#### **ALGORITHMS FOR THE DETERMINATION OF RADIAL DISTANCE, SLOPE, AND CURVATURE**

#### *Radial distance determination using translations instead of rotations*

Radial object distance is defined as the Euclidean distance between the base of a vibrissa and the point at which it makes contact with an object. Previous work has shown that an object's contour can be extracted by continuous rotation of a vibrissa against an object (Kaneko et al., 1998; Solomon and Hartmann, 2010). As the vibrissa "sweeps" against the object through a rotation, measurement of the bending moment (torque) at the vibrissa base permits the radial object distance to be continuously computed, and the object contour thereby inferred. This technique was empirically validated with sweeps of a vibrissa past three differently shaped objects (Solomon and Hartmann, 2010). The rotational sweep models vibrissal motion during typical "whisking" behavior, in which the rat rotates its vibrissae at their base.

During wall-following behavior, however, rats do not always whisk, but sometimes keep their vibrissae protracted against the wall. This behavior is better modeled as a relative translation between vibrissa and object, rather than as a rotation. Although Solomon and Hartmann (2010) developed equations to describe the translational sweep of a vibrissa past an object, the equations were not validated in hardware.

In the present study, we experimentally validate the translation equations developed in Solomon and Hartmann (2010) and use them to extract object contours. Equations 1–3 below provide only a brief overview of the differences between radial distance extractions during rotation versus translation. A more complete description of radial distance extraction based on rotation is provided in Solomon and Hartmann (2010).

As shown in previous work (Kaneko et al., 1998; Solomon and Hartmann, 2006, 2008, 2010, 2011; Birdwell et al., 2007) the radial distance *r*0 to the initial contact point on an object can be calculated using

$$r\_0 = k \frac{\phi\_0}{M\_0} \tag{1}$$

where *k* = 3*EI*, *E* is the elastic modulus of the vibrissa, *I* is the area moment of inertia, φ0 is a small pushing angle beyond initial contact (typically about 3◦), and *M*0 is the bending moment at the vibrissa base.

Once the radial distance *r*0 to the initial contact point is calculated, the whisker undergoes either a small rotation (*d*φ) or translation (*dL*). Radial distance at the current time step (*ri*) can then be calculated based on estimates of radial distance at previous time steps.

The procedures for calculating *ri* after a translation or rotation are quite similar, but they differ in the calculation of the magnitude of the vector δ*i*−1, perpendicular to the longitudinal axis of the vibrissa at the contact point. With the x-axis defined to be parallel to the vibrissa at *t* = *t*0, Solomon and Hartmann (2010) define δ*i*−1 for rotation as:

$$\vec{\delta}\_{i-1} = -r\_{i-1} \cdot d\phi \cdot \begin{bmatrix} \sin \phi\_{i-1} \\ \cos \phi\_{i-1} \end{bmatrix} . \tag{2}$$

where *ri*−1 is the radial distance at the previous time step (before rotation), *d*φ is the incremental rotation, and φ*i*−1 is the angle between the vibrissa base at time *t*0 and the base at time *ti*−1 (after rotation).

In contrast, δ*i*−1 for translation is defined by:

$$
\overline{8}\_{i-1} = \begin{bmatrix} 0 \\ dL \end{bmatrix}. \tag{3}
$$

where *dL* is the incremental linear movement between time steps. It is important to note that the base of each vibrissa is fixed perpendicular to the direction of linear movement for this definition of δ*i*−1.

The value of δ*i*−1 is then used to calculate the next radial distance *ri*, regardless of whether the whisker was translated or rotated. For a detailed explanation of this calculation, please refer to **Figure 4** of Solomon and Hartmann (2010) and Equations 6–11 also in Solomon and Hartmann (2010).

#### *Computing object slope and curvature from radial distances*

The local slope (μ) and the local curvature (κ) of an object can be calculated based on radial distance measurements. To calculate either μ or κ, we first define our coordinate system as shown in **Figure 1**. The distance *dsi* is the distance the base-point of the whisker has translated in timestep *i*, and the direction of *ds* is always coincident with the rat's velocity. Radial distance estimates at multiple time points, the angle *d*θ, and the distance *ds* can then

**FIGURE 1 | Coordinate system for the determination of slope and curvature.** At each time step, radial distance is calculated relative to the object. The rat's velocity is indicated by *V*rat, and the arc *s* is always parallel to *V*rat. The distance *dsi* is the distance in the direction of *s* between discrete time measurements *ti* and *ti*−1. *dsi*−1 is the distance in the direction of *s* between discrete time measurements *ti*−1 and *ti*−2. *Ri* is the radial distance measure at the *i th* time step, *d*θ*i* is an angle measure of the change in direction of *V*rat, and μ*i* represents the local slope of the object.

be used to find the slope μ*i* at each time step:

$$\mu\_i = \frac{R\_i - R\_{i-1} + ds\_i (d\theta\_{i-1} + d\theta\_i)}{ds\_i - R\_i (d\theta\_{i-1} + d\theta\_i) + R\_{i-1} d\theta\_{i-1}} \tag{4}$$

Slope μi is defined by the differences between radial distance measures (*Ri* and *Ri*−1) divided by the distance traveled (*dsi*). Each *d*θ term accounts for changes in the rat velocity vector between radial distance measurements.

Calculation of curvature is similar, but requires us to first calculate the slope μ*i*−1 obtained in timestep *i* − 1:

$$\mu\_{i-1} = \frac{R\_{i-1} - R\_{i-2} + ds\_{i-1}d\theta\_{i-1}}{ds\_{i-1} - R\_{i-1}d\theta\_{i-1}}\tag{5}$$

Curvature is then defined using both μ*i* and μ*i*−1:

$$\kappa\_i = \frac{\left| \frac{\mu\_i - \mu\_{i-1}}{ds\_i} \right|}{\left[ 1 + (\mu\_i)^2 \right]^{\frac{3}{2}}} \tag{6}$$

This equation is simply the discretized version of the definition of curvature measured in 2D Cartesian coordinates. Notably, the commonly-used simplification of Equation 6 that assumes μi is small compared to unity (and therefore reduces the equation to just the numerator) was empirically found to be inaccurate. The small slope assumption does not always hold in our calculations.

When the rat's heading does not change significantly between timesteps, *d*θ*i* and *d*θ*i*−1 are both zero and Equation 4 can be simplified to:

$$
\mu\_i = \frac{R\_i - R\_{i-1}}{ds\_i} \tag{7}
$$

where *R* and *s* still represent differences between radial distance measures (*Ri* and *Ri*−1) over the distance traveled (*dsi*). Equation 5 can be similarly simplified and Equation 6 remains the same.

#### *Calculating slope and curvature based on radial distances measured with multiple whiskers*

So far, we have described the sensory data obtained by a single whisker over time. Equations 4–7 apply equally well, however, to multiple whiskers making simultaneous contact with an object, at different locations on the object. In this case, the subscript *i* in Equations 4–7 should be interpreted to mean the *i* th whisker, instead of the *i* th timestep. **Figure 2** illustrates this idea. **Figure 2A** illustrates that calculation of slope and curvature with a single vibrissa requires memory of the data acquired on the previous timestep. In contrast, **Figure 2B** illustrates that calculation of slope and curvature with multiple vibrissa can be achieved within a single time step. This second method is equivalent to the rat integrating information across vibrissae within the array. In practice, calculation of slope and curvature can occur over both a single whisker (with memory) and multiple whiskers (at different spatial locations) simultaneously.

#### **PREDICTION OF FUTURE POINTS OF CONTACT AND FUTURE CURVATURE OVER TWO DIFFERENT SPATIAL SCALES**