A Simulation Study on Selective Stimulation of C-Fiber Nerves for Chronic Pain Relief

It has been reported that stimulating nociceptive unmyelinated C nerves (<inline-formula> <tex-math notation="LaTeX">$C$ </tex-math></inline-formula>) near the sarcolemma could induce the secretion of endogenous opioids that relieve chronic pain. However, a substantial concern remains: concomitant stimulation might cause acute pain to nociceptive myelinated <inline-formula> <tex-math notation="LaTeX">$\text{A}\delta $ </tex-math></inline-formula> nerves (<inline-formula> <tex-math notation="LaTeX">$A\delta$ </tex-math></inline-formula>), which generally have a lower activation threshold than the <inline-formula> <tex-math notation="LaTeX">$C$ </tex-math></inline-formula>. However, few studies have reported on <inline-formula> <tex-math notation="LaTeX">$C$ </tex-math></inline-formula> selectivity over <inline-formula> <tex-math notation="LaTeX">$A\delta $ </tex-math></inline-formula> (<inline-formula> <tex-math notation="LaTeX">$C$ </tex-math></inline-formula>-selectivity). In this study, the <inline-formula> <tex-math notation="LaTeX">$C$ </tex-math></inline-formula> and <inline-formula> <tex-math notation="LaTeX">$A\delta $ </tex-math></inline-formula> nerves were modeled using the Hodgkin-Huxley (HH) and McIntyre-Richardson-Grill (MRG) models, respectively. Two potential stimulation schemes, including bipolar square waves and burst modulated alternating current, together with a new stimulation scheme named inhibit-<inline-formula> <tex-math notation="LaTeX">$A\delta $ </tex-math></inline-formula> (i-<inline-formula> <tex-math notation="LaTeX">$A\delta$ </tex-math></inline-formula>) that inhibits the excitability of the <inline-formula> <tex-math notation="LaTeX">$A\delta $ </tex-math></inline-formula>, were systematically investigated. Their stimuli were given to the <inline-formula> <tex-math notation="LaTeX">$C$ </tex-math></inline-formula> and <inline-formula> <tex-math notation="LaTeX">$A\delta $ </tex-math></inline-formula> nerves through point electrodes located near the nerve fibers. Simulation results revealed that i-<inline-formula> <tex-math notation="LaTeX">$A\delta $ </tex-math></inline-formula> presented the highest <inline-formula> <tex-math notation="LaTeX">$C$ </tex-math></inline-formula>-selectivity, which provides a basis for non-invasive effective chronic-pain relief.


