Supraspinal control of motoneurons after paralysis enabled by spinal cord stimulation

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.


INTRODUCTION
In the last two decades, a series of animal and clinical studies showed that epidural electrical spinal cord stimulation (SCS) can restore some degree of motor control of lower and upper limbs in people with spinal cord injury and stroke [1][2][3][4][5][6][7] .The most striking observation that emerged from clinical trials is the evidence that, during SCS, humans with motor paralysis appear able to voluntarily control their previously paralyzed limbs 2,8 .In other words, SCS does not generate movement per se, but rather facilitates volitional movement.In an effort to explain these results, mechanistic investigations on SCS focused either on the identi cation of the neural pathways that are directly recruited by SCS [9][10][11][12] or the involvement of movement-generating spinal circuits such as re exes [10][11][12] , synergies 13 and central pattern generators [14][15][16] .However, these studies never considered the contribution of residual supraspinal inputs, thereby failing to explain the facilitation of voluntary movements.This poses a problem because the current design of SCS protocols is focused on the recruitment of spinal neural structures 5,11,17 , without considering the role of volitional inputs.Understanding how volitional inputs can shape motor output during SCS could lead to more effective stimulation protocols targeted to maximize residual volitional input.
One of the most popular views into how SCS produces this facilitation is that electrical stimulation exerts a neuromodulatory effect on spinal circuits 2,8,18 .More speci cally, SCS would increase the excitability of spinal circuits and, particularly, spinal motoneurons, making them more responsive to volitional residual supraspinal inputs 2,8 .In other words, through this neuromodulatory mechanism, SCS would enable descending bers that survive the injury to regain control of motoneurons below the lesion.Unfortunately, this conjecture only provides a phenomenological description of the effects of SCS but fails to propose an experimentally veri able mechanism of action by which SCS restores motor function.
In order to transform this conjecture into an experimentally veri able hypothesis, we started from identi ed properties and mechanisms of SCS.It is well-established that SCS directly activates sensory afferents in the dorsal roots and dorsal columns [9][10][11]19,20 . Imporantly, these afferents form mono-and polysynaptic excitatory connections to spinal motoneurons and convey SCS-mediated inputs to motoneurons.Thus, stimulation of the sensory afferents could provide motoneurons with additional excitatory drive after paralysis.Indeed, studies in monkeys showed that SCS pulses potentiate monosynaptic corticospinal motor-evoked potentials 21 .Here, we started from this body of electrophysiological evidence to study the effects of excitatory sensory inputs elicited by SCS on volitional control of spinal motoneurons.
We built on established biophysics of motoneuron membrane dynamics [22][23][24][25] to calculate their membrane potential while receiving SCS-mediated excitatory drive from sensory afferents and excitatory inputs from residual supraspinal bers.We found that motoneurons can re during voluntary drive by transforming subthreshold sensory inputs generated by SCS pulses into suprathreshold events, thereby mimicking voluntary activity.However, the SCS parameter space at which this phenomenon occurred shrunk as a function of lesion severity.We validated these predictions with electrophysiology in two monkeys, two persons with mild stroke and one person with severe spinal cord injury, who were implanted with SCS leads.Our results suggest that SCS does not facilitate voluntary drive.Instead, SCS directly generates motoneuron action potentials, but only when residual supraspinal drive is present.

Fundamental representation of the biophysics of SCS
To identify the membrane mechanisms underlying motoneuron recruitment during SCS (Fig. 1a), we started from four assumptions grounded on experimental evidence (Fig. 1b).First, SCS directly recruits large-diameter sensory afferents, particularly, Ia afferents 10 .Second, electrical stimulation causes simultaneous depolarization in all recruited afferents, thereby producing non-natural volleys of action potentials.Third, in contrast to this non-natural input, supraspinal drive can be modeled as a stochastic process of non-synchronized axon ring 26,27 .Fourth, a lesion reduces the number of available excitatory supraspinal bers, thus impacting the ability of the cortex to produce action potentials in spinal motoneurons.
We then designed the simplest known circuit involved in voluntary motor control: motoneurons, Ia afferents and supraspinal bers (Fig. 1c).In this three-element model, the membrane potential of spinal motoneurons is described with a modi ed Hodgkin-Huxley biophysical representation 25 .Instead, Ia afferents and supraspinal bers are treated as point-processes with given ring rates connected to synaptic boutons, triggering excitatory postsynaptic potentials (EPSPs) in the motoneurons.
We modeled the stimulation intensity as the number of Ia afferents activated by SCS, allowing us to vary between no activation at all and recruitment of the entire Ia population 24 .Recruited Ia afferents elicit synchronous volleys that re at a given SCS frequency, a phenomenon that holds for stimulation frequencies up to 500 Hz 24 .Alternatively, voluntary commands are conceptualized by time-varying ring rates in the supraspinal Poisson inputs.Variability is introduced at every element of this circuit (see methods, Extended Data Fig. 1a-c).Finally, we modeled a lesioned system by modifying ion conductances of the spinal motoneurons membrane (see methods) and by reducing the number of available supraspinal bers, thereby simulating different degrees of lesion severities.
With this model, we calculated the changes in membrane potentials of spinal motoneurons caused by simultaneous supraspinal and SCS inputs.We then used it to explore the effects of SCS parameters on the ability of residual supraspinal inputs to produce ring rates in spinal motoneurons, i.e., the ability to recover voluntary motor control.

