If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Department of Biomedical Engineering, Duke University, Durham, NC, USADepartment of Electrical and Computer Engineering, Duke University, Durham, NC, USADepartment of Neurobiology, Duke University School of Medicine, Durham, NC, USADepartment of Neurosurgery, Duke University School of Medicine, Durham, NC, USA
Novel method for computational design of non-regular temporal patterns of stimulation.
Novel patterns show equivalent efficacy compared to conventional 50 Hz stimulation and improve efficiency.
Patterns evaluated in both healthy and nerve-injured animals.
Temporal patterns of stimulation represent a novel dimension for improving the efficacy of spinal cord stimulation to treat chronic neuropathic pain.
We hypothesized that nonregular temporal patterns of stimulation designed using a computational model would be superior to conventional stimulation at constant frequencies or completely random patterns of stimulation.
Using a computational model of the dorsal horn network and an optimization algorithm based on biological evolution, we designed an optimized pattern of spinal cord stimulation with comparable efficacy and increased efficiency relative to constant frequency (CF) stimulation. We evaluated the effect of different temporal patterns on individual neurons recorded in the dorsal horn of urethane-anesthetized rats.
The optimized pattern and 50 Hz CF stimulation produced greater inhibition of spontaneously firing neurons recorded in vivo than random 50 Hz stimulation or a pattern designed intentionally with poor fitness. Spinal Cord Stimulation (SCS) led to significant changes in the firing patterns of recorded units, and stimulation patterns that generated significant inhibition also tended to reduce entropy and regularize the firing patterns of units, suggesting that patterns of dorsal horn neuron activity may be important for pain perception in addition to the firing rate.
These results demonstrate that the computational model can be used as a tool for optimizing stimulation parameters and suggest that optimized temporal patterns may increase the efficacy of spinal cord stimulation.
]. Spinal cord stimulation (SCS) is an effective and cost-effective option for treating chronic pain conditions that are refractory to other management options such as opioids, physical therapy, surgical interventions, or radio-frequency ablation [
]. The large parameter space limits the number of parameter combinations that can be tested during a session, and programmers typically limit selection of stimulation frequency and pulse durations to commonly used combinations [
Conventional SCS uses constant frequency (CF) stimulation (i.e., regular temporal patterns), but non-regular temporal patterns of stimulation may be more effective and more efficient than CF stimulation [
]. Mechanosensory neural circuits respond to both the rate and timing of sensory inputs, with each component coding unique characteristics. For example, mechanoreceptive afferents encode features such as the pressure, orientation, or vibration through the specific timing of spikes [
]; therefore, it is plausible that the timing of electrical pulses may contribute as much as the pulse rate to the efficacy of SCS. A recent pre-clinical study demonstrated that stochastically varying SCS pulse times could achieve comparable effects to CF SCS [
]. However, commercial development of novel modalities of SCS has focused primarily on higher stimulation frequencies. For example, burst (500 Hz intraburst, 40 Hz interburst) and kilohertz frequency stimulation have shown similar or superior reductions in pain compared to conventional stimulation [
], suggesting that further improvement of non-regular SCS waveforms is possible.
Given the importance of temporal dynamics in the processing of sensory inputs, we hypothesized that optimized, non-regular temporal patterns of SCS would have equivalent or superior efficacy while potentially improving energy efficiency compared to CF SCS. We used a modified genetic algorithm [
], an optimization approach that mimics biological evolution, to design temporal patterns of SCS in a validated computational model of the dorsal horn network. We quantified the effects of individual model-optimized patterns on the activity of dorsal horn neurons recorded in vivo in healthy rats and in a rat model of neuropathic pain. An optimized temporal pattern of SCS produced effects on net inhibition and excitation of dorsal horn neurons comparable to 50 Hz SCS. However, the optimized pattern improved waveform energy efficiency and may reduce the variance in responses to SCS across heterogeneous pain states. Additionally, prior studies have focused on changes in the firing rate of dorsal horn neurons, but we also identified changes in the patterns of activity in neurons inhibited by SCS.
2. Materials and methods
2.1 Computational optimization
The computational model used for optimization has been described previously [
]. Briefly, the model consisted of three types of neurons: an inhibitory interneuron (IN), excitatory interneuron (EX), and wide-dynamic range (WDR) projection neuron. The firing rate of the WDR neuron represented the model analog for pain because WDR firing rate is correlated with the amount of pain [
]. Each node within the model represented a distinct sensory receptive field area and contained one IN, one EX, and one WDR neuron, and the inputs to each node were the spike times from 15 Aβ fibers, 15 Aδ fibers, and 30 C fibers. Ongoing baseline activity in afferent fibers was simulated by generating spike times from a homogeneous Poisson distribution with mean firing rates based on recordings following peripheral nerve injury (Aβ & Aδ fibers: μ = 2.2 Hz, C: μ = 1.5 Hz) [
]. The synaptic connections within an individual node were based on dorsal horn recordings and tuned so that model activity matched observed dorsal horn neuron behavior such as wind-up and A-fiber mediated inhibition [
]; in their recordings, they identified a “center” receptive field area wherein low-threshold peripheral stimulation produced excitation of recorded deep dorsal horn neurons and a “surround” receptive field area that wherein stimulation produced mixed excitation and inhibition of recorded neurons. The synaptic connections in the model were optimized to replicate this center-surround architecture with zone 1 representing the center receptive field area and zones 2 and 3 representing the surround receptive field area. Excitatory and inhibitory connections between nodes came from interneurons to projection neurons in the surrounding nodes.
Neuropathic pain leads to disinhibition in the dorsal horn through a variety of mechanisms, and progression of pain may contribute to decreased efficacy of SCS over time [
]. The default model state was the computational network model with all synaptic conductances left unchanged, and we generated model pain states by varying connections influencing GABAergic inhibition. First, we reduced the weight of synaptic inputs from all Aβ fibers in the model by 50%. Second, we reduced the weight of all GABAergic synapses from IN neurons to WDR neurons by 50%. When we implemented reduced GABAergic inhibition or reduced conductance from Aβ fiber inputs in the dorsal horn, the responses of model WDR neurons to CF stimulation shifted compared to the default network state, and these shifts matched prior model responses [
]. Genetic algorithms function analogously to natural selection in a biological system: stimulation patterns are scored according to “fitness” quantified by a cost function, and the most effective patterns are used as “parents” to generate offspring patterns. We used all improvements to the algorithm described in Cassar et al. [
] to enable rapid convergence and maximize performance. The algorithm created individual stimulation patterns with one bit per millisecond (i.e., “1” for a pulse and “0” for no pulse at each millisecond). Each node in the model received the temporal pattern of stimulation through the Aβ fiber inputs including a delay representing the propagation time from the stimulation site. However, we varied the percentage of fibers that received the pattern to match the percentage of C and Aδ fibers that were active in each node. 100% of fibers were active for zone 1 nodes, 75% were active for zone 2 nodes, and 50% were active for zone 3 nodes. The distribution of active fibers was estimated based on the clinical assumption of optimal pain relief with complete paresthesia overlap during SCS [
]. Temporal patterns were evaluated over a 10 s window in the model following 8 s with no SCS. Wind-up in the WDR neurons occurred within this initial 8 s window and increased the number of spikes of the WDR neuron from ∼45 to ∼55 spikes/second (Hz). We initialized the GA with 50 different randomly generated temporal patterns (“organisms”).
In the computational model, the primary metric for scoring patterns was a proxy for stimulation efficacy: how much the pattern reduced model WDR neuron firing rate during stimulation, because WDR firing rate is correlated with pain [
]. A secondary scoring measure was efficiency, quantified as the average pulse frequency of the pattern, because decreased pulse repetition rate corresponds to less energy consumption, which, in turn, corresponds to improved device battery life. Third, we directly compared the response of each temporal pattern to the performance of CF stimulation with the same average pulse frequency. Further, we implemented multiple neuropathic pain states in the computational model representing different indications of SCS [
] and scored stimulation patterns favorably if they reduced the variance of WDR firing rates across the pain states. For pain states in the model, we used multiple instantiations of the computational model with varying synaptic strengths to identify patterns that performed well across heterogeneous conditions. The fitness of each pattern was evaluated using a customized scoring function (Equation (1)) for each of nine different model states: the baseline model state and then a continuum of eight model states representing reduced GABAergic inhibition and Aβ inputs as described above. The net fitness of each organism was calculated as the mean rank of the fitness across the simulated model states (default and pain states). We used rank ordering because scores were not normally distributed across model states although individual components were normally distributed within states. Each new generation consisted of four types of patterns: the best patterns from the previous generation (“elites”, 5 organisms, based on highest net fitness), patterns with mixed components of the best patterns from the previous generations (“children”, 25 organisms), new random patterns (“random immigrants”, 5 organisms), and newly designed patterns based on the temporal features of the previous generation (“competitive immigrants”, 15 organisms). The entropy of stimulation patterns was calculated using binned interpulse intervals with logarithmically sampled bins [
]. All animal procedures were approved by the Duke University Institutional Animal Care and Use Committee. Briefly, we recorded populations of individual dorsal horn (DH) neurons from urethane-anesthetized male Sprague-Dawley rats weighing between 300 and 500 g. We recorded from a total of eight rats, four animals that had previously undergone a spared nerve injury procedure, and four non-injured animals. For the spared nerve injury procedure, animals were anesthetized using a combination of isoflurane (4% for induction, 2% for maintenance, inhaled, Abbott Laboratories) and carprofen (subcutaneous at incision site, 5.0 mg/kg of 50 mg/ml solution, Sigma Aldrich). After ensuring rats were sufficiently anesthetized by checking for a pinch response, we cleaned the incision site with iodine and 70% isopropyl alcohol, exposed the sciatic nerve at the terminal branch point, ligated the tibial and peroneal branch using 5-0 silk suture, and resected a ∼3 mm section distal to the ligation [
]. We sutured the muscle and skin closed around the incision and monitored animals to ensure recovery. We performed behavioral assessments using von Frey filaments to quantify paw withdrawal thresholds prior to the surgery and on day 3, 5, 7, and 14 following the surgery [
]. Rats were placed in a chamber and allowed to acclimate to the chamber for 20 min prior to behavioral testing.
For acute neural recordings, we induced anesthesia with isoflurane (3%, inhaled) and maintained anesthesia with urethane (1.2 g/kg subcutaneous injection in the shoulder, Sigma Aldrich). We ensured rats were anesthetized by checking for a withdrawal reflex to pinching the hindpaw and applied a second (0.4 g/kg, i.p.) and third (0.1 g/kg i.p.) dose of urethane as necessary. We performed a tracheotomy to maintain respiration throughout the procedure. Animals were placed on top of a heating blanket and fixed using a stereotaxic frame. We monitored heart rate, temperature, and respiration throughout the procedure (Gaymar T/Pump and PhysioSuite; Kent Scientific, Torrington, CT). We exposed the sciatic nerve and placed a cuff electrode around the nerve proximal to the terminal branch point. We performed a laminectomy from the L1 to T13 spinal segments and resected the dura over the exposed spinal cord. We placed a custom four contact electrode with 0.6 mm × 0.6 mm platinum iridium contacts underneath the vertebra rostral to the exposed spinal cord between T10 and T11 (Dyconex, Zurich Switzerland). We determined motor threshold using constant 50 Hz, 200 μs stimulation through the ipsilateral contact pairs by finding the minimum stimulation amplitude that led to a perceptible motor response. After determining motor thresholds, we injected gallamine triethiodide (1 ml/h at 0.02 g/ml i.p. infusion, Santa Cruz Biotechnology, Dallas, TX) to paralyze animals and waited 1 h prior to beginning recordings.
We recorded single unit activity using 16- or 32- contact single shank electrodes (NeuroNexus, Ann Arbor, MI). Electrodes were inserted along the exposed spinal cord medial to the dorsal root entry zone at a 30–45° angle with respect to the coronal plane. Individual units were identified and recorded using a 32-channel MAP data acquisition system (Plexon, Dallas, TX). We used brush as a search stimulus to identify units within the receptive field. Units were identified online and then confirmed in offline analysis (Offline Sorter v3, Plexon, Dallas, TX) using various features (principal components, waveform energy, waveform amplitude). We evaluated responses to brushing (light pressure with a paintbrush), pressing (mild pressure with forceps), pinching (moderate pressure with an arterial clip), and crushing (heavy pressure with forceps) before and after electrical stimulation. Neurons were classified as wide-dynamic range (WDR) if they exhibited an increase in activity to each successive mechanical input. We evaluated responses to 30 s of SCS for each temporal pattern applied at 60–70% motor threshold and randomized the order patterns were applied across trials. We waited at least 30 s between each new stimulation pattern. For some blocks, we stimulated the sciatic nerve at C-fiber thresholds (50-60x motor threshold) to drive activity in dorsal horn neurons, and then simultaneously stimulated the dorsal columns with different SCS patterns. Following experiments, animals were euthanized with Euthasol (0.5 mL intraperitoneal injection, Virbac Corp.) followed by a bilateral thoracotomy.
2.3 Data analysis
Statistical tests, including Chi-squared tests, ANOVA, and individual pairwise post-hoc comparisons were performed using MATLAB and JMP Pro 14 (SAS, Cary, NC). P-values of p < 0.05 were deemed statistically significant unless otherwise stated. We identified neurons that were responders to stimulation (either excited or inhibited) by comparing activity during stimulation to spontaneous activity or activity during sciatic nerve stimulation alone. We calculated post-stimulus time histograms (PSTHs) for the response to the first pulse of each stimulus (i.e., assuming a frequency of 1 Hz because each pattern repeated every 1 s) after blanking the stimulus artifact. We determined the optimal number of bins based on an optimization function calculated from the spike count and spike variance based on a prior study [
]. We sampled spontaneous activity over a 30 s window preceding the trial using the same bin width and artificially blanked windows for the stimulus artifact. Neurons were labeled as responders if their activity was significantly different than spontaneous activity (z > 1.96) for three consecutive bins.
Entropy and unit entropy were calculated according to previously described methods [
]. Entropy was calculated with one dimensional ISIs (Equation (7)) and unit entropy was calculated according to the two-dimensional ISI probability distribution (Equation (8)). Correlations between neurons were calculated using the spike time tiling coefficient, a bounded measure of correlation that is independent of firing rate [
3.1 Designing temporal patterns of spinal cord stimulation using a computational model
We used a validated computational model of the dorsal horn network to quantify the effects of temporal patterns of stimulation on model neuron activity (Fig. 1a and b). A variety of optimization strategies exist for highly complex, multidimensional problems, and we used a genetic algorithm implemented specifically for design of temporal patterns of stimulation (Fig. 1c) [
]. The fitness of each pattern was evaluated using a cost function that combined four components of performance (see Equation (1), Methods). Across multiple iterations of the GA to test the robustness, the patterns with the highest fitness had a higher frequency component with instantaneous frequencies between 80 and 100 Hz and a lower frequency component with instantaneous frequencies between 10 and 20 Hz. We saw dramatic improvement in the fitness of each component of the cost function over the first five generations (Fig. 2b), and we did not observe substantial alterations in the GA pattern after the 50th generation (Supplemental Fig. 1a). The GA converged on temporal patterns with similar distributions of IPIs whether the pre-defined length of the temporal pattern had been set to 500 ms or 2000 ms (Supplemental Figs. 1a–c). While the IPI distributions within the patterns were similar, the lengths of the high frequency and low frequency components varied with the length of the pattern indicating that the balance of high frequency and low frequency components mattered more than the specific length of either component.
We scored the temporal patterns across the model pain states, in addition to the baseline model state, to identify patterns that performed well across heterogeneous conditions. CF stimulation above 50 Hz completely suppressed model WDR neuron activity across each of the individual model nodes before the implementation of model variants to represent pain states (Fig. 2a1). Therefore, the performance improvements we observed between the temporal patterns and CF stimulation were due to improvements in efficiency that accompanied equivalent efficacy. After reducing GABAergic inhibition, optimized temporal patterns improved efficacy relative to CF stimulation because CF stimulation did not completely inhibit model WDR activity. In computational model states with reduced conductance of Aβ fiber inputs to inhibitory interneurons CF stimulation did completely inhibit model WDR neuron activity, albeit requiring higher frequencies than in the baseline network state (Fig. 2a3). Nevertheless, the non-regular temporal patterns applied in this network state had efficacy scores comparable to CF stimulation and improved efficiency. Overall, the average frequency of the optimized patterns usually converged to between 30 and 40 Hz. The best pattern from the GA optimization (pattern G for “Good”) also reduced the variance of neuron responses across pain states compared to 50 Hz CF stimulation. Across multiple pain states, the optimized stimulation patterns with high fitness had similar temporal properties, including entropy, entropy per pulse, average frequency, and IPI CV (Fig. 2c), indicating that the algorithm converged to specific temporal features that improved performance. Overall, the algorithm identified patterns based on specific temporal features that performed better than the best CF stimulation across a variety of metrics in the computational model.
We designed four irregular temporal patterns to test in vivo. First, we selected the best pattern from the GA optimized across multiple pain states (pattern G). Like other iterations of the GA, this pattern exhibited two primary frequency bands: the high frequency section had IPIs between 8 and 12 ms and lasted for 230 ms while the low frequency section had IPIs between 50 and 70 ms and lasted for 770 ms (Fig. 3b). We also generated patterns with random IPIs drawn from a Poisson distribution with an average frequency of 30 Hz (pattern R1 for “Random 1”) or 50 Hz (pattern R2 for “Random 2”, Fig. 3b). Finally, presuming that the computational model could also predict temporal patterns that would perform worse than CF stimulation, we selected the pattern with the worst fitness with the same cost function with an average frequency less than 50 Hz (pattern B for “Bad”). Patterns R1, R2, and B had a much wider distribution of IPIs than pattern G, and they did not suppress the model WDR neuron as effectively as pattern G suggesting that neither of the random patterns nor pattern B would improve performance in vivo.
3.2 Evaluating the effect of temporal patterns on neural responses in vivo
We quantified the effect of multiple patterns of SCS on single unit activity recorded in the lumbar dorsal horn of urethane-anesthetized rats. We recorded a total of 101 units in four healthy male Sprague-Dawley rats and 117 units in four rats with spared nerve injury (SNI). The SNI rats all developed hypersensitivity in the two weeks following the nerve injury procedure prior to the neural recordings (Fig. 4b). Across the entire population, the changes in firing rate in response to SCS were heterogeneous, but stimulation patterns produced significantly different responses (Kruskal-Wallis test, healthy: χ2 = 30.1, p < 0.001 SNI: χ2 = 12.4, p < 0.05). The heterogeneous responses made it difficult to distinguish the effects of different stimulation patterns, and we focused on the population of neurons most likely to represent the model WDR neurons. Model WDR neurons were spontaneously firing in neuropathic model instantiations and were generally inhibited by stimulation, so we classified neurons based on their spontaneous firing rates and whether they were inhibited or excited by SCS (Fig. 5a1 ). 41% of units recorded in healthy animals and 46% of units recorded in SNI animals exhibited spontaneous firing rates of at least 5 Hz, and we focused on responses of these units to different temporal patterns of SCS (Supplemental Fig. 2b).
We also evaluated the effect of individual stimulation patterns on a subgroup of spontaneously firing neurons (13 in healthy animals and 13 in SNI animals) inhibited by SCS that were classified as WDR neurons in their responses to mechanical stimulation. These units showed a monotonic increase in their response to successive brush, press, pinch, and crush inputs on the ipsilateral hindpaw (Supplemental Fig. 2a) and represent the most similar population to the WDR model neurons. These WDR neurons showed a similar distribution of firing rates overall between healthy and SNI animals, but the shift in neurons with high spontaneous firing rates went from 39% to 51% in SNI animals (Fig. 5a2), indicating that these WDR neurons may have been sensitized following neuropathic pain induction. While mechanical classification is a common way to identify pain responses in rats [
], neurons with biphasic waveforms are more likely to be inhibitory, and neurons with monophasic waveforms are more likely to be excitatory. We found that stimulation was more likely to evoke responses in biphasic neurons in naïve animals (Supplemental Fig. 6b1), and furthermore these increases in activity were more pronounced for 50 Hz CF stimulation, pattern G, and pattern R1.
Responses to temporal patterns of SCS differed between healthy and SNI animals. SCS inhibited a significantly higher proportion of neurons in SNI animals than in healthy animals (56% vs. 28%, χ2 = 39.9, p < 0.001), justifying the use of model pain states in addition to the model default state to design patterns. 50 Hz, pattern G, pattern R1, and pattern R2 all inhibited a higher percentage of neurons in SNI animals than in healthy animals, while pattern B excited a higher percentage in SNI animals (Fig. 5b and c). We did not detect a difference in the percent of neurons inhibited by different patterns of SCS in healthy animals (χ2 = 2.1, p = 0.7), but there was a difference in SNI animals (χ2 = 13.4, p < 0.01) and both pattern G (χ2 = 10.3) and pattern R2 (χ2 = 10.3) inhibited a larger percentage of neurons compared to pattern B (p < 0.05, Bonferroni corrected). These results support our hypothesis that pattern G performed better in vivo than pattern B.
We evaluated the mean change in firing rate of individual neurons that were inhibited or excited by SCS separately in healthy and SNI animals and found that 50 Hz stimulation and pattern G produced the largest changes in firing rates for both neurons inhibited and excited by SCS. In neurons inhibited by SCS, net inhibition was significantly greater during 50 Hz SCS and pattern G than pattern B or pattern R2 (Fig. 5d1, Aligned Rank Transform (ART) ANOVA [
], p < 0.05, post-hoc Tukey's Test, p < 0.05). In neurons excited by SCS there was a significant effect for pain state, stimulation pattern, and the interaction between them (ART ANOVA pain state p < 0.05, stimulation pattern p < 0.001, pain state∗stimulation pattern, p < 0.001). These results support our previous finding that these patterns increase the inhibitory drive in the dorsal horn; however, differences were not as pronounced for SNI animals (Supplemental Fig. 6b2): In healthy animals, 50 Hz CF, pattern G, and pattern R1 all increased “biphasic” neural firing rates more than pattern B, but we did not detect a difference in average increase in firing rates across patterns in SNI animals. These results indicate that 50 Hz CF, pattern G, and pattern R1 may be better at driving putatively inhibitory neural activity despite each pattern having a different average frequency.
In the subgroup of spontaneously firing neurons (13 in healthy animals and 13 in SNI animals) inhibited by SCS that were classified as WDR neurons, changes in firing rate of aligned with the trends we observed in the broader population of neurons, with greater inhibition during 50 Hz and pattern G compared to pattern B and pattern R2 (Fig. 5e, ART ANOVA, p < 0.05). Overall, these results indicate that 50 Hz CF stimulation and pattern G had larger effects on dorsal horn neuron firing rates than other patterns of SCS. Sciatic nerve stimulation at intensities sufficient to activate C-fibers elevated the activity of neurons, but response distribution and firing rate differences between stimulation patterns for neurons inhibited or excited by SCS were not significant (Supplemental Fig. 3). However, sciatic nerve stimulation experiments were intended only to assess effects following generally elevated firing rates and did not necessarily replicate the conditions of the model-based optimization, while experiments with spontaneously firing neurons following neuropathic pain were intended to more closely replicate the conditions tested in the model.
We also quantified changes in the temporal properties of the responses of individual units during stimulation with different temporal patterns. We calculated the firing pattern entropy (Hx), a measure of the randomness within the spike train, and calculated correlations between pairs of neurons recorded simultaneously using the spike time tiling coefficient (STTC), a bounded correlation measure independent of firing rate [
]. Correlations between units may be a biomarker of SCS efficacy, as prior experiments demonstrated an increase in correlations between putatively inhibitory neurons during SCS and reduced correlations between putatively excitatory and inhibitory neurons during effective SCS [
]. We hypothesized that 50 Hz CF stimulation and pattern G would reduce the firing pattern entropy in neural responses since stimulation patterns with higher fitness generally had lower entropy. We first assessed the effect of stimulation patterns on WDR neurons in the computational model and found that 50 Hz CF stimulation and pattern G reduced entropy compared to baseline activity across 10 simulations (Fig. 6a1 ). Entropy at baseline was similar for all model WDR neurons. Interestingly, model WDR neurons showed a decrease in correlations with 50 Hz CF stimulation, but increased correlations during all temporal patterns of stimulation (Fig. 6a2), indicating that CF stimulation may be more effective at reducing correlations between neurons. We also calculated unit entropy (Hxunit), a two-dimensional estimate of entropy based on the interspike interval (ISI) distribution (see Methods), for each neuron. Without stimulation, the ISI distributions of units recorded in healthy animals tended to have higher peaks, indicating regular firing patterns, while ISI distributions of units recorded in SNI animals tended to be broader, indicating less regular firing (Supplemental Fig. 4). The entropy and unit entropy of neurons recorded in SNI animals were higher than in healthy animals (Fig. 6b), indicating that neuropathic pain may increase the randomness of information conveyed through neuron spiking and alter the temporal code in addition to rate codes.
SCS led to significant changes in the activity patterns of recorded units, and both neurons inhibited by SCS and neurons excited by SCS showed more narrow ISI distributions in healthy animals than in SNI animals, indicating more regular firing. 50 Hz SCS reduced entropy in neurons inhibited by SCS compared to no stimulation and relative to pattern B, which increased entropy (ANOVA, p < 0.001, post-hoc Tukey's, p < 0.05, Fig. 6c1). Only pattern B failed to reduce the unit entropy relative to no stimulation (ANOVA, p < 0.001, post-hoc Tukey's, p < 0.05, Fig. 6c2), and all patterns reduced unit correlations compared to no stimulation (ANOVA, p < 0.001, post-hoc Tukey's, p < 0.05, Fig. 6c3). Correlations between spontaneously firing units were not significantly different between healthy and SNI animals without stimulation. While ISI distributions were consistent across stimulation patterns in units excited by SCS, pattern B and pattern R2 also led to a much broader distribution of ISIs than the other stimulation patterns (Supplemental Fig. 4). Stimulation pattern did not have a significant effect on entropy in excited neurons (p = 0.07), but patterns B and pattern R2 reduced the unit entropy in excited neurons (ANOVA, p < 0.001, post-hoc Tukey's, p < 0.05, Fig. 6b2), and all patterns significantly increased correlations (ANOVA, p < 0.001, post-hoc Tukey's, p < 0.05, Fig. 6b3). Overall, these results demonstrate that the stimulation patterns that drove significant inhibition (50 Hz and pattern G) also tended to reduce entropy and regularize the firing patterns of units inhibited by stimulation, i.e., they more consistently modulated the rate coding of DH neurons (both inhibitory and excitatory) compared to pattern B and pattern R2.
Chronic neuropathic pain disrupts inhibitory mechanisms in the dorsal horn and alters the balance of excitation and inhibition. SCS is thought to counteract some of the disinhibition and hyperexcitability by driving activity in inhibitory interneurons to restore healthy processing of somatosensory inputs. We combined model-based optimization and recordings of dorsal horn unit activity to investigate the efficacy of temporal patterns of stimulation for treating chronic pain. We demonstrated that a temporal pattern of SCS designed using computational evolution had comparable effects as 50 Hz stimulation, the clinical standard [
], but with a lower average frequency, thereby improving the energy efficiency of stimulation, and achieving more consistent performance across diverse pain states.
We designed novel temporal patterns of SCS using a computational model and an evolutionary algorithm to identify patterns that were effective across multiple modeled neuropathic pain states. Model-based optimization enabled efficient exploration of the parameter space while optimizing for multiple components of the fitness function. We used a genetic algorithm developed specifically to combine principles of neural stimulation with elements of biological evolution [
]. The genetic algorithm incorporated specific functions to predict characteristics of the stimulation patterns that would be more effective, and this may have led, in part, to convergence to a solution with two primary frequencies. The resulting pattern of SCS had efficacy comparable to 50 Hz CF stimulation while improving energy efficiency by 20%, an important consideration for SCS devices because improved battery life increases the time between surgeries for battery replacement [
] or increases the intervals between recharging. We also used a fitness function that accounted for the resilience of temporal patterns across multiple modeled neuropathic pain states including reduced GABAergic inhibition and reduced Aβ inputs within the dorsal horn network model [
]. Nevertheless, the pattern we identified had comparable or superior efficacy and comparable or superior efficiency compared to the best CF stimulation pattern even when extended across multiple pain states (Fig. 2). However, neuropathic pain-related changes in the dorsal horn network are quite diverse and incorporating patient-specific models in future work may provide better representation of the variance in responses across a population [
The optimal pattern designed with the genetic algorithm had effects on dorsal horn neurons recorded in vivo similar to those of 50 Hz stimulation. We focused on spontaneously active neurons because we evaluated the effect of temporal patterns on spontaneously active WDR neurons in the model. However, this population likely included both excitatory and inhibitory interneurons in addition to projection neurons. Prior recordings of projection neuron responses indicate that there is substantial heterogeneity in responses to CF stimulation, even within nociceptive-specific and WDR projection neurons [
], and thus we focused our analysis on the magnitude of inhibition and excitation during different stimulation patterns. 50 Hz stimulation and pattern G, the optimal GA pattern, had the strongest effects on spontaneously active dorsal horn neurons, including both greater inhibition and greater excitation, than other patterns. 50 Hz stimulation and pattern G also led to greater net inhibition of spontaneously firing mechanically classified WDR neurons than other patterns (Fig. 5e). Further, pattern B, the worst pattern from the GA, had the weakest net effects on dorsal horn neurons. While we cannot evaluate effects on neuropathic pain without behavioral tests, strong activation of the dorsal horn network indicates an increased drive through dorsal column axons to the dorsal horn network. Random stimulation at 30 Hz had larger net effects than random stimulation at 50 Hz, indicating that the higher frequency dorsal column stimulation does not necessarily translate to stronger network effects.
Both neuropathic pain and individual stimulation patterns altered the patterns of activity of individual units. Furthermore, trends in the changes to firing patterns in response to each stimulation pattern were similar between model WDR neurons and neurons recorded in vivo that were inhibited by stimulation. Neurons recorded in nerve-injured animals showed greater randomness in their firing patterns, or increased unit entropy (Fig. 6, Supplemental Fig. 4). However, this could be due to differences among the units we happened to sample with our recording electrode. Nevertheless, increased entropy in neuropathic pain is an interesting finding because it parallels observations from other neurological disorders. For example, increased entropy in the GP and SNr is also associated with transitions to parkinsonian states [
]. Similarly, we hypothesize that increased entropy in chronic pain decreases information transmission from dorsal column inputs to projection neurons, but additional experiments are necessary to investigate this hypothesis. Although our metrics of entropy are independent of firing rate, the changes in entropy paralleled changes in firing rates in spontaneously firing neurons inhibited by SCS. As well, stimulation reduced correlations between neurons inhibited by SCS and increased correlations between neurons excited by SCS. The correlation metric was independent of firing rate, and the decrease of correlation between units inhibited by SCS may indicate a restoration towards the baseline state through decreased transmission between excitatory neurons. However, decreased entropy leading to decreased correlations between units is counterintuitive because more regular firing could often lead to increased correlations. In this case, increased correlations may correspond to mixed somatosensory processing, or spreading excitatory inputs to other receptive field areas or sensory processing circuits. Decreasing correlations may limit the spread of excitation across dorsal horn circuits responsible for processing different sensory modalities [
The results represent an initial test of optimized temporal patterns of SCS with several limitations. First, the fitness function used for optimization combined four elements of performance, and a fitness function that increased the relative importance of efficacy may have resulted in a pattern that produced greater neuronal inhibition than 50 Hz. We used a population of models to represent multiple pain states, but these pain states were limited to changes in GABAergic inhibition and reduced Aβ inputs, and future studies may consider other forms of disinhibition as well [
]. Additionally, SCS-mediated inhibition of WDR neurons in the model was dependent on the number of fibers activated, i.e., stimulation amplitude, but this was outside of the parameters we considered for this study. Also, we did not consider several important dimensions within the stimulation parameter space to focus on temporal pattern optimization. We assumed that each pulse in the temporal patterns would activate a fixed number of dorsal column inputs in the computational model, but this ignores important interactions between stimulation frequency, pulse width, and amplitude with axonal activation. A more realistic depiction of the effects of stimulation patterns on dorsal column axons would strengthen our understanding of the effects of individual temporal patterns, albeit with a substantially higher computational cost. We also applied each pattern at conventional paresthesia-based stimulation amplitudes (60–70% of motor threshold), but recent clinical studies indicated that stimulation at subthreshold amplitudes with conventional frequencies may be more effective than stimulation at higher amplitudes that produce paresthesias [
]. We did not incorporate supraspinal effects or non-neuronal effects into the model as these effects are difficult to quantify and incorporate into a computational model. Finally, while efficiency is an important consideration for battery life, we did not consider programming time with temporal patterns or interactions between temporal patterns and active contacts that may impact programming.
Within our in vivo studies, we focused our analyses on individual dorsal horn neurons, but future experiments should evaluate behavioral effects of stimulation in awake animals as a next step to potential clinical use. We also did not classify individual units as projection neurons and focused on classifying individual neurons based on their spontaneous firing rates and responses to mechanical inputs. We did attempt to classify neurons based on their waveform shape [
], and found significant differences between these classification groups (Supplemental Fig. 6). However, differences were not as clear, in part because we were still averaging responses across a diverse population to find the net response. Classifying neurons based on their responses to stimulation may oversimplify some aspects of the analysis, but it also reveals important functional differences between the stimulation patterns. Differences between stimulation patterns were also less pronounced in the SNI animals than naïve animals, perhaps indicating that the model did not accurately capture all the relevant changes that occur in neuropathic pain. Conversely, the search strategy of minimizing response variation (i.e. maximizing pattern resilience) across model states may have precluded the discovery of a pattern optimized to the specific dorsal horn changes that occur following SNI (e.g. if such a pattern proved suboptimal for a different modeled state). However, the differences observed in naïve animals are encouraging and provide some validation for the model, and future work may involve the customized design of patterns tailored to specific pain states. While previous studies have characterized temporal properties of afferent inputs, the temporal properties of dorsal horn neuron activity and the relationship to pain behaviors remains unknown. Our results showing that SCS led to significant changes in the firing patterns of recorded units highlight the importance of considering the temporal patterns of neural activity, as opposed to only quantifying the effects on average firing rates. We also did not perform assessments on awake, behaving animals in this study, although the firing rate of dorsal horn neurons is strongly correlated with pain sensation [
]. However, previous experiments that evaluated random patterns of stimulation did not find significant differences in behavioral effects between random patterns and CF stimulation until 45–60 min after stimulation onset, but we only evaluated immediate effects of stimulation that occurred within the first 30 s [
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
John E Gilbert is an inventor on intellectual property on temporal patterns of spinal cord stimulation, owned by Duke University, for which he receives royalty payments.
Tianhe C. Zhang is an inventor on intellectual property on temporal patterns of spinal cord stimulation, owned by Duke University, for which he receives royalty payments. Tianhe Zhang is a full-time employee of Boston Scientific Neuromodulation, and Tianhe Zhang owns shares of Boston Scientific stock.
Rosana Esteller is an employee of Boston Scientific and holds stocks with this company and also has stocks with NeuroPace, Inc.
Warren M Grill is an inventor on intellectual property on temporal patterns of spinal cord stimulation, owned by Duke University, for which he receives royalty payments. Warren M Grill is Principal Investigator on Research Grants to Duke University sponsored by Boston Scientific . Warren M Grill serves as a Paid Consultant for Boston Scientific.
We acknowledge the Duke Compute Cluster for their assistance with our simulations. We thank Danielle Degoski for her assistance with experiments. This work was supported by a grant to Duke University from Boston Scientific Corporation . This work was also supported by the Robert Plonsey Fellowship from the Duke University Department of Biomedical Engineering . T.Z. and R.E. are employees of Boston Scientific Corporation . J.E.G, N.D.T, T.Z. and W.M.G. have received royalty payments from Boston Scientific Corporation for licensed I.P. W.M.G received compensation from Boston Scientific as a member of the Neuromodulation Scientific Advisory Board. T.Z. and R.E. own Boston Scientific stock and R.E owns stock in NeuroPace.
Appendix A. Supplementary data
The following is the Supplementary data to this article: