PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Neurosci Methods. Author manuscript; available in PMC 2010 December 14.
Published in final edited form as:
PMCID: PMC2904617
NIHMSID: NIHMS209682

Automatic classification of motor unit potentials in surface EMG recorded from thenar muscles paralyzed by spinal cord injury

Jeffrey Winslow, PhD,a,b Marina Dididze, MD, PhD,a and Christine K Thomas, PhDa,c

Abstract

Involuntary electromyographic (EMG) activity has only been analyzed in the paralyzed thenar muscles of spinal cord injured (SCI) subjects for several minutes. It is unknown if this motor unit activity is ongoing. Longer duration EMG recordings can investigate the biological significance of this activity. Since no software is currently capable of classifying 24 hours of EMG data at a single motor unit level, the goal of this research was to devise an algorithm that would automatically classify motor unit potentials by tracking the firing behavior of motor units over 24-hours. Two-channels of thenar muscle surface EMG were recorded over 24-hours from 7 SCI subjects with a chronic cervical level injury using a custom data logging device with custom software. The automatic motor unit classification algorithm developed here employed multiple passes through these 24-hour EMG recordings to segment, cluster, form global templates and classify motor unit potentials, including superimposed potentials. The classification algorithm was able to track an average of 19 global classes in 7 24-hour recordings with a mean (± SE) accuracy of 89.9 % (± 0.98%) and classify potentials from these individual motor units with a mean accuracy of 90.3% (± 0.97%). The algorithm could analyze 24 hours of data in 2–3 weeks with minimal input from a person, while a human operator was estimated to take more than 2 years. This automatic method could be applied clinically to investigate the fasciculation potentials often found in motoneuron disorders such as amyotrophic lateral sclerosis.

Keywords: EMG, long term recording, motor unit classification, spinal cord injury

1.0 Introduction

Spinal cord injured (SCI) individuals often have no voluntary control over skeletal muscles that are innervated from spinal segments below the lesion. However, these paralyzed muscles can contract involuntarily (spasm) in response to various stimuli such as vibration, cutaneous inputs, and changes in temperature and body position (Kawamura et al., 1989; Little et al., 1989; Thomas and Ross, 1997). In other situations, weak involuntary contractions have been recorded for up to 30 minutes, particularly in hand muscles (Stein et al. 1990; Zijdewind and Thomas 2001; Zijdewind et al. 2004). Whether this spontaneous motor unit activity takes place over extended time periods is unknown. This activity appears to be spontaneous in that no obvious stimulus generates the motor unit activity. These different types of muscle activity can be recorded as electromyographic activity (EMG).

Long-term EMG recordings, for example over 24 hours, are important because they can provide detailed information about motor unit activity and muscle behavior. For example, they can be used to investigate the concept that only small amounts of neuromuscular activity are required to sustain muscle and motor unit properties (Kernell et al., 1987; Kim et al., 2007). They could also be used to evaluate whether medications dampen activity. Since the orderly recruitment of motoneurons will largely determine which neurons are activated, and medications may only influence certain motoneurons, it is important to examine these kinds of questions at the motor unit level.

Spontaneous EMG recorded from thenar muscles paralyzed by spinal cord injury largely involve single motor unit activity. Thus, these data provide a unique opportunity to non-invasively investigate motor unit behavior using surface electrodes. To analyze this single motor unit behavior, it is critical to identify the action potentials that belong to a single motor unit and to follow that unit activity over time. A computerized method is required to accomplish this analysis because it is not practical to do these tasks manually (Zennaro et al, 2003). This process is known as EMG decomposition (Mambrito and De Luca, 1984). Many algorithms have been derived to decompose EMG arising from both surface and intramuscular recordings (Florestal et al., 2009), but there have been no previous attempts to classify 24-hours of EMG.

Many investigators have constructed software packages with various levels of user interaction capable of decomposing intramuscular EMG from able-bodied subjects into its constituent motor unit potentials. The Zoom (Edin et al, 1988) and Spike 2 (Cambridge Electronic Design, Cambridge, England) software packages can classify individual motor unit potentials when the user manually sets specific waveform amplitude and duration features for which to search, but both require extensive operator assistance. McGill et al. (2005) composed EMGLAB, an interactive EMG decomposition tool written in Matlab that is template-based and classifies one or more EMG channels by using a novel similarity metric (Florestal et al., 2006). Zennaro et al.’s (2003) EMG-LODEC program uses wavelet-based features and a template-based method to classify motor unit potentials in longer recordings ranging from three minutes to several hours.

To aid the decomposition of surface EMG into motor unit potentials, investigators have used multi-channel arrays to provide spatial information about the motor unit (De Luca et al., 2006), sometimes in the form of bipolar and Laplacian montages (Kleine et al., 2007). Algorithms have used hierarchical clustering (Kleine et al. 2007), a fuzzy rule set based on differences in firing rates, peak to peak values of motor unit potentials, and motor unit based waveform correlations (Chauvet et al., 2003) or convolution kernel compensation (CKC), a non-template based method of classification (Holobar et al., 2009).

In general, decomposition of surface EMG into single motor unit potentials has largely been completed for weak, isometric voluntary contractions of limited duration. The relative stability of this recording situation means that classification algorithms have been able to depend on the shapes of motor unit waveforms and relatively well documented patterns of motor unit activity to identify the potentials that belong to different single motor units. Manual editing has also been used to improve classification accuracy. However, quantification of algorithm performance has been limited because user-defined classification of motor unit potentials is time consuming.

The aim of our study was to devise a multi-level, multi-step template-based classification procedure to automatically classify and analyze involuntary single motor unit activity in EMG recordings that are 24 hours in length. Classification of motor unit potentials in long-term surface EMG recordings during involuntary contractions presents different issues and challenges that existing classification software packages are ill-equipped to handle. For example, spontaneous motor unit activity that occurs without any obvious trigger may involve unpredictable changes in contraction strength, muscle length, and motor unit firing rates. The limited recordings of spontaneous motor unit activity from paralyzed muscles show that sustained motor unit activity can involve regular or highly variable firing rates. Sporadic potentials are common, as are doublets (Thomas and Ross 1997; Zidewind and Thomas, 2001; Zijdewind et al. 2001). Until the firing behavior of these spontaneously active motor units is analyzed more extensively, it is difficult to use low or high frequency cutoffs, activity patterns, or firing variability to constrain an algorithm. Thus, at present, an algorithm that strives to classify spontaneous motor unit activity accurately has to rely largely on motor unit potential shape. Moreover, to capture a representative picture of slow, irregular motor unit activity, like the fasciculation potentials often seen in muscles of people with motoneuron disorders (Drost et al. 2007), it is necessary to make long recordings.

