Temporal coherence of mechanical stimulation modulates tactile form perception

The human hand can detect both the form and texture of the contact surface. The detection of skin displacement (sustained stimulus) and change in skin displacement (transient stimulus) are thought to be mediated in different tactile channels; however, tactile form perception may use both channels to reconstruct an uniﬁed tactile imagery of contact surface in the brain. Here, we studied whether the temporal frequency and temporal coherence information of tactile stimuli encoded in sensory neurons could be used to recognize the form of the contact surface. We used a known tactile phenomenon, the ﬁshbone tactile illusion (FTI), as a probe for tactile perception in humans. It uses an archetypal surface geometry that has a smooth central bar and textures (ridges and grooves) in its adjacent areas. By stroking the central bar back and forth with a ﬁngertip, an observer typically perceives an indented surface geometry even though the bar is physically ﬂat. We used a passive high-density pin matrix to extract only the vertical information of the contact surface, while excluding any tangential force caused by rubbing the surface. Participants in the psychological experiment reported the indented surface geometry by tracing the FTI textures with pin matrices of different spatial densities (1.0- and 2.0- mm pin intervals). Participants reported a steep decrease in the relative magnitude of the indented surface geometry when pins in the adjacent areas vibrated synchronously. To investigate possible mechanisms for tactile perception in the FTI, we developed a computational model of sensory neurons to estimate temporal patterns of action potentials from tactile receptive ﬁelds. Our computational data suggest that i) the temporal asynchrony of sensory neuron responses correlated with the relative magnitude of indented geometric perception, and ii) the temporal change of displacements in tactile stimuli correlated with the asynchrony of simulated sensory neuron responses for the ﬁshbone surface patterns. Based on these results, we suggest that both the temporal frequency and asynchrony of activity in sensory neurons could produce tactile form perception.


