On the Reuse of Motor Unit Filters in High Density Surface Electromyograms Recorded at Different Contraction Levels

We analyzed the efficiency of Motor Unit (MU) tracking across different experimentally recorded contraction levels of Biceps Brachii (BB), Tibialis Anterior (TA), First Dorsal Interosseous (FDI) and Abductor Digiti Minimi (ADM) muscles and in simulated conditions. We used the Convolution Kernel Compensation (CKC) algorithm to estimate the MU filters from high-density electromyograms (HDEMGs) at contraction levels ranging from 5 % to 70 % of Maximum Voluntary Contraction (MVC), and applied these MU filters to all the recorded contractions levels of the same muscle. For each MU filter estimation-application pair we assessed the number of identified MUs, Pulse-to-Noise Ratio (PNR), Mean Discharge Rate (MDR) and MUAP peak-to-peak (P2P) amplitude, normalized by HDEMG P2P amplitude. In simulated conditions, we also calculated the Precision ( $Pr$ ) and Sensitivity ( $Se$ ) of MU firing identification. We also reported the number of jointly identified MUs for all possible contraction level pairs. The results in simulated conditions confirmed the efficiency of MU filter transfer from low to high contraction levels. The large majority of MUs identified at lower contraction levels were tracked successfully at higher contraction levels, though many of them were not directly identified at higher contraction levels. Sensitivities of identified MU firings were above 97 %, decreasing only by about 1 % when transferring the MU filters from low to high contraction levels. In the experimental conditions the efficiency of the MU tracking depended significantly on the investigated muscle. The highest efficiency was demonstrated in the FDI muscle, where about 90 % and 70 % of MUs identified at the contraction level of MU filter estimation were also tracked at 10 % and 20 % higher contraction level, respectively. These figures decreased to 47 % and 20 % of MUs in TA, and to 21 % and 18 % of MUs in the BB muscle. The observed differences between the MU filter transfer efficiencies are likely caused by the differences in the size and anatomy of the investigated muscles.