Since short-term surface EMG recordings from thenar muscles paralyzed chronically by cervical SCI show a prevalence of spontaneous EMG, we recorded data from these muscles to address some of the issues presented by long-term involuntary muscle contractions. To provide more spatial information for classification of any given motor unit potential, two channels of EMG were recorded simultaneously using surface electrodes. The duration of the recordings also required that the analysis be largely automatic and time efficient. Thus, important steps like resolution of superpositions were constrained to preserve accuracy and time. For the data to provide meaningful information about involuntary motor unit activity in muscles paralyzed by SCI, classification accuracy is imperative. Thus, we employed a multiple-pass segmentation and classification method, rather than a sequential technique, to use a broad view of the data for decision making and robustly track potentials over 24 hours. Furthermore, both the dataset and the combination of methods that form the overall algorithm are novel. Thus, extensive manual analysis was used to evaluate the performance of the automated procedure. Comparisons were made between the results from the algorithm and a person, two people, the two surface EMG sources, as well as intramuscular and surface EMG.

2.0 Materials and methods

2.1 Subject pool

The subject pool consisted of 7 subjects, aged 27–60, with chronic (>1 yr) cervical level injury at C4 (n=2), C5 (n=2), or C6 (n=3). Six subjects were male while one was female. Three subjects (subjects 3, 6 and 7) took no medication to dampen muscle spasms. The thenar muscles (abductor pollicis brevis, flexor pollicis brevis, and opponens pollicis) were not under voluntary control in either hand as confirmed by manual muscle examination (Thomas, 1997). All of the procedures were approved by the University of Miami Institutional Review Board. All subjects signed consent forms prior to participation. One 24-hour recording was made per subject.

2.2 Recording Setup

Two channels of EMG were recorded from the thenar muscles in order to provide additional motor unit potential features to aid in classification. For a 2-channel recording, five self-adhesive electrodes (Superior Silver electrodes, Uni-Patch, MN) were trimmed to approximately 1 cm × 2.5 cm using scissors so that they could be positioned across the thenar muscles. The distal electrode was positioned across the metacarpal-phalangeal joint, the proximal electrode at the base of the thenar muscles, and two closely spaced electrodes across the muscle bellies at their midpoint (Fig. 1). Another electrode was placed just proximal to the wrist as a ground. The electrodes were then secured to the thumb skin overlying the thenar muscles using a layer of hypafix medical tape, followed by several wrappings of Co-Flex adhesive bandage (Andover Healthcare Inc., Salisbury, MA) to make sure that the electrodes maintained the same position on the skin during the EMG recording. The median nerve was stimulated to evoke maximal compound action potentials (M-waves) both immediately before and after 24-hour recordings. The peak-to-peak amplitude of the distal (and proximal) pre-and post-M-waves were compared in order to objectively measure electrode stability. Neither the distal (mean ± SE pre: 8.4 ± 1.3 mV, post: 8.8 ± 1.3 mV, p=0.12) nor the proximal amplitudes (pre: 6.0 ± 1.2 mV, post: 6.2 ± 1.2 mV, p=0.46) differed significantly across time (paired t-tests).

Fig. 1
Electrode configurations for recording surface EMG from the thenar muscles of the right hand. A wrist electrode was utilized as the common for both preamplifier channels. Two channels of thenar EMG were recorded, one more distally located (distal channel) ...

The electrodes were connected to the inputs of two differential preamplifiers (manufactured by Motion Control, Salt Lake City, UT) that had gains of approximately 400. The preamplifier outputs were fed to a preprocessing unit that was responsible for scaling the signal to an approximate 2 V baseline and further amplifying the signal to fill the input range of the data logger. The preprocessing unit was equipped with analog bandpass filters with cutoff frequencies of 10 Hz and 1000 Hz. A portable, battery powered data logging device recorded the output of the preprocessing unit. The data logging device consisted of the Tattletale 8 Logger, manufactured by Onset Computer Corporation, driven by specially designed software. It had a 12-bit analog-to-digital converter that accepted input signals from 0 to 4.096 V. Data was stored on a 1.0 Gigabyte type I compact flashcard for later analysis. The analog-to-digital converter used a sampling frequency of 3.2 KHz, the same sampling frequency used in the laboratory. The data logger was calibrated prior to and after the 24 hour recording with a 100 Hz, 1 mV peak to peak sine wave obtained from a signal generator.

2.3 Multi-level, multi-step classification algorithm

All data were processed using Matlab 7.0 and 7.4 (The Mathworks, Natick, Massachusetts), a programming suite whose capability ranges from simple processing tasks to complex computations. The 24-hour recording was broken into 480, 3-minute portions for each channel. All processing was done independently for each channel. Three minutes of data were found to offer the best recording length for clustering of potentials and provided feasible lengths for manual classification. Data were further processed by applying a zero phase 2nd order 60 Hz notch filter, a zero phase 2nd order bandpass filter with a passband of 30 to 500 Hz, and a simple moving average filter.

The motor unit action potential classification algorithm utilized two channels of EMG (distal and proximal) and consisted of four processing passes through the 24-hour recording (Fig. 2). These passes included: 1) segmentation of motor unit potentials. This involved threshold determination, in which a global (24-hour) threshold was computed in order to isolate (segment) individual motor unit potentials from the baseline; 2) clustering, in which similar potentials were grouped into their respective motor unit classes. Representative templates were also calculated for each class of potential; 3) global template formation where the most stable, commonly occurring motor unit classes were found throughout the 24-hour recording; 4) final classification, where all unclassified potentials were compared to the global motor unit templates in order to achieve a classification with or without the need for superposition resolution. The subsequent sections describe these four passes in greater detail.

Fig. 2
Block diagram for motor unit potential classification. The classification algorithm consisted of 4 major passes: 1) segmentation, in which individual potentials were isolated from the continuous EMG; 2) clustering, in which similar potentials were grouped ...

2.3.1 Segmentation