Introduction
The human hand can detect both the form and texture information of a contact surface 1 ; tactile perception of surface edges are produced when a tactile stimulus is presented onto the finger pad. Human observers tend to actively scan the surface when they examine the roughness of a surface texture 2,3 . This shows that the human perceiver can use both skin displacement and change in skin displacement to detect and recognize the outside world through touch.
Previous studies have shown that tactile form perception and coarse roughness are mediated by slowly adapting type I afferents, which detect skin displacement 4,5 . Meanwhile, roughness perception is also mediated by rapidly adapting afferents 6 , which detect temporal changes in skin displacement (e.g., flutter vibration 7 ) . Detection of skin displacement and changes in skin displacement are thought to be conveyed in different tactile channels; however, tactile form perception may use both types of information. For example, the fishbone tactile illusion (FTI) is a perceptual phenomenon in which a person feels an indented geometry even though the surface being felt is geometrically flat 8 (Fig. 1a). This phenomenon occurs because the flat smooth surface is surrounded by a rough texture. Previous research has shown that it occurs when relative motion between the contact surface and the fingertip exists, suggesting that not only the skin displacement but also the change in skin displacement can produce tactile form perception 9,10 .
In order to explore the contribution of change in skin displacement in tactile form perception, we first conducted a psychological study to test whether we could eliminate the tangential displacement from tactile stimulus using an apparatus that can extract only the vertical displacement of the contact surface at specific spatial resolutions. We found that human participants experienced the low occurrence of indented form perception in the FTI (the attenuation of the tactile illusion) under certain conditions, where the tactile stimulus was temporally coherent. Second, in order to explain the attenuation of the tactile illusion, we developed a sensory neuron model by combining the ion conductance model of mechanosensitive channels 11 and Hodgkin-Huxley model 12 to simulate the response of the receptive field that detects the change in skin displacement. Using the developed mathematical model, we tested whether the temporally coherent tactile stimulus could produce synchronized responses of sensory neurons. By combining the results of the psychological experiment and computational analysis, we propose that human tactile form perception in the fingertip is produced by both skin displacement information and the change of skin displacement encoded in tactile sensory neurons. maximum values at the 0.4 -1.4 mm rib intervals (Fig. 2). The probability decreased with an increasing rib interval above 1.8 mm.
The same data were also analysed using pairwise comparisons 13 . Supplementary Fig 1 was plotted by polling all of data from twelve participants who conducted the comparison task for each pair. The probability of answering "more indented" is equivalent to the relative magnitude of indentation, so in the following we used the relative magnitude of indentation to describe the strength of indented form perception.
The relative magnitude of the indentation in the FTI can be modulated using a pin matrix. This apparatus has multiple straight pins that move passively up and down according to the height of the geometric surface beneath the pin matrix ( Fig.   1b) 8 . The advantage of using a pin matrix is that the perceiver can change the spatial frequency of the tactile stimulus arbitrarily by changing the interval between adjacent pins (Fig. 1c). Moreover, we could eliminate the tangential displacement of the tactile stimulus caused by the mechanical constraint of passive pins. We arranged two different kinds of spacings between the pins (1.0 and 2.0 mm) and two different pin diameters (0.8 and 1.8 mm), without changing the size of contact area to the fingertip. We conducted the same experiment in using a bare finger (see Fig. 1c for configurations). By using a pin matrix with a spacing of 1.0 mm (PM1), the relative depth of the indented geometry was exaggerated. Strikingly, the effect of modulation produced by pin matrices was even larger when a perceiver used a pin matrix with a spacing of 2.0 mm (PM2 and PM3).
We observed two phenomena when using pin matrices. First, the spatial resolution of the pin matrix significantly modulated the geometric tactile perception (Fig. 2, see also Table 1). Second, the probability of determining that a 1.0 mm rib interval was "more indented" decreased from the surface with pin matrices of 2.0 mm spatial interval (PM2 and PM3) to that with a 1.0 mm spatial interval (PM1).
The psychological experiment raised two important questions regarding tactile geometric perception in humans. (i) Is the temporal frequency of sensory neuron responses produced by the influence of fishbone ribs and/or (ii) does temporal asynchrony between sensory neurons affect the relative magnitude of an indented geometric tactile perception in the FTI?
For the first question, we could predict that the narrower the rib interval becomes, the larger the responses of sensory neurons. Therefore, we can assume that the temporal frequency of sensory neuron responses increases with a change in the interval of pin spacing in the pin matrices.
For the second question, using the 2.0 mm pin matrices (PM2 and PM3), all the pins in the rib areas would go up and down simultaneously when a perceiver touches the surface with rib intervals of integers (1.0, 2.0, 3.0, and 4.0-mm). We could assume that the temporal asynchrony of the responses in tactile sensory neurons corresponds to the magnitude of indented geometric perception in the FTI.
To address these questions, we developed a computational model for a tactile sensory neuron based on previous literature.
With this model, we calculated both the temporal frequency and asynchrony of the responses of tactile sensory neurons.

Computational experiment
We developed a simple computational model that captures sensory neuron responses in the peripheral nervous system. Fig. 3 shows the computational model and how the receptive fields respond to mechanical input. This mathematical model has a hierarchical structure composed of multiple stages of touch reception. Briefly, the mechanical deformation produced by touching the pins with a fingertip leads to the activation of mechanosensitive channels (Fig. 3a, b). This assumption is based on the fact that physical quantities related to local membrane stretch can predict the responses of sensory neurons 14,15 . The responses of multiple mechanoreceptors are integrated into one sensory neuron response at a spike initiation site (Fig. 3c,   d). For the simulation, we assumed that the same type of receptive-field mechanoreceptor system is distributed uniformly in the fingertip (Fig. 3e, Materials and Methods). We found that modelling the sensory neuron responses from the behaviours of mechanosensitive ion channels was critical because we assumed that a corresponding timing of mechanical input and the timing of sensory neuron responses was required to elucidate the attenuated indentation of tactile form perception in pin-matrix conditions (Fig. 2).
We performed numerical experiments corresponding to the PM1 condition with different rib intervals. A representative result of the numerical simulation is shown in Fig. 4. As surface texture moved, some of mechanoreceptors responded, leading to responses in sensory neurons (see also Supplementary Fig. 4 and Supplementary Video 1-5).
Spatial event plots revealed that the time for the occurrence of action potentials was more synchronised for the 1.0-mm rib interval than for the 1.

