Computational Analysis of Kilohertz Frequency Spinal Cord Stimulation for Chronic Pain Management

SPINAL cord stimulation (SCS) is the most common neurostimulation therapy to treat neuropathic pain conditions that are refractory to conventional medical management. Clinical SCS was first tested in 1967, and its demand has dramatically increased over the years with approximately 35,000 systems sold in 2008 alone.^1,2  SCS is a Food and Drug Administration–approved therapy, typically considered as a final treatment option, with a primary indication for refractory neuropathic limb pain. In the United States, SCS is primarily used to manage failed back surgery syndrome and complex regional pain syndrome.^3–5 

Although conventional SCS applied at a rate between 40 and 80 Hz has been a widely used clinical therapy for decades, it has a limited success rate (approximately 50% of patients receive ≥50% reduction in pain).⁵  There has been recent interest in the use of much higher frequencies in an attempt to improve the clinical results with SCS. Kilohertz stimulation frequencies have shown the ability to generate rapid and reversible conduction block in peripheral nerve models^6–8  and have gained significant attention in recent years. Initial clinical data with a novel device capable of delivering kilohertz frequency SCS (KHFSCS) suggest promising clinical benefits and paresthesia-free effects.^9–12  These studies have also shown a patient preference for KHFSCS over conventional SCS and the ability of KHFSCS to provide pain relief in patients who failed conventional SCS. However, results with KHFSCS remain inconsistent. A recent double-blind, placebo-controlled crossover trial concluded that KHFSCS was not better than sham treatment.¹³ 

The limited and contradictory clinical data available for KHFSCS emphasize the need for a detailed and systematic characterization of its therapeutic mechanisms that currently remain unclear. Several potential therapeutic mechanisms have been suggested: direct activation, conduction block, pseudospontaneous activation, transmission delays, or conduction failure at branch points.^8,14,15  Direct conduction block of action potentials is often considered the most logical mechanism because therapeutic KHFSCS does not produce paresthesias.^9–12 

To determine the pain relief mechanisms of KHFSCS, we must understand the electric fields generated by the stimulation waveform and its direct effects on the neural elements of the spinal cord. This knowledge can be difficult to gain experimentally, and, in the past, several groups have used computational models to study conventional SCS.^16–23  The goal of this study was to use similar theoretical techniques to investigate the effects of KHFSCS on the spinal cord. Our model infrastructure consisted of a finite element model (FEM) of an SCS lead implanted in the epidural space along with multicompartment cable models of dorsal root (DR) and dorsal column (DC) fibers in the spinal cord. This approach permitted systematic characterization of numerous variables and their influence on the direct effects of KHFSCS: waveform shape, dorsal cerebrospinal fluid (dCSF) thickness, lead location, fiber collateralization, and fiber size. The data indicate that direct activation of the spinal cord elements may be possible with KHFSCS; however, it is unlikely that clinical KHFSCS generates direct conduction block within the spinal cord.

We used computer models to investigate the direct response of spinal cord axons to KHFSCS used in chronic pain management. The computer models had two main components: (1) FEM of an SCS lead implanted in the epidural space and (2) electrical models of spinal cord axons. We performed the model analysis in three steps: (1) from the FEM, we calculated the extracellular voltages generated in the spinal cord and surrounding tissue during KHFSCS; (2) we generated electrical models of axons within the spinal cord; and (3) we assessed the direct axonal response to KHFSCS by applying the extracellular voltages (step 1) to the axon models (step 2). The text below provides an overview of these modeling procedures (see the appendix for a detailed description of the modeling parameters).

The first step in our model analysis was to estimate the extracellular voltages generated in the spinal cord during KHFSCS. We performed this estimation using finite element analysis. Finite element analysis is a computerized method for predicting how an object (i.e., the spinal cord and surrounding tissue) reacts to various forces (i.e., the electric fields generated during KHFSCS).²⁴  In finite element analysis, the object is represented by thousands-to-millions of geometrical shapes or finite elements, such as tetrahedrons. Mathematical equations containing information (i.e., electrical conductivity) connecting each point in the object are then used to estimate the response (i.e., extracellular voltage) of each finite element. A computer then sums the response of each individual element to estimate the response of the complete object.

To perform this analysis, we created a three-dimensional FEM of the lower thoracic spinal cord and its surrounding anatomy (fig. 1A). The FEM consisted of the gray and white matter of the spinal cord, surrounding cerebrospinal fluid, dura, epidural fat, vertebral bone, and a surrounding general thorax layer. The dimensions of the spinal cord and the white and gray matter boundaries were defined by human cadaver samples of the lower thoracic spinal cord.²⁵  The FEM also contained an explicit representation of an eight-electrode percutaneous lead implanted in the epidural fat dorsal to the spinal cord. The electrode was placed on the dorsal surface of the dura along the spinal cord midline. Electrical conductivities were assigned to each domain based on experimental data available in the literature (Table 1).^21,22,26  In this study, all simulations were performed for bipolar stimulation with a separation of 8 mm center-to-center between active electrodes (fig. 1A).

