Mechano-neurophysiological model of the fingertip for spatiotemporal firing pattern prediction of SA1 mechanoreceptors during embossed letter scanning

Numerical simulation can be used to observe spatiotemporal firing responses of tactile receptors during dynamic tactile exploration, and it provides a more understanding of the mechanism of tactile perception. In this study, we developed an improved mechano-neurophysiological model of the fingertip that employs a realistic fingertip structure and accurate contact mechanics while scanning embossed letters. To confirm the potential of the model, we simulated the spatiotemporal firing patterns of slowly adapting type-1 (SA1) mechanoreceptors while scanning the embossed letter “G” and compared the simulation result with the existing experimental data in neurophysiology. Although the experimental data were reconstructed from a single nerve fiber, the simulation simultaneously observed the responses of multiple SA1 receptors, which resulted in a more obscure “G” spatiotemporal firing pattern than that in the previous experiment. This result supports existing data from another psychophysical experiment that demonstrates that it is harder to recognize embossed letter “G” accurately during letter scanning. This finding suggests that the spatiotemporal firing pattern from multiple SA1 receptors may show an obscure “G” pattern while scanning the embossed letter “G”.


Introduction
Human beings perceive the shape and roughness of an object by grazing their fingers over textured surfaces. This tactile perception is provided by neural signals from the responses of tactile mechanoreceptors during the dynamic exploration. The numerical simulation of this spatiotemporal firing response leads to a more quantitative understanding of the underlying mechanisms of tactile perception. The quantitative understanding will contribute to the development of advanced haptic technology and bionic medical prostheses.
Previous studies [1][2][3][4][5] have developed numerical simulation models to predict the responses of tactile receptors for a quantitative understanding of the mechanisms of tactile perception. In one of the approaches for the simulation, mechanical states (strain Toyota Central R&D Labs., Inc., Aichi, Japan and stress) at the receptors were calculated using a finite element (FE) model of a fingertip and then converted into the firing rate of receptors, based on the insight that the mechanical states are correlated with the receptor responses. 1,2 This approach reasonably enables the prediction of receptor responses. However, it cannot predict the timing of the firing response, that is, the action potential. Hence, Gerling et al. 3 developed a FE fingertip model coupled with a neurophysiological model (i.e. mechano-neurophysiological model) to predict the timing of the firing response. The FE model simulated the mechanical states at the tactile receptors during the finger skin stimulation. The neurophysiological model simulated the transformation of mechanical states into the current and membrane potential of the receptors. This model succeeded in predicting the sphere discrimination of human beings, while that simulated the condition where the finger pad was indented vertically by stimulators. In contrast, some studies simulated the firing responses during the ''dynamic tactile exploration,'' in which the finger pad is stimulated by moving an object vertically and laterally across the skin. Kim et al. 4 and Saal et al. 5 simulated the firing responses of tactile receptors during the movement of a surface with embossed letters across a finger pad. For this simulation, they also developed a mechano-neurophysiological model, and this model enabled the prediction of the firing responses of tactile receptors during the dynamic tactile exploration.
However, previous models 4, 5 have several limitations. Although the mechanical part of the mechanoneurophysiological model in those studies 4,5 was a continuum mechanical model, it did not represent a realistic fingertip structure. Moreover, the point-load distribution was shifted as a substitution for emulating the transient dynamic contact mechanics. These factors may hinder the accurate prediction of spatiotemporal firing responses of tactile receptors that could provide important information about objects.
If the mechano-neurophysiological model includes the mechanical part that treats a realistic fingertip structure and accurate contact mechanics during dynamic tactile exploration, it can predict the more natural spatiotemporal firing responses against the previous models. [1][2][3][4][5] In addition, the model can simultaneously observe the responses from multiple tactile receptors. These advantages will expand the existing knowledges about the spatiotemporal firing responses of tactile receptors. For example, although Johnson and and Phillips 6 have measured the spatiotemporal firing pattern on stimulation of the finger pad by the embossed letter, the pattern was reconstructed from a single nerve fiber. We hypothesized that the improved model would observe different responses from existing data. 6 Therefore, through this study, we aim to develop an improved mechano-neurophysiological model of the fingertip and to simulate the spatiotemporal firing patterns of tactile receptors during dynamic tactile exploration. The mechanical part of the model considered the realistic fingertip structure and accurate contact mechanics during the dynamic tactile exploration; the neurophysiological part of the model used fewer input parameters. We then simulated the spatiotemporal firing patterns of the slowly adapting type 1 (SA1) mechanoreceptors during dynamic tactile exploration (i.e. scanning of an embossed letter), to confirm the potential and advantage of the developed model. The result of the simulation was then compared with experimental data from another study. 6