I. INTRODUCTION
Chronic pain is one of the significant factors that reduce patients' quality of life (QOL) because the condition is poorly understood and minimally treatable using existing therapies [1]. The most notable feature of chronic pain is its uncontrolled continuous pain sensation even without external stimulation. Moreover, a survey showed that the retention rate of chronic pain increases with age [2]. For countries facing rapidly aging populations, the treatment of chronic pain is an essential medical and even social issue [3].
Psychological therapy can improve pain treatment by facilitating positive subjective attitudes, such as treatment willingness and good interpersonal relationships [4]. However, its total effect still needs to be verified [3]. Two main methods are available for chronic pain relief: 1) stimulation therapy (such as massage, physical therapy, and acupuncture) and 2) drug therapy [3]. For mild pain, stimulation therapy can The associate editor coordinating the review of this manuscript and approving it for publication was Yizhang Jiang . help evacuate substances that cause pain in the affected areas. For more severe pain, stimulation therapy can improve the inhibition of projection neurons that transmit noxious information by stimulating large-diameter non-nociceptive afferent nerves, such as Aβ, thereby blocking the signaling of nociceptive nerves and relieving the pain [5]. However, this treatment based on gate control theory cannot achieve a longlasting analgesic effect. For chronic pain, generally, it is more effective to take medication and use non-drug therapy to assist the medication. Opioid analgesics are used for very serious chronic pain, and these drugs are highly addictive and likely to produce a broad spectrum of adverse effects [6].
It was reported in [7] that a pain immunity mechanism could inhibit the transmission of strong pain signals into the brain, causing spontaneous secretion of opioid substances that bind to opioid receptors. Endogenous opioids have the same analgesic effect as synthetic opioid drugs but lack addictive properties [7]. Later, researchers found that electroacupuncture stimulation could release endogenous opioids [8] and that needle electrodes could be used to stimulate peripheral nociceptive afferent nerves [9], prompting the nervous system to secrete opioids that act on the µ-receptor for analgesia [10]. Since then, treatments using endogenous opioid theory for sustained analgesia, such as using laser [11], heat [12], and electrical stimulation [13] to induce endogenous opioids, have been developed. The following two facts are well documented.
1) Different types of pain induce different endogenous opioid substances to bind to their corresponding receptors [14]. Moreover, in some needle stimulation studies, researchers found that different types of opioids were induced by different frequencies of stimulation. Especially, the µ-receptor with the strongest analgesic effect corresponds to low-frequency stimulation; thus, stimulation to C fibers (C) produces the most significant analgesic effects [14].
2) Compared with peripheral nerve endings in the skin surface, nociceptive afferent nerves in deep tissues (such as tissues near the sarcolemma) play a decisive role in pain suppression [15].
Based on these two facts, the research on using needle electrodes for the stimulation of deep tissue nerve fibers for endogenous opioid has progressed further ahead of other means of stimuli, e.g. using laser or heat. Using needle electrodes to stimulate deep tissue nerve fibers, such as the C nerve near the sarcolemma, has been explored to achieve long-lasting analgesic effects (upon induction with µ-related opioids) [16]. However, trauma caused by persistent acute pain elicited by concomitant stimulation of Aδ fibers (Aδ) represents a concern for patients receiving treatment [16]. Moreover, the needle stimulation for the treatment requires specific technical skills that can only be acqiuired with proper training. Thus, it is difficult to push this pain relief technology into daily living use. New techniques in this direction need to be developed to effectively and safely stimulate deep tissue nerve fibers for endogenous opioid to have an accessible chronic pain relief. The research and progress of using needle electrodes will show insights and guidance to the development of the techniques using other kinds of stimuli.
We aim to realize endogenous-opioid based effective chronic pain relief through surface stimulations. This requires stimulation schemes with the following features: 1) high C-selectivity over Aδ (C-selectivity), consideration of multiple adversarial factors in future implementations of surface stimulations, and the requirement for C-selectivity to be as high as possible; and 2) conduction of the stimuli from the skin surface to the vicinity of nociceptive fibers. For the first requirement, i.e., C-selectivity, it is challenging because C fibers are normally 0.2-1.5 µm in diameter and are thinner than Aδ fibers that have a diameter of 1.0-5.0 µm. On the other hand, generally, thicker myelinated nerves are much easily excited than the thin unmyelinated nerves [17] since the bigger membrane area reduces the equivalent resistance. Consequently, to activate thin unmyelinated C fibers before activating thicker myelinated Aδ fibers, the C-selectivity represent a significant difficulty. In this study, we focused on the first requirement/issue to explore the possibility of selective stimulation of C. Moreover, the findings that are useful for the second requirement will also be addressed.
Several papers have reported selective stimulation for motor function restoration [18]. Selective stimulation of afferent nerves has been studied for a long time [19]- [21]. Injecting drugs to realize C-selectivity was addressed in [22]. Other than drugs, electrical stimulation has been used to stimulate different afferent nerves selectively. For electrical stimulations, the effects of different waveforms and stimulation parameters (frequency bands, inter-simulation interval (ISI)) alone and in combination on selectivity and alternation of the activation threshold of certain types of afferent nerves have been investigated [21], [23]. Neurometer (Neurotron, Baltimore, MD, USA) uses sine waves with frequencies of 2000, 250, and 5 Hz to stimulate Aβ, Aδ, and C fibers, respectively [24]. However, as has been noted, low-frequency sine waves cannot guarantee the selective stimulation of unmyelinated C fibers [21,25]. In the low-frequency band, square waves exhibit better stimulating effects on unmyelinated nerves compared with sine waves with long pulse duration (>8 ms) [26]. In [26], low-frequency bipolar square waves (BSW) are reported as an option for C-selectivity stimulation. It is essential to explore the possibility of BSW with asymmetric cathodal and anodal stimuli, which has not been a focus of study in the research area of selective stimulation.
Moreover, high-frequency bipolar stimuli make it difficult to stimulate C fibers. This technique is not suitable for inducing endogenous opioids [24]. However, the influence of low-frequency wave modulated high-frequency waves, such as ''Russian electrical stimulation'' [27], which is more recently known as burst-modulated alternating current (BMAC), on C-selectivity is unclear. This influence is attributed to the fact that high-frequency carrier waves may improve the penetration of stimulus to human skin [28]. Thus, carrier waves might solve another problem for real pain-relief applications via surface stimulation to C-fiber nerves near the sarcolemma. Furthermore, for BSW and BMAC, the effect of the polarity asymmetry ratio of cathodal and anodal stimuli, an important parameter of any charge-balanced waves, on C-selectivity has not been investigated. In summary, a systematic investigation of the effect of BSW and BMAC on C-selectivity is strongly required.
On the other hand, effective endogenous opioid secretion could be realized by not only decreasing the activation threshold of C, but also modifying the activation threshold of, or even blocking the conduction of Aδ. Pre-pulses could be used to modify the activation threshold of afferent nerves. Grill et al. proposed pre-pulse stimulation to change the neural threshold of an unmyelinated nerve model (implemented using the Hodgkin-Huxley (HH) model), which showed that cathodal stimulation suppresses and anodal stimulation promotes nerve excitability [29].
Bostock et al. proposed a method to diagnose neurology diseases according to the neuronal excitability change in response to pre-pulse stimulation to myelinated nerves [30]. More importantly, their results are contrary to the stimulating VOLUME 8, 2020 effects of the HH model; thus, myelinated and unmyelinated nerves might exhibit relatively different characteristics in pre-pulse stimulation.
Animal experiments have demonstrated that nonpharmaceutical nerve conduction block can be achieved in mammals [31]. However, the experiments, which used plier electrodes to stimulate nerves directly by high-frequency high-strength waves, are not suitable for human therapy. Furthermore, a simulation study [32] showed that high-frequency stimulation could gradually open potassium channels of myelinated nerves and increase its threshold, thus temporarily reducing the excitability of the myelinated nerves. This phenomenon can occur even without causing a nerve conduction block [32]. Based on this observation, a new scheme for improving C-selectivity called inhibit-Aδ (i-Aδ) was proposed and compared with the other potential schemes.
In this study, the C and Aδ nerves were modeled using the HH [33], [34] and McIntyre-Richardson-Grill (MRG) model [35], respectively. The stimuli of the three stimulation schemes, BSW, BMAC, and i-Aδ, were given to C and Aδ fibers through point electrodes located near the fibers, and their effect on C-selectivity was systematically compared. Parameters were investigated in a series of experiments designed based on our understanding of ion channel dynamics of the nerve models. To allow a quantitative comparison between different schemes, C-selectivity, which is the selectivity of unmyelinated C fibers over myelinated Aδ fibers, was defined as the ratio (R th ) of the excitation threshold of the two fibers and used as the major evaluation criterion.
Moreover, although the threshold ratio R th is a useful performance index to evaluate the C-selectivity of a stimulation waveform, it is difficult to understand the mechanisms involved in the dynamic changes of each ion channel. Therefore, the phase portrait analysis [36] was employed to disclose the dynamic process involved in the high C-selectivity.
The remainder of the paper is organized as follows. Section 2 introduces models that simulate C and Aδ afferent fibers and also describes the evaluation criterion, constraints of simulation, and the stimulation wave details used in experiments. Section 3 introduces the experimental results of C fiber selective stimulation. Section 4 presents a discussion of experiments and models, the contributions of this study, and future directions. Section 5 presents the conclusions.

