Advertisement

Evaluating optimized temporal patterns of spinal cord stimulation (SCS)

  • John E. Gilbert
    Affiliations
    Department of Biomedical Engineering, Duke University, Durham, NC, USA
    Search for articles by this author
  • Tianhe Zhang
    Affiliations
    Neuromodulation Research and Advanced Concepts, Boston Scientific Neuromodulation, Valencia, CA, USA
    Search for articles by this author
  • Rosana Esteller
    Affiliations
    Neuromodulation Research and Advanced Concepts, Boston Scientific Neuromodulation, Valencia, CA, USA
    Search for articles by this author
  • Warren M. Grill
    Correspondence
    Corresponding author. Duke University, Department of Biomedical Engineering, Fitzpatrick CIEMAS, Room 1427, 101 Science Drive, Box 90281, Durham, NC, 27708-0281, USA.
    Affiliations
    Department of Biomedical Engineering, Duke University, Durham, NC, USA

    Department of Electrical and Computer Engineering, Duke University, Durham, NC, USA

    Department of Neurobiology, Duke University School of Medicine, Durham, NC, USA

    Department of Neurosurgery, Duke University School of Medicine, Durham, NC, USA
    Search for articles by this author
Open AccessPublished:July 30, 2022DOI:https://doi.org/10.1016/j.brs.2022.07.046

      Highlights

      • 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.

      Abstract

      Background

      Temporal patterns of stimulation represent a novel dimension for improving the efficacy of spinal cord stimulation to treat chronic neuropathic pain.

      Objective

      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.

      Methods

      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.

      Results

      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.

      Conclusions

      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.

      Graphical abstract

      Keywords

      Abbreviations:

      Dorsal horn (DH), Dorsal columns (DC), Spinal cord stimulation (SCS), Motor threshold (MT), Perception threshold (PT), Activation threshold (AT), Inter-spike interval (ISI)

      1. Introduction

      Chronic neuropathic pain is a frequent and compelling reason for seeking medical attention and significantly decreases quality of life [
      • Elliott A.M.
      • Smith B.H.
      • Penny K.I.
      • Cairns Smith W.
      • Alastair Chambers W.
      The epidemiology of chronic pain in the community.
      ,
      • Meyr A.J.
      • Saffran B.
      The pathophysiology of the chronic pain cycle.
      ,
      • Nahin R.L.
      Estimates of pain prevalence and severity in adults: United States, 2012.
      ]. 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 [
      • Kemler M.A.
      • Raphael J.H.
      • Bentley A.
      • Taylor R.S.
      The cost-effectiveness of spinal cord stimulation for complex regional pain syndrome.
      ,
      • Manca A.
      • Kumar K.
      • Taylor R.S.
      • Jacques L.
      • Eldabe S.
      • Meglio M.
      • Molet J.
      • Thomson S.
      • O'Callaghan J.
      • Eisenberg E.
      • Milbouw G.
      • Buchser E.
      • Fortini G.
      • Richardson J.
      • Taylor R.J.
      • Goeree R.
      • Sculpher M.J.
      Quality of life, resource consumption and costs of spinal cord stimulation versus conventional medical management in neuropathic pain patients with failed back surgery syndrome (PROCESS trial).
      ,
      • Rosenblatt R.A.
      • Catlin M.
      Opioids for chronic pain: first do No harm.
      ,
      • Rosenblum A.
      • Marsch L.A.
      • Joseph H.
      • Portenoy R.K.
      Opioids and the treatment of chronic pain: controversies, current status, and future directions.
      ]. SCS efficacy is dependent on stimulation parameters, including the pulse amplitude, pulse duration, and pulse repetition rate (frequency) [
      • Yearwood T.L.
      • Hershey B.
      • Bradley K.
      • Lee D.C.
      Pulse width programming in spinal cord stimulation: a clinical study.
      ]. 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 [
      • Miller J.P.
      • Eldabe S.
      • Buchser E.
      • Johanek L.M.
      • Guan Y.
      • Linderoth B.
      Parameters of spinal cord stimulation and their role in electrical charge delivery: a review.
      ].
      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 [
      • Grill W.M.
      Temporal pattern of electrical stimulation is a new dimension of therapeutic innovation.
      ]. 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 [
      • Birznieks I.
      • Vickery R.M.
      Spike timing matters in novel neuronal code involved in vibrotactile frequency perception.
      ,
      • Johansson R.S.
      • Birznieks I.
      First spikes in ensembles of human tactile afferents code complex spatial fingertip events.
      ,
      • Saal H.P.
      • Wang X.
      • Bensmaia S.J.
      Importance of spike timing in touch: an analogy with hearing?.
      ], while the firing rate of neurons that project to higher order brain centers in primates is correlated with the amount of pain reported by human subjects [
      • Simone D.A.
      • Sorkin L.S.
      • Oh U.
      • Chung J.M.
      • Owens C.
      • LaMotte R.H.
      • Willis W.D.
      Neurogenic hyperalgesia: central neural correlates in responses of spinothalamic tract neurons.
      ]. According to the gate control theory, SCS drives activity in large diameter afferents that in turn drive inhibition in dorsal horn networks [
      • Melzack R.
      • Wall P.D.
      Pain mechanisms: a new theory.
      ]; 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 [
      • Edhi M.M.
      • Heijmans L.
      • Vanent K.N.
      • Bloye K.
      • Baanante A.
      • Jeong K.S.
      • Leung J.
      • Zhu C.
      • Esteller R.
      • Saab C.Y.
      Time-dynamic pulse modulation of spinal cord stimulation reduces mechanical hypersensitivity and spontaneous pain in rats.
      ]. 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 [
      • Al-Kaisy A.
      • Van Buyten J.-P.
      • Smet I.
      • Palmisani S.
      • Pang D.
      • Smith T.
      Sustained effectiveness of 10 kHz high-frequency spinal cord stimulation for patients with chronic, low back pain: 24-month results of a prospective multicenter study.
      ], but they are energy-inefficient when delivered continuously, resulting in rapid battery depletion or frequent recharging requirements [
      • Paz-Solis J.
      • Thomson S.
      • Jain R.
      • Chen L.
      • Huertas I.
      • Doan Q.
      Exploration of high and low frequency options for subperception spinal cord stimulation using neural dosing parameter relationships: the HALO study.
      ,
      • Vesper J.
      • Slotty P.
      • Schu S.
      • Poeggel-Kraemer K.
      • Littges H.
      • Van Looy P.
      • Agnesi F.
      • Venkatesan L.
      • Van Havenbergh T.
      Burst SCS microdosing is as efficacious as standard burst SCS in treating chronic back and leg pain: results from a randomized controlled trial.
      ] and potentially therapeutically suboptimal as well [
      • Hagedorn J.M.
      • Engle A.M.
      • Ghosh P.
      T RD. Device profile of the Proclaim XR neurostimulation system for the treatment of chronic pain: an overview of its safety and efficacy.
      ], 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 [
      • Leitner J.
      • Westerholz S.
      • Heinke B.
      • Forsthuber L.
      • Wunderbaldinger G.
      • Jager T.
      • Gruber-Schoffnegger D.
      • Braun K.
      • Sandkuhler J.
      Impaired excitatory drive to spinal GABAergic neurons of neuropathic mice.
      ], 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 [
      • Gilbert J.E.
      • Zhang T.
      • Esteller R.
      • Grill W.M.
      Distributed network model of effects of spinal cord stimulation (SCS) on wide-dynamic range (WDR) neurons.
      ]. 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 [
      • Simone D.A.
      • Sorkin L.S.
      • Oh U.
      • Chung J.M.
      • Owens C.
      • LaMotte R.H.
      • Willis W.D.
      Neurogenic hyperalgesia: central neural correlates in responses of spinothalamic tract neurons.
      ]. 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) [
      • Kajander K.C.
      • Bennett G.J.
      Onset of a painful peripheral neuropathy in rat: a partial and differential deafferentation and spontaneous discharge in Aβ and Aδ primary afferent neurons.
      ,
      • Liu X.
      • Eschenfelder S.
      • Blenk K.-H.
      • Jänig W.
      • Häbler H.-J.
      Spontaneous activity of axotomized afferent neurons after L5 spinal nerve injury in rats.
      ,
      • Wall P.D.
      • Gutnick M.
      Ongoing activity in peripheral nerves: the physiology and pharmacology of impulses originating from a neuroma.
      ]. 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 [
      • Zhang T.C.
      • Janik J.J.
      • Grill W.M.
      Modeling effects of spinal cord stimulation on wide-dynamic range dorsal horn neurons: influence of stimulation frequency and GABAergic inhibition.
      ]. Synaptic connections between nodes were based on the spatial distribution of inputs described by Hillman and Wall [
      • Hillman P.
      • Wall P.D.
      Inhibitory and excitatory factors influencing the receptive fields of lamina 5 spinal cord cells.
      ]; 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 [
      • Kumar K.
      • Taylor R.S.
      • Jacques L.
      • Eldabe S.
      • Meglio M.
      • Molet J.
      • Thomson S.
      • O'Callaghan J.
      • Eisenberg E.
      • Milbouw G.
      • Buchser E.
      • Fortini G.
      • Richardson J.
      • North R.B.
      Spinal cord stimulation versus conventional medical management for neuropathic pain: a multicentre randomised controlled trial in patients with failed back surgery syndrome.
      ,
      • Sandkuhler J.
      The role of inhibition in the generation and amplification of pain.
      ]. For example, neuropathic pain alters the balance of GABAergic inhibition and reduces the drive onto inhibitory interneurons from Aβ fibers [
      • Leitner J.
      • Westerholz S.
      • Heinke B.
      • Forsthuber L.
      • Wunderbaldinger G.
      • Jager T.
      • Gruber-Schoffnegger D.
      • Braun K.
      • Sandkuhler J.
      Impaired excitatory drive to spinal GABAergic neurons of neuropathic mice.
      ,
      • Moore K.A.
      • Kohno T.
      • Karchewski L.A.
      • Scholz J.
      • Baba H.
      • Woolf C.J.
      Partial peripheral nerve injury promotes a selective loss of GABAergic inhibition in the superficial dorsal horn of the spinal cord.
      ]. 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 [
      • Zhang T.C.
      • Janik J.J.
      • Grill W.M.
      Modeling effects of spinal cord stimulation on wide-dynamic range dorsal horn neurons: influence of stimulation frequency and GABAergic inhibition.
      ]. All network simulations were completed using Neuron v7.7 [
      • Hines M.
      • Carnevale N.
      The NEURON simulation environment.
      ].
      We used a genetic algorithm (GA) designed specifically for optimizing temporal patterns of neural stimulation [
      • Cassar I.R.
      • Titus N.D.
      • Grill W.M.
      An improved genetic algorithm for designing optimal temporal patterns of neural stimulation.
      ]. 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. [
      • Cassar I.R.
      • Titus N.D.
      • Grill W.M.
      An improved genetic algorithm for designing optimal temporal patterns of neural stimulation.
      ] 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 [
      • North R.B.
      • Ewend M.G.
      • Lawton M.T.
      • Piantadosi S.
      Spinal cord stimulation for chronic, intractable pain: superiority of "multi-channel" devices.
      ]. 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 [
      • Simone D.A.
      • Sorkin L.S.
      • Oh U.
      • Chung J.M.
      • Owens C.
      • LaMotte R.H.
      • Willis W.D.
      Neurogenic hyperalgesia: central neural correlates in responses of spinothalamic tract neurons.
      ]. 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 [
      • D'Mello R.
      • Dickenson A.H.
      Spinal cord mechanisms of pain.
      ] 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 [
      • Dorval A.D.
      Estimating neuronal information: logarithmic binning of neuronal inter-spike intervals.
      ]. For these simulations, we used a weight matrix of [2,2,2,1,1] for W1⋯W5. GA code was implemented using MATLAB (R2020a, Mathworks, Natick MA).
      Fitnessnet=W1zscore(Effect)+W2zscore(Control)+W3zscore(Surround)+W4zscore(Efficiency)+W5zscore(Variance)
      (1)


      Effect=FR(WDRz1)
      (2)


      Control=FR(WDRz1)FR(WDRz1CF)
      (3)


      Surround=0.66WDRZ2+0.33WDRZ3
      (4)


      Efficiency=#pulses
      (5)


      Effect=FR(WDRz1)
      (6)


      2.2 Neural recordings

      The protocol for neural recordings has been described elsewhere [
      • Gilbert J.E.
      • Zhang T.C.
      • Esteller R.
      • Grill W.M.
      Pharmacological disinhibition degrades neural population responses to Aβ-fiber electrical stimulation (Aβ-ES).
      ,
      • Zhang T.C.
      • Janik J.J.
      • Peters R.V.
      • Chen G.
      • Ji R.R.
      • Grill W.M.
      Spinal sensory projection neuron responses to spinal cord stimulation are mediated by circuits beyond gate control.
      ]. 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 [
      • Decosterd I.
      • Woolf C.J.
      Spared nerve injury: an animal model of persistent peripheral neuropathic pain.
      ]. 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 [
      • Bonin R.P.
      • Bories C.
      • De Koninck Y.
      A simplified up-down method (SUDO) for measuring mechanical nociception in rodents using von Frey filaments.
      ]. 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 [
      • Shimazaki H.
      • Shinomoto S.
      A method for selecting the bin size of a time histogram.
      ]. 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 [
      • Dorval A.D.
      • Grill W.M.
      Deep brain stimulation of the subthalamic nucleus reestablishes neuronal information transmission in the 6-OHDA rat model of parkinsonism.
      ]. Interspike interval distributions (ISIs) were binned logarithmically with 10 bins per decade [
      • Dorval A.D.
      Estimating neuronal information: logarithmic binning of neuronal inter-spike intervals.
      ,
      • Dorval A.D.
      • Grill W.M.
      Deep brain stimulation of the subthalamic nucleus reestablishes neuronal information transmission in the 6-OHDA rat model of parkinsonism.
      ]. 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 [
      • Cutts C.S.
      • Eglen S.J.
      Detecting pairwise correlations in spike trains: an objective comparison of methods and application to the study of retinal waves.
      ].
      Hx=aP(ISIa)log2P(ISIa)
      (7)


      Hunit=baP(ISIa,ISIb)log2P(ISIa|ISIb)
      (8)


      3. Results

      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) [
      • Cassar I.R.
      • Titus N.D.
      • Grill W.M.
      An improved genetic algorithm for designing optimal temporal patterns of neural stimulation.
      ]. 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.
      Fig. 1
      Fig. 1Designing temporal patterns of spinal cord stimulation using a computational model. a, Example inputs and outputs for a single node within the dorsal horn network. The model contains primary afferents (Aβ, Aδ, and C) that transmit inputs to local inhibitory (IN) and excitatory (EX) interneurons and a wide-dynamic range (WDR) projection neuron. Constant frequency (CF) or temporal patterns of SCS are represented as inputs to the dorsal columns along the Aβ fibers. The voltage traces of the output neurons without any stimulation inputs are shown on the right. Scale bar, 100 ms and 30 mV. b, Distributed multinodal model architecture. Each circle represents a node of the model with unique inputs. The center node (Zone 1) receives excitatory inputs from zone 2 nodes and inhibitory inputs from both zone 2 and zone 3 nodes. All connections between nodes are from interneurons to WDR neurons. Each zone receives pattern inputs from dorsal columns/Aβ fibers. c, Temporal patterns of stimulation are the inputs to the model. The algorithm evaluates the fitness of the pattern based on how much they reduce pain in the model, the efficiency of the pattern, the reduction in variance of the pain reduction across multiple trials, and the efficacy of the pattern compared to CF controls. The genetic algorithm operates similar to biological evolution and designs new patterns of stimulation based on the patterns with the highest fitness quantified via a cost function (Equations ((1), (6))) from the previous round and repeats the process for multiple generations.
      Fig. 2
      Fig. 2Example genetic algorithm evolution across pain states. a1, Default CF response in the model and GA optimization with temporal patterns. Each dot represents an individual pattern and patterns for all 50 generations are shown here. Dot color shows the fitness of each individual pattern. Patterns are scored individually within pain states, normalized, and then combined by rank within generations to evaluate the best performing patterns across pains states. Black line represents the default model response to CF stimulation from 1 to 150 Hz. In the default network state, patterns did not improve efficacy compared to CF stimulation but did improve efficiency. a2, Modeled neuropathic pain state with a 50% reduction in the strength of Aβ inputs to IN model neurons. a3, Modeled neuropathic pain state with a 50% reduction in GABAergic inhibition. b, Evolution of the fitness function across 50 generations. The fitness rapidly improves across the first 5 generations, and then slowly improves thereafter. Here fitness was calculated using the z-score for all patterns (i.e., prior to ranking patterns) to show the improvement that occurs across generations. Inputs to the network model are randomly seeded for each generation, which leads to variance in the scores across generations for the same temporal pattern. c, Properties of optimized patterns plotted across generations for patterns from the GA run in (a). Entropy (c1), entropy per pulse (c2), mean frequency (c3), and CV of IPIs (c4) for all patterns. GA converges to an entropy of approximately two bits and a frequency of between 40 and 60 Hz within 50 generations. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
      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.
      Fig. 3
      Fig. 3Temporal patterns of spinal cord stimulation. a, Evolution of pattern with the highest fitness every five generations for the first 50 generations. The final pattern represents the output of the optimization algorithm. Patterns are 1000 ms long and are repeated twice in the figure. b, Temporal patterns tested in vivo: 50 Hz CF stimulation, the pattern with the highest fitness from the GA (G), the pattern with the worst fitness with an average frequency below 50 Hz (B), and patterns with randomly sampled spike times with a mean pulse rate of 30 Hz (R1) and mean pulse rate of 50 Hz (R2). c, Frequency prevalence within each type of pattern. We binned interpulse intervals and calculated the frequency prevalence as the number of interpulse intervals in each bin divided by the average frequency of the bin to balance low and high frequencies within each GA pattern.

      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).
      Fig. 4
      Fig. 4Measuring the effects of temporal patterns of spinal cord stimulation on single unit activity in vivo in urethane-anesthetized rats a, Depiction of experimental setup. For nerve injured animals, we exposed the sciatic nerve and ligated the tibial and peroneal branch of the nerve. In the acute recording setup in both healthy and nerve-injured animals we cuffed the sciatic nerve to stimulate at C-fiber amplitudes. We placed a 4-contact paddle electrode over the dorsal columns between T10 and T11 and inserted a 16- or 32-contact multichannel recording electrode into the dorsal horn between T13 and L1.
      Fig. 5
      Fig. 5Effect of temporal patterns of stimulation on neurons recorded in vivo. a1, Binned spontaneous firing rates of individual units recorded in vivo in healthy and SNI animals. Units with greater than 5 Hz spontaneous firing rate were included in subsequent analysis. a2, Same as a1, but for neurons classified as wide-dynamic range neurons based on their responses to brush, press, pinch, crush stimulation on the ipsilateral hindpaw (see ). The firing rate of neurons recorded in SNI animals was significantly larger than the firing rate of neurons recorded in healthy animals for all neurons and for WDR neurons (Two-sample Kolmogorov-Smirnov test, p < 0.001). b, Percent of neurons that were inhibited (blue), excited (red), or neither (grey). c, Raw (c1) and normalized (c2) changes in firing rate for each stimulation pattern for neurons with a spontaneous firing rate greater than 5 Hz. Within each stimulation pattern, the left bar displays neurons recorded in healthy animals and the right bar shows neurons recorded in SNI animals. Neurons are sorted within each group. d1, Mean change in firing rate for neurons that were inhibited by stimulation. Error bars represent standard error. Lines connecting groups show a significant difference in the mean normalized firing rate. d2, Same as (d1), but for neurons excited by SCS. e, Same as (d1), but for neurons classified as wide-dynamic range neurons that were inhibited by SCS. Black bars represent significant differences across groups (stimulation pattern or pain model). Grey bars represent significant differences within groups where there was a significant interaction between stimulation pattern and pain model (Aligned Rank Transform ANOVA [58], p < 0.05, post-hoc Tukey's test). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
      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 [
      • Deuis J.R.
      • Dvorakova L.S.
      • Vetter I.
      Methods used to evaluate pain behaviors in rodents.
      ], other methods have been proposed to identify functional differences between populations of dorsal horn neurons based on their waveform shape [
      • Lee K.Y.
      • Ratté S.
      • Prescott S.A.
      Excitatory neurons are more disinhibited than inhibitory neurons by chloride dysregulation in the spinal dorsal horn.
      ,
      • Snyder A.C.
      • Morais M.J.
      • Smith M.A.
      Dynamics of excitatory and inhibitory networks are differentially altered by selective attention.
      ]. We used similar methods to classify dorsal horn neurons based on their average waveform shape: based on prior experiments [
      • Lee K.Y.
      • Ratté S.
      • Prescott S.A.
      Excitatory neurons are more disinhibited than inhibitory neurons by chloride dysregulation in the spinal dorsal horn.
      ], 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 [
      • Wobbrock J.O.
      • Findlater L.
      • Gergle D.
      • Higgins J.J.
      The aligned rank transform for nonparametric factorial analyses using only anova procedures.
      ], 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 [
      • Cutts C.S.
      • Eglen S.J.
      Detecting pairwise correlations in spike trains: an objective comparison of methods and application to the study of retinal waves.
      ]. 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 [
      • Gilbert J.E.
      • Zhang T.C.
      • Esteller R.
      • Grill W.M.
      Pharmacological disinhibition degrades neural population responses to Aβ-fiber electrical stimulation (Aβ-ES).
      ]. 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.
      Fig. 6
      Fig. 6Effects of temporal pattern of SCS on patterns of neural activity. a1, Change in average entropy of model WDR neurons across 10 simulations for each stimulation pattern. a2, Change in correlations between model WDR neurons across 10 simulations for each stimulation pattern. b1, Change in the average entropy of individual units excited by stimulation for each pattern. Error bars show standard error. b2, Change in unit entropy of individual units. b3, Change in correlations between units excited by stimulation. c, Same as (b), but for units inhibited by SCS. Bars represent significant differences between stimulation patterns. Bars in legend represent significant difference between healthy and SNI animals (ANOVA p < 0.05, post-hoc Tukey's Test p < 0.05).
      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.

      4. Discussion

      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 [
      • North R.B.
      • Kidd D.H.
      • Zahurak M.
      • James C.S.
      • Long D.M.
      Spinal cord stimulation for chronic, intractable pain: experience over two decades.
      ], 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 [
      • Cassar I.R.
      • Titus N.D.
      • Grill W.M.
      An improved genetic algorithm for designing optimal temporal patterns of neural stimulation.
      ]. 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 [
      • Howell B.
      • Lad S.P.
      • Grill W.M.
      Evaluation of intradural stimulation efficiency and selectivity in a computational model of spinal cord stimulation.
      ] 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 [
      • Rigoard P.
      • Basu S.
      • Desai M.
      • Taylor R.
      • Annemans L.
      • Tan Y.
      • Johnson M.J.
      • Van den Abeele C.
      • North R.
      • Group P.S.
      Multicolumn spinal cord stimulation for predominant back pain in failed back surgery syndrome patients: a multicenter randomized controlled trial. pain.
      ]. The neuropathic pain states shifted the response of model WDR neurons to CF stimulation (Fig. 2), demonstrating that optimal stimulation parameters may shift across pain conditions [
      • Zhang T.C.
      • Janik J.J.
      • Grill W.M.
      Modeling effects of spinal cord stimulation on wide-dynamic range dorsal horn neurons: influence of stimulation frequency and GABAergic inhibition.
      ]. 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 [
      • Lempka S.F.
      • Zander H.J.
      • Anaya C.J.
      • Wyant A.
      • Ozinga JGt
      • Machado A.G.
      Patient-specific analysis of neural activation during spinal cord stimulation for pain.
      ].
      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 [
      • Zhang T.C.
      • Janik J.J.
      • Peters R.V.
      • Chen G.
      • Ji R.R.
      • Grill W.M.
      Spinal sensory projection neuron responses to spinal cord stimulation are mediated by circuits beyond gate control.
      ], 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 [
      • Dorval A.D.
      • Grill W.M.
      Deep brain stimulation of the subthalamic nucleus reestablishes neuronal information transmission in the 6-OHDA rat model of parkinsonism.
      ], and corresponding decrease in information transmission to the thalamus [
      • Guo Y.
      • Rubin J.E.
      • McIntyre C.C.
      • Vitek J.L.
      • Terman D.
      Thalamocortical relay fidelity varies across subtahlamic nucleus deep brain stimulation protocols in a data-driven computational model.
      ]. 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 [
      • Sandkuhler J.
      The role of inhibition in the generation and amplification of pain.
      ,
      • Prescott S.A.
      • Ratte S.
      Pain processing by spinal microcircuits: afferent combinatorics.
      ,
      • Seal R.P.
      Illuminating the gap: neuronal cross-talk within sensory ganglia and persistent pain.
      ]. Regardless, the strong relationship observed between changes in firing rates, entropy, and correlations from these studies support their use as metrics for evaluating SCS parameters [
      • Gilbert J.E.
      • Zhang T.C.
      • Esteller R.
      • Grill W.M.
      Pharmacological disinhibition degrades neural population responses to Aβ-fiber electrical stimulation (Aβ-ES).
      ].
      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 [
      • Coull J.A.M.
      • Boudreau D.
      • Bachand K.
      • Prescott S.A.
      • Nault F.
      • Sík A.
      • De Koninck P.
      • De Koninck Y.
      Trans-synaptic shift in anion gradient in spinal lamina I neurons as a mechanism of neuropathic pain.
      ,
      • Peirs C.
      • Seal R.P.
      Neural circuits for pain: recent advances and current views.
      ]. 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 [
      • Metzger C.S.
      • Hammond M.B.
      • Paz-Solis J.F.
      • Newton W.J.
      • Thomson S.J.
      • Pei Y.
      • Jain R.
      • Moffitt M.
      • Annecchino L.
      • Doan Q.
      A novel fast-acting sub-perception spinal cord stimulation therapy enables rapid onset of analgesia in patients with chronic pain.
      ]. 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 [
      • Lee K.Y.
      • Ratté S.
      • Prescott S.A.
      Excitatory neurons are more disinhibited than inhibitory neurons by chloride dysregulation in the spinal dorsal horn.
      ], 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 [
      • Simone D.A.
      • Sorkin L.S.
      • Oh U.
      • Chung J.M.
      • Owens C.
      • LaMotte R.H.
      • Willis W.D.
      Neurogenic hyperalgesia: central neural correlates in responses of spinothalamic tract neurons.
      ]. 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 [
      • Edhi M.M.
      • Heijmans L.
      • Vanent K.N.
      • Bloye K.
      • Baanante A.
      • Jeong K.S.
      • Leung J.
      • Zhu C.
      • Esteller R.
      • Saab C.Y.
      Time-dynamic pulse modulation of spinal cord stimulation reduces mechanical hypersensitivity and spontaneous pain in rats.
      ].

      CRediT authorship contribution statement

      John E. Gilbert: Investigation, Formal analysis, Methodology, Software, Writing – original draft. Tianhe Zhang: Conceptualization, Project administration, Writing – review & editing. Rosana Esteller: Conceptualization, Project administration, Writing – review & editing. Warren M. Grill: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – review & editing.

      Declaration of competing interest

      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.

      Acknowledgements

      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:

      References

        • Elliott A.M.
        • Smith B.H.
        • Penny K.I.
        • Cairns Smith W.
        • Alastair Chambers W.
        The epidemiology of chronic pain in the community.
        Lancet. 1999; 354: 1248-1252
        • Meyr A.J.
        • Saffran B.
        The pathophysiology of the chronic pain cycle.
        Clin Podiatr Med Surg. 2008; 25: 327-346
        • Nahin R.L.
        Estimates of pain prevalence and severity in adults: United States, 2012.
        J Pain. 2015; 16: 769-780
        • Kemler M.A.
        • Raphael J.H.
        • Bentley A.
        • Taylor R.S.
        The cost-effectiveness of spinal cord stimulation for complex regional pain syndrome.
        Value Health. 2010; 13: 735-742
        • Manca A.
        • Kumar K.
        • Taylor R.S.
        • Jacques L.
        • Eldabe S.
        • Meglio M.
        • Molet J.
        • Thomson S.
        • O'Callaghan J.
        • Eisenberg E.
        • Milbouw G.
        • Buchser E.
        • Fortini G.
        • Richardson J.
        • Taylor R.J.
        • Goeree R.
        • Sculpher M.J.
        Quality of life, resource consumption and costs of spinal cord stimulation versus conventional medical management in neuropathic pain patients with failed back surgery syndrome (PROCESS trial).
        Eur J Pain. 2008; 12: 1047-1158
        • Rosenblatt R.A.
        • Catlin M.
        Opioids for chronic pain: first do No harm.
        Ann Fam Med. 2012; 10: 300-301
        • Rosenblum A.
        • Marsch L.A.
        • Joseph H.
        • Portenoy R.K.
        Opioids and the treatment of chronic pain: controversies, current status, and future directions.
        Exper Clin Pychopharmacol. 2008; 16: 405-416
        • Yearwood T.L.
        • Hershey B.
        • Bradley K.
        • Lee D.C.
        Pulse width programming in spinal cord stimulation: a clinical study.
        Pain Physician. 2010; 13: 321-335
        • Miller J.P.
        • Eldabe S.
        • Buchser E.
        • Johanek L.M.
        • Guan Y.
        • Linderoth B.
        Parameters of spinal cord stimulation and their role in electrical charge delivery: a review.
        Neuromodulation. 2016; 19: 373-384
        • Grill W.M.
        Temporal pattern of electrical stimulation is a new dimension of therapeutic innovation.
        Curr Opin Biomed Eng. 2018; 8: 1-6
        • Birznieks I.
        • Vickery R.M.
        Spike timing matters in novel neuronal code involved in vibrotactile frequency perception.
        Curr Biol. 2017; 27 (e1482): 1485-1490
        • Johansson R.S.
        • Birznieks I.
        First spikes in ensembles of human tactile afferents code complex spatial fingertip events.
        Nat Neurosci. 2004; 7: 170-177
        • Saal H.P.
        • Wang X.
        • Bensmaia S.J.
        Importance of spike timing in touch: an analogy with hearing?.
        Curr Opin Neurobiol. 2016; 40: 142-149
        • Simone D.A.
        • Sorkin L.S.
        • Oh U.
        • Chung J.M.
        • Owens C.
        • LaMotte R.H.
        • Willis W.D.
        Neurogenic hyperalgesia: central neural correlates in responses of spinothalamic tract neurons.
        J Neurophysiol. 1991; 66: 228-246
        • Melzack R.
        • Wall P.D.
        Pain mechanisms: a new theory.
        Science. 1965; 150: 971-979
        • Edhi M.M.
        • Heijmans L.
        • Vanent K.N.
        • Bloye K.
        • Baanante A.
        • Jeong K.S.
        • Leung J.
        • Zhu C.
        • Esteller R.
        • Saab C.Y.
        Time-dynamic pulse modulation of spinal cord stimulation reduces mechanical hypersensitivity and spontaneous pain in rats.
        Sci Rep. 2020; 1020358
        • Al-Kaisy A.
        • Van Buyten J.-P.
        • Smet I.
        • Palmisani S.
        • Pang D.
        • Smith T.
        Sustained effectiveness of 10 kHz high-frequency spinal cord stimulation for patients with chronic, low back pain: 24-month results of a prospective multicenter study.
        Pain Med. 2014; 15: 347-354
        • Paz-Solis J.
        • Thomson S.
        • Jain R.
        • Chen L.
        • Huertas I.
        • Doan Q.
        Exploration of high and low frequency options for subperception spinal cord stimulation using neural dosing parameter relationships: the HALO study.
        Neuromodulation. 2021;
        • Vesper J.
        • Slotty P.
        • Schu S.
        • Poeggel-Kraemer K.
        • Littges H.
        • Van Looy P.
        • Agnesi F.
        • Venkatesan L.
        • Van Havenbergh T.
        Burst SCS microdosing is as efficacious as standard burst SCS in treating chronic back and leg pain: results from a randomized controlled trial.
        Neuromodulation: Technol Neural Interface. 2019; 22: 190-193
        • Hagedorn J.M.
        • Engle A.M.
        • Ghosh P.
        T RD. Device profile of the Proclaim XR neurostimulation system for the treatment of chronic pain: an overview of its safety and efficacy.
        Expet Rev Med Dev. 2020; 17: 499-505
        • Leitner J.
        • Westerholz S.
        • Heinke B.
        • Forsthuber L.
        • Wunderbaldinger G.
        • Jager T.
        • Gruber-Schoffnegger D.
        • Braun K.
        • Sandkuhler J.
        Impaired excitatory drive to spinal GABAergic neurons of neuropathic mice.
        PLoS One. 2013; 8e73370
        • Gilbert J.E.
        • Zhang T.
        • Esteller R.
        • Grill W.M.
        Distributed network model of effects of spinal cord stimulation (SCS) on wide-dynamic range (WDR) neurons.
        Proc Soc Neurosci. 2018; ([Abstract 216.11/X7])
        • Kajander K.C.
        • Bennett G.J.
        Onset of a painful peripheral neuropathy in rat: a partial and differential deafferentation and spontaneous discharge in Aβ and Aδ primary afferent neurons.
        J Neurophysiol. 1992; 68: 734-744
        • Liu X.
        • Eschenfelder S.
        • Blenk K.-H.
        • Jänig W.
        • Häbler H.-J.
        Spontaneous activity of axotomized afferent neurons after L5 spinal nerve injury in rats.
        Pain. 2000; 84: 309-318
        • Wall P.D.
        • Gutnick M.
        Ongoing activity in peripheral nerves: the physiology and pharmacology of impulses originating from a neuroma.
        Exp Neurol. 1974; 43: 580-593
        • Zhang T.C.
        • Janik J.J.
        • Grill W.M.
        Modeling effects of spinal cord stimulation on wide-dynamic range dorsal horn neurons: influence of stimulation frequency and GABAergic inhibition.
        J Neurophysiol. 2014; 112: 552-567
        • Hillman P.
        • Wall P.D.
        Inhibitory and excitatory factors influencing the receptive fields of lamina 5 spinal cord cells.
        Exp Brain Res. 1969; 9: 284-306
        • Kumar K.
        • Taylor R.S.
        • Jacques L.
        • Eldabe S.
        • Meglio M.
        • Molet J.
        • Thomson S.
        • O'Callaghan J.
        • Eisenberg E.
        • Milbouw G.
        • Buchser E.
        • Fortini G.
        • Richardson J.
        • North R.B.
        Spinal cord stimulation versus conventional medical management for neuropathic pain: a multicentre randomised controlled trial in patients with failed back surgery syndrome.
        Pain. 2007; 132: 179-188
        • Sandkuhler J.
        The role of inhibition in the generation and amplification of pain.
        in: Castro-Lopes J. Current topics in pain: 12th world congress on pain seattle. IASP Press, WA, USA2009
        • Moore K.A.
        • Kohno T.
        • Karchewski L.A.
        • Scholz J.
        • Baba H.
        • Woolf C.J.
        Partial peripheral nerve injury promotes a selective loss of GABAergic inhibition in the superficial dorsal horn of the spinal cord.
        J Neurosci. 2002; 22: 6724-6731
        • Hines M.
        • Carnevale N.
        The NEURON simulation environment.
        Neural Comput. 1997; 9: 1179-1209
        • Cassar I.R.
        • Titus N.D.
        • Grill W.M.
        An improved genetic algorithm for designing optimal temporal patterns of neural stimulation.
        J Neural Eng. 2017; 14
        • North R.B.
        • Ewend M.G.
        • Lawton M.T.
        • Piantadosi S.
        Spinal cord stimulation for chronic, intractable pain: superiority of "multi-channel" devices.
        Pain. 1991; 44: 119-130
        • D'Mello R.
        • Dickenson A.H.
        Spinal cord mechanisms of pain.
        Br J Anaesth. 2008; 101: 8-16
        • Dorval A.D.
        Estimating neuronal information: logarithmic binning of neuronal inter-spike intervals.
        Entropy. 2011; 13: 485-501
        • Gilbert J.E.
        • Zhang T.C.
        • Esteller R.
        • Grill W.M.
        Pharmacological disinhibition degrades neural population responses to Aβ-fiber electrical stimulation (Aβ-ES).
        in: Proceedings of the neural interfaces 2021: the NANS-NIC joint virtual meeting. 2021
        • Zhang T.C.
        • Janik J.J.
        • Peters R.V.
        • Chen G.
        • Ji R.R.
        • Grill W.M.
        Spinal sensory projection neuron responses to spinal cord stimulation are mediated by circuits beyond gate control.
        J Neurophysiol. 2015; 114: 284-300
        • Decosterd I.
        • Woolf C.J.
        Spared nerve injury: an animal model of persistent peripheral neuropathic pain.
        Pain. 2000; 87: 149-158
        • Bonin R.P.
        • Bories C.
        • De Koninck Y.
        A simplified up-down method (SUDO) for measuring mechanical nociception in rodents using von Frey filaments.
        Mol Pain. 2014; 10: 26
        • Shimazaki H.
        • Shinomoto S.
        A method for selecting the bin size of a time histogram.
        Neural Comput. 2007; 19: 1503-1527
        • Dorval A.D.
        • Grill W.M.
        Deep brain stimulation of the subthalamic nucleus reestablishes neuronal information transmission in the 6-OHDA rat model of parkinsonism.
        J Neurophysiol. 2014; 111: 1949-1959
        • Cutts C.S.
        • Eglen S.J.
        Detecting pairwise correlations in spike trains: an objective comparison of methods and application to the study of retinal waves.
        J Neurosci. 2014; 34: 14288-14303
        • Deuis J.R.
        • Dvorakova L.S.
        • Vetter I.
        Methods used to evaluate pain behaviors in rodents.
        Front Mol Neurosci. 2017; 10: 284
        • Lee K.Y.
        • Ratté S.
        • Prescott S.A.
        Excitatory neurons are more disinhibited than inhibitory neurons by chloride dysregulation in the spinal dorsal horn.
        Elife. 2019; 8e49753
        • Snyder A.C.
        • Morais M.J.
        • Smith M.A.
        Dynamics of excitatory and inhibitory networks are differentially altered by selective attention.
        J Neurophysiol. 2016; 116: 1807-1820
        • Wobbrock J.O.
        • Findlater L.
        • Gergle D.
        • Higgins J.J.
        The aligned rank transform for nonparametric factorial analyses using only anova procedures.
        in: Proceedings of the SIGCHI conference on human factors in computing systems. 2011: 143-146
        • North R.B.
        • Kidd D.H.
        • Zahurak M.
        • James C.S.
        • Long D.M.
        Spinal cord stimulation for chronic, intractable pain: experience over two decades.
        Neurosurgery. 1993; 32: 384-395
        • Howell B.
        • Lad S.P.
        • Grill W.M.
        Evaluation of intradural stimulation efficiency and selectivity in a computational model of spinal cord stimulation.
        PLoS One. 2014; 9: 1-25
        • Rigoard P.
        • Basu S.
        • Desai M.
        • Taylor R.
        • Annemans L.
        • Tan Y.
        • Johnson M.J.
        • Van den Abeele C.
        • North R.
        • Group P.S.
        Multicolumn spinal cord stimulation for predominant back pain in failed back surgery syndrome patients: a multicenter randomized controlled trial. pain.
        2019
        • Lempka S.F.
        • Zander H.J.
        • Anaya C.J.
        • Wyant A.
        • Ozinga JGt
        • Machado A.G.
        Patient-specific analysis of neural activation during spinal cord stimulation for pain.
        Neuromodulation. 2020; 23: 572-581
        • Guo Y.
        • Rubin J.E.
        • McIntyre C.C.
        • Vitek J.L.
        • Terman D.
        Thalamocortical relay fidelity varies across subtahlamic nucleus deep brain stimulation protocols in a data-driven computational model.
        J Neurophysiol. 2008; 99: 1477-1492
        • Prescott S.A.
        • Ratte S.
        Pain processing by spinal microcircuits: afferent combinatorics.
        Curr Opin Neurobiol. 2012; 22: 631-639
        • Seal R.P.
        Illuminating the gap: neuronal cross-talk within sensory ganglia and persistent pain.
        Neuron. 2016; 91: 950-951
        • Coull J.A.M.
        • Boudreau D.
        • Bachand K.
        • Prescott S.A.
        • Nault F.
        • Sík A.
        • De Koninck P.
        • De Koninck Y.
        Trans-synaptic shift in anion gradient in spinal lamina I neurons as a mechanism of neuropathic pain.
        Nature. 2003; 424: 938-942
        • Peirs C.
        • Seal R.P.
        Neural circuits for pain: recent advances and current views.
        Pain Research. 2016; 354: 578-584
        • Metzger C.S.
        • Hammond M.B.
        • Paz-Solis J.F.
        • Newton W.J.
        • Thomson S.J.
        • Pei Y.
        • Jain R.
        • Moffitt M.
        • Annecchino L.
        • Doan Q.
        A novel fast-acting sub-perception spinal cord stimulation therapy enables rapid onset of analgesia in patients with chronic pain.
        Expet Rev Med Dev. 2021; 18: 299-306