Real-time control of hind limb functional electrical stimulation using feedback from dorsal root ganglia recordings

The University of Pittsburgh Institutional Animal Care and Use Committee approved all procedures performed for this study, which utilized four adult, male cats. Each animal was initially anesthetized with a ketamine (10 mg kg⁻¹)–acepromazine (0.15 mg kg⁻¹) mixture. The cat was subsequently intubated or a tracheotomy was performed to facilitate maintenance of anesthesia using isoflurane (1–2.5%). Respiratory rate, end-tidal CO₂, rectal temperature, O₂ perfusion in the tongue or ear, heart rate and blood pressure were monitored continuously on a vitals monitor (Advisor, SurgiVet) and maintained within normal ranges. Intravenous lines were inserted into the cephalic veins bilaterally for fluid infusion (3–7 mL kg⁻¹ hr⁻¹). One limb received a one-to-one ratio of lactated Ringer's solution and 5% Dextrose while the other limb received Plasmalyte. A single-lumen polypropylene catheter (3.5 Fr, Sovereign) was inserted into the urethra to drain the bladder.

After completion of all experimental objectives, as detailed below, each cat was euthanized with an intravenous dose of potassium chloride (10 mL of 2 mEq mL⁻¹) while deeply anesthetized.

One to three stimulating electrodes were inserted into each of the primary flexion and extension muscle groups of the left hind limb joints: anterior sartorius (hip flexion), semimembranosus (hip extension), biceps femoralis (knee flexion), vastus lateralis (knee extension), tibialis anterior (ankle flexion) and lateral and/or medial gastrocnemius (ankle extension). See figure 1(a) for a diagram of the overall setup and the major muscle groups targeted for stimulation. Each muscle was exposed during surgery and electrodes were placed on or in the muscle belly. When possible, we used electrical stimulation of the muscle surface to locate the motor point, characterized by having a low threshold and strong recruitment. Most muscle electrodes were patch electrodes (epimysial) anchored to the fascia via 3-0 silk suture. Intramuscular electrodes were used to access deeper muscles. The patch and intramuscular electrodes were made in-house using multi-stranded stainless steel wire (AS 632, Cooner Wire Company). In one experiment commercial intramuscular electrodes (Permaloc, Synapse Biomedical) were inserted in each muscle with a 16 gauge needle. A needle inserted percutaneously in the abdomen or ipsilateral foot served as the return electrode.

Zoom In Zoom Out Reset image size

The placement of each muscle electrode was inexact. Non-target muscles may have been recruited when targeting ankle extension (such as soleus), ankle flexion (such as extensor digitorum longus), knee flexion (such as gracilis) and hip flexion (such as psoas). Our goal in these experiments was not to precisely replicate normal limb movements during locomotion but to obtain sufficient muscle recruitment such that the limb could be moved through the normal range of motion in a stepping-like pattern. These electrode placements and stimulation protocols were sufficient to meet that objective.

After revealing the spinal laminae by reflecting the paraspinal muscles overlying the L5-S1 vertebrae, a laminectomy was performed to expose the lumbar DRG (L5-S1) on the left side. The cat was placed in a custom frame, which supported the torso, spine and pelvis while allowing the hind limb to move freely. Microelectrode arrays (90 channel MultiPort Utah Array, Blackrock Microsystems) were inserted into the L6 and L7 DRG (10 × 4 and 10 × 5, respectively) with a pneumatic inserter (Blackrock) and the wire bundles were secured to the dura with suture. A bone screw in the iliac crest was used as a ground electrode, and a recording reference wire was placed near the spinal cord.

Neural signals from the microelectrode arrays were recorded with a signal processing system (RZ-2, Tucker Davis Technologies) at 25 kHz and band-pass filtered (300–3000 Hz). An amplitude threshold was set on each electrode channel above the noise floor and a spike event was stored each time this threshold was crossed. A spike event consisted of a time stamp and a snippet of voltage data representing the spike waveform (0.7–1.2 ms duration). Spike sorting of the waveforms on each channel was performed on the RZ-2 processor using principal components analysis and K-means clustering (K set to 2). Before attempting FES, the clusters were manually verified and combined into a single cluster when the clusters represented the same afferent or both clusters contained multi-unit activity.