A. COMPUTATIONAL MODELS OF C AND Aδ
In this section, the applicability of two nerve models, namely the HH and MRG models, as the computational models for C and Aδ, respectively, and model implementation details are described.

1) THE HH MODEL FOR C
The HH model can be used to describe the behaviors of sodium and potassium currents of the giant squid [33], [37]. This model has been used to model the behaviors of nerves, muscles, and endocrine cells in various organisms ranging from invertebrates to mammals [32], [38], [39]. The model has been applied to mammalian unmyelinated nerves [40]. A model of C that takes action potential propagation into consideration has been established to reproduce experimental data from patients and animals in earlier experiments [41]. In our study, the basic portion of the model was the same as that reported in [41] and expressed by differential equations (1) to (3). These equations were adapted to human C fibers by reducing the resistivity, fiber diameter, etc. according to the experimental results in [42], [43]. Equation (1) shows the change in transmembrane current by external stimulation. Equation (2) presents an expression of the ionic transmembrane current I ion , and equation (3) shows the differential equation of the variables.
In equation (1), C m is the membrane capacitance, V m is the intracellular potential, t is time, I ion is the ionic transmembrane current, and I ext is the externally applied stimulation current. In equation (2), m, h, and n define the open status of those ion channel gates as the opening or closing of sodium ion channels and the opening of potassium ion channels, respectively. In addition, g Na , g K , and g Lk represent the maximal conductance of sodium and potassium ion channels and the current leakage, respectively. E represents the resting potential of each channel. In equation (

2) THE MRG MODEL FOR Aδ
Research efforts have been made to develop models of myelinated nerve fibers [44], [45]. Among them, the MRG (McIntyre, Richardson, and Grill) model, which was originally a geometrically and electrically accurate model of mammalian motor nerve fibers [35], has been used to simulate myelinated afferent nerves [46]- [48], and its validity has been demonstrated by a couple of mammalian motor nerve experiments [49], [50].
The MRG model is a double-cable axon model consisting of nodes separated by inter-medullary segments wrapped by myelin ( Fig. 1. B). The internode is divided into ten segments: two node myelin sheath attachment segments (MYSA), two fluted region segments (FLUT), and six stereotyped internode segments (STIN) ( Fig. 1: adapted from [35]). Mathematically, the MRG model can be expressed as equation (1) based on the study of the classical HH model. However, for the different ion channels, the ionic transmembrane current I ion The nodal segment of the Ranvier segment (i.e., nodal circuit on the right) contains fast sodium (Na), persistent sodium (Na p ), and slow potassium (K ) channels, and leakage resistance (Lk) with nodal capacitance (C n ). The inter-nodal segment contains resistance and capacitance of the myelin sheath (G m and C m ) and inter-nodal double-cable structures (G i and C i ) [35].
in equation (1) should be expressed as follows: where g Naf is the nodal membrane conductance of fast sodium channel, g Nap is the conductance of the persistent sodium channel, g Ks is the conductance of the slow potassium channel and g Lk refers to the linear leakage conductance. m, h, p, and s define the open probability of each ion channel in the MRG model. Among them, m and h represent the opening and closing of fast sodium ion channels, respectively. In addition, p defines the opening of sodium persistent ion channels, and s defines the opening of slow potassium ion channels.

B. EVALUATION CRITERIA
The activation threshold of C fibers (Th HH ) calculated from the HH model, the activation threshold of Aδ fibers (Th MRG ) calculated from the MRG model, and their ratio R th were used to assess C-selectivity. The threshold ratio, R th , is expressed as follows: Based on equation (5), a lower Th HH and a higher Th MRG lead to a lower R th , which indicates increased C-selectivity.