Supraspinal inputs transform subthreshold SCS-mediated EPSPs into action potentials
Using our model, we inspected the membrane potential of motoneurons before and after a lesion.Prior to a lesion, supraspinal inputs generated regular and sustained action potentials in the motoneuron membrane (Fig. 2a, top and middle panels).In our model, the slow calcium dynamics of spinal motoneurons expectedly limited their ring rate to 40 Hz 27 (Fig. 2a, bottom panel).We then represented a lesion by suppressing the activity of about 70% of supraspinal bers (Fig. 2b, white left).Incoming EPSPs produced a mean tonic depolarization in the membrane of 5.6 mV, ramping up within 100 ms after the start of the simulated supraspinal activity.This tonic depolarization effectively set the motoneuron membrane closer to ring threshold (around − 60 mV).However, it was insu cient to produce action potentials, thus leading to simulated paralysis.We then simulated a train of subthreshold SCS that alone (i.e., without supraspinal inputs) did not su ce to directly produce action potentials in motoneurons (Fig. 2b, blue middle).As opposed to the tonic membrane depolarization produced by residual natural inputs, subthreshold SCS induced a transient depolarization.This transient depolarization occurred through discrete EPSPs events in consequence of the synchronized recruitment of Ia afferents following each stimulation pulse.
Next, we combined subthreshold SCS with residual supraspinal inputs, which led to a robust production of action potentials that seemed aligned to the stimulation pulses (Fig. 2b, gray right).Indeed, when calculating the interspike interval (ISI) of the motoneuron pool, we found that it coincided with the stimulation interpulse interval (Fig. 2e).In other words, for each stimulation pulse, motoneurons produced an action potential.Subsequently, we increased the stimulation frequency to further explore this phenomenon and found that, for SCS frequencies above 19 Hz (still subthreshold), motoneurons were not able to generate action potentials for each SCS pulse because of the time constraint in their slow calcium dynamics (Fig. 2c,f; Extended Data Fig. 1d; Extended Data Fig. 2).Speci cally, motoneurons membrane skipped those SCS pulses that were elicited within the calcium dynamics following a motoneuron action potential (Fig. 2f light arrows in the inset).In turn, the number of skipped pulses increased as a function of the SCS frequency (Extended Data Fig. 3c-g).However, action potentials still occurred only in consequence of SCS pulses and, despite this decoupling above 19 Hz SCS, higher stimulation frequencies increased non-monotonically motoneuron pool ring rates (Fig. 2d).
In summary, we observed that tonic depolarization from residual supraspinal inputs transforms subthreshold SCS-mediated transient depolarizations into suprathreshold events in the motoneuron membrane.Yet, motoneurons were not able to follow these suprathreshold inputs for all frequencies, skipping some stimulation pulses before producing action potentials.This resulted in characteristic peaks in motoneuron ISI distributions at integers of SCS interpulse intervals.

Voluntary vs involuntary motoneuron activity during SCS
Given the in uence of SCS frequency on motoneuron ring rates, we extended our analysis to study the effects of SCS parameters on the production of volitional motoneuron activity.To this end, we quanti ed the ring activity in a simulated pool of motoneurons during SCS alone (Fig. 3a) and during concurrent SCS and residual input (Fig. 3b), while varying SCS amplitude (i.e., percentage of Ia afferents recruited, 10-100%) and frequency (8-256 Hz).We found that, for certain parameters (a nonlinear trade-off above 30-70% and 12-124 Hz), SCS alone directly produced sustained motoneuron ring rate (> 8 Hz) 28 , even without supraspinal input (Fig. 3a).This well-known phenomenon has been thoroughly described in animal and human studies as suprathreshold SCS 3,17,29 .This combination of parameters sharply differed from the subthreshold parameters (below 30-70% and 12-124 Hz, Fig. 3a), where SCS alone did not produce ring rates in the motoneuron pool, being effectively "sub-threshold" 8 .
When we combined SCS with residual supraspinal inputs, we observed a general increase in motoneuron ring rates (Fig. 3b).Interestingly, most of the subthreshold SCS parameters that led to sparse or no ring rate showed sustained production of motoneuron ring rates.This pattern presented two large "regimes'' in the SCS-parameter space that we termed voluntary and involuntary movement regimes (Fig. 3b; Extended Data Fig. 4).In the voluntary movement regime, SCS-mediated EPSPs became suprathreshold (only with concurrent inputs from residual supraspinal bers) and motoneurons showed sustained ring rates (8-20 Hz).Instead, in the involuntary movement regime, SCS dominated motoneuron ring rates (20-30 Hz), with little in uence from incoming residual supraspinal inputs.
In summary, we showed that there exists a well-de ned range of stimulation parameters in which residual supraspinal inputs transform subthreshold SCS pulses into suprathreshold events.

Residual supraspinal input can modulate motoneuron ring rate
The fact that residual supraspinal inputs can transform subthreshold SCS pulses into action potentials could imply that residual supraspinal inputs solely act as a gating mechanism.In this case, motoneuron ring rates would be fully controlled by SCS parameters and volitional modulation of forces would be impossible.Since this would contradict experimental evidence in humans, we explored the role of residual supraspinal inputs in the ring rates of motoneurons during SCS.We simulated motoneuron output during 45 Hz SCS for different strengths of residual supraspinal inputs by changing supraspinal mean ring.We observed that the motoneuron pool ring rates linearly increased with residual supraspinal bers ring rates during SCS amplitude within the voluntary regime (i.e., 731.7% modulation of the motoneuron ring rate at 40% SCS amplitude).Instead, in the involuntary regime, motoneuron ring rates were not affected by residual supraspinal inputs (i.e., 11.5% modulation at 80% SCS amplitude, Fig. 4a).We hypothesized that this modulation could occur by two mechanisms: 1) controlling the number of motoneurons that respond to SCS and 2) controlling the number of skipped pulses by recruited motoneurons.
We inspected simulated membrane behavior to test for the rst mechanism and found that, if residual supraspinal inputs was su ciently strong (i.e., higher residual supraspinal bers ring rates), it could induce stronger depolarization of motoneurons membrane, thereby bringing motoneurons closer to ring threshold and increasing the number of recruited motoneurons (see an example of a motoneuron that was recruited by SCS with high residual supraspinal inputs (dark blue), but not low (light blue), Fig. 4b).Extending these results over the motoneuron pool, the number of motoneurons that responded to SCS increased with higher residual supraspinal bers ring rates (i.e., recruitment probability increase of 82.7% from supraspinal bers ring at 5 Hz vs 60 Hz, Fig. 4d).This means that residual supraspinal inputs can actively control the number of motoneurons that respond to SCS.
Second, we found that, as the ring rate of the residual supraspinal bers increased, motoneuron depolarization was not only stronger but brought motoneurons to threshold faster.This increase in depolarization speed during continuous SCS pulses effectively decreased the number of SCS pulses skipped (i.e., SCS pulses that were not transformed into suprathreshold events).For instance, the motoneuron in Fig. 4c skipped 4 to 5 SCS pulses with residual supraspinal bers ring rate of 20 Hz, but it only skipped 3 with a residual supraspinal bers ring rate of 45 Hz.This means that the residual supraspinal inputs can also control the number of skipped SCS pulses for motoneurons that respond to SCS.Indeed, ISI histograms normalized by the SCS interpulse interval were shifted towards smaller values for stronger residual supraspinal inputs (Fig. 4e-g).
To illustrate the impact of SCS on ne motor control, we used a task where residual supraspinal inputs followed a sinusoidal function and we estimated the force produced by a pool of motoneurons 30 (see methods, Fig. 4h,i).We estimated the error by comparing the ring rate of the pool of motoneurons with the prelesion case.(Fig. 4j), thereby assessing SCS parameters e cacy.We found that SCS parameters in the voluntary movement regime assisted residual supraspinal inputs to follow the task (Fig. 4h,j).Instead, SCS parameters in the involuntary movement regime impeded residual supraspinal inputs modulation of motoneuron activity, yielding high task errors (Fig. 4i,j).
In summary, our biophysical model predicted that, while every action potential is triggered by an SCS pulse, residual supraspinal inputs can still volitionally modulate motoneuron activity by changing the number of motoneurons that respond to SCS and the number of SCS pulses that they skip effectively, affecting recruitment and ring rates of motoneurons.

