Effect of 3D microstructure of dermal papillae on SED concentration at a mechanoreceptor location

The feeling of touch is an essential human sensation. Four types of mechanoreceptors (i.e., FA-I, SA-I, FA-II, and SA-II) in human skin signalize physical properties, such as shape, size, and texture, of an object that is touched and transmit the signal to the brain. Previous studies attempted to investigate the mechanical properties of skin microstructure and their effect on mechanoreceptors by using finite element modeling. However, very few studies have focused on the three-dimensional microstructure of dermal papillae, and this is related to that of FA-I receptors. A gap exists between conventional 2D models of dermal papillae and the natural configuration, which corresponds to a complex and uneven structure with depth. In this study, the three-dimensional microstructure of dermal papillae is modeled, and the differences between two-dimensional and three-dimensional aspects of dermal papillae on the strain energy density at receptor positions are examined. The three-dimensional microstructure has a focalizing effect and a localizing effect. Results also reveal the potential usefulness of these effects for tactile sensor design, and this may improve edge discrimination.


Introduction
There are four types of mechanoreceptors, namely Pacinian corpuscles, Ruffini endings, Merkel cells, and Meissner corpuscles, which work as haptic sensors for humans. Pacinian corpuscles, Ruffini endings, Merkel cells, and Meissner corpuscles are the end organs of rapidly adapting type-II mechanoreceptors (FA-II) and slowly adapting type-II mechanoreceptors (SA-II) that correspond to SA-I mechanoreceptors and FA-I mechanoreceptors, respectively. They are positioned at different depths of the skin (Fig 1a). Specifically, SA-I and FA-I receptors share the intermediate region between the epidermis and dermis. As shown in Fig 1,  to fingerprints, and dermal papillae that correspond to an aspect parallel to fingerprints. The anatomical components of skin and geometry of ridges are described in [1] and [2].
There is a paucity of research on the mechanotransduction mechanism by which a tactile stimulus is transformed into a neural signal. Finite-element modeling (FEM) is a reliable method for investigating the mechanical properties of skin components and the manner in which they affect mechanoreceptors. Various modeling techniques (from macrostructures to microstructures) were proposed in the biomechanics literature.
In term of macrostructures, a significant model corresponds to the three-dimensional (3D) model of a whole fingertip that was developed based on real human and monkey fingertip geometry by Dandekar et al. [4]. The strain energy density (SED) calculated at the SA-I receptor location of this model was matched with neurophysiological data reported by Phillips and Johnson [5]. A model of the same scale was developed by Gerling et al. [6] with extensive transduction and neural sub-models that convert the SED into neural spikes. The model matched surface deflection in human experiments conducted by Srinivasan [7] and exhibited a good correlation with single afferent responses obtained by Johnson [5] as well as psychophysical results derived by Goodwin [8]. However, fingerprints and other microstructures of skin, such as intermediate ridges, dermal papillae, and fibril structures, were not included in the forementioned macro-structured 3D models.
Results indicate that the microstructures strongly affect the deformation of skin. Hence, an important challenge involves the modeling of skin microstructure to investigate its effect on mechanoreceptors. Maeno et al. [9] developed a two-dimensional (2D) plane strain model based on cross sections of a human fingertip as opposed to using a whole-fingertip model. Their "intermediate ridges" and fingerprints (also called "papillary ridges") are based on a 2D model developed by Srinivasan [10]. Maeno [9] barely exhibited any effect of the intermediate ridge with respect to the stress concentration at the positions of the FA-I receptors. A potential explanation corresponds to the lack of anatomically detailed dermal papillae in the model.

The microstructure of skin in three dimensions
A gap exists between the natural geometry and conventional FE models because 2D models assume identical geometry along the depth direction (Fig 1b). In contrast, Cauna's observation [2] [3] revealed the convex-concave structure of dermal papillae. Additionally, FA-I receptors are located at common apexes of both dermal papillae and the intermediate ridge (Fig 1c), and thus, a 3D model with both intermediate ridges and dermal papillae (3D microstructure) is required to investigate their mechanotransduction.
Vodlak [12] developed a full-scale fingertip with a single representative volume of dermal papillae (containing one FA-I receptor) that focused on the anatomical detail of the FA-I receptor and its transduction method. This model is potentially useful in investigating 3D microstructure although modeling the 3D microstructure at the macro-scale is complex and computationally expensive.
This paper presents a small-scale model that includes the 3D microstructure of dermal papillae. The proposed model is validated with a standard line-load method. This is followed by comparing the 3D model with a model with a 2D microstructure to demonstrate how the 3D microstructure affects the SED concentration at mechanoreceptor positions, i.e., at the