During FES, stimulation artifacts were removed using a dual approach. Most artifacts were automatically rejected by the clustering algorithm. In addition, a simple artifact rejection algorithm was implemented in the RZ2 processor, as described in [34]. Briefly, if the sum of spike-threshold crossing events across all channels during a 400 µs detection window exceeded a set limit (54 channels; 60%) then all events in a corresponding 2 ms rejection window were excluded from the spike count vector. Since active electrodes were stimulated synchronously at 30 Hz, this rejection window retained over 90% of time available for neural decoding.

For each sorted unit, after artifact rejection, the RZ-2 processor calculated spike counts (C_i) in 50 ms bins. A spike count vector (SC_i) for each unit was generated by smoothing counts over prior time points, as indicated by equation (1):

$$\begin{equation} SC_i \left( v \right) = \frac{{C_i \left( v \right)}}{2} + \frac{{C_i \left( {v - 1} \right)}}{3} + \frac{{C_i \left( {v - 2} \right)}}{6}.\end{equation} \tag{ 1 }$$

In equation (1), v is a time index corresponding to 50 ms time bins for the spike count calculation of neuron i. This smoothing approach was selected to reduce jitter in the spike count vector.

Active LED markers were placed over the left iliac crest and the hip, knee, ankle joints and toes. The marker locations were recorded with a 6-camera motion capture system (Impulse, PhaseSpace Motion Capture), sampled at 120 samples per second. A custom LabVIEW (National Instruments) program was used to calculate the hip, knee and ankle joint angles from the marker locations in real-time. As in previous studies, the limb position was represented in two reference frames (figure 1(b)) [29, 36]. Intersegmental angles were used to measure the angular position of the hip (φ_H), knee (φ_K) and ankle (φ_A) joints. Segment lengths for the femur, shank and foot were measured and combined with the joint angle measurements to determine the position of the toe relative to the hip. We also measured the toe position relative to the hip expressed as a vector in polar coordinates. The length of the vector was the hip to toe distance (R) and the orientation of the vector (θ) was measured with respect to the horizontal. For simplicity, we included the three joint and two end point position variables in a single state vector (X_m; m = 1, 2, 3 for hip, knee and ankle angles, m = 4, 5 for hip-to-toe distance and orientation angle).

Both kinematic representations of the hind limb and the smoothed neural spike count vectors were streamed via User Datagram Protocol to the finite state controller (see below) at 50 ms intervals. A haptic robot (Phantom Premium 1.5HF, Sensable Technologies Inc.) was attached to the plantar surface of the left foot and programmed to create a virtual floor, rendering ground reaction forces during the stance phase of a step cycle (transition from state 4 to state 1 in figure 1(b)).

A multivariate linear regression model was used to estimate the limb position as a function of the ensemble neural spike counts [28, 29]. All observed DRG units were included in the regression model. The estimated hind limb joint angles ($X_m^*$; m = 1, 2, 3 for hip, knee, ankle) and the hip-to-toe distance and angle with the horizon (m = 4, 5) were obtained from equation (2):

$$\begin{equation} X_m^* = \hat \beta _{m0} + \mathop \sum \limits_{i \in S_N } \hat \beta _{mi} SC_i .\end{equation} \tag{ 2 }$$

In equation (2), S_N refers to the set of N units that were sorted in real-time by the clustering algorithm and $\hat \beta _{mx}$ are the least-squares estimates of regression coefficients.

A CWE stimulator (FNS16) was used to drive muscle electrodes in a monopolar configuration. This 16-channel stimulator allowed stimulus parameter updates in 'real-time' using an RS-232 interface. We developed a custom finite state controller using LabVIEW that modified the stimulation parameters based on the location of the hind limb end point. The controller used the actual toe position (X_m), measured by the motion capture system, or the estimated toe position ($X_m^*$) based on the neural recording to update the stimulation parameters. In addition, the controller could be programmed to run in an open-loop mode with each stimulation state operating for a fixed duration.

The controller operated in one of four distinct states based on the current position of the toe. The four states corresponded to different phases of the step cycle. In state 1, the controller stimulated the ankle and knee flexor muscles to lift the leg. In state 2, the hip and ankle flexors were activated to move the foot forward. In state 3, the hip flexors remained active while the ankle and knee extensor muscles were stimulated to plant the foot on the virtual floor rendered by the haptic robot. Finally, in state 4, the hip, and ankle extensors were stimulated to pull the leg backwards (against the virtual ground).

