HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via ISI Web of Science (3) Citing Articles via Google Scholar Google Scholar Articles by De Luca, C. J. Articles by Nawab, S. H. Search for Related Content PubMed PubMed Citation Articles by De Luca, C. J. Articles by Nawab, S. H.

INNOVATIVE METHODOLOGY

Decomposition of Surface EMG Signals

¹NeuroMuscular Research Center, ²Department of Electrical and Computer Engineering, and ³Department of Biomedical Engineering, Boston University, Boston, Massachusetts

Submitted 4 January 2006; accepted in final form 25 May 2006

	ABSTRACT

TOP ABSTRACT INTRODUCTION BACKGROUND METHODS RESULTS DISCUSSION APPENDIX: SIGNAL PROCESSING... GRANTS ACKNOWLEDGMENTS REFERENCES

 
This report describes an early version of a technique for decomposing surface electromyographic (sEMG) signals into the constituent motor unit (MU) action potential trains. A surface sensor array is used to collect four channels of differentially amplified EMG signals. The decomposition is achieved by a set of algorithms that uses a specially developed knowledge-based Artificial Intelligence framework. In the automatic mode the accuracy ranges from 75 to 91%. An Interactive Editor is used to increase the accuracy to >97% in signal epochs of about 30-s duration. The accuracy was verified by comparing the firings of action potentials from the EMG signals detected simultaneously by the surface sensor array and by a needle sensor. We have decomposed up to six MU action potential trains from the sEMG signal detected from the orbicularis oculi, platysma, and tibialis anterior muscles. However, the yield is generally low, with typically 5 MUs per contraction. Both the accuracy and the yield should increase as the algorithms are developed further. With this technique it is possible to investigate the behavior of MUs in muscles that are not easily studied by needle sensors. We found that the inverse relationship between the recruitment threshold and the firing rate previously reported for muscles innervated by spinal nerves is also present in the orbicularis oculi and the platysma, which are innervated by cranial nerves. However, these two muscles were found to have greater and more widespread values of firing rates than those of large limb muscles.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION BACKGROUND METHODS RESULTS DISCUSSION APPENDIX: SIGNAL PROCESSING... GRANTS ACKNOWLEDGMENTS REFERENCES

 
The electromyographic (EMG) signal is composed of the action potentials from groups of muscle fibers organized into functional units called motor units (MUs). This signal can be detected with sensors placed on the surface of the skin or with needle or wire sensors introduced into the muscle tissue. When only two or three MUs in the vicinity of the sensors are active, it is usually possible to visually identify most of the individual MU action potentials because the incidence of superposition among the individual MU action potentials is relatively low. However, when the EMG signal contains the activity of four or more MUs the individual action potentials become, in large part, indistinguishable to the naked eye because the incidence of superposition among two or more MU action potentials becomes numerous and the shapes of the MU action potentials may approach in similarity.

In many cases it is desirable to study and/or use the information contained in the timing of the discharges of individual motor units, such as in investigations for furthering the understanding of how motor units are controlled by the CNS in generating force or for assessing the degree of dysfunction in upper motoneuron diseases such as cerebral palsy, Parkinson's disease, amyotrophic lateral sclerosis (ALS), stroke, and other disorders. This may be achieved by "decomposing" the EMG signal. The concept is depicted in Fig. 1. A decomposed EMG signal provides all the information available in the EMG signal. The timing information provides a complete description of the interpulse interval, firing rate, and synchronization characteristics. The morphology of the shapes of the MU action potentials provides information concerning the anatomy and health of the muscle fibers.

View larger version (40K): [in this window] [in a new window]   FIG. 1. Pictorial outline of the decomposition of the surface EMG signal into its constituent motor unit action potentials. (Adapted from De Luca et al. 1982a.)

There have been numerous and varied approaches to extracting action potentials from neural and muscle activity over the past three decades. The mid-to-late 1960s produced a flurry of computer-based activity directed at identifying the individual action potentials and discharge times of neural activity by shape discrimination. Dominant among these pioneering attempts were the works of Gerstein and Clark (1964), Simon (1965), Keehn (1966), and Glaser and Marks (1966). Applications to separate the EMG signals did not appear until a full decade later (Andreassen 1983; De Figueiredo and Gerber 1983; De Luca and Forrest 1972; De Luca et al. 1982a,b; Guiheneuc et al. 1989; LeFever and De Luca 1978; Mambrito and De Luca 1984; McGill et al. 1985). This initial flurry was followed by sustained interest during the past two decades (Broman 1988; Christodoulou and Pattichis 1999; Fang et al., 1999; Hochstein et al. 2002; Iani et al. 1994; Jongen et al. 1996; Loudon et al. 1992; Nawab et al. 2004a,b; Stashuk and de Bruin 1988; Türker et al. 1989; Zennaro et al. 2001; among many others). To date all the above approaches and techniques that have been able to identify individual MU action potentials in the superimposed EMG signal and provide useful physiological information have used indwelling sensors to detect the signal. These types of sensors have obvious disadvantages arising from their invasive nature.

In this report we describe a successful attempt at accurately decomposing EMG signals detected from surface sensors.

	BACKGROUND

