Reliability of spike triggered averaging of the ... - Wiley Online Library

RELIABILITY OF SPIKE TRIGGERED AVERAGING OF THE SURFACE ELECTROMYOGRAM FOR MOTOR UNIT ACTION POTENTIAL ESTIMATION XIAOGANG HU, PhD,1 WILLIAM Z. RYMER, MD, PhD,1,2,3,4 and NINA L. SURESH, PhD1 1

Sensory Motor Performance Program, Rehabilitation Institute of Chicago, 345 E Superior Street, Room 1378, Chicago, Illinois, 60611, USA 2 Department of Physical Medicine and Rehabilitation, Feinberg School of Medicine, Northwestern University, Chicago, Illinois, USA 3 Department of Physiology, Northwestern University, Chicago, Illinois, USA 4 Department of Biomedical Engineering, Northwestern University, Chicago, Illinois, USA Accepted 13 February 2013 ABSTRACT: Introduction: The reliability of estimated motor unit parameters using spike triggered averaging (STA) of the surface electromyogram (sEMG) has not been tested thoroughly. We investigated factors that may induce amplitude bias in estimated motor unit action potentials (MUAPs) and shape variations. Methods: An sEMG record was simulated. MUAPs were then estimated from the STA of the simulated EMG. Results: Variations in MUAP duration led to under-estimation of real MUAP amplitude, while synchronized firing led to over-estimation of amplitude. Spurious firing resulted in over-estimation of the amplitude of small motor units but under-estimation of the amplitude of large ones. Variability in amplitude and high firing rates had minimal influence on amplitude estimation. High firing rates and variation in MUAP duration led to large variations in MUAP shape. Estimation errors also correlated with shape variations. Conclusions: Recommendations to enhance the accuracy of the STA estimates have been proposed. Muscle Nerve 48: 557–570, 2013

Spike triggered averaging (STA) is a system-identification technique that has been used extensively to understand neuromuscular properties by focusing typically on motor unit electrical or contractile features. The STA procedure aligns segments of a recorded electrical or mechanical signal relative to time-locked triggering events and can attenuate non–time-locked activities from the complex signal and extract information related uniquely to the triggering events. The method is usually used to align discrete neural events, such as a neuronal discharge, with a continuous signal such as the electromyogram or a muscle force record. STA of electromyographic (EMG) activity has been used to identify functional connectivity between cortical neurons with motoneurons innervating muscle fibers.1–3 Motor unit (MU) properties such as the number of MUs,4,5 MU conduction velocity,6,7 and MU size distribution8,9 have also Abbreviations: CV, coefficient of variation; EMG, electromyogram; ISI, inter-spike interval; MU, motor unit; MUAPs, motor unit action potentials; P–P, peak–peak; sEMG, surface EMG; STA, spike triggered averaging Key words: action potential amplitude variation; action potential classification; duration variation; firing rate; Spike triggered average; synchronized firing Correspondence to: X. Hu; e-mail: [email protected] This study was supported by NIH R24 HD50821-07 awarded to W.Z.R. C 2013 Wiley Periodicals, Inc. V

Published online 20 February 2013 in Wiley (wileyonlinelibrary.com). DOI 10.1002/mus.23819

Reliability Assessment of STA of sEMG

Online

Library

been estimated using STA of surface recorded or intramuscular EMG. In clinical neurophysiology, MU properties derived from the STA of EMG have been applied to assess neuronal diseases such as amyotrophic lateral sclerosis (ALS) and spinal muscular atrophy.10,11 STA can also be used to examine neuronal degenerative changes (e.g., changes in action potential amplitude and morphological properties) due to aging4,9,12 or accompanying a cerebral stroke.13–16 Although stroke is a central nervous system disorder, there are known changes in motoneuron and muscle fiber morphology and function that affect the features of the motor unit action potential and its firing patterns. Despite the widespread use of STA, this technique has not been tested thoroughly to determine the influence of various concurrent factors on the reliability of STA estimation, including amplitude bias in the estimated motor unit action potential (MUAP). The objective of this study was to examine the quantitative characteristics of estimated MUAP shape derived from STA in relation to the known MUAP and the sEMG, and to provide guidelines regarding the derivation of a secure STA estimate from such sEMG signals. Motor unit parameters estimated by STA may be biased by various factors. During sEMG recordings, the action potentials composing the EMG interference pattern may change, which may bias the averaged action potential derived from STA. Superposition of action potentials from different MUs can also distort the STA-derived MUAP estimate. Decomposition errors in MU firing event detection may bias the action potential estimation. In particular, spurious firings that belong to other MUs with different action potential shapes may distort the STA estimates. It has been shown that spurious triggering events can also bias motor unit number estimation.17 In light of these uncertainties, we examined several factors that can influence the reliability of STA derived MUAP. Five variables were examined: amplitude variations within a MUAP train, varying duration of a MUAP train, action potential superposition due to high firing rates, synchronized MUSCLE & NERVE

October 2013

557

firing effects, and spurious event classification during firing event discrimination. To quantify the influence of these factors, sEMG was simulated from experimentally derived action potential shapes and firing properties. Overall, we identified possible sources of systematic errors in STA-estimated MUAPs, and our observations may help provide a means for assessing the reliability of the STA estimates. We recommend that certain procedures be followed in signal acquisition and processing to minimize the MUAP estimation bias during the application of STA techniques. METHODS Derivation of MU Discharges and Action Potentials.