I. INTRODUCTION
In humans, skeletal muscles consist of muscle fibers that are grouped into Motor Units (MUs). MUs differ in size and many other physiological parameters. Each of them is composed of several tens to several hundreds of muscle fibers, and the motor neuron that innervates all the fibers within that MU.
The associate editor coordinating the review of this manuscript and approving it for publication was Ganesh Naik .
Each electrical pulse in the motor neuron triggers the Single Fiber Action Potentials (SFAPs) in the innervated muscle fibers and, consequently, their contraction. Since the SFAPs of the fibers constituting the MU appear synchronously in time, they can be summed up conceptually into so-called Motor Unit Action Potential (MUAP). The MUAPs are large enough to be detected within a few centimeters from their origin, also on the surface of the skin above the investigated muscle, constituting a surface electromyogram (EMG). VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ The human motor system increases the muscle force by increasing the number of active MUs and the frequency of their activation, the so-called MU firing rate [1], [2]. In recent years, decomposition of a surface EMG into contributions of individual MUs enabled analysis of neural coding in relatively large populations of MUs, and, consequentially, many new insights into the neurophysiology of the motor system [1], [3]- [8]. However identification of individual MUs from a surface EMG is not a trivial task. Due to their low selectivity, surface electrodes detect overlapped contributions of a relatively high number of MUs. This results in highly complex EMG signals. Thus, computationally intensive computer algorithms are required to separate the contributions of MUs in the surface EMG and identify their firing patterns [9], [10]. Importantly, several tens of MUs are active in the detection volume of the surface electrodes, and not all of them are identified by decomposition algorithms [3]. Small and deep MUs contribute small MUAPs that are difficult to discriminate, and, therefore, constitute physiological noise in the process of MU identification [3], [4].
Two families of surface EMG decomposition algorithms have been introduced, namely template matching [5] and blind source separation (BSS) [4], [9]- [12]. While template matching algorithms aim to estimate the MUAP templates and match them to EMG waveforms, the BSS approaches invert the EMG mixing model and estimate the so-called MU filters. Each individual MU has a unique filter, which defines the weights in the linear spatio-temporal combination of EMG channels that yields the estimation of the individual MU spike train. For this reason, the BSS approaches require multichannel EMG recordings, the so-called High-Density EMG (HDEMG) recordings, but, compared to template matching algorithms, gain computational attractiveness when long EMG recordings need to be decomposed [3], [4]. Namely, under the constraint of isometric and nonfatiguing muscle contractions, MU filters can be estimated from the approximately first 10 s of the EMG recordings, and, afterwards, applied to several minutes long isometric recordings, yielding the complete MU firing pattern. This procedure of MU filter application is computationally very efficient, and may be conducted in real time [12], [13]. On the other hand, the process of MU filter estimation requires a few tens of seconds of processing time per MU on modern computers [9].
The contraction levels in the MU filter application phase may differ from the contraction levels in the MU filter estimation phase [3], [4]. This is highly beneficial for online MU tracking as well as for offline MU analysis. For example, short muscle contractions do not allow for reliable MU filter estimation. By using the MU filters that are pre-estimated on calibration contractions, we may benefit from statistical power of MU filters extracted from longer calibration contractions and transfer the filters to shorter contractions. This also greatly speeds up identification of MU spike trains at different contraction levels or in different experimental conditions (e.g., pre and post tested intervention) as searching for and refining of the MU filters requires more than 99 % of HDEMG decomposition time. Reuse of previously estimated filters may therefore greatly speed up the HDEMG decomposition process and neural code analysis. By combining the MU filters estimated at different contraction levels, we may also potentially increase the number of identified MUs from a given contraction level, and, therefore, improve the analysis of MU population coding.
However, there are several factors that can hinder the efficiency of such MU filter application. First, the length of muscle fibers may change with contraction level in isometric conditions [14] and the effect of this change on MU filter application likely depends on the muscle anatomy. Second, the number of unidentifiable MUs, i. e. the level of physiological noise, increases with the contraction level. In addition, following the Henemann's size principle [2] the size of MUs increases progressively with the force. For this reason, the lately activated MUs contribute the MUAPs that are significantly larger than the MUAPs of early recruited MUs, hindering the identification of small MUs at high contraction levels. Therefore, the efficiency of MU filter application across different muscle contraction levels may be muscle specific.
In this study, we first investigated the efficiency of MU filter estimation and application to different muscle contraction levels in simulated HDEMGs with different levels of noise. Afterwards, we analyzed the efficiency of MU filter application in experimental HDEMGs of four different muscles of different sizes and anatomies, namely, First Dorsal Interosseous (FDI), Abductor Digiti Minimi (ADM), Biceps Brachii (BB) and Tibialis Anterior (TA). Noteworthy, MU filter transfer analyzed in this study does not apply to fast fatiguing or dynamic muscle contractions with the extensive changes of MUAP shapes.

A. HDEMG MODEL AND MU FILTERS
In isometric contractions of skeletal muscles, HDEMG signals can be modeled by the following convolutive model [9]: where stands for time-wise extended vector of HDEMG signals, with extension factor F set between 10 and 20 [9], is the time-wise extended noise vector, and is similarly extended vector of MU spike trains. Here, the spike train of the j-th MU is defined as where δ() denotes the Delta function and k-th spike appears at time moment τ j (k). The mixing matrix H comprises L samples long MUAPs of J active MUs, as detected by M uptake electrodes [9].
The Convolution Kernel Compensation (CKC) method [9] estimates the MU filter iteratively aŝ where α determines the speed of convergence, E() stands for mathematical expectation, g(t) is a nonlinear weighting function, e. g. g(t) = log(1 + t 2 ) and C y = E(y(n)y T (n)) represents the correlation matrix of extended HDEMG measurements. After each iteration of Eqs. 6 and 7, the estimate of MU spike train gets updated MU filterf j may be estimated on one set of HDEMG signals (the MU filter estimation phase), and applied to another set of HDEMG signals (the MU filter application phase). Detailed descriptions of the CKC method are provided in [9] and [15].