TOP ABSTRACT INTRODUCTION BACKGROUND METHODS RESULTS DISCUSSION APPENDIX: SIGNAL PROCESSING... GRANTS ACKNOWLEDGMENTS REFERENCES

 
For the past three decades we have been improving a technique that we have termed Precision Decomposition I (PD I). The evolution of this technique, specifically designed to enable physiological experiments, has been described by LeFever and De Luca (1978), LeFever and De Luca (1982), Mambrito and De Luca (1984), Broman (1988), and De Luca and Adam (1999). The technique has been useful for decomposing indwelling EMG (iEMG) signals detected by needle sensors during isometric contractions and has been used in numerous physiological studies (see, e.g., Adam and De Luca 2003, 2005; Adam et al. 1998; De Luca et al. 1982a,b; Erim et al. 1999; Finsen et al. 2006; Kamen and De Luca 1989; Makasado et al. 1995; Masuda and De Luca 1991; Roark et al. 2002; Sogaard 1995; Westad and Westgaard 2005). Briefly, the technique consists of identifying action potentials in the iEMG signal and assigning them to specific motor units by classifying the shapes and amplitudes of the action potentials. The assignments of the action potentials are made on the basis of template matching and the probability of firing of the individual motor units being tracked. Slight modifications in the shape of the action potential can be tracked by adapting the templates to these modifications. Superpositions of action potentials are resolved. Artificial Intelligence techniques are used to accomplish all these tasks. These algorithms were able to decompose no more than eight MU action potential trains in the iEMG signal with an accuracy of 60 to 70% in an automatic mode, requiring 12-min processing time for 1-s acquisition time (De Luca et al. 1982a,b). The technique has provisions for using an Interactive Editor program to obtain decomposition accuracies of 100% in some records of EMG signals.

Throughout this communication, the decomposition accuracy for the ith decomposable motor unit train is defined as

where N_FIR(i) is the number of true firings of the MU and N_FN(i) and N_FP(i) are respectively the number of false negatives and the number of false positives produced by the decomposition algorithm for that MU. The term "true firings" refers to the firings that were obtained by an expert operator using the Interactive Editor on automatically decomposed data. The overall decomposition accuracy for a signal with N decomposable MU trains is then obtained as

The rationale behind this unweighted average is that the accurate decomposition of any MU train is of the same significance as that of any other MU train regardless of its duration, number of firings, and so forth.

Recently, we made significant improvements to the algorithms, rendering a considerably more powerful version that we have labeled Precision Decomposition II (PD II). The new algorithms are based on our own Artificial Intelligence knowledge-based approach specifically designed to manage signal-processing algorithms that perform two main categories of functions: 1) they identify differences in shapes and track changes in the shapes of the action potentials under a variety of conditions and 2) they resolve complex superpositions. The algorithms have been described in publications by Nawab and Lesser (1992), Winograd and Nawab (1995), Hochstein et al. (2002), and Nawab et al. (2004a,b). We succeeded in automatically decomposing iEMG signals from the quadrifilar needle sensor (De Luca et al. 1982a) and later with quadrifilar wire sensors (De Luca and Adam 1999) with typically 10 MU action potential trains with an accuracy of typically 85% with a processing time eightfold that of the acquisition time. In one extraordinary case we were able to decompose automatically an iEMG signal containing 14 MU action potential trains, of which the most significant eight could be detected with an accuracy of 98% (Nawab et al. 2004a).

The quadrifilar needle sensor has the advantage that it can be repositioned after an insertion or, if necessary, relocated, thereby increasing the probability of obtaining a quality signal that can be decomposed. The fine-wire version may be placed in deep muscles located under an overlying layer of muscle and it usually provides no discomfort once inserted; however, it has some disadvantages. Once inserted it cannot be precisely relocated within the muscle. One can pull the wire out fractions of a millimeter, although this procedure can be done only once or twice and with little control over the precise placement of the sensor.

Both versions of sensors have inherent limitations:

• They must be inserted into the muscle. This requires a clinical preparation involving sterilization of the sensors and the environment where the insertion is to be made.
  
• They render difficult or, in some cases, preclude the investigation of MU firing properties in muscles located in sensitive and dangerous areas, such as the muscles of the lips, eyelids, tongue, and most facial muscles, as well as other regions.
  
• They carry the risk—albeit low—of infection and disease transmission.
  
• They cause minor damage to the muscle tissue from which they are detecting the signal, which in turn influences the shapes of the action potentials.
  
• They are not well tolerated by individuals who have needle aversion, such as children.
  
• Once these sensors are inserted, the subject must remain very steady. A minor movement of 0.1 mm may cause the shapes of the MU action potentials to change, thus precluding the continued identification of a specific unit and generally incapacitating the decomposition algorithms from identifying actions potentials in the remainder of the contraction.
  

Many of the disadvantages and limitations of indwelling sensors can be mitigated by developing technology that uses a sensor to detect the surface EMG (sEMG) signal and algorithms that can decompose the sEMG signal. Such an approach presents difficult challenges because the sEMG signal is more complex than that detected by the relatively more selective indwelling sensors.

