Computational approach to understand temporal and spatial tactile transmission processes from mechanical stimuli of the index fingertip to the primary somatosensory cortex

Mechanisms of information transmission using tactile sense are one of major concerns in producing simulated experience in virtual or augmented reality as well as in compensating elderly or impaired people with diminished tactile sensory function. However, important mechanism of the difference of peak latency in the primary somatosensory cortex (SI) between electrical and mechanical stimulations of finger skin is not fully understood. We propose a computational approach to fuse a computational model to simulate temporal and spatial transmission processes from mechanical stimuli to the SI and experimental method using a magnetoencephalograph (MEG). In our model, a tactile model that combined a three-dimensional mechanical model of fingertip skin and a neurophysiological model of a slowly adapting type 1 (SA1) mechanoreceptor was integrated with a somatosensory evoked field (SEF) response model. Electrical and mechanical stimulations were applied to the same locations of the right or left index fingertips of three subjects using a MEG. By identifying parameters of the SEF response model using the electrical stimulation test data, predicted first peak latency due to a mechanical stimulus was identical to its average value obtained from the mechanical stimulation test data, while the spatial map predicted at the multiple SA1 receptors qualitatively corresponded to the MEG image map in the timings of peak latency. This suggests that mechanical change in the skin and neurophysiological responses generate the difference of peak latency in SI between electrical and mechanical stimulations. The computational approach has the potential for detailed investigation of mechanisms of tactile information transmission.


Introduction
The perception of touch entails mechanical interactions with the external physical world, including contact with objects and tactile scans of textured surfaces. Perceived tactile information allows individuals to control forces generated by our hands and fingers as well as their motions when detecting contact, when grasping and manipulating objects delicately, and when evaluating textures. Some patients who suffer strokes have impaired ability to grasp and lift objects because of tactile sensory or motor deficits (Blennerhassett et al., 2007;Hermsdörfer et al., 2003;Carey et al., 2011). Likewise, elderly people tend to have less fine texture discrimination ability and hand strength because of diminished tactile sensory and motor function (Skedung et al., 2018;Carmeli et al., 2003). Therefore, it is critical to understand the mechanisms of tactile perception from mechanical stimuli to the primary somatosensory cortex (SI) in producing simulated experience in virtual or augmented reality as well as compensating for tactile sensation in stroke patients as well as in elderly people with diminished tactile sensation.
In the mechanism of tactile perception, the latency from the stimuli to the central nervous system is one of the most important functions of somatosensory perception. A few experimental studies (Hari et al., 1993;Yamashiro et al., 2009;Hesse et al., 2010;Onishi et al., 2010) used magnetoencephalographs (MEG) to investigate a somatosensory evoked field (SEF) waveform and the peak latency of the SI when the skin pad of the hands or fingers are stimulated electrically or mechanically. Yamashiro et al. (2009) stimulated the dorsum of the right hand electrically with a train of current-constant square wave pulses (pulse duration, 0.5 ms) and the peak latency of the early SI was about 34 ms in average, and Hesse et al. (2010) also stimulated the lateral dorsum of the left hand mechanically by being tapped once every 2-3 s and the peak latency of the SI was 55 ms in average. Hari et al. (1993) stimulated distal volar phalanxes of the fingers electrically with 0.2 or 0.3 ms constant current pulses and the peak latency of the SI was 30 ms in average, and Onishi et al. (2010) also stimulated the tip of the right index finger mechanically with a constant stimulus of 1000 ms duration and the peak latency of the SI was about 58 ms in average. Comparing Yamashiro et al. (2009) and Hesse et al. (2010) which stimulated the dorsum of the hand, the difference of peak latency between electrical and mechanical stimulations is estimated to be 21 ms by using difference of the average values. Comparing Hari et al. (1993) and Onishi et al. (2010) which stimulated the fingers, the difference is estimated to be 28 ms. Thus, the difference of peak latency between electrical and mechanical stimulation is estimated to be approximately 25 ms in average. However, different experimental setups and different stimulation methods were used in these studies. There are no experimental studies that stimulated the same location on a fingertip electrically or mechanically under the same experimental conditions using a MEG measurement system. Therefore, the mechanism of the differences in the peak latency of the SI between electrical and mechanical stimulations are still unclear.
Previous studies Sripati et al., 2006;Kim et al., 2010;Saal et al., 2015Saal et al., , 2017 suggested that the mechanism of tactile perception involves tactile data produced by mechanical stimulation that are acquired through mechanoreceptors within the skin, encoded as trains of action potentials or spikes with unique firing properties, which are transmitted to the central nervous system, especially to the somatosensory cortex. Some computational models (Sripati et al., 2006;Kim et al., 2010;Maeno et al., 1998;Kobayashi et al., 2012;Gerling et al., 2014;Pham et al., 2017;Hamasaki et al., 2018) were proposed to represent the mechanism of tactile perception; computational models can effectively represent some specific mechanical states (e.g., stress or stain) that are difficult to measure experimentally. Computational models are also used to determine the effect of age and gender as well as the shape, roughness, and stiffness of objects in the context of tactile perception. Maeno et al. (1998) developed a two-dimensional (2D) finite element (FE) model of human fingertip skin and predicted the stress distribution at the mechanoreceptors of Merkel cells, Meissner's corpuscles, Ruffini endings, and Pacinian corpuscles when a frictional force was applied. Kobayashi et al. (2012) developed an electronic circuit model to obtain SEF functions using electrical stimulations of the right thumb. They simulated the peak latencies of the SEFs at the primary somatosensory cortex (SI) that were similar to data measured using a 39-channel superconducting quantum interference device (SQUID) system. Gerling et al. (2014) developed a computational model of cutaneous skin and tactile neurons. They combined three-dimensional (3D) FE representations of the fingertip to mimic skin mechanics, a biphasic transduction function to convert stresses and strains in simulated skin into transmembrane current of slowly adapting type 1 (SA1), and a leaky-integrate-and-fire (LIF) model to generate the timing of the spikes. Then, they validated the model against mechanoreceptor population responses to psychophysical discrimination of the spheres. Nevertheless, there are no computational models to simulate spatial and temporal transmission processes from mechanical stimuli to the SI. Therefore, previous computational models are not used for better understanding of the mechanism of the difference in the peak latency of SI between electrical and mechanical stimulations of the index fingertip.
The purpose of this study was to propose a computational approach for better understanding of the mechanism of the difference in the peak latency of SI between electrical and mechanical stimulations. For this purpose, we developed a computational model for simulating temporal and spatial transmission processes from mechanical stimuli of the index fingertip to the SI, and conducted a series of experimental tests in which stimulations were electrically or mechanically applied to the same location of the right or left index fingertip of each subject using a MEG measurement system. Then we discussed the mechanism of the difference of the peak latency from comparisons between simulation results and test data.