State-switch regions were defined as boxes in the four quadrants of the expected range of motion of the hind limb. As needed, the dimensions of each state region box were adjusted to match the range of motion that could be achieved during each phase of the stimulation cycle. At each time point (50 ms increments) the position of the toe was evaluated and a state change was triggered when the toe position entered one of the pre-determined state-switch regions (see figure 1(b)). This closed-loop control scheme cycled through different stimulation states, resulting in a stepping motion.

To configure the stimulation parameters, we first determined the threshold current for muscle twitch recruitment for each electrode. Next, groups of muscles were stimulated simultaneously to generate coordinated limb movements for each phase of a stepping cycle. Stimulation amplitudes (0.5–20 mA) were adjusted manually to achieve the desired movement. If necessary, additional wire electrodes were inserted percutaneously to improve recruitment of specific muscles. The stimulation waveforms were charge-balanced cathodal-leading stimuli with a half-amplitude recovery phase at a fixed stimulation frequency (30 Hz) and pulse width (200 µs). Current levels for each stimulating channel were controlled individually, in multiples of the respective threshold current.

During a training phase, the controller used the actual limb state (obtained by the motion capture system) to update the finite state controller. The stimulation parameter sequence drove the hind limb through each of the four identified gait states. This resulted in closed-loop FES control where the feedback was provided using the camera setup (see alternate feedback path in figure 1(a)). Estimates of the regression coefficients ($\hat \beta _{mx}$) were updated every 200 ms throughout the training phase using the smoothed spike count vectors (SC_i) and the recorded limb kinematics (X_m). The training phase was terminated when the estimated coefficients stabilized, which typically took 1–2 min.

Next, we used the predictions of the limb state ($X_m^*$) based on DRG afferent activity as feedback for the finite state controller. Here, the estimated toe position was inferred from the observed spike counts (SC_i) by equation (2) at 50 ms intervals. In the first experiment (cat I), a joint angle ($X_{1,2,3}^*$) reference frame was used to estimate limb position from the afferent firing rates, and used to calculate toe position. In the last experiment (cat K), the toe position in polar coordinates ($X_{4,5}^*$) was estimated directly from the afferent firing rates. For the purpose of decoding toe position during stepping, both reference frames can be used somewhat interchangeably as prior work by our group [29–31, 36] has shown that the firing rates of many primary afferent neurons correlate strongly with limb movements represented in both reference frames.

Before activating the closed-loop FES, the hind limb would be in resting position, typically hanging between states 4 and 1. When the controller was turned on, the initial target state was 1. Subsequently, the controller drove the limb from state to state to generate a complete step cycle (duration ∼2–5 s per step; see table 1). Approximately 10–15 step cycles were performed in each trial, with 5–15 min rest periods between trials.

Subsequently, we performed multiple perturbation trials in which the limb was obstructed to prevent movement to the next state transition zone. In one experiment, the foot was held manually for ≤2 s before being released. In another experiment, the haptic robot was programmed to generate a force simulating a virtual barrier for 1 s once the limb entered state 2. After that fixed time period, the virtual barrier was removed. The period of time that the barrier actually impeded the leg movement varied, as it was located between states 2 and 3.

An additional goal of these experiments was to examine the effects of FES on the neural signals that are recorded in the DRG. FES can produce two distinct types of artifact that can contaminate the limb-state information conveyed in the records of neural activity. The most common stimulation artifact is that produced by volume conduction of stimulation currents from the muscle to the recording electrodes. These appear simultaneously with the stimulation pulses applied to the muscle. A second type of stimulation artifact may be produced by sensory fibers recruited directly by electrical stimulation applied to the muscle. This type of artifact appears at a short latency (∼2–5 ms) corresponding to the propagation delay from the stimulation site to the DRG. Since this neural activity is produced by stimulation, it does not necessarily carry information about the actual movement of the leg. Although we were successful at generating closed-loop FES control using DRG recordings and artifact rejection, we performed this analysis to achieve a more thorough assessment of neural signal quality in the presence of muscle stimulation.