The notion of decomposing the sEMG signal has been considered by many for two decades (e.g., Masuda et al. 1985) for the obvious reason that it can be used by nonclinical researchers interested in motor control. Kleine et al. (2000) were able to decompose up to five concurrently active MUs obtained in a multistep, highly user interactive template-matching procedure from surface EMG signals detected with a 128-channel sensor during a low-level (5%) maximal voluntary contraction (MVC). In the past two years, there has been increased interest in the problem of decomposing the sEMG signal. Östlund et al. (2004) and Holobar and Zazula (2004) investigated the problem with various signal-processing techniques, but the few results presented were restricted to contraction levels of 10% MVC, included no direct measure of accuracy, and assumed stationary MU action potential shapes and uncorrelated discharge behavior. Gazzoni et al. (2004) described a technique that automatically decomposes 10 simulated MU trains and detects, in select cases, two to three MU trains from real sEMG signals, but makes no attempt at decomposing superimposed action potentials. The PD I technique that we used two decades ago (LeFever and De Luca 1978) was similarly tested, with perfect results, on simulated EMG signals, but failed miserably when applied to real EMG signals. We argued the essential importance of testing decomposition techniques on real EMG signals in Mambrito and De Luca (1984). Most recently, Garcia et al. (2005) presented a semiautomatic decomposition algorithm to extract the average firing rates of one to two MUs from EMG signals detected with a surface electrode array. No accuracy measures for these results were given, but the authors provide a success rate of 19/26 complete firing trains for a set of synthetic EMG signals composed of 10 MU action potential trains.

In this report we describe a technology that is capable of decomposing several simultaneous MU action potential trains from real sEMG signals. Although it is in the early stages of development, we can demonstrate its efficacy. In some cases the present version of this new technology yields an accuracy similar to that obtained by our early algorithms of PD I using indwelling sensors. We refer to our new technique for the decomposition of sEMG signals as Precision Decomposition III (PD III). When fully developed, PD III should find applications in neurology, ergonomics, sports medicine, physical medicine, and space medicine.

	METHODS

TOP ABSTRACT INTRODUCTION BACKGROUND METHODS RESULTS DISCUSSION APPENDIX: SIGNAL PROCESSING... GRANTS ACKNOWLEDGMENTS REFERENCES

 
The surface sensor array

We chose a sensor design that was similar to that of the original needle sensor, but with larger dimensions. As shown in Fig. 2A, the sensor consists of four cylindrical probes (0.5 mm diameter) with blunted ends that protrude from the housing so that when pressed against the skin they make a surface indentation. The probes are placed at the corners of a 3.6 x 3.6-mm square. The diameter of the probes was chosen to be as small as possible without piercing the skin when pressed forcefully. As in the intramuscular sensors mentioned above, each probe provides a detection surface. The leads are connected to the inputs of differential amplifiers. The four detection surfaces provide six differential combinations or channels. Signals from four of these channels are amplified and filtered with a bandwidth of 250 Hz to 2.0 kHz to remove any movement artifact at the low end of the spectrum and any excessively long tails of the action potentials. The signals are then stored. This arrangement is shown in Fig. 2B. We found that for the sEMG signals the use of four differential channels improved the accuracy of the decomposition in PD III, whereas for iEMG signal decomposition in PD II, three differential channels were sufficient.

View larger version (19K): [in this window] [in a new window]   FIG.2. A: surface electromyographic (sEMG) sensor array containing the 4 pins that detect the sEMG signals. B: differential combinations that produce 4 channels of sEMG signals. Detected signals from each channel are band-pass filtered before being decomposed by the algorithms.

Acquisition of sEMG signal

Three volunteers with no known neurological disorders participated in the tests. All subjects read and signed an informed consent form approved by the Institutional Review Board of Boston University. To detect the signal, the surface sensor array was pressed against the skin above the muscle of interest. No skin preparation or conductive gel was required. Sufficient pressure was applied to provide good electrical contact as evidenced by the best signal-to-noise ratio of the detected signals. (In our system, this is accomplished by viewing the detected signal on a computer screen in real time.) We detected signals from various muscles from the face, neck, and in limbs. Here we report on two muscles of the head and neck, the orbicularis oculi and the platysma, and on one limb muscle, the tibialis anterior. The head and neck muscles were chosen because they are not conveniently or comfortably assayed with a needle sensor and because their architecture differs from that of the more commonly studied limb muscles that have unified origins and insertions. The tibialis anterior was chosen as a typical representative of limb muscles. The sensor array was placed in the following locations. In the tibialis anterior, it was placed in the distal third (lengthwise) of the muscle; in the orbicularis oculi, it was placed on the lateral aspect of the superior eye lid; and in the platysma, it was placed on the lateral aspect of the neck, halfway between mandible and clavicle.

In the tibialis anterior, we recorded the torque output of the muscle and provided it as feedback to the subject so that a constant-force isometric contraction was generated during the detection of the sEMG signal to be decomposed. This was done by placing the limb in an apparatus that restrained the movement of the limb and was instrumented with a high-stiffness force gauge (model MB-250; Interface, Scottsdale, AZ). The torque output of the muscle generated during the contraction was expressed as a percentage of the maximal torque generated before the contraction. The subjects were asked to maintain a constant force at 20–50% MVC. In the muscles of the head and neck, where it was not possible to record the torque, the feedback to the subject was provided by the root-mean-squared (RMS) value of the EMG signal detected with a standard sEMG differential sensor (DE 2.1; Delsys, Boston, MA) placed adjacent to the array. This differential sensor consists of a parallel-bar configuration having detection surfaces that are 1 mm wide, 10 mm long, and spaced 10 mm apart. This sensor detects a global sEMG signal that is representative of the overall activity of the muscle being monitored. The global sEMG signal was detected with a bandwidth of 20–450 Hz and the RMS was calculated and smoothed using hardware circuitry. When processed in this fashion, the time-varying RMS value provides a reasonable analog to the force produced by the muscle (Basmajian and De Luca 1985; Lawrence and De Luca 1983). The time-varying value of the RMS of the global sEMG signal, or the estimated force, was expressed as a percentage of the value obtained during a maximal contraction. In the contraction of the orbicularis oculi, the subject was requested to squint with one eye and to attempt to maintain the global RMS constant. In the platysma, the subject was requested to tense the neck musculature and to maintain the global RMS constant.

