SUPRASPINAL CONTROL OF MOTONEURONS AFTER PARALYSIS ENABLED BY SPINAL CORD STIMULATION

Abstract

Spinal cord stimulation (SCS) restores motor control after spinal cord injury (SCI) and stroke. This evidence led to the hypothesis that SCS facilitates residual supraspinal inputs to spinal motoneurons. Instead, here we show that SCS does not facilitate residual supraspinal inputs but directly triggers motoneurons action potentials. However, supraspinal inputs can shape SCS-mediated activity, mimicking volitional control of motoneuron firing. Specifically, by combining simulations, intraspinal electrophysiology in monkeys and single motor unit recordings in humans with motor paralysis, we found that residual supraspinal inputs transform subthreshold SCS-induced excitatory postsynaptic potentials into suprathreshold events. We then demonstrated that only a restricted set of stimulation parameters enables volitional control of motoneuron firing and that lesion severity further restricts the set of effective parameters. Our results explain the facilitation of voluntary motor control during SCS while predicting the limitations of this neurotechnology in cases of severe loss of supraspinal axons.

Extended Data Fig. 1 |. Motoneuron firing is robust to parameter changes after a lesion.

a-c, Changes up to ±10% in key parameters of the model (mean synaptic weight from residual supraspinal fibers to motoneurons (a), mean synaptic weight from sensory fibers to motoneurons (b) and the conductivity of the N-type Ca2+ ion channel (c)) do not produce dramatic changes in motoneuron firing rate, thereby indicating that the model is robust to changes in key parameters. d, Normalized ISI probability distribution for the different modifications of the key parameters in a-c. The transformation of subthreshold SCS-pulses in suprathreshold events do not depend on the fine tuning of these parameters.

Extended Data Fig. 2 |. Motoneuron ion channel dynamics in our biophysical model.

a, Ion currents of a single motoneuron in response to different SCS frequencies (same amplitude, 50%) along with the motoneuron membrane potential during concurrent innervation of residual supraspinal inputs (gray). Ionic currents: L-type Ca2+ (Ical), Ca2+-activated K (Ikca), delayed rectifier K+ (Ikrect), N-type Ca2+ (Ican), and nonlinear fast Na+ (Ina). At low stimulation frequencies (<19 Hz), motoneurons generate action potentials with concurrent supraspinal input for each stimulation pulse. At higher frequencies, motoneurons skip some stimulation pulses due to N-type Ca2+ slow time dynamics. b, Ion currents of the same motoneuron at 45 Hz SCS (50% amplitude) for different conductances of the N-type Ca2+ ion channel. Setting the conductance close to 0 relaxes N-type Ca2+ dynamics constraints. Consequently, less rectifier K+ ions leave the membrane to stabilize the motoneuron membrane potential, leading to a shorter refractory period in the overall motoneuron membrane and producing more action potentials following each stimulation pulse. c, Ion current of the N-type Ca2+ ion channel varying its conductance and its resulting membrane potential. When membrane dynamics allow the entrance of N-type Ca2+ ions, a cascade of Ca2+-activated K ions respond proportionally to counteract this entrance, leading to exaggerated refractory periods in the motoneuron membrane. Arrows indicate stimulation pulses.

Extended Data Fig. 3 |. Motoneuron firing rate increases non-monotonically with SCS frequency.