The extracellular voltages generated during KHFSCS were the output of this first step in the model analysis. To calculate the voltage at each point in the model tissue, we placed current sources at the appropriate stimulating electrodes, set the outer model surface to ground (i.e., 0 V) and then solved the Poisson equation. These tissue voltages were then interpolated onto the spinal cord axon models described below (step 2).

The second step of our model analysis was to define computer models of spinal cord axons. The fundamental purpose of SCS is to modulate neural activity with electric fields. Theoretical and experimental studies have demonstrated that axonal activation is the principal effect of stimulation within the central nervous system.²⁷  With regard to SCS, studies have shown that the two axon types most likely affected by SCS are the large-diameter myelinated DR fibers and Aβ fibers within the DCs.^21,28  Therefore, we included computer models of both DR and DC fibers in our analysis (fig. 1B). DC and DR fibers were represented by a previously published compartmental model of a mammalian axon.²⁹  In this axon model, the nodes of Ranvier contained active (i.e., voltage-gated fast Na⁺, persistent Na⁺, and slow K⁺ ion channel conductances) and passive (i.e., leak conductance, capacitance) membrane properties.

Dorsal column fibers, running longitudinally along the rostrocaudal axis, were placed on a regular grid (200 μm for the mediolateral direction and 100 μm for the dorsoventral direction) within the white matter boundary of the DC defined by the FEM.²² 

Dorsal root fibers had a three-dimensional axon trajectory in which they entered the spinal cord at a 45-degree angle with respect to the transverse plane and approximated the anatomy of the dorsal rootlets in the lower thoracic spinal cord.^20,30  DR fibers were placed in 1-mm intervals along the rostrocaudal axis (fig. 1B). Near the dorsal horn of the spinal cord, the DR fiber branched into a daughter fiber that traveled along the rostrocaudal axis within the DC.²² 

The third step in our model analysis was to assess the direct axonal response to KHFSCS. We performed this final step by applying the extracellular voltages (step 1) to the axon models of DC and DR fibers (step 2). We then calculated the stimulation amplitudes required for activation and conduction block under a variety of model parameters.

We calculated the minimum stimulation amplitudes required to generate action potentials (i.e., activation threshold) and to block the conduction or propagation of action potentials along the axon (i.e., block threshold) for an individual fiber. We defined thresholds as the peak amplitude for a single phase of the symmetric biphasic KHFSCS waveforms (fig. 2C).

The activation threshold was defined as the minimum peak amplitude required to generate one or more action potentials in a particular fiber (fig. 2C).

It is important to note that, unlike conventional (approximately 50 Hz) SCS, KHFSCS has the potential to block action potential conduction along an axon.⁸  The block threshold was defined as the minimum peak amplitude required to block the conduction of an action potential from one end of an axon to the opposite end. To calculate the block threshold, an internal stimulating electrode was placed at one end of an axon (i.e., most caudal node of a DC fiber) and was used to apply a test pulse after a KHFSCS waveform had been applied for 40 ms. Successful conduction block was defined as the condition in which no action potential propagated from the caudal end of the axon to the rostral end (fig. 2C).

We varied a number of model parameters to investigate their potential significance in KHFSCS. These parameters included waveform shape, dCSF layer thickness, lead location, fiber collateralization, and fiber size (Table 2). To examine the potential influence of each parameter, we calculated the activation and conduction block thresholds for DC and DR fibers for each parameter set.

We considered two KHFSCS waveforms in this study. The first KHFSCS waveform was a continuous 10-kHz sinusoidal waveform. The second KHFSCS waveform was a rectangular waveform with symmetric cathodic and anodic phases with a pulse width of 30 μs and an interpulse interval of 20 μs applied at a rate of 10 kHz. This rectangular KHFSCS waveform closely resembled the reported waveform parameters of the Senza device manufactured by Nevro, Inc., USA.¹⁰  We also calculated activation thresholds for a conventional SCS waveform with a pulse width of 210 μs applied at a rate of 50 Hz, which represented common clinical stimulation parameters.^31–33 

In this study, the dCSF layer thickness represented the distance between the dorsal surface of the spinal cord and the dural sac. KHFSCS thresholds were calculated for dCSF thicknesses of 2.0, 3.2, and 4.4 mm.

We defined lead location as the position of the lead relative to the dural surface and the spinal cord midline. To examine the effects of dorsoventral lead position on KHFSCS thresholds, we initially placed the lead on the dural surface at the spinal cord midline and moved it in the dorsal direction, that is, away from the dura. We calculated thresholds for distances of 0.0, 0.2, 0.4, 0.8, and 1.2 mm between the dura and the lead. To examine the effects of mediolateral lead position on KHFSCS thresholds, we again placed the lead at the spinal cord midline and moved it laterally along the surface of the dura. We calculated thresholds for lead positions offset 0.0, 1.0, 2.0, and 3.0 mm relative to the spinal cord midline.

