Modeling of skin thermal pain: A preliminary study Applied Mathematics and Computation 205 (2008) 37–46 Contents lists available at ScienceDirect Applied Mathematics and Computation journal homepage: www.elsevier .com/ locate/amc Modeling of skin thermal pain: A preliminary study Feng Xu a, Ting Wen a, Keith Seffen a, Tianjian Lu b,* a Engineering Department, Cambridge University, Cambridge CB2 1PZ, UK b MOE Key Laboratory of Strength and Vibration, School of Aerospace, Xi’an Jiaotong University, Xi’an 710049, PR China a r t i c l e i n f o a b s t r a c t Keywords: Skin tissue Thermal pain Holistic model Thermomechanical model 0096-3003/$ - see front matter � 2008 Elsevier Inc doi:10.1016/j.amc.2008.05.045 * Corresponding author. E-mail address: tjlu@mail.xjtu.edu.cn (T. Lu). Skin thermal pain is one of the most common problems for humans in everyday life. The understanding of the underlying physical mechanisms is still not clear, and modeling of them is very limited. In this paper, a holistic mathematical model for quantifying skin ther- mal pain sensation is developed. The model is composed of three interconnected parts: peripheral modulation of noxious stimuli, which converts the energy from a noxious ther- mal stimulus into electrical energy via nerve impulses; transmission, which transports these neural signals from the site of transduction in the skin to the spinal cord and brain; and modulation and perception in the spinal cord and brain. With this model, a direct rela- tionship is built between the level of thermal pain sensation and the character of noxious stimuli. � 2008 Elsevier Inc. All rights reserved. 1. Introduction Skin is the largest single organ of the body. It consists of several layers and plays a variety of important roles including sensory, thermoregulatory and host defense etc. However, in an extreme environment, an uncomfortable feeling or pain sen- sation is evoked due to heat or cold: the skin fails to protect the body when the temperature lies outside of the normal phys- iological range. In medicine, with advances in laser, microwave and similar technologies, various thermal therapeutic methods have been developed recently and widely used to cure disease/injury involving skin tissue; but the corresponding problem of pain relief has limited the further application and development of these thermal treatments. As one of the most important sensations, pain has been long studied over a range of scales, from the molecular level (ion channel) to the entire human neural system level. Thermal stimulation, as one of the three main pain stimulators (thermal, mechanical and chemical), has been widely used in the study of pain, as in the examination of tissue injury and sensitisation mechanisms, and in the quantification of therapeutic effects due to pharmacological, physical, and psychological interven- tions. However, the understanding of the underlying mechanisms is still far from clear. One reason is that pain is influenced by many factors, including physiological, such as stimulus intensity and duration, and psychological factors such as attention and empathy. The utilization of computational models in the field of pain has been very limited, with previous attempts at pain mod- eling focusing mainly on acute pain. There are nonetheless strong arguments for mathematical modeling [1,2], it can handle extremely complex theories and can be used to predict behaviors which have perhaps previously gone unnoticed. The meth- od is also non-invasive. Several mathematical models have been developed at different levels: at the molecular level and cel- lular level [1,3–9]; and at the level of neuron network [10,11]. None of these models have correlated the external stimulus parameters directly with the pain sensation level, and no transmission process has been considered. For example, the ‘‘gate control theory” has been used to extrapolate the relevant features of pain and has been translated into a mathematical model . All rights reserved. mailto:tjlu@mail.xjtu.edu.cn http://www.sciencedirect.com/science/journal/00963003 http://www.elsevier.com/locate/amc 38 F. Xu et al. / Applied Mathematics and Computation 205 (2008) 37–46 by Britton and Skevington [1,7–9], but it only considers the perception and modulation processes of pain sensation at spinal cord and brain, without considering the transduction process occurring within the skin tissue. These issues are addressed in this study. The focus of this paper is placed upon the superficial nociceptive acute pain. Neuropathic pain and chronic pain as well as psychological factors that influence the pain are not modeled at present stage. In the following sections, the mechanisms of skin thermal pain are first described with emphasis on the transduction function of skin nociceptors; the development of the model is presented next; finally, a simple case of surface heating is used for the preliminary analysis of thermal pain. 2. Skin thermal pain The physiology of pain has been studied extensively and a detailed description of the current understanding on this sub- ject can be found in [12–15]. The pain pathway is schematically shown in Fig. 1 and can be simply described as: (1) trans- duction: when a stimulus is applied to the skin, the nociceptors located there will trigger action potentials by converting the physical energy from a noxious thermal, mechanical, or chemical stimulus into electrochemical energy; (2) transmission: the signals are then transmitted in the form of action potential chains (similar to pulse trains) via nerve fibers from the site of transduction (periphery) to the dorsal horn of the spinal cord, which then activates the associated transmission neuron; (3) perception: the appreciation of signals arriving in higher structures as pain; (4) modulation: descending inhibitory and facil- itory input from the brain that influences (modulates) nociceptive transmission at the level of the spinal cord. 2.1. Nociceptors in skin tissue Nociceptors, the special receptor for pain sensation, are the first cells in the series of neurons leading to the sensation of nociceptive pain. Nociceptors transduct mechanical, chemical and/or thermal energy to ionic current (noxious stimuli into depolarizations that trigger action potentials), conduct the action potentials from the peripheral sensory site to the synapse in the central nervous system, and convert the action potentials into neurotransmitter release at the presynaptic terminal (frequency modulation) [16]. Nociceptors in skin are generally one of three kinds of peripheral nerves: myelinated afferent Ad fibers and Aa fibers, as well as unmyelinated C afferent fibers, while thermal pain sensations are mediated mainly by both thin myelinated Ad- and unmyelinated C-fibers [17]. Many unmyelinated C-fibers can be traced far into the epidermal layer. The study of myelinated mechanical nociceptor endings in cat hairy skin showed that the free nerve endings of pain recep- tors exist around the depth of 50 lm [18], while the receptor depth in the hairy skin of monkeys estimated from responses to ramped stimuli ranged from 20 to 570 lm with a mean of 201 lm [19,20]. Fig. 1. Schematic of the holistic skin thermal pain model. F. Xu et al. / Applied Mathematics and Computation 205 (2008) 37–46 39 2.2. Ion channels of nociceptors All the essential functions of nociceptors depend on ion channels [12,16], including heat activated channels, capsaicin receptors, ATP-gated channels, proton-gated channels or acid-sensing channels, nociceptor-specific voltage-gated Na+ chan- nels, mechanosensitive channels, etc. In spite of these different channels, they are generally gated by three different meth- ods, namely, thermal (hot or cold), mechanical, and chemical, with thresholds of �43 �C and �0.2 MPa for the first two, respectively. When a noxious stimulus is applied to a nociceptor (larger than the threshold), the corresponding ion channels will be opened, resulting in a current, thereby triggering the action potential. 3. Model development According to the above description of pain pathway, a holistic model of skin thermal pain is developed in this section that is composed of three sub-models: model of transduction, model of transmission, and model of perception and modulation. The schematic of the holistic model is shown in Fig. 1. The details of each sub-model are given below. 3.1. Model of transduction 3.1.1. Thermomechanical model of skin tissue It is now accepted that thermal pain is determined by the temperature at the location of nociceptors, not at the skin sur- face. When the skin is heated to a temperature beyond the threshold (�43 �C), thermal damage of skin tissue begins accu- mulate. The thermal damage causes cells to break down, releasing a number of tissue byproducts and mediators, which will activate and sensitize the nociceptors. Furthermore, thermally induced stresses due to non-uniform temperature distribu- tions may also lead to the sensation of thermal pain, in addition to the pain generated by heating alone. Supporting evidence shows that, for the same level of nociceptor activity, a heat stimulus evokes more pain than a mechanical stimulus, and that tissue deformation due to heating and cooling may explain the origins of pain [21,22]. A thermomechanical model of skin tissue has been developed in our previous study [23]. The skin is treated as a layered – ‘‘laminated” – material according to its real structure, whose overall properties are assembled in a composite manner, as shown in Fig. 2. Its thermomechanical behavior is simplified as a ‘sequentially coupled’ problem, in other words, the mechanical and thermal properties are inde- pendent of each other. Solutions of the temperature, thermal damage and thermal stress fields are given next. 3.1.1.1. Temperature profile. The temperature field can be obtained by solving the bioheat transfer model of skin tissue. For example, for a one-layer skin model under surface contact heating at temperature Ts, the closed-form solution of tempera- ture at the location of nociceptor, Tn, can be obtained as [23]: Tnðzn; tÞ ¼ T0ðznÞ þ 2a H Ts � k dT0ðznÞ dzn ���� z¼0 � �X1 m¼1 bm sinðbmðzn þ HÞÞ 1 ab2m þ -bqbcb qc 1� e�ab 2 mt� -bqbcb qc t � � ; ð1Þ where T0ðzÞ is initial temperature field in the tissue with z ¼ 0 corresponding to skin surface; q, c, k are the density, specific heat and thermal conductivity of skin tissue, respectively; qb, cb are the density and specific heat of blood, -b is the blood Fig. 2. Skin structure and corresponding idealized skin model. 40 F. Xu et al. / Applied Mathematics and Computation 205 (2008) 37–46 perfusion rate; a ¼ k=qc is thermal diffusivity of skin tissue; H is total thickness of skin tissue; and bm ¼ mp=H;m ¼ 1;2;3; . . . 3.1.1.2. Thermal damage Z t � � X ¼ 0 A exp � Ea RT dt; ð2Þ DegðtÞ ¼ 1� expð�XðtÞÞ: ð3Þ X is the dimensionless indicator of damage, Deg is the degree of thermal damage (Deg ¼ 0 equals no damage while Deg ¼ 1 equals fully damaged), A is a material parameter (frequency factor), Ea is the activation energy, and R ¼ 8:314 J=mol K is the universal gas constant. 3.1.1.3. Thermal stress. rk ¼ Ek ��kkDT þ C1ð1þ mkÞ PM i¼1 R zi zi�1 Ei�kiDT dz � � þ C2ð1þ mkÞ PM i¼1 R zi zi�1 Ei�kiDTzdz � �� � þzð1þ mkÞ C2 PM i¼1 R zi zi�1 Ei�kiDT dz � � þ C3ð1þ mkÞ PM i¼1 R zi zi�1 Ei�kiDTzdz � �� � 8>>>< >>>: 9>>>= >>>; : ð4Þ r is the in-plane stress; subscript ‘‘k” denotes the kth layer of skin; E ¼ E=ð1� m2Þ, �k ¼ ð1þ mÞk, where E is the Young’s mod- ulus, m is the Poisson’s ratio, k is thermal expansion coefficient; C1, C2, C3 are constants depending on the relative thickness of each layer of skin tissue. 3.1.2. Current As previously discussed, the pain signal starts from the current induced by the opening of ion channels in nociceptors. Since ion channels are generally gated by three different methods (thermal, mechanical and chemical), there are correspond- ingly three different currents. The total current can be calculated as: Ist ¼ Iheat þ Ichem þ Imech: ð5Þ The heat current, Iheat, is assumed to be a function of nociceptor temperature, Tn, and thermal pain threshold, Tt , given as: Iheat ¼ fhðTn; TtÞ: ð6Þ Here, Tt is assumed to be 43 �C [24,25]. The chemical current, Ichem, is assumed to depend on the thermal damage degree of skin tissue, Deg: Ichem ¼ fcðDegÞ: ð7Þ The mechanical current, Imech, is assumed to be a function of the stress at the location of nociceptor, rn, and mechanical pain threshold, rt: Imech ¼ fmðrn;rtÞ; ð8Þ where rt is assumed to be 0.2 MPa [26]. 3.1.3. Frequency modulation When the evocated current overpasses the threshold, an action potential is generated. It is now accepted that the inten- sity of external stimulation is carried through the frequency of these impulses, fs, not the exact magnitude or shape of the signal, which can be calculated as: fs ¼ ffmðIstÞ: ð9Þ Unfortunately, there exists no study to quantify the current–frequency relation. We need, therefore, to choose a model for the generation of an action potential. Although there has been no relative analysis on nociceptor kinetics, all neurons have been found to behave qualitatively similar to that described by the Hodgkin–Huxley (H–H) model of nerve excitation [27]. Therefore, the H–H model is adopted here to model the frequency response of skin nociceptors. Mathematically, the H–H model can be described as: Cmem dVmem dt ¼ Ist þ gNam 3hðENa � VmemÞ þ gKn 4ðEKNa � VmemÞ þ gLðEL � VmemÞ; ð10Þ where Vmem is the membrane potential (mV), which is positive when the membrane is depolarized and negative when the membrane is hyperpolarized; ENa, EK and EL are the corresponding reversal potentials for sodium, potassium and leakage current components, given by the Nernst equation (all in mV), respectively; gNa, gK and gL are the maximal ionic conductance through sodium, potassium and leakage current components (all in mS/cm2), respectively; t is time (ms); Cmem is membrane capacity per unit area (lF/cm2); Ist is the stimulation-induced current density, as defined by Eq. (5), which is positive F. Xu et al. / Applied Mathematics and Computation 205 (2008) 37–46 41 when the current is outward (lA/cm2); INa, IK and IL are the sodium, potassium and leakage current components (lA/cm 2), respectively. Each of the three ionic currents (INa, IK and IL) is driven by an electrical potential difference and a permeability coefficient, or conductance [27]. The conductances of the ionic currents are regulated by voltage dependent activation and inactivation variables, known as gating variables (m, n and h), which are given as [27]: dx dt ¼ axð1� xÞ � bxx: ð11Þ Here, x can be either m, n and h, which are voltage- and time-dependent parameters; ax and bx are rate constants (in s �1). ax and bx can be approximated from experiments [27] as: an ¼ �0:01ðVmem þ 50Þ exp½�ðVmem þ 50Þ=10� � 1 ; bn ¼ 0:125 exp �ðVmem þ 60Þ 80 � � ; ð12Þ am ¼ �0:1ðVmem þ 35Þ exp½�ðVmem þ 35Þ=10� � 1 ; bm ¼ 4 exp �ðVmem þ 60Þ 18 � � ; ð13Þ ah ¼ 0:07 exp �ðVmem þ 60Þ 20 � � ; bh ¼ 1 exp½�ðVmem þ 30Þ=10� þ 1: ð14Þ Here, the various constants are obtained empirically from experimental data [27], given as: Cm ¼ 1:0 lF=cm2, ENa ¼ 55 mV, EK ¼ �72 mV, xfac ¼ 1, gNa ¼ 120 mS=cm2, gK ¼ 36 mS=cm2, gL ¼ 0:3 mS=cm2. 3.2. Model of transmission This model simulates the transmission of noxious stimulus triggered neural signals from the skin along the respective fibers to spinal cord and brain. Since the stimulation intensity is carried through the firing rate, or frequency, of these im- pulses, the actual shape, amplitude and duration of a single spike are not taken into account. The time for this transmission, tt , is obtained once the conduction velocity and the corresponding nerve length are determined. 3.2.1. Nerve length ðLnÞ For simplicity, the present study only considers the case that the skin is subjected to a thermal loading restricted to a small area (i.e., localized stimuli). Consequently, all fibers connected with the skin nociceptors are postulated to have the same length. Here, a length of 1 m is chosen, approximately the distance between a finger and the spinal cord [28]: Ln ¼ 1 m: ð15Þ 3.2.2. Conduction velocity ðvcÞ The conduction velocity ðvcÞ is found to be directly related to fiber diameter ðDÞ [13]. A linear relationship between vc and D is assumed: vc ¼ vcN ¼ cdD; ð16Þ where vcN is the conduction velocity under normal condition and cd is the coefficient between diameter and velocity. In addition to fiber diameter, temperature also influences the nerve conduction velocity. For example, it has been found that the conduction velocity is reduced in a cold environment and the rate of conduction velocity reduction is identical for slow and fast fibers [29,30]. Thus, the conduction velocity, considering temperature effect, can be calculated as: vc ¼ vcT ¼ cTvcN; ð17Þ where vcT is the velocity when the fiber is in an environment of temperature T, cT is the influence factor of temperature on the conduction velocity: cT is given by Paintal [29,30], which was obtained from experiments using cat. Here it is assumed that the resulting rules are also applicable to human peripheral nerves in view of their structural and functional similarities. 3.2.3. Transmission time ðttÞ Once the nerve length and conduction velocity are known, the transmission time can be calculated as: tt ¼ Ln=vc: ð18Þ 3.3. Model of modulation and perception The theory of gate control [31] is employed here to describe the modulation and perception process of skin thermal pain. A schematic description of the theory is presented in Fig. 3. In general, it is the small fibers ðc;AdÞ that carry information about noxious stimuli whilst the role of large fibers ðAbÞ is to carry information about less intense mechanical stimuli. As the signal from the ðc;AdÞ fibers is routed through substantia gelatinosa (SG) to central transmission (T) cells and onwards, the double inhibition (indicated by the minus signs) actually strengthens the signal. The signal from the Ab fibers, however, is diminished in strength when routed through SG. Fig. 3. Schematic of the model of modulation and perception. 42 F. Xu et al. / Applied Mathematics and Computation 205 (2008) 37–46 Gate control theory shown in Fig. 3 has been translated into a mathematical model by Britton and Skevington [1,7–9], given as: Table 1 Thermo Parame Blood s Blood s Arteria Core te si _V i ¼ �ðV i � V i0Þ þ gliðxlÞ þ gmiðxmÞ ð19Þ 0:7 _V i ¼ �ðV i þ 70Þ þ 60 tanhðhlixlÞ þ 40 tanhðfmðVmÞÞ ð20Þ se _Ve ¼ �ðVe � Ve0Þ þ gseðxs;VeÞ ð21Þ 0:7 _Ve ¼ �ðVe þ 70Þ þ 40 tanhðhsexsÞ½1þ 3 tanhð4f eðVeÞÞ� ð22Þ where subscripts i, e, t and m stand for inhibitory SG cell, excitory SG cell, T-cell and midbrain, respectively; sj is time con- stant, Vj is membrane potential; Vj0 is initial membrane potential; xj is firing frequency; xl and xs are signals (frequency) from large and small fibers, respectively; functions gjk represent the effects of inputs ðjÞ to a cell ðkÞ on its steady state slow po- tential; tanhð4f eðVeÞÞ is the NMDA (N-methyl-D-aspartic acid) component of the equation. The NMDA receptor has been ar- gued to be responsible for phenomena in pain sensation such as wind-up, where the response of a neuron increases progressively as the neuron is repeatedly stimulated. hli and hse are the proportions of the inputs that pass through interneu- rons in the SG and correspondingly ð1� hliÞ and ð1� hseÞ are the proportions passing through to the T-cell. The output from the T-cell, Vt , is taken to be directly related to the pain level, such that if the T-cell exceeds its firing threshold of �55 mV then the noxious signal is transmitted to the next relay point. If the noxious signals reach the cortex, they are perceived as pain [1]. 4. Results and discussion: case study For illustration, the applicability of the developed holistic thermal pain model is illustrated using a simple surface heating case. The skin surface initially at normal temperature of 37 �C is suddenly taken into contact with a hot source of constant temperature (50 �C) for 20 s. Using the previous skin thermomechanical model, the temperature history of nociceptor is ob- tained first, which is then used as the input for the neural model. The nociceptors are assumed to be C fibers with a conduc- tion velocity of 1 m/s. The one-dimensional, three-layer skin model is used where effect of blood perfusion is only considered in the dermis layer. The relevant parameters of skin tissue used are summarized in Tables 1 and 2. 4.1. Validation of the model The predicted membrane potential and frequency responses under stimulus of different nociceptor temperatures ðTnÞ are presented in Fig. 4b. Compared with experimental observations of nociceptors [25] shown in Fig. 4a, a good agreement has physical properties of blood ters References ensity (kg/m3) 1060.0 [42] pecific heat (J/kg K) 3770.0 [43] l blood temperature (�C) 37 mperature (�C) 37 Table 2 Thermophysical properties of skin tissue for three-layer model Parameters Value References Thermal expansion coefficient (�10�4/K) Epidermis 1 Assumption Dermis 1 Assumption Subcutaneous fat 1 Assumption Poisson’s ratio (–) Epidermis 0.48 [44] Dermis 0.48 [44] Subcutaneous fat 0.48 [44] Young’s modulus (MPa) Epidermis 102 [45] Dermis 10.2 [45] Subcutaneous fat 0.0102 [45] Skin density (kg/m3) Epidermis 1190.0 [42] Dermis 1116.0 [42] Subcutaneous fat 971.0 [42] Skin thermal conductivity (W/mK) Epidermis 0.235 [46] Dermis 0.445 [46] Subcutaneous fat 0.185 [46] Skin specific heat (J/kg K) Epidermis 3600.0 [47] Dermis 3300.0 [47] Subcutaneous fat 2700.0 [47] Metabolic heat generation (W/m3) Epidermis 368.1 [48] Dermis 368.1 [48] Subcutaneous fat 368.3 [48] Thickness (mm) Epidermis 0.1 [49] Dermis 1.5 [50] Subcutaneous fat 4.4 Assumption -80 -40 0 40 0.0 4.03.02.01.0 T n = 45 O C V m ( m V ) t (s) -80 -40 0 40 Tn = 51 O C a b Fig. 4. (a) Experimental observations of the responses of C nociceptors in mouse glabrous skin evoked by heat stimuli of 45 and 51 �C (response threshold to heat was 43 �C) [25]; (b) predicted membrane potential from our model. F. Xu et al. / Applied Mathematics and Computation 205 (2008) 37–46 43 been achieved for the intensity–response relationship although. In particular, with the model, the predicted neural impulse rate is comparable to that of the actual nociceptors. 4.2. Influence of nociceptor location The nociceptors at different depths (20, 50, 100, 200 lm) are considered. The predicted temperature history of the noci- ceptor is shown in Fig. 5a, which is then used as the input for the pain model. The corresponding pain level is shown in Fig. 5b. It is seen from Fig. 5 that, at the same stimulus intensity, the pain level is higher if the nociceptor is located closer to the surface of skin. This may be used to explain why different pain thresholds were obtained by different studies for the same stimulus [19,20]. 0 5 10 15 20 39 42 45 48 51 T s = 50 O C z n (μm) 20 50 100 200 T n O C t (s) Threshold of nociceptor 43 O C 0 5 10 15 20 -70 -60 -50 -40 -30 -20 -10 0 V t ( m v) t (s) Ts= 50 O C z n (μm) 20 50 100 200 Threshold of pain -55 mV (a) Temperature history (b) Pain level Fig. 5. Influence of nociceptor location on (a) temperature and (b) perceived pain level. 0 5 10 15 20 39 42 45 48 51 zn = 50 μm T n O C t (s) Ts OC 46 48 50 Threshold of nociceptor 43 O C 0 5 10 15 20 -70 -60 -50 -40 -30 -20 -10 0 46 48 50 zn = 50 μm V t ( m v) t (s) Threshold of pain -55 mV (a) Temperature history (b) Pain level Ts OC Fig. 6. (a) Nociceptor temperature and (b) pain level for selected skin temperatures and fixed nociceptor location. 42 45 48 51 54 -40 -30 -20 -10 0 10 20 Model results Linear fit V t ( m v) Ts OC z n = 50 μm Fig. 7. Influence of stimulus intensity on pain level. 44 F. Xu et al. / Applied Mathematics and Computation 205 (2008) 37–46 F. Xu et al. / Applied Mathematics and Computation 205 (2008) 37–46 45 4.3. Influence of stimulus intensity Naturally, it is expected that the magnitude of skin surface temperature has a great effect on pain. Fig. 6 presents the noci- ceptor temperature and pain level as functions of time for selected values of skin surface temperature. The corresponding relationship between stimulus intensity and pain sensation is plotted in Fig. 7. A nearly linear relation stimulus intensity and pain level exists in the temperature range studied (45–50 �C; Fig. 7). Similar results have been observed in literature. For example, LaMotte and Campbell [32] found that heat pain increases monotonically with stimulus intensity between 40 and 50 �C whlist Torebjork et al. [33] observed a linear relationship between the mean responses of CMHs recorded in conscious humans and the median ratings of pain over the temperature range of 39–51 �C. 5. Summary Preliminary results from theoretical modeling of skin thermal pain are presented. By considering the underlying biophys- ical and neurological mechanisms of pain sensation, a holistic mathematical model for quantifying skin thermal pain is developed. The model is composed of three interconnected parts: peripheral modulation of noxious stimuli, which converts the energy from a noxious thermal stimulus into electrical energy via nerve impulses; transmission, which transports these neural signals from the site of transduction in the skin to the spinal cord and brain; and modulation and perception in the spinal cord and brain. In the model of modulation, a thermomechanical model of skin tissue is employed for the calculation of stimulus intensity at the location of nociceptor. With this model, the intensity of thermal pain can be predicted directly in terms of the character of noxious stimuli: as a corollary, for a given thermal treatment profile, a given medical procedure can be determined to be painful or not. There are nonetheless several limitations of the model, as described below. (1) In the model, the neurons are assumed to be one point, whose numbers and real shape are not taken into account. We will therefore need to consider the actual mor- phology of the neuron. (2) The skin is assumed to be a thermo-elastic material although viscoelasticity has been shown to play an important role in the mechanical behavior of skin [34–37] and touch sensation [38–41]. (3) The parameters used in our model are not initially obtained from experiments on nociceptors, which is due to the fact that there is no appropriate experimental or theoretical data about skin nociceptors. These deficiencies of the model will be addressed in our future the- oretical and experimental studies. Acknowledgements This work is supported by the Overseas Research Studentship (ORS) and Overseas Trust Scholarship of Cambridge Univer- sity, the Xi’an Jiaotong University Research Funding, the National Natural Science Foundation of China (10572111, 10632060), the National 111 Project of China (B06024), and the National Basic Research Program of China (2006CB601202). References [1] N.F. Britton, S.M. Skevington, On the mathematical modelling of pain, Neurochemical Research 21 (9) (1996) 1133–1140. [2] P.D. Picton, J.A. Campbell, S.J. Turner, Modelling chronic pain: an initial survey, in: Eighth International Conference on Neural Information Processing, Shanghai, China, 2001, pp. 1267–1270. [3] U. Fors, M.L. Ahlquist, R. Skagerwall, L.G. Edwall, G.A. Haegerstam, Relation between intradental nerve activity and estimated pain in man – a mathematical model, Pain 18 (4) (1984) 397–408. [4] U.G. Fors, M.L. Ahlquist, L.G. Edwall, G.A. Haegerstam, Evaluation of a mathematical model analysing the relation between intradental nerve impulse activity and perceived pain in man, International Journal of Bio-medical Computing 19 (3–4) (1986) 261–277. [5] U.G. Fors, L.G. Edwall, G.A. Haegerstam, The ability of a mathematical model to evaluate the effects of two pain modulating procedures on pulpal pain in man, Pain 33 (2) (1988) 253–264. [6] U.G. Fors, H.H. Sandberg, L.G. Edwall, G.A. Haegerstam, A comparison between different models of the relation between recorded intradental nerve impulse activity and reported pain in man, International Journal of Bio-medical Computing 24 (1) (1989) 17–28. [7] N.F. Britton, M.A. Chaplain, S.M. Skevington, The role of N-methyl-D-aspartate (NMDA) receptors in wind-up: a mathematical model, IMA Journal of Mathematics Applied in Medicine and Biology 13 (3) (1996) 193–205. [8] N.F. Britton, S.M. Skevington, A mathematical model of the gate control theory of pain, Journal of Theoretical Biology 137 (1) (1989) 91–105. [9] N.F. Britton, S.M. Skevington, M.A.J. Chaplain, Mathematical modeling of acute pain, Journal of Theoretical Biology 3 (4) (1995) 1119–1124. [10] H. Minamitani, N. Hagita, A neural network model of pain mechanisms computer simulation of the central neural activities essential for the pain and touch sensations, IEEE Transactions on Systems, Man and Cybernetics 11 (1981) 481–493. [11] M. Haeri, D. Asemani, S. Gharibzadeh, Modeling of pain using artificial neural networks, Journal of Theoretical Biology 7 (220(3)) (2003) 277–284. [12] M.J. Caterina, D. Julius, Sense and specificity: a molecular identity for nociceptors, Current Opinion in Neurobiology 9 (1999) 525–530. [13] D. Julius, A.I. Basbaum, Molecular mechanisms of nociception, Nature 413 (6852) (2001) 203–210. [14] M.J. Millan, The induction of pain: an integrative review, Progress in Neurobiology 57 (1999) 1–164. [15] J. Brooks, I. Tracey, From nociception to pain perception: imaging the spinal and supraspinal pathways, Journal of Anatomy 1 (2005) 19–33. [16] E.W. McCleskey, M.S. Gold, Ion channels of nociception, Annual Reviews in Physiology 61 (1999) 835–856. [17] R.A. Meyer, J.N. Campbell, S.N. Raja, Peripheral neural mechanisms of nociception, see Wall & Melzack, 1994, pp. 13–44. [18] L. Kruger, E.R. Perl, M.J. Sedivec, Fine structure of myelinated mechanical nociceptor endings in cat hairy skin, The Journal of Comparative Neurology 198 (1) (1981) 137–154. [19] D.B. Tillman, R.D. Treede, R.A. Meyer, J.N. Campbell, Response of C fibre nociceptors in the anaesthetized monkey to heat stimuli: correlation with pain threshold in humans, Journal of Physiology 485 (Pt. 3) (1995) 767–774. [20] D.B. Tillman, R.D. Treede, R.A. Meyer, J.N. Campbell, Response of C fibre nociceptors in the anaesthetized monkey to heat stimuli: estimates of receptor depth and threshold, The Journal of Physiology 485 (Pt. 3) (1995) 753–765. 46 F. Xu et al. / Applied Mathematics and Computation 205 (2008) 37–46 [21] J. Van Hees, J. Gybels, C nociceptor activity in human nerve during painful and non painful skin stimulation, Journal of Neurology, Neurosurgery, and Psychiatry 44 (7) (1981) 600–607. [22] A.V.S. Reuck, J. Knight, Touch, Heat and Pain, Churchill Ltd, 1966. [23] F. Xu, T. Wen, K.A. Seffen, T.J. Lu, Biothermomechanics of skin tissue, Journal of the Mechanics and Physics of Solids 56 (5) (2008) 1852–1884. [24] A. Patapoutian, A.M. Peier, G.M. Story, V. Viswanath, Thermo TRP channels and beyond: mechanisms of temperature sensation, Natural Review in Neuroscience 4 (8) (2003) 529–539. [25] D.M. Cain, S.G. Khasabov, D.A. Simone, Response properties of mechanoreceptors and nociceptors in mouse glabrous skin: an in vivo study, Journal of Neurophysiology 85 (4) (2001) 1561–1574. [26] N.C. James, A.M. Richard, Neurobiology of Nociceptors, Oxford University Press, Oxford, 1996. [27] A.L. Hodgkin, A.F. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve, The Journal of Physiology 117 (1952) 500–544. [28] L. de Medinaceli, J. Hurpeau, M. Merle, H. Begorre, Cold and post-traumatic pain: modeling of the peripheral nerve message, Biosystems 43 (3) (1997) 145–167. [29] A.S. Paintal, Effects of temperature on conduction in single vagal and saphenous myelinated nerve fibres of the cat, The Journal of Physiology (London) 180 (1965) 20–49. [30] A.S. Paintal, The influence of diameter of medullated fibres of cats on the rising and falling phases of the spike and its recovery, Journal of Physiology (London) 184 (1966) 791–811. [31] R. Melzack, P.D. Wall, Pain mechanisms: a new theory, Science (New York, NY) 150 (699) (1965) 971–979. [32] R.H. LaMotte, J.N. Campbell, Comparison of responses of warm and nociceptive C-fiber afferents in monkey with human judgments of thermal pain, Journal of Neurophysiology 41 (2) (1978) 509–528. [33] H.E. Torebjork, R.H. LaMotte, C.J. Robinson, Peripheral neural correlates of magnitude of cutaneous pain and hyperalgesia: simultaneous recordings in humans of sensory judgments of pain and evoked responses in nociceptors with C-fibers, Journal of Neurophysiology 51 (2) (1984) 325–339. [34] S.J. Kirkpatrick, I. Chang, D.D. Duncan, Viscoelastic anisotropy in porcine skin: acousto-optical and mechanical measurements, International Society for Optical Engineering, Bellingham, WA 98227-0010, United States, Saratov, Russian Federation, 2005, pp. 174–183. [35] F. Khatyr, C. Imberdis, P. Vescovo, D. Varchon, J.M. Lagarde, Model of the viscoelastic behaviour of skin in vivo and study of anisotropy, Skin Research Technology 10 (2) (2004) 96–103. [36] J.Z. Wu, R.G. Dong, W.P. Smutz, A.W. Schopper, Non-linear and viscoelastic characteristics of skin under compression: experiment and analysis, Bio- medical Materials and Engineering 13 (4) (2003) 373–385. [37] F.H. Silver, J.W. Freeman, D. DeVore, Viscoelastic properties of human skin and processed dermis, Skin Research and Technology 7 (1) (2001) 18–23. [38] G. Moy, U. Singh, E. Tan, R. Fearing, Human psychophysics for teletaction system design, Haptics-e 1 (3) (2000). [39] J.Z. Wu, K. Krajnak, D.E. Welcome, R.G. Dong, Analysis of the dynamic strains in a fingertip exposed to vibrations: correlation to the mechanical stimuli on mechanoreceptors, Journal of Biomechanics 39 (13) (2006) 2445–2456. [40] R.J. Gulati, M.A. Srinivasan, Human Fingerpad under Indentation I: Static and Dynamic force response, ASME, New York, NY, USA, Beever Creek, CO, USA, 1995. pp. 261–262. [41] D.T.V. Pawluk, A viscoelastic model of the human fingerpad and a holistic model of human touch, Harvard University, 1997. [42] F.A. Duck, Physical properties of tissue: A comprehensive reference book, Academic Press, 1990. [43] D.A. Torvi, J.D. Dale, A finite element model of skin subjected to a flash fire, Journal of Biomechanical Engineering 116 (3) (1994) 250–255. [44] A. Delalleau, G. Josse, J.M. Lagarde, H. Zahouani, J.M. Bergheau, Characterization of the mechanical properties of skin by inverse analysis combined with the indentation test, Journal of Biomechanics 39 (9) (2006) 1603–1610. [45] F.M. Hendriks, D. Brokken, C.W. Oomens, D.L. Bader, F.P. Baaijens, The relative contributions of different skin layers to the mechanical behavior of human skin in vivo using suction experiments, Medical Engineering & Physics 28 (3) (2006) 259–266. [46] W. Elkins, J.G. Thomson, Instrumented Thermal Manikin, Acurex Corporation, Aerotherm Division Report AD-781, 1973, pp. 176. [47] F.C. Henriques, A.R. Moritz, Studies of thermal injury, 1. The conduction of heat to and through skin and the temperatures attained therein. A theoretical and an experimental investigation, American Journal of Pathology 23 (1947) 531–549. [48] W. Roetzel, Y. Xuan, Transient response of the human limb to an external stimulus, International Journal of Heat and Mass Transfer 41 (1) (1998) 229– 239. [49] P. Sejrsen, Measurement of cutaneous blood flow by freely diffusible radioactive isotopes, Danish Medical Bulletin 18 (1972) 1–38. [50] S. Dahan, J.M. Lagarde, V. Turlier, L. Courrech, S. Mordon, Treatment of neck lines and forehead rhytids with a nonablative 1540-nm Er:glass laser: a controlled clinical study combined with the measurement of the thickness and the mechanical properties of the skin, Dermatologic Surgery 30 (6) (2004) 872–879. Modeling of skin thermal pain: A preliminary study Introduction Skin thermal pain Nociceptors in skin tissue Ion channels of nociceptors Model development Model of transduction Thermomechanical model of skin tissue Temperature profile Thermal damage Thermal stress Current Frequency modulation Model of transmission Nerve length ({L}_{{\rm n}}) Conduction velocity ({v}_{{\rm c}}) Transmission time ({t}_{t}) Model of modulation and perception Results and discussion: case study Validation of the model Influence of nociceptor location Influence of stimulus intensity Summary Acknowledgements References