Computational model
We developed a computational model combining three existing computational models: a mechanical 3D finite element model of the human index fingertip skin to simulate skin deformation by mechanical stimuli; a neurophysiological model including a transduction model to convert stresses or strains in the simulated skin into the transmembrane current at the mechanoreceptors of slowly adapting type 1 (SA1) and a model to simulate the membrane potential of SA1 receptors by the transformed current; and a SEF response model of the SI using an electronic circuit model to simulate the transmission of neural signals generated by the membrane potential to the SI. Fig. 1 shows the outline of the developed computational model that simulates a series of flows in tactile perception from the mechanical stimulation of the SI.

Mechanical model
We developed a 3D FE model of the human index fingertip skin and validated it against the reaction force during its contact with a flat surface (Hamasaki and Iwamoto, 2019). The external shape and position of the bone in the index fingertip skin model were created by referring to a model constructed from computed tomography images (Shimawaki and Sakai, 2007). The skin of the fingertip model was composed of four tissue layers: the stratum corneum, epidermis, dermis, and subcutaneous tissue, with thicknesses of 0.17 mm (Egawa et al., 2007), 0.65 mm (Maeno et al., 1998;Shao et al., 2010), 0.8 mm (Maeno et al., 1998;Shao et al., 2010), and 0.7 mm (Shao et al., 2010), respectively. We focused on SA1 mechanoreceptors, which are evaluated by tactile receptors. The SA1 receptor is located at the boundary space between the epidermis and dermis. This 3D model was discretized into finite elements using a 10-node tetrahedron element type. There were approximately 50,600 elements with 74,200 nodes. The stimulator had a sharp tip, the shape of which was similar to the tip of a toothpick in the mechanical device used in the experimental test described in Fig. 6, and stimulated the index fingertip pad of the model. The mechanical model in the proposed computational model was used to simulate contact between the fingertip skin and the stimulator, and the mechanical states at which the tactile receptors were calculated. Fig. 2 shows the fingertip and stimulator.
The material models and properties applied to the FE model shown in Fig. 2 were determined referring on those in the literature (Shao et al., 2010;Wu et al., 2006). A first-order polynomial hyperelastic model (neo-Hookean model) was adopted for the stratum corneum and epidermis, and a second-order polynomial hyperelastic model was adopted for the dermis and subcutaneous tissues. Bone and nail were assumed to be linear elastic materials. 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. In the case of N = 1(i = 1, j = 0), a first-order polynomial model is applied, whereas in case of N = 2, a second-order polynomial model is applied. Table 1 shows the material properties adopted for each material. The stimulator was defined as a rigid body. We developed a 3D FE model of the human index fingertip skin and validated it against the reaction force during its contact with a flat surface (Hamasaki and Iwamoto, 2019).

Neurophysiological model
The neurophysiological model was divided into two parts: 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. We used a mathematical model proposed by Gerling et al. (2014) as the transduction model. However because the current of SA1 attenuates exponentially when the stimulation intensity is constant (Ikeda et al., 2014), we modified the transduction model to represent the attenuation of SA1 responses. We tried to use two types of membrane potential models that are widely used for simulating the behavior of neural