Correlation analysis between psychological data and computational result
In order to explain our psychological results, we conducted a correlation analysis to evaluate the effect of temporal frequency and temporal asynchrony of sensory neuron responses on the psychological data. To evaluate the effect of temporal frequency of sensory neuron responses, we adopted an averaged number of responses from sensory neurons in adjacent areas (fish rib) to gauge the contribution to psychological results. To evaluate the effect of temporal asynchrony, we employed the Shannon entropy for the probability of action potential occurrence. If the Shannon entropy has a small value, the temporal asynchrony of the action potential occurrence is low, implying that multiple tactile receptive fields respond simultaneously. (2.0-mm pin spacing) at different rib intervals. Fig. 5a shows the averaged frequency of action potentials for all seventy-two receptive fields. Using the simulation data, we obtained the maximum frequency of action potentials from sensory neurons was obtained at a rib interval of 0.4 mm, which coincides with the rib interval at which the maximum relative magnitude of the indented geometry was observed (Fig. 2). The same data, however, does not correspond to the sudden dips of the indentation We found that the Shannon entropy corresponded to the psychological result. The Kendall and Spearman's rank coefficients significantly correlated in the Shannon entropy and perceived indentation ranking for all pin matrix conditions (Table 2).
However, these correlated with the frequency of action potentials and perceived indentation only for the PM1 condition. This correspondence between the perceptual phenomenon and the computational model (frequency and asynchrony of sensory neuron responses) may explain the lower probability of indented geometric perception at the 1.0-and 2.0-mm rib intervals, as shown in Fig. 2. In particular, the lower the Shannon entropy values, the lower the relative magnitude to perceive indented forms.

Discussion
We reported that the relative magnitude of depth perception in the fishbone pattern decreased at specific ridge intervals and pin matrix conditions (Fig. 2). The attenuation of the relative magnitude was not found in the bare-finger condition; this suggested that a vertical displacement of spatially discrete tactile stimuli could produce this phenomenon. To explain what we observed in the psychological experiment, we developed a simple computational model of sensory neuron responses by integrating the ion conductance model of mechanosensitive channels in the cell membrane 11 and the Hodgkin-Huxley model 12 . Based on previous studies on anatomical structures of peripheral nerves in the finger pad, we assumed that four mechanoreceptors were connected to a spike initiation zone where the ion currents are integrated ( Supplementary Fig. 2) 16 .
A recent study on computational models of slowly adapting touch receptors has showed that the generator current coming from mechanoreceptors is integrated at a spike initiate zone 17,18 . Lesniak and Marshall et al. proposed that neural spikes are generated at a spike initiation zone, where the generator currents from mechanoreceptors are integrated. In the follow-up study of the same group, Gerling et al. formulated a function for a generator current by combining slowly, rapidly, and ultra-slowly inactivating currents from a mechanoreceptor unit 18 . Their computational model used a leaky integrate-and-fire model to simulate neural dynamics. We adopted the Hodgkin-Huxley model in the computational analysis to generate action potentials at the spike initiation zone in each receptive field without assuming the spike initiation threshold.
We uniformly arranged 72 low-threshold mechanoreceptive fields to calculate the total sensory neuron responses according to the mechanical displacements of the surface pattern presented through the pin matrices (Fig. 3, Supplementary Fig. 5, 6). We calculated the temporal frequency and asynchrony of the computed sensory neuron responses (Fig. 5), and found that temporal asynchrony of sensory neurons could explain our psychological experiment ( Table 2). Even though our mathematical model does not fully depend on neurophysiological data as in a previous study 15 , we could capture the human nature of tactile form perception, at least, of the fishbone tactile illusion.