Mechanical model
The mechanical model was used to simulate contact between the finger skin and an object during dynamic tactile exploration, and mechanical states at tactile receptors were calculated. This study focuses on the SA1 tactile receptor because of its prominent role in spatial perception. 7 A three-dimensional mechanical model representing the fingertip was developed using the finite element method ( Figure 1). The external shape and position of the bone corresponded to those of the model constructed from computed tomography images. 8 The skin of the fingertip model was composed of four tissue layers; the stratum corneum, epidermis, dermis, and subcutaneous tissue, whose thicknesses were 0.17 mm, 9 0.65 mm, 10,11 0.8 mm, 10,11 and 0.7 mm, 11 respectively. The SA1 receptors were located at the space bordering the epidermis and dermis. This three-dimensional model was discretized into finite elements, that is, 10node tetrahedron elements. There were approximately 81,000 elements with 121,000 degrees of freedom. Figure 2 shows the structure of the stimulator, an embossed letter ''G,'' in contact with the finger pad. The width of ''G'' was 6.5 mm, and the height of the embossed part was 0.6 mm. 6 In the literature, 6 the firing pattern of SA1 receptors when the finger pad was stimulated by the embossed letter ''G'' is clearly shown as experimental data.
The material models and properties applied to the FE model were determined with reference to those in the literature. 11,12 A first-order polynomial hyperelastic model (neo-Hookean model) was adopted to the stratum corneum and epidermis, and a second-order polynomial hyperelastic model to the dermis and subcutaneous tissues. The bone and nail were assumed to be linearly elastic. The polynomial models were defined using a function of strain energy density per unit volume, W , as follows, where I 1 and I 2 are the first and second deviatoric strain invariants, respectively; J is the elastic volume ratio; and C ij and D i are the material parameters. C ij is related to the shear modulus and D i is related to the bulk modulus. N =1 (i=1, j=0) is applied to a first-order polynomial model and N =2 to a second-order polynomial model. Table 1 shows each adopted material property. For the stratum corneum and epidermis, the Young's modulus and the Poisson's ratio equivalent to the properties are also shown. The stimulator was defined as a rigid body. Figure 3(a) and (b) demonstrate the validation of the developed mechanical model. Figure 3(a) shows that the mechanical model can account for the relevant reaction force of an actual finger contacting a flat plate. 10 In addition, Figure 3(b) shows that the model can account for skin surface deflection. 13 In the test, the finger pad was indented by the stimulator with a sharp edge (with a 0.025-mm radius of curvature), and the indented displacement was 1.0 mm. 13

Neurophysiological model
The neurophysiological model was divided into two models-the transduction model, which transforms the mechanical states at SA1 receptors obtained by the mechanical model into the current of SA1 receptors; and the membrane potential model, which simulates the membrane potential of SA1 receptors by the transformed current. The transduction model adopted the model suggested by Gerling et al., 3 since their model requires minimal input parameters to simulate the firing responses of SA1 receptors. In addition, the current of SA1 is confirmed to have decayed exponentially when the stimulation intensities become constant. 14 Therefore, the transduction model was modified to represent this decay of SA1 responses from the model suggested by Gerling et al. 3 The membrane potential model adopted the leaky integrate-and-fire (LIF) model, which is one of the most widely used models for simulating the behavior of neural systems. The transduction model (equations (2), (3), and (4)) and the membrane potential model (equation (5)) used in this study are shown below.   Table 1.
where I t ð Þ is the current of an SA1 receptor, U t ð Þ is the strain energy density (SED) at the SA1 receptor calculated by the mechanical model, I peak is the current obtained from equation (2) immediately before dU t ð Þ=dt ł 0, t peak is the time immediately before dU t ð Þ=dt ł 0, t is the time constant, b, K s , and K d are the parameters for transduction, and Dt is the time increment. In addition, V t ð Þ, R, and C are the membrane potential, membrane resistance, and membrane capacity of that SA1 receptor, respectively. Equation (5) is solved using a fourth-order Runge-Kutta numerical integration method.