We added collaterals to DC fibers to assess the effects of fiber collateralization or branching on KHFSCS thresholds. The branched DC fibers consisted of a parent fiber projecting along the rostrocaudal axis and daughter branches connected near the center of the parent fiber. The daughter collaterals were oriented perpendicular to the parent DC fiber and projected in the ventral direction. We calculated KHFSCS activation thresholds for DC fibers with a single collateral as well as DC fibers with multiple collaterals placed at adjacent nodes of Ranvier. We also examined multiple parent fiber-to-collateral diameter ratios (1.0, 1.6, and 2.0), multiple fiber diameters (5.7 to 16.0 μm), and KHFSCS with a monopolar stimulation configuration. The membrane surface area of branching nodes was increased to 150% relative to nonbranching nodes.¹⁹  The DC parent fibers had a length of approximately 120 mm with node spacing dependent on the axon diameter.²⁹  The collateral length and node spacing were also diameter dependent with lengths of 4.5, 4.5, and 3.8 mm for collateral diameters of 5.7, 7.3, and 11.5 μm, respectively.

We also varied the diameter of DC fibers to investigate the effects of fiber size on KHFSCS thresholds. We calculated activation and conduction block thresholds for three fiber diameters: 7.3, 11.5, and 15.0 μm.

Several clinical and experimental studies have investigated the effects of kilohertz frequency electrical stimulation on the nervous system with stimulation waveforms having different parameters (e.g., sinusoidal vs. rectangular; monophasic vs. biphasic; continuous vs. discontinuous; voltage controlled vs. current controlled).^{6–8,10,13–15,34} We elected to focus our analysis on the waveforms that we believed were most relevant to current clinical applications of kilohertz frequency stimulation for chronic pain applications. Therefore, we considered continuous current-controlled sinusoidal and rectangular KHFSCS waveforms to examine the potential effects of waveform shape on activation and block thresholds (fig. 3). Both KHFSCS waveforms were applied at a rate of 10 kHz. We also considered a conventional SCS waveform applied at a rate of 50 Hz.

The activation and conduction block thresholds for both KHFSCS waveforms and the activation thresholds for the conventional SCS waveform were calculated for all DC and DR fibers (fig. 3). At the activation threshold, axons did not fire continuously but typically only generated a number of action potentials during the first 30 ms or less of KHFSCS (fig. 2C). However, it was possible to generate continuous firing at a rate of several hundred hertz with a range of stimulus amplitudes between the activation threshold and block threshold (data not shown). At the block threshold, KHFSCS also produced transient spiking, an effect referred to as the “onset response” (fig. 2C).⁸ 

Figure 3A shows contour plots of the activation and block thresholds for DC fibers as a function of position within the DC. Figure 3B shows the activation thresholds for DR fibers as function of rostrocaudal level of the DR and daughter DC fiber branch point relative to the center of the cathode. The sinusoidal waveform produced a minimum activation threshold and a block threshold of 4.7 and 8.4 mA, respectively, for DC fibers. For DR fibers, the sinusoidal waveform produced minimum activation and block thresholds of 4.6 and 9.2 mA, respectively. The rectangular waveform produced minimum activation and block thresholds of 4.6 and 8.0 mA, respectively, for DC fibers. For DR fibers, the rectangular waveform generated minimum activation and block thresholds of 4.3 and 8.5 mA, respectively. The conventional SCS waveform produced minimum excitation thresholds of 1.2 and 1.6 mA for DC and DR fibers, respectively.

Both of the KHFSCS waveforms had significantly higher activation thresholds relative to the conventional SCS waveform. However, the rectangular KHFSCS waveform required slightly lower amplitudes for activation and conduction block relative to the sinusoidal KHFSCS waveform. For a given amplitude, the rectangular KHFSCS waveform with a 30-μs pulse width also injected less charge per phase relative to the continuous sinusoidal KHFSCS waveform (i.e., 6% difference). Therefore, we only present the rectangular KHFSCS waveform for the remaining analyses.

The dorsal aspect of the spinal cord is separated from the dura by a layer of cerebrospinal fluid (fig. 4A). The thickness of this dCSF layer is believed to be one of the most important variables in determining therapeutic stimulation amplitudes in SCS as well as the ability to stimulate the desired neural targets.³⁵  The thickness of the dCSF varies significantly as a function of spinal level and is also highly variable between patients.³⁶  The middle and lower thoracic spinal levels have a thicker dCSF layer than other spinal levels, and higher amplitudes are required for therapeutic stimulation.³⁵  The thicker dCSF layer also makes it more difficult to activate the DC fibers without generating unwanted side effects due to activation of DR fibers.³⁷  Therefore, we examined the effects of dCSF thickness on activation and block thresholds in KHFSCS (fig. 4). Three thicknesses (2.0, 3.2, and 4.4 mm) were examined that largely covered the estimated range of the dCSF layer for the lower thoracic spinal cord (i.e., 3.6 ± 1.6 mm at the T11 spinal level).³⁶ 