To derive physiologically realistic properties of a typical sEMG record, the MUAP and respective discharge properties were estimated from raw sEMG data collected from a healthy subject. The subject performed isometric index finger abduction at 30% of maximum force. sEMG of the first dorsal interosseous muscle was recorded using a surface sensor array (Delsys Inc.). Single MU firing events and MUAP templates were then extracted from a high-yield decomposition system (version 1.0.0.28). Detailed information for the decomposition algorithms has been described previously.18 Representative MU templates and firing events were derived from the output of the decomposition algorithm. The signal-derived MUAP and firing events from the decomposition algorithm of the 13 MUs are shown in Figure 1A and 1B. Reconstruction of EMG. The firing time statistics (i.e., mean firing rate and coefficient of variation [CV] of firing rate) and MUAPs of the 13 units were used to construct a simulated sEMG signal. For the reconstructed EMG, the firings and MUAP templates used here represent the “true” properties for each selected MU. In the subsequent sections, the MUAPs used to derive the simulated sEMG are termed true MUAPs. The EMG simulation diagram is shown in Figure 2A. First, the MUAP train of each MU was constructed from the manipulated MUAP template and firing events (The detailed manipulation at each component is described in the next section). The MUAP trains of the 13 MUs were then summed linearly, and baseline noise was added to form the interference EMG signal. Baseline noise was simulated from a Gaussian distribution with a zero mean and a standard deviation (do), being the same as the original EMG baseline noise (calculated from the first 5 s quiescent period where no MU activity was detected). Spike Triggered Averaging. Spike triggered averaging (STA) was performed on the simulated EMG signal. The different variations in the EMG 558


FIGURE 1. Derivation of motor unit action potential (MUAP) templates and firing properties from real surface electromyogram (sEMG). A: 13 representative MUAP templates. B: Firing properties of the 13 MUs.

simulation are described in the next section. The firing times of individual MUs were used as triggers for the STA calculation. The time interval used to derive the action potential estimate was set as beginning 10 ms before and ending 10 ms after the firing time. The firing time was then placed at the algebraic center of the time window. Two separate STA features were calculated. These features were designed to assess the variability of the waveform shape and the relation with the amplitude bias in STA estimates. 1. The STA estimated MUAP was calculated based on a window (EMG segment) length of 4 s, and the MUAP variation across different segments was then estimated. This 4-s epoch yielded approximately 50–100 firing events in each window, and the window was then shifted over the sEMG signal using a step size of 0.5 s. The coefficient of variation (CV: standard deviation normalized by the mean) of the peak–peak (P–P) amplitude of the MUAP, for each step change, was calculated as a measure of the variability of the waveform shape. One exemplar MUAP derived from STA of recorded sEMG over different windows is shown in Figure 2B, and the different sliding window is denoted by window index. The P–P amplitude was the difference between the maximal positive and negative peaks of the MUAP. The P–P amplitude and duration over different windows is shown in Figure 2D and 2E. For this particular EMG MUSCLE & NERVE

October 2013

FIGURE 2. Factors influencing spike triggered averaging (STA) estimates and STA features. A: Reconstruction of electromyographic (EMG) and simulation of the 5 factors influencing the STA estimates. Gaussian noise was added to the motor unit action potential (MUAP) amplitude or duration to simulate shape variations. ISIs were scaled to simulate different firing rates, and temporal synchrony was imposed on the firing events. Baseline Gaussian noise was also added to the constructed EMG signals. STA of constructed EMG was performed, and spurious firing events were added to simulate firing detection errors. B: Exemplar MUAP estimates over the sliding window. The coefficient of variation (CV) of P–P amplitude was calculated over different window segments. Window index denotes different window segments over time. C: Correlation estimate between STA-derived MUAP and “true” MUAP. D,E: P–P amplitude and duration of the 13 MUs estimated over different window segments.

signal, the CV of amplitude ranged from 0.052 to 0.200, and the CV of duration ranged from 0.043 to 0.291. 2. The MUAP estimations derived using the STA method were compared with the “true” MUAPs. The maximum linear correlation coefficient between the STA estimate (calculated over the entire trial duration) and the “true” MUAPs was computed as a second measure of the reliability of the waveform average. The STA-derived template was shifted 10 ms backward and forward relative to the “true” MUAPs to find the maximum correlation coefficient. This correlation quantifies the shape similarity between STA and “true” MUAPs and signifies the reliability of the STA waveform average. One exemplar MUAP shape correlation is shown in Figure 2C. Possible

Sources

of

Bias

and

Variation

in

STA

The MUAP shape variation can arise from both amplitude and duration changes during the recordings. Superposition of action potentials between units can also come from high firing rates and from synchronized firing between units. In addition, human or machine errors in firing event detection can also influence the reliability of STA estimates. These parameters were simulated 1 at a time during EMG reconstruction.

Estimate.

Amplitude Noise. In the simulation, the amplitude variation across spikes for a particular MU was modeled as a Gaussian distribution with a zero mean and a standard deviation (dAmp) scaled with the P–P Reliability Assessment of STA of sEMG