The origin of temporal asynchrony in computed sensory afferent responses
Temporal asynchrony of computed sensory afferent responses originates from the temporal incoherence of the tactile stimulus.
We calculated the timing of the maximum displacement of the tactile stimulus presented in different pin matrices (Supplementary 5/24 Fig. 6), and calculated the temporal incoherences of the tactile stimulus using the same metric (the Shannon entropy, equations 13 ) as temporal asynchrony of computed sensory afferent responses. We computed the temporal incoherence of the displacement and displacement changes of pins in the pin matrices ( Fig. 6 a and b). We found that the temporal incoherence of displacement changes showed similar trends as the temporal asynchrony of the computed sensory afferent responses (Fig. 6 c). This result implies that the temporal incoherences of displacement changes presented on the skin may explain the relative magnitude of indentation in the FTI.
A previous study reported that human participants showed a lower detection rate for incoherent tactile stimulus after exposure to a vibrotactile stimulus 19 . This implies the existence of a temporal asynchrony (and a motion) detector in humans. Kuroki et al. 2018 19 assumed that the detection of asynchrony in tactile stimuli is used to perceive a moving tactile stimulus. In addition, our data imply that the detection of asynchrony in tactile stimuli is also used in tactile form perception. Based on our psychological results, we propose a model for tactile form perception in humans (Fig. 7). The perceived magnitude of indentation rank correlated with that of the temporal frequency of simulated neurons under the PM1 touch condition, but not under PM2 and PM3 conditions. On the other hand, the magnitude of indentation correlated with the Shannon entropy ( Table   2). Based on these observations, our working model proposes a hypothesis that both the temporal frequency and temporal asynchrony of action potentials in sensory neurons may play a role in tactile form perception.

Computational model predicts that tactile sensory afferent system that can detect the change of displacement may play a role in tactile form perception
Correlation analysis showed that the responses of tactile sensory channels that detect changes in displacements can mediate the information concerning tactile form perception. In general, all four kinds of tactile sensory afferents in the glabrous skin of the human hand respond to the onset of mechanical input, which is equivalent to the change in displacement 20 . For scanned surface patterns, Blake et al. recorded the responses of SA I and RA mechanoreceptive afferents innervating the glabrous skin of the rhesus monkey, and reported that both sensory afferents responded at the edges of scanned raised surfaces 21 . This physiological data suggests that both SA I and RA channels could contribute tactile form perception in fingertip, and this is consistent with our psychological computational data.
A recent study of surface feature detection in first-order tactile neurons has proposed the involvement of the temporal structure of a sensory afferent neuron. Johansson and Birznieks reported that first spikes of human tactile afferents can encode the direction of applied force 22 . Pruszynski and Johansson reported that both the temporal frequency (intensity) and temporal structure of two types of first-order tactile neurons (SA I and RA) could explain the edge orientations of tactile form scanned across a human fingertip 23 . Similarly, neurophysiological data and computational data of modelled sensory neuron responses support the view that millisecond precision spike timing is utilised to perceive the skin vibration 15,24 .
Based on our results, we assume that RA and SA I can affect tactile form perception. Even though a group of studies support that tactile form perception is mediated by SA I 1 , and surface texture is encoded in both SA I and RA I neurons 45 .
Weber et al. 6 also reported that RA could encode both spatial and temporal variations of natural textures.

6/24
Based on the current computational model, we cannot conclude that SA I and/or RA medicate tactile form perception produced by the FTI. However, we propose that RA I can encode the temporal incoherence of tactile stimulus. In addition, the density of receptive fields simulated in our computational model is similar to the neurophysiological arrangement of receptive fields in RA I 16 . This view is worth addressing by conducting neurophysiological recordings 20 or simulating with established computational models 15 by presenting the equivalent tactile stimuli to both SA I and RA channels in the future.

Limitation of this study
In our model, we considered only the displacement of pins as mechanical input, but did not consider the viscoelastic characteristics of the skin. Mechanoreceptors that form tactile receptive fields are embedded in the skin and consist of multiple layers of different mechanical properties. If we could take the viscoelastic effect of the skin into account in our model, we may be able to establish a mathematical model that reconstitutes the observed neural responses to indented 14,25 and/or scanned forms 21 observed in previous neurophysiological recordings. We consider that our computational model could explain, at least, the decreased magnitude of indented perception in the FTI through pin matrices.

Conclusion
In summary, this study shows that frequency of action potentials in sensory neurons alone cannot fully explain the tactile form perception in the FTI. Temporal synchronicity of the tactile stimuli corresponds to the rank of the indented form perception.
We predict that sensory neuron responses (indexed by the Shannon entropy) in addition to the spatial-rate coding theory, may explain the diminished perception in the FTI. Our mathematical model explains the form tactile perception, at least, in the range of small mechanical deformation up to hundreds micro meter based on surface asperity of the texture used in the current study. Several studies that had challenged to model the viscoelastic behaviors of the skin and its mechanical filtering effect, and showed the correspondence between mechanical quantities and neural responses 6,14,18,26,27 . Future study should address whether we can control the relative magnitude of perceived indentation by changing the temporal incoherence of tactile stimuli, and see the correspondence between the psychological result (perceived magnitude of geometric indentation) and temporal incoherence as an index of tactile form cues.  Table 1. Statistical analysis of the effect of rib intervals in each touch conditions (Bare finger, PM1, 2, and 3) (see Fig. 1a). All main effects of touch conditions were significant (p = 0.000). *** and C.I. denote p < 0.001 and confidence interval.