Segmentation represents the method of breaking up the recording into individual motor unit potentials. This is done by extracting the data points that comprise a potential in a window for later classification. Segmentation occurred in three phases and on both the distal and proximal channels independently. The three steps involved estimation of the baseline noise, determination of the minimal signal required for segmentation of potentials, followed by segmentation. In the first phase (Fig. 3), baseline signal estimates for each of the 480 3-minute records were computed. The signal was windowed every 35 data points (the approximate duration of two phases of a potential) to produce peak to peak values per window. The minimum values taken every 150 data points (a regional EMG baseline estimate from the approximate duration of 4 potentials) were averaged to represent the baseline signal for a given 3-minute record. Thus, there were 480 different baseline signals, one for each 3-minute record. In the second phase, a global segmentation amplitude threshold was determined. This threshold represented the quiet signal or the smallest amplitude at which motor unit potentials could be observed with respect to the noise. The means and standard deviations of all 480 of the 3-minute baseline signal estimates were computed and the standard deviations were ordered from least to greatest. The 3-minute record whose standard deviation was twice that of the overall smallest standard deviation was utilized as the standard deviation of the baseline for the entire recording. The corresponding mean for this selected 3-minute record was used as the mean baseline for the entire recording. The doubling of the baseline variability with respect to its minimum value during the recording served as an estimate of the overall baseline. This estimate performed well in all seven recordings that were analyzed, providing good baseline estimates. Finally, the global segmentation amplitude threshold was the mean baseline plus 2 standard deviations. The third and final phase actually segmented candidate motor unit potentials using the global segmentation threshold. To ensure that the alignment of the windows used to compute peak to peak values did not affect segmentation, the original EMG signal was windowed again. Peak to peak values were computed just as in phase one, but this time it was repeated at an offset of half a window length. All regions where peak to peak values exceeded the global segmentation threshold were flagged and the minimum of the region was taken as the fiducial alignment point. Individual motor unit potentials were represented by 35 data points to either side of the fiducial alignment point, equivalent to the approximate duration of 2 motor unit potentials, enabling the resolution of any superpositions.

Fig. 3
Segmentation phase one. A peak to peak filter determined the peak to peak EMG value (dotted trace) in consecutive 32 point windows (10 ms). The baseline signal estimate (black trace) was computed by finding the minimum value in fixed 150 point (47 ms) ...

A noise rejection method complemented the global segmentation amplitude threshold in preventing noise from being segmented as candidate potentials. Noise was eliminated based on frequency content and was distinguishable from motor unit potentials based on the fact that actual potentials had most spectral content at the lower frequencies.

Muscle spasms (large EMG bursts) were also problematic and contaminated the spontaneous single motor unit activity that was to be classified. The automated classification algorithm was only designed to classify discrete visible potentials and superpositions consisting of a small number of motor unit potentials. Spasms were manually labeled by an expert classifier and were excluded from all results. On average, 25.4 minutes (range: 1.6–60.1 min) were removed from the recordings due to spasms.

2.3.2 Clustering and representative template formation

Transitive clustering, a novel clustering technique based on the transitive property, was used to cluster individual potentials independently for each channel and in each 3-minute EMG record. The Euclidean distances (the vector version of the simple Euclidean distance formula) between all possible combinations of all potentials were computed along with their corresponding correlation coefficients. The correlation coefficient is simply the computed covariance between the potentials divided by the product of their standard deviations. Potentials were paired if the Euclidean distance between them was less than 0.25 and their correlation coefficient was greater than 0.95. These thresholds were found to be the optimum values to achieve the most compact and separate clustering. Transitive clustering performs according to the following example with 3 different potentials labeled A, B, and C. If potentials A and B were paired and potentials B and C were also paired, then potentials A, B, and C should be clustered according to the transitive property. A cluster had to have at least 3 member potentials per 3-minute record to be retained. This limit insured that true, re-occurring potentials would be followed.

Since both channels were clustered independently and were given different labels, clusters were then matched between channels to yield the final motor unit classes in each 3-minute EMG record (Fig. 4). Since thenar fibers probably go from one end of the muscle to the other (Westling et al., 1990), and the EMG electrodes lay across the bellies of the muscles in a distal/proximal configuration, motor unit potentials were recorded simultaneously in the two channels. Thus, each channel had a different representation of the same single motor unit potential, additional information used to ensure accurate clustering of potentials from any given motor unit. Even though most clusters in one channel had a direct match with a corresponding cluster in the opposite channel, those with no direct matches were still retained. New corresponding clusters were formed in the channel with no direct matches, thereby increasing the number of clusters found by the algorithm. Template potentials in each channel were computed by taking the mean of all member potentials in a given cluster. Each cluster was represented by two template potentials for subsequent classification procedures, one in each channel.

Fig. 4
Cross channel template matching. Each channel was clustered separately, yielding different labels (Clusters A and B in the proximal channel and Clusters X, Y, and Z in the distal channel). Each potential has a unique representation in each channel, and ...

2.2.3 Global template formation

Composite templates were formed for each cluster in each 3-minute record by joining the motor unit class templates from the distal and proximal channels to form one template of double length. Thus two unique signatures of each cluster were combined. The composite templates were then clustered transitively across 3-minute records to form the global motor unit classes to be tracked throughout the 24-hour recording. The term global motor unit class refers to a motor unit whose activity has been tracked by the classification algorithm. Since many motor units may be firing during a 24-hour recording only the most active motor units were tracked in each 24-hour recording, while those that were only occasionally active were ignored. For a motor unit to be tracked (referred to as a global motor unit class), it had to consist of at least five member templates, or be present in 15 minutes of the recording (3 potentials per 3-minute recording in 5 3-minute recordings).

After transitive clustering, there were some global motor unit classes that required merging or elimination. Since high correlation and Euclidean distance thresholds were utilized to avoid tracking mistakes, a graphical user interface (GUI) was provided for a human operator to perform the required merging. Some global motor unit classes were in fact noise because some noise was periodic and clustered as a legitimate class. A human operator also removed these noise classes using the GUI.

There were some remaining template potentials throughout the 24-hour recording that were unclassified after the global classes had been formed. Attempts to assign global template labels to these unclassified templates were made by means of fuzzy membership values applied in a nearest neighbor fashion. The representative global motor unit templates for each global motor unit class were those classified composite templates that occurred closest to the 3-minute EMG record in question. Fuzzy membership values were used to compare composite member potentials from each local motor unit class to the representative templates from each global class. Fuzzy membership functions were utilized at this stage, rather than the Euclidean distance metric, because they provided both flexibility in differentiating between similar shapes and a robust normalization factor that standardized the “fuzzy” distances between the target motor unit potentials and template shapes as given by the summation in the denominator of equation (2).

Fuzzy membership values were computed using equations (1) and (2).

d2(XjVi)=n=1M(Xj(n)Vi(n))2
(1)

μij=[1d2(Xj,Vi)]1q1i=1K[1d2(Xj,Vi)]1q1
(2)