C. PHASE PORTRAIT ANALYSIS
To understand and interpret the mechanism behind high C-selectivity, phase portrait analysis of membrane potential and other important ion channel variables was employed. Phase portrait analysis has been used to analyze the behavior of the nonlinear equations of the HH model [51]. Fig. 2 shows the change of the time constant of each variable as the membrane potential V m increases in the HH and MRG models. It is clear that the time constant of variable m is considerably reduced compared with h and n, which means the m exerts the greatest effect on system behavior. Based on this observation, it is reasonable to set h and n constant at their resting value for the next step in phase portrait analysis. Three singular points are present in the phase portrait of (V m , m): a stable resting point, a threshold saddle point, and an exciting stable point. In general, the closer the threshold saddle point to the stable resting point, the more susceptible the nerve is to excitation. The change in threshold saddle point also indicates changes in model excitability.
Although h and n were assumed to be constants, they are also variables, and changes in h and n cause changes in membrane potential as well as the threshold of the saddle point and exciting stable point in the phase portrait. For the convenience of discussion and given their relatively small  Fig. 3 C: inhibit-Aδ wave (i-Aδ) contains a high-frequency symmetric sine wave as its pre-pulse for Aδ inhibition. With frequency, intensity, and duration as the parameters, and BSW or BMAC with optimal parameters identified through a previous investigation as its main pulse. Fig. 3 D: An illustration of the change in membrane potential and potassium ion channel permeability variable s by high-frequency sine wave stimulation. Marks a, b, and c were used to express the three moments after the high-frequency pre-pulse stimulation, a: 0 ms, b 10 ms, and c 30 ms. effect on membrane potential, only constant h and n cases were analyzed and discussed.
The phase portrait analysis has been applied to myelinated nerve models, such as the Frankenhaeuser-Huxley (FH) model [52]. However, few studies have assessed the dynamics of the MRG model in the phase portrait, which may lead to stimulation schemes favoring C-selectivity over Aδ if compared with the dynamic behavior of C. Therefore, in this study, phase portrait analysis of the HH and the MRG model responding to stimuli of different stimulation schemes with different parameters was performed to clearly define the mechanism underlying high C-selectivity.

D. POTENTIAL STIMULATION SCHEMES TO REALIZE SELECTIVE STIMULATION FOR LONG-TERM PAIN RELIEF
In this section, the waveforms, parameters, and features of the three potential stimulation schemes, including BSW, BMAC, and i-Aδ, are schematically explained.  Introduction section. Each period contains an anodal stimulus and a cathodal stimulus. The polarity asymmetry ratio is the quotient of the amplitude of the cathodal stimulus and the anodal stimulus, which has not been investigated in depth in the literature. Square waves with different stimulus frequencies (from 5 to 500 Hz) and polarity asymmetry ratio values (from 1:6 to 6:1) were investigated in terms of C-selectivity.
2) BMAC Fig. 3. B shows the stimulus waveform of the BMAC. Besides burst frequency and carrier frequency, the ratio of interstimulus interval (ISI) in carrier waves are decisive parameters. The burst frequency investigated ranges from 5 to 500 Hz, while the carrier waves investigated range from 100 to 10 kHz. In addition, their interaction with burst waves and subsequent effects on C-selectivity were investigated. Moreover, the C-selectivity of the BMAC with different polarity asymmetry ratio values was compared with that of the BSW.

3) I-Aδ
The inhibit-Aδ wave (i-Aδ) ( Fig. 3. C) contains a high-frequency sine wave as its pre-pulse for Aδ inhibition; frequency, intensity, and duration as parameters; and BSW or BMAC with optimal parameters identified through the previous investigation as its main pulse. The stimulation intensity, pulse duration and frequency of the pre-pulse were investigated in the study. Table 1 provides the assessed parameters. For the main pulse, BSW and BMAC waves with a burst frequency of 50 Hz and duration of 20 ms were tested. Fig. 3. D presents an illustration of the change in membrane potential and potassium ion channel permeability variable s based on high-frequency sine wave stimulation. Three important moments a, b, and c are defined to explain the role of variable s in the Discussion section.

E. SIMULATION SETUP
For both nerve fibers, only a single fiber case was considered. Each fiber consisted of 40 nodes. The cathode electrode was set 1 mm away from the middle node along the radial direction, and the anode electrode was set far enough away from the object fiber along the radial direction. Stimuli with a specific current density (i.e., thresholds of each model) were applied to each nerve fiber model through the cathode electrode. The neural behaviors of each modeled nerve fiber were recorded and analyzed.
For the analysis, the unmyelinated and myelinated nerve models were first confirmed using past animal experimental data and then used to evaluate and compare three different stimulation schemes targeting high C-selectivity, which were evaluated using the C-selectivity index R th .
Both models were implemented on MATLAB 2013a (The MathWorks, Inc., Natick, Massachusetts, USA). The 4th-order Runge-Kutta method was employed to solve differential equations of the HH and MRG models numerically. A time step of 0.1 µs was used in the numerical computation. All the computer simulations were run on a desktop computer with a quad-core, 3.4 GHz processor. It takes 240 s for the Matlab code to compute one axon.