Increases in dCSF thickness required higher KHFSCS amplitudes for activation and block in both DC and DR fibers in agreement with the clinical experience with conventional (approximately 50 Hz) SCS (fig. 4). For DC fibers, the minimum activation thresholds were 2.2, 4.6, and 8.2 mA and block thresholds were 3.9, 8.0, and 14.3 mA for dCSF thicknesses of 2.0, 3.2, and 4.4 mm (fig. 4C). For DR fibers, the minimum activation thresholds were 2.7, 4.3, and 6.4 mA and block thresholds were 5.5, 8.5, and 11.2 mA for dCSF thicknesses of 2.0, 3.2, and 4.4 mm (fig. 4C). For a thin dCSF layer (2.0 mm), there was a significant increase in the likelihood of DC and DR fiber activation, and it was possible to achieve conduction block of large-diameter DC fibers within the clinically relevant range (0.5 to 5 mA).¹⁰ 

This study considered KHFSCS with cylindrical leads that are typically implanted percutaneously with a Tuohy-style needle. When implanting these leads, much care is taken to place the leads at the desired rostrocaudal and mediolateral locations. However, it is difficult to control the dorsoventral position of the lead. In addition to the dorsoventral variability in lead location, it is also possible for the lead to be placed or migrate lateral to the spinal cord midline.

Because of the potential difficulties in controlling lead position, we examined the effects of electrode location on the activation thresholds for DC and DR fibers (fig. 5). To examine the effects of the dorsoventral position of the lead, we calculated the activation thresholds at a range of distances between the dural surface and the lead surface (fig. 5A). For DC fibers, minimum activation thresholds were 4.6, 5.6, 6.3, 7.6, and 8.8 mA for distances of 0.0, 0.2, 0.4, 0.8, and 1.2 mm, respectively. For DR fibers, minimum activation thresholds were 4.3, 4.9, 5.3, 6.0, and 6.7 mA for distances of 0.0, 0.2, 0.4, 0.8, and 1.2 mm, respectively.

We also examined the effects of lateral lead migration by calculating the activation thresholds for a range of lateral lead offsets relative to the spinal cord midline (fig. 5B). For DC fibers, the minimum activation thresholds were 4.6, 4.7, 5.1, and 5.8 mA for distances of 0.0, 1.0, 2.0, and 3.0 mm. For DR fibers, the minimum activation thresholds were 4.3, 3.6, 3.3, and 3.3 mA at distances of 0.0, 1.0, 2.0, and 3.0 mm.

This analysis demonstrated that dorsoventral movement of the lead away from the dural surface resulted in an exponential increase in the activation thresholds for both DC and DR fibers. However, activation thresholds increased more rapidly for DC fibers relative to DR fibers. Lateral movement of the lead produced an exponential increase in the activation thresholds for DC fibers and an exponential decrease in the activation thresholds for DR fibers. Therefore, any movement of the stimulating lead away from the surface of the dura and/or the spinal cord midline increased the selective stimulation of DR fibers over DC fibers. These trends matched previous clinical and modeling results of conventional SCS.^37–39 

Dorsal column fibers are not simple straight axons, but they send out several collaterals at regularly spaced intervals that project into the gray matter of the spinal cord where they freely arborize.^19,40  These projections occur at the nodes of Ranvier and can occur at an average spacing of approximately 1 mm, depending on the fiber type, fiber diameter, and distance from the corresponding DR fiber.^19,40  Modeling studies have shown that fiber collaterals lower conventional SCS activation thresholds by as much as 50%.¹⁹  Therefore, we examined the effects of DC fiber collateralization on the activation thresholds during KHFSCS. We incorporated DC fiber collaterals that were oriented perpendicular to the rostrocaudal orientation of the parent DC fiber and projected in the ventral direction.