Participants
Twelve participants (age range 19-22 years, six females; all right handed) participated in the psychological experiment. The

Tactile stimulus
We used a fishbone-like pattern 10 of various intervals between adjacent ribs (Fig. 1a). The fishbone surface pattern was described by three parameters: (i) the width of the central strip, (ii) width of a rib, and (iii) the distance between adjacent ribs.
The central strip width was 4 mm, which is within the range of optimal widths for perceiving the contrast in stimulus intensity Mean estimates of probability are plotted, and error bars indicate the 95% confidence interval.

9/24
The surface of fingertip Receptive field  Figure 3. Mathematical model. a A pin of the pin matrix produces skin deformation that leads to b opening a mechanosensitive ion channel in a mechanoreceptor. c The change in receptor potential due to the mechanoreceptor's physiological response is conveyed electrically to a spike initiation site d where four mechanoreceptors are connected. This mechanoreceptor-sensory neuron system works as e a responsive receptive field distributed uniformly in the fingertip. The concept of mechanoreception was developed by combining previously proposed models in a mechanoreceptor 17,28,29 .    We used a passive pin matrix (Fig. 1b, see Nakatani et al., 2005) to present vibrotactile spatial patterns, as well as rubbing the surface directly with a bare finger. This was done in an apparatus to extract height information from the contact surface beneath the pin matrix. We used pins of 0.8 mm or 1.8 mm in diameter; therefore, this pin-matrix could overcome the height up to 0.4 or 0.9 mm. The supporting plate contained holes of 0.9 mm in diameter so that each pin could move up and down smoothly, according to the contact surface height beneath the pin matrix. The advantage of using the pin matrix is that it eliminates most of the tangential displacement of the finger skin.

Procedures
Each participant entered the room ten minutes before the experiment so that the participant's hands would be at room temperature during the experiment. Participants had a comfortable chair and were asked to wear an eye mask during the experiment. The experimenter instructed the participant to touch various tactile stimuli with either a bare finger or the specified pin matrix provided by the experimenter. During the experiment, participants conducted a pair-wise comparison task. In a trial, participants examined two different tactile stimuli in a trial, and were asked to state which one they felt had a deeper indentation (concave

12/24
shape) in its centre in a two alternative forced-choice fashion. Thirty-six possible tactile stimulus pairs were examined in one session, and two sessions were examined in a block. Within a block, the experimenter asked each participant to touch the presented tactile stimuli under one of four touch conditions (bare finger, with PM1, PM2, or PM3). The order of the presented tactile stimulus pair was randomised between blocks for each participant. The order of the touch conditions was also randomised between participants. The participants rested for five min every 30 min to avoid sensory fatigue. The experiment lasted two hours on a given day and was completed over four days.

Data analysis
We analysed the participants' responses to non-musical and musical sounds separately. We conducted six trials for each touch condition, but we excluded the first two trials of data because the participants took time to get used to the experimental protocol.
The mean probability of answering more indentation was used in the one-way repeated analysis of variance (ANOVA) to observe the within-participant factor (rib interval). When the sphericity was violated, we used Greenhouse-Geisser correction to calculate the number of degrees of freedom of the F distributions. We also calculated the relative magnitude of indentation in the fishbone pattern by employing Thurstone's case V model 13 . Kendall and Spearman rank correlation coefficients were calculated between the psychological data (the rank of mean probability of answering more indented) and numerical results (the rank of temporal frequency or Shannon entropy of sensory neuron responses from receptive fields). Statistical analysis was conducted using IBM SPSS Statistics (Version 25).

