Abstract
Various computational approaches have been applied to predict aspects of animal behavior from the recorded activity of populations of neurons. Here we invert this process to predict the requisite neuromuscular activity associated with specified motor behaviors. A probabilistic method based on Bayes' theorem was used to predict the patterns of muscular activity needed to produce various types of desired finger movements. The profiles of predicted activity were then used to drive frequency-modulated muscle stimulators to evoke multijoint finger movements. Comparison of movements generated by electrical stimulation with desired movements yielded root mean squared errors between ∼18 and 26%. This reasonable correspondence between desired and evoked movements suggests that this approach might serve as a useful strategy to control neuroprosthetic systems that aim to restore movement to paralyzed individuals.
- bayesian statistics
- electromyography
- kinematics
- neuroprosthetics
- functional electrical stimulation
- motor control
Fundamental insights into how arrays of neurons encode motor or sensory variables can be gained from computational methods that attempt to reconstruct or predict aspects of animal behavior or sensory stimuli from the recorded activity of neural populations (Georgopoulos et al., 1986, 1988; Schwartz, 1993; Wilson and McNaughton, 1993; Deadwyler and Hampson, 1997; Rieke et al., 1997;Brown et al., 1998; Nicolelis et al., 1998; Zhang et al., 1998;Wessberg et al., 2000). The accuracy with which a behavior such as the direction of limb movement or the path of an animal navigating a maze can be reconstructed provides an estimate of the amount of behaviorally relevant information represented in the discharge of the recorded neurons.
It should also be possible to invert this process to predict neural activity from behavior. One application of such an approach would be to identify the patterns of neuromuscular activity across a population of muscles needed to elicit desired movements in paralyzed individuals using functional electrical stimulation. Functional electrical stimulation involves artificial activation of paralyzed muscles with implanted electrodes (Keith et al., 1988; Hoshimiya et al., 1989;Kilgore et al., 1989; Smith et al., 1998) and has been used successfully to improve the ability of quadriplegics to perform activities for daily living (Mulcahey et al., 1997). The range of motor behaviors that can be generated by functional electrical stimulation, however, is limited to a relatively small set of preprogrammed movements, such as hand grasp and lateral and palmer pinch (Triolo et al., 1996).
A broader range of movements has not been implemented primarily because of the substantial challenge associated with identification of the patterns of muscle stimulation needed to elicit specified movements. Most limb movements, even those involving a single digit, require intricate coordination among multiple muscles that act across several joints (Schieber, 1995; Rose et al., 1999). Such complex mechanical systems do not readily lend themselves to deterministic solutions. Although electromyographic (EMG) signals recorded from able-bodied subjects can be used to identify patterns of muscle activity associated with a particular movement (Hoshimiya et al., 1989), this painstaking method yields control signals appropriate only for the motor task from which the EMG signals were originally recorded.
In an attempt to overcome this limitation, we have used a probabilistic method called Bayes' theorem to predict the patterns of muscle stimulation needed to produce, in theory, an unlimited set of movements across multiple joints. Our use of Bayes' theorem was based on previous studies that used this method to reconstruct various forms of motor behavior from recorded neural activity (Brown et al., 1998; Zhang et al., 1998; Tresch and Kiehn, 2000). The bidirectionality of Bayes' theorem facilitated the inverse prediction of neuromuscular activity from behavior required for the present investigation (Rieke et al., 1997). The aim of this study, therefore, was to determine whether implementation of Bayes' theorem was an effective method for predicting the muscle stimulation patterns needed to artificially evoke a variety of finger movements. A reasonable correspondence between desired and evoked movements was observed in this study, indicating that this approach might provide a flexible means to control functional electrical stimulation and thereby expand the repertoire of motor functions available to paralyzed individuals. An abstract of this work has been published previously (Seifert et al., 2001).
MATERIALS AND METHODS
Overview. The general approach taken in this study involved two stages, as outlined in Figure1. In the first stage, EMG and joint kinematic signals were recorded during a variety of finger movements in one subject. These signals were then used as inputs to a computer algorithm that characterized the relationship between muscle activity and kinematics using a probabilistic method known as Bayes' theorem. In the second stage, the probabilistic relationship between muscle activity and kinematics identified in the first stage was used to predict muscle activity associated with a new set of intended or desired movements of the finger. The predicted patterns of muscle activity were then transformed into frequency-modulated trains of pulses that were used to control a set of muscle stimulators to evoke finger movements in other subjects. The accuracy of the method was evaluated by comparing evoked movements with the corresponding desired movements. Details of the procedures are given in the following sections. The Institutional Human Investigation Committee approved the procedures, and all subjects gave their informed consent to participate in the study.
Diagram depicting two main stages of the experimental procedures. Stage 1, Muscle activity and kinematic signals recorded during a variety of finger movements in one subject were used to determine the probabilistic relationship between muscle activity and movement using Bayes' rule. Stage 2, The relationship established by application of Bayes' rule in Stage 1 was used to predict muscle activity associated with a new set of desired kinematics. Predicted muscle activity was transformed into frequency-modulated trains of pulses, which were used to control a set of muscle stimulators to evoke finger movements in other subjects. Evoked movements were compared with the corresponding desired movements to evaluate the accuracy of the method.
Joint angle and EMG acquisition. A healthy human subject sat in a dental chair with his or her forearm supported on a platform and stabilized in a mid-supinated position between two foam-padded rods as shown in Figure 2. Three flexible strain gauge transducers (Biopac, Santa Barbara, CA) were used to record joint angles from the metacarpalphalangal (MCP) joint, the proximal interphalangeal (PIP) joint, and the distal interphlangeal (DIP) joint of the third digit. This digit was used because fewer muscles insert onto it compared with the thumb, index finger, or little finger and because of its greater independence of movement compared with digit four (Robinson and Fuglevand, 1999; Häger-Ross and Schieber, 2000). The joint angle transducers were attached with double-sided tape across each joint after the subject had donned a vinyl glove. The glove was worn to improve adhesion of the transducers. A plastic extension was glued to the glove over the fingernail to lengthen the distal segment and thereby allow the transducer to be fixed across the DIP joint. Once attached to the subject, each of the joint angle transducers was calibrated using a metal frame that held the joints at specified angles. Angular position was measured with respect to a neutral (fully extended) orientation of the joints with positive angles referring to flexion and negative angles indicating hyperextension. Joint angle signals were amplified (gain of 1000; World Precision Instruments, Sarasota, FL) and sampled with a computerized data acquisition system (Spike 2; Cambridge Electronics Design, Cambridge, UK) at ∼2000 Hz.
Experimental setup for recording joint angles and muscle activity and for stimulating muscle. The subject's arm was supported on a horizontal platform, and the wrist was secured in a mid-supinated position between two padded rods. Strain gauge transducers that measure joint angle were fixed over the MCP, PIP, and DIP joints of the middle finger. Tungsten microelectrodes, inserted through the skin, were used to record muscle activity from or to stimulate the middle finger (digit 3) compartments of the FDP3, FDS3, and ED3.
Tungsten microelectrodes were used to record EMG signals from the main muscles that control flexion and extension of the third digit [i.e., the digit 3 compartments of the flexor digitorum profundus (FDP3), flexor digitorum superficialus (FDS3), and extensor digitorum (ED3)]. The tungsten microelectrodes (1–5 μm tip diameter, ∼3 mm of insulation removed from the tip, 250 μm shaft diameter; Frederick Haer Co., Bowdoinham, ME) were inserted through the skin and directed toward the target muscle. Low-intensity constant current pulses (∼0.4 mA, 1 msec duration, and 1 pulse/sec) were delivered via a stimulator coupled to a stimulus isolation unit (S88 and SIU7; Grass Instruments, West Warwick, RI), while the intramuscular electrode position was adjusted manually until a site was found that elicited motor responses in one of the target muscles. Activation of FDP was distinguished from that of FDS by the presence of evoked movements in the distal phalanx. Once the placement of the electrodes in the target muscles had been verified by electrical stimulation, the electrodes were then connected to alternating current-coupled differential amplifiers (model 12; Grass Instruments). Surface electrodes (Ag-AgCl, 4 mm diameter) attached to the skin over the distal radius served as reference electrodes. EMG signals were amplified with a gain of 1000, bandpass filtered (30–1000 Hz), and digitally sampled at ∼2000 Hz.
Training data. Once the position transducers and electrodes were in place, the subject was asked to perform a variety of unrestrained flexion–extension movements of the middle finger in which contact was not made with external surfaces. Some movement of the other fingers also occurred inadvertently. However, EMG and joint angle data were recorded only for the middle finger movements, which were used for subsequent training of the Bayes' algorithm to yield the probabilistic relationships between muscle activity and joint kinematics. The movements were designed to cover much of the joint space associated with relatively natural movements. The duration of the training set was 60 sec.
Desired movements. Next, the subject was instructed to make a sequence of movements from which a set of desired movements was extracted from the recorded joint angle trajectories. These movements consisted of repeated tapping motions similar to key presses, pushing movements involving simultaneous extension of the PIP and DIP joints and flexion of the MCP joint of the middle finger in a motion similar to that which occurs when sliding a small object away from the hand across a flat surface, and pulling movements involving flexion of the PIP and DIP joints and extension of the MCP joint as if sliding a small object toward the hand. The three types of movements were performed repetitively for ∼10 sec each. The subject was instructed to make movements at a comfortable pace but to vary the duration of the movement from one cycle to the next. The entire 30 sec sequence was performed twice. From this record, five 10 sec segments were extracted that were used to represent different types of desired movements: tapping, pushing, pulling, transition from pushing to tapping movements, and transition from tapping into pulling movements. EMG data were collected during the desired movements and used for comparison with the predicted patterns of EMG.
Signal processing. In off-line digital analysis of the training set, EMG signals were full-wave rectified and low-pass filtered at 2 Hz. Joint angular velocities were calculated for each joint by digital differentiation of the joint angle data. Positive values for joint angular velocity indicated flexion movements, whereas negative values indicated extension movements. Joint angle, joint angular velocity, and EMG signals were all resampled at ∼200 Hz/signal. EMG magnitude was normalized to a percentage of the peak EMG within the training set and rounded to the nearest 1% increment. Joint angles and joint angular velocities were rounded into intervals of 1° and 1°/sec, respectively.
Bayesian reconstruction algorithm. Bayes' theorem is a technique that uses conditional probabilities to predict the likelihood of an outcome given that a particular event or set of events has occurred. The basic form of Bayes' theorem can be written as follows:
Equation 1where P(A‖B) is the probability that variable A takes on a particular value given different levels of variable B. In neurophysiology experiments, A is often a controlled parameter related to a sensory stimulus or a behavior, and B typically is an index of neural activity. P(B‖A) is the probability that variable B attains a specific value given different levels of A. P(A) is the distribution representing the probabilities for observing different levels of A. In practice, the denominator termP(B) is treated as a normalization constant that represents the sum of probabilities across all levels ofA for the distribution indicated in the numerator of Equation 1, namely:
Equation 2This normalization simply ensures that total probability represented by Equation 1 is equal to 1.
In the present case, the variables of interest were joint kinematics (θ) and muscle activity (EMG). In contrast to previous studies, in which neural activity has been used to predict some aspect of behavior (Georgopoulos et al., 1986, 1988; Schwartz, 1993; Wilson and McNaughton, 1993; Tresch and Kiehn, 2000) or features of sensory stimuli (Rieke et al., 1997; Nicolelis et al., 1998), our goal was to predict the requisite neuromuscular activity needed to generate a particular motor behavior. Consequently, the general form of equation 1 became:
Equation 3In our application of Bayes' theorem, six kinematic variables [three joint angle trajectories (Θj) and the three associated joint angular velocities (
] were used to predict activity in a muscle (EMGi). Equation 3 was first applied individually for each of the six kinematic parameters. Then, under the simplifying assumption of independence among kinematic parameters, the probability of EMG given values for all six kinematic parameters was given by the product of the individual probabilities, namely:
Equation 4The assumption of independence for angles and angular velocities across the joints of a finger is not altogether valid. To account for relationships among the kinematic parameters, however, would have required a substantially more complex process for the computations. Consequently, we opted for a more tractable form for predicting EMG represented by Equation 4 at the possible expense of some loss in accuracy and theoretical rigor. Furthermore, although other nonprobabilistic methods could have been used to predict EMG from kinematics, Zhang et al. (1998) have shown that the Bayesian approach provides accurate reconstruction of spatial motor behavior from the activity of neurons also when implementing the simplifying assumption of independence among neurons. Moreover, Zhang et al. (1998) showed that reconstruction using the Bayes' method was somewhat better than nonprobabilistic methods, such as population vector coding and template matching.
Equation 4 was applied separately to determine the muscle activity pattern for each of the three muscles (i.e., for EMG1, EMG2, and EMG3). Once the probabilistic relationships between joint kinematics and muscle activity had been established by application of Bayes' theorem on the training data, a set of new joint angles and angular velocities could be entered into the algorithm to predict the associated patterns of muscle activity (Fig. 1). For convenience, in the present case, the new set of kinematic data was obtained from the same subject from whom the training data were obtained. In theory, however, any set of desired joint trajectories could be used to predict muscle activity patterns. Ten second segments of kinematic data recorded during a variety of finger movements (but not used in the training of the algorithm) were used as inputs to the Bayes' algorithm for prediction of muscle activity. Longer segments were not used because of limitations in memory associated with generation of timing files needed to control the muscle stimulators.
Muscle stimulation. The predicted patterns of muscle activity that were based on the desired movement trajectories were converted into frequency-modulated trains of constant current pulses. To reduce computation time, successive 100 msec epochs of predicted muscle activity were consolidated into a single average value. Stimulus frequency was then linearly related to the amplitude of the average muscle activity and held constant over the 100 msec period. Stimulus frequencies ranged from 10 to 50 Hz for predicted EMG values from 20 to 100% of the peak EMG obtained in the training set. Stimulus frequencies between 10 and 50 Hz approximately correspond to the range of firing rates recorded in human motor units during voluntary contraction (Bellemare et al., 1983). The long time constant associated with the low-pass filtering of the rectified EMG led to a relatively slow decay of the EMG after a burst, such that the filtered EMG often did not reach baseline levels before the onset of a subsequent burst. Therefore, to avoid continuous stimulation of muscle because of this filter-induced prolongation of EMG, activity levels below an arbitrarily chosen threshold value of 20% of the peak EMG were not converted into stimulus pulses.
In separate sessions on five subjects (one of whom, subject A, was the subject from whom the training data were obtained), tungsten microelectrodes with ∼3 mm of insulation removed from the tip were placed into the same muscles recorded from during the training session. Joint angle transducers were applied in the same way as described above. The electrodes were connected to three independent stimulators and isolation units. The amplitude of the current pulses (1 msec in duration) was then adjusted independently for each stimulator. These adjustments were made while delivering 1 sec trains at 30 Hz to each muscle. Once a stimulus intensity was found that evoked, based on subjective criteria, a moderately brisk movement that spanned ∼50–75% of the joint range of motion, the stimulus intensity was then maintained at that level for the remainder of the experiment.
Three channels of a digital–analog converter (Spike2; Cambridge Electronics Design) were then used to deliver pulse sequences associated with the desired movements to trigger the three stimulators. During these trials, subjects were encouraged to relax the hand and not to resist the movements generated by the stimulators. The subjects were not informed of the specific type of finger movements that were to be evoked. The initial resting configuration of the finger was not specified. Five trials were evoked for each of five types of movements (tapping, pulling, pushing, tapping followed by pulling, and pushing followed by tapping). Each trial consisted of a 10 sec set of three pulse sequences delivered simultaneously to the three muscles.
Data analysis. The resulting evoked movements were recorded using three position transducers as described previously. The joint angle trajectories of both the evoked and desired movements were normalized such that the maximum joint angle within each trial was set to 100%, and the minimum joint angle was assigned a value of 0%. The evoked joint angles were then compared with the desired joint angles by calculating the root mean square (rms) difference over the 10 sec trial. The average rms error for each subject was calculated over five trials for each joint and movement. Statistical analysis of rms error was performed using a two-way repeated-measures ANOVA with joint and desired-movement type as factors. Post hoc assessment of significant differences across levels within a factor was performed using a Tukey test. Differences among means were considered to be significant for p < 0.05.
RESULTS
Figure 3 shows a segment of the training data consisting of the unprocessed EMG signals and the corresponding rectified and smoothed EMG (RS-EMG) signals from the three muscles, the three joint angle trajectories, and the three joint angular velocities obtained while the subject (subject A) performed unrestrained movements of the middle finger. In this example, as was often the case, the kinematic pattern for the PIP and DIP joints was very similar.
Short segment of the data recorded in one subject (subject A) during unrestrained movements of the middle finger that was used as input to Bayes' theorem to establish probabilistic relationships between muscle activity and joint kinematics. The bottom six traces show the unprocessed EMG signals and associated full-wave RS-EMG signals recorded from digit 3 compartments of ED3, FDS3, and FDP3. The top six traces show the joint angle and joint angular velocities for the MCP, PIP, and DIP joints of digit 3.
The method by which conditional probability distributions were constructed from data such as those shown in Figure 3 is shown as a diagram in Figure 4. For example, the joint angle values for joint 1 (θ1) associated with an activity level of 20% of the peak EMG in muscle 1 (Fig. 4,blue arrows) or 30% of the peak EMG (Fig. 4, red arrows) were used to generate the conditional probability distributionsP(θ1‖EMG1 = 20%) andP(θ1‖EMG1 = 30%) shown in Figure 4, B and C, respectively. To aid in visualization, these distributions were then depicted as strips of colored elements, with hot colors indicating high probability and cool colors representing low probability. The resulting distributions were plotted on the joint probability distribution,P(θ, EMG), as vertical bands shown in Figure4D. This process was repeated for each 1% increment in EMG amplitude to fill the entire space defined by the joint probability distribution in Figure 4D. Once completed, the color at any location on this plot indicated the probability that muscle 1 attained a particular value of EMG when joint 1 was at the specified angle.
Diagram depicting method by which joint probability distributions were generated for each combination of kinematic parameter and EMG signal. A, For each increment in magnitude of EMG activity (e.g., 20% of peak EMG,blue horizontal line; 30% of peak EMG, red horizontal line) for muscle 1, the corresponding joint angle (θ) values (vertical arrows) from joint 1 were used to construct the conditional probability distributions,P(θ1 EMG1 = 20%), shown in B, andP(θ1 EMG1 = 30%), shown in C. The histogram representations shown inB and C were transposed into bands of colored elements based on the scale depicted on the ordinate. These colored representations are depicted in D as vertical bands for which the color of any element indicates the probability that joint 1 passed through the angle specified by the ordinate given that the EMG in muscle 1 attained a value of 20 or 30% of the peak EMG. This process was repeated for each increment in EMG magnitude to fill in the surface shown in D representing the joint probability distribution, P(θ1, EMG1).
A set of six joint probability distributions (one for each kinematic parameter) was generated for each of the three EMG signals. From these joint probability distributions, it was possible to predict the pattern of EMG activity given a new set of desired movements. The process by which this was done is illustrated in Figure5. For clarity, only two of the six kinematic parameters representing the desired movement are shown. Figure 5A depicts a section of the desired angular trajectory for the MCP joint and the corresponding desired angular velocity of the MCP joint. The vertical line on the desired movement trajectories represents the time at which a prediction of the EMG was to be made. The horizontal arrows in Figure 5A indicate the specific values of the desired MCP angle and MCP angular velocity at that instant.
Conversion of desired movements into predicted EMG by application of Bayes' theorem. For clarity, only two of the six kinematic parameters used in the prediction of EMG in one muscle (ED3) are shown. A, Short time segment of desired angular trajectory for the MCP joint and associated angular velocity (Vel.) for the MCP joint. Positive angular velocities represent flexion, and negative angular velocities indicate extension movements. At the time instant indicated by the vertical line, the corresponding desired values for MCP angle and angular velocity were ∼25° and −200°/sec, respectively. These values were then used to select from the joint probability distributions, P(θ, EMG) (B, derived from training data recorded in one subject), the conditional probability associated with the desired joint angle, and joint angular velocity (thin rectangles superimposed on thecolor plots). EMG values are represented as a percentage of the peak EMG value recorded during the training set, and only EMG values above the threshold level for converting to stimulus pulses (20% of peak EMG) are shown. The specific conditional probabilities indicated on the color plots are redrawn as histograms (C) and represent the likelihood that a kinematic parameter θ, such as MCP angle, will attain a specific value,y (e.g., 25°), given different levels of EMG. These histograms were then multiplied by the overall probability of observing different levels of EMG, P(EMG) (D). The resultant histograms are shown inE. The resultant histograms were then divided by the total probability in the histogram (Σ) to yield the normalized histograms shown in F. The normalized histograms were then multiplied together to yield the conditional probability distribution shown in G,P(EMG θ1,
), which represents the likelihood of obtaining different levels of EMG given that the MCP joint angle is 25° and the MCP joint angular velocity is −200°/sec. The average value of that distribution (large arrowhead) was used as the best estimate of the EMG given the specified values of the kinematic parameters. In the present application, six kinematic parameters (angles and angular velocities for MCP, PIP, and DIP joints) were actually used to predict the most likely value of the EMG. This process was repeated for each increment in time over the entire trial.
Figure 5B shows the joint probability distributions for MCP angle versus ED3-EMG (Fig. 5B, top) and for MCP angular velocity versus ED3-EMG (Fig. 5B,bottom). The conditional probabilities associated with the specific values of the desired kinematics at the time instant in question (i.e., the regions of the color plots within thethin rectangles) are redrawn as histograms in Figure5C. These histograms represent the conditional probability that a kinematic parameter, θ, will attain a specified value,y, given different levels of EMG, namely,P(θ = y‖EMG). Then, in accordance with Bayes' theorem, these histograms were multiplied by the overall probability of encountering different levels of EMG during the training trial, i.e., P(EMG) (Fig. 5D). In this case, the probability distribution, P(EMG), was relatively uniform over different values of EMG. Consequently, the shapes of the resulting distributions, P(θ = y‖EMG) ×P(EMG), were similar to the original P(θ =y‖EMG) distributions (Fig. 5E). For normalization purposes, these resultant probability distributions were then divided by the total probability in that distribution, Σ (see Eq. 3). The normalized distributions, P(θ =y‖EMG) × P(EMG)/Σ, are shown in Figure5F. The total probability across each of these normalized distributions has a value of 1. This ensured that each kinematic parameter provided equal weight in the prediction of EMG.
Bayes' theorem specifies that the normalized distributions in Figure5F are equivalent to the conditional probability of observing various levels of EMG given the specified value of the kinematic parameter, namely, P(EMG‖θ =y) (Eq. 3). Therefore, to estimate the most likely value of the EMG given the simultaneous occurrence of the specific values of MCP angle and MCP angular velocity, the normalized histograms in Figure5F were multiplied together (Eq. 4). The outcome of that multiplication is shown in Figure 5G. A measure of the central tendency of that distribution was calculated (mean,arrowhead), and that value was then used as the predicted value of the EMG for that time instant (Rieke et al., 1997). In the present application of Bayes' theorem, six kinematic parameters (angle and angular velocity for each of three joints) were actually used to predict the most likely value of the EMG. This process was then repeated at each increment in time over the duration of the desired movement trial. The same procedure was performed separately to predict the EMG activity for each of the three muscles.
The landscape of the joint probability distribution for MCP joint angle versus EMG was relatively uniform over most of the surface shown in Figure 5B (top). This was also the case for the PIP and DIP joints. Consequently, the shapes of the conditional probability distributions associated with different values of joint angle were similar. Hence, the ability of joint angle, by itself, to predict different levels of EMG was poor. This was in contrast to joint angular velocity, in which a systematic change in the conditional probability distribution occurred for different values of angular velocity. For example, in Figure 5B (bottom), as extension angular velocity increased (i.e., increasing negative values), there was a progressive shift in the probability density toward higher values of ED3 EMG. Therefore, the inclusion of joint angular velocity was important for the prediction of EMG.
Once EMG activity for a desired set of kinematic data was predicted, the EMG signal was converted into a frequency-modulated pulse pattern. Figure 6A shows a typical predicted EMG signal derived from the process outlined in Figure 5 for a push–tap movement. For comparison, Figure6A also shows the actual ED3 EMG signal recorded (but not used in the prediction) during the trial to obtain the desired movements. The correspondence between actual and predicted EMG was quite good, with an average rms error (±SD) across all muscles and movements of 12.1 ± 3.2% (range, 6.9–16.2%). The dashed horizontal line indicates the threshold level for converting predicted EMG into stimulus pulses. Figure 6B shows the stimulus pulse pattern resulting from the predicted EMG signal shown in Figure 6A. Stimulus frequency was a linear function of the EMG amplitude such that the greater the amplitude of the EMG, the higher the frequency of pulses.
An example of predicted EMG and the corresponding stimulus pulse pattern. A, Predicted EMG (thick trace) for the ED3 muscle using Bayes' theorem based on a set of desired kinematics. Superimposed on this trace is the actual ED3 EMG (thin trace) recorded (but not used in the prediction) during the trial used to obtain desired movements. The rms error between actual and predicted EMG was 8.8% of peak amplitude of the actual EMG. The dashed line indicates the threshold below which conversion to stimulus pulses was not performed.B, Timing of stimulus pulses derived from the predicted EMG shown in A. Stimulus frequency was linearly related to the amplitude of muscle activity. Frequencies ranged from 10 to 50 Hz corresponding to 20–100% peak value of the EMG obtained in the training set.
The stimulation pulse patterns derived from the predicted EMG signals for the three muscles were used to trigger three separate stimulators. In general, the resulting evoked movements were highly consistent over repeated trials. For example, Figure7A shows the angular displacement of the MCP joint superimposed for five repeated trials of evoked tapping movements in one subject (subject C). The reproducibility of the movement was good, indicating that factors such as fatigue or electrode movement did not noticeably affect the evoked responses over the course of the experiment. Indeed, in only one case in a different subject (subject A) did the pattern of evoked movement change markedly over the course of five trials for one type of movement. This was probably attributable to migration of one of the electrodes outside the target muscle. In this case, only the first two trials were used in the analysis. For all other movement types and subjects, all five trials were included in the analysis.
Example of evoked tapping movement in one subject (subject C). A, Five superimposed tracesof the evoked angular displacement of the MCP joint.B, Desired angular displacement of the MCP joint (thin trace) and one trial of evoked movements of MCP (thick trace). Differences in joint angle bias and differences in magnitude of movements led to a relatively large rms error between two traces of 37.9% of maximum angular displacement of the desired movement. C, After normalization of thetraces in B, the correspondence in the pattern of motion at the MCP joint between desired and evoked trials is good (rms error, 13%).
Figure 7B compares the joint angle trajectory for one of the trials of evoked movement shown in Figure 7A with the desired trajectory. A relatively constant bias in joint angle between desired and evoked movements led to a large rms error (38% of the maximum angular displacement of the desired movement) despite similarities in the underlying pattern of motion. After normalization (see Materials and Methods), the correspondence between desired and evoked movements was good, as reflected in the comparatively low rms error value of 13% (Fig. 7C).
Figure 8 shows an example of normalized evoked and desired trajectories for the three joints of the finger in one subject (subject B) for a movement that involved a transition from tapping into pulling at ∼4.5 sec into the trial. During the initial tapping portion of the trial, all three joints moved more or less in phase, flexing and extending together. However, little angular displacement was evoked at the DIP joint during this phase, because the predicted level of EMG activity for the FDP (the only muscle that acts to flex the DIP joint) was less than the threshold level (20% of the peak EMG) set for conversion into stimulus pulses. At the time of transition from tapping to pulling (∼4.5 sec), a brisk extension of the MCP joint altered the phase relationship among the joints such that extension of the MCP joint then occurred while the PIP and MCP joints were flexing. This subtle transition in the phase relationship among the joints was reproduced with good fidelity in the evoked movement. During the latter pulling phase of the movement, the greatest discrepancy between evoked and desired movement was in the PIP joint, where the depth of flexion was shallower for the evoked compared with the desired movements. However, a good match between evoked and desired movements occurred during this phase for the MCP and DIP joints. Over the trial, rms errors were 15.4, 19.1, and 15.4% for the MCP, PIP, and DIP joints, respectively.
Normalized angular trajectories for MCP, PIP, and DIP joints during a trial involving a transition from tapping to pulling movement at ∼4.5 sec into the trial for one subject (subject B). Thick traces indicate evoked movements; thin traces indicate desired movements. The rms errors ranged from 15 to 19% in this trial. The drawing of the hand was adapted fromHepp-Reymond et al. (1996).
An issue of particular interest in the present study was whether prediction of muscle activity based on kinematic and EMG measurements taken from one subject could be used to generate desired movements in other subjects. In most cases, the pattern of evoked movements was similar across subjects. For example, Figure9 shows evoked trajectories of the MCP joint for the five subjects during trials involving tapping motion. Also, superimposed on these traces is the desired trajectory. Qualitatively, there was a good correspondence in the pattern of evoked movements across the five subjects, particularly in the timing of transitions from flexion to extension. However, some minor differences existed in the relative magnitudes of the movements across subjects at different phases of the trial. Quantitatively, the match between the desired and evoked movements in these trials was quite good for all subjects, with rms errors ranging from 12 to 17%. Furthermore, the correspondence between evoked and desired movements was no better for the subject in whom the original training data were obtained (subject A) than for the other subjects. This was true across all movements, as revealed by a one-way ANOVA in which no statistical difference in the magnitude of rms errors was detected across subjects for all evoked movements. Therefore, patterns of finger muscle activity predicted from data obtained in one subject can be used as templates to generate finger movements in other subjects that are reasonably close to desired movements. Whether this holds for more complex movements involving more muscles and joints is yet to be determined.
Normalized angular displacement of the MCP joint during tapping motion evoked in five subjects. Eachtrace was obtained in a different subject. Thedashed trace indicates the desired trajectory. Overall, there was a high degree of consistency in the evoked movements across subjects, which corresponded well to the desired movements. rms errors: subject A, 16.5%; subject B, 11.9%; subject C, 12.8%; subject D, 15.4%; subject E, 15.7%. Only data from subject A were used to train the Bayes' algorithm and predict the muscle activation patterns used to evoke movements in all subjects.
To evaluate the overall performance of the Bayes' stimulation technique, the rms errors between the desired and evoked movements for all five subjects are summarized in Figure10. For each subject, the average rms error over the five trials of each 10 sec movement sequence was calculated for each joint. Figure 10A shows the mean rms error and SD for the five types of movements tested (tap, push, pull, push tap, and tap pull) averaged across the three joints for all subjects. The normalized rms errors range from 17.8 to 26.5%. ANOVA indicated there was a significant difference among the mean rms error values across the different types of movements. Post hocanalysis revealed that the only significant difference in rms error among movements was between tapping and pushing movements. The lower error in tapping may have been caused by the relative simplicity of this movement, which involved alternating flexion and extension of all three joints together, whereas the pushing task required a more complex coordination with the MCP joint, flexing and extending out of phase with the other two joints. Figure 10B shows the normalized error for the different joints (MCP, PIP, and DIP) across all movement conditions. The average errors ranged from 21.8 to 23.8%, with no statistical difference in the amount of error measured for different joints. Overall, the errors were relatively modest, suggesting that the evoked movements corresponded reasonably well to the desired movements.
Mean and SD rms errors between desired and evoked responses across different movements and different joints.A, rms error across five different movements: tapping, pushing, pulling, pushing to tapping, and tapping to pulling. Errors range from 18 to 26%. The rms error for tapping was significantly less than that for pulling. B, rms error across three different joints (MCP, PIP, and DIP). There was no significant difference in rms error across joints.
DISCUSSION
Here we have shown that it is feasible to estimate the patterns of neuromuscular activity associated with a range of multijoint finger movements and to use those patterns to evoke desired movements with good fidelity using electrical stimulation. The foundation of our approach was based on previous studies that have used Bayes' theorem to reconstruct features of motor behavior from the activity of neural populations (Brown et al., 1998; Zhang et al., 1998; Tresch and Kiehn, 2000). An implicit aim in those studies was to evaluate the amount of information contained within the activity of neural ensembles related to the behavior under study. In the present investigation, we inverted this approach and used Bayes' theorem to estimate the activity in a small ensemble of muscles based on the motion of a multijoint system. We have also shown that the probabilistic relationship between EMG and kinematics derived from one individual can be used to predict patterns of activity appropriate to control muscles in other individuals. The practical importance of this finding is that existing functional electrical stimulation systems using chronically implanted electrodes (Kilgore et al., 1989; Smith et al., 1998) but implementing a probabilistic control strategy as described here could be trained on an able-bodied subject and then be deployed in paralyzed individuals.
The relationship between muscle activity and joint kinematics has been explored previously using other analytical techniques. One approach has been to predict muscle force from EMG activity using Hill-type models of muscle dynamics (Hof and Van den Berg, 1981; Olney and Winter, 1985;Winters and Stark, 1987; Soechting and Flanders, 1997). Predicted muscle forces are then used as inputs to a linked-segment model of a joint system to estimate joint kinematics using classical equations of motion (Zajac and Gordon, 1989; Kashima et al., 2000). Although this type of approach has provided an important framework for understanding the control of limb movements, such analytical methods are extremely complex, even for one or two joint systems (Winters and Stark, 1987;Zajac and Gordon, 1989), and are susceptible to several sources of error (Soechting and Flanders, 1997).
Recently, artificial neural networks have been implemented in an attempt to predict limb trajectory from muscle activity (Cheron et al., 1996; Au and Kirsch, 2000). From a practical standpoint, the advantage of this approach is that there is no need to specify an explicit algorithm that represents the complex set of interactions by which activity in several muscles is transformed into movement of a limb possessing multiple degrees of freedom. Instead, the interconnected elements that comprise the artificial neural network learn relationships among a set of input and output variables when trained with example data. Such neural networks have been shown to yield excellent predictions of complex arm movements based on EMG activity recorded from several muscles (Cheron et al., 1996; Au and Kirsch, 2000).
We followed a similar approach to that involving artificial neural networks in that no attempt was made to represent the internal mechanisms that underlie the relationship between muscle activity and movement. However, unlike the studies of Cheron et al. (1996) and Au and Kirsch (2000), our goal was to predict muscle activity patterns from movements rather than to predict movement from muscle activity. For this purpose, we used Bayes' theorem to ascertain the most likely value of muscle activity given a set of kinematic variables recorded from multiple joints of a finger. Furthermore, we extended the work ofCheron et al. (1996) and Au and Kirsch (2000) in that we used predicted muscle activity associated with desired movements to drive muscle stimulators to artificially elicit finger movements. The evoked movements were reasonably similar to the desired movements (mean rms error ranged from 18 to 26%), suggesting that this approach ultimately might serve as a useful strategy in attempts to restore movement in paralyzed individuals. The desired movements used in the present study, however, were relatively simple and repetitive. Consequently, it remains to be determined whether more elaborate and episodic movements involving more of the degrees of freedom of the hand can be reproduced using the approach described here.
Disparity between desired and evoked trajectories in the present study arose primarily because of two categories of error: (1) those associated with prediction of muscle activity from joint kinematics and (2) those related to transformation of predicted muscle activity into actual muscle activity through electrical stimulation. Overall, errors associated with prediction of EMG patterns from joint kinematics were modest (Fig. 6A). One probable cause for errors associated with this first category was that we did not obtain EMG recordings from some of the intrinsic muscles of the hand (such as second and third dorsal interossei), which can assist in flexing the MCP joint of the middle finger. Their role was likely to be particularly important during movements in which the MCP joint was flexed while the PIP and DIP joints were extended, such as occurs during different phases of pushing and pulling tasks. These movements usually were those associated with the largest rms errors (Fig.10A). Therefore, in the training data, there were some movements that could not be readily accounted for in the activity of the muscles from which we did record. Inclusion of additional muscles in the training set should help to further reduce errors in the kinematic-based prediction of muscle activity patterns.
Another factor that likely contributed to EMG prediction errors was the simplifying assumption of independence among the kinematic parameters used in our implementation of Bayes' theorem. Such independence clearly was not the case, as is evident in the congruity of the joint angle trajectories depicted in Figure 3. Consequently, the prediction of EMG likely would have improved had we accounted for correlation among kinematic parameters in our representation of Bayes' theorem.
The other category of error in the present study was related to the attempt to artificially recreate the active state of the muscle developed during voluntary contraction using a frequency-modulated pulse train delivered to the muscle through a single electrode. A number of simplifying assumptions and approximations were required to implement such a transfer function. In natural muscle contraction, the force exerted by a muscle is dependent on the number of muscle fibers recruited and on the rate of action potentials imposed by the motor neurons on the active fibers (Fuglevand et al., 1993). The muscle fibers are organized into motor units, which are groups of spatially dispersed fibers innervated by branches of the same motor axon. Variation in the strength of contraction is brought about by concurrent change in both recruitment and rate coding of motor units. The intensity of the electromyogram detected with large surface area electrodes is also influenced by both recruitment and rate coding such that a fixed (and practically linear) relationship exists between muscle force and EMG (Fuglevand et al., 1993). Accordingly, the magnitude of the EMG provides an index of the active state of muscle, which in turn is related to its mechanical (force) output.
The conversion of a predicted level of EMG into muscle activation in the present study involved delivery of current pulses through intramuscular electrodes. Because only one electrode was placed in each muscle, and because the magnitude of the stimulus current was held at a fixed level, variations in muscle activity were brought about entirely through changes in rate coding. Although this method was relatively simple to implement, it did not emulate the actual process by which muscle activity is modulated. The inclusion of a means to concurrently vary recruitment and rate coding, for example, by altering both the amplitude and frequency of the delivered current pulses, would likely improve the reproduction of the active state of the muscle and thereby enhance the match between desired and evoked movements.
Nevertheless, the overall performance of the present approach was satisfactory in reproducing desired movements and would seem to justify additional exploration and improvement of the Bayes' stimulation method. One future direction would be to include contact force signals, perhaps mediated through tactile sensors, that together with kinematic signals could be used to predict EMG activity associated with tasks that involve interaction with external objects. Another logical extension of the current method would be to expand the number of muscles included in the algorithm to predict muscle activity associated with a wide range of movements of an entire limb. However, a major obstacle to the practical implementation of such a system relates to how a paralyzed individual would supply the desired movement trajectory as input to the trained Bayes' algorithm. One possible solution would be to provide a menu of stored desired movements from which the patient could select using nonparalyzed muscles. This approach, although feasible, would not take advantage of the flexible nature of the Bayes' method.
A promising alternative would be to decipher the desired movement trajectory directly from ensembles of neurons in the cerebral cortex (Nicolelis, 2001). Previous work in nonhuman primates has indicated that the activity of populations of neurons in the primary motor, premotor, and parietal cortices can be used to predict the intended direction of hand motion during reaching movements toward targets distributed in extrapersonal space (Georgopoulos et al., 1986, 1988;Kalaska et al., 1990; Schwartz, 1993; Kakei et al., 1999, 2001;Wessberg et al., 2000). Moreover, Chapin et al. (1999) and Wessberg et al. (2000) have shown that it is possible to interpret the cortical code for a desired or actual movement and to use that brain-derived signal to control a robotic device in real time. Therefore, desired trajectories extracted from cortical recordings conceivably could be used as inputs to the Bayes' stimulation method to produce movements in an arm and hand instead of in a robot (Hoffer et al., 1996). Such an integrated system would restore movement and independence to paralyzed individuals.
Footnotes
This work was supported by National Institutes of Health Grant NS 39489 (A.J.F.). We thank Drs. Bruce McNaughton, Richard Zemel, and Gail Koshland for advice and helpful suggestions and Claudia Stanescu for technical assistance.
Correspondence should be addressed to Dr. Andrew J. Fuglevand, Department of Physiology, College of Medicine, P.O. Box 210093, University of Arizona, Tucson, AZ 85721-0093. E-mail:fuglevan{at}u.arizona.edu.