Experimental validation in monkeys
We conducted animal experiments to validate the underlying model assumptions and predictions.We designed an experimental setup in anesthetized macaque monkeys (n = 2; Mk-JC, Mk-Ka), who share with humans monosynaptic corticospinal innervations of spinal motoneurons 31,32 .Brie y, to replicate our simulations in a controlled fashion, we applied SCS at C6-C8 cervical spinal segments, while manipulating supraspinal input by means of deep-brain stimulation of the internal capsule 33 , which is a deep-brain structure hosting corticospinal bers.In parallel, we measured motoneuron output directly, via intraspinal neural recordings, and indirectly, via evoked intramuscular electromyographic (EMG) in the hand muscles (Fig. 5a).
Our main model prediction was that supraspinal bers can modulate motoneuron membrane and reduce their threshold for incoming SCS-mediated EPSPs.If this is true, conditioning a pulse of SCS with a single pulse of internal capsule stimulation (ICS) should facilitate SCS and produce larger SCS-triggered spinal re exes in targeted muscles.Moreover, this facilitation should occur at latencies that are compatible with the time coincidence of monosynaptic Ia-afferent mediated EPSPs and cortex-mediated excitatory inputs 21 .
We thus designed a stimulation protocol to condition SCS with a pulse of ICS at different delays (5-50 ms, Fig. 5b).In both monkeys, mean peak-to-peak of SCS-induced re ex amplitudes signi cantly increased when SCS was conditioned by a pulse of ICS, demonstrating that corticospinal inputs can facilitate SCS-induced re exes (Mk-JC: 111.3%, 17.2%, 18.9%, 7.2% in APB; 285.8%, 590.9%, 793.3%, 557.7% in FDM.Mk-Ka: 19.6%, -5.12%, -1.8%, -18.3% in FDM; 44.6%, 8.5%, 22.4%, -7.7% in APB, Fig. 5c,d).Speci cally, the greatest re ex peak-to-peak amplitude was obtained at 5 ms delay in Mk-Ka, when both EPSPs coincided in the hand muscles (Extended Data Fig. 5a).Another important feature of our model is that both Ia afferents and supraspinal bers excite motoneurons through fast glutamatergic synapses, resulting in a similar effect on the motoneuron membrane.Thus, SCS-induced EPSPs should modulate corticospinal motor-evoked potentials.Therefore, we repeated the same experiments but now conditioning ICS with a pulse of SCS at different delays.Indeed, corticospinal motor-evoked potentials mean peak-to-peak amplitudes signi cantly increased in comparison to control (i.e., SCS alone) in both monkeys, in contrast to previous reported results 21 (Fig. 5e,f).These results show that both inputs can depolarize membrane potentials of motoneurons and facilitate a subsequent EPSP.
Second, we validated the assumption that SCS recruits sensory afferents eliciting synchronous volleys at a given stimulation frequency (Fig. 1b).To do so, we inspected postsynaptic motoneuron ring when receiving SCS-elicited sensory inputs alone.As shown in our simulations (Extended Data Fig. 3b), if SCS produces synchronous volleys in the recruited afferents, motoneurons should produce spikes aligned to these sensory inputs.We delivered SCS at a high frequency (100 Hz) in both monkeys in the same anesthetized procedure, thereby isolating motoneuron ring from alternative sources of modulation other than SCS.We then identi ed intraspinal single motor units from the ventral channels in the linear probe (Fig. 5g), corresponding to the location of motoneuron pools in lamina IX 34,35 .We computed ISI distributions and detected spectral peaks at exact integers of SCS pulses intervals (Fig. 5h).In both monkeys, motoneuron action potentials occurred near SCS pulses.
Third, we sought to con rm the link between SCS parameters and motoneuron ring modulation.Simulations revealed that the number of skipped SCS pulses increased with the SCS frequency (Extended Data Fig. 3c-g).Therefore, we computed the ISI cumulative distributions, which re ected this pattern (Fig. 5j).Indeed, for low SCS frequencies (30 Hz), sorted motor units skipped less SCS pulses in both monkeys and, in turn, the cumulative distribution of the ISI increased faster for less stimulation pulses in comparison to 100 Hz.We expanded this analysis for medium frequencies (60 Hz) and found consistent results in Mk-Ka and less pronounced in Mk-JC.Lastly, we computed the overall ring rate for the pool of sorted motor units and found that maximum ring activity was obtained at higher stimulation frequencies (i.e, 100 Hz in Mk-Ka and 60 Hz in Mk-JC, Fig. 5i).Hence, as predicted by the model (Fig. 2d), despite pulse-skipping at these high stimulation frequencies, motoneurons ring increased.
These experiments validated some of the model assumptions and results.However, being the monkeys anesthetized, we could not demonstrate our central results concerning combined volitional input and SCS and, speci cally, the transformation of subthreshold SCS-induced EPSP into suprathreshold events by the residual supraspinal inputs.