amplitude of the MUAP. At each level of amplitude noise, a constant coefficient of variation (CV) of P–P amplitude was specified across all obtained units. Different noise levels (calculated based on the CV of P– P amplitude) ranging from 0 to 0.4, with an increment of 0.05 were simulated. We have simulated CV of amplitude up to 0.4, which is expected to be much higher than that in real experimental settings, because a MUAP with 40% amplitude change will likely to be classified as a different MU. This also holds for other simulated factors. The original firing events from the decomposition system routinely display a certain level of synchrony between units (see the Results section) as reflected in the realistic firing properties of MUs.19,20 This synchrony may also induce bias and variations in the STA estimate. To avoid the confounding influences of amplitude variability and synchronized firing on STA estimate errors, simulated firing events were obtained to form action potential trains. Specifically, the inter-spike intervals (ISI) of each unit were modeled to follow a Gaussian distribution with the mean and CV calculated from the original firing times derived from the steady state firing stage. Simulated sEMG was then obtained using the MUAPs with varying amplitude and nonsynchronized spike trains. An STA analysis of the simulated EMG was then performed using the simulated MU d firing times. The STA estimated MUAP (MUSP) from the reconstructed EMG was compared with the true MUAP used to simulate the EMG. The nord Amp 2MUAPAmp Þ= malized amplitude error [ðMUSP MUSCLE & NERVE

October 2013

559

MUAP Amp 3100] of each MU was calculated to quantify the STA estimation bias. The 2 STA features (CV of P–P amplitude across segments of sEMG signal and correlation d and –UAP) were calcucoefficient between MUSP lated. The changes in these 2 parameters as a function of added amplitude noise were then computed to quantify the estimation variability. The correlation of the 2 features with the normalized amplitude error (calculated in absolute term without considering over- or under-estimation error) was also calculated to examine whether the STA features could predict the amplitude error, because the amplitude error from the STA estimate is typically unknown in experimental EMG data. Duration Noise. Similar to the amplitude noise simulation, the duration variation across spikes for a particular MU was modeled as a Gaussian distribution with a zero mean and a standard deviation (dDuration) that was scaled with the duration of the MUAP. At each level of duration noise, a constant CV of P–P duration was specified across all MUs. Different noise levels (based on the CV of P–P duration) ranging from 0 to 0.4 with an increment of 0.05 were simulated. The same simulated MU firing events used in the amplitude variation were used for the duration variability simulation. Simulated EMG was obtained from these firing times and action potentials with varying duration. STA analysis of the simulated EMG was then performed; amplitude estimation error and STA features (CV of P–P amplitude and correlation coefficient) were calculated to evaluate the influence of duration noise on the STA estimate. Adjustment of Firing Rate. High firing rates increase the degree of action potential superposition between MUs, which can potentially compromise the STA estimation. To test the influence of firing rate on the STA estimates, different interspike interval (ISI) ranges were simulated. Specifically, the ISIs of each MU followed a Gaussian distribution with a mean that was scaled based on the original firing events. The scaling of ISI ranged from 0.5 to 4 with an increment of 0.5. A scaling factor of 0.5 means that the mean ISI was half of the original ISIs, and a scaling factor of 4 means the mean ISI was 4 times the original ISIs. The CV of ISIs in all simulated conditions was maintained at the same level as the original ISIs. The sEMG was constructed using the simulated firing times and MUAP templates. The STA analysis of the simulated EMG was performed, and amplitude error and MUAP features (CV of amplitude and correlation coefficient) were calculated to evaluate the influence of firing rate on the STA 560


estimate. With the changes of ISI, the number of firing events in each 4-s window also changes. To make the reliability of STA estimation comparable between different scaling factors, the moving window used to calculate the CV of P–P amplitude was also adjusted in the same scale as the ISI; namely, the window length was 2 s (4 s 3 0.5) in the scaling factor of 0.5 condition, and the window length was set at 16 s (4 s 3 4) when the scaling factor was 4. Accordingly, the duration of the simulated EMG signal was also adjusted such that the duration was 4 times longer than the original when the scaling factor was 4. Synchronization. MUAP superposition can also arise from synchronized firing across different MU trains. To determine the influence of such synchronization on the STA estimate, the timing of independent firing events was adjusted to impose temporal correlation between some of the action potentials of different MUs. The synchronization was simulated using a previously described method.21,22 Different levels of synchronization were simulated. The percentage of firing events selected as reference firing and the percentage of MUs selected as synchrony MUs varied from 0% to 40%. The actual level of synchronization from the decomposition-derived MUAP trains was estimated using the synchronization index.23 First, the crosscorrelation histogram24 between MU pairs was calculated. Then the peak around time 0 was estimated using the synchronization index [(peak area – mean area) / (total area/2)] as a measure of synchrony between MU pairs. Detection Errors in Firing Events. During identification of spike times, there can be missed spikes and/or added spurious firing errors. A firing omission does not influence the STA estimate, because that action potential is missed and will not be taken into account during the averaging calculation. In contrast, a false positive error (i.e., a spurious firing) may influence the reliability of the STA estimates, because an action potential of 1 unit is misclassified into another unit. Therefore, only spurious firings were injected during STA calculation. Specifically, firing times were simulated from a Gaussian distribution with the mean and standard deviation derived from the original firing events. No synchronization was imposed between units. sEMG was obtained using simulated firing times and invariant MUAP templates across firing events. Spurious firings were then injected into the firing times. The injected firing error satisfied the following constraints: first, the minimum ISI after error injection is not smaller than the original minimum ISI of the unit, and second, the injected firing times come from existing firings of other units. MUSCLE & NERVE

October 2013

A proportion of firing intervals of 1 unit was selected randomly, and existing firings from other MUs were also selected randomly. The amount of injected spurious firing was calculated as the error injection rate (the number of injected firing events divided by the total number of firing events). STA analysis of sEMG was then performed based on the firing events with spurious firings added. The influence of spurious firing detection on STA estimates was evaluated through the amplitude error and waveform shape variations. RESULTS