Parameters for transduction model
The parameters b, K s , and K d for the transduction model are identified in this section. The values of C, R, and V (the threshold for LIF model) were adopted from a previous study. 3 Experimental data of a previous study 15 were used for the parameter identification. The experiment 15 measured the firing responses of an SA1 receptor when the finger pad was indented by cylindrical bars. The firing rate data 15 under four conditions, when the cylinder curvatures were 16, 8, 4, 2 [1/ in], were used for the parameter identification. Figure 4 shows the experimental data. 15 First, b and K s were identified. The current required for generating the firing rate from the experimental data was calculated using equation (5). Since b and K s are parameters related to the steady-state (500 ms to 2400 ms in Figure 4), the average value of the firing rate of the steady-state in Figure 4 was used. Table 2 shows the experimental data used for the calculation, the current obtained by the calculation, and the reconstructed firing rate calculated Table 1. Material properties of the mechanical model. 11,12 Young's modulus, E (MPa) Poisson's ratio, n (-)   Table 2. Experimental data for the computation, the current obtained via the calculation, and the reconstructed firing rate obtained using the current obtained and equation (5). using the obtained current and equation (5). Second, we illustrate the relationship between the obtained current and the SED at the SA1 receptor when the cylindrical bars contact the finger pad, as shown in Figure 5. The SED was calculated using the FE model described above. The results were linearized to derive the approximate expression. In this expression, the intercept was set to b and the slope was set to K s . Finally, the b and K s obtained were applied to determine K d to minimize the squared error of the difference in the peaks of the firing rate between the simulation result and experimental data, with a curvature of 2 [1/in], during the transient state (0 ms to 500 ms in Figure 4). In addition, t was also determined to minimize the squared error of the difference in firing rate between the simulation result and experimental data from 100 ms to 2400 ms with a curvature of 2 [1/in] in Figure 4. Table 3 lists the identified parameters. Figure 6 illustrates the simulation results representing the firing pattern of SA1 similar to the experiment 15 by using the mechanical and neurophysiological models in this study.

Simulation
In the previous experiment 6 that was compared with the present simulation results, the stimulator made contact with the fixed finger and moved across the finger   pad to represent the dynamic tactile exploration. To emulate this experimental condition, the stimulator was displaced by 1.97 mm in the Z-direction, which corresponds to a force of 0.6 N at the contact point. The stimulator was then displaced at a speed of 20 mm/s in the Y-direction across the finger pad, with the bone in a fixed position. The frictional coefficient of the contact surface between the finger pad and the surface was set to 0.3. This value was equivalent to that between the finger pad and a nylon surface, 16 which was used as the stimulus surface in the experiment. 6 Figure 7 shows the locations of the SA1 receptors for simulating the firing responses and the initial positional relationship in the simulation. When the stimulator was at rest (0 mm), the SA1 receptor was located halfway between the embossed letters ''F'' and ''G''. When the stimulator was displaced by 10 mm, the SA1 receptor was located halfway between the letters ''G'' and ''H''. The simulation was performed using a commercial FE code Abaqus (SIMULIA, Dassault Systemes, France).  its local maximum (Figure 8(a) and top of (b)). In addition, the bundle of high firing rate patterns occurred thrice, when the stimulator was displaced from 1.5 mm to 2.3 mm, 5.2 mm to 6.0 mm, and 7.3 mm to 8.5 mm (top of Figure 8(b)). Figure 9(a) illustrates the raster plots of the spatiotemporal firing pattern during the scanning of the embossed letter ''G''. Each color (each row of the raster plots) indicates the firing responses of each SA1 receptor in Figure 7. Figure 9(b) shows the previous experimental data 6 to be compared. When the stimulator was displaced from 0 mm to 5 mm, the SA1 receptors at the center of the finger pad fired first. As the stimulator moved farther, the receptors at the outer side of the finger pad fired gradually after the firing of the center receptor (Figure 9(a)). In contrast, when the stimulator was displaced from 5 mm to 10 mm, the SA1 receptors at the outer side of the finger pad were the first to exhibit high firing rates (Figure 9(a)).