FE models
The model was considered as a cut-away cube of a fingertip and contained six fingerprints along with 12 intermediate ridges with exterior measurements as shown in Fig 2. The configuration of the limiting ridges is excluded. The model is relatively small, and the fingerprints are arranged evenly. The depth of skin calculated from the surface of the epidermis to the center of bone typically approximately corresponds to 11.55 mm while the convex of the surface is less than 0.121 mm. The difference ratio approximately corresponds to 1.04%, and thus it is reasonable to assume that the surface is flat in this case. This assumption is also supported by anatomical observations [13]. A thin layer of subcutaneous tissue is used as a buffer, and its thickness is determined through a line-load validation as described in the next section. This modeling method reduces the calculation cost while allowing the model to have a sufficient thickness for deformation in three dimensions.
Two models were created for the simulation, namely 2D and 3D ridged models. The 2D ridged model includes intermediate ridges, which are concave-convex structures with trapezoid apexes. The dermal papillae are considered to possess even geometry (Fig 3a). The 3D ridged model includes the same intermediate ridges and the concave-convex dermal papillae along the depth (144 tips in total). The dermal papillae are simplified into trapezoid apexes and bases with the same width as those of the intermediate ridge. The length of dermal papillae is set at 0.15 mm, and this is equivalent to the average length of an FA-I receptor in humans (Fig 3b).
The models are multi-layered and include the epidermis, dermis, and subcutaneous tissue (Fig 2) with elastic properties as adapted from previous studies [9]. The Young's modulus corresponds to 0.136 MPa for the epidermis, 0.08 MPa for the dermis, and 0.034 MPa for the subcutaneous tissue. The Poisson's ratio corresponds to 0.48 for each layer [9] [11] [14]. In this study, we focus on the mechanical response at FA-I location under small displacement loads, and thus the model is assumed as linearly viscoelastic. Sophisticated hyperelastic models [15] are excluded as [16] [17] suggested that a non-linearity model is necessary for large displacement loads. As suggested by [17], the outer most layer is linear elastic and the others are viscoelastic. The total stress is as follows: where t denotes time, g(t) denotes the stress relaxation function, and σ 0 (t) denotes the instantaneous stress. We define g(t) as a two-term Prony series.
where g i and τ i denote the stress relaxation parameters that are obtained by fitting the net force response of the model to the experimental data as adopted from [17]. The fitting program is developed and implemented in Python with the support of the Scipy (v0. 19.0) library. The finite element software ANSYS release 16.0 (ANSYS Co.) is used to mesh and analyze the responses. The mesh uses 20-node solid brick elements for surface-to-surface contact. The coefficient of friction between the models and indenters is assumed as zero due to small indentations. The numbers of nodes and elements in 3D ridged models correspond to 267,156 and 88,482, respectively. The numbers of nodes and elements in 2D ridged models correspond to 114,879 and 34,131, respectively. The differences in numbers are due to the complexity of

Viscous parameters calibration for dynamic experiments
We simulate the indentation experiments with a cylindrical indenter as performed in [17]. The indenter is circular with a diameter of 0.5 mm. The indenter is ramped up to 0.2 mm (0.5 mm/ s) and held for 3 s and then completely retracted. We calculate the contact force at the surface of the fingerprints as a function of time.