To examine the effects of DC fiber collateralization on activation thresholds during KHFSCS, the location of the lowest threshold DC fiber (i.e., superficial DC fiber at the spinal cord midline) was varied along the rostrocaudal axis relative to the stimulating electrodes (fig. 6A). For an 11.5-μm diameter DC fiber, the collateral diameter was varied among 5.7, 7.3, and 11.5 μm (i.e., fiber-to-collateral diameter ratios of 2.0, 1.6, and 1.0). The average fiber diameter-to-collateral diameter ratio has been estimated to be on the order of 3.1 ± 0.7 in Ia and Ib primary afferents originating from the cat hind limb muscles.^19,41  Smaller ratios (i.e., larger collateral diameters) will produce a larger decrease in activation thresholds. Therefore, the range examined in this study represents an upper limit on activation threshold reduction due to DC fiber collateralization. Collateralization produced a maximum decrease in activation threshold of 0.0, 4.2, and 29.8% for fiber collateral diameters of 5.7, 7.3, and 11.5 μm, respectively (fig. 6A). Although Figure 6A shows that there was a significant reduction in the activation thresholds for a fiber-to-collateral diameter ratio of 1.0, the more physiologically relevant ratio of 2.0 (i.e., collateral diameter of 5.7 μm) produced no decrease in activation threshold for a bipolar stimulation configuration.

The minimum thresholds occurred when the branching nodes were centered at the rostrocaudal level of the stimulating electrodes. There was an increase in activation threshold when the branching node was between the rostrocaudal level of the anode and cathode where the depolarizing force of the extracellular stimulation (i.e., second spatial derivative of the extracellular voltages or activating function⁴² ) was close to zero (figs. 2A and 6A). Although clinical KHFSCS is performed with bipolar stimulation, we performed the same analysis for a monopolar stimulation configuration with a distant ground. Even under monopolar stimulation configurations, fiber collateralization produced a maximum threshold reduction from 5.7 to 5.3 mA (7.0%) for a fiber-to-collateral diameter ratio of 2.0 (data not shown).

The preceding analysis within the current section, “Fiber Collateralization,” only considered a DC fiber with a single collateral; however, the presence of multiple collaterals can further decrease the activation thresholds.^19,43  Therefore, we also calculated the activation threshold during KHFSCS for a (11.5 μm diameter) DC parent fiber with 11 (5.7 μm diameter) collaterals at adjacent nodes of Ranvier. Multiple collaterals produced trends similar to those shown in figure 6A. For bipolar stimulation, there was no decrease in the activation thresholds, but a significant increase when the branching nodes were near the rostrocaudal levels of the stimulating electrodes (data not shown). For monopolar stimulation, the presence of these 11 collaterals produced a maximum threshold reduction from 5.7 to 5.2 mA (8.8%) (data not shown).

These small reductions in KHFSCS activation thresholds do not match previous conventional SCS modeling results that demonstrated a major reduction in DC activation threshold due to fiber collateralization for small-diameter fibers (6 μm) with smaller diameter collaterals (2 μm).¹⁹  However, studies have shown that conventional SCS directly affects large-diameter fibers due to their low activation thresholds relative to small-diameter fibers.²¹  Therefore, we investigated the effect of fiber collateralization for a range of fiber diameters (5.7 to 16.0 μm) with a fiber-to-collateral diameter ratio of 1.0 (i.e., fiber diameter = collateral diameter). A single fiber collateral produced a major reduction in activation threshold for small-diameter fibers, but this effect decreased with increasing fiber diameter (fig. 6, B and C).

Extracellular electrical stimulation can excite myelinated axons by generating action potentials at the nodes of Ranvier. For myelinated fibers, the activation threshold is largely determined by the spacing between adjacent nodes of Ranvier.⁴⁴  This internodal spacing increases as a function of axon diameter, and therefore, large-diameter fibers have a lower threshold than smaller fibers. Previous studies suggest that conventional SCS functions through direct activation of large-diameter myelinated axons within the DC.²¹  Therefore, fiber diameter is an important variable to consider in KHFSCS because the DC of the human spinal cord consists of axons with a wide range of diameters (average fiber diameter approximately 5.0 μm and a maximum diameter of 16.0 μm at lower thoracic levels).²⁸  The density and distribution of fiber sizes are also a function of mediolateral position within the DC, with higher densities of medium-diameter and large-diameter fibers in the lateral DCs.²⁸  Therefore, activation and block thresholds were calculated for three DC fiber diameters (7.3, 11.5, and 15.0 μm) that represented a wide range of fiber diameters found within the human spinal cord and included large-diameter fibers that are most likely to be affected by SCS (fig. 7).^21,28  The results showed the expected trend of large-diameter fibers having the lowest activation and block thresholds. The minimum activation thresholds were 10.8, 4.6, and 3.1 mA for DC fiber diameters of 7.3, 11.5, and 15.0 μm, respectively (fig. 7B). The minimum block thresholds were 19.1, 8.0, and 5.5 mA for DC fiber diameters of 7.3, 11.5, and 15.0 μm, respectively (fig. 7B). The minimum DR fiber activation and block thresholds of 4.3 and 8.5 mA are also shown in figure 7B. This analysis shows that it is possible to have significant direct activation of large-diameter fibers (≥11.5 μm) within the current clinical range of stimulation amplitudes (0.5 to 5 mA).¹⁰  However, even for the largest diameter fibers (15.0 μm), conduction block was not possible within the clinical range.