a, Normalized ISI probability distribution of a motoneuron pool prior to the lesion (i.e., without SCS) showing no peaks at integer numbers. For comparison purposes with panel b, the time is set in SCS pulses units with a stimulation frequency of 45 Hz. b, Normalized ISI probability distribution of the same motoneuron pool during SCS alone (45 Hz, 50% amplitude). Spikes reliably follow SCS pulses. c, Firing rates of a pool of motoneurons as a function of SCS frequencies (constant 60% amplitude). Motoneurons are concurrently receiving excitatory inputs from residual supraspinal fibers. d-g, Normalized ISI probability distribution of a motoneuron pool in SCS pulses units for different SCS frequencies: 19 (d), 29 (e), 45 (f) and 60 (g) Hz (note arrows in c). The different peaks in the Normalized ISI probability distribution is a consequence of transforming SCS pulses into action potentials in the motoneuron membrane. In general, motoneuron firing rate increases with SCS frequency (c). However, when motoneurons are unable to produce an action potential for each SCS pulse, the firing rate drops (see arrows e and g in panel c). Given the membrane dynamics, following an action potential, there is a refractory period that sets the ISI greater than the interpulse interval (IPI) and motoneurons skip SCS pulses. This skipping effectively decreases the firing rate of the motoneuron pool. Hence, in the local maximum marked with the d arrow, the ISI is the same as the IPI (d) (i.e., all SCS pulses trigger an action potential) while the minimum coincides where almost all motoneurons skip one pulse (e). However, the IPI decreases with higher SCS frequencies, producing an increase in firing rate between point (e) and (g). There is another local maximum in firing rate when motoneurons skip two SCS pulses (f), while again the minimum coincides where most motoneurons skip two pulses (g). Finally, at increasing SCS frequencies, the firing rate of the pool increases consequently.

Extended Data Fig. 4 |. Simulations of SCS parameter regimes in an isometric task.

a, Simulation of a pool of motoneuron pool during a maximum voluntary contraction task during phasic activation of supraspinal fibers (2 seconds active, 2 seconds inactive). Traces represent the normalized force (top), motoneuron pool firing rate (middle), supraspinal input count (bottom) in a pre- (gray) and postlesion (black) system. b, Same as c during continuous SCS at parameters belonging to the voluntary movement regime (SCS frequency 60 Hz, amplitude 30%). Motoneurons fire only coinciding with supraspinal activation, achieving a firing rate only 19% smaller than the prelesion case. Bottom trace indicates the continuous SCS amplitude. c, Same as b during continuous SCS at parameters belonging to the involuntary movement regime (SCS frequency 60 Hz, amplitude 80%). Motoneurons fire also when supraspinal fibers are silent.

Extended Data Fig. 5 |. Monkey EMG response latencies and intraspinal ISI distributions.

a, EMG response latencies for each monkey evoked by a single SCS pulse (H-reflex, blue) and ICS pulse (motor-evoked potential (MEP), red). b, Histograms of ISI normalized probability distribution of the sorted motor units during 30 and 60 Hz SCS in both monkeys. Square and error bars indicate the mean distribution of the data with 95% confidence interval.

Extended Data Fig. 6 |. Lesion characterization and SCS implantation in humans with paralysis.

a, Sagittal, coronal, and axial T1-weighted magnetic-resonance imaging (MRI) 2D projections for SCS03 and SCS04 using a 3T Prisma system (Siemens) using a 64-channel head and neck coil (T1-weighted structural image was captured using a magnetization-prepared rapid gradient echo sequence: repetition time = 2,300 ms; echo time = 2.9 ms; field of view = 256 × 256 mm2; 192 slices, slice thickness = 1.0 mm, in-plane resolution = 1.0 × 1.0 mm). The segmented lesion is shown in red for both participants and performed as described in. R indicates the Right hemisphere. b, X-rays for SCS03 and SCS04 showing the position of the spinal leads. The red lines mark the same anatomical location across the X-rays to facilitate interpretation. Minimal displacement occurred after initial implantation. c-d, X-rays for STIM01 showing the position of the spinal leads and the lesion. The red lines mark the same anatomical location across the X-rays to facilitate interpretation. Minimal displacement occurred after initial implantation.

Extended Data Fig. 7 |. Action potentials triggered by SCS pulses are not stimulation artifacts in humans with paralysis.

a, For the experiments with STIM01, we recorded the stimulation artifact with a surface EMG placed on the dorsal. We defined the stimulation timestamps through EMG threshold-crossing above 0.05 mV (dashed line). b, For the experiments with the stroke participants, traces were recorded from high-density EMG: in gray showing the presence of stimulation artifacts and in black after removing the stimulation artifacts (showing an example of SCS03, 60 Hz SCS). c, For the experiments with the stroke participants, stimulation artifact removal was achieved by applying a series of notch filters at harmonics of SCS frequency to remove peaks in the Fast Fourier transform (shown in gray) of the bandpass filtered signal (see methods). d, Examples of the probability distribution of the time between the most recent SCS pulse and an action potential for each motor unit (MU; 40 Hz SCS, 1 ms bin) using surface EMG in STIM01. If the action potentials detected by the acquisition software were stimulation artifacts, they would be aligned to the SCS pulse, thereby showing histograms with a high peak at time 0. Given that the times from the most recent SCS pulse and the action potentials are distributed with a mean between 12.5 and 25 ms, the action potential detected can not be artifacts of SCS. e, Same as d using the high-density EMG in the stroke participants (SCS03, SCS04; 60 Hz SCS, 1 ms bin).