Equations (1) and (2) were taken from the fuzzy clustering procedure described by Zouridakis and Tam (2000). In equation (1), X is the vector representing an individual potential while V is the vector representing the nearest neighbor template. In the denominator of equation (2), K represents the total number of global templates. The quantity q is a fuzzy parameter that controls the fuzziness of the membership values. Values of q closer to 1 make the fuzzy membership values either 0 or 1, while larger values of q make the values more distributed in-between 0 and 1. A q of 1.8 was used in this case because other authors found success using this value (Zouridakis and Tam, 2000).

In order for a representative template from a 3-minute record to be assigned a global motor unit class label, at least 2 of its member potentials required fuzzy membership values larger than 0.7. This threshold was similarly found to be optimal by manually tracking global classes in each 24-hour recording. An example of the templates for the global classes from the subject 4 recording is shown in Fig. 5. Thus, all these steps ensured that the most commonly occurring motor units were followed over time and that they had templates that were stable over time.

Fig. 5
Global classes in one 24 hr recording. Eight global motor unit classes, each with a distal and proximal channel template, were found in the subject 4 recording during the global template formation stage of the overall classification algorithm.

2.3.4 Final classification

The aim of the final stage of the algorithm was to classify all individual potentials that remained unclassified. All unclassified potentials in each 3-minute EMG record were compared with nearest neighbor global motor unit class templates in order to obtain a classification. Individual potentials had to meet Euclidean distance (0.5) and correlation coefficient (0.94) thresholds when compared with the target global motor unit class templates in both channels. These threshold values were found to optimally reduce the number of misclassifications while maximizing those correctly classified. If a potential was still unclassified after this comparison, the superposition resolution routine attempted to resolve the potential into its constituent global motor unit classes. If there was still no classification after all of these steps, then the process was repeated after applying an attenuation windowing technique. In attenuation windowing, a Kaiser window was multiplied by the unclassified potential to attenuate any noise or unknown potentials at the periphery of the individual potential’s window.

In order to resolve the superimposed potentials into their constituent global motor unit classes, a modified dual channel peel-off method (Stashuk, 2001) was employed. Cross-correlation was used to find the correct time shifts at which to add motor unit template class potentials in order to best reconstruct the target superimposed potential. The correlation coefficients were then computed between each target potential and each template in both channels using these time shifts. The highest correlating global motor unit class template overall was selected as the constituent template irrespective of the channel. The chosen constituent global motor unit class template was then subtracted from the target potential, with the residual being retained. The entire process was repeated using the residual to build a superimposed potential to match the target potential. This process was stopped if the constructed superimposed potential exceeded the Euclidean distance and correlation thresholds in both channels or after a total of three constituent global motor unit classes were attempted and no match was found. This upper constituent limit was enforced because any waveform can be constructed with a large number of template shapes and misclassifications could result.

2.4 Performance Evaluation

The classification algorithm was evaluated in 3 ways. The segmentation performance (the ability of the algorithm to isolate potentials correctly and minimize the extraction of noise) was determined by comparing the potentials that were segmented by the algorithm against those segmented by an expert, considered the gold standard. The segmentation performance was based on five 3-minute records from one recording.

Classification performance (the ability of the algorithm to correctly classify individual motor unit potentials) was primarily based on 12 minutes of EMG data manually classified by an expert in each of seven recordings for a total of 84 minutes of data. When there were less than 10 member potentials for a given global motor unit class in the 12 minutes of data, additional 3-minute records were analyzed to provide at least 20 potentials for any given global motor unit class. A second person independently classified the potentials for 7 motor unit classes in the same 12 minutes of data. Intra-class correlations (ICC) were calculated in SPSS to determine algorithm-to-person and person-to-person reliability. Paired t-tests were used to evaluate the classification agreement between the algorithm and a person versus the agreement between two people. Classification accuracy was also determined using the two-source method (Mambrito and DeLuca 1984). Two motor unit classes from subject 1 were classified automatically and manually when intramuscular EMG and surface EMG were recorded simultaneously for 6 min. Furthermore, the importance of sampling rate was evaluated by reclassifying 25 motor unit classes from 5 subjects in 6 3-min recordings after halving the sampling rate (1.6 KHz). The classification accuracy was computed using gold standards derived from manual classification. The importance of using two channels of EMG for motor unit identification was evaluated by investigating the motor units that were identified in common between channels within the manually classified EMG data and those motor units that were uniquely found by only one channel.

The tracking performance was evaluated in two ways: 1) Tracking ability (the number of motor units the classification algorithm was able to track out of the estimated total number of active motor units during each recording), and 2) tracking accuracy (how well the algorithm could correctly classify representative templates as global motor unit classes over time), given that there could be in excess of 100 active units during the recordings (Yang and Stein, 1990) and the total number of active motor units actually recorded could only be estimated. The tracking accuracy was evaluated by comparing the algorithm’s template classifications over the entire 24-hour recording with those of an expert obtained through manual classification.

3.0 Results

3.1 Segmentation performance

The segmentation routine was evaluated by comparing the potentials that were segmented by the algorithm against those segmented by an expert in EMG recognition. In five 3-minute sections from one 24-hour recording containing 11,158 potentials, 93.7% of the motor unit potentials (n=10,455) were correctly segmented by the algorithm and 6.3% of potentials were missed (n=699). Only 250 potentials (2.3%) in these five 3-minute sections were false positives, or noise. This low false positive rate demonstrates the success of the noise rejection methods. All of the potentials that were missed were actually of too small a magnitude (less than 15 µV) to be adequately classified by the algorithm, or a person. If the segmentation performance were evaluated just on the basis of correctly finding all classifiable potentials, then all of the potentials were correctly segmented.

3.2 Classification performance

The classification performance of the algorithm was evaluated by comparing the results of the algorithm to the classification results of an expert, achieved manually by using a custom designed GUI-based interface. The number of potentials evaluated manually for each tracked unit ranged from 20, for the most sporadically active unit, to 2420, for the most frequently firing motor unit. Classification performance for each global motor unit class was based on 5 parameters: 1) correct classification (the percentage of manually classified potentials that were correctly labeled by the classification algorithm); 2) missed potentials (the percentage of the manually classified potentials that were left unclassified by the algorithm); 3) misclassified potentials (the percentage of potentials that were assigned an incorrect global motor unit class label); 4) false positives (the percentage of potentials that were assigned a global motor unit class label when they instead should have been left unclassified); 5) accuracy (Acc; the percentage of potentials that were correctly classified with respect to all potentials assigned a given motor unit class label using equation (3)).

%Acc.=ncorrectncorrect+nmissed+nmisclassified+nfalsepositives
(3)

