Biophysical modeling
We implemented a biophysical model of motoneurons, Ia afferents and supraspinal fibers in Python 3.8 using NEURON 7.853. The following number of neurons and their properties have been validated in similar models10,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):
$$C\frac{dV}{dt}={\sum }_{x}^{ions}{g}_{x}\left({V-V}_{x}\right) + I$$
1
where C is the membrane capacity, \(I\) is current artificially injected to the neuron and \({g}_{x}\) is the conductance of the ion \(x\), 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 rectifier 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 defined a model of postlesion motoneuron with modified ionic properties. Specifically, we reduced to 60% the Ca2+-activated K ion channel conductance54. 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 fibers. The number of supraspinal fibers was defined by their synaptic weight and the number of pulses per second to each motoneuron55. In contrast to motoneurons, both EPSPs were modeled as point processes characterized by their conduction velocity and muscle fiber length. These synapses followed an exponential function with a reversal potential of \({E}_{syn}\) = 0.5 mV and a decay time constant \(\tau\) = 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 fibers point processes were activated, an EPSP was delivered in the motoneuron dendritic tree. To account for the variability in supraspinal fiber 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 firing rate, we used the model proposed in30. Here, we described the main equations (Equations 2–5). The twitch force generated for a spike from neuron \(i\) at time \(t\) is modeled as:
$${f}_{i}\left(t\right)={g}_{i}\frac{{P}_{i}t}{{T}_{i}}exp(1-t/{T}_{i})$$
2
.
The peak twitch is defined as \({P}_{i}=exp(b\cdot i)\) and the contraction time is \({T}_{i}=90(1/{P}_{i}{)}^{1/4.2}\). Finally, the gain, \({g}_{i}\) was modeled depending on the normalized mean interspike interval of neuron \(i\) (\(IS{I}_{i})\) and the contraction time \({T}_{i}\):
For \({T}_{i}/IS{I}_{i}<0.4\):
,
otherwise:
$${g}_{i}=\frac{S({T}_{i}/IS{I}_{i})}{{T}_{i}/IS{I}_{i}}$$
4
,
where \(S\left(x\right)\) is the sigmoidal function\(S\left(x\right)=1-exp\left(2{x}^{3}\right).\)
Finally, the total normalized force is:
$$F\left(t\right)={\sum }_{i=1}^{N}{\sum }_{j=1}^{{N}_{i}}{f}_{i,j}(t-{t}_{i,j})$$
5
,
where \(N\) is the number of motoneurons and \({N}_{i}\) is the number of spikes in neuron \(i\). 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 in56. Second, in a terminal procedure, certified 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), flexor digiti minimi (FDM)). The reference was inserted in the lower back or the tail with a needle electrode. Subsequently, we fixed the monkey’s head in a stereotactic frame (Kopf, Model 1530, Tujunga, CA, USA) with the cervical spine fixed 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 data56. The probe positioning was contralateral to the muscle implants. We precisely drilled a penetration hole in the skull for the implantation of a fixation bolt. Next we inserted a radiofrequency electrode (RF-SE-10 Reusable Stainless Steel Electrode, Abbott). We confirmed 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 confirmed 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 flow 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-profile 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). Specific stimulation configurations and parameters are defined in each figure. After data amplification (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 filter 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 filtered between 30 and 800 Hz with a 2nd order Butterworth filter. We then identified the peak-to-peak amplitude of reflex-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 offline 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 filter of the raw signal to remove artifacts at 60 Hz. We then applied a 3rd order Butterworth digital filter (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 Offline Sorter, Dallas, Texas), we identified 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 first 5 seconds of spiking data after the start of the stimulation to perform our analysis on stationary conditions. This analysis was performed offline in MATLAB (version R2019b).
Analysis of torque activity. Mean torque signal was computed as the module of all three measured directions and filtered by applying a 15 kHz low bandpass filter with a 3rd order Butterworth filter. 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 offline 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 efficacy of SCS to improve motor control in people with chronic post-stroke upper-limb hemiparesis6. A detailed description of the trial, including inclusion criteria and study design, has been described in6 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 in6.
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 flexion where he was asked to follow a predefined 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-profile 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 (specifically, 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 firing rates were identified from high-density signals by the spike trains of motor units, which were identified by the Convolution Kernel Compensation (CKC) technique57. The high-density EMG signals were first bandpass filtered between 20 and 500 Hz with a 3rd order Butterworth filter. Additionally, SCS artifact was removed by notch filtering at the SCS frequency. SCS harmonics were also removed using notch filtering by identifying peaks in the power spectrum associated with SCS. The filtered 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 firing rates. During this step, physiologically irregular motor unit firings (<5 Hz and >50 Hz) were discarded60,61. The irregular motor unit firings were identified by following the guidelines in literature and accuracy of identified 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 first 60 ms. Moreover, given the fluctuations 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 filter and baseline correction with the mean resting torque for the first 3.33 seconds. This analysis was performed offline in MATLAB (version R2019b).
Analysis of firing rate activity. Only for SCS03 and SCS04, we computed the firing 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 identified firing 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 figures correspond to all the repetitions under the same conditions and session. This analysis was performed offline in MATLAB (version R2019b).
Statistical analyses
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 classified as outliers and discarded. For illustration purposes, the mean distribution of the data with 95% confidence interval was indicated for each distribution. This analysis was performed offline in MATLAB (version R2019b).