A. SIMULATION VALIDATION
The implemented HH and MRG models were first validated from two different aspects: strength-duration (S-D) curves and conduction velocity. Fig. 4. A shows the S-D curve of the Aδ model (red dotted line) with the curves measured from the human median and ulnar nerves [53] for sensory axons (gray zone). Fig. 4. B shows the S-D curve of the C model (red dotted line) with the curves measured from cat nerves (gray zone) [54]. Both models matched experimental results in the literature.
The conduction velocity of the C model (diameter: 1 µm) and Aδ model (diameter: 5 µm) were 0.98 m/s and 27.5 m/s, respectively, which is consistent with measurement experimental results of 0.5-2 m/s and 5-40 m/s for C and Aδ, respectively, using the same diameters as simulation settings as reported in [55], [56].  5. A shows the acquired thresholds of the C and Aδ models, namely, the HH and MRG, respectively, stimulated by sine and square waves at a different frequency. As shown in the figure, using square stimulation, C has a lower activation threshold than Aδ when the stimulation frequency is less than 100 Hz. However, using sine stimulation, C showed a lower activation threshold than Aδ only for the frequency range of 50-100 Hz. Moreover, the thresholds of the sine stimulus are always greater than those for square stimulation at the same frequency. Fig. 5. B shows the threshold of the C and Aδ models when the polarity asymmetry ratio was changed but the charge balance was maintained for stimulation safety, thereby resulting in an equal integral value of the anodal and cathodal stimuli. The total duration of the stimulus is 20 ms. The trend of the threshold change caused by the sine stimulus and the square stimulus is similar, but the threshold values of the sine stimulus are always greater than that of the square stimulus. From a ratio value 1:6 to 2:1, C showed a lower threshold than Aδ. There is a crossover between 2:1 and 3:1. After the crossover, it is more difficult to activate C than Aδ.

2) BMAC
Considering the advantage of high-frequency waves in deep stimulation, it is necessary to verify the effect of the carrier waves on the selective stimulation. Fig. 6 shows the threshold change when stimulated by pulses containing carrier waves with different frequencies (100 Hz-50 kHz) and burst waves of 5, 20, and 50 Hz. At lower carrier frequencies, the C and Aδ have a similar activation threshold. However, the threshold of the Aδ increases more than that of C when the carrier frequency increases. Fig. 7 shows the activation threshold of the two models when they were stimulated by burst waves with different frequencies and polarity asymmetry ratio values. The carrier frequency was 10 kHz. The results were identical to those of the BSW shown in Fig. 5. Fig. 8 shows the activation threshold of the two models when they were stimulated by carrier waves with different ISI. The burst frequency of BMAC was 50 Hz, and the carrier frequency was 10 kHz. The activation threshold is approximately linear as the stimulus duration increases. However, the charge required for stimulation remains the same. Fig. 9 shows the effect of different parameters of the prepulse: frequency, duration, and intensity on the potassium ion channel permeability variable s. The opening of the potassium ion channel leads to an increase in the activation threshold.

3) I-Aδ
As shown in Fig. 9 A1, B1, and C1, stronger and more long-lasting stimuli with a lower frequency result in a higher variable s of the myelinated nerves. BMAC is used as the stimulus with duration of 20 ms, carrier frequency of 10 kHz, and polarity asymmetry ratio of 1:1. Using the strength-s and strength-R th graphs ( Fig. 9. A1, A2) as examples, the higher strength of high-frequency pre-pulse stimulation could increase the value attained by variable s of the myelinated nerves. When the variable s reached a certain value, the nerve could be activated. Therefore, the strength-s curve shows a plateau phase at high stimulation strength, indicating  activation of the myelinated nerves by the pre-pulse stimulation, which should be avoided considering that the purpose of introducing the pre-pulse is to improve C-selectivity. As the duration of pre-pulses increases ( Fig. 9. B1, B2), the variable s first increases rapidly, then increases slowly, and eventually saturates at a certain value depending on the intensity and frequency of the stimulus. The change of the threshold ratio R th is consistent with the change in variable s. The frequency-s and frequency-R th graphs ( Fig. 9. C1, C2) exhibit an inverse tendency compared with graphs of strength and duration. Fig. 9 A2, B2, and C2 present the threshold ratio R th with BMAC or BSW as the main pulses, and pre-pulse parameters that are the same as those explored in Fig. 9 A1, B1, and C1, respectively. BMAC is used as a stimulus with duration of 20 ms, 10-kHz carrier frequency, and 1:1 polarity asymmetry ratio. Fig. 10 shows the change in R th when stimulated with a pre-pulse wave with different strength values. The other parameters of the pre-pulse waves include duration of 20 ms and frequency of 100 kHz followed by an asymmetric BMAC wave with burst duration of 20 ms and different polarity asymmetry ratio values (1:1, 1:3, and 1:6). When the strength was 0, there was no pre-pulse wave before the BMAC. The results are similar to the results shown in Fig. 5. B and Fig. 7. B. Polarity asymmetry ratios of 1:3 and 1:6 exhibited reduced R th values (0.3931 and 0.3638, respectively) compared with that of 1:1 (0.6075). However, the three curves showed the same tendency.