Our general results show that variations in action potential duration led to under-estimation of the action potential amplitude and that synchronized firing between units led to over-estimation of the amplitude. Spurious firing errors resulted in over-estimated amplitudes of small motor units but under-estimated amplitudes of large motor units. Variations in action potential amplitude and high firing rates had little influence on the action potential amplitude estimates. Variation in action potential duration and high firing rate led to large variations (increased CV of amplitude and reduced correlation) in action potential shape. The STA-derived action potential features (CV of amplitude and correlation) could also predict the estimation bias. The results of the 5 simulated variables are now described in more detail. Amplitude Noise. The STA estimated P–P amplitude depicted as a function of the true amplitude at different amplitude noise levels is shown in Figure 3A. Overall, the STA derived amplitude correlated with the true amplitude very well. The added amplitude noise had a minimal influence on the amplitude estimate, and there was no consistent shift of amplitude estimates with added noise, regardless of the amplitude of the MUAP. Figure 3B illustrates the normalized amplitude error, which was approximately 10% with no noise added and increased slightly to approximately 16% with increasing noise level. The change of CV of amplitude (estimated across a moving window of 4 s of the EMG segment) due to added noise for each unit is shown in Figure 3C. There was no consistent change of CV of amplitude with increasing noise level as shown for individual MUs (thin lines). When the change of CV was averaged across units (thick line), the maximum increase of CV was only 0.011. The change of correlation is shown in Figure 3D. The reduction of correlation was small for most of the MUs as the noise level increases. The average reduction of correlation was also small (up to 0.07) across units when the highest noise level was added. Reliability Assessment of STA of sEMG

To predict amplitude estimation errors, the STA derived features (including CV of amplitude and correlation) were correlated with the normalized amplitude error. As shown in Figure 3E and 3F, the units with larger amplitude errors tended to have varying action potential shapes (i.e., higher CV of amplitude and lower correlation). The STA estimated P–P amplitude, presented as a function of the true amplitude at different duration noise levels is shown in Figure 4A. Unlike the amplitude noise effect, with increasing duration of added noise, the STA estimated amplitudes tend to decrease consistently, and the magnitude of reduction is stronger (in an absolute sense) in larger MUAPs. The relative amplitude errors (Fig. 4B) of the 13 MUs also reduced with increasing noise level, and the degree of reduction was similar across different units. The change of CV of amplitude for each unit is shown in Figure 4C. Most units showed increased CV of amplitude with increasing noise level. The averaged change of CV across units (thick line) due to added noise increased substantially (up to 0.023). Similarly, a reduction of correlation was evident in all the MUs (thin lines in Fig. 4D), and the average reduction of correlation was up to 0.12. However, there was no obvious association between the normalized amplitude error and the STA features as shown in Figure 4E and 4F.

Duration Noise.

Superposition Due to High Firing Rate. Different mean firing rates were simulated to examine the influence of superposition between MUs on STA estimates. ISIs of all the units were scaled to adjust the firing rates. The STA-derived amplitude as a function of the true amplitude at different ISI scales is illustrated in Figure 5A. To our surprise, there was no systematic bias in the amplitude estimates with increasing firing rates. Similarly, the normalized amplitude error (Fig. 5B) did not show a systematic shift with increasing firing rates (reduction of ISI scales). The change of CV of amplitude due to high firing rate for each unit is shown in Figure 5C. All the MUs showed increased CV of amplitude with increasing firing rates (i.e., reduction of ISI scales). The averaged increase of CV of amplitude (thick line) reached 0.13 at the scale of 0.5 where the firing rates were twice as high as the original firing rates. The scale of ISI showed a minimal influence on the correlation feature both for individual MUs (thin lines in Fig. 5D) and for averaged correlation across units (thick line). As shown in Figure 5E and 5F, the normalized amplitude error correlated well with the STA features. MUSCLE & NERVE

October 2013

561

FIGURE 3. Influence of amplitude variation. A: Amplitude of spike triggered averaging (STA) motor unit action potentials (MUAPs) as a function of true MUAPs. The embedded exemplar action potential (in thick lines) is the true MUAP, and variations of amplitude due to added amplitude noise are shown in thin lines. The color of symbols changes from black to gray as the amplitude noise level increases. B: Normalized amplitude errors of the STA estimates. The size of symbols increases with the size of the MUAP. C: Change of coefficient of variation (CV) of amplitude of individual motor units (MUs) (thin lines) and averaged change of CV of amplitude across MUs (thick line). D: Change of correlation of individual MUs (thin lines) and averaged change of correlation. E: Normalized amplitude error (un-signed) as a function of the CV of amplitude from STA of EMG segments. F: Normalized amplitude error (un-signed) as a function of the correlation. Superposition Due to Synchronization. Different levels of synchrony between MUs were also simulated to examine the influence of superposition on STA estimates. The STA-derived amplitude as a function of the true amplitude is illustrated in Figure 6A. With increasing synchrony, STA tended to over-estimate the amplitude substantially. The relative amplitude error (Fig. 6B) tended to be positive for most MUs, and the amplitude error reached over 100% for some small MUs. The change of CV of amplitude due to synchrony is shown in Figure 6C. There was no consistent increment of CV with the level of synchronization. The averaged increment of CV due to synchronization was also small. The 562


reduction of correlation due to synchrony was also small for most MUs (Fig. 6D), and the averaged correlation reduced by up to 0.07 with increasing synchronization level (thick line). The associations between amplitude errors and STA features are shown in Figure 6E and 6F. MUs with larger amplitude errors tended to have larger CV of amplitudes and lower correlations. The synchronization levels between the original firing spike trains from the decomposition were estimated using the synchronization index.23 The synchronization index matrices of all MU pairs are illustrated in Figure 7A and 7B for 2 separate recorded EMG signals from 2 subjects. The synchronization index with the MU itself was set at MUSCLE & NERVE

October 2013