The overall performance statistics are shown in Table 1. Of all 132 global motor unit classes identified over 7 recordings, the percent accuracy ranged from 38.1 to 100.0 % (mean ± SE, 90.3 ± 0.97%). More than three-quarters (n=106) of all global classes tracked had accuracies greater than 85% shown as a cumulative distribution for all global motor unit classes tracked throughout all recordings (Fig. 6) and by subject (Fig. 7). Only a small percentage of potentials were missed, misclassified, or false positives.

Fig. 6
Cumulative distribution of classification performance. The accuracy for all global motor unit classes for all recordings were sorted from low to high then plotted. This order was then applied when plotting the false positive or misclassified potentials. ...
Fig. 7
Percent accuracy with and without superpositions. The data for each global class of each recording are shown with different symbols against a line of unity (black line).
Table 1
Individual classification performance (n=132 global motor unit classes)

Two different expert classifiers created 12 min gold standard classifications for 7 global motor unit classes in the subject 2, and 3 recordings. On average, the two expert classifiers agreed on classifications 95% of the time. The automatic classification algorithm classified these same classes with an accuracy of 92 %. Algorithm-to-person (ICC: 0.998, α=0.999) and person-to-person agreement (ICC: 0.999, α=1.0) were high and not differ significantly (p=0.313).

The results were further validated by a two sources approach. The surface and intramuscular EMG recordings were in 93.9% agreement for two motor units from subject 1. The algorithm accuracy with respect to the expert classifier’s gold standard was 96.8%.

Potentials that were superpositions were more difficult for the algorithm to classify accurately. The classification accuracy was better for almost all global motor unit classes when superpositions were excluded (Fig. 7; data above the unity line) thus explaining some of the reduction in performance. Overall, the average reduction in accuracy when superpositions were included was 5.1%, ranging in-between −49.0% to 40.5% (standard error = 0.08%). In general, recordings with a higher number of global motor unit classes (Subjects 1, 3 and 7) showed greater increases in accuracy when superpositions were removed because the presence of more classes meant a better opportunity for superpositions.

The number of motor unit classes that could be identified at 3200 kHz was reduced to 58 % by halving the sampling rate (from 43 to 25 units). The potentials with peak to peak amplitudes of less than 50 µV were not stable enough to be clustered and identified at the 1600 Hz sampling rate. The classification accuracy was not significantly different (p = 0.1076) for the 25 motor units identified at both sampling rates.

To examine the importance of a two-channel recording in unit identification prior to global class formation, the number of motor units found in the distal and proximal channels separately, and in both channels, was examined over 32 3-minute records. Only 100 motor units were found in both channels. 168 motor units were found only in the distal channel, and 137 units were unique to the proximal channel (n=405 total). Thus, if only one (distal) channel were used for unit identification, 34% (137 out of 405) of the motor units would not be identified in these 3-minute records.

3.3 Tracking performance

The tracking accuracy was computed by comparing the classification of global class templates by the algorithm over time against the results of an expert in EMG recognition (Fig. 8, Table 2). For all 132 global motor unit classes, accuracy averaged 89.9% (SE: 0.98%, range 50.0% to 100.0%). More than 80% of all global motor unit classes had tracking accuracies (n=105) that were better than 85%. Reductions in tracking accuracy largely reflected false positives. Fig. 9 shows templates and overlaid potentials from two global motor unit classes with different accuracies. The global motor unit class in A had a classification accuracy of 95.0% and a tracking accuracy of 98.2%. In B, the corresponding data were 80.0% and 90.6%.

Fig. 8
Cumulative distribution of tracking performances. Tracking accuracy was sorted for all global motor unit classes for all recordings then plotted. The same order was then applied when plotting the false positive misclassified data.
Fig. 9
Global motor unit class tracking and clustering. Global motor unit class templates (upper) and overlaid potentials (lower) from the distal and proximal channels for two units at the start of the recording, at hour 14, and the last hour of the recording ...
Table 2
Tracking accuracy statistics (n=132 global motor unit classes)

Since there could be in excess of 100 motor units active during any 24-hour EMG recording, neither the classification algorithm, nor an EMG classification expert could track all of the units over time. The total number of motor units that were not tracked was estimated by determining the average number of templates found per class in the 3-minute records. Each global motor unit class could have at most 1 template in each 3-minute record, for a maximum of 480 templates in a 24-hour recording. The average number of missed global motor unit classes was computed by dividing the total number of unclassified templates in a recording by the average number of templates per global motor unit class. For example, for subject 1 there were 6372 global motor unit class templates for all 22 of the global motor unit classes that were tracked, making an average of 290 global motor unit class templates per motor unit class (Table 3). The average global motor unit classes missed is simply computed by dividing the 626 unclassified templates by the 290 average templates per class to yield approximately 2 motor unit classes that were missed. Thus, the estimated percentage of global motor unit classes that were tracked for the subject 1 recording was 91% (22 tracked out of an estimated 24 total global motor unit classes).

Table 3
Tracking ability: Estimates of the percentage of the total number of motor units tracked

The classification algorithm was able to track 132 motor units out of the total estimate of 164 active motor units in 7 recordings (76%). The tracking ability ranged from a minimum of 55% to a maximum of 92% (standard error = 2%). On a per recording basis, the classification algorithm was able to track the activity of an average of 19 motor units ranging from a minimum of 7 motor units to a maximum of 48 motor units (standard error = 2%).

3.4 Execution time

The automated EMG classification program required just over two weeks on average to process a 24-hour recording, but manual classification times were much longer, including the stages that must be manually accomplished by a human operator) on a 2.8 GHz Intel processor based machine (Table 4).

Table 4
Classification program approximate execution times

The global template formation and manual tracking verification steps are the same in the case of automatic and manual classification because these steps are required in manual classification as well. A classification expert required a complete 40 hour week (5 days) to manually classify 12 minutes of EMG with assistance from software. If extrapolated to the entire 24-hour recording, the time required skyrockets to about 2 and a half years.

4.0 Discussion

A multi-step classification algorithm was devised to automatically classify single motor unit activity in 24-hour surface EMG recordings from paralyzed thenar muscles. The algorithm was able to classify potentials as belonging to global motor unit classes and track these classes over 24 hours as accurately as a person. The major differences were autonomy and time. The algorithm was able to classify entire 24 hour records with minimal input from a user and was much faster than a person.

4.1 Classification Performance