Challenges of sEMG decomposition

Any approach for decomposing EMG signals must be able to deal with four major complexities that occur within the signal. These complexities are shown in Fig. 3: 1) Superposition of action potentials from different MUs, 2) Large dynamic range of the amplitudes among the action potentials of different MUs of interest, 3) Shape changes across the different action potentials of each MU (arising from slight movement between the sensor and muscle fibers and/or intracellular processes), and 4) Similarity of shape at various times among the action potentials of different MUs. These phenomena may also act in concert with each other to make the decomposition task all the more difficult.

View larger version (25K): [in this window] [in a new window]   FIG. 3. Stylized examples of the various challenges presented by the realistic behavior of EMG signals detected with indwelling sensors with small detection volume and susceptible to movement.

View larger version (18K): [in this window] [in a new window]   FIG. 4. Top: one channel of signal detected from the surface sensor array shown in Fig. 2A. Bottom: one channel of signal detected simultaneously from a needle quadrifilar sensor. Both signals were detected during constant 30% maximal voluntary contraction (MVC) isometric contractions from the tibialis anterior muscle.

The PD III algorithms used to decompose the sEMG signal are modified versions of those previously developed by us for the PD II system (Hochstein et al., 2002; Nawab et al. 2002a,b, 2004a,b). The algorithms in both PD II and PD III are designed to address the difficulties inherent in trying to separate overlapping action potentials at various points in the EMG signal by using an Artificial Intelligence framework called IPUS ("Integrated Processing and Understanding of Signals") Lesser et al. 1995; Nawab and Lesser 1992; Winograd and Nawab 1995). Although other signal-processing approaches have also been investigated for performing EMG decompositions, the IPUS approach has provided significantly more accurate results on intramuscular EMG data (Nawab et al. 2004a,b). The PD III system represents our initial attempt at incorporating the IPUS approach into the decomposition of surface EMG signals.

As illustrated in Fig. 5, the PD III system consists of four processing stages. It takes as input an sEMG signal, x(t), and produces as its final output a set {y_M^(j) (t)|j = 1, 2, ... , N_M} of N action potential trains. PD III first applies a digital IIR band-pass filter (see the APPENDIX for details) to the input EMG signal, x(t), and passes the result, x_B(t), through a maximum a posteriori probability (MAP) receiver (LeFever and De Luca 1982) designed to avoid false positives in the classification of action potential detections (see the APPENDIX for details). The MAP receiver may also be viewed as producing estimates {y_M^(j)(t)|j = 1, 2, ... , N_M} of N_M action potential trains. The corresponding decomposition error is evaluated as the square of the error signal, e_M(t) = x_B(t) – _j=1^N_M y_M^(j)(t). This error arises primarily because in trying to avoid false positives the MAP receiver 1) often splits off a motor unit train into two or more trains because of shape differences among various action potentials of the original train and 2) often fails to find a classification for signal data that is composed of complex superpositions of several different action potentials. Subsequent stages of processing are aimed at reducing these two types of decomposition error.

View larger version (14K): [in this window] [in a new window]   FIG. 5. Block diagram presenting the most important components of the Precision Decomposition III (PD III) algorithms.

The fourth PD III stage in Fig. 5 is designed to compensate for the ineffectiveness of the MAP receiver in signal regions involving complex superpositions. In such regions, significant interference between two or more action potentials makes the resultant data very different from any other data encountered in the input signal. In such cases, the MAP receiver declares the data as belonging to a new MU train but never finds a matching action potential later in the signal. Such "degenerate" trains are abandoned and their corresponding data are reanalyzed by PD III's fourth stage. The reanalysis is aimed at identifying the nondegenerate trains from the MAP receiver that are consistent with the data in the superposition regions. This reanalysis begins with an iterative correlation strategy (see the APPENDIX for details) to assign probabilities of occurrence to each nondegenerate train in each superposition region. These iterative correlation results then undergo a utility maximization process (Von Neumann and Morgenstern 1944) to convert the probabilities of occurrence into decisions about which MU trains are actually consistent with the superposition data.

Interactive Editor

The results of the automatic decomposition could be improved by using the Interactive Editor to identify and correct wrongly assigned action potentials. Figure 6 shows an actual screen-shot of the pop-up window that appears in response to a click on a gap in the bar plot on the top, where each bar indicates a firing time of a motor unit. The pop-up window and its associated functions allow the user to compare the signal data in that gap with the various templates to determine whether there is sufficient evidence to indicate that the program made a mistake. In this particular instance, the operator judged that the signal data matched the visual characteristics of the motor unit with the gap in its bar plot. Consequently, the user was able to modify the decomposition results as indicated by the same MU's bar plot at the bottom of Fig. 6. By using the interactive editor, it is possible to obtain decomposition accuracies that approach 100%.