We examined the effects of stimulation on DRG unit responses during passive movements of the leg generated by a robotic arm (VS-6556E/GM, DENSO Robotics) attached to the base of the foot. DRG neural responses were recorded in two conditions: robot and robot­ + FES. During the robot + FES condition, electrical stimulation was applied continuously to the muscle electrodes in the same pattern and amplitudes used during closed-loop FES while the robot moved the leg. The robot arm produced the same leg movements during the robot and robot + FES conditions and allowed for a direct comparison of the neural activity generated by movement alone and movements produced in the presence of FES. The passive movements consisted of ramp-and-hold flexion–extension movements (e.g. between states 1 and 3 in figure 1(b)) or complete step cycles (e.g. all four states in figure 1(b)).

A multi-step process was followed to examine the effect of electrical stimulation on observed DRG unit responses with respect to the position and velocity of the toe relative to the hip. Only the hip-to-toe distance was used in this analysis as this kinematic parameter is highly correlated with the joint angles during these simple movements, providing a simple, yet characteristic representation of the whole-leg movements. The DRG neural spike activity was sorted offline in OpenSorter (TDT, Inc.) after each experiment. Robot + FES trials were sorted based on units identified during sorting of robot trials. Stimulation artifacts and DRG units that did not appear during robot trials were marked as outliers in robot + FES trials and excluded. DRG unit firing rates were calculated by convolving spike times with a 150 ms triangular window, as in [29].

We define the terms 'position gain' and 'velocity gain' to represent the response of a DRG neuron to changes in leg position and leg velocity. The position gain of each DRG unit was estimated by the difference in the mean firing rate between the flex and extend positions of movement, when velocity was near zero, divided by the mean change in the toe position. The velocity gain of each DRG unit was estimated by the difference in the mean firing rate between states when the limb was moving (|velocity| > 0) and hold states (velocity ≈ 0) divided by the mean difference in velocity. Figure 2 shows example limb movements and firing rates that were used for each gain calculation. Units without consistent polarity in gain across trials or an average correlation coefficient (ρ) between firing rate and the kinematic parameter less than 0.3 were excluded. These steps were performed to remove units that were not related to limb movement from the analysis.

Zoom In Zoom Out Reset image size

Position and velocity gain values obtained from the robot and robot + FES conditions for each unit were compared using the 'gain normalization' metric computed with equation (3):

$$\begin{equation} {\rm GN} = \frac{{{\rm Gain}\left( {{\rm robot} + {\rm FES}} \right) - {\rm Gain}({\rm robot})}}{{{\rm Gain}\left( {{\rm robot} + {\rm FES}} \right) + {\rm Gain}({\rm robot})}}.\end{equation} \tag{ 3 }$$

This gain normalization value (GN) yields one of three relationships: units with GN near 0 did not have a change in the unit response during electrical stimulation; units with GN near +1 had a higher effective gain during stimulation; and units with GN near –1 had a higher effective gain without stimulation. Subscripts of 'p' and 'v' denote GN for position and velocity, respectively. For the example units of figure 2, each GN was just below zero, indicating a slightly lower gain during stimulation.

The previous analysis examined the effects of stimulation on individual neurons. We also examined the global effects of stimulation on neural population decoding. Repeated multivariate linear regressions were performed on each robot + FES and robot trial in Matlab (Mathworks) to estimate the leg position from the recorded activity, for different numbers of DRG units as an input. Sorted units in each trial were ranked by their individual correlation coefficient to the hip-to-toe distance. Starting with the best unit and increasing up to all units in a trial, linear regression was performed to fit unit firing rates to the position. For robot + FES trials, regression estimates of hip-to-toe distance were repeated for several alternative unit sorting methods (all snippets including stimulus artifacts in a single unit sort per channel; all non-stimulus artifacts in a single unit sort per channel; and only stimulus artifacts in a single unit sort per channel) and compared to the limb position regression estimate for the ideal unit sorting approach, with stimulus artifacts excluded and all single units sorted appropriately.

The root mean square error (RMSE), as calculated in equation (4), was used to determine the error in limb position estimates:

$$\begin{equation} {\rm RMSE} = \sqrt {\frac{{\mathop \sum \nolimits_{u = 1}^t \left( {P_u - P_u^* } \right)^2 }}{t}} .\end{equation} \tag{ 4 }$$