Experiments
Model validation. The displacement at the surface of model (in response to 50 μm line load indenter) is compared to experimental data from [7] to validate the mechanical response of the model. The line-load indenter results in a displacement of 1 mm that is perpendicular from the surface of the models (at the rate of 1 mm/s) and is held for 2 s. The surface deflection is computed at t = 2 s to allow the model to reach its stable condition. Three restraint conditions are tested to determine the boundary condition that best approximates the observed displacement behavior, namely full-restrained (bottom and surrounds), bottom-faces-restrained, and bottom-edges-restrained conditions (Fig 5a). The models used for the three boundary conditions are denoted as 1x020s, 1x020f, and 1x020e.
Another boundary test focuses on the thickness of the subcutaneous layer. The thickness of the subcutaneous layer varies in range of 0.2 mm, 0.4 mm, 0.6 mm, 0.7 mm, 0.8 mm, and 1 mm and the six conditions are as denoted 1x020e, 1x040e, 1x060e, 1x070e, 1x080e, and 1x100e, respectively). All the models were bottom-edges-restrained.
Ramp-and-hold experiments. In order to investigate the effect of a 3D microstructure, a 2.0×2.0 mm 2 bar provides displacement to a depth of 0.5 mm from the model surface at the rate of 5 mm/s and is held for 2 s. The displayed time point is the first time that the indenter The strain rate is constant during the ramp phase, and thus the stress-strain curves are non-linear. For the purpose of convenience, the concept of equivalent strain energy density is used [18]. The convex-concave shape of dermal papillae is considered, and the stress-strain behavior at their tips is similar to the elastic-plastic behavior at a notch tip. The Glinka's concept relates equivalent strain energy density for an elastic fictitious material as well the real elastic-plastic material during a low strain. Hence, the SED is estimated as follows where σ e denotes equivalent stress and e denotes equivalent strain. The equivalent stress and equivalent strain is related to the principal stresses and strains as given by the following equations where ν denotes the material Poisson's ratio.
The SED distribution at 10 specific x-axis-crossed surfaces, namely eight in-range and two out-of-range from the indenter for comparison, are recorded and then superimposed.
In order to examine whether the 3D microstructure can help in discerning a similar indenter, we additionally analyze the normalized SED distribution of the following three indenters: a gap indenter (comprising two 0.5×2 mm 2 indenters with a 1.0×2.0 mm 2 gap between inside edges), a 90-degree-rotated gap indenter, and a solid circle indenter (with a diameter of 2 mm). The normalized SED (ε i ) at the i th sample of data is as follows where A represents the dataset of SED, and A i denotes the absolute value of SED at the i th sample. The normalized SED is contoured into a 10 × 10 table with each cell relating to a specific tip of dermal papillae (i.e., 100 tips in total). Only "hot-spots" that correspond to the positions that fall into the right side of the scale bar are plotted; i.e., the threshold corresponds to 0.6. The size is adjusted for better visualization of the distribution as opposed to considering the absolute value. Vibration experiments. The FA-I is sensitive to stimulus in range of 20 to 50 Hz, and thus indentation experiments with four vibrating stimuli (with different frequencies corresponding to 20 Hz, 30 Hz, 50 Hz, and 80 Hz) are conducted. The indenter is a 2×2 mm 2 plane. The amplitude is maintained at 0.5 mm. Previous studies [19] indicated the highest spike rates of FA-I when indenters reached their deepest position. Thus, the frequency-dependent SED response is plotted at a single FA-I near the edge of the indenter at the following time point (e.g., t = 0.005 s for 50 Hz stimuli). The thickness of the subcutaneous layer appears to strongly affect the result of line-load validation. The models with a thickness of the subcutaneous layer ranging from 0.6 mm to 1 mm well fit both the experimental data [7] and the linear-elastic model from [11] as shown in Fig  5c. Each set of 2D and 3D ridged models were treated equally, and thus, the elastic characteristics are similar in both conditions. The 1x100e model (that includes original dimensions with a 1-mm-thick subcutaneous layer; bottom-edge restrained) fit experimental data [7] the most. Hence, the 1x100e model was used in ramp-and-hold experiments and vibration experiments. and thus, the observed surface is shifted along the x-axis to the cross section where the SED is highest. The SED is high at the corners of indenters in the 2D ridged model (Fig 6a) although it diffuses separately into the dermal papillae in the 3D ridged model (Fig 6b). A high SED is observed at each tip of the dermal papillae.