B. SIMULATED HDEMG
To test the efficiency of MU filter transfer (estimation of MU filters at certain muscle contraction level and its application to other contraction levels), we simulated HDEMG signals of a 130 mm long, 60 mm wide and 30 mm deep fusiform muscle with an elliptical cross-section. 200 MUs, each comprising from 24 to 2048 fibers, were distributed randomly in cross-section area of the simulated muscle, with fiber density of 20 fibers/mm 2 [16]. The normally distributed MUAP conduction velocity of 4 ± 0.3 m/s was simulated, and the innervation zone was spread 5 mm around the fiber center. SFAPs were generated by the cylindrical volume conductor model [17] with the bone radius set to 25 mm, and with the skin and subcutaneous fat layer thicknesses set to 1 mm and 4 mm, respectively. Individual MUAPs were generated as a sum of corresponding SFAPs. Circular electrodes with diameter of 1 mm were simulated as being arranged into a 10 × 9 array. The electrode array was centred over the simulated innervation zone. The inter-electrode distance was set to 5 mm.
The MU firing pattern was generated as per [18], with the parameters adapted to the BB muscle. The recruitment firing rate was 8 pulses-per-second (pps) and increased with force up to 36 pps at maximum muscle excitation. The recruitment of MUs followed the Henneman's size principle [2].
The length of simulated contractions was set to 20 seconds, with the sampling rate of 2048 Hz. Six different constant muscle excitation levels were simulated, namely 10 %, 20 %, 30 %, 50 % 60 % and 70 %, yielding 105, 136, 155, 178, 186 and 193 active MUs. Ten muscles were simulated for  each contraction level. The MUs were distributed randomly within the muscle tissue in each simulated muscle.
Coloured Gaussian noise (frequency band of 20-500 Hz) with Signal-to-Noise Ratio (SNR) of ∞ dB, 20 dB and 10 dB was added to the HDEMG in three separate simulation runs.

C. EXPERIMENTAL HDEMG
Approximately 20 second long HDEMG signals were recorded in different studies [15], [19], [20] from neurologically intact young subjects performing isometric contractions of the FDI, ADM, TA and BB muscles at different Maximum Voluntary Contraction (MVC) levels (Table 1). Visual feedback on force was provided to the participants. All the experiments were conducted in accordance with the Declaration of Helsinki, and were approved by the local Ethical Committee. The participants gave informed consent prior to participation in the study. The details on the experimental protocol are provided in Table 1 and in [15], [19], [20]. Experimental HDEMG recordings are exemplified in Fig. 1.

D. DATA ANALYSIS
In both simulated and recorded HDEMG signals, we used the CKC method [9] to estimate the MU filters. This yielded several MU filters per each HDEMG, constituting one filter VOLUME 9, 2021 bank per each simulated run/participant, contraction level and SNR/muscle.
In some cases we had more than one trial available per participant and contraction level. In such cases, the trial was selected that yielded the highest number of MUs, with Pulseto-Noise Ratio (PNR) [15] above or equal to 25 dB and Mean Discharge Rate (MDR) above or equal to 7 Hz. This PNR threshold was selected as no manual editing of MU identification or MU tracking results was applied. To our experiences, the spike trains of MUs with PNR > 25 dB may often be manually edited to yield the PNR > 30 dB, as recommended in [15].
Each estimated MU filter bank was then applied to the HDEMG of each contraction level in the same muscle, yielding the corresponding MU spike trains (the MU filter application phase). Also in the application phase, the MU trains with either the PNR below 25 dB or the MDR below 7 Hz were discarded. Fig. 2 depicts an example of MU identification in the MU filter estimation and application phase, along with the HDEMG signals. Fig. 3 depicts representative examples of individual MU spike trains as automatically identified by the filter that was estimated from the HDEMG recordings at 30 % of MVC, and applied to HDEMG recordings at 30 %, 40 % and 50 % of MVC.
For each SNR/muscle, simulated run/participant and estimation-application MVC % pair we calculated the following metrics: The number of identified MUs, PNR, MDR and MUAP peak-to-peak (P2P) amplitude, normalized by HDEMG P2P amplitude. Additionally, we calculated the number of jointly identified MUs for all possible contraction levels pairs in each simulated SNR/muscle and in simulated HDEMGs, we also calculated the precision (Eq. 9) and sensitivity (Eq. 10) of MU firing identification, with MU firing detection tolerance set to 0.5 ms.
Here, TP, FP and FN stand for the number of true positives, false positives and true negatives, respectively.
The Shapiro-Wilk test with significance level α = 0.05 rejected the normal distribution of performance metrics in 45 % of the cases. Thus, we opted for the paired non-parametric Friedman test with significance level α = 0.05 to check for significant differences in the results. For each MU contraction level used for MU filter estimation, we compared the results across the filter application contraction levels. Only the cases with at least two identified MUs were taken into account for the PNR, MDR and normalized MUAP P2P amplitude comparisons. Furthermore, only the MUs jointly identified at both compared estimation-application contraction level pairs were tested due to the pairwise comparison.