1) EFFECT OF BIPOLAR STIMULUS
C-selectivity could be interpreted and confirmed by phase portrait analysis. Fig. 11 shows the phase portrait of the two nerve models for the sodium ion channel variable m and the membrane potential V m , in their resting state, during the anodal stimuli, and responding to BMAC with duration  of 20 ms, carrier frequency of 10 kHz and an polarity asymmetry ratio of 1:1. In Fig. 11. C., points d and e express the resting state and the peak of an action potential, respectively. Point e is the saddle point, expressing the threshold  , p: 0.0914, s: 0.007), which means that it has a lower threshold and is more likely to be excited. On the other hand, the HH model exhibits greater variability (>4 µA/cm 2 ). Thus, when given an anodal stimulus with sufficient intensity, its equilibrium points that indicate the threshold of the nerve fiber disappear in the phase portrait. In contrast, the change in the MRG model is minimal, and its equilibrium points show minimal changes, too. When the anodal stimulus is not sufficient to cause unmyelinated threshold reduction, the equilibrium points of the phase portrait did not disappear, as shown Fig. 11. E. These findings demonstrate the behavior of the HH model in response to a 4 µA/cm 2 anodal stimulus. The variability in the HH model indicates that it is possible to change its activation threshold by controlling the stimulation parameters. Fig. 12 shows the phase portrait of the two nerve fiber models after applying a pre-pulse with 100 kHz and 29 mA/cm 2 sine waves for 20 ms, as illustrated in Fig. 3. C. Three states are introduced in Fig. 3. D, which are indicated as a, b, and c in Fig. 12. These points represent the moment when the pre-pulse stimulation ended, 10 ms after sine stimulation ended, and 30 ms after sine stimulation ended and reached its resting state, respectively. The change of the threshold saddle point (from c to a) of the MRG model ( V 1 + V2 = 8.0 mV) is greater than that of the HH model ( V = 1.9 mV) with the same high-frequency sine stimulus. Regarding the recovery of the saddle point (expressing threshold) to its resting state, the MRG model ( t 2 is approximately 35 ms) is slower than that of the HH model ( t 1 is approximately 8 ms). Thus, the high-frequency sine wave stimulation has larger and longer effect on the change of the Aδ threshold than that of C. Phase portrait analysis disclosed the possibility that pre-pulses improve C-selectivity. In other words, pre-pulses may promote endogenous opioid secretion.

IV. DISCUSSION
Neural models have been used to study the behavior of unmyelinated and myelinated sensory nerves [33], [35]. In section 3.1, the strength and duration results of the HH model for unmyelinated C fibers and the MRG model for myelinated Aδ fibers were compared with data recorded in human and animal measurement experiments. As shown in Fig. 4, these experiments were well matched, verifying the validity of HH and MRG as models of the two sensory nerve fibers.
In the previous research, the effect of the pulse duration on different individual nerves was studied [53], [54]. The threshold of one single unmyelinated nerve was greater than that of one single myelinated nerve at any pulse duration, and their threshold values were similar when stimulated with longer duration [53], [54]. This difficulty is experienced in the selective stimulation of C. Other decisive factors are needed to reverse the threshold of C to Aδ. Through the systematic investigation of the effect of the three stimulation schemes, we expected to identify such potential important factors for improving C-selectivity to determine the stimulus waveform that effectively facilitates the secretion of endogenous opioids.

A. EFFECTS OF THE THREE STIMULATION SCHEMES 1) BSW
BSW was compared with sine waves and monopolar pulses, and the effect of the polarity asymmetry ratio was assessed.
As shown in Fig. 5. A, when the frequency is increased to greater than 100 Hz, which means that the duration of the bipolar square stimulus is less than 5 ms, the threshold of the HH model is greater than that of the MRG model. This requirement of long-duration low-frequency stimulation is consistent with a previous report [24] that highlighted the effect of a sine wave with very low frequency (5 Hz) on C. However, sine waves have a reduced stimulation effect compared with square waves. C activation relies on the difference between the values of the variables m and h, reflecting the opening and closing of the fast sodium ion channels, respectively (Fig. 2. A). Given its slowly increasing stimulation intensity, the low-frequency sine wave allows the variable h to change for a longer period of time, thus resulting in a smaller difference between the m and h. Thus, it is more difficult to activate C than square waves. Therefore, the square wave stimulus is more suitable for high C-selectivity.
Some researchers denoted that pre-pulses with different strength and duration have various effects on nerve dynamics [29], [30], [33]. Our simulation results in Fig. 5. A showed that the anodal stimulus has minimal inhibitory effects on the MRG model (SQUARE.MRG greater than MONO.MRG) but a significant excitation effect on the HH model (SQUARE.HH lower than MONO.HH). The anodal stimulus can reduce the threshold of the C, which can be further made clear by phase plane analysis in section 4.2. As shown in the figure, a greater effect is achieved by an anodal stimulus with increased strength and longer duration. Therefore, bipolar pulses are more effectively stimulate C fibers than monopolar pulses.
More importantly, through our systematical investigation of the BSW, the role of the polarity asymmetry ratio in C-selectivity was clarified; this is a necessary consequence of the following arguments.
First, considering the safety of the stimulus waveforms, the stimuli bipolarity is a necessity for both maintain the charge (the product of strength and duration) balance between the anodal and cathodal part of their stimulation waveforms [57] and improving nerve excitability [29]. Secondly, in the case of symmetric polarity balance, a strong anodal stimulus required the matching of a strong cathodal stimulus that is about to activate the Aδ and lead to unnecessary sharp pain. Third, as shown in Fig. 5. B, the strength of stimulation has a greater influence on the threshold of C. When the strength of the anodal stimulus is strong and its duration is short (i.e., as noted with a polarity asymmetry ratio of 1:6, as shown in Fig. 5. B), the threshold of the HH model is lower than the stimulation with weak strength but long duration. This simulation result indicated that the anodal stimuli with stronger and shorter duration resulted in a high polarity asymmetry ratio, benefiting C-selectivity.
Since the polarity asymmetry ratio has not been in terms of C-selectivity, it is a new dimension for designing effective selective stimulation. However, it is not preferable to further increase the polarity asymmetry ratio because it results in either a stronger anodal stimulus, which might be limited by the highest stimulation strength for human subjects, or weaker cathodal stimulus, which could not provide sufficient power to activate C. Thus, it could be optimized based on the details of the selective stimulation, which is the aim of our future study.

