Mathematical Model for Skin Pain Sensation under Local Distributed Mechanical Compression for Electronic Skin Applications

Skin pain resulting from mechanical compression is one of the most common pains in daily life and the indispensable information for electronic skin to perceive external signals. The external mechanical stimuli are transduced into impulses and transmitted via nerve fiber, and finally, the sensation is perceived via the procession of the nerve system. However, the mathematical mechanism for pain sensation due to mechanical stimuli remains unclear. In this paper, a mathematical model for skin pain sensation under compression is established, in which the Flament solution, the revised Hodgkin–Huxley model, and the mathematical model gate control theory are considered simultaneously. The proposed model includes three parts: a mechanical model of skin compression, a model of transduction, and a model of modulation and perception. It is demonstrated that the pain sensation degree increases with the compression amplitude and decreases with deeper nociceptor location in the skin. With the help of the proposed model, the quantitative relationship between compression pain sensation and external mechanical stimuli is revealed, which has a significant benefit in promoting the design and mechanism research of electronic skin with pain perception function.


Introduction
Skin, the largest human organ, wraps the surface of the body and is in direct touch with the external environment [1]. It has many essential functions including protecting, excreting, and thermoregulation [2], and one of the most important functions is sensing external stimuli, including mechanical [3,4] and thermal stimuli [5][6][7]. As illustrated in Figure 1, external mechanical stimuli could be perceived as the sensation of touching with the transduction of skin and the modulation of the nerve system [8]. However, extreme mechanical stimuli will cause unpleasant pain sensations, or even damage the skin [9,10].
There were two main types of theories explaining the pain sensation mechanism before the 1960s. The first one is the specific theory, which indicates that each kind of pain owns a specific pathway to the brain [11,12]. The second one is the pattern theory, which proposes that the pain information transmitted into the brain is coded via the spatiotemporal pattern of the impulses [13,14]. However, both these theories have their shortcomings: the specific pain theory cannot explain the fact that pain can be relieved by rubbing the injured skin, and the decoding mechanism for the pattern theory remains ambiguous. The proposition of gate control theory (GCT) [15] makes it possible to explain many experimental phenomena, but GCT cannot elucidate the relationship between external stimuli and the nerve impulse. There were two main types of theories explaining the pain sensation me before the 1960s. The first one is the specific theory, which indicates that each pain owns a specific pathway to the brain [11,12]. The second one is the patter which proposes that the pain information transmitted into the brain is coded spatiotemporal pattern of the impulses [13,14]. However, both these theories h shortcomings: the specific pain theory cannot explain the fact that pain can be rel rubbing the injured skin, and the decoding mechanism for the pattern theory ambiguous. The proposition of gate control theory (GCT) [15] makes it possib plain many experimental phenomena, but GCT cannot elucidate the relation tween external stimuli and the nerve impulse.
Although pain has been studied from the molecular level to the level of t nervous system for a long period, the computation models of pain sensation limited [16,17]. Researchers have developed mathematical models at differen molecular [18,19], cellular [20], and neuron network levels [11,21,22]. Britton e established a mathematical model for GCT and used it to explain the quality of et al. [5,8,16] have systematically studied the thermal pain sensation. By comb transduction, transmission, modulation and perception process, the temperatu mal stress, and chemical burn are all taken into consideration in thermal pain Based on Xu's theory, Yin et al. [9] studied the skin pain sensation under the h epidermal electronic devices.
There are three main pain stimuli: thermal, mechanical, and chemical sti [16]. Among these, mechanical stimulation has drawn little attention in pain there are few pieces of research about the evaluation of the human sensation u tactile sensing of electrical skin. The aforementioned issues motivated us to s superficial acute pain sensation due to mechanical stimuli. In this paper, firstly chanical model of compression on the skin is investigated and the theoretical s tribution is obtained based on the Flament solution. Then, setting the stress as the pain sensation is evaluated based on the modified Hodgkin-Huxley mode gate control theory. The theoretical stress distribution is verified by finite eleme sis (FEA). The influences of compression magnitude and nociceptor location ar vestigated. Although pain has been studied from the molecular level to the level of the entire nervous system for a long period, the computation models of pain sensation are still limited [16,17]. Researchers have developed mathematical models at different levels: molecular [18,19], cellular [20], and neuron network levels [11,21,22]. Britton et al. [20] established a mathematical model for GCT and used it to explain the quality of pain. Xu et al. [5,8,16] have systematically studied the thermal pain sensation. By combining the transduction, transmission, modulation and perception process, the temperature, thermal stress, and chemical burn are all taken into consideration in thermal pain analysis. Based on Xu's theory, Yin et al. [9] studied the skin pain sensation under the heating of epidermal electronic devices.
There are three main pain stimuli: thermal, mechanical, and chemical stimulation [16]. Among these, mechanical stimulation has drawn little attention in pain analysis; there are few pieces of research about the evaluation of the human sensation under the tactile sensing of electrical skin. The aforementioned issues motivated us to study the superficial acute pain sensation due to mechanical stimuli. In this paper, firstly, the mechanical model of compression on the skin is investigated and the theoretical stress distribution is obtained based on the Flament solution. Then, setting the stress as an input, the pain sensation is evaluated based on the modified Hodgkin-Huxley model and the gate control theory. The theoretical stress distribution is verified by finite element analysis (FEA). The influences of compression magnitude and nociceptor location are also investigated.