A. SIMULATED HDEMG
The results on simulated HDEMG signals are depicted in Figs. 4, 6 and 7. Each row contains the results of MU filters' estimation at specified contraction levels and application to different contraction levels. In each cell, the MU filter application levels with significant differences from the cell's application level are listed in the brackets.
Both the precision and the sensitivity metrics were very high, averaging above 99 % in most of the cases. In a few cases, the precision and sensitivity of MU identification was significantly higher in the estimation than the application phase. In those cases, the precision and sensitivity metrics were still above 98 % and above 97 %, even with an SNR of 10 dB. For comparison, the number of MUs that were jointly identified by directly applying CKC to all possible contraction level pairs (no MU filter transfer between the contraction levels) are reported in Fig. 5. The individual cells display the number of MUs that were jointly identified by CKC at two different contraction levels. The number of jointly identified MUs decreased relatively fast with the increasing difference in muscle contraction levels, suggesting that the MU filters of early recruited MUs are difficult to identify at higher contraction levels. As revealed by the Friedman statistical test (α = 0.05), the identification of MU filter at lower contraction levels and its application to higher contraction levels (Fig. 4) identified significantly more low-threshold MUs at higher contraction levels that direct MU identification by CKC (Fig. 5). On the other hand, the transfer of MU filters from high to low contraction levels did not always identify more MUs than direct MU identification by CKC (Figs. 4  and 5).
The MUAP P2P amplitude of identified MUs, as normalized by HDEMG P2P amplitude is depicted in Fig. 8. Whereas in the MU filter estimation phase P2P MUAP amplitudes of identified MUs accounted for ∼33 % of HDEMG P2P amplitude, this ratio decreased significantly with the increasing contraction level in the MU filter application phase (Fig. 8).

B. EXPERIMENTAL HDEMG
The results for experimental HDEMG signals are depicted in Figs. 9 -13. The diagonal entries of Fig. 9 represent the cases where the application and estimation contraction levels are the same. In general, MU filter estimation yielded more MUs in FDI and TA than in the BB muscle. As indicated by the non-diagonal entries of the appropriate panels of Figs. 9 and 11, the FDI also had good MU filter transfer efficiency. As expected, MDR increased with contraction level in all the muscles studied (Fig. 12).
The number of jointly identified MUs by direct CKC application are reported in Fig. 10, with cells displaying the number of MUs that were jointly identified at two different contraction levels. Comparison by the Friedman statistical test (α = 0.05) revealed that the identification of MU filter at lower contraction levels and its application to higher contraction levels (Fig. 9) identified significantly more low-threshold MUs at higher contraction levels that direct MU identification by CKC (Fig. 10) in FDI and ADM, but this was not always the case in BB and TA muscle. Similar results were observed also when transferring MU filters from high to low contraction levels ( Figs. 9 and 10).
The MUAP P2P amplitudes of identified MUs, as normalized by HDEMG P2P amplitudes are depicted in Fig. 13. In comparison with simulated HDEMG (Fig. 8), the relative MUAP P2P amplitudes were significantly lower in the MU filter estimation phase. Nevertheless, the smallest MUAP P2P amplitudes were observed in FDI muscle, which demonstrated relatively high MU tracking efficiency across different contraction levels. This demonstrates that the efficiency of MU tracking cannot easily be assessed from MUAP amplitudes alone.