Table 1
Material properties of the mechanical model (Shao et al., 2010;Wu et al., 2006 systems, that is, (1) the leaky integrate-and-fire (LIF) model, which is the most widely used because of high computational efficiency but a kind of rudimentary model, and (2) the Izhikevich (IZH) model (Izhikevich, 2006), which is widely adopted to emulate SA1 firing responses (Wei et al., 2019;Zhengkun and Yilei, 2017) with more biologically accurate features like Hodgkin-Huxley (HH) model. The transduction model (Eqs. (2)-(5)) used in this study are presented below: 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 Eq.n (3) immediately before dU(t)/ dt ≤ 0, t peak is the time immediately before dU(t)/dt ≤ 0, τ is the time constant, β, K S , and K d are the parameters for transduction, and Δt is the time increment. The membrane potential models in this study are presented below. In the LIF model, V(t), R, and C are the membrane potential, membrane resistance, and membrane capacity of the SA1 receptor, respectively. In the LZH model, u(t), k, a, b, V r , and V t are the recovery current, neuron's rheobase, recovery time constant, resonance, resting membrane potential, and instantaneous threshold potential, respectively. Eqs. (6)-(8) are solved using a fourth-order Runge-Kutta numerical integration method. Fig. 3 shows the SED at the tactile receptor and firing responses of the receptor during the mechanical stimulation of the fingertip pad similar to that used in the experimental tests described in Fig. 6, which were simulated using the developed neurophysiological model integrated with the mechanical model of the human index fingertip skin. The first firing occurred after a delay of 10 ms from the onset of the stimulation (0 ms). The five sequential firing responses occurred by 25 ms when the stimulation turned to unloading. In contrast, firing responses did not occur during the unloading of stimulation. The neurophysiological model integrated with the mechanical model of the human index fingertip skin was also validated against the experimental test data from (Srinivasan and LaMotte, 1991) and reproduced the firing responses of SA1 receptor when the fingertip pad was indented by cylindrical bars with curvatures of 16, 8, 4, and 2 [1/in] (Fig. 4). Although the values of C, R, and V, which are the thresholds for the LIF model (when V > V, V is reset to zero) in Eq. (6), were adopted from a previous study , the parameters β, K S , and K d of the transduction model (Eqs. (2) and (3)) were re-identified to represent the experimental data (Srinivasan and LaMotte, 1991). The parameters R and V were identified to 66.11 × 10 6 Ω and 50.10 mV, respectively. In addition, the parameters β, K S , and K d of the transduction model (Eqs. (2) and (3)) for the IZH model were identified to represent the experimental data (Srinivasan and LaMotte, 1991). Table 2 lists the identified parameters for LIF and IZH models. For IZH models, the parameters k, a, b, V r and V t were identified to 0.7 pA/mV 2 , 0.03 m s − 1 , − 2, − 60 mV, and − 40 mV, respectively. Fig. 4 shows the comparisons between the simulation results predicted by the mechanical model and neurophysiological model (LIF or IZH model) as well as the experimental test data obtained from Srinivasan and LaMotte (1991). The comparisons found that the responses from the two models were comparable. Previous study (Valadez-Godínez et al., 2017) demonstrated that LIF model has the potential to show the same performance as HH model regarding spike-timing. Therefore, we adopted the LIF model with high computational efficiency to the following simulations.

SEF response model
The SEF response model was developed to predict the peak latency of SI based on an electronic circuit model (Fig. 5) and a transfer function (9) proposed by Kobayashi et al. (2012). L is a coil, C is a capacitor, R 1 , R 2 , and R 3 are resistances, and D is a delay element to describe latency. Note that C and R, shown in Eq. (9), are not related to those in Eq. (6). The SEF response model receives neural signals generated by the membrane potential from the neurophysiological model integrated with the mechanical model of the human index fingertip skin as inputs and outputs of the SEF responses at the SI: Fig. 3. The SED at a tactile receptor and firing responses of the receptor during mechanical stimulation of the fingertip pad.

Measurement of brain activity elicited by tactile stimulation
To validate the abovementioned computational model and to simulate the spatial and temporal transmission processes from mechanical stimuli to the SI, some experimental data on the latency from the stimuli to the central nervous system are needed. A few experimental studies (Hari et al., 1993;Yamashiro et al., 2009;Jousmäki et al., 2007;Onishi et al., 2010) used magnetoencephalographs (MEG) to determine the SEF waveform, and the peak latency of the SI when the skin pads of the hands or fingers are stimulated electrically or mechanically. Nevertheless, there are no experimental studies on the peak latency of the SI that stimulates the same location on a fingertip electrically or mechanically under the same experimental conditions using a MEG measurement system. Therefore, the differences in the peak latency of the SI between electrical and mechanical stimulations remain unclear.
We conducted a series of tests using a MEG measurement system to obtain the peak latency of SI and brain activity when the right or left index fingertip skin pads of three subjects were stimulated either electrically or mechanically. Test procedures and data processing using the MEG system are described below.

MEG recordings
The SEF were recorded with a whole-head 306-channel sensor array (TRI-UX, ELEKTA Neuromag, Stockholm, Sweden), which consisted of 102 identical triple sensor elements, in a magnetically shielded room at Kizawa Memorial Hospital. Each sensor element consisted of two orthogonal planar gradiometers and one magnetometer coupled to a multi-SQUID, thereby providing three independent measurements of magnetic fields. We analyzed the MEG signals recorded by 204-channel planar-type gradiometers. A band-pass filter for the recording was set to 0.01-330 Hz, and the sampling rate was 1 kHz. After completion of the recording, movement compensation was performed by temporal extension signal-space separation (tSSS) using Maxfilter 2.2 software (Neuromag; Elekta, Stockholm, Sweden) and applied off-line to the recorded raw data to reduce artefactual signals arising from the outside of the sensor array and to correct for the changes in the head position and associated movement-related artifacts (Medvedovsky et al., 2007;Taulu et al., 2004Taulu et al., , 2005. Before MEG recordings, three anatomical landmarks (nasion and bilateral preauricular points) and five head position indicator (HPI) coils were attached to specific sites on each subject's head, and a threedimensional (3D) digitizer was used to measure the anatomical landmarks of the head with respect to the HPI coils. The precise location of the head with respect to the sensory array was determined using the HPI coils.  (Srinivasan and LaMotte, 1991).

Table 2
Parameters determined for the neurophysiological model. The anatomical landmarks provide the spatial information necessary for the integration of magnetic resonance imaging (MRI) and MEG data. MRI scans were obtained from all subjects using a 3.0-T Philips Achieva scanner (Philips Medical Systems, Best, NL). A high-resolution T1weighted anatomical image was acquired using a 3D fast field echo sequence with the following parameters: TR = 6.9 ms; TE = 4.6 ms; flip angle 8 • ; field of view: 240 mm × 240 mm; matrix: 240 × 240; slice thickness: 1 mm with no gap; in-plane resolution of 1.0 mm × 1.0 mm. The structural volumes were used for MEG co-registration and spatial normalization.
Three volunteers from our department (two males and one female, aged 29-38 years) participated in this study. All participants were righthanded (Edinburgh Handedness Test;Oldfield, 1971) and had no known neurological disorders. Written informed consent was obtained from all participants. The study was approved in advance by the ethics committee of Kizawa Memorial Hospital and Toyota Central R&D Labs., Inc.

Test procedures
Each subject was seated comfortably inside a magnetically shielded room with the head firmly positioned inside the MEG system. Then, two types of tactile stimulations, electrical or mechanical, were applied to the index fingertip pulp of the subject's right or left hand. Fig. 6 shows the electrical and mechanical stimulation devices used in the study. For electrical stimulation, a two-point stimulus was delivered using a striptype Ag electrode (1 mm in width and 0.1 mm in thickness, anode 2 mm from cathode) was applied to the index fingertip pulp surface of the right hand with a constant stimulus of 2 ms. The two-point stimulus distance was determined as two conditions (2 and 4 mm) based on the two-point discrimination of human fingers described in the literature (e.g., Koo et al., 2016). Because we focused on one-point stimulation, time differences between two-point stimuli were set to zero. For mechanical stimulation, air pressure was applied to the output portion (air-puff stimulator) of the mechanical stimulation device inside the shielded room through a tube from the outside of the room. Because the pressure was weak (less than 10 gf), the subjects could not easily perceive it using their tactile sensors. Using a mechanical stimulation device with parallel leaf springs and a toothpick, infinitesimal deformation of the air-puff stimulator with a pneumatic load can be converted to a stroke to be perceived by each subject's index fingertip. One-point stimulation was applied to the index fingertip pulp of the subject's hand with a stimulus of 50 ms duration, and the MEG was recorded. Because changes in the SEF are very small, especially in mechanical stimulations, the data obtained on the SEF were processed using signal averaging. For mechanical stimulations, signal averaging was applied over 150 trials, while it was applied over 100 trials in electrical stimulations to improve the rejection of the background noise not related to the stimulations. The interstimulus interval (ISI) was set to 4.0 s because Jousmäki et al. (2007) reported that prolongation of the ISI from 2 to 4-5 s elicited additional activation in the second somatosensory (SII) cortex despite the fact that the typical ISI is 1.5-2.0 s. During the recording, subjects were instructed to wear eye masks to reduce the number of eye movements and earplugs so as not to hear the stimulation timing identified by the sounds from the pneumatic drive in the mechanical stimulations. There were four experimental test conditions: Case Nos. 1 and 2 were electrical stimulations with distances of 2 and 4 mm, respectively. Case Nos. 3 and 4 were mechanical stimulations of the right and left hands, respectively. In the electrical stimulations, the onset timing of the stimulation was 50 ms after the onset time of the MEG recording while the air was applied for approximately 50 ms from the onset time of the MEG recording in the mechanical stimulations.

Data analysis
Data analysis was performed using SPM12 (statistical parametric mapping software; Wellcom Department of Imaging Neuroscience, London, UK) in MATLAB (Math Works, Natick, MA). To analyze SEF, the band-pass filter was set between 2 and 100 Hz and the band-stop filter was set between 55 and 65 Hz; then, these filters were applied to SEF processes using signal averaging. In SPM12, to identify the sources of the evoked activities, variational Bayesian equivalent current dipoles (VB-ECDs) are used to circumvent the drawbacks of the classical ECD using a simple best-fitting optimization method with least square error criteria (Litvak et al., 2011). A template head model with meshes of normal sizes was adjusted to the individual head shape obtained from the MRI data measured for each subject by identifying the left and right preauricular points as well as the nasion. To create the electromagnetic forward model based on Maxwell's equation, a single-shell sphere model, which is aligned to the inner skull surface, was selected, and then inverse reconstruction was performed. The reconstruction of MEG sources is an ill-posed problem because the solution does not depend continuously on the data and there is no unique solution in the absence of prior information or constraints. Thus, the empirical Bayesian framework based on hierarchical linear models was used to obtain a solution with several priors or multiple constraints in SPM12 (Mattout et al., 2006). Once the inversion was complete, the time course of the source with maximal activity was obtained. The time course of the source was converted into the peristimulus time course, in which the stimulation timing was set at 0 ms. Fig. 7 shows the posterior probability maps (PPMs) and the magnetic field time history at SI, which was identified as a source of evoked maximum activation obtained from the MEG data using SPM12. The left side ( Fig. 7(a)) shows the time course of the source with maximal activity at 94 ms after MEG recording; i.e., the time when the maximum activation was generated at SI; the right side ( Fig. 7(b)) shows the maximum intensity projections (MIPs) of the PPMs at the same time. The data were obtained from Case No. 1 of Subject 1. The MIPs were projected onto a glass brain with three orthogonal planes. The PPMs in Fig. 7(b) are the images of the probability or confidence that an activation exceeds a specific threshold in which Bayesian inference is performed with prior variances estimated from the experimental data of this study based on Friston and Penny (2003), Friston et al. (2002) and Friston et al. (2002). In Fig. 7(b), the voxels with 96% confidence that the effect sizes at voxels are greater than 1% of the global mean are represented in black. Fig. 7(a) indicates the time course of the source with maximal activity; i.e., the position of the red circle in Fig. 7(b) at 94 ms, which corresponds to the index finger area of the contralateral SI according to the functional anatomy map of the somatosensory and motor cortex (Penfield and Boldrey, 1937). The time course in Fig. 7(a) is supposed to be a SEF waveform, that is, a magnetic field evoked on the index finger area of SI; it is totally different from the SEFs waveforms of SI in previous studies (Hari et al., 1993;Yamashiro et al., 2009;Jousmäki et al., 2007;Onishi et al., 2010), which were directly measured from the MEG sensors attached to a subject's head surface. The SEF waveforms of Case No.1 or 2 of Subjects 1 or 2 in electrical stimulation tests and those of Case No. 4 of Subjects 1, 2, or 3 were obtained from the PPMs (Figs. 9, 10 and S1). Then, averaged curves were obtained for the SEF waveforms in electrical and mechanical stimulations (Figs. 9 and 10). In addition, time-frequency analyses were applied to the averaged SEF curves of electrical and mechanical stimulations to investigate the strength of SI activity. The time-frequency representation maps were estimated for the averaged SEF curves of electrical and mechanical stimulations by computing the spectrogram distribution. The spectrogram distribution corresponds to the squared modulus of the short-time Fourier transform (STFT) and a spectral energy density of the locally windowed signal x(u)h(u − t) is obtained as follows: The window h of the STFT is the smoothing window function (here we used Gaussian window of 23 points). x(u) is an averaged SEF curve of electrical or mechanical stimulation. i is the imaginary unit and ν is the frequency. The analysis was performed in MATLAB by using the Time-Frequency Toolbox (http://tftb.nongnu.org/).

Simulation setup
We used the developed computational model described in Fig. 1 and simulated the spatial and temporal tactile transmission processes from one-point mechanical stimulation to the SI using the measured test data described in Fig. 6. The simulation conditions were set to reproduce the test conditions.
In electrical stimulation, we only used the SEF response model to simulate the temporal tactile transmission process from one-point electrical stimulation to the SI because somatosensory perception by electrical stimulation is not perceived at receptors inside the skin, but rather it is generated by the elicitation of nerve firing. Then, the median nerve directly transmits the neural signal to the SI via the spinal cord and thalamus when the index fingertip is stimulated electrically. The parameters were re-identified using the four SEF waveforms in electrical stimulations of the study described in Fig. 6 (see Fig. S1). We used Case No. 1 or 2 of Subjects 1 or 2 for re-identification) ( Table 3). Test data for Case No. 1 or 2 of Subject 3 were not available because we could not obtain the PPMs of SI for Subject 3. This is probably because the visual cortex or other functional parts in the brain were more activated than the somatosensory cortex during the experimental tests of Subject 3. The parameter identification method is shown as follows: First, L, C, R 1 , and D were determined using the back-calculating Eq. (9) to represent the SEF waveform measured in this study. Second, using the calculated parameters as the initial values, optimization was performed to minimize the difference between the measured SEF waveform and the waveform obtained from Eq. (9). The optimization was performed using the fminsearch function of MATLAB. In the optimization, the SEF waveforms were normalized because the units of the waveforms were different between the experiment and Eq. (9). Table 3 lists the initial and identified parameters for an average curve of the four SEF waveforms in the electrical stimulations. Fig. 8(A) shows an image of the simulation results using the SEF response model under the condition of electrical stimulation. The peak latency of SI in the electrical stimulation was obtained from the peak time of the simulated waveform predicted by the SEF response model (Fig. 8(A)(1)).
In the mechanical stimulation, we used a computational model including the three component models to simulate the spatial and temporal tactile transmission process from one-point mechanical stimulation to the SI. In the mechanical stimulation experiment, the stimulator was pressed vertically against the index fingertip pad with the force that the subject felt weak, and it took 50 ms to load and unload the stimulation. To simulate this condition, a forced displacement of 0.5 mm was applied to the stimulator of the mechanical model for 25 ms, and then the displacement was returned to zero for another 25 ms, as shown in Fig. 3. This mechanical simulation using the 3D FE model of the human index fingertip skin was performed by using the commercial FE code Abaqus (SIMULIA, Dassault Systemes, Vélizy-Villacoublay, France). After the mechanical simulation, the firing responses of SA1 were simulated using the calculated SED at SA1 and the neurophysiological model (Fig. 3), and the obtained firing responses were inputted into the SEF response model to simulate the SEF waveform in the mechanical stimulation. Fig. 8(B) shows an image of the simulation results using the computational model under the condition of mechanical stimulation. As mentioned previously, the first firing occurred after a delay of 10 ms from the onset of the stimulation (0 ms), which corresponds to the latency from the mechanical stimulation to the timing of the first spike at SA1 (Fig. 8(B)(2)) and the six firing responses of SA1, including the first firing and the sequential five firing responses correspond to the latency due to the mechanical stimulus duration of 25 ms calculated by the SEF response model (Fig. 8(B)(3)). Therefore, the peak latency of the SI in the mechanical stimulation was obtained from the sum of the latency from the mechanical stimulation to the first firing of SA1 ( Fig. 8(B)(2)) and the latencies due to the sequential five firing responses ( Fig. 8(B)(3)), in which each firing response was assumed to generate the latency to the SI shown in Fig. 8(A)(1). In this study, the number of SA1 receptors was set to 99 (red points in Fig. 12) around the mechanical stimulation point. We obtained a SEF response curve for a mechanical stimulation using signal averaging of the SEF response curves predicted for the 99 receptors. Because the magnitude of the firing pulse calculated from the neurophysiological model is a dimensionless quantity, the magnitude of the firing pulse was set as 1 V according to Kobayashi et al. (2012) when the firing responses were inputted into the SEF response model.
In addition, multiple SA1 receptors in the 3D FE model of the human index fingertip skin respond to mechanical stimulation. We assumed nodes in a part of the 3D FE model as the SA1 receptors, in which the number was set to 99, and the firing response of each SA1 receptor was transmitted to the corresponding cortical neuron in the SI by referring to scatter maps illustrating somatotopic and converging representations of fingers in the SI (Iwamura et al., 1980;Duncan and Boynton, 2007). Based on the simulations of the electrical or mechanical stimulations, we predicted the SEF waveform at the SI and spatial distribution of the SEF responses around the SI (hereafter, called the SEF map) in the mechanical stimulation. . 9 shows a comparison of SEF waveforms at the SI in the electrical stimulation between the simulation results predicted by the SEF response model with re-identified parameters (Table 3) and the average curve with the standard deviation of the experimental data of the four cases (Case No. 1 or 2 of Subjects 1 or 2) obtained from this study. The figure also includes the SEF response predicted by the computational model for mechanical stimulation. The SEF response model with the reidentified parameters almost reproduced the first peak latency of the SI (45 ms) obtained from the average curve of the experimental test data. Although the simulated waveform and the peak latency can be compared with the experimental test data, the amplitude of the waveform has no meaning and cannot be compared with the test data. In addition, the predicted latency between electrical and mechanical stimulations was 25 ms, and then the predicted latency from the mechanical stimulation to the SI was 70 ms. Fig. 10 shows a comparison of the SEF waveform at the SI in the mechanical stimulation between the simulation results predicted by the proposed computational model and the average curve with the standard deviation of the experimental test data of the three cases (Case No. 4 of Subjects 1, 2, or 3) obtained from this study. The figure includes the SEF waveforms of the experimental data of the three cases. The first peak latency of SI (70 ms) predicted by the computational model was identical to the average value of the first peak latencies from the three cases (70 ms), in which the first peak latencies were 88, 59, and 64 ms for Subjects 1, 2, and 3, respectively. The pattern of the SEF time history curve predicted by the computational model was different from that of the average curve of the experimental data. However, the predicted SEF curve had a positive value in the initial peak and a negative value in the second peak, similarly to the SEF curves obtained from the experimental data of the three subjects. Fig. 11 shows the time-frequency representation map for the averaged SEF waveform of the test data with electrical stimulation and that of the test data with mechanical stimulation. The timing of the peak latency of SI was 56 ms after electrical stimulation with the frequency range 0-50 Hz and 97 ms after mechanical stimulation with frequency range 0-50 Hz. Fig. 12shows the comparisons of the SEF map at the SI in the mechanical stimulation between the computational model and experimental test data of Case No. 4 of Subject 1, of which the SEF curve pattern (deep blue line) was relatively similar to that predicted using the computational model (red line) in Fig. 10. The red points in the 3D FE model of the index fingertip skin of Fig. 12 show SA1 receptors. The SEF map predicted by the computational model at points A and B shown in Fig. 10, that is, at 88 and 212 ms were compared with the MEG data of Case No. 4 of Subject 1, which were obtained using SPM12. A SEF map was depicted by plotting the SEF value of each receptor at each time. We assumed that the SEF maps correspond to magnetic fields at the cortical neurons around the SI obtained from the MEG data. The SEF map predicted by the model corresponds to an upper right quarter of a large rectangle shown in the MEG data of Fig. 12. The SEF map predicted by the model was similar to the MEG data at the contralateral SI, although the SEF responses of parts other than the contralateral SI were not simulated in this study.

Discussion
We proposed a computational approach to fuse a computational Table 3 Parameters applied to the SEF response model. model to simulate temporal and spatial transmission processes from mechanical stimuli to the SI and experimental method using a MEG, which is critical to determine how to produce simulated experience in virtual or augmented reality as well as compensate for elderly people or impaired humans with diminished tactile sensation. First, we discuss the validity of the peak latency of SI obtained from a series of experiments conducted in this study. The first peak latency of the SI for electrical stimuli in previous studies (Hari et al., 1993;Yamashiro et al., 2009) was found from 20 to 40 ms, which were obtained using a MEG measurement system similar to the one used in the present study. In our study, the average value of the first peak latency of SI was 45 ms (see Fig. 9) in the electrical stimulation tests. The first peak latency of SI for mechanical stimuli in previous studies (Jousmäki et al., 2007;Hesse et al., 2010;Onishi et al., 2010) was 52-69 ms, also obtained using a MEG measurement system similar to the one we used here. In the present study, the average value of the first peak latency of the SI was 70 ms (see Fig. 10) in the mechanical stimulation tests. Therefore, the peak latencies of the SI in both electrical and mechanical stimuli obtained in this study were almost identical to the upper values of previous studies. As mentioned in Introduction, from previous studies (Hari et al., 1993;Yamashiro et al., 2009;Hesse et al., 2010;Onishi et al., 2010), the difference of peak latency between electrical and mechanical stimulation is estimated to be approximately 25 ms in average, although different experimental setups and different stimulation methods were used in these studies. In our study, the peak latencies for electrical stimuli and mechanical stimuli were obtained at the same location of the right or left index fingertip using the same experimental setups. The average value of the first peak latency of SI was 45 ms (see Fig. 9) in the electrical stimulation tests while the average value of the first peak latency of the SI was 70 ms (see Fig. 10) in the mechanical stimulation tests. The difference of peak latency of the SI between electrical and mechanical stimulations obtained in this study was 25 ms and is identical to the value estimated from previous studies. This indicates the validity of the peak latency of SI obtained from this study. We obtained the first peak latency of SEF waveforms generated from PPMs that have the potential to identify the brain activity at the exact anatomical location of various functional brain areas by referring to the functional anatomy map (Penfield and Boldrey, 1937), for instance, with accuracy in the anatomical location to identify the difference between the index finger and ring finger, which cannot be directly obtained from the conventional outputs of the SEF waveforms of 204-channels planar-type gradiometers used in previous studies. Further experimental studies with more subjects are necessary to investigate the      M. Iwamoto et al. difference between our method using PPMs and conventional methods to obtain the peak latency of SI activity.
In addition, the time-frequency analyses demonstrated that the timing of the peak latency of SI was 56 ms after electrical stimulation with the frequency range 0-50 Hz and 97 ms after mechanical stimulation with frequency range 0-50 Hz (Fig. 11). In previous study (Papadelis et al., 2011), the timing of the peak latency of SI is 21 ms after stimulation with the frequency range 50-150 Hz. There are two differences between this study and the previous study (Papadelis et al., 2011). First, the frequency range of this study was smaller with low band (0-50 Hz) than that of previous study with medium band (50-150 Hz). Second, the timing of peak latency of SI for electrical stimulation in this study was 35 ms delayed from that in previous study. The reasons of these differences are because the previous study stimulated the median nerve whereas this study stimulated the index fingertip. The median nerve stimulation produces strong and well-defined MEG waveforms and the short-latency responses are found less than 40 ms after stimulus onset (Papadelis et al., 2011). On the other hand, the electrical stimulation to the fingers produces relatively weaker MEG waveforms and the short-latency responses occurred at 36-38, 16-18 ms later than the prominent N20m deflection in the response to stimulation of the median nerve at the wrist (Hari et al., 1993). In addition, the timing of peak latency of SI was 45 ms in average for electrical stimulation and 70 ms in average for mechanical stimulation in this study. However, from time-frequency analysis performed in this study, the timing was 56 ms for electrical stimulation and 97 ms for mechanical stimulation. The reason of differences in the timings of peak latency of SI is due to the data processing method. The values of 45 and 70 ms were the average values for the first peak of SEF waveform of each subject in the electrical and mechanical stimulations, respectively while the values of 56 and 97 ms were the first peaks of an averaged curve for SEF waveforms of three subjects in the electrical and mechanical stimulations, which can be obtained from the first peaks of black solid lines in Figs. 9 and 10, respectively. Since the stimulations to the fingers produce weaker SEF waveforms and the waveforms have different first peak latencies among subjects, further experimental studies with more subjects are needed to determine the first peak timing of SI from the SEF waveforms obtained using PPMs.
Second, with regard to the validity of the proposed computational model to predict the peak latency of SI in tactile perception by mechanical stimulation of the index fingertip, as shown in Fig. 1, the proposed model includes three components: a mechanical model, a neurophysiological model, and a SEF response model. The mechanical model is used to simulate the contacts between the fingertip skin and the stimulator, and to predict the mechanical states of stresses or strains at the tactile receptor, SA1. The mechanical model was validated against the reaction force during its contact with a flat surface (Hamasaki and Iwamoto, 2019) and showed good agreement with experimental test data (Maeno et al., 1998;Seria et al., 1997;Shao et al., 2009). The mechanical model was also validated against the mean impulse rate of SA1 during dynamic scanning of textured surfaces with various dot spacings (Hamasaki and Iwamoto, 2019) and demonstrated that the temporal history of stress at SA1 predicted by the model, which was assumed as temporal changes in the number of impulses of SA1, strongly correlated with the experimental text data (R 2 = 0.93) obtained from Phillips et al. (1992). The neurophysiological model integrated with the mechanical model of the human index fingertip skin was validated against the firing responses of SA1 when the fingertip pad was indented by cylindrical bars with different curvatures and showed good agreement with experimental test data from (Srinivasan and LaMotte, 1991) (Fig. 4). The proposed model includes the SEF response model differently from previous computational models including mechanical and neurophysiological models (Sripati et al., 2006;Gerling et al., 2014;Saal et al., 2017). The SEF response model is used to simulate the temporal tactile transmission process from the SA1 to the SI, in which we assumed that the median nerve directly transmits the neural signal to the SI via the spinal cord and thalamus, similarly to the transmission process when the index fingertip is stimulated electrically. Thus, the parameters of the SEF response model were identified to reproduce the feature of an average curve of the SEF waveforms obtained from the electrical stimulation tests. Using the proposed model, we performed a simulation of the temporal transmission process from one-point mechanical stimulation of the index fingertip skin to the SI and predicted the first peak latency from a mechanical stimulus to the SI as 70 ms, which was identical to the average value of the first peak latencies of the three subjects in the mechanical stimulation tests. In addition, predicted SEF curve had a positive value in the initial peak and a negative value in the second peak, similarly to the experimental test data, for example, which correspond to A and B of Subject 1, respectively (Fig. 10). This result suggests that the proposed model has the potential to simulate the temporal transmission process from a mechanical stimulus to the SI.
Third, we discuss the potential of the proposed model to simulate the spatial tactile transmission process. The proposed model calculated firing responses at multiple SA1 receptors based on changes of those strains due to the skin deformation when one-point mechanical stimulation was applied to the index fingertip skin and converted the calculated firing nerve signal at each SA1 receptor to brain activity at the SI cortical neurons using the SEF response model. We obtained the brain activity map (SEF map) with reconfiguration according to the scatter maps, illustrating the somatotopic and converging representation of fingers in the SI (Iwamura et al., 1980;Phillips et al., 1988;Duncan and Boynton, 2007) shown in Fig. 12. The predicted SEF maps were compared with the magnetic field map obtained from the MEG data of Case No. 4 of Subject 1, and brain activity predicted by the proposed model qualitatively corresponded to those of the test data at 88 and 212 ms after a mechanical stimulus. This result suggests that the proposed model has the potential to simulate the spatial transmission process from a mechanical stimulus to the SI.
As described above, the proposed model has the potential to simulate temporal and spatial tactile transmission processes when mechanical stimuli are applied to human finger skin. From the assumption that the temporal tactile transmission process from the SA1 to the SI in mechanical stimulations is identical to the process from electrical stimuli to the SI and the SEF response model produces the latency at the SI in electrical stimulations, we found that mechanical change in the skin and neurophysiological responses generate the difference of peak latency in SI between electrical and mechanical stimulations. Using the mechanical model of the proposed model, we investigated the aging effect on tactile perceptions by increasing skin stiffness with aging (Hamasaki et al., 2018). Increasing skin stiffness with aging could change the pattern of stress or strain at mechanoreceptors such as SA1 and decrease the firing responses and the resulting S1 activity. The spatial acuity of skin at the fingertip deteriorates noticeably with age as assessed by two-point threshold measurement, whereas there is a decline in sensory nerve conduction velocity and in the amplitude of the sensory action potential with age (Wickremaratchi and Llewelyn, 2006). Using the neurophysiological model of the proposed model, we can investigate the effects of aging or impairment of the firing responses of the mechanoreceptors by changing the time constant τ or the parameters for transduction β, K S , and K d used in the neurophysiological model. In addition, elderly subjects have a higher threshold for perception of electrical stimuli compared with younger subjects (Wickremaratchi and Llewelyn, 2006). Using the SEF response of the proposed model, we can also investigate the effects of aging or impairment on the transmission speed of neural signals from the mechanoreceptors to the SI by changing the parameters of the delay element to describe latency D or resistances R 1 , R 2 , and R 3 . Additional experimental tests with elderly people using a MEG would be helpful to determine these parameters of the SEF response model. Therefore, the proposed computational approach can help estimate individual differences in tactile perception for mechanical stimulations and understand the mechanisms of tactile perception from mechanical stimuli to the SI in producing simulated experience in virtual or augmented reality as well as compensating for tactile sensation in stroke patients as well as in elderly people with diminished tactile sensation.
This study has some limitations. First, the proposed model did not reproduce the pattern and vertical axis value of the SEF time history curves obtained from the experimental data. Nevertheless, the predicted pattern had the tendency of the pattern of the measured SEF curves. Second, a single SA1 tactile receptor was evaluated to simulate the response of the tactile receptors, although multiple tactile receptors other than SA1, the SA2 (Ruffini endings), the RA1 (Meissner corpuscles), and the RA2 (Pacinian corpuscles) responded. However, the SA1 responses significantly affect temporal and spatial transmission processes when human finger skin is stimulated in a mechanical manner with sustained and slow compressive loadings adopted in this study. Further studies are necessary to investigate the difference in the vertical axis value of the SEF curves between the magnetic field measured using MEG and the electrical activity of cortical neurons at the predicted SI using the proposed computational model and to understand how the SEF waveforms of SI are affected by the firing response of the other types of tactile receptors, including Meissner corpuscles, Pacinian corpuscles, and Ruffini endings elicited by mechanical stimuli with shear or tensile loadings.

Conclusion
In conclusion, our findings suggest that the proposed computational approach has the potential to simulate both the temporal and spatial tactile transmission processes when mechanical stimuli are applied to human finger skin and investigate the mechanisms of tactile perception from mechanical stimuli to the SI. Using the proposed computational approach, we found that mechanical change in the skin and neurophysiological responses generate the difference of peak latency in SI between electrical and mechanical stimulations.
In a future study, we intend to expand the proposed model to a handsomatosensory cortex model combining the 3D fingertip skin model including multiple mechanoreceptors of not only SA1 but also SA2, RA1, and RA2, according to the densities and the corresponding neurophysiological and SEF response models for simulating grasping objects or scanning textured surfaces to produce simulated experience in virtual or augmented reality and improve the tactile ability of patients or elderly people with impaired tactile sensation.

Data availability statement
A human index fingertip skin model can be calculated using the commercial finite element (FE) code Abaqus (https://www.3ds.co m/products-services/simulia/products/abaqus/). The human index fingertip skin FE model and MATLAB codes for a neurophysiological model and a SEF response model are available with documents including how to use the models by sending an e-mail to the corresponding author. Four SEF waveforms of electrical stimulation tests (Case No. 1 or 2 of Subjects 1 or 2) used for re-identification are included on Supporting information. If readers are interested in experimental data other than those shown in this manuscript, they are also available by sending an email to the corresponding author.