Spatial sampling using pin matrix
In the mathematical model, the fishbone pattern used in the experiment (Fig. 1a) is represented as a two-dimensional (2D) region in which the backbone lies along the y axis and the ribs along the x axis. The backbone has a width L w , and the ribs are repeated infinitely in the y direction with rib thickness L 1 and rib-rib spacing L 2 , which makes the spatial rib period L 1 + L 2 .
Then, the height of the fishbone pattern at (x, y), denoted by Φ(x, y), is given by where h 0 is the height of the fishbone pattern and k takes all integer values.
Each pin in the pin matrix moves in the z direction when scanning along the fishbone pattern in the y direction. Here, we choose a coordinate that moves with the pin matrix, which is equivalent to considering the stationary pin matrix on the moving fishbone pattern Φ(x, y −V t) with the scanning speed V . Given a 2D pin coordinate (x p , y p ), the vertical displacement h(x p , y p ) of the pin is calculated by taking into account the spherical shape of the pin tip, as follows: When Φ(x p , y p ) = h 0 , the displacement coincides with the height of the fishbone pattern h(x p , y p ) = h 0 . When Φ(x p , y p ) = 0, the displacement depends on the distance ξ between the pin tip (x p , y p ) and the nearest point of the fishbone pattern, denoted by (x * , y * ). In the present 13/24 fishbone pattern, there are only four possibilities: (x * , y * ) = (x p ± ξ , y p ), (x p , y p ± ξ ). Then, the displacement is determined by where r p = d 2 /2 is the radius of curvature of the pin tip, and ξ * = 2r p h 0 − h 2 0 is the distance within which the pin is in contact with the edge of the printed pattern. Note that we assume r p ≥ h 0 ; if r p < h 0 , then the pin cannot overcome the bump.
Thus, the pin displacement h(x, y) depends on not only the fishbone pattern Φ(x, y), but also the pin radius r p .

Mechanical displacement of the skin
The pin displacement is given to the fingertip and converted to mechanoreceptor deformation. Because the mechanoreceptor lies beneath the epidermis and close to the surface (∼ 0.4 mm), we simplify the detailed elastic effects as follows: The fingertip that is in contact with a pin with displacement h(t) at time t is assumed to have the same amount of vertical displacement h(t), which then causes the mechanoreceptor to deform. We assume a linear relationship between the fingertip displacement h(t) and the mechanoreceptor channel displacement σ (t), namely σ (t) = χh(t), where χ is the conversion ratio of the input stimulus and the channel displacement, chosen appropriately as a free model parameter.

Response of mechanosensitive channels in a mechanoreceptor
Mechanical stimulation of mechanoreceptors can produce receptor currents that pass through the cell membrane of mechanoreceptors. Receptor currents are achieved by the function of mechanosensitive channels located at the membrane of mechanoreceptors 30 . As a mechanosensitive channel, Piezo 2 is widely recognised in mammalian mechanoreception 30 , and has been observed in representative mechanoreceptors [31][32][33] . The piezo 2 channel has a rapidly adapting mechanosensitive current with a time constant of less than 10 ms 30,32 . Inactivation and adaptation properties in mechanosensitive currents have been pointed out as two of the main properties of receptor currents 34 : Inactivation means that the receptor currents diminish when two subsequent stimuli are given, whereas adaptation means that the stimulus threshold, above which a receptor current is triggered, becomes greater when a continuous stimulus is applied.
Prešern et al. proposed a mathematical model of a mechanoreceptor's conductance that reflects these two properties 11 : They considered the mechanoreceptor conductance g(t) in terms of the channel open probability p(t) and the channel inactivation probability q(t), namely, where g max is the maximal conductance of the receptor current. The probabilities p(t) and q(t) are governed by where x is the strength of the mechanical stimulus, and affects both the channel open and inactivation probabilities, while τ p , τ q , x p , x q , k p , k q , and α p are constants. The stimulus applied to the mechanoreceptor is represented by σ (t), which in this model is given by the channel displacement σ (t).
In summary, a stimulus from the fingertip deformation is converted into an input σ (t) = χh(t), which induces an electrical response from a mechanoreceptor as the conductance g(t). Supplementary Fig. 3 shows the typical behavior of the channel state probabilities p(t) and q(t) and the resulting mechanoreceptor conductance g(t) when a rectangular input signal of h(t) is provided.