In (4), P_u and $P_u^*$ refer to the measured and estimated hip-to-toe position at time u across the t data points in a trial. Confidence intervals were calculated to determine statistical differences. A p-value ≤0.05 (for a 95% confidence interval) was considered significant. Where relevant, data is reported as average ± one standard deviation.

Real-time decoding of lumbar DRG activity was performed in two of the experiments (cats I, K). In experiment I, 129 DRG units were used in the real-time model, which estimated the angular position of the hip, knee and ankle joints ($X_{1,2,3}^*$). The estimated joint angles and measured segment lengths for the femur and shank were used to calculate an estimate for the toe position. At the end of the 61 s training period, which had 29 complete step cycles, the limb position estimate had an RMSE of 6.4 mm. In experiment K, 124 DRG units were used in the real-time model, which used polar coordinates ($X_{4,5}^*$) to estimate the toe position. In this experiment, a 73 s training period consisting of 28 complete limb step cycles led to a final limb position estimate RMSE of 4.8 mm. The first 20 s of the training period is recreated for cat K in figure 3.

Zoom In Zoom Out Reset image size

Closed-loop FES control of the left hind leg was attempted in three cats (cats I, J and K). In cat J, insufficient muscle recruitment was achieved with the FES electrodes made in-house, which prevented an evaluation of closed-loop control. This failure led to the purchase of the commercial intramuscular electrodes utilized in experiment K. The results discussed here are from the two successful experiments.

Example step cycles and electrode combinations from cat K are shown in figure 4. The electrode locations and stimulation sequence used to control the limb position are indicated below the kinematic trace. Four distinct stimulation channel combinations were used for transitioning between the four states of the step cycle. As the foot was transitioned out of state 4, there were small oscillations resulting from contact with the haptic robot simulated floor.

Zoom In Zoom Out Reset image size

In experiment I, a total of 74 complete FES-controlled stepping cycles were performed using DRG signals as feedback, of which 57 cycles were not impeded with a perturbation. Across six trials, 12.3 ± 9.2 cycles were performed for 64.4 ± 20.4 s per trial. Rest periods between closed-loop trials were 3.6 ± 1.8 min. In cat K, a total of 94 complete FES-controlled stepping cycles were performed, of which 70 cycles were not impeded. Across seven trials, 13.4 ± 6.2 cycles were performed for 77.6 ± 8.2 s per trial. Rest periods between trials were 16.8 ± 8.4 min. Table 1 summarizes the effective linear stepping distance and time per experiment. Partial stepping cycles, at the beginning of a trial when the limb moved to state 1 or at the end of a trial if recording was stopped before the limb returned to state 1 are not included in this summary.

A total of 20 manual perturbations were performed in cat I (1.01 ± 0.47 s average duration) and 24 automated (i.e. robot-controlled) perturbations were performed in cat K (0.76 ± 0.11 s average duration). Table 1 summarizes the effective linear stepping distance and time during perturbation cycles, per experiment. After each perturbation, the controller correctly identified that the position of the limb had not yet reached the next state. This shows that the estimates of the limb state were not influenced by the stimulation itself and were based purely on sensory neural activity associated with the actual limb position. Once the perturbation ceased, the hind limb continued to advance toward the desired state and usually completed the step cycle successfully. After 8 of the perturbations, the controller first had to force the limb to re-enter the perturbed state before continuing due to movement that occurred during the perturbation. Two successful stepping cycles with perturbation are shown in figure 4.

Across the two experiments, 3 and 20 step errors, respectively, occurred in which the limb either re-entered a state (18) or skipped a state (5) in the step cycle. These errors occurred during 18 of the full step cycles, representing an error rate of 10.7%. Eight of these step errors occurred after a perturbation, as described above. After each error the controller was always able to complete a step cycle, although in two cases there was more than 1 error before the cycle finished.

Closed-loop FES stepping cycles were performed in several trials in each experiment, with multiple cycles per trial (12.3 ± 9.2 and 13.4 ± 6.2 cycles per trial, for I and K). Within each trial, regardless of the number of cycles, a trend was typically observed for the total path length that the toe followed: cycle path lengths increased to a maximum near the mid-point of the trial before decreasing. In figure 5, the mean cycle path lengths in five equal duration periods (20% intervals) of each trial were normalized to the maximum period mean cycle path length. Path lengths for individual trials within each experiment and across all trials are indicated. Path lengths during the first periods of all trials were significantly lower than during periods 2–4 (p < 0.05). Mean path lengths during the final period were lower than but not significantly different than during the middle period (p = 0.07). These results suggest that the flexibility and recruitment of the limb increased after the start of a trial and that fatigue may have occurred later in each limb movement sequence.