FIGURE 4. Influence of duration variation. A: Amplitude of spike triggered averaging (STA) motor unit action potentials (MUAPs) as a function of true MUAPs. The embedded exemplar action potential (in thick lines) is the true MUAP, and variations (increase or decrease) of duration due to added duration noise are shown in thin lines. B: Normalized amplitude errors of the STA estimates. C: Change of coefficient of variation (CV) of amplitude of individual MUs (thin lines), and averaged change of CV of amplitude across motor units (MUs) (thick line). D: Change of correlation of individual and averaged across MUs. E: Normalized amplitude error as a function of the CV of amplitude. F: Normalized amplitude error as a function of the correlation.

zero. The cross-correlation histogram used to calculate the synchronization index of 2 exemplar motor unit pairs: MU 7 versus 6 (high synchrony) and MU 7 versus 8 (low synchrony) are displayed in Figure 7C. As shown in Figure 7A, the synchronization indices were above 0.1 for most MU pairs and tended to be higher for larger MU pairs; the maximum synchronization index reached 0.46 for subject 1. The synchronization indices (Fig. 7B) were also above 0.1 for most MU pairs for subject 2 (note that the color coding is different between subjects 1 and 2). Compared with subject 1, the synchronization indices of MU pairs from subject 2 were still nonuniformly distributed but were below 0.18 for most MU pairs. sEMG signals were constructed using the original firing events from these Reliability Assessment of STA of sEMG

2 subjects. STA analysis of the simulated sEMG was performed, and the STA-estimated MUAP amplitudes as a function of the true amplitude for both subjects are shown in Figure 7D. The STA analysis tended to over-estimate the amplitude of MU with higher synchronization indices (e.g., the larger MUs of subject 1). Firing Event Detection Errors. Firing event detection errors (i.e., spurious firing events) were simulated during the STA calculation. One example of spurious firings is shown in Figure 8. The spike trains represent the original firing events used to construct the sEMG, and the crosses represent the randomly injected spurious firing errors. The original firings combined with these spurious firings MUSCLE & NERVE

October 2013

563

FIGURE 5. Influence of superposition due to high firing rate. A: Amplitude of spike triggered averaging (STA) motor unit action potentials (MUAPs) as a function of true MUAPs. B: Normalized amplitude errors of the STA estimates. C: Change of coefficient of variation (CV) of amplitude of individual and averaged across motor units (MUs). D: Change of correlation of individual and averaged across MUs. E,F: Normalized amplitude error.

were used as the triggering signal for the STA calculation. The STA-derived amplitude, presented as a function of true amplitude at different rates of error injection (the number of injected firing events divided by the total number of firing events) is illustrated in Figure 9A. The STA calculation tended to over-estimate the MUAP amplitude of small MUs and tended to under-estimate the amplitude of large MUs. The biases became larger with higher rate of error injections. Similarly, the normalized amplitude errors (Fig. 9B) tended to be positive (i.e., over-estimation) for small MUs, and the errors reached up to 38%. The amplitude errors tended to be negative (i.e., under-estimation) for larger MUs, and the errors reached up to 29%. The change of CV of amplitude with inclusion of spurious events is shown in Figure 9C. There 564


was no consistent increment of CV with increasing rate of error injection. The averaged change of CV across units (thick line) reached up to 0.022 when 20% of the firing events were spurious errors. Similarly, the reduction of correlation (Fig. 9D) due to spurious error was small for most units, and the averaged reduction of correlation (thick line) was 0.028 when 25% of the firings were spurious errors. The correlations between the amplitude error and STA estimation features are shown in Figure 9E and 9F. A larger amplitude error tends to correspond to a higher CV of amplitude and a lower correlation coefficient. DISCUSSION

This study investigated several possible factors that can influence the reliability of STA estimated MUAPs. The 5 variables examined included MUAP MUSCLE & NERVE

October 2013

FIGURE 6. Influence of superposition due to synchronized firing. A: Amplitude of spike triggered averaging (STA) motor unit action potentials (MUAPs) as a function of true MUAPs. B: Normalized amplitude errors of the STA estimation. C: Change of coefficient of variation (CV) of amplitude of individual and averaged across motor units (MUs). D: Change of correlation of individual and averaged across MUs. E,F: Normalized amplitude error.

amplitude variation, MUAP duration variation, superposition of units due to high firing rates, synchronized firing of MUs, and spurious errors in action potential discrimination. Regarding the estimation bias of the STAderived amplitude, the results show that the amplitude variations had little influence on the estimation errors in the STA-derived MUAP amplitude; however, the duration variations led to systematic under-estimation of amplitude, which was unexpected. When the degree of MUAP superposition between units was varied using 2 different methods, high firing rates and synchronization, the superposition induced by high firing rates had an unexpected minimal influence on the accuracy of the STA-derived MUAP amplitude estimation; however, synchronized firing led to over-estimation of the amplitude. Spurious action potential Reliability Assessment of STA of sEMG

discrimination also led to estimation biases. The direction of bias depended on the MUAP amplitude in that under-estimated amplitudes were observed in large MUAPs, and over-estimated amplitudes were observed in small MUAPs. Regarding the shape variation in the STA-estimated MUAPs, the duration variations and high firing rates led to a higher degree of shape variation. Namely, a substantial increase in amplitude variation (CV of amplitude) and a consistent reduction of the correlation between the STAderived action potentials and the true action potentials were observed when duration noise was added to the MUAPs and firing rates were increased. The results also showed that the amplitude errors correlate with the shape variations (STA features), even though a factor that led to obvious amplitude errors may not necessarily lead MUSCLE & NERVE

October 2013

565

FIGURE 7. Influence of synchronized firing in realistic decomposed motor units (MUs). A: Synchronization index of all possible MU pairs of subject 1. Synchronization with MU itself was set at zero. B: Synchronization index of all possible MU pairs of subject 2. C: Cross-correlation histogram between MU 7 and 6 (top panel) and between MU 7 and 8 (bottom panel) of subject 1. D: Amplitude of spike triggered averaging (STA) motor unit action potentials (MUAPs) as a function of true MUAPs of the 2 subjects. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