Ramp-and-hold experiment results
The 3D distribution of SED at the tips of dermal papillae after pressing by the solid bar indenter is shown in Fig 7. The vertical axis presents the magnitude of SED, the horizontal axis presents the distance along the z-axis from the center, and the depth axis presents the distance along the x-axis from the center. The peak value of SED in the 3D ridged model exceeds that in 2D ridged model. A valley is observed at the center of the apex and indicates an interchange zone (peak-valley-peak) in both models. The contour of the normalized SED in the case of solid bar indenter reveals that the interchange zone of the 3D ridged model along the z-direction is more significant than that of the 2D ridged model (Fig 8). Along the x-direction, the interchange zone is equally noticeable in the two models.
A summary of the normalized dataset is shown in Fig 8. In the 2D ridged model, the SED is mostly high at the four corners of the indenters (for example, solid bar indenter and 90degree-rotated gap indenter) while it mostly high at the edges of the indenters (and especially along the z-axis) in the 3D ridged model. The 3D ridged model produces distributions of the normalized SED at the FA-I receptor location that differentiate the indenters. Furthermore, the 2D ridged model shows little difference between the normalized SED distributions of the solid bar indenter and gap indenter. It is evident that the 3D microstructure significantly influences encoding spatial differences among the indenters. In order to examine the change of SED when the frequency increases, the SED at a FA-I location (with the maximum SED value) is plotted and compared with the following position in the 2D ridged model (Fig 10). In the 3D ridged model, the SED increases when the frequency increases to 80 Hz while it remains constant in the 2D ridged model.

Discussion
Previous studies showed the effect of microstructure on stress/strain behavior at the position of mechanoreceptors (SA-I and FA-I) where the microstructure amplifies the SED concentration at the SA-I receptor's location. However, most studies were conducted in two dimensions since they lack the depth configuration of dermal papillae and overlooked the effect of the FA-I receptor positions. The present study involved proposing and developing a model that presents the 3D microstructure of skin in significant detail (i.e., a 3D ridged model). The boundary conditions were determined through a validation experiment with a standard line load method (Fig 5). The viscous parameters were determined by fitting the normalized contact force on the circular indenter to the experimental data reported by [17]. The approach reduces the calculation cost while ensuring that the performance of the model is comparable to that of other models. The limitations of the model will be discussed in the final part of this section.
Our main finding indicates that the 3D microstructure appears to modify the distribution of SED at FA-I receptor's positions such that it helps in discerning the spatial configuration of similar indenters (Fig 8). This finding agrees and interprets the conclusion of a recent study by Harih et al. [20] in which the 2D and 3D fingertip FE models were justified during static contact simulation. Although the results of the 2D and 3D model in [20] were similar (possibly due to the lack of microstructure), the 3D model provides additional insights into the third dimension. The results suggest that the 3D microstructure provides focalizing effects and especially localizing effects as described below.