Kilohertz frequency spinal cord stimulation for chronic pain management is a new and promising technology; however, recent clinical studies have presented conflicting results. Although it is possible that outcome inconsistency stems from clinical trial design and patient selection, it is also possible that clinical outcomes may have been affected by technical limitations, such as lead positioning and stimulation programming choices, and differences in the stimulation waveform parameters (e.g., pulse width, frequency). Addressing such limitations can be challenging, particularly for KHFSCS, because it does not produce paresthesias.^9,11,12  Understanding the mechanisms responsible for the clinical benefits (or failure) of KHFSCS is critical to reduce variability in how the therapy is delivered to the target neural elements. This knowledge will allow for future clinical studies to be more reproducible and to definitively answer efficacy questions. The present data represent a first step in this process and may be useful in KHFSCS therapy optimization.

This study used a computer model to characterize the effects of KHFSCS on spinal cord axons. Direct activation of spinal cord axons is one of the most obvious potential therapeutic mechanisms of KHFSCS. Therefore, we computed the KHFSCS thresholds for action potential generation under a number of conditions. In general, the results of this study suggest that KHFSCS requires significantly higher amplitudes for excitation relative to conventional approximately 50 Hz SCS (fig. 3). These amplitudes were often in the upper end or outside of the range currently used in clinical practice (i.e., 0.5 to 5 mA).^10,12 

Regarding DC fibers, it was possible to excite large-diameter fibers (≥11.5 μm) within the clinical range, especially for a thin dCSF layer (i.e., 2 mm) (fig. 4). However, any movement of the lead from the “ideal” location (i.e., spinal cord midline and adjacent to the dural surface) resulted in an increase in the activation threshold (fig. 5). DC fiber collaterals have been shown to produce significant reductions in the stimulation amplitudes required for the activation of small-diameter fibers with conventional SCS.¹⁹  In the current study, fiber collateralization of large-diameter DC fibers only produced a minor reduction in KHFSCS activation thresholds and increased activation thresholds at certain rostrocaudal levels for bipolar stimulation (fig. 6).

With respect to DR fibers, it was possible to excite large-diameter (15 μm) DR fibers within the clinical range of stimulus amplitudes, especially with a thin dCSF layer (2 mm) (fig. 4). As with DC fibers, dorsoventral movement of the lead away from the dural surface led to an increase in the DR fiber activation thresholds outside of the clinical range (fig. 5A). However, lateral displacement of the lead produced a moderate reduction in the DR fiber activation threshold (fig. 5B) and higher selectivity for stimulating DR fibers over DC fibers.

The data show that within the clinical range (0.5 to 5 mA), it is possible that KHFSCS leads to direct activation of large-diameter (≥11.5 μm) fibers in the superficial layer of the DC as well as large-diameter DR fibers close to the cathode or anode. Therefore, direct activation of spinal cord fibers may be a potential mechanism of action of KHFSCS. However, it is important to consider that the minimum activation thresholds (i.e., most excitable fibers) were within the upper limits of this clinical range and above the clinically effective amplitudes reported in a small cohort of patients.⁴⁵  These coupled results reduce the likelihood that the clinical benefit thus far reported for KHFSCS^9–12  stems from direct activation of large-diameter spinal cord fibers. A lack of activation is also consistent with clinical evidence showing that, within the clinically effective amplitude range of 0.5 to 5 mA, patients do not experience the paresthesias that are a hallmark of DC or DR fiber activation.^10–12  Paresthesias are only generated at significantly higher amplitudes.⁴⁵  In turn, our theoretical predictions are in line with prior experimental evidence showing that KHFSCS does not modulate the firing rate of gracile nucleus neurons in a rat model of neuropathic pain.¹⁵ 

Direct conduction block of spinal cord fibers may be the most commonly assumed mechanism to mediate the effects of KHFSCS. This notion stems from prior studies that demonstrate conduction block generated by kilohertz frequency electrical stimulation in peripheral nerve models^6–8  as well as the lack of paresthesias reported in clinical studies.^9–12  We tested the thresholds for conduction block in DC and DR fibers under a number of different conditions. Conduction block thresholds were always higher than activation thresholds and this relation matched previous studies in the peripheral nervous system.⁴⁶  Conduction block thresholds were almost always outside of the clinical amplitude range, the only exception being large-diameter DC fibers when a very thin dCSF layer was assumed (2 mm) (fig. 3). However, the main current clinical application of KHFSCS involves stimulation at the lower thoracic (T8 to T11) levels which typically present a thicker dCSF layer.³⁵  The prediction that KHFSCS does not generate direct conduction block in the spinal cord matches prior observations. Before axonal conduction block is generated with kilohertz frequency stimulation, there is an initial increase in action potential firing in the target axons, called the “onset response.”^8,46  This onset response can be observed experimentally during KHFSCS by recording increased activity in wide dynamic range neurons³⁴  and manifests behaviorally as several signs of discomfort for the first few minutes of stimulation that are different from the signs of sensory threshold.¹⁴  Although this onset response has been observed in animal models of KHFSCS, no paresthesias or other subjective perceptions have been reported during clinically effective KHFSCS in human patients.¹⁰ 