FIGURE 8. Exemplar segment of firing events (spikes) and injected spurious firings (crosses). Only motor unit (MU) 1, 3, …, 13 are shown. 566


MUSCLE & NERVE

October 2013

FIGURE 9. Influence of spurious firings. A: spike triggered averaging (STA) -estimated amplitude of motor unit action potentials (MUAPs) as a function of the amplitude of the true MUAPs. B: Normalized amplitude errors of the STA-estimated MUAPs at different firing rate errors. C: Change of coefficient of variation (CV) of individual and averaged across motor units (MUs). D: Change of correlation of individual and averaged across MUs. E,F: Normalized amplitude error.

to substantial action potential shape variations. Overall, this study provides information about the factors that influence the reliability of the STA-estimated MUAPs and offers potential calculation procedures that can assess the accuracy of the STA estimates and possibly exclude the unreliably estimated MUAPs. Amplitude and Duration Variations. The amplitude and duration of an MUAP train can vary due to multiple factors. During EMG recording, the variation may arise from the position shift of the recording electrodes with respect to the muscle fibers. There may also be impedance changes of Reliability Assessment of STA of sEMG

the electrode-electrolyte interface, possibly due to sweat or electrical contact changes.25 Additional variation may also arise from the fact that not all muscle fibers within an MU are activated simultaneously by the motoneuron, and that the muscle fiber conduction velocity can fluctuate.26 Focal temperature changes can also modulate the MUAP duration27,28 and sometimes can induce MUAP amplitude drift.29 During the STA calculation, approximately 300–500 firing events for the entire simulation duration and approximately 60–100 firings for a segment of EMG were used to estimate the MUAP. Given that the added amplitude noise had a zero MUSCLE & NERVE

October 2013

567

mean, the fluctuation of amplitude can be canceled out when a sufficient number of firing events were used to estimate the action potential shape. Therefore, it seems that the STA can tolerate a spike-by-spike amplitude variation up to 40% and still derive a consistent and accurate action potential shape if the trial duration is at least 10 seconds long. Given that duration variation influences the degree of superposition between action potentials and that the superposition effect can be canceled out when a large number of firing events are used, one would expect that duration noise has a minimal influence on the STA estimates. However, duration variation leads to systematic underestimation of MUAP amplitude. Essentially, when the multi-phasic action potentials of a particular MU are aligned based on the firing, the duration noise can lead to phase misalignment between spikes of the same MU, and this is the case in both shortened and lengthened durations. Specifically, as shown in the exemplar action potentials embedded in Figure 4A, duration change can lead to summations of positive and negative phases of action potentials between spikes and result in amplitude cancelation. Because mismatched phases occur in both directions (lengthened or shortened), the STA-derived amplitude bias (under-estimates) increases consistently with the level of duration variation. We have simulated up to 40% of variation in amplitude and duration, which is expected to be much higher than what occurs in real experimental settings, because a MUAP with 40% change will likely to be classified as a different MU. Our results suggest that the estimation bias due to amplitude variation can be reduced by including a large number of firing spikes during the STA calculation. However, the control and assessment of estimation bias due to duration variation are not straightforward. An indirect way of assessing the biases is proposed later (see discussion of potential procedures to reduce bias – below). MUAP Superposition. The simulated mean interspike intervals are sometimes set to be as short as approximately 20 to 30 ms. Given that the recorded action potential durations are approximately 10 ms, substantial superposition of adjacent action potentials between different MUs is expected. However, the results show that the high firing rate has a minimal influence on the accuracy of STA estimates. It is possible that the superimposed action potentials may have different phases partially overlapped. When a large number of firing events were used during the STA calculation, the different phases of interfering action potentials 568


can evidently be canceled out. On the other hand, when a segment (4 s) of EMG was used to estimate the action potential shape, the limited number (60 to 100) of firing events in the segment of EMG is not sufficient to diminish the superposition effect; as a result, the variation of amplitude across segments increases substantially with firing rates. The findings suggest that an average of a large number of action potentials is necessary to cancel out the superposition effect and to derive reliable estimates of action potentials. Therefore, to reduce estimation bias, a long EMG recording duration (approximately 10 s) is recommended to minimize the influence of high firing rates. Unlike high firing rates, synchronized firings lead to substantial over-estimation errors. The absolute estimation errors were generally the same across MUs, so the relative errors of the small MUs were much larger than the relative errors of the large MUs (Fig. 6B). This is likely due to the fact that an action potential of a large size is superimposed on a small action potential, because a few large action potentials mixed with several small action potentials are likely to have more averaging bias than when a few small action potentials are mixed with several large action potentials. Therefore, interference is more severe on small MUs than on large ones. To reduce contamination of the STA due to synchronized firing, one could compute event correlations across MUs and eliminate the correlated firing events, provided that there are still a large number of firings left for the STA calculation. Spurious Action Potential Identification. An earlier study has shown that spurious action potentials can inflate the amplitude of small MUAPs.17 Our simulation extends these earlier findings and shows that the spurious spikes bias the estimate of the MUAP amplitude in different ways depending on the action potential size. Namely, the STA provides over-estimated amplitude in small MUs and under-estimated amplitude in large MUs. The over-estimation of small MUs is likely due to the fact that smaller MUs have a higher chance of including a larger action potential during the average calculation, because a small MU sits at the lower end of the rank-ordered MUAP amplitudes. By the same token, the inclusion of small action potentials during the STA estimation of large MUAPs leads to under-estimate amplitude. In this simulation, the spurious firing events of 1 MU are drawn from the firings of other MUs with equal probability regardless of the size of the MUAPs. Typically, the spurious action potential assignment occurs when 2 action potentials have similar shapes and amplitudes. The spurious MUSCLE & NERVE