Zoom In Zoom Out Reset image size

A comparison of DRG unit decoding during electrical stimulation for robot-controlled limb movements was performed in three cats (H, I, K). In cats H and I, two-state ramp-and-hold limb movements (whole-limb flexion and extension) were performed. In cat K, four-state stepping was performed. Across the three experiments, spike snippets from the robot-only trials were sorted on 222 channels yielding 353 non-artifact units. Four channels in experiment I and 44 channels in experiment K did not yield neural data. In the robot + FES trials, we discarded units that were not present in the robot-only condition. The discarded units included waveforms corresponding to stimulation artifacts and afferent activity produced by electrical stimulation. The discarded units are referred to as artifacts in analyses described below.

From the identified units, 242 were used in this analysis (60 units in both position and velocity analyses) and 111 units were excluded. The excluded units comprised 90 units which had a consistent gain polarity but ρ < 0.3 for either the position or velocity analysis and 20 units which failed both criteria in both the position and velocity analyses. Among the excluded units were 46 clear single units and 32 multi-units which were not correlated to the limb movement and 33 noise units on channels where the snippet threshold was set too low.

Figure 2 shows the position and velocity responses of example units during trials with and without stimulation. The GN for each unit is −0.05, indicating that a lower gain was observed for each unit during the trials with stimulation. Figure 6(a) summarizes the GN_p and GN_v across all non-excluded units, per experiment and in total. In experiment I (p < 0.001) and for all units (p < 0.02), the mean GN_p was significantly lower than 0, indicating that the presence of electrical stimulation led to a lower observed mean unit gain. The mean GN_v was significantly lower than zero (p < 0.029) in all comparisons.

Zoom In Zoom Out Reset image size

The lower mean unit gains during electrical stimulation did not significantly affect the mean linear regression estimate of the hip-to-toe position. In figure 6(b) the mean RMSE of the regression estimate, for robot + FES and robot trials, is given for up to 30 rank-sorted units in the model. There was no significant difference between the means at any number of units in the model (p > 0.19). The standard deviation of the position estimate errors was larger for robot + FES trials (≤2.1 times larger), but not significantly (p > 0.07). Beyond 30 units in either trial type, the mean position estimate RMSE was between 4 and 6 mm.

We compared the accuracy of the estimated position in robot + FES trials achieved with spike data obtained with different offline unit sorting methods. First, estimates of the toe position were compared between sorted DRG units and unsorted DRG units including artifacts. The units were rank-ordered by their individual correlation with toe position and the difference in the RMSE between the two approaches was calculated (figure 7(a)). For 2–9 units, identifying individual units per channel and removing artifacts yielded significantly or near significantly better estimates (0.02 < p < 0.14) than not performing any sorting on each channel. For greater than 20 units, the RMSE difference approached 1.0 mm and was not significantly different (p > 0.14).

Zoom In Zoom Out Reset image size

To further investigate the effects of sorting, robot + FES trials were sorted in two additional means: artifacts only in a single sort code and non-artifacts only in a single sort code. Figure 7(b) shows the RMSE of the toe position estimates for the different sorting approaches using 20 units. For non-artifacts only in a single sort code in the robot + FES trials, the regression estimates were not different than when all threshold crossings were combined as a single unit per channel. These two approaches were consistently less accurate than when units are individually sorted with artifacts rejected across the various number of units in the model (1.8 ± 2.5 mm average difference for 1–30 units in each model). As indicated in the figure, end point position estimates were very inaccurate when only artifacts were in the model (9.0 ± 5.8 mm average difference against other approaches for 1–30 units in models). The end point estimates for the robot trials showed a similar trend, with a model based on unsorted activity only slightly less accurate than a model using well-sorted units (p > 0.50). In summary, artifacts or the use of unsorted units led to slightly less accurate model estimates and reduced unit gains with greater variability; however the effect was minimal for beyond 10 units in the models.