Although it is well accepted that conventional SCS functions through direct activation of spinal cord fibers,^47,48  the potential mechanisms of action of KHFSCS are unknown. The results of this study suggest that KHFSCS may not function explicitly through direct activation or conduction block of spinal cord fibers and may function through more complex or subtle mechanisms. There are a number of other potential mechanisms through which KHFSCS could provide pain relief. Some of these mechanisms include pseudospontaneous activation, transmission delays, or conduction failure at branch points. Pseudospontaneous activity is related to the concept that an individual neuron’s response to a near-threshold stimulus may be subthreshold or suprathreshold due to the stochastic nature of ion channel gating.⁴⁹  Therefore, KHFSCS at subthreshold or near-threshold amplitudes may lead to asynchronous or pseudospontaneous activation of afferent fibers. KHFSCS may also alter the axonal response to incoming signals¹⁴  or produce action potential conduction failure at axon branch points.^{14,19,50–52}  KHFSCS could potentially induce ion accumulation in the extracellular or periaxonal space that could affect activation and/or conduction block thresholds.^53,54  Last, paresthesia-free analgesia may not require KHFSCS, as it can also be achieved with conventional SCS waveforms applied in short bursts at 40 Hz.⁵⁵ 

Although the KHFSCS model presented in this study was based on standard computational modeling principles, it was subject to a number of potential limitations, and its accuracy needs to be confirmed with future experimental and clinical studies. We only considered frequency-independent electrostatic solutions to the electric field generated during KHFSCS. It is possible that the frequency-dependent properties of the electrode–electrolyte interface of the stimulating electrodes as well as the surrounding biological tissue could affect the electric field generated during KHFSCS. We did not account for the potential effects of electrode encapsulation that could alter the electric field generated during KHFSCS.^22,56  However, the FEM did produce monopolar electrode impedances (430 Ω) that were within the clinically observed range for chronically implanted percutaneous SCS leads and activation thresholds with conventional 50 Hz SCS that were similar to a previously published computer model of SCS.²² 

We used a nonlinear axon model derived from mammalian motor axons²⁹  that may have a limited ability to represent the behavior of sensory neurons in the spinal cord. The axon model also was not explicitly parameterized for electrical stimulation rates in the kilohertz range. However, this axon model has shown a high degree of accuracy in reproducing the in vivo axonal response to kilohertz frequency electrical stimulation in the peripheral nervous system.⁴⁶ 

Finally, this study represents a first step to describe the mechanisms that are more likely to mediate the clinical effects of KHFSCS. Future studies will be needed, and greater complexity will be added to include detailed axon branching, synaptic terminals, cell bodies, and dendrites. Future work may examine the network effects of KHFSCS.^57,58  Future modeling studies could also be used to investigate optimized electrode designs, stimulation configurations, and waveform parameters that lower excitation/block thresholds and/or improve the ability to affect specific neural targets.

The results presented in this study represent significant strides toward theoretical characterization of the potential pain relief mechanisms of KHFSCS. Although a number of variables were considered, the activation and block thresholds for the most excitable fibers were in the upper limits or outside of the current clinical range of KHFSCS. Therefore, this study suggests that clinically effective KHFSCS may not function explicitly through direct activation or conduction block of spinal cord fibers and alternative concepts should be explored and evaluated.

The authors thank John T. Gale, Ph.D. (Neurosciences, Lerner Research Institute, Cleveland Clinic, Cleveland, Ohio), Hyun-Joo Park, Ph.D. (Center for Neurological Restoration, Cleveland Clinic), and Aaron J. Fleischman, Ph.D. (Biomedical Engineering, Lerner Research Institute, Cleveland Clinic), for their assistance with the finite element modeling. The authors also thank Niloy Bhadra, M.D., Ph.D. (Biomedical Engineering, Case Western Reserve University, Cleveland, Ohio), and Narendra Bhadra, Ph.D. (Neural Engineering Center, Case Western Reserve University), for their helpful discussions with kilohertz frequency electrical stimulation of neural tissue, and Kabilar Gunalan, M.S. (Biomedical Engineering, Case Western Reserve University), for his assistance running the NEURON simulations. The authors also thank Michael A. Moffitt, Ph.D. (Neuromodulation Research and Advanced Concepts, Boston Scientific Neuromodulation, Valencia, California), for his helpful comments and discussion.

This work was supported by the Louis Stokes Cleveland Veterans Affairs Medical Center, Cleveland, Ohio.