October 2013

Clinical

examiner to extract single action potentials. As a result, the recorded MUAP amplitude is limited to small sizes. Additionally, multiple needle insertions are often required to obtain a representative sample of the motor unit pool, and these multiple insertions are painful. The application of our decompositionenhanced STA of surface EMG can overcome these limitations and extract a relatively large range of MUAPs nonrinvasively, even during high levels of muscle contractions. For example, in a recent preliminary study done in hemispheric stroke survivors using STA analysis of sEMG signals, we have found that the MUAP amplitude is reduced in paretic muscles of chronic stroke survivors (Fig. 3A).14 Indeed, the novel decomposition techniques combined with STA analysis can provide a powerful tool for the diagnosis and assessment of neural degenerative disorders. Nevertheless, caution should be taken during use of the STA technique in the estimation of MUAP morphological parameters. In particular, the potential confounding factors simulated in this study should be considered to reduce estimation bias.

Motor unit number estimation is typically used to assess motor unit loss in neural disorders (e.g., ALS5,11,30 and stroke13,31,32). In this technique, the STA of the EMG is routinely calculated to estimate mean MUAP size, and the maximum M-wave is evoked by supramaximal stimulation of the motor nerve.4 The motor unit number can then be obtained from the division of the M-wave by the mean MUAP size. If there is a systematic bias in either the MUAP size estimate or in M-wave measurement, the motor unit number estimation may be misleading. For example, the degree of synchronized firing between units is often elevated in these pathological populations.33 As shown in this study, synchrony leads to over-estimated MUAP amplitude, which, in turn, can lead to under-estimated motor unit number. Therefore, the potential bias from synchrony should be taken into account (by eliminating synchronized firings) to derive an accurate estimate. Alternatively, MUAP morphological parameters such as amplitude, duration, and polyphasia have also been used to assess neuromuscular changes in ALS and post stroke. For instance, reduced MUAP amplitude in paretic muscles has been used to signify motoneuron or muscle fiber loss,34,35 and increased MUAP amplitude or polyphasia of MUAPs have been used to assess compensatory reinnervation of muscle.15 In these cases, single MUAPs are extracted directly from intramuscular recorded EMG; however, this intramuscular method has critical limitations in that the force is necessarily constrained to low levels to allow the

Potential Procedures to Reduce Bias. Estimation bias due to superposition can be reduced readily by increasing recording duration and by eliminating the correlated firings from the triggering signal. However, the reduction of biases due to variations in action potential duration and spurious action potential detection is more uncertain, because the spike-by-spike variation level and the frequency of spurious firing errors are not known routinely. One indirect way of assessing the estimation bias is to examine the STA features systematically. For example, the amplitude and duration variability of MUAPs can be compared across differing segment lengths of the EMG record. Or the correlation coefficients of estimated MUAPs with the real MUAP signal can be assessed using different approaches to evaluate estimation bias. Our results show that the estimation errors correlate with the STA features (amplitude variation and action potential shape correlation) and provide ways to exclude the STA estimates with large biases. For example, selection criteria can be set on the STA features to extract only stable waveforms, and the MUAPs with larger estimation biases may be removed. Alternatively, as shown in the plots of estimation error vs. STA features, the MUAPs form clusters at low estimation errors (e.g., Fig. 6E and 6F). Classification analyses, such as cluster analysis, can be applied to the multi-dimensional STA features to select reliably estimated MUAPs.

firings could be drawn in a higher probability from MUs with similar size and shape; however, the specific parameters of the probability (e.g., the standard deviation of the distribution) are unknown, and an equal probability was used as an alternative. In a sense, the simulation generates the worst possible spurious errors that a decomposition algorithm or a human can make. Therefore, the observed estimation bias may not be as high as in the estimation due to realistic spurious firing errors. Even in this simulation procedure, the amplitude errors were more or less the same at error rates from 0 to 0.1 (Fig. 9B), beyond which the estimation errors start to show a systematic increment. The results suggest that the STA can tolerate approximately 10% spurious action potential errors without inducing estimation biases. Note that the 10% spurious firing errors represent the worst scenario, because there may be errors in missed firings that had no influence on the STA estimates. Implications

for

the

Use

of

Applications.


STA

in

MUSCLE & NERVE

October 2013

569

CONCLUSIONS