View larger version (56K): [in this window] [in a new window]   FIG. 6. Interactive Editor program. Motor unit (MU) train at the top contains an apparent error—the gap in the oval. By examining the underlying 3-channel EMG data it is possible to identify the unlabeled waveform as belonging to MU 4. Interactive Editor facilitates the decision-making process by providing local firing statistics for each motor unit, the ability to overlap a motor unit template with the signal waveform, and graphical representation of the residual signal energy after subtraction of the template.

	RESULTS

TOP ABSTRACT INTRODUCTION BACKGROUND METHODS RESULTS DISCUSSION APPENDIX: SIGNAL PROCESSING... GRANTS ACKNOWLEDGMENTS REFERENCES

View larger version (29K): [in this window] [in a new window]   FIG. 7. Top: examples of the 4 channels of the sEMG signals detected by the surface array sensor from the tibialis anterior muscle during a 50% MVC contraction. Dotted vertical lines indicate the location of identified motor unit action potentials. Bottom: firing times of 6 motor units extracted from the decomposed sEMG signals above. Note that the location of the action potentials indicated by the bars corresponds to the identifiable pulses in each channel of the sEMG signals. Algorithms were able to identify the recruitment of MU 5 and MU 6, events that are not visually recognizable in the sEMG signals. Also, the algorithms resolved 7 cases of superpositions that are not obvious from the appearance of the sEMG signals.

The remainder of the decomposition of the trapezoid contraction is presented in Fig. 8. This is the result of the edited version, which required 12 h of expert user interaction. The automatic decomposition was obtained in 8 min with an accuracy of 90.4%. Figure 8A presents all the firings of six MUs throughout the complete duration of the 25-s-long contraction. (We refer to this as the Bar Plot.) The dark line represents the force produced by the subject. The arrows on the force axis indicate the recruitment threshold of the MUs. Figure 8B presents the interpulse interval, plotted vertically for each firing of all the MUs. (We refer to this as the Dot Plot.) This presentation is useful for visually identifying the errors in the decomposition, which are seen as dots that fall either noticeably below or above the average firing interval. Figure 8C presents the time-varying mean firing rates of each MU. The time-varying mean firing rates are obtained by passing the impulse trains corresponding to the firing trains of the top plot through a unit-area Hanning filter of 1.0-s duration. The decreasing firing rate value as a function of time during the constant force contraction is consistent with similar previous observations by De Luca et al. (1996). The shapes of the action potentials of the MUs for all four channels are presented in Fig. 8D. The shapes fluctuated throughout the contraction. The thick line represents the mean value of the amplitude of the action potentials and the thin, dashed lines above and below represent the SD.

View larger version (45K): [in this window] [in a new window]   FIG. 8. Results of a decomposition of the sEMG signal detected from the tibialis anterior muscle during a 50% MVC contraction. This automatic decomposition yielded an accuracy of 90.4% of the presented data, which was edited with Interactive Editor. A: bar plot of the individual firings of 6 concurrently active motor units. Continuous trace represents the isometric force generated by the subject. B: dot plot of the interpulse intervals of each motor unit. Vertical displacement indicates the time between adjacent firings. C: plot of the time-varying mean firing rates of each motor unit. D: shapes of the action potentials of 6 motor units as recovered from the raw sEMG signal through spike-triggered averaging. Mean (solid line) ± SD (dashed line) are plotted for each action potential.

View larger version (39K): [in this window] [in a new window]   FIG. 9. Results of sEMG decomposition from the orbicularis oculi and platysma. Data are examples of motor unit firing patterns during a 30 and a 50% MVC contraction from one subject. Plots on the left present the firing times of motor units; plots on the right represent the time-varying mean firing rates of the motor units. Solid black lines represent the root-mean-squared value of the EMG signal detected with standard surface sensors, which provides an indication of the force exerted by the muscle (see text for details).

Proof of accuracy and consistency

When one proposes to decompose a signal whose composition is not known a priori, it is incumbent on the proponent to prove that the decomposition is performed accurately. In our case it is also necessary to provide evidence that the rudimentary PD III decomposition technique provides information similar to that of the more established PD II technique, which decomposes the iEMG signal detected with indwelling sensors. One obvious approach for proving accuracy would be to construct mathematically a signal of known components (MU action potential trains) and then proceed to decompose the synthesized signal into the constituent action potentials and check the veracity of the outcome. Over two decades ago Mambrito and De Luca (1984) showed that if the shapes of the action potentials remain invariant throughout the contraction, this approach is not sufficiently challenging, even when white noise having an amplitude 20% of the average amplitude of the action potentials was added to the signal. The challenges presented by real signals and identified in the preceding section represent important complexities that must be dealt with. Additionally, the algorithms must be able to cope with other unforeseen signal components that are occasionally presented by real data as a result of abrupt displacements of the sensor with respect to the source of the action potentials. Consequently, we used the "two-source" technique that we first introduced in 1984.