Dr. Kilgore has equity ownership in Neuros Medical, Inc. (Willoughby Hills, Ohio). Dr. McIntyre is a paid consultant for Boston Scientific Neuromodulation (Valencia, California) and is a shareholder in the following companies: Surgical Information Sciences, Inc. (Minneapolis, Minnesota); Autonomic Technologies, Inc. (Redwood City, California); Cardionomics, Inc. (Forest Lake, Minnesota); Neuros Medical, Inc.; and Enspire DBS Therapy, Inc. (Wilmington, Delaware). Dr. Machado has the following disclosures: distribution rights from intellectual property for Enspire DBS Therapy, Inc., Autonomic Technologies, Inc., and Cardionomics, Inc.; consultant for Icahn School of Medicine at Mount Sinai (New York, New York), Spinal Modulation (Menlo Park, California), and Functional Neuromodulation (St. Paul, Minnesota); and fellowship program support from Medtronic, Inc. (Minneapolis, Minnesota). Dr. Lempka declares no competing interests.

Unless specified otherwise, the dorsal cerebrospinal fluid layer had a thickness of 3.2 mm, a value within the range clinically observed at the lower thoracic levels.³⁶  The dura thickness was set to 300 μm, and its dorsal surface was flattened for computational simplicity.^22,59  Each electrode contact had a diameter of 1.25 mm and length of 3 mm and were separated from adjacent electrodes by electrode insulation 1 mm in length.

The model geometry was defined and meshed with the 3-matic Module within the Mimics Innovation Suite (Materialise, Belgium). We specified higher mesh densities near the electrode array as well as within a 35-mm long region of interest surrounding the electrode array. The total model length was 201 mm with a diameter of 70 mm. The FEM consisted of more than 12.7 million first-order tetrahedral elements. After the model geometry and mesh were generated, it was exported to the FEM software package, COMSOL Multiphysics (COMSOL, Inc., U.S.A.). Within COMSOL, electrical conductivities were assigned to each domain based on experimental data available in the literature (Table 1).^21,22,26 

In this study, all simulations were performed for bipolar stimulation applied through two active contacts separated by an inactive contact (i.e., separation of 8 mm center-to-center) (fig. 1A). We selected a distance of 8 mm between the anode and cathode because it represents common clinical programming selections in spinal cord stimulation. To calculate the electric fields generated by bipolar KHFSCS, we placed unit current sources (i.e., 1 A) of opposite polarity at the cathode and anode and set the outer surfaces of the general thorax layer to ground (i.e., 0 V) (fig 1A). We then calculated the voltage distributions (Φ) generated in the tissue (stiffness matrix, σ) based on the specified current sources (I) by solving the Poisson equation:

We calculated electrostatic FEM solutions for these unit current sources with an iterative equation solver using the conjugate gradient method. We refined the mesh density until further increasing the mesh density produced a maximum of less than 2% difference in the activation thresholds calculated for the neural elements considered in this study. Doubling the total volume of the FEM produced a maximum of less than 2% difference in the predicted activation thresholds.

We represented dorsal column (DC) and dorsal root (DR) fibers with a previously published compartmental model of a mammalian motor axon.²⁹  This model reproduces experimental data by accurately representing the ion channels at the nodes of Ranvier as well as matching the geometry of the paranode, internode, and myelin to measured morphology. This model incorporates a double-layer cable model that accounts for the finite impedance of the myelin sheath. The nodes of Ranvier contain fast Na⁺, persistent Na⁺, and slow K⁺ nonlinear conductances as well as the linear leakage conductance and the membrane capacitance. All equations and parameters for the axon model were defined in the study by McIntyre et al.²⁹  Unless specified otherwise, each DC fiber was 159 mm long with a fiber diameter of 11.5 μm and 128 nodes of Ranvier. For the DR fibers, the mother fiber had a diameter of 15.0 μm and a length of 43.5 mm with 31 nodes of Ranvier. The daughter DC fiber had a diameter of 11.5 μm and a length of 118 mm with 95 nodes of Ranvier.

The lengths of the DC and DR fibers were sufficient to ignore potential edge effects that can occur due to the wide mean membrane depolarization that occurs during kilohertz frequency stimulation.⁴⁶  Simulations were performed with the software package, NEURON, within the Python programming environment.⁶⁰  Model solutions were calculated using backward Euler implicit integration with a time step of 0.002 ms.

To assess the direct neural response to KHFSCS, the voltage distributions calculated in the FEM were ported to the Python programming environment and directly applied to the axon models of the DC and DR fibers. Because the bulk conductivity (σ) is linear, the voltage distributions generated by different current source waveforms or magnitudes were scaled versions of the original FEM solutions with unit currents. The scaled voltage distributions were interpolated onto the model neurons described above using the extracellular mechanism within NEURON. The activation and block thresholds were determined using a bisection algorithm (error < 0.1 mA).