Extended Data Fig. 8 |. Experimental validation in humans with spinal cord injury and stroke.

a, Torque traces (top, overlapped with a trace from another repetition), firing rate (middle) and raster plot of the motor units (bottom) followed by histogram of ISI normalized probability distribution of the extracted motor units in STIM01 for all repetitions during the maximum voluntary contraction task during SCS (40 Hz) at increasing amplitudes (7 mA, 10.8 mA, 12.5 mA; 2 repetitions for each stimulation amplitude). At 7 mA SCS, no motor units were detected. b, Histogram of ISI normalized probability distribution of the extracted motor units in STIM01 during 40 Hz SCS (above, 2 repetitions) and 60 Hz SCS (below, 2 repetitions) recorded in other sessions when performing the isometric task at 35% of maximum force. c, Histograms of ISI probability distribution of the extracted motor units in SCS04 and SCS03 with and without stimulation for all force levels for all trials (18 repetitions for SCS04 for all for levels, 9 repetitions in SCS03 for 5% isometric force level, 7 for 20%, 7 for 35% and 4 for 50%).

Extended Data Fig. 9 |. Simulation of SCS biophysics in a severe lesion.

a, Root-mean square (RMS) of the difference between the motoneuron firing rate pre- vs postlesion for all simulated SCS parameters in a severe lesion (same as in Fig. 4j, but with a system with 11 supraspinal fibers). As with isometric force simulations in Fig. 7b, there is a shrinking of the voluntary movement regime where small task errors occur for a very restrained combination of SCS parameters. The gray lines in the heatmap indicate the voluntary movement regime for this severe lesion (Fig. 7b). b, Simulation of the normalized force and motoneuron pool firing rate during phasic activation of supraspinal fibers (2 seconds active, 2 seconds inactive) and continuous SCS at parameters belonging to the involuntary movement regime. Residual supraspinal fibers sculpt SCS drive ro produce functional motoneuron activity. Motoneurons firing is inhibited by the phasic activation of residual supraspinal fibers (red). c, Same as b, but SCS initiates 500 ms after the initiation of supraspinal fibers (red). Motoneurons now fire only coinciding with temporal activation of SCS.

Extended Data Fig. 10 |. EMG responses during concurrent stimulation in monkeys.

a, Root-mean square (RMS) of the EMG responses during phasic ICS and continuous SCS at increasing SCS amplitudes in addition to Fig. 7e (APB in Mk-Ka, FDM in Mk-JC). Statistical quantification of the RMS EMG responses (***p<0.001; **p<0.01; *p<0.05; Kruskal-Wallis test with 21, 14 points for 500 uA SCS at 0.7 mA ICS and at 0.5 mA ICS, respectively, in APB in Mk-Ka; 15, 18 points for 600 uA in Mk-Ka; 16, 21 points for 700 uA in Mk-Ka; 8, 19 points for 850 uA in Mk-Ka; 4, 8 points for 80 uA, 60 Hz SCS in Mk-JC; 11, 11 points for 110 uA, 60 Hz SCS in Mk-JC; 9, 8 points for 140 uA, 60 Hz SCS in Mk-JC; 8, 5 points for 300 uA, 60 Hz SCS in Mk-JC; 9, 10, 10 points for 80 uA, 30 Hz SCS in APB in Mk-JC; 8, 10, 8 points for 80 uA, 30 Hz SCS in FDM Mk-JC). b, RMS of the EMG responses during continuous SCS alone at increasing SCS amplitudes of the muscles shown in a (above) and Fig. 7e (below). Statistical quantification of the RMS EMG responses (***p<0.001; **p<0.01; *p<0.05; Kruskal-Wallis test with 27, 21, 19, 17 and 9 points for no stim, 500, 600, 700, 850 uA SCS, respectively, in APB in Mk-Ka; 10, 5, 11, 9, 8 points in FDM in Mk-JC; 27, 21, 18, 15, 10 points in FDM in Mk-Ka; 10, 5, 11, 9, 7 points in APB in Mk-JC). RMS EMG responses appeared to be significant in Mk-Ka, however, these responses were still below the lowest RMS value measured during concurrent phasic ICS.