Modeling and Analysis
There are three main types of pain: nociceptive pain, inflammatory pain, and neuropathic pain [16]. Nociceptive pain has been studied from different aspects before. Figure 1 is a schematic of the pain sensation process, including (1) transduction: the different kinds of outside stimuli are converted into nerve impulses by receptors (nociceptors) in the skin; (2) transmission: the nerve impulses, which carry the information of stimuli, transmit to the dorsal horn through nerve fibers; (3) modulation: downward inhibition and facilitation of nociceptive transmission; (4) perception: the evaluation of signals received in higher-order structures of the nervous system. Thus, a mathematical model of compression pain is established in this paper, which consists of three parts: the mechanical model of skin compression, the model of transduction, and the model of modulation and perception. The time delay caused by the transmission is omitted in this paper. In the following parts of this section, the details of the model will be demonstrated.

Mechanical Model of Skin Compression
Similar to thermal pain, which has been studied by Yin [9] before, compression pain is determined by stress at the location of the nociceptors, instead of the surface of the skin. As compression is applied on the surface of the skin (Figure 2), the sensation of touching will first be generated. With the growth of the compression amplitude, the stress on the tissue will increase and reach a critical value [9], which may lead to the sensation of pain. Mechanical pain is not only an unpleasant feeling but also potentially damaging to the skin and nervous system [10]. Thus, the mechanical model of skin is developed in this part, where the skin is regarded as a two-dimensional semi-infinite solid for simplicity, as demonstrated in Figure 2. A uniform compression with amplitude q and length 2a is applied on the top surface of the skin. To obtain the stress field of skin, the Flament answer to the semi-infinite plane strain problem is utilized here [23]:

Mechanical Model of Skin Compression
Similar to thermal pain, which has been studied by Yin [9] before, compress is determined by stress at the location of the nociceptors, instead of the surfac skin. As compression is applied on the surface of the skin (Figure 2), the sens touching will first be generated. With the growth of the compression amplitu stress on the tissue will increase and reach a critical value [9], which may lea sensation of pain. Mechanical pain is not only an unpleasant feeling but also po damaging to the skin and nervous system [10]. Thus, the mechanical model o developed in this part, where the skin is regarded as a two-dimensional semi solid for simplicity, as demonstrated in Figure 2. A uniform compression with am q and length 2a is applied on the top surface of the skin. To obtain the stress field the Flament answer to the semi-infinite plane strain problem is utilized here [23]: The mechanical model of compression on the skin.
The strain field of skin can be acquired by integrating ξ in Equation (1) from When y belongs to different intervals: (−a, a) or (a, ∞)∪(−∞, −a), the integration res distinct, noting that it is compression pain that is being investigated, and y stand position of the nociceptor on the skin [19]. Hence, it is reasonable to assume that ciceptor locates at the x-axis in compression pain, for the convenience of the a leading the integration in (1)−(3) to: The strain field of skin can be acquired by integrating ξ in Equation (1) from −a to a. When y belongs to different intervals: (−a, a) or (a, ∞)∪(−∞, −a), the integration results are distinct, noting that it is compression pain that is being investigated, and y stands for the position of the nociceptor on the skin [19]. Hence, it is reasonable to assume that the nociceptor locates at the x-axis in compression pain, for the convenience of the analysis, leading the integration in Equations (1)-(3) to: Since there is no shear stress, σ x and σ y are the two principal stresses at the location of the nociceptors with σ x < σ y < 0 where σ x is the minor principal stress and σ y is the major principal stress on the x-axis.