We placed a surface sensor array (Fig. 2A) on the skin above the tibialis anterior muscle and a quadrifilar needle sensor was inserted directly beneath the array sensor, approximately 1 cm into the muscle. It was expected that some of the motor units would contribute signals to both sensors and others to only one sensor, depending on the proximity of the muscle fibers to the separate sensors. The EMG signal from each sensor was automatically decomposed by the respective decomposition algorithm and then independently further processed with the assistance of the Interactive Editor. The firing times of the motor units from both sensors are shown in Fig. 10. Note that three of the five MUs decomposed from the iEMG signal were also found in the sEMG. The expanded time sequence in the middle indicates that the MUs are firing at precisely the same time. This was so for 97.6%, or 996 sEMG-detected out of 1,021 iEMG-detected firings of the three MU action potential trains. This proves that the rudimentary PD III can obtain the same information from the sEMG signal as the PDII obtains from the iEMG signal. The precise coincidence of all the action potentials in three trains extracted from the two signals also proves that both the PD II and the PD III decomposition algorithms are able to accurately decompose the signals because the probability that both algorithms make exactly the same errors in both signals is extremely small.

View larger version (16K): [in this window] [in a new window]   FIG. 10. Comparison of decomposition results of simultaneously recorded indwelling EMG (iEMG) and sEMG signals. Open spaces on the sEMG decomposition indicate some of the errors made by the algorithms. Note the perfect match for the motor unit firings in the expanded time interval.

	DISCUSSION

TOP ABSTRACT INTRODUCTION BACKGROUND METHODS RESULTS DISCUSSION APPENDIX: SIGNAL PROCESSING... GRANTS ACKNOWLEDGMENTS REFERENCES

In addition to demonstrating a technical success, the decomposed data reveal two interesting physiological implications. The first concerns the behavior of the firing rates. It is noteworthy, that during a contraction the hierarchical organization of the firing rates of the orbicularis oculi and platysma, which are innervated by the VII cranial nerve, appears similar to that of the tibialis anterior, which is innervated by the L4, L5, and S1 spinal nerves. In particular, values of the mean firing rates of MUs are inversely proportional to their recruitment thresholds, such that earlier-recruited, low-threshold MUs maintain higher firing rates than later-recruited, higher-threshold units. (The occasional crossover between adjacent firing rates seen in Figs. 8 and 9 can be attributed to noise components in the common drive to the motoneuron pool.) The layered appearance of the plots of the firing rates versus force of successively recruited MUs led us to refer to this phenomenon as "onion skin" (De Luca and Erim 1994; De Luca et al. 1982a). This characteristic can be seen in data reported by other investigators who have signal processing knowledge (Freund et al. 1975; Hoffer et al. 1987; Masakado et al. 1995; McGill et al. 2005; Person and Kudina 1972; Rose and McGill 2001; Stashuk and De Bruin 1988). Some (Gydikov and Kosarov 1974; Moritz et al. 2005; Taylor and Enoka 2004) have argued against the hierarchical arrangement of the firing rates based on simulation studies or by combining the firing rate versus time curves of several subjects and/or contractions and emphasizing the lack of the hierarchical arrangement. Clearly the latter approach introduces intrasubject variations and cannot be used to describe the relative behavior of firing rates present within a contraction.

The second implication is that the dynamic range of the firing rate of MUs innervated from the cranial nerves (i.e., the orbicularis oculi and the platysma) is similar to that of the small muscles of the hand, which are innervated by the cervical nerves. De Luca et al. (1982a) reported that the firing rates of the first dorsal interosseous reached values of 41 pps; Kukulka and Clamann (1981) and Bigland-Ritchie et al. (1983) reported those of the adductor pollicis to reach 32 and 35 pps, respectively. These firing rate values are in contrast to those observed in large limb muscles, as exemplified by the tibialis anterior in this study and various other muscles in previously reported studies (Adam and De Luca 2005; De Luca and Erim 1994; De Luca et al. 1982a,b). In the two examples provided in Fig. 9, values of the firing rates of the orbicularis oculi and the platysma have a considerably greater spread among MUs of different recruitment threshold compared with the behavior of the MUs in the tibialis anterior (Fig. 8). Also the firing rates reach greater values in the orbicularis oculi and the platysma. For example, during the 50% MVC contractions, the firing rates of the platysma reach values of 35 pps, whereas those in the tibialis anterior reach about 17 pps. Although speculative, the higher dynamic range and absolute values of the firing rates of the orbicularis oculi and the platysma might be explained by a lack of recurrent inhibition by the Renshaw system. Its absence has been reported by Windhorst (1996) in the nearby oculomotor and masticatory muscles.

	APPENDIX: SIGNAL PROCESSING DETAILS

TOP ABSTRACT INTRODUCTION BACKGROUND METHODS RESULTS DISCUSSION APPENDIX: SIGNAL PROCESSING... GRANTS ACKNOWLEDGMENTS REFERENCES

 
Signal acquisition in PD III

• Number of signal channels: 4
  
• Sampling rate of input signal: 20,000 Hz
  
• Analog antialiasing filter: fourth-order Butterworth filter; 24 dB/octave slope; 9,500 Hz cutoff
  
• Analog high pass filter: second-order Butterworth filter; 12 dB/octave slope; 100 Hz cutoff
  

First stage of PD III

• Digital filter type: eighth-order Butterworth band-pass
  
• Digital filter's lower cutoff: 24 dB/octave <250 Hz
  
• Digital filter's upper cutoff: 24 dB/octave >2,000 Hz
  

Second stage of PD III

INPUT TO MAP RECEIVER. The input to the MAP receiver stage is a segmented version of x_B(t), the filtered input signal. Each signal segment is established by virtue of its signal amplitudes being, on average, above a threshold that is determined adaptively by the decomposition system in accordance with a set of dynamic range criteria for the amplitudes of the decomposable MU trains. The segmented signal is then used to initiate and update action potential templates for each hypothesized MU. These templates are, in turn, used by the MAP receiver to help classify signal components that are detected by amplitude peaks in the segmented signal. The task of the MAP receiver is to associate a hypothesized MU with each detected signal component.