Volitional control of spinal motoneurons in humans during SCS
To verify our main ndings concerning concurrent natural activity from supraspinal bers and SCS, we conducted experiments in humans that received spinal cord implants.Speci cally, we leveraged our existing Institutional Review Board protocols in two participants with poststroke arm hemiparesis (SCS03, SCS04), who participated in a pilot clinical trial exploring the effects of cervical SCS on motor control of the hand and arm 6 (Extended Data Fig. 6a,b), and an observational study in a participant with severe spinal cord injury (STIM01, ASIA A), who received an SCS implant on the lumbosacral spinal cord to reduce neuropathic pain in the lower limbs (Extended Data Fig. 6c,d).By testing our hypothesis in humans with stroke and spinal cord injury (SCI), we aimed at demonstrating that there are fundamental principles underlying the mechanisms of SCS irrespective of the speci c disease that caused the lesion to the descending pathways.
We started by testing the effects of subthreshold SCS on motoneuron ring rates in STIM01, who had complete leg paralysis in consequence of a SCI; hence, a severe lesion.We asked the participant to contract his paralyzed leg in an isometric task at the maximum voluntary strength with and without SCS.
Concurrently, we extracted motor units using surface EMG 36 (Fig. 6a,b; Extended Data Fig. 7a).While the participant could not produce any force without SCS, SCS led to substantial force production (Fig. 6c; Extended Data Fig. 8a).However, the stimulation amplitude required to enable this voluntary production was just within the suprathreshold region.Indeed, we found signi cant motor responses even in moments in which the leg was supposed to be at rest (Extended Data Fig. 8a, brown area).When we inspected the ISI distribution of motoneurons during the isometric contraction task, we found clear peaks aligned to SCS pulses (Fig. 6d; Extended Data Fig. 8a) that closely resembled our simulations (Fig. 2e,f), which could not be explained by stimulation artifacts (Extended Data Fig. 7d).
We hypothesized that similar mechanisms should also apply for the control of upper limbs, which are known to be more dependent on direct corticospinal input 6 .Thus, we tested our predictions in two participants with stroke (SCS03, SCS04) that received a cervical spinal cord implant.Our participants with stroke retained the ability to produce force also without SCS representing a milder lesion severity.We exploited this capability to test the validity of our model assumption that describes the ring of supraspinal bers as stochastic (Fig. 1b) in natural conditions.As depicted by our simulations (Extended Data Fig. 3a), if supraspinal bers re stochastically, motoneuron ring should show non-synchronous ring patterns during volitional movement alone 37,38 from the computed ISI distributions from single motor units from the paretic arm muscle (Extended Data Fig. 7b,c,e).Without SCS, as expected, the ISI distributions did not depict peaks of activity across time (plotted as a function of the number of SCS pulses for comparison purposes, Fig. 6f).Instead, during stimulation, once again the ISI distribution matched our ndings with our simulations and in the SCI participant (Fig. 6f; Extended Data Fig. 8b).We could not test the involuntary SCS parameter regime as high stimulation amplitudes induced pain in these participants.
Lastly, we tested our model prediction that residual supraspinal inputs can modulate motoneuron ring rate and, thus, force by controlling the number of SCS pulses skipped.Thus, we sought to measure the motor unit ring for different force levels, which correspond to different levels of supraspinal input strength (Fig. 6e,g).In particular, we parameterized the lower and upper force boundaries of our forcemodulated task, which represent a percentage of the maximum force, and asked the participants to follow a target trace between 5-20%, 20-35% and 35-50% of their maximum force (Fig. 6g).The ISI cumulative distributions obtained from this task showed that both SCS03 and SCS04 skipped less SCS pulses when producing high forces, in contrast to low forces (Fig. 6h; Extended Data Fig. 8c).
In summary, our results in humans show that motor unit ring is aligned to stimulation pulses, regardless of motor diseases and lesion severity.We also demonstrated that, when present, residual supraspinal inputs can modulate motor unit ring induced by SCS to grade force production by reducing pulseskipping.

Subthreshold SCS e cacy may depend on lesion severity
In our human data we observed that, while stroke participants retained the ability to largely modulate SCS output, our SCI participant required high-stimulation amplitudes to enable movement that were already slightly suprathreshold.Therefore, we sought to exploit our biophysical model to study the theoretical impact that lesion severity could have on the range of effective parameters of SCS.We replicated in our model an isometric contraction task for a range of lesion severities (73-96%) and found that increasing the lesion severity to the supraspinal bers reduced the combinations of subthreshold SCS parameters that generated sustained motoneuron ring (> 8 Hz, Fig. 7a) with concurrent activation of the residual supraspinal inputs.In turn, this pattern was translated into the shrinking of the voluntary movement regime to smaller sizes (from 87 to 15 combinations of SCS parameters, Fig. 7b,c; Extended Data Fig. 9a).
We tested these results again in monkeys Mk-JC and Mk-Ka, where we could causally control the amount of inputs from the corticospinal tract via ICS.We superimposed subthreshold continuous SCS to bursts of ICS (Fig. 7d) and found that EMG responses produced during ICS bursts were potentiated as a function of the SCS amplitude.We then repeated the same stimulation protocol but at a very low ICS amplitude, thereby mimicking a case with very limited residual corticospinal inputs (i.e., a severe lesion).In this case, the root-mean square of the EMG responses decreased but it could be restored at higher SCS amplitudes, as in the case of our human patient with SCI (e.g., in Mk-JC, 8.4 mV during 140 uA SCS at 0.7 mA ICS vs 8.8 mV during 300 uA at 0.5 mA ICS, Fig. 7e; Extended Data Fig. 10), which was also clear in the mean torque produced by the hand grip (Fig. 7d,f).
Together, these results show that the useful parameter space at which SCS can facilitate production of volitional muscle activity shrinks with the reduction of the residual supraspinal inputs, thereby suggesting that the e cacy of SCS in facilitating motor control is inversely proportional to the severity of the damage to supraspinal pathways.