Model of Transduction Current and Frequency Modulation
When the skin is stimulated by the compression, the nociceptors trigger potential action as a result of ion transportation across the cell membrane of the neuron [19]. The membrane potential of the nerve excitation can be described by the revised Hodgkin-Huxley (H-H) model [18,24,25], as illustrated in Figure 3. The total current consists of the transportation of ions across the corresponding ion channel and the current that charges the membrane capacity, which gives: where V m (mV) stands for the membrane potential of the nociceptors, and t (ms) represents the time. C m = 1 µF/cm 2 is the membrane capacitance per unit area [9]. I Na , I K , I A , and I L correspond to the current induced by the sodium ions (Na + ), potassium ions (K + ), fast transient K + , and leakage current component. I mech is the current induced by the mechanical stimuli, which is revealed to be a function of the stress at the location of the nociceptor: When the skin is stimulated by the compression, the nociceptors trigger po action as a result of ion transportation across the cell membrane of the neuron [1 membrane potential of the nerve excitation can be described by the revised Ho Huxley (H-H) model [18,24,25], as illustrated in Figure 3. The total current consist transportation of ions across the corresponding ion channel and the current that c the membrane capacity, which gives: a   C  I  V  t  I  I  I  I  I  +  −  +  +  =  + where Vm (mV) stands for the membrane potential of the nociceptors, and t (ms) sents the time. Cm = 1 μF/cm 2 is the membrane capacitance per unit area [9]. INa, IK, IL correspond to the current induced by the sodium ions (Na + ), potassium ions (K transient K + , and leakage current component. Imech is the current induced by the m ical stimuli, which is revealed to be a function of the stress at the location of th ceptor: Here, σt = 20 kPa [26] is the mechanical threshold and H(x) is the Heaviside fu Cm1 = 2 μA/cm 2 , Cm2 = 2, Cm3 = −1 μA/cm 2 are the mechanical stimuli-related co and Ishift = 8.1 mA is the current to guarantee the action potential when σ > σt [27].  Here, σ t = 20 kPa [26] is the mechanical threshold and H(x) is the Heaviside function. C m1 = 2 µA/cm 2 , C m2 = 2, C m3 = −1 µA/cm 2 are the mechanical stimuli-related constants and I shift = 8.1 mA is the current to guarantee the action potential when σ > σ t [27].
The ionic current I Na , I K , I A , and I L driven by the membrane potential could be described as: where E Na = 55 mV, E K = −72 mV, E A = −75 mV, and E L = −17.5 mV are the reversal potentials as demonstrated in Figure 3. κ i represents different ionic conductance here and κ Na = 120 mS/cm 2 , κ K = 20 mS/cm 2 , κ A = 47.7 mS/cm 2 , and κ L = mS/cm 2 [9]. The m, n, h, A, and B are the gating variables satisfy: where x represents the gating variables. For x = m, h, n, the expressions of τ x and x ∞ are given as follows: V m can be solved by substituting Equations (9)-(20) into Equation (8). Instead of the amplitude, the frequency of V m represents the intensity of mechanical stimuli, which is utilized in the analysis of modulation and perception in the following section.

Model of Modulation and Perception
There are various theories explaining the relationship between neural excitation and pain sensation. Among them, the gate control theory, which is proposed by Melzack [15], is one of the most successful theories, which precisely describes the modulation and perception procedure of skin pain sensation. Britton et al. [20] give the mathematical model of GCT, which is schematically illustrated in Figure 4. The mathematical descriptions are given as: where the V i and V e are the potentials of the inhibitory substantia gelatinosa (SG) and excitatory SG cells, and V t and V b stand for the potential of central transmission cell (T-cell) and midbrain. The subscript l and s represent the large fibers (C, Aδ) and small fibers (Aβ), respectively. θ li and θ se represent the proportion of excitation transmitted to inhibitory and excitatory SG cells through large and small nerve fibers. Both θ li and θ se are taken as 0.8 in this paper. x i denotes the frequency of the membrane potential transmitted on the fiber, and the definition of function f j (V j ) is Here, the subscript j represents i, e, t, b in Equations (21)- (24). V j0 = −70 mV is the initial membrane potential, and V thr = −55 mV is the threshold potential value for the pain sensation [9]. The output of T-cell V t is in direct relation to the pain sensation because the noxious signal is transmitted to the cortex when V t exceeds the threshold value. Note, it has been demonstrated that intense noxious stimuli information is transmitted through small fibers while the slight stimuli information is carried by large fibers. Since compression pain is analyzed in this paper, the x l and x s are taken to be zero and the frequency of V m , respectively [9].
Here, the subscript j represents i, e, t, b in Equations (21)− (24). Vj0 = −70 itial membrane potential, and Vthr = −55 mV is the threshold potential valu sensation [9]. The output of T-cell Vt is in direct relation to the pain sensatio noxious signal is transmitted to the cortex when Vt exceeds the threshold v has been demonstrated that intense noxious stimuli information is transm small fibers while the slight stimuli information is carried by large fiber pression pain is analyzed in this paper, the xl and xs are taken to be zero quency of Vm, respectively [9].

Stress Distributions in the Skin
The strain distributions of the skin in Equations (4)−(6) were verified ment analysis. A two-dimensional rectangular solid with 10 m × 5 m was ABAQUS software, which was large enough to simulate semi-infinity. A p kPa with radius a = 10 mm was applied on the top surface of the skin and faces of the skin were fixed. Young's modulus and Poisson's ratio were se kPa and ν = 0.3 [6,28]. The six-node triangular CPS6M element was chosen

Stress Distributions in the Skin
The strain distributions of the skin in Equations (4)-(6) were verified by finite element analysis. A two-dimensional rectangular solid with 10 m × 5 m was established in ABAQUS software, which was large enough to simulate semi-infinity. A pressure q = 20 kPa with radius a = 10 mm was applied on the top surface of the skin and the other surfaces of the skin were fixed. Young's modulus and Poisson's ratio were set to be E = 20 kPa and ν = 0.3 [6,28]. The six-node triangular CPS6M element was chosen to discretize the model. The size of the elements was controlled in the range from 0.005 m to 0.2 m with a fine mesh around the compressed region. The number of the elements was 3524, which was large enough to guarantee the convergence of the FEA model. Figure 5a demonstrates the Mises stress distribution obtained from the FEA results. It can be observed that the Mises stress reached the maximum around the compressed region and dramatically decreased to zero when moving away. Figure 5b shows the comparison of σ x (blue) and σ y (red) along the red path in Figure 5a. The solid lines and the dots represent the theoretical solutions and FEA results, respectively. The perfect agreement verifies the theoretical solutions. Besides, it shows that the absolute value of σ x is greater than σ y . Thus, the minor principal stress σ x was utilized in the calculation of membrane potential.
Micromachines 2022, 13, x FOR PEER REVIEW 7 of 10 the model. The size of the elements was controlled in the range from 0.005 m to 0.2 m with a fine mesh around the compressed region. The number of the elements was 3524, which was large enough to guarantee the convergence of the FEA model. Figure 5a demonstrates the Mises stress distribution obtained from the FEA results. It can be observed that the Mises stress reached the maximum around the compressed region and dramatically decreased to zero when moving away. Figure 5b shows the comparison of σx (blue) and σy (red) along the red path in Figure 5a. The solid lines and the dots represent the theoretical solutions and FEA results, respectively. The perfect agreement verifies the theoretical solutions. Besides, it shows that the absolute value of σx is greater than σy. Thus, the minor principal stress σx was utilized in the calculation of membrane potential.

Influence of the Compression Amplitude
Obviously, the level of compression pain is related to the pressure amplitude. Thus, the compression skin pain sensation with different amplitudes is studied in this section, as illustrated in Figure 6. The compression radius was fixed at a = 10 mm and the location

Influence of the Compression Amplitude
Obviously, the level of compression pain is related to the pressure amplitude. Thus, the compression skin pain sensation with different amplitudes is studied in this section, as illustrated in Figure 6. The compression radius was fixed at a = 10 mm and the location of the nociceptor was assumed to be at x = 1.6 mm, y = 0, corresponding with the approximate epidermis thickness of the skin [1]. The nonlinear differential equations of the revised H-H model and the mathematical GCT were solved based on the four-order Runge-Kutta method. Figure 6a demonstrates the membrane potential V m under different compression amplitude values (q = 15 kPa, q = 25 kPa, q = 35 kPa). It was demonstrated that the membrane potential V m showed a quasi-periodical variation with time and the stimuli intensity increased with the growing compression amplitude. Figure 6b demonstrates the frequency of V m was positively related to compression amplitude q after surpassing the mechanical threshold σ t = 20 kPa. The membrane frequency here was determined by fast-Fourier transformation (FFT). The output of T-cell V t under different compression values is illustrated in Figure 6c. It was found that V t increased and reached a plateau finally. The outputs of T-cells caused by q = 25 kPa and q = 35 kPa both exceeded the threshold value (−55 mV) for pain sensation.

Influence of the Nociceptor Location
It can be observed from Figure 5b that the stress on the skin varied w Thus, it was necessary to study the influence of the location of the nocicept pain sensation, as illustrated in Figure 7. The nociceptor was set to locate with varying depths to the surface of the skin (1 mm, 10 mm, and 50 mm). T the compression amplitude were taken as a = 10 mm and q = 25 kPa. It was s frequency of Vm decreased with increasing depth, shown in Figure 7a with th mm), the red line (10 mm), and the blue line (50 mm), respectively. The frequ stable when the depth continued to increase, corresponding to a stress lev (Figure 7b). The output of T-cell Vt also decreased with the increasing dept ceptor as shown in Figure 7c. With the same compression amplitude of q

Influence of the Nociceptor Location
It can be observed from Figure 5b that the stress on the skin varied with the depth. Thus, it was necessary to study the influence of the location of the nociceptor on the skin pain sensation, as illustrated in Figure 7. The nociceptor was set to locate on the y-axis with varying depths to the surface of the skin (1 mm, 10 mm, and 50 mm). The radius and the compression amplitude were taken as a = 10 mm and q = 25 kPa. It was shown that the frequency of V m decreased with increasing depth, shown in Figure 7a with the gray line (1 mm), the red line (10 mm), and the blue line (50 mm), respectively. The frequency became stable when the depth continued to increase, corresponding to a stress level less than σ t (Figure 7b). The output of T-cell V t also decreased with the increasing depth of the nociceptor as shown in Figure 7c. With the same compression amplitude of q = 25 kPa, the output of T-cell V t became less than the threshold value when the depth reached 50 mm compared with the other shallower locations (10 mm and 50 mm). In other words, skin with a thinner epidermis is more likely to experience pain sensations. According to our life experience, our hands are indeed more sensitive to pain until the skin grows thick calluses.
( Figure 7b). The output of T-cell Vt also decreased with the increasing dep ceptor as shown in Figure 7c. With the same compression amplitude of q output of T-cell Vt became less than the threshold value when the depth re compared with the other shallower locations (10 mm and 50 mm). In othe with a thinner epidermis is more likely to experience pain sensations. Acc life experience, our hands are indeed more sensitive to pain until the ski calluses.

Conclusions
In this paper, a model for compression pain sensation of skin is developed, which includes: (1) a model of skin compression, where the stress distribution on the skin is obtained theoretically based on the Flament solution; (2) a model of transduction, where the electrical signal converted from external mechanical stimuli is analyzed with the revised H-H model; (3) a model of modulation and perception, where the correlation between the nerve impulse and the pain sensation is revealed through the mathematical model of GCT theory. The skin compression model was verified by the perfect agreement between theoretical stress distribution and FEA results. Factors that influence pain sensation were also investigated, including compression amplitude and location of the nociceptor. However, there are some limitations to the model developed in this paper, which still need further investigation. (1) The skin was modeled as a semi-infinite elastic solid in this paper with static loading on the surface, although it has been demonstrated that the viscoelasticity of skin has to be considered in dynamical analysis. (2) The mechanical stimuli compression was uniform in this study, and different kinds of mechanical loadings need to be investigated to account for the complexity of the real situation. (3) The parameters in the model of transduction and the model of modulation and perception have not been compared with experimental results. It is hoped that the parameters could be determined by comparing the mathematical pain sensation model with experimental results.