IV. DISCUSSION
The results on synthetic HDEMG confirmed the efficiency of the MU filter transfer from low to high contraction levels ( Figs. 4 and 6). The large majority of MUs identified at lower contraction levels were identified successfully at higher contraction levels (see above diagonal entries in Fig. 4), regardless of the simulated noise level. The PNR decreased with the contraction level (Fig. 6), but stayed well above 30 dB, suggesting an accuracy > 90 % in MU firing identification [15]. Indeed, in synthetic HDEMG, both the precision and the sensitivity of identified MU firing patterns were always above 97 %, decreasing only by about 1 % when transferring the MU filter from low to high contraction levels (e. g. at SNR of 20 dB, the precision of MU filter identification decreased from 99.9 ± 0.3 % to 99 ± 0.9 % when MU filters, estimated at 10 % MVC level, were applied to 70 % MVC level). In accordance with the simulated MDR, the MDR of the identified MUs increased with the muscle contraction level.
Transfer of MU filters allowed for identification of significantly more early recruited MUs at higher contraction levels than by direct CKC-based MU identification (Figs. 4 and 5). Due to the orderly recruitment [2] all the MUs identified at lower contraction levels should also be active at higher contraction levels, whereas vice versa is not true. The relatively rapid decrease of the number of identified MUs when  The number of jointly identified MUs by direct CKC application to simulated HDEMG with different levels of noise. Each cell depicts the number of MUs (mean ± standard deviation) that were jointly identified from a simulated HDEMG contraction level pair. In each cell, the contraction levels pairs with significant differences from the cell's contraction level pair are listed in the brackets. FIGURE 6. PNR in dB (mean ± standard deviation) in simulated HDEMG with different levels of noise. In each cell, the MU filter application levels with significant differences from the cell's application level are listed in the brackets. Empty (white) cells correspond to cases with fewer than two identified MUs.
transferring the MU filter from high to low contraction levels demonstrates the bias in the MU identification. At higher contraction levels, mostly high-recruitment threshold MUs were identified. Thus, in simulated conditions, MU filter estimation at lower contraction levels appears to be essential for identification of low-recruitment threshold MUs at high contraction levels.
The results on experimental HDEMG confirmed our hypothesis of variability in MU filter transfer efficiency across different muscles. The highest efficiency was achieved in the FDI muscle, where about 90 % and 70 % of MUs, identified at a contraction level of C % MVC was also identified at C + 10 % MVC and at C + 20 % MVC, respectively ( Fig. 9). At higher contraction level differences, the portion of identified MUs decreased significantly (∼45 % at C + 30 % MVC, ∼31 % at C + 40 % MVC and ∼7 % at C + 50 % MVC). PNR also decreased when transferring the MU filter from lower to higher contraction levels (above-diagonal entries in Fig. 11). On the other hand, when transferring the MU filter from higher to lower contraction levels, the PNR of identified MUs seldom decreased significantly (belowdiagonal entries in Fig. 11).   In the TA muscle, the efficiency of MU filter transfer was lower than in the FDI. About 47 % of MUs identified at contraction level of C % MVC were identified also at C + 10 % MVC, and about 20 % of MUs were identified also at C + 20 % MVC. At C + 30 % MVC and C + 40 % MVC this percentage decreased to 10 % and 2 %, respectively. VOLUME 9, 2021 FIGURE 10. The number of jointly identified MUs by direct CKC application to experimental HDEMG of different muscles. Each cell depicts the number of MUs (mean ± standard deviation) that were jointly identified from an experimental HDEMG contraction level pair. In each cell, the contraction levels pairs with significant differences from the cell's contraction level pair are listed in the brackets. FIGURE 11. PNR in dB (mean ± standard deviation) in experimental HDEMG of different muscles. In each cell, the MU filter application levels with significant differences from the cell's application level are listed in the brackets. Empty (white) cells correspond to cases with fewer than two identified MUs.
In accordance with these results, the PNRs of the identified MU spike trains were much lower in TA than in the FDI muscle, further demonstrating the increased complexity of MU identification in TA compared to the FDI muscle.
The lowest efficiency of MU transfer was demonstrated in the BB muscle. On average, only 21 % of MUs identified at a contraction level of C % MVC were also identified at C + 10 % MVC, about 18 % of MUs at C + 20 % MVC and 17 % of MUs at C + 30 % MVC. MU spike trains identified from the BB had the lowest PNR among all the investigated muscles (Fig. 11).
In ADM, the efficiency of the MU filter transfer was investigated on a relatively small range of contraction levels. Thus, the results can only be compared to similar experimental conditions in the TA and BB muscles (the lower three panels of Fig. 9). Higher efficiency of MU transfer was observed in ADM than in TA and BB. These intermuscular differences are likely caused by the differences in the size and anatomy of FIGURE 12. MDR in Hz (mean ± standard deviation) in experimental HDEMG of different muscles. In each cell, the MU filter application levels with significant differences from the cell's application level are listed in the brackets. Empty (white) cells correspond to cases with fewer than two identified MUs.