2) BMAC
When stimulated by carrier waves modulated by burst waves, i.e., BSW, as shown in Fig. 6, the thresholds of both myelinated and unmyelinated nerves increased as the carrier frequency increases; however, the threshold of the unmyelinated nerves increased and saturates more rapidly compared with myelinated nerves. The threshold of C is considerably reduced compared with Aδ. The only exception is carriers with frequencies of 100 Hz to 200 Hz, at which point the values were similar. As reported in [58], low-frequency burst-modulated waveforms (burst frequency: 20 Hz with 2-5 pulses in each burst period) can reduce the activation of Aδ even if the equivalent carrier frequency was quite low [58]. Our results in Fig. 6 showed that the effect on selective stimulation of C increased significantly when the carrier frequency increased from 200 Hz to 1 kHz. After 1 kHz, the stimulation effect increases very slowly, and the effect is saturated at 10 kHz. Therefore, based on our simulation results, the carrier waves with higher frequency (>10 kHz) favor much more increased C-selectivity.
Comparing the results of BMAC in Fig. 7 with those of BSW shown in Fig. 5, the strength value of the activation threshold of both nerves to BMAC (carrier frequency 10 kHz) is greater due to the reduced charge as a consequence of adding carrier waves. However, the shape of the curves in Fig. 7 is similar to that of the BSW-only case (Fig. 5), indicating that a high-frequency carrier (>10 kHz) has no further contribution to myelinated or unmyelinated nerve fibers simulated by the MRG and HH model, respectively.
This finding partially agrees with Grill et al.'s results regarding the fact that a high-frequency carrier (>2 kHz) does not further contribute to the activation of myelinated fibers [59]. However, for the HH model, when simulating unmyelinated sensory nerves, there have not been any reports on their responses to high-frequency carriers to the best of our knowledge. Through this experiment, it is clear their activation threshold is not affected by the high-frequency carrier waves when their frequency is greater than a certain hertz (approximately 1 kHz), which is similar to that of the MRG model.
Although the increase in duration causes a reduction in the thresholds, as shown in Fig. 8, the overall required charge did not change. Concretely, the charge of BSW at 50 Hz to activate the HH model is 47.4 ms × µA/cm 2 . This value is exactly the same as that of BMAC with a burst frequency of 50 Hz and carrier frequency of 10 kHz. This finding is consistent with that reported in a previous study [59], revealing that the frequency of carrier waves does not contribute to the activation of the high-frequency band because it depends on the charge required for stimulating the MRG model. This experiment demonstrates that this notion is also true for the HH model. In summary, carrier waves in the carrier frequency range (2-10 kHz) favor C-selectivity. In addition, increasing the carrier frequency beyond 10 kHz does not further improve C-selectivity.
On the other hand, as shown in Gabriel's research, the electrical properties of human tissues are frequency dependent [28]. For percutaneous electrical stimulation, high-frequency stimulation could be effective to lower the conductivity of the epidermis to reduce the current consumption at the epidermis and enhance stimulation to targeted deep layers. Therefore, although the use of BMAC to improve C-selectivity saturated after a carrier frequency of 10 kHz, the findings could be helpful to select the correct carrier waves to realize the stimulation to C in deep layers while maintain the high C-selectivity that resulted from BSW. Future studies are needed to confirm this point.

3) I-Aδ
The feasibility of high-frequency stimulation for nerve conduction block was first confirmed in mammalian peripheral motor nerve experiments [16], [31]. The underlying mechanism is that the open-close function of the myelinated slow potassium ion channel, denoted by variable s, is slower than that of the other channels at rest potential as shown in Fig. 2. B, and this rate slowly returns to normal after the high-frequency stimulation ends [32]. This mechanism was first investigated for nerve selective stimulation in this study, harnessing the differential effect of high-frequency stimulation on myelinated and unmyelinated nerves to increase the threshold gap between the two nerves and thereby improving the R th .
As shown in Fig. 9, stimuli with stronger intensity, longer duration, and lower frequency increased the value of the variable s. However, the stimulation parameters should be carefully selected to prevent activating the myelinated nerves by the pre-pulses of the i-Aδ. Moreover, as noted in Fig. 10, it is clear that the pre-pulses could work with the main pulses, especially the asymmetric polarity of BSW and BMAC, to contribute to C-selectivity. If no clear interaction occurred between different factors contributing to C-selectivity, the final R th might reflect the product of multiple factors, such as those resulting from the pre-pulse and the polarity asymmetry ratio. This notion could be revealed by the results shown in Fig. 10. When there was no pre-pulse (the strength of the pre-pulse was 0, only-BMAC case), the quotient of the R th of the stimulation with a polarity asymmetry ratio of 1:1 and 1:6 was 1.669 (0.6074:0.3638). When the strength of the pre-pulse was 20 and 25 mA/cm 2 (pre-pulse with BMAC cases), the quotients were 1.688 (0.5412:0.3206) and 1.657 (0.5000:0.3018), respectively. In other words, i-Aδ and BMAC waveform influence C-selectivity independently. However, this should be further investigated by exploring the different methods to combine the pre-pulses and main pulses, e.g., beginning main pulse stimulation at different levels of the variable s, connecting them with the anodal or cathodal stimulus at different intervals, etc.