BASICS OF MAP RECEIVER. Denoting the ith MU by the symbol u_i, the MAP receiver assigns a particular MU, say, u_k, to the signal component such that

where is the data vector for the detected pulse. In the presence of white Gaussian noise with variance ², the MAP decision criterion may be stated as

where s_i is the data vector for the template of the ith MU, and P(u_i) is the a priori probability that the detected pulse belongs to the pulse train u_i. This a priori probability is estimated by a hazard function calculation (LeFever and De Luca 1982), assuming that the first-order interpulse intervals have a Gaussian density for each MU.

TEMPLATE INITIATION. An EMG signal segment d is declared as an initiation of a new motor unit train if for each preexisting template s, the value that minimizes |d – s|² is such that |1 – | > 0.5 or |d – s|²/|s|² > 0.5.

TEMPLATE UPDATE. When the MAP receiver determines a match between segment data d and a preexisting template s(old), the segment data are used to update the template in accordance with the recursive relation: s(new) = s(old) + d. The initial template and all of its subsequent updates are saved for later use by the decomposition system.

Third stage of PD III

TEMPLATE SELECTION. When comparing two trains to determine whether they are both split off from a single train, the decomposition system selects the pair of templates (one from each train) that are closest in time to each other.

PROBABILITY ESTIMATION. In comparing template s_j with template s_k we estimate the probability that they are both from the same train as P_j,k = (s_j · s_k)/(|s_j||s_k|), where the · symbol denotes the dot-product operation; = for 0 1, 1/ for > 1, = 0 for < 0, and = (s_j · s_k)/(s_j · s_j).

Fourth stage of PD III

The basic strategy is that a representative template for each MU train is correlated against the entire filtered signal. The representative template is selected as the one that was most recently updated. The goal is to determine locations where the "maximal amplitude" MU, i.e., the one whose representative template attains the greatest amplitude from among all the motor units, may be present. To achieve this goal, the receiver first determines the locations where the correlation function of the maximal MU has a local maximum that is greater than the correlation values at that point of all the other MUs. At each of the determined locations, an estimate of the probability that the maximal MU's action potential actually occurred at the point is obtained as

where the template for the MU is represented by the vector p, the corresponding signal data are represented by the vector d, is the scale factor that minimizes the value of |d – p|², and = if 0 1; = 1/ if > 1; and = 0 if < 0.

Conceptually, represents the degree to which d and p are co-linear and (d – p) represents the orthogonal component of the modeling error. A probability threshold is also established by the decomposition system as the largest probability value obtained by correlating the maximal-amplitude MU's template with the templates of the remaining MUs. At all the locations where the estimated probability is above the probability threshold, the template of the maximal MU (appropriately scaled) is subtracted from the sEMG signal to obtain the input signal for the next iteration. The next iteration is the same as the previous one except for the removal from consideration of the maximal MU of the previous iteration, thus leading to the declaration of another MU as the maximal MU for that iteration. Once all the iterations are complete, a utility maximization procedure is used for converting probabilities of occurrence into decisions about which MU may be actually responsible for which action potential contribution in the signal data.

	GRANTS

TOP ABSTRACT INTRODUCTION BACKGROUND METHODS RESULTS DISCUSSION APPENDIX: SIGNAL PROCESSING... GRANTS ACKNOWLEDGMENTS REFERENCES

 
This work was supported by National Institutes of Health Bioengineering Research Partnership Grant 1R24HD-38585 from the National Center for Medical Rehabilitation Research of the National Institute of Child Health and Human Development, the National Aeronautics and Space Administration Grant NCC 9–127 and Delsys Inc.

	ACKNOWLEDGMENTS

TOP ABSTRACT INTRODUCTION BACKGROUND METHODS RESULTS DISCUSSION APPENDIX: SIGNAL PROCESSING... GRANTS ACKNOWLEDGMENTS REFERENCES

 
The authors thank A. Morgan and A. Trongnetrpunya for assistance in processing the data.

	FOOTNOTES

 
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

Address for reprint requests and other correspondence: C. J. De Luca, NeuroMuscular Research Center, 19 Deerfield Street, Boston, MA 02215 (E-mail: cjd{at}bu.edu)

	REFERENCES

TOP ABSTRACT INTRODUCTION BACKGROUND METHODS RESULTS DISCUSSION APPENDIX: SIGNAL PROCESSING... GRANTS ACKNOWLEDGMENTS REFERENCES

 
Adam A and De Luca CJ. Recruitment order of motor units in human vastus lateralis muscle is maintained during fatiguing contractions. J Neurophysiol 90: 2919–2927, 2003.[Abstract/Free Full Text]

Adam A and De Luca CJ. Firing rates of motor units in human vastus lateralis muscle during fatiguing isometric contractions. J Appl Physiol 99: 268–280, 2005.[Abstract/Free Full Text]

Adam A, De Luca CJ, and Erim Z. Hand dominance and motor unit firing behavior. J Neurophysiol 80: 1373–1382, 1998.[Abstract/Free Full Text]