There is no feasible way to ascertain the entire performance of the algorithm over a full 24-hour EMG recording. This is due to the inherent need for an algorithm of this type in the first place. It would take a human operator more than 2.5 years to manually classify an entire 24-hour EMG recording, more than 157 times that needed by the algorithm to accomplish the same feat. This time estimate is conservative due to the fact that humans can make errors and make inconsistent decisions across time.

The classification algorithm can only be expected to perform as well as an expert EMG classifier. Since two EMG classification experts agreed with each other 95% of the time, on average, and the classification algorithm performed similarly on those same 7 classes (92% accuracy), these performance statistics provide a reasonable estimate of the ability of the algorithm to correctly classify single motor unit potentials over 24 hours. Detailed performance comparisons for all of the 132 units classified over 24 hours show that, without any manual editing, the algorithm correctly identified an average of 90.3% of the potentials classified by a person. The algorithm results were validated further by a two sources approach. When surface EMG plus intramuscular signals were both processed by the classification algorithm separately, the classification agreement was 93.9% for two motor units. The algorithm accuracy for these same units was 96.8% compared to an expert classifier. Together, all of these verifications strongly suggest that the algorithm performance was robust because there is a low probability that motor unit potentials could be misclassified so consistently by two people, two independent sources of EMG, as well as the algorithm and a person.

4.2 Sources of error in classification

There were several error sources affecting automatic motor unit classification and tracking. A low sampling rate caused some classification inaccuracies. For example, global motor unit templates that had w-shaped waveforms had poorer accuracies in both classification and tracking (Fig. 10). About half (n=13) the global motor unit classes with accuracies less than 85% had w-shaped templates in one channel. W-shapes are problematic because the relative positions of the peaks in potentials could not be consistently recorded by the data logger when sampling at 3.2 KHz. Alignment of these potentials was problematic because either of the negative points on the “w” could be the most negative point in the potential at different times, even within the same 3-minute record of EMG. Recording with a higher sampling rate would help ensure more consistent recording of the peak positions of the potentials.

Fig. 10
Hard to track motor unit class. Three representative template potentials from one global class recorded at times that differ by more than a few hours. The “W” shapes in these potential representations make this motor unit class hard to ...

Superpositions decreased the classification accuracy of global motor unit classes an average of 5.1% (Fig. 7), performance decrements that were influenced to the level of muscle activity. Accuracy for the subject 2 recording accuracy only decreased an average of 1.1% due to superpositions because the least number of global classes (n=7) were tracked in this recording and only 26,180 potentials were classified. Classification accuracy was decreased most by superpositions in the subject 3 recording (11.2%). More than twice as many units (n=16) were tracked and this record had more than four times as many classified potentials (138,585).

The peel-off superposition resolution method employed here (Stashuk, 2001) is most successful with partial superpositions rather than fully destructive superposition where the sum out of phase waveforms completely cancel each other. In order for this situation to be problematic for the dual channel peel off method employed here, destructive superposition would have to occur in both channels simultaneously. The probability that fully destructive superpositions occurred in our data is very low because all of the motor unit template shapes observed throughout all the recordings were unique and none were completely out of phase with each other for a given channel. To alleviate destructive superpositions, one could either employ more EMG channels, so that an additional channel is more likely to have a partial superposition than a destructive one, or utilize timing information if the superposition involves a periodic motor unit. The addition of timing information could help identify when a motor unit could fire and locate where potentials could occur in the record. Unfortunately, most motor units observed in this study fired sporadically (119 out of 132, or 90.1%), precluding the use of timing information as a classification feature.

The peak to peak amplitude of the global classes recorded in either the distal or proximal channels was unrelated to classification accuracy (distal and proximal EMG both r2≤ 0.005). Amplitude resolution did impact classification and tracking performances, however. The amplitude of the average global motor unit class potential only occupied about 25% of the full scale range of the data logger although most of the range was filled when muscle spasms were triggered (Fig. 11). Those potentials with recorded peak to peak values of less than 18 µV (0.3% of the logger range) were too small to be accurately classified by an expert. After smoothing and filtering these motor unit potentials were visible, legitimate motor unit potentials, but when recorded with low amplitude resolution, these smaller potentials did not always assume the same waveform shape, making their classification unreliable. If the gain were increased, the EMG would have been recorded with a higher amplitude resolution and smaller motor unit potential waveforms would have been more stable. When more classifiable potentials are detectable, more global motor unit classes can be formed and superpositions can also be added. Thus, overall algorithm performance is, in part, a balance between the amount of motor unit activity and its magnitude.

Fig. 11
Amplitude resolution. Approximately 25 % of the available voltage range of the data logger was used to record the motor unit potentials. The potentials that could be classified had raw voltage peak to peak values of approximately 300 µV (A, B). ...

4.3 Motor unit tracking ability

Motor unit tracking and detection were enhanced by using two channels of EMG. The recording configuration used in this study positioned electrodes longitudinally along the muscle fiber, and transversely, since the electrodes were wide enough to overlap many fibers simultaneously. The presence of two channels enabled 405 motor units to be identified over 32 3-minute records. Only 268 motor units would have been identified if only the distal channel was used, or 237 motor units would have been found if only the proximal channel were used. On average, adding a second channel boosted the identifiable motor units by 51% (405 versus 268 units).

Even though the tracking ability could only be estimated, the percentage of motor units that could be identified by classification software was comparable to that found by Farina et al. (2008). These authors recorded surface EMG from the abductor digiti minimi muscle using different sized square grids of circular electrodes (diameter 1 mm spaced at 5 mm distances) in bipolar and laplacian configurations. More motor units were detected when more channels, and therefore more signal features, were used. The estimated tracking ability in our study (76%) is on the order of that obtained using a 5 × 5 grid of electrodes in a bipolar or laplacian configuration (Farina et al., 2008). However, the tracking ability in our study is even more notable when considering that the motor units were followed by an algorithm for 24 hours, rather than a few minutes. In future work, additional EMG channels could be incorporated into the classification algorithm and this may boost both classification and tracking accuracy.

The ability of the algorithm to track different motor units over time was also influenced by the number of active motor units. When more units were active, fewer motor units could be tracked due to the presence of superpositions. If a given motor unit is more active, more superpositions will result as well. For example, the algorithm was able to track an estimated 55% (n=48) of the global motor unit classes in the subject 7 recording. The large number of unclassified templates (4384) was responsible for the poor tracking ability. These unclassified templates may have been complex superpositions that the classification program could not resolve.