Integration of mechanoreceptor response into sensory neuron responses
We consider a neural response only in a spike initiation zone 35 . We also assumed that the receptor potentials (or generator potential in 17,35,36 ) from branched sensory neurons can be integrated at a spike initiation site, where signals from the receptors are integrated into the input current that, when sufficiently strong, initiates a nerve pulse as an output signal in the sensory neurons. We do not consider the spatial degree of freedom needed to describe the travelling of pulses. We describe the neural activity using the Hodgkin-Huxley model 12,37 to simulate an action potential in a sensory neuron. In this model, neural dynamics are described by four variables: the action potential of the membrane (v), the activation variables of the potassium channel (m) and sodium channel (n). and the inactivation variable of the sodium channel (h). Their dynamics are governed by where C m is the membrane capacitance; v K , v Na , and v L are the equilibrium potentials of potassium, sodium, and leak channels, respectively; and g K , g Na , and g L are the conductances of these channels. The stimulus current I stim represents the input from the receptors connected to the neuron. Here, we consider the simplified situation in which the mechanoreceptor channels in individual receptors are connected directly to the neuron, and we model the stimulus current I stim as the sum of the receptor 15/24 currents I j from N r receptors: where g j (t) is the conductance of receptor j calculated using (3), and v eq is the equilibrium potential of the membrane. The offset potential v offset is introduced to ensure that the change in the ion conductance from zero is sufficient to initiate an action potential when v = v eq . The functions α m , α n , α h , β m , β n , and β h are given as follows 12 : which were determined through extensive parameter searches (details are shown in Supplementary Fig. 7).

Overview of ensemble of sensory neuron responses
The fingertip is considered as a 10 mm × 20 mm rectangular region, on which N n = 72 neurons are placed in a square lattice with distance L x = L y = 10/6 mm. Each neuron has a circular receptor field with radius r = 1 mm, in which N r = 4 receptors are randomly distributed. See Supplementary Fig. 2.

Numerical experiment
We performed numerical experiments corresponding to PM1, PM2, and PM3. For each experimental set, we prepared twenty scanning speed V is set to 50 mm/s, which corresponds to the psychological experiment conditions. See Supplementary Fig. 4 and Supplementary Movie.

Data analysis
Each numerical experiment yields output data as a set of time series produced by all the neural activities represented by the neural output v(t) from t = 0 to t = T = 800 ms. In each neural activity, a firing event is defined as a moment when v(t) passes a threshold value v * from below. We set v * = 40.0 mV.
We divide the entire simulation time T into n time intervals. In each interval ∆t = T /n, we count the number of firing events for each neuron. Let N i be the sum of the firing events from all neurons in [(i − 1)∆t, i∆t] (i = 1, . . . , n) and the fraction of firing events. Then, the Shannon entropy H, which measures the degree of randomness, is defined as Note that H depends on ∆t, which we choose as ∆t = 4 ms and the mean firing rateF per neuron, averaged over N n neurons, is given bȳ Mechanoreceptors were randomly distributed in a plane, and a receptive field consisting of four mechanoreceptors that were closest to the centre of the receptive field (left). We assumed that the length of separation between the centre of the receptive fields was L x = L y = 1.0 mm in the numerical simulation.  Fig. 3 d) where the sum of ion currents is integrated. If the sum of the ion currents is over the threshold, the spike initiation site could produce an action potential (b bottom). Displacements and conductances of four mechanoreceptors (red, green, blue cyan) are plotted with slight shifts for clarity. An overview of the phase diagram for determining free parameters χ and x q that govern the characteristics of neural responses. a. For three different rib intervals and two different pin conditions, numerical simulations were performed by varying χ and x q . For each parameter set, we simulated the response of a single neuron with four receptors, which were on the same pin and received identical signals. The time profiles of the membrane potentials were calculated and classified into six distinct patterns according to the mean number of responses during one period of input signals, denoted by p: p = 0 (cross), 0 < p < 1 (inverted triangle), p = 1 (circle), 1 < p < 2 (diamond), p = 2 (square), and p > 2 (triangle). Typical profiles are shown in b. We chose the parameter set from the parameter region in which (i) membrane potentials that were in accordance with the pin displacement (circle) were observed in as many cases as possible, and (ii) membrane potentials were preserved even in the case of rib interval = 0.2. Based on this observation, we set the free parameters as χ = 0.07 and x q = 6.0 µm.