Andreassen S. Computerized analysis of motor unit firing. In: Computer-Aided Electromyography, edited by Desmedt JE. Basel: Karger, 1983, p. 150–163.

Basmajian J and De Luca CJ. Muscles Alive (5th ed.). Baltimore, MD: Williams & Wilkins, 1985.

Bigland-Ritchie B, Johansson R, Lippold OCJ, Smith S, and Woods JJ. Changes in motoneuron firing rates during sustained maximal voluntary contractions. J Physiol 340: 335–346, 1983.[Abstract/Free Full Text]

Broman H. Knowledge-based signal processing in the decomposition of myoelectric signals. IEEE Trans Eng Med Biol 7: 24–28, 1988.[CrossRef]

Castañon DA. Efficient algorithms for finding the K best paths through a trellis. IEEE Trans Aerospace Electronic Syst 26: 405–410, 1990.[CrossRef]

Christodoulou CI and Pattichis CS. Unsupervised pattern recognition for the classification of EMG signals. IEEE Trans Biomed Eng 46: 169–178, 1999.[CrossRef][ISI][Medline]

De Figueiredo RJP and Gerber A. Separation of superimposed signals by a cross-correlation method. IEEE Trans Acoust Speech Signal Process 31: 1084–1089, 1983.[CrossRef]

De Luca CJ and Adam A. Decomposition and analysis of intramuscular electromyographic signals. In: Modern Techniques in Neuroscience Research, edited by Windhorst U and Johansson H. Heidelberg, Germany: Springer-Verlag, 1999, p. 757–776.

De Luca CJ and Erim Z. Common drive of motor units in regulation of muscle force. Trends Neurosci 17: 299–305, 1994.[CrossRef][ISI][Medline]

De Luca CJ, Foley PJ, and Erim Z. Motor unit control properties in constant-force isometric contractions. J Neurophysiol 76: 1503–1516, 1996.[Abstract/Free Full Text]

De Luca CJ and Forrest WJ. An electrode for recording single motor unit activity during strong muscle contractions. IEEE Trans Biomed Eng 19: 367–372, 1972.[Medline]

De Luca CJ, LeFever RS, McCue MP, and Xenakis AP. Behaviour of human motor units in different muscles during linearly varying contractions. J Physiol 329: 113–128, 1982a.[Abstract/Free Full Text]

De Luca CJ, LeFever RS, McCue MP, and Xenakis AP. Control scheme governing concurrently active human motor units during voluntary contractions. J Physiol 329: 129–142, 1982b.

Erim Z, Beg MF, Burke D, and De Luca CJ. Effects of aging on motor-unit control properties. J Neurophysiol 82: 2081–2091, 1999.[Abstract/Free Full Text]

Fang J, Agarwal GC, and Shahani BT. Decomposition of multi-unit electromyographic signals. IEEE Trans Biomed Eng 46: 685–697, 1999.[CrossRef][ISI][Medline]

Finsen F, Søgaard K, Graven-Nielsen T, and Christensen H. Activity patterns of wrist extensor muscles during wrist extensions and deviations. Muscle Nerve 31: 242–251, 2005.[CrossRef][ISI][Medline]

Freund HJ, Budingen HJ, and Dietz V. Activity of single motor units from human forearm muscles during voluntary isometric contractions. J Neurophysiol 38: 933–946, 1975.[Abstract/Free Full Text]

Garcia GA, Okuno R, and Akazawa K. Technological developments in Japan—a decomposition algorithm for surface electrode-array electromyogram: a noninvasive, three-step approach to analyze surface EMG signals. Eng Med Biol Mag IEEE 24: 63–72, 2005.[CrossRef]

Gazzoni M, Farina D, and Merletti R. A new method for the extraction and classification of single motor unit action potentials from surface EMG signals. J Neurosci Methods 136: 165–177, 2004.[CrossRef][ISI][Medline]

Gerstein GL and Clark WA. Simultaneous studies of firing patterns in several neurons. Science 143: 1325–1327, 1964.[Abstract/Free Full Text]

Glaser EM and Marks WB. The on-line separation of inter-leaved neuronal pulse sequences. Rochester Conf. on Data Acquisition in Biology and Medicine, 1966, p. 137–156.

Guiheneuc P, Doncarli C, and Le Bastard M. Concomitant multitrain automatic analysis of motor unit potentials and macroEMG. In: Computer-Aided Electromyography and Expert Systems, edited by Desmedt JE. Amsterdam: Elsevier, 1989, p. 83–90.

Gydikov A and Kosarov D. Some features of different motor units in human biceps brachii. Pfluegers Arch 347: 75–88, 1974.[CrossRef][ISI][Medline]

Hochstein L, Nawab SH, and Wotiz R. An AI-based software architecture for a biomedical application. SCI-2002. Proc 6th World Multiconf Systemics Cybernetics, Informatics, Orlando, FL, July 2002, vol. XI, p. 60–64.

Hoffer JA, Sugano N, Loeb GE, Marks WB, O'Donovan MJ, and Pratt CA. Cat hindlimb motoneurons during locomotion. II. Normal activity patterns. J Neurophysiol 57: 530–553, 1987.[Abstract/Free Full Text]

Holobar A and Zazula D. Correlation-based decomposition of surface electromyograms at low contraction forces. Med Biol Eng Comput 42: 487–495, 2004.[CrossRef][ISI][Medline]

Iani C, Stalberg E, Falck B, and Bishof