Discussion
In this study, an improved mechano-neurophysiological model of the fingertip was developed. The mechanical part of the model considers the realistic fingertip structure and accurate contact mechanics during dynamic tactile exploration, resulting in the neurophysiological part of the model requiring fewer input parameters than existing models. The developed model then simulated the spatiotemporal firing patterns of SA1 receptors during dynamic tactile exploration, that is, embossed letter scanning, to confirm the potential and advantage of the model.
We succeeded in simulating the characteristics of the spatiotemporal firing pattern in experimental data, 6 which is considered to play an important role in the perception of the shape of the embossed letter ''G''. The bundle of high firing rate pattern from the SA1 receptor at the center of the finger pad was observed thrice during the scanning of the embossed letter ''G'' (top of Figure 8(b)). The tendency of this firing pattern was also found in the experimental data 6 (bottom of Figure  8(b)). This firing pattern is considered to play a role in the perception of the shape of the letter ''G'', that is, the far-left edge, the subsequent horizontal bar ''-'' shape, and the clearance between them. In addition, the firing pattern of Figure 9(a) is similar to the C-shaped part of the letter ''G'', and, hence, the curvature shape and central clearance of ''G'' can be perceived. This tendency of the spatiotemporal firing pattern can also be found in the experimental data 6 (Figure 9(b)).
The present model is helpful to discuss the firing responses while observing the skin cross-section deformation and its mechanical states. Figure 8(a), the top of Figure 8(b) and (c) reveal that the SA1 receptor fired up until the letter ''G'' passed underneath the receptor. This result suggests that the stimulus was already being perceived up until the letter ''G'' passed underneath the receptor, although it might seem reasonable that it was perceived for the first time when ''G'' was located underneath the receptor. In the in vivo experiment, it is difficult to discuss this firing response while observing the skin cross-section deformation and its mechanical states. However, the simulation model developed has no such difficulty. This advantage will further help clarify the tactile perception mechanism in dynamic tactile exploration.
The model developed enables simulating the natural spatiotemporal firing patterns of the receptors during the scanning. The simulation results showed that the firing pattern of SA1 receptors at the outer side of the finger pad was sparse, revealing a rather obscure ''G'' in the overall appearance of the firing pattern ( Figure  9(a)), while the experimental results revealed an obvious ''G'' in the firing pattern (Figure 9(b)). Under the simulation conditions in this study, the center of the finger was always on the center-line of the stimulator (see Figure 7). In addition, since the finger pad has a curvature in the finger-width direction (X-direction), the contact force between the pad and the embossed letter ''G'' was at its maximum in the center of the finger pad, and it decreased toward the outer side of the finger pad. Therefore, the SA1 receptors at the outer side showed sparser firing patterns than those at the center of the finger pad. Thus, the overall appearance of the firing pattern revealed an obscure ''G''. In contrast, each row of the raster plot in the experimental data 6 (Figure 9(b)) exhibits the firing pattern from the same single SA1 receptor, because these firing responses were recorded in each trial while repeatedly shifting the stimulator by 200 mm in the X-direction. Therefore, without being affected by the contact force, such as the simulation condition, the overall appearance of the firing pattern showed an obvious ''G'' in the experimental data, 6 unlike the simulation result. The spatial firing data of the SA1 receptor 6 was constructed from a single fiber by repeatedly scanning the letters while shifting the stimulator, possibly because it was difficult to accurately identify all the target SA1 afferents from an enormous number of nerve fibers. In the present simulation model, there is no such difficulty, and it is possible to observe the natural spatiotemporal firing patterns of SA1 receptors unlike in the experiment. 6 This advantage will further clarify the mechanism of tactile perception during dynamic tactile exploration.
It is thought that the appearance of the firing pattern in the present simulation result can be an obscure ''G'' firing pattern. Although it was surprising to note the high resolution with which the firing pattern of the experimental data 6 was able to map ''G'', whether human beings, for example, can actually perceive the detailed shapes of ''G'' with a width of 6.5 mm remains to be seen. Indeed, the experimental data 17 of letter recognition by tactile exploration showed that subjects were not able to recognize ''G'' accurately. It is assumed that the obscure ''G'' firing pattern in the present simulation (Figure 9(a)) may explain this incorrect perception. In addition, according to Figure 9(a) and (b), when present above the right-hand gap of ''G'', the SA1 receptors in the simulation result fired strongly. In the experimental data, 6 however, the receptors from the same position fired weakly (see arrows I in Figure 9(a) and (b)). This indicates that the firing pattern of the simulation result cannot perceive the gap, while that of the experimental data 6 can. Considering that the twopoint discrimination ability in the finger pad for human beings is approximately 2 mm, 18 human beings might not be able to perceive the gap of ''G'' because the width of the gap is at most 1 mm in the X-direction.
The present model has some advantages over the previous ones. [4][5][6] Kim et al. 4 and Saal et al. 5 simulated the firing responses of tactile receptors when the stimulator along with the embossed letters came in contact with the finger pad, by a combination of a mechanical model and neurophysiological model. However, their mechanical model was a simplified continuum model, and it did not represent a realistic fingertip structure. Furthermore, although these works 4,5 seem to treat the transient dynamic contact during embossed letters scanning, the contact was substituted by shifting the point-load distribution on the skin surface. Although Johnson and Phillips 6 measured the spatiotemporal firing pattern on stimulation of the finger pad by the embossed letter, the pattern was reconstructed from a single nerve fiber. The present study, however, treated the realistic fingertip structure and contact mechanics during the embossed letter scanning by using a threedimensional FE model. In addition, the present model was able to simultaneously observe the spatiotemporal firing pattern from multiple SA1 receptors. We believe that these advantages of our model will provide a different approach to expand the existing knowledge of the spatiotemporal firing responses of tactile receptors. We hypothesized that our model would observe different responses from existing data. 6 In this study, we observed a more obscure ''G'' spatiotemporal firing pattern than that in experimental data 6 from previous studies. This result, as mentioned above, supports the data from the earlier psychophysical experiment 17 that demonstrates that it is harder to recognize embossed letter ''G'' accurately during letter scanning. This finding suggests that the spatiotemporal firing pattern from multiple SA1 receptors may show an obscure ''G'' firing pattern while scanning the embossed letter ''G''. Since the previous models 4,5 and measurements 6 did not offer such a new suggestion, these are original and beneficial ideas introduced in this study.
In addition, the present model can accurately treat contact mechanics, including frictional behavior, unlike previous models. 4,5 Doi and Fujimoto 19 showed that the frictional coefficient between a finger pad and a textured surface affects the recognition of the textured pattern. The present model, which can treat frictional behavior, will contribute to a more quantitative understanding of the underlying mechanisms of that tactile perception. This is one of our main future works. The present model has several limitations, which should be investigated in the future. First, we did not treat the human skin as a viscoelastic material. It is known that viscoelastic properties give the finger skin its dynamic behavior. 12 However, the relationship between these properties and the firing of tactile receptors is not clear. Previous studies 20,21 have shown that the firing decay of SA1 receptors (e.g. Figure 6) may be related to stress relaxation due to the viscoelastic properties of the skin. In contrast, the decay was able to be confirmed by stimulating an SA1 receptor alone. 14 In addition, it has been reported that stress relaxation alone could not account for the decay. 22 For these reasons, the present study did not consider the viscoelastic property aimed at representing the decay. In fact, the actual firing responses, including the decay, can be reconstructed without considering the skin as a viscoelastic material, by treating the decay characteristic of the SA1 receptor in this study ( Figure 6). However, since the relationship between the firing responses of the tactile receptor and the skin's viscoelastic property is not clear, it should be investigated in the future. The second limitation is that the present model did not account for the frequency response characteristics of the tactile receptor, which would be required to simulate the receptor responses during flutter and vibration in the skin. Thirdly, the present fingertip model did not include fingerprints, to reduce the computational cost. However, the fingerprint structure may be important in cases where simulating the contact between fingerprints and the asperity of an object is needed. To reduce the computational cost while studying the influence of fingerprints on the firing of tactile receptors, a twodimensional fingertip model should be adopted. In addition, adjusting or modifying the present model may need more comparison with other experimental data. For example, data on widely deformed finger pad and responses of tactile receptors under various stimulations. In future works, improving these aspects may give better simulation results and findings.

Conclusion
In this study, we developed an improved mechanoneurophysiological model of the fingertip that employs a realistic fingertip structure and accurate contact mechanics while scanning embossed letters. To confirm the potential of the model, we simulated the spatiotemporal firing patterns of slowly adapting type-1 (SA1) mechanoreceptors while scanning the embossed letter ''G'' and compared the simulation result with existing experimental data. 6 Although the experimental data 6 was reconstructed from a single nerve fiber, the simulation simultaneously observed the responses of multiple SA1 receptors. As a result, we observed a more obscure ''G'' spatiotemporal firing pattern than that in the previous experiment. 6 This result supports the existing data from another psychophysical experiment 17 that demonstrates that it is harder to recognize the embossed letter ''G'' accurately during scanning. This finding suggests that the spatiotemporal firing pattern of multiple SA1 receptors may show an obscure ''G'' firing pattern while scanning the embossed letter ''G''.