DISCUSSION
In this study, we combined biophysical modeling with animal and human experiments to study the neural mechanisms underlying the recovery of voluntary control of spinal motoneurons during SCS after paralysis.We found that, after a lesion, residual supraspinal inputs can tonically depolarize spinal motoneurons bringing transient EPSPs induced by subthreshold SCS into action potentials.This transformation causes motoneuron ring to be aligned to SCS pulses, but supraspinal input can still in uence motoneuron activity by controlling the number of active motoneurons and skipped SCS pulses.
However, supraspinal control of motoneuron ring only occurs for a restricted range of subthreshold SCS parameters and the space of effective parameters shrinks for more severe supraspinal axon loss.Our results imply that the e cacy of SCS depends on the amount of spared supraspinal input.

On the validity of the concept of neuromodulation
It is a common tendency to de ne the enabling effect that we describe here as "neuromodulation".However, our ndings exclude the possibility that SCS, delivered at frequencies below 250 Hz, is modulating the membrane of neurons in the strict sense of the term "modulation", i.e., by biasing membrane excitability.In fact, it is the residual natural activity that can exert this biasing action by means of fast neurotransmitters such as glutamate, but also by the slower serotonergic and noradrenergic systems 39 .This natural biasing of cell excitability will then transform subthreshold, discrete SCS events into suprathreshold events (Fig. 2).While from behavioral observations, both interpretations are indistinguishable (since movement can only occur when intention is present), the speci c mechanism does imply important differences on motoneurons behavior.For example, our model explains why stimulation frequency is related to proportional changes in strength and kinematics observed in our experiments (Fig. 6) and others 3,6,29,40,41 : SCS determines motoneuron spike times, forcing them to re coherently with stimulation pulses.However, at the same time, residual supraspinal inputs can control the number of SCS pulses transformed into action potentials, thereby modulating motoneuron ring rates and producing controllable forces.This is in contrast to other neurotechnologies that directly drive motoneurons 42,43 , which would instead produce highly synchronized motoneuron ring and, in turn, result in poorly controllable forces.In other words, this is one of those cases in which the simple description of behavior does not offer satisfactory explanatory power of a phenomenon that requires analysis of neural output to be correctly interpreted 44 .In summary, while our results did con rm the "enabling" effects of SCS, they also suggest that the term "neuromodulation" should be used with caution, particularly, when studying the effects at the postsynaptic level, and suggest that similar effects may occur with other neurostimulation technologies like SCS for pain or deep-brain stimulation (DBS).

Supraspinal control of movement after paralysis
Through the analysis of the membrane mechanisms of the response of spinal motoneurons to SCS pulses, we studied the computational principles that govern the integration of SCS mediated inputs to motoneurons and residual supraspinal activity.However, beyond spinal motoneurons, these mechanisms are likely to hold also for spinal interneurons that receive signi cant inputs from Ia afferents activated by SCS 45,46 .For example, central pattern generators 47,48 , inhibitory interneurons 49 and excitatory premotor interneurons involved in muscle synergies 50 , all receive strong sensory inputs from proprioceptive afferents.Therefore, similarly to what we observed in motoneurons, residual supraspinal modulation can bring subthreshold SCS events into suprathreshold action potentials also in these network elements.In other words, the divergent connectivity of Ia afferents towards large parts of the sensorimotor network guarantees the brain access to many crucial components.This interpretation could explain why, despite the well-known low speci city of SCS 11,17 , we can observe repeatable, strong clinical effects 4,5 .In fact, we can speculate that a more selective technology that only enables control of motoneurons, or other speci c sets of neurons, would fail at restoring motor control.The success of SCS does not stem from the activation of speci c neural elements, but rather from the widespread occurrence of the membrane mechanism that we identi ed here in many spinal network nodes via the projection of activated sensory afferents, thus enabling complex, volitional motor control.

Clinical implications
One of our most critical ndings is that the well-de ned SCS parameter regime that enables volitional control of spinal motoneurons shrinks as the number of residual supraspinal inputs decreases, i.e., as the severity of the lesion increases (Fig. 7).In other words, in a mild lesion (say 70% damage), where su cient residual supraspinal bers survived, the space of effective SCS parameters can be more than 7 times larger than in a severe lesion (say 95% damage).This is caused by the fact that, to compensate for severe loss of excitatory drive and to produce su cient ring to enable movement, SCS amplitude must be set to values that rapidly enter the "suprathreshold" region, which we termed "involuntary movement" regime, meaning that activation and, likely, muscle co-contraction, will occur at rest.This intrinsic limitation of SCS emerges when inspecting aggregated experimental human data 51 that show a clear trade-off between e cacy and lesion severity in SCI.In this scenario, the residual ability of people with paralysis to inhibit unwanted SCS inputs, or the possibility to acquire it with physical training, becomes more important than their capacity to excite spinal circuits (Extended Data Fig. 9b).In fact, it is conceivable that a person with severe paralysis could learn to sculpt suprathreshold SCS inputs either by presynaptic inhibition of the Ia afferents or postsynaptic inhibition of the motoneurons, thereby suppressing unwanted activity and directing SCS excitatory drive towards movement-relevant neural pathways.It is possible that this is one of the effects of physical training, which, by strengthening synaptic connectivity, can enable voluntary control of SCS 2,8,18 .Finally, technology can also improve the e cacy of SCS, for example, closed-loop, spatiotemporal modulation of SCS 3,29,41,52 can help support the sculpting of SCS and facilitate people with severe injuries to produce complex movements (Extended Data Fig. 9c).While these alternatives can help increase the applicability of SCS, it is important to highlight that intrinsic limitations will eventually reduce SCS applicability and effectiveness for people with severe paralysis irrespective of the speci c condition causing paralysis.It is conceivable that these limitations will be even more relevant for the recovery of upper limb movements given the strong dependence of hand control on corticospinal drive.