FIGURE 13
. MUAP P2P amplitude of identified MUs, normalized by HDEMG P2P amplitude (%, mean ± standard deviation) in experimental HDEMG signals of different muscles. In each cell, the MU filter application levels with significant differences from the cell's application level are listed in the brackets. Empty (white) cells correspond to cases with fewer than two identified MUs. the muscles. Fusiform muscles exhibit much lower diversity of MUAP shapes at the surface of the skin than pennate muscles. Furthermore, although not measured in our study, the subcutaneous tissue thickness was likely significantly lower in FDI and ADM than in the TA and BB muscles. Thus, the thickness of subcutaneous tissue could also be a major factor in MU filter transfer efficiency.
On average, direct CKC-based MU identification identified 2.5 ± 1.3, 2.1 ± 0.1, 0.2 ± 0.2 and 0.7 ± 0.4 fewer early recruited MUs at high contraction levels than the MU filter transfer in FDI, ADM, TA and BB muscle, respectively. Comparison of simulated and experimental conditions revealed significant discrepancies in the MU tracking performance. Although still the subject of further investigations, these discrepancies likely originate from the differences in complexity of simulated and real muscles, as well as from the differences in the MU recruitment strategies. With increase in muscle excitation, larger MUs get recruited [2] and the amplitude of their MUAPs at the surface of the skin depends on their location in the muscle tissue. Differences in distribution of larger and smaller MUs could partially explain the observed variations in MU tracking efficiency among different muscles as well as between simulated and experimental conditions. But as demonstrated by the results in Figs. 8 and 13 the relative P2P MUAP amplitudes alone cannot fully explain the observed discrepancies. Noteworthy, in experimental conditions the muscle geometry may change slightly with the contraction level [14], causing the changes in MUAP shapes on the surface of the skin. These changes are likely more pronounced in big than in small, and in fusiform than in pennate muscles, as suggested by the results of this study. Further investigations are required to confirm these claims.
In conclusion, the MU tracking across different muscle contraction levels is feasible, supporting the estimation of MU filters at low contraction levels and application of these filters to higher contraction levels. This may lead to a higher number of identified MUs at higher contraction levels, and may decrease the bias of the CKC method towards high-recruitment threshold MUs. But the efficiency of such MU filter transfer depends significantly on the anatomy of the investigated muscle, and likely also on the thickness of subcutaneous tissue layers. In dynamic and fatiguing muscle contractions, the identification and reuse of MU filters becomes much more challenging, due to the changes in the geometry of the muscle. Therefore further investigations of MU filter transfer efficiency are required. His research interests include biomedical signal processing and machine learning.
ALEŠ HOLOBAR (Member, IEEE) received the Ph.D. degree in computer science.
He is currently a Professor and the Head of the Institute of Computer Science at the Faculty of Electrical Engineering and Computer Science, University of Maribor, Slovenia. He has coauthored 75 articles in peer-reviewed journals, four book chapters, and more than 100 conference contributions on digital signal processing, human-machine interfaces, biomedical signal processing, and rehabilitation engineering.