Fig. 1 |. Biophysical model with underlying assumptions.

a, Schematic of the spinal circuit involved in the recovery of voluntary movement during SCS composed of the monosynaptic Ia-to-motoneuron reflex circuit and supraspinal monosynaptic connections. Spinal motoneurons receive excitatory inputs from residual supraspinal fibers and Ia afferents recruited by SCS. b, Four assumptions grounded on experimental evidence that constitute the biophysical model of spinal motoneurons. c, Three-element biophysical model: motoneuron, Ia afferents inputs and supraspinal inputs. Biophysical model of spinal alpha motoneuron encompasses an electronic-equivalent dendritic tree, multicompartment soma and myelinated axon, connected through an initial axon segment. Motoneuron membrane dynamics follow dedicated Hodgkin–Huxley equations for multiple ion channels: L-type Ca2+, Ca2+-activated K, delayed rectifier K+, N-type Ca2+, and nonlinear fast Na+.

Fig. 2 |. Supraspinal input facilitates SCS.

a, Simulation of the membrane potential of a single motoneuron driven by supraspinal fibers prior to the lesion (top), raster plot of each supraspinal fiber firing stochastically (middle), ionic current of N-type Ca2+ of the motoneuron (bottom). Motoneuron firing rate is limited by the slow calcium dynamics (yellow). b, Simulation of the membrane potential of a single motoneuron after the lesion (white background), during SCS alone (16 Hz, 50% amplitude; blue) and concurrent residual supraspinal inputs and SCS (gray, top), raster plot of each residual supraspinal fiber firing stochastically after the lesion (middle), raster plot of each Ia afferent firing synchronously after each SCS pulse (bottom). c, Same as b at a higher SCS frequency (45 Hz, 50% amplitude). d, Mean motoneuron pool firing rate as a function of SCS frequencies, coupled up to 19 Hz (i.e., motoneuron firing rate equals stimulation frequency up to 19 Hz SCS represented by the dashed gray line) and effectively decoupled at higher stimulation frequencies while still increasing motoneuron firing. e, Probability distribution of the motoneuron pool ISI normalized by the stimulation interpulse interval (16 Hz SCS), indicating the number of SCS pulses between action potentials (in this case, 1 SCS pulse). Inset indicates the membrane potential of a motoneuron during concurrent supraspinal inputs and SCS, highlighting in yellow the time constraint of the slow calcium dynamics. f, Same as e at a higher SCS frequency (45 Hz). Arrows indicate stimulation pulses. Light arrows indicate skipped pulses elicited within the slow calcium dynamics in inset in f.

Fig. 3 |. SCS parameter regimes.

a, Simulated heatmap of mean firing rates of a motoneuron pool during SCS alone (5 second simulation) for a range of frequencies and amplitudes. Solid line separates subthreshold from suprathreshold SCS parameters. b, Same as a during concurrent SCS and residual supraspinal inputs innervating the pool of motoneurons. Dashed line separates the combinations of parameters that produce sustained motoneuron firing rates (>8 Hz). The parameters between the grey lines represent the combination of SCS parameters that produce sustained firing rate only with concurrent SCS and residual supraspinal inputs (i.e., voluntary movement regime).

Fig. 4 |. Residual supraspinal inputs can modulate motoneuron firing rate during SCS.