Limitations
In this study, we only considered monosynaptic connections to spinal motoneurons coming from supraspinal axons as well as Ia afferents.Obviously, many more components that play a role in regulating excitability in the cord are stimulated by SCS, for example, large cutaneous mechanoreceptors, Ib afferents and group II afferents 10,11 .Nevertheless, as we discussed above, the basic membrane mechanism that regulates the changes in the membrane are likely to hold when considering these other inputs and, importantly, for other types of neurons that receive inputs from these afferents, such as Iainhibitory interneurons, central pattern generators and excitatory premotor interneurons implementing muscle synergies.Finally, our model only focuses on immediate effects of SCS and cannot be used to directly predict long-term changes that occur after physical training.Future work should be directed towards modeling long-term circuit changes that represent a major driver of SCS clinical e cacy.
In conclusion, we believe that detailed understanding of the neural mechanisms that drive the recovery of voluntary motor control will enable the optimization of SCS and other neurotechnologies to maximize e cacy of the residual brain-spinal architecture and control movement after paralysis, thereby promoting true motor recovery.
We implemented a biophysical model of motoneurons, Ia afferents and supraspinal bers in Python 3.8 using NEURON 7.8 53 .The following number of neurons and their properties have been validated in similar models 10,11,24,25,27 .
Motoneuron model.We simulated a pool of 169 motoneurons.Each motoneuron comprised an integrated multicompartmental soma connected to an electronic-equivalent dendritic tree and myelinated axon, through an initial axon segment.The established approach to describe the dynamics of the motoneuron membrane potential as a function of the ion channels is the Hodgkin-Huxley model.Eq. 1 describes the membrane potential (V): where C is the membrane capacity, is current arti cially injected to the neuron and is the conductance of the ion , which depends on the membrane potential and ion concentrations.Thus, we implemented dedicated Hodgkin-Huxley equations for multiple ion channels (L-type Ca2+, Ca2+-activated K, delayed recti er K+, N-type Ca2+, and nonlinear fast Na+) to study the integration of postsynaptic potentials from multiple presynaptic sources of modulation in the motoneuron membrane.Provided the increase in excitability in the motoneuron membrane after a lesion, we de ned a model of postlesion motoneuron with modi ed ionic properties.Speci cally, we reduced to 60% the Ca2+-activated K ion channel conductance 54 .We assumed that the motor cortex is unaltered after a lesion.
Excitatory synapses.Each motoneuron received excitatory inputs stochastically distributed along the dendritic tree from 60 Ia afferents, which were synchronously activated by each SCS pulse, and supraspinal bers.The number of supraspinal bers was de ned by their synaptic weight and the number of pulses per second to each motoneuron 55 .In contrast to motoneurons, both EPSPs were modeled as point processes characterized by their conduction velocity and muscle ber length.These synapses followed an exponential function with a reversal potential of = 0.5 mV and a decay time constant = 0.5 ms.
Both EPSPs represented fast, glutamatergic synapses with similar dynamics to drive motoneuron excitation.Their synaptic weight followed a normal distribution (mean 0.047, standard deviation 0.0094) 24 .When supraspinal bers point processes were activated, an EPSP was delivered in the motoneuron dendritic tree.To account for the variability in supraspinal ber diameters, synaptic transmission delays followed a normal distribution (mean 0.25 ms, standard deviation 0.05 ms) plus a period of 1 ms.Instead, Ia afferents travel times were modeled as a latency between SCS pulses and elicited EPSPs that followed a log-normally distributed stochastic jitter (mean 1.0 ms, standard deviation 0.5 ms) 11 , thereby considering the difference in afferent diameter.Ia-afferent recruitment occurred when SCS pulses arrived outside the membrane refractory period.EPSP activation dynamics from both sources of modulation was modelled equally and followed a normal distribution (mean 1.6 ms, standard deviation 0.16 ms).
Force model.To estimate the force from the motoneuron pool ring rate, we used the model proposed in 30 .Here, we described the main equations (Equations 2-5).The twitch force generated for a spike from neuron at time is modeled as:
The peak twitch is de ned as and the contraction time is .Finally, the gain, was modeled depending on the normalized mean interspike interval of neuron ( and the contraction time : where is the number of motoneurons and is the number of spikes in neuron .We always reported the force normalized with the maximum prelesion force. ) .
Experimental procedures in monkeys Animals.All procedures were approved by the University of Pittsburgh Animal Research Protections and Institutional Animal Care and Use Committee (ISOOO17081).We conducted experiments in 2 anesthetized Macaca Fascicularis (Mk-JC: age 6, 7.5 kg, male; Mk-Ka: age 5, 6.8 kg, female).The animals were housed in the primate facility at the Division of laboratory Animal Resources at the University of Pittsburgh.The animals had unrestricted access to water and food as well as daily enrichments.
Surgical procedures.Animals underwent two surgical procedures.First, in a survival procedure, we acquired in vivo structural magnetic-resonance imaging and computed tomography for each animal.
Collected imaging data was used to plan the trajectories of the deep-brain probe that stimulated the posterior limb of the internal capsule, which innervates the hand muscles.We further describe this procedure in 56 .Second, in a terminal procedure, certi ed neurosurgeons performed muscle implantation, robotic deep probe implantation and spinal probe implantation (described below in detail) while the monkeys had continuous ventilation support and were under full anesthesia induced with ketamine (10 mg/kg, i.m.) and under continuous intravenous infusion of propofol (1.8-5.4 ml/kg/h) and fentanyl (0.2-1.7 ml/kg/h).Anesthesia conditions were maintained until the conclusion of the experiments, where animals were euthanized with a single injection of pentobarbital (86 mg/kg, i.v.).
We inserted a pair of EMG needle electrodes (disposable single 7mm Subdermal needle electrode, Rhythmlink) into the hand muscles (abductor pollicis brevis (APB), exor digiti minimi (FDM)).The reference was inserted in the lower back or the tail with a needle electrode.Subsequently, we xed the monkey's head in a stereotactic frame (Kopf, Model 1530, Tujunga, CA, USA) with the cervical spine xed in a prone position.
We proceeded with the robotic deep probe implantation (ROSA robot) 56 .To accurately implant a probe in the hand area of the internal capsule (AC-PC level), a frontal insertion trajectory was guided by the robot software (ROSA One(R) Robot Assistance Platform).Trajectory planning was achieved due to previously collected imaging data 56 .The probe positioning was contralateral to the muscle implants.We precisely drilled a penetration hole in the skull for the implantation of a xation bolt.Next we inserted a radiofrequency electrode (RF-SE-10 Reusable Stainless Steel Electrode, Abbott).We con rmed the correct position of the probe by inspecting the hand muscles EMG activity in response to ICS (2 Hz, 0.8-2 mA) that should produce monosynaptic activation of the cervical motoneurons.
Lastly, we performed a laminectomy from C3 to T1 vertebrae and exposed the cervical spinal cord to the nerves that innervated targeted muscles.After opening the dura matter, we placed a stimulator (bipolar, 73602-190/10, Ambu Neuroline, Denmark) on the cervical dorsal roots (Mk-JC: C8, Mk-Ka: C6).We con rmed its positioning by verifying EMG responses in the target muscles.Furthermore, we implanted a 32-channel linear probe (A1x32-15mm-100-177-CM32 Linear Probe with 32 pin Omnetics Connector, NeuroNexus, Ann Arbor, MI, USA) at the cervical spinal segment rostral to the stimulation site (Mk-JC: C7, Mk-Ka: C7) using Kopf micromanipulators (Kopf, Model 1760, Tujunga, CA, USA).To do so, we pierced the pia using a surgical needle to ensure safety probe penetration.Its placement was ~ 1.2 mm parallel to midline (ipsilateral to the muscle implants) due to the ow of sensory information of Ia afferent in the dorsoventral direction consistent with the placement of the probe.Its length (3.1 mm) is designed to cover the gray matter of the monkey's spinal cord and record extracellular potentials as well as intraspinal spiking activity.The reference was inserted under the skin next to the probe insertion site.
Grasp force data acquisition.We placed a force transducer in the anesthetized monkey's hand (ipsilateral to the muscle implants).Grasp force data was captured with a 6-axis low-pro le force and torque (F/T) sensor (Mini40, ATI Industrial Automation, North Carolina).The sensor was powered using a multi-axis (F/T) transducer system (ATI DAQ F/T) that also calibrated the force data.The sensor measured three degrees of force (+/-810 N for Fxy, +/-2400 N Fz) and torque (+/-19Nm Txy, +/-20Nm Tz) simultaneously.A small rod was mounted to the xy-plane of the sensor to be a grip handle for the animal.
Data acquisition and electrophysiology.We employed an AM stimulator (model 2100 A-M Systems, Sequim, WA, USA) for the stimulation of the dorsal roots and internal capsule as biphasic symmetric square pulses in the form of single or train of pulses and detailed pulse widths (Mk-JC: 0.40 ms for SCS, 1 ms for ICS; Mk-Ka: 0.45 ms for SCS, 0.8 ms for ICS).Speci c stimulation con gurations and parameters are de ned in each gure.After data ampli cation (Nano2-HV headstage) and digital processing, neural recordings (EMG and intraspinal spiking data) were recorded (Ripple Neuro Grapevine Neural Interface 1.14.3) at 30 kHz sampling frequency along with the force/torque data (same sampling frequency).Through the data acquisition platform (Trellis), we developed an in-house interface that computed the EMG stimulation-trigger averaged responses in real-time to determine the stimulation amplitude for motor threshold responses for each muscle.Moreover, intraspinal recordings were treated with a broadband lter between 300 Hz and 3 kHz and a voltage threshold (3 root-mean square) to determine neural events.
Analysis of muscle activity.EMG signal was calculated for each muscle differentially between the two needle electrodes and was bandpass ltered between 30 and 800 Hz with a 2nd order Butterworth lter.
We then identi ed the peak-to-peak amplitude of re ex-mediated responses evoked by the paired stimulation in a window of 25 ms after the SCS pulse.Moreover, we computed the root-mean square of the response during phasic ICS alone and concurrent phasic ICS and continuous SCS.This analysis was performed o ine in MATLAB (version R2019b).
Analysis of single unit activity.We sorted spikes from the 8 most ventral channels of the 32-channel linear probe.Previously, we applied a comb lter of the raw signal to remove artifacts at 60 Hz.We then applied a 3rd order Butterworth digital lter (800-5000 Hz).Spikes were detected by waveform threshold crossing.All intraspinal recordings were concatenated by channel for comparison purposes across performed experiments.Using a spike sorting software (Plexon O ine Sorter, Dallas, Texas), we identi ed single motor units by channel with semi-automatic k-means clustering.Regarding the histograms of the normalized ISI probability distribution of the sorted units, we used a bin width of 0.1 ms.Previously, we noticed that the application of continuous SCS resulted in an initial short hyperexcitable period, thus, we discarded the rst 5 seconds of spiking data after the start of the stimulation to perform our analysis on stationary conditions.This analysis was performed o ine in MATLAB (version R2019b).
Analysis of torque activity.Mean torque signal was computed as the module of all three measured directions and ltered by applying a 15 kHz low bandpass lter with a 3rd order Butterworth lter.
Provided that monkeys were anesthetized, we obtained the baseline value by subtracting the mean torque response during the no stimulation period during ICS alone for each amplitude (0.7 mA and 0.5 mA).This analysis was performed o ine in MATLAB (version R2019b).
Experimental procedures in humans Trials information.All experimental protocols were approved by the University of Pittsburgh Institutional Review Board (IRB): STUDY19090210, for the stroke participants (SCS03 and SCS04), and STUDY22020031, for the spinal cord injury participant (STIM01).All participants provided informed consent according to the procedure approved by the IRB of the University of Pittsburgh and participants were compensated for each day of the trial and for travel and lodging during the study period.
The stroke participants (SCS03 and SCS04) were recruited to participate in an exploratory clinical trial to obtain preliminary evidence of safety and e cacy of SCS to improve motor control in people with chronic post-stroke upper-limb hemiparesis 6 .A detailed description of the trial, including inclusion criteria and study design, has been described in 6 and can be found on ClinicalTrials.gov(NCT04512690).The spinal cord injury participant (STIM01) was recruited to participate in a basic research study enrolling people already living with an epidural spinal cord stimulator (STUDY22020031).During the experimental sessions in STIM01, we changed the stimulation parameters to perform the different tasks using the Medtronic Intellis Platform.At the end of each session, we reset the stimulation parameters to their clinical values.
Participants Information.STIM01 (Black male, 24 years old) was diagnosed as ASIA A after a spinal cord injury at T6 vertebra level 5 years ago.He was implanted with a spinal cord stimulation to reduce pain and spasticity 4 years ago.Since the implantation, he has combined SCS and physical therapy to regain control of his lower limbs.When he enrolled in the study, he was able to voluntarily move his left leg during SCS.STIM01 was implanted with a Medtronic paddle array from the bottom of T11 to the top L1 vertebrae.SCS03 (White male, 45 years old) suffered a hemorrhagic stroke followed by craniotomy 7 years before participating in the study and was diagnosed with a Fugl-Meyer motor score of 28.SCS04 (Black male, 43 years old) suffered an ischemic stroke in the anterior frontal and temporal lobe 5 years before participating in the study and was diagnosed with a Fugl-Meyer motor score of 19.To minimize risks in SCS03 and SCS04, participants were temporally implanted for 29 days, after which electrodes were explanted.
Torque experiments.We performed the torque experiments in all subjects using a robotic torque dynamometer (HUMAC NORM, CSMi) that can be set to measure torque in different joints.These experiments and data acquisition are described in detail in 6 .
For STIM01, we measured his maximum torque during knee extension for 6 repetitions of 2 seconds (Extended Data Fig. 8a).We also asked STIM01 to hold a torque level of 35% of his maximum torque during 30 seconds (Fig. 6c).In both tasks, we placed a 4-channel EMG sensor (Trigno Galileo Sensor, Delsys) to record EMG activity on the Vastus Lateralis and Rectus Femoris.
SCS03 performed an isometric force-modulated task for elbow exion where he was asked to follow a prede ned force trace ranging from 5-20%, 20-35%, and 35-50% of his maximum voluntary contraction force with and without SCS (Fig. 6g; Extended Data Fig. 8c).This experiment was repeated with and without SCS.A 64-channel high-density EMG grid (TMSi, The Netherlands) was placed over Bicep Brachii to record muscle activation during the experiment.
Because of SCS04's inability to perform the force-modulated task at the elbow, he performed a similar task of grasp force-modulated task at different force levels, ranging from 5-20%, 20-35%, and 35-50% of his maximum torque, both with and without SCS (Fig. 6g; Extended Data Fig. 8c).The grasp force data was measured with a 6-axis low-pro le force and torque (F/T) sensor (Mini40, ATI Industrial Automation, North Carolina).A 32-channel high-density EMG grid (TMSi, The Netherlands) was placed over the anterior forearm (speci cally, Flexor Carpi Radialis) to record muscle activity during the task.
Motor unit decomposition.In STIM01, we used the Neuromap software (Delsys) to decompose the EMG activity in single motoneuron action potentials.For the experiments with SCS03 and SCS04, single motor unit ring rates were identi ed from high-density signals by the spike trains of motor units, which were identi ed by the Convolution Kernel Compensation (CKC) technique 57 .The high-density EMG signals were rst bandpass ltered between 20 and 500 Hz with a 3rd order Butterworth lter.Additionally, SCS artifact was removed by notch ltering at the SCS frequency.SCS harmonics were also removed using notch ltering by identifying peaks in the power spectrum associated with SCS.The ltered high-density EMG signals were decomposed by the DEMUSE software (v6.1;The University of Maribor, Slovenia) 58, 59 .
The decomposed motor unit spike trains were processed to identify single motor units and calculate their ring rates.During this step, physiologically irregular motor unit rings (<5 Hz and >50 Hz) were discarded 60,61 .The irregular motor unit rings were identi ed by following the guidelines in literature and accuracy of identi ed motor units was determined by the Pulse-to-Noise Ratio (>25 dB) 62,63 (Extended Data Fig. 7b,c,e).
Analysis of torque activity.To account for the weight of the paretic limb in STIM01, torque responses were baseline corrected with the mean resting torque for the rst 60 ms.Moreover, given the uctuations in the responses, we smoothed the data applying a moving window of 100 ms.Lastly, we computed the maximum torque during the maximum voluntary contraction periods for each repetition.Only for SCS04, given that the force-modulated task was performed with the same data acquisition system as in the monkeys, torque responses were analogously treated by applying the same low bandpass lter and baseline correction with the mean resting torque for the rst 3.33 seconds.This analysis was performed o ine in MATLAB (version R2019b).
Analysis of ring rate activity.Only for SCS03 and SCS04, we computed the ring rate in a 50 ms window with a sliding window of 5 ms over the entire duration of the trial and zeropadded for 45 ms and smoothed it by applying a moving window of 100 ms.We identi ed ring rates for each unit during a constant force level, i.e., a constant supraspinal innervation (replicating what we did in our biophysical model).Like in the monkey normalized ISI probability distributions, we used a bin width of 0.1 ms.As indicated in the captions, data included in the gures correspond to all the repetitions under the same conditions and session.This analysis was performed o ine in MATLAB (version R2019b).

Statistical
All statistical comparisons were performed using the Kruskal-Wallis test, which is a nonparametric approach to the one-way ANOVA.Subsequently, followup multiple comparison tests on pairs of sample medians were performed.Data points outside 1.5 times the interquartile range plus the 25th/75th percentile were classi ed as outliers and discarded.For illustration purposes, the mean distribution of the data with 95% con dence interval was indicated for each distribution.This analysis was performed o ine in MATLAB (version R2019b). Figures

Supplementary Files
This is a list supplementary les associated with this preprint.Click to download. ExtendedDataFigures.docx i