Tracking errors may have resulted from template changes over time. These changes may have resulted from alterations in the conductive properties of the electrodes or electrode movement. Electrodes were taped onto the skin to ensure that they did not move. The maximal M-waves that were recorded before and after the 24 hour recording differed by an average of 7% for the distal EMG and by 2% in the proximal channel. These data illustrate that electrodes do not necessarily behave like matched pairs over time and that the signals recorded by both channels can change by varying amounts over time. Even so, only a 7 % change in potential amplitude can result from electrode issues. Changes in global template amplitude that exceed this level must arise from other sources.

4.4 Multiple-pass versus sequential method

Motor unit potentials were classified over several passes in this study, instead of being classified sequentially as they were segmented as in most other EMG decomposition packages (McGill et al., 2005; Zennaro et al., 2003). Our multiple-pass procedure for classification provided a broad view of the data both for decision making and for robustly tracking potentials over 24 hours. One problem with sequential algorithms is that motor unit potentials that fire sporadically at first and then regularly later on in the recording will have some potentials missed if the earlier portions of the recording are not re-analyzed. The classification method presented here analyzes the entire recording after all motor units template shapes are found to make sure that no motor unit firings are missed. The average motor unit tracked by the classification algorithm across all seven recordings started firing 2 hours and 18 minutes from the start of the experiment. On average, 3.8% (6,145 out of 163,169 per recording) of all motor unit potentials that the classification algorithm was able to classify would be missed had the classification method been sequential and not re-analyzed.

Our non-sequential classification technique also provides for more robust tracking of motor unit potentials over 24 hours. EMG-LODEC (Zennaro et al., 2003) accomplishes motor unit tracking using a weighted averaging technique with a forgetting factor to update template waveforms. This technique gives more weight to recently occurring templates. In contrast, our classification algorithm tracks motor unit potentials over 24 hours by uniting the motor units that have been individually classified within a series of small records. In effect, the weighted averaging method was replaced with an additional classification stage – the classification of motor unit potential templates found in each of the smaller records over the entire recording. Our tracking method is more robust because when tracking is done over the entire recording a posteriori, more information is available – individual templates and member potentials – to make more informed classifications for tracking. With a weighted average, the current template waveform is updated from previous templates, with each successive template being computed using more information that the previous one.

4.5 Execution Time

The automated EMG classification program required just over two weeks on average to process a 24-hour recording on a 2.8 GHz Intel processor based machine. These execution times seem lengthy, but they were yielded by average computers. The execution times could be reduced to a day or less if more powerful computers were used.

Execution times could also be reduced by reducing the sampling rate, thus reducing the data to process, but the tracking ability of the algorithm would be negatively impacted. The algorithm could only identify about 58% (25 out of 43) of the motor units in the data sampled at half the rate. Smaller potentials (peak to peak amplitudes < 50 µV) were not stable enough to be clustered and identified at the 1600 Hz sampling rate, leading to further bias against motor units with smaller potentials, either because they are weak motor units or deep below the skin. Halving the sampling rate did not alter the classification accuracies, however.

The superposition resolution procedures for the classification algorithm presented here also limited computational complexity and therefore reduced execution time. Small reductions in the execution times of the superposition resolution procedure can drastically reduce the overall execution time for a 24-hour recording because, on average, 844,000 candidate motor unit potentials had to be inspected per recording. Instead of trying all possible combinations, only the highest correlating templates in each channel at each peel off stage were selected. In this way, superimposed potentials were resolved more optimally and at an increased speed because the peel off was done on each channel simultaneously.

Furthermore, time was spared by not using firing rate information to classify motor unit potentials. Holobar et al. (2009) utilized convolution kernel compensation, a non-template based method of classification, which did not require motor units to fire periodically either. However, this method did involve the complexity of computing the inverse of a correlation matrix, limiting the execution speed of the algorithm and the amount of data that can be processed at once.

Finally, while algorithm execution times can be reduced, the time required for an expert to manually classify EMG remains the same. It was estimated that a classification expert would require about 2 and a half years to complete the same analysis. This extraordinary amount of time not only makes the manual classification of an entire 24-hour recording virtually impossible, but also demonstrates the importance of developing this automated classification software.

4.6 Applications

The classification algorithm described here enables a detailed analysis of long-term single motor unit activity, including the computation of motor unit firing patterns. At the gross level (Fig. 12), the analysis yields the times the tracked motor units were active. This can be expanded to investigate activity differences between awake and asleep time, to investigate the amounts of neuromuscular activity required to sustain muscle and motor unit properties (Kernell et al., 1987; Kim et al., 2007) or to evaluate whether medications dampen activity. Muscle spasms are a quality of life issue for spinal cord injured individuals, for example by disrupting sleep. Many SCI individuals take anti-spasm medication such as baclofen to dampen the responsible involuntary muscle activity (Taricco et al., 2006). These drugs may suppress the activity of some motor units but not others. For example, Subject 2 administered a baclofen/clonazepam dosage at 3 different times during the 24 hour recording (Fig. 12). The approximate times of the drug administration are shown as vertical lines along with the dose number (as in number of administrations). Using the classification algorithm, we can perhaps answer why medication is more effective in some individuals than others since only some motor units may be affected by anti-spasm medications.

Fig. 12
Motor unit activity over 24 hours. Global motor unit class template shapes are shown on the left side and the motor unit activity is displayed on the right for global motor unit classes 1, 4, and 6 for the subject 2 recording. For a motor unit class to ...

At a finer level, the firing rate histograms produced by the algorithm can be used to identify possible sources of spontaneous motor unit activity. There were two main types of firing behavior during the 24-hour recordings, sustained and sporadic (Fig. 13). The motor units that fired tonically often had average firing rates near 4–6 Hz, similar to the firing rates observed in other people (Zijdewind and Thomas, 2001; Zijdewind et al., 2004) and may reflect the activation of persistent inward currents in the motoneuron (Heckman et al., 2005). In contrast, the occasional firing of a motor unit at irregular rates, often below 1Hz, may reflect activation of the motoneuron by synaptic noise (Matthews, 1996), and represent the fasciculation potentials often seen in motoneuron disorders like amyotrophic lateral sclerosis (Kleine et al., 2008). The algorithm described here can classify these kinds of long duration motor unit recordings and help further our understanding of their biological importance.

Fig. 13
Motor unit firing rate histograms. The y-axis presents the motor unit potentials that fired for a given global class as a percentage. Global class 21 shows sporadic firing behavior because most of the potentials are firing below 1 Hz. Global class 5 exhibits ...

Acknowledgments