a, Motoneuron pool firing rate (5 seconds simulations) as a function of residual supraspinal inputs firing rate in the voluntary movement regime (SCS parameters: 45 Hz and 30% amplitude) and involuntary regime (SCS parameters: 45 Hz and 80% amplitude). b, Membrane potential of the same motoneuron with residual supraspinal inputs at 20 and 45 Hz. This motoneuron is recruited at 45 Hz but not at 20 Hz. c Membrane potential of another motoneuron that skips 3 SCS pulses at 20 Hz and 4 or 5 at 45 Hz. SCS parameters: 69 Hz and 30% amplitude. d, Percentage of motoneurons recruited (with a firing rate >8 Hz) as a function of the residual supraspinal fibers firing rate e,f, Two examples of ISI normalized probability distributions for residual supraspinal firing rates at 20 and 45 Hz during SCS (69 Hz, 30%). g, ISI normalized cumulative distributions show that residual supraspinal inputs can modulate motoneuron firing rate by controlling the number of skipped SCS pulses. h-i, Fine movement task in the voluntary (69 Hz, 30% amplitude) and involuntary (69 Hz, 80% amplitude) movement regimes. Firing rates are computed in 100 ms bin windows. From bottom to top: supraspinal inputs that follow a sinusoidal wave, motoneuron firing rate and normalized force. Blue area indicates the difference in firing rate from the prelesion case (i.e., error). j, Root-mean square (RMS) of the difference between the motoneuron firing rate pre- vs postlesion (dark blue regions in h and i) for all simulated SCS parameters. The gray lines in the heatmap indicate the voluntary movement regime defined in the isometric task (Fig. 3b).

Fig. 5 |. Experimental validation of model assumptions and predictions in monkeys.