Overall, this study identified the sources of STA biases that can arise from physiological properties of the motor unit pool and from signal recording and processing procedures. The results suggest that duration variation in MUAP, synchronized firings, and spurious action potential detections can bias the estimates systematically. In general, the STA can still provide valid estimation if appropriate procedures are taken to eliminate unreliable estimates. REFERENCES 1. Fetz EE, Cheney PD. Postspike facilitation of forelimb muscle activity by primate corticomotoneuronal cells. J Neurophysiol 1980;44: 751–772. 2. Muir RB, Lemon RN. Corticospinal neurons with a special role in precision grip. Brain Res 1983;261:312–316. 3. Poliakov AV, Schieber MH. Multiple fragment statistical analysis of post-spike effects in spike-triggered averages of rectified EMG. J Neurosci Methods 1998;79:143–150. 4. Brown WF, Strong MJ, Snow R. Methods for estimating numbers of motor units in biceps-brachialis muscles and losses of motor units with aging. Muscle Nerve 1988;11:423–432. 5. McComas AJ. Invited review: motor unit estimation: methods, results, and present status. Muscle Nerve 1991;14:585–597. 6. Stalberg E. Propagation velocity in human muscle fibers in situ. Acta Physiol Scand Suppl 1966;287:1–112. 7. Farina D, Arendt-Nielsen L, Merletti R, Graven-Nielsen T. Assessment of single motor unit conduction velocity during sustained contractions of the tibialis anterior muscle with advanced spike triggered averaging. J Neurosci Methods 2002;115:1–12. 8. Milner-Brown HS, Stein RB, Yemm R. The orderly recruitment of human motor units during voluntary isometric contractions. J Physiol (Lond) 1973;230:359–370. 9. Roeleveld K, Stegeman DF, Falck B, Stalberg EV. Motor unit size estimation: confrontation of surface EMG with macro EMG. Electroencephalogr Clin Neurophysiol 1997;105:181–188. 10. Bromberg MB, Swoboda KJ, Lawson VH. Counting motor units in chronic motor neuropathies. Exp Neurol 2003;184(Suppl 1):S53– S57. 11. Bromberg MB, Forshew DA, Nau KL, Bromberg J, Simmons Z, Fries TJ. Motor unit number estimation, isometric strength, and electromyographic measures in amyotrophic lateral sclerosis. Muscle Nerve 1993;16:1213–1219. 12. Campbell MJ, McComas AJ, Petito F. Physiological changes in ageing muscles. J Neurol Neurosurg Psychiatry 1973;36:174–182. 13. McComas AJ, Sica RE, Upton AR, Aguilera N. Functional changes in motoneurones of hemiparetic patients. J Neurol Neurosurg Psychiatry 1973;36:183–193. 14. Hu X, Suresh AK, Li X, Rymer WZ, Suresh NL. Impaired motor unit control in paretic muscle post stroke assessed using surface

570


15. 16. 17. 18. 19. 20. 21. 22. 23. 24.

25. 26. 27. 28. 29. 30. 31. 32. 33.

34. 35.

electromyography: a preliminary report. San Diego, CA: IEEE; 2012. p 4116–4119. Lukacs M. Electrophysiological signs of changes in motor units after ischaemic stroke. Clin Neurophysiol 2005;116:1566–1570. Kallenberg LA, Hermens HJ. Motor unit properties of biceps brachii in chronic stroke patients assessed with high-density surface EMG. Muscle Nerve 2009;39:177–185. Bromberg MB, Abrams JL. Sources of error in the spike-triggered averaging method of motor unit number estimation (MUNE). Muscle Nerve 1995;18:1139–1146. Nawab SH, Chang SS, De Luca CJ. High-yield decomposition of surface EMG signals. Clin Neurophysiol 2010;121:1602–1615. Sears TA, Stagg D. Short-term synchronization of intercostal motoneurone activity. J Physiol (Lond) 1976;263:357–381. Datta AK, Stephens JA. Synchronization of motor unit activity during voluntary contraction in man. J Physiol 1990;422:397–419. Taylor AM, Steege JW, Enoka RM. Motor-unit synchronization alters spike-triggered average force in simulated contractions. J Neurophysiol 2002;88:265–276. Yao W, Fuglevand RJ, Enoka RM. Motor-unit synchronization increases EMG amplitude and decreases force steadiness of simulated contractions. J Neurophysiol 2000;83:441–452. De Luca CJ, Roy AM, Erim Z. Synchronization of motor-unit firings in several human muscles. J Neurophysiol 1993;70:2010–2023. Dietz V, Bischofberger E, Wita C, Freund HJ. Correlation between the dischanges of two simultaneously recorded motor units and physiological tremor. Electroencephalogr Clin Neurophysiol 1976;40:97– 105. Roy SH, De Luca G, Cheng MS, Johansson A, Gilmore LD, De Luca CJ. Electro-mechanical stability of surface EMG sensors. Med Biol Eng Comput 2007;45:447–457. Zwarts MJ. Evaluation of the estimation of muscle fiber conduction velocity. Surface versus needle method. Electroencephalogr Clin Neurophysiol 1989;73:544–548. Buchthal F, Guld C, Rosenfalck P. Action potential parameters in normal human muscle and their dependence on physical variables. Acta Physiol Scand 1954;32:200–218. Bertram MF, Nishida T, Minieka MM, Janssen I, Levy CE. Effects of temperature on motor unit action potentials during isometric contraction. Muscle Nerve 1995;18:1443–1446. Denys EH. The influence of temperature in clinical neurophysiology. Muscle Nerve 1991;14:795–811. Boe SG, Stashuk DW, Doherty TJ. Motor unit number estimates and quantitative motor unit analysis in healthy subjects and patients with amyotrophic lateral sclerosis. Muscle Nerve 2007;36:62–70. Li X, Wang YC, Suresh NL, Rymer WZ, Zhou P. Motor unit number reductions in paretic muscles of stroke survivors. IEEE Trans Inf Technol Biomed 2011;15:505–512. Hara Y, Masakado Y, Chino N. The physiological functional loss of single thenar motor units in the stroke patients: when does it occur? Does it progress? Clin Neurophysiol 2004;115:97–103. Enoka RM, Christou EA, Hunter SK, Kornatz KW, Semmler JG, Taylor AM, et al. Mechanisms that contribute to differences in motor performance between young and old adults. J Electromyogr Kinesiol 2003;13:1–12. Lukacs M, Vecsei L, Beniczky S. Large motor units are selectively affected following a stroke. Clin Neurophysiol 2008;119:2555–2558. Krarup C. Lower motor neuron involvement examined by quantitative electromyography in amyotrophic lateral sclerosis. Clin Neurophysiol 2011;122:414–422.

MUSCLE & NERVE

October 2013