The diffusion of SED and the focalizing effect
The focalizing effect describes the trend of SED concentration at the tips of dermal papillae (FA-I receptor's positions). The effect occurs when the contact surface exhibits uneven 3D microstructure and SED concentration at a mechanoreceptor location geometry and the stiffness of the epidermis and dermis are significantly different. Gerling [21] demonstrated the lensing effect of ridged geometry in a 2D model of a microstructure and concluded that stress is more concentrated at the bases of intermediate ridges. In the 3D model of the microstructure, we confirmed a similar effect given that a high concentration of SED was present. at the tips of the dermal papillae in Fig 6. Localizing effect of 3D microstructure and its potential roles in spatial discrimination The localizing effect describes the peaks of the 3D distribution of SED located at the tips of dermal papillae near the edges of the indenter (Fig 7). The contrast between these positions and others produced a specific pattern that reflects the shape and the size of the indenter (Fig 8). The localizing effect was confirmed in both ramp-and-hold experiments (Fig 8 and S1 Fig)  The 3D ridged model appears to be sensitive to the edges of the indenter with a valley of SED values at the center. In the case of the solid bar indenter, the interchange zone in the xaxis (peak-valley-peak) is more noticeable than that in the z-axis (Fig 8). This behavior may be related to the presence of fingerprints along the x-axis given that the fingerprints are believed to affect the stress/strain concentration at the dermal papillae as shown in [9].
Interestingly, the normalized SED distributions at the FA-I receptor location were sensitive to the edges of indenters in most experiments (see Fig 8). This finding does not contradict the hypothesis that the SA-I receptors are more sensitive than FA-I receptors in terms of shape discrimination (for details, see [5,22]). However, it indicates that the mechanical response at FA-I receptor location might provide complementary information for better discrimination under dynamic conditions.

Effect of 3D microstructure in vibratory conditions
In the 3D model of the microstructure, the SED response at a FA-I location increased when the vibratory frequency increased from 20 Hz to 80 Hz (Fig 10), irrespective of the widelyknown maximum sensitive frequency of the FA-I receptor (50 Hz). The increase is understandable since the viscous force is proportional to velocity. In this study, the amplitude was kept constant in all experiments, and thus the higher frequency stimulus is displaced at a higher velocity.
Interestingly, the SED response at a FA-I location in the 2D ridged model remained unchanged (Fig 10). Evidently, either the increase was extremely small in the 2D ridged model or the energy was further concentrated at other positions, e.g. SA-I locations. Meanwhile, in the 3D ridged model, the increase in SED response at a FA-I location was high due to the focalizing effect. This result indicates the requirement for a 3D model of microstructure for FA-I related analysis in vibratory conditions. An electrophysiological report published by Bensmaia [23] showed that the spike rates also gradually increased from 20 Hz to 80 Hz, and this was similar to the SED behavior observed in our study (Fig 11). This result is in agreement with the results obtained by previous studies that suggested a relationship between the dynamic spike rate and local SED at a mechanoreceptor's location [17] (in this case, this is at the FA-I location).
The potential use of 3D microstructure in tactile sensor design In terms of tactile sensor design, the elastic cover is highly essential to protect the sensor from being damaged by external forces [24]. However, Shimojo [25] demonstrated that the elastic cover significantly reduced the spatial resolution of the sensor even if the cover thickness is only 0.2 mm. Various studies mimicked the natural geometry of human finger microstructure in sensor design to improve tactile sensing ability while retaining protection such as studies that indicated that artificial epidermal ridges enhanced strain gage sensitivity by 1.8 times when compared with tactile sensors without ridges [26] or enhanced tactile shape discrimination [27,28]. An optical tactile sensor developed by Chorley et al. [29] with an array of identical urethane pins as artificial dermal papillae were capable of encoding edge information.
Here, we propose the potential use of a 3D microstructure in tactile sensor design. Given the assumption that the sensing elements respond equally, placing the same at different positions of 3D microstructure allows us to achieve different effects intentionally. The positions at the tips of dermal papillae enhance the spatial discrimination of the tactile sensor (our finding) while the positions at the bases of the intermediate ridges enhance the sensitivity (as suggested in [11]). The contour of the 3D ridged model in Fig 8 implies the tactile image obtained by embedding sensing elements at the tips of dermal papillae under the following circumstances.

The limitations of proposed model and future works
The limitation of the model is that the model is small, and this restricts the size of the indenters. The current model cannot be used to show the response at sides (left-right) of the fingertip. It is necessary to develop a large-scale model with respect to the complexity and computational cost to facilitate comparisons with physiological or electrophysiological experiments. Physiological or electrophysiological experiments are usually conducted with indenters larger than 2 mm. In this study, the neuronal responses of FA-I mechanoreceptors were not examined due to the lack of a neuronal model. Additionally, the neuronal response of FA-I mechanoreceptor is suggested as a combination of population activities [23,30]. It is unknown as to precisely how the FA-I mechanoreceptor modifies responses during the combination. Hence, the determination of the manner in which the FA-I mechanoreceptor processes the population mechanical responses into neural responses is an interesting challenge for future work. To investigate the neuronal responses of FA-I mechanoreceptors, one could use a transduction function to transform the SED into neural current. The transduction function was assumed as a sigmoidal stimulus-current curves in [31], and as a two linear function in [6]. Then the spike times could be obtained by using the leaky integrate-and-fire model.

Conclusion
We developed a model with the 3D microstructure of skin and compared it with conventional 2D microstructure to examine the effect of natural configuration on the SED distribution of mechanoreceptors. The analysis was conducted with various types of indenters under dynamic conditions. The results indicate that the 3D microstructure modifies the SED at FA-I mechanoreceptors such that it helps in discerning similar indenters. The 3D microstructure enhanced the SED response at a FA-I mechanoreceptor when the stimulus frequency increased.
The cover design of available tactile sensors would benefit from our findings. The model produced a specific response pattern of SEDs corresponding to each indenter. This feature along with the currently popular neural sub-model may allow us to effectively discriminate between similar indenters.