The authors would like to thank Drs. Dejan Tepavac and Chuck Cone. The completion of this project would not have been possible without the guidance of Dejan Tepavac and the oscilloscope probe that was graciously provided by Chuck Cone, Chair of the Electrical Engineering Department of Embry Riddle Aeronautical University, Prescott, Arizona. This research was funded by USPHS grant NS-30226, State of Florida Appropriation #677, and The Miami Project to Cure Paralysis.

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References

  • Chauvet E, Fokapu O, Hogrel JY, Gamet D, Duchene J. Automatic identification of motor unit action potential trains from electromyographic signals using fuzzy techniques. Med Biol Eng Comput. 2003;41(6):646–653. [PubMed]
  • De Luca CJ, Adam A, Wotiz R, Gilmore LD, Nawab SH. Decomposition of surface EMG signals. J Neurophysiol. 2006;96(3):1646–1657. [PubMed]
  • Drost G, Kleine BU, Stegeman DF, Van Engelen BG, Zwarts MJ. Fasciculation potentials in high-density surface EMG. J Clin Neurophysiol. 2007;24(2):301–307. [PubMed]
  • Edin BB, Backstrom A, Backstrom LO. Single unit retrieval in microneuorgraphy: A microprocessor-based device controlled by an operator. J Neurosci Methods. 1988;24(2):137–144. [PubMed]
  • Farina D, Francesco N, Gazzoni M, Enoka R. Detecting the unique representation of motor unit action potentials in the surface electromyogram. J Neurophysiol. 2008;100(3) 1223-123. [PubMed]
  • Florestal JR, Mathieu PA, Malanda A. Automated decomposition of intramuscular electromyographic signals. IEEE Trans Biomed Eng. 2006;53(5):832–839. [PubMed]
  • Florestal JR, Mathieu PA, McGill KC. Automatic decomposition of multichannel intramuscular EMG signals. J Electromyogr Kinesiol. 2009;19(1):1–9. [PubMed]
  • Heckman CJ, Gorassini M, Bennett DJ. Persistent inward currents in motoneuron dendrites: Implications for motor output. Muscle Nerve. 2005;31(2):135–156. [PubMed]
  • Holobar A, Farina D, Gazzoni M, Merletti R, Zazula D. Estimating motor unit discharge patterns from high-density surface electromyogram. Clin Neurophysiol. 2009;120(2):551–562. [PubMed]
  • Kawamura J, Ise M, Tagami M. The clinical features of spasms in patients with a cervical cord injury. Paraplegia. 1989;27(3):222–226. [PubMed]
  • Kernell D, Eerbeek O, Verhey BA, Donselaar Y. Effects of physiological amounts of high-and low-rate chronic stimulation on fast-twitch muscle of the cat hindlimb. I. Speed- and force-related properties. Journal of Neurophysiol. 1987;58(3):598–613. [PubMed]
  • Kim SJ, Roy RR, Zhong H, Suzuki H, Ambartsumyan L, Haddad F, Balwin KM, Edgerton VR. Electromechanical stimulation ameliorates inactivity-induced adaptations in the medial gastrocnemius of adult rats. J Appl Physiol. 2007;103(1):195–205. [PubMed]
  • Kleine BU, Stegeman DF, Schelhaas HJ, Zwarts MJ. Firing pattern of fasciculations in ALS: Evidence for axonal and neural origin. Neurology. 2008;70(5):353–359. [PubMed]
  • Kleine BU, Van Dijk JP, Lapatki BG, Zwarts MJ, Stegeman DF. Using two-dimensional spatial information in decomposition of surface EMG signals. J Electromyogr Kinesiol. 2007;17(5):535–548. [PubMed]
  • Little JW, Micklesen P, Umlauf RC, Britell C. Lower extremity manifestations of spasticity in chronic spinal cord injury. Am J Phys Med Rehabil. 1989;68(1):32–36. [PubMed]
  • Mambrito B, De Luca CJ. A technique for the detection, decomposition, and analysis of the EMG signal. Electroencephalogr Clin Neurohyisol. 1984;58(2):175–188. [PubMed]
  • Matthews PBC. Relationship of firing intervals of human motor units to the trajectory of post-spike after-hyperpolarization and synaptic noise. J Physiol. 1996;492(2):597–628. [PubMed]
  • McGill K, Lateva ZC, Marateb HR. EMGLAB: An interactive EMG decomposition program. J Neurosci Methods. 2005;149(2):121–133. [PubMed]
  • Stashuk DW. EMG signal decomposition: How can it be accomplished and used? J Electromyogr and Kinesiol. 2001;vol. 11(3) 151-1732001. [PubMed]
  • Stein RB, Brucker BS, Ayyar DR. Motor units in incomplete spinal cord injury: Electrical activity, contractile properties, and the effects of biofeedback. J Neurol Neurosurg Psychiatry. 1990;53(10):880–885. [PMC free article] [PubMed]
  • Taricco M, Pagliacci MC, Telabo E, Adone R. Pharmacological interventions for spasticity following spinal cord injury: Results of a cochrane systematic review. Eura Medicophys. 2006;42(1):5–15. [PubMed]
  • Thomas CK. Contractile properties of human thenar muscles Paralyzed by Spinal Cord Injury. Muscle Nerve. 1997;20(7):788–799. [PubMed]
  • Thomas CK, Ross BH. Distinct patterns of motor unit behavior during muscle spasms in spinal cord injured subjects. J Neurophysiol. 1997;77(5):2847–2850. [PubMed]
  • Westling G, Johansson RS, Thomas CK, Bigland-Ritchie B. Measurement of contractile and electrical properties of single human thenar Motor units in Response to intraneural motor-axon stimulation. J Neurophysiol. 1990;64(4):1331–1338. [PubMed]
  • Yang JF, Stein RB. Methods for estimating the number of motor units in human muscles. Ann Neurol. 1990;28(4):487–495. [PubMed]
  • Zennaro D, Wellig P, Koch VM, Moschytz GS, Laubli T. A software package for the decomposition of long term multichannel EMG signals using wavelet coefficients. IEEE Trans Biomed Eng. 2003;50(1):58–69. [PubMed]
  • Zijdewind I, Thomas CK. Spontaneous motor unit behavior in human thenar muscles after spinal cord injury. Muscle Nerve. 2001;24(7):952–962. [PubMed]
  • Zijdewind I, Peterson L, Thomas CK. Firing patterns of spontaneously active motor units in spinal cord injured subjects. presented at the Active Dendrites in Motor Neurons Conference; Jun. 2004; Boulder, Colorado. [PubMed]
  • Zouridakis G, Tam DC. Identification of reliable spike templates in multi-unit extracellular recordings using fuzzy clustering. Comput Methods Programs Biomed. 2000;61(2):91–98. [PubMed]