a, Experimental setup in anesthetized monkeys and schematic of the spinal circuits involved during dorsal root stimulation and activation of corticospinal fibers via deep-brain stimulation of the internal capsule. In response to the stimulation, EMG hand muscles responses were recorded and intraspinal spikes captured. b, Schematic of the paired stimulation protocol. A single pulse of SCS (top) is facilitated with ICS conditioning at different delays (bottom) and vice versa, yielding greater peak-to-peak amplitudes of the EMG traces of the conditioned pulse. c, Stimulation-triggered average of the EMG responses when conditioning with ICS a single pulse of SCS at different delays in FDM muscle in Mk-JC. d, Peak-to-peak amplitudes of the EMG responses in the hand muscles (FDM, APB) conditioned with ICS in both monkeys. e, Same as c when conditioning with SCS a single pulse of ICS. f, Same as d conditioning with SCS. Statistical quantification of the peak-to-peak amplitudes (***p<0.001; **p<0.01; *p<0.05; Kruskal-Wallis test with 117, 116, 119, 113, 121 points for FDM muscle for control, 5 ms, 10 ms, 50 ms, 100 ms delay during ICS conditioning in Mk-JC; 118, 96, 113, 100, 119 points for APB muscle during ICS conditioning in Mk-JC; 58, 64, 61, 53, 58 points for FDM muscle during iCS conditioning in Mk-Ka; 71, 60, 60, 49, 61 points for APB muscle during iCS conditioning in Mk-Ka. 118, 119, 116, 118, 117 points for FDM muscle for control, 5 ms, 10 ms, 50 ms, 100 ms delay during SCS conditioning in Mk-JC; 117, 97, 119, 98, 114 points for APB muscle during SCS conditioning in Mk-JC; 59, 60, 59, 67, 58 points for FDM muscle during SCS conditioning in Mk-Ka; 58, 60, 60, 71, 58 for APB muscle during SCS conditioning in Mk-Ka). g, Samples of a raster plot corresponding to sorted motor unit waveforms recorded from the ventral channels of the linear probe (motor units sorted from a 30-second trial). h, Histograms of normalized ISI probability distribution of the sorted motor units during 100 Hz SCS in both monkeys. i, Firing rates of the sorted motor units for different frequencies in Mk-JC (left) and Mk-Ka (right). Statistical quantification of the firing rates (***p<0.001; **p<0.01; *p<0.05; Kruskal-Wallis test with 22, 26 and 21 points for 30 Hz (0.14 mA), 60, Hz (0.14 mA) and 100 Hz (0.11 mA) in Mk-JC. 14, 13 and 12 points for in Mk-Ka (0.7 mA SCS). j, Normalized ISI cumulative distribution of the sorted motor units during continuous subthreshold SCS at 30, 60 and 100 Hz. Square and error bars indicate the mean distribution of the data with 95% confidence interval.

Fig. 6 |. Experimental validation in humans with paralysis.

a, Schematic of the experimental setup in participants with upper-limb paralysis performing a force task with their paretic arm (same setup was adapted for the participant with lower-limb paralysis). b, Schematic of motor unit decomposition using surface or high-density EMG performing a voluntary movement along with a sample raster plot for a detected motor unit waveform. c, Example of torque traces (top, overlapped with a trace from another repetition), firing rate (middle) and raster plot of the motor units (bottom) of both extensor muscles (Vastus Lateralis, Rectus Distal, sorted by average firing rate) during the isometric force task (35% of maximum force) in STIM01 during 40 Hz SCS. d, Histogram of ISI normalized probability distribution of the extracted motor units in STIM01 during 40 Hz SCS (above, 2 repetitions) and 60 Hz SCS (below, 2 repetitions). e, Example of torque traces (top, overlapped with trace from all repetitions), firing rate (middle) and raster plot of the motor units (bottom) during the force-modulated task (35–50% of the maximum force) in SCS04 with (60 Hz SCS) and without stimulation. f, Histograms of ISI normalized probability distribution of the extracted motor units in SCS04 (top, 18 repetitions) and SCS03 (bottom, 4 repetitions) with and without stimulation during 50% force-level period. g, Schematic of the parametrized isometric force-modulated task between pairs of force levels (force level is represented as percent of maximum force). Its lower and upper-limits were modified to perform a 5–20% maximum force task, 20–35% and 35–50%. h, Normalized ISI cumulative distribution of the extracted motor units during SCS in SCS04 (left) and SCS03 (right) during SCS for all force levels for all trials (18 repetitions for SCS04 for all for levels, 9 repetitions in SCS03 for 5% isometric force level, 7 for 20%, 7 for 35% and 4 for 50%).

Fig. 7 |. Lesion severity reduces voluntary movement regime.

a, Simulation of the membrane potential of a single motoneuron driven by residual supraspinal fibers and SCS (16 Hz, 30%) in a mild (30 supraspinal fibers, top) and severe lesion (11 supraspinal fibers, bottom). b, Simulated heatmap indicating the combinations of SCS parameters within the voluntary movement regime for different degrees of lesion severity (mildest lesion in light red, the most severe lesion in dark red; 5 seconds simulation). c, Number of SCS parameter combinations within the voluntary movement regime for different degrees of lesion severity. d, Raw EMG and torque traces to illustrate the stimulation protocol employing phasic ICS (47 Hz) during continuous subthreshold SCS (60 Hz) at increasing SCS amplitudes (in Mk-JC, APB muscle). e, Root-mean square (RMS) of the EMG responses during phasic ICS and continuous SCS at increasing SCS amplitudes. Statistical quantification of the RMS EMG responses (***p<0.001; **p<0.01; Kruskal-Wallis test with 5, 5 points for 80 uA for 0.7 mA ICS and 0.5 mA ICS in Mk-JC; 10, 11 points for 110 uA; 9, 9 points for 140 uA; 7, 7 points for 300 uA; 21, 16 points for 500 uA for 0.7 mA ICS and 0.5 mA ICS in Mk-Ka; 18, 11 points for 600 uA; 17, 21 points for 700 uA; 9, 21 points for 850 uA). f, Mean torque during phasic ICS and continuous SCS at increasing SCS amplitudes. Statistical quantification of the mean torque (***p<0.001; **p<0.01; *p<0.05; Kruskal-Wallis test with 5, 9 points for 80 uA for 0.7 mA ICS and 0.5 mA ICS in Mk-JC; 10, 11 points for 110 uA; 9, 9 points for 140 uA; 8, 5 points for 300 uA; 18, 15 points for 500 uA for 0.7 mA ICS and 0.5 mA ICS in Mk-Ka; 18, 20 points for 600 uA; 16, 20 points for 700 uA; 9, 20 points for 850 uA). Square and error bars indicate the mean distribution of the data with 95% confidence interval.