B. NERVE MODEL DYNAMICS
To explain the phenomenon whereby the activation threshold of unmyelinated nerves is less than that of myelinated nerves, phase portrait analysis was used based on the investigation of the time constant of different ion channel variables (Fig. 2). The phase portrait analysis of the HH model has been described in detail in the study of Fitzhugh et al. [51]. However, there have been few reports on its application to the MRG myelinated fiber model, especially in terms of selective nerve stimulation.
The phase portrait of the HH and MRG models in Fig. 11 demonstrate that given the anodal stimulus, the polarization equilibrium point e of the HH model moves downwards close to the resting static equilibrium point d until the two equilibrium points disappear and only the hyperpolarization equilibrium point f remains in the phase portrait. In this situation, the state of the nerve model should have moved towards the hyperpolarization equilibrium point, causing the nerve to activate. However, given the influence of the external anodal stimulus, the neural state remains in its current state. When the external anodal stimulation is complete, the nerve terminates its equilibrium state and moves to the only equilibrium point left. Moreover, since the closing of the sodium ion channels variable h shown in Fig. 2 is greater than the opening variable m, C requires a longer time to close the sodium ion channel for activation.
As shown in Fig. 12, high-frequency sine pre-pulse waves cause the resting state isocline of dV m /dt = 0 (blue curve in Fig. 11) to increase from state c to state a. This effect increased the threshold (i.e., point e, in Fig. 11) of both the C and Aδ. The change in threshold of the Aδ is considerably increased compared with C ( V 1 + V 2 = 8.0 mV vs. V = 1.9 mV). In addition, the recovery time of the Aδ, which is the period of time for the threshold at state a to return to state c, is longer than that of C (35 ms vs. 8 ms). This difference in the phase portrait shows that the threshold of both C and Aδ increases after receiving the high-frequency sine stimulus, while the threshold of Aδ is more affected and longer. This feature greatly benefits C-selectivity, providing a powerful influencing factor that promotes the secretion of endogenous opioids.

C. CONTRIBUTION AND LIMITATIONS
The ultimate goal of our study is to selectively stimulate unmyelinated C near the sarcolemma with surface electrodes. In this study, potential stimulation waveforms were explored by investigating their threshold compared with myelinated Aδ and the dynamic behavior of the unmyelinated and myelinated nerve models responding to the stimuli from a point source near the two nerves. The major contributions of the study are summarized as follows.
1. For BSW and BMAC, the polarity asymmetry ratio was first identified and explored as a new dimension to design stimulation waveforms for C-selectivity.
2. For BMAC, a frequency band of carrier waves that improve C-selectivity was identified thorough investigations of the different carrier frequency-threshold characteristics of C and Aδ. The result is partially consistent with relevant studies on MRG in the literature. Regarding the effect of carrier frequency on C, our study is the first attempt to investigate these effects compared with Aδ to the best of our knowledge.
3. The feasibility of high-frequency stimulation, which was first confirmed by effective nerve conduction block in a mammalian peripheral motor nerve experiment [31], was first investigated for selective nerve stimulation in this study. It was demonstrated that high-frequency pre-pulse stimulation inhibits the activation of Aδ by increasing the value of the variable s, which plays a key role in the MRG model and can be used to design pre-pulses and its junction with main pulses. The mechanism was demonstrated by applying phase portrait analysis to both the HH and MRG models.
However, the effect of wave conduction from the skin surface to deep layers and multiple factors affecting nerve activation and selectivity in actual stimulation cases, such as the distribution of free endings and nerve fibers, nerve fiber radius, and tissue morphology near the stimulation area, have not been investigated. In our future studies, a conduction model will be introduced to verify the effect of the stimulation waveforms proposed and explored in this study on C-selectivity to ensure the possibility of surface stimulation in promoting the secretion of endogenous opioids. Moreover, the safety issue should be reconsidered. Studies have demonstrated that when a mammalian unmyelinated nerve is stimulated with an intensity several times its threshold, the nerve will lose its response [60], which clearly revealed the upper limit of unmyelinated nerve stimulation. The stimulation parameters should be optimized under this constraint.

V. CONCLUSION
In this study, to investigate the increased selectivity of C fibers over Aδ fibers and subsequent secretion of endogenous opioids, we compared three stimulation waveforms, two canonical waveforms (BSW and BMAC), and one original scheme proposed based on the nerve block mechanism using an unmyelinated nerve fiber model and a myelinated nerve fiber model that was implemented based on HH equations and MRG equations, respectively. We verified the action potential propagation and their behavior in the phase plane. The ratio between the activation threshold of C and Aδ, namely R th , was used as the C-selectivity index. The results of the computational experiments showed that the polarity asymmetry ratio of BSW, the carrier frequency of BMAC, and the pre-pulses represent important factors when designing stimulation schemes for C-selectivity. Moreover, i-Aδ takes advantage of all three factors to achieve high C-selectivity. This study is an important step towards effective noninvasive chronic-pain relief.

APPENDIX I
See Table 2.

APPENDIX II
See Table 3.