Motor Unit Identification in the M Waves Recorded by High-Density Electromyogram

Objective: We describe and test the methodology supporting the identification of individual motor unit (MU) firings in the motor response (M wave) to percutaneous nerve stimulation recorded by surface high-density electromyography (HD-EMG) on synthetic and experimental data. Methods: A set of simulated voluntary contractions followed by 100 simulated M waves with a normal distribution (MU mean firing latency: 10 ms, Standard Deviation - SDLAT 0.1-1.3 ms) constituted the synthetic signals. In experimental condition, at least 52 progressively increasing M waves were elicited in the soleus muscle of 12 males, at rest (REST), and at 10% (C10) and 20% (C20) of maximal voluntary contraction (MVC). The MU decomposition filters were identified from 15-20 s long isometric plantar flexions performed at 10-70% of MVC and, afterwards, applied to M waves. Results: Synthetic signal analysis demonstrated high accuracy of MU identification in M waves (precision ≥ 85%). In experimental conditions 42.6 ± 11.2 MUs per participant were identified from voluntary contractions. When the MU filters were applied to the M wave recordings, 28.4 ± 14.3, 23.7 ± 14.9 and 20.2 ± 13.5 MU firings were identified in the maximal M waves, with individual MU firing latencies of 10.0 ± 2.8 (SDLAT: 1.2 ± 1.2), 9.6 ± 3.0 (SDLAT: 1.5 ± 1.3) and 10.1 ± 3.7 (SDLAT: 1.7 ± 1.6) ms in REST, C10 and C20 conditions, respectively. Conclusion and significance: We present evidence that supports the feasibility of identifying MU firings in M waves recorded by HD-EMG.


I. INTRODUCTION
T HE M wave is a short-latency compound muscle action potential (CMAP) usually evoked by electrical stimulation of the nerve trunk. It can be recorded by electromyographic (EMG) electrodes placed on the skin over the muscle of interest. The M wave is represented by a typical shape, where the duration, peak-to-peak amplitude and area under the curve are usually assessed [1], [2], [3], [4], and is employed in many research fields, from Sports Sciences to clinical practice [5], [6], [7], [8]. The M wave is sensitive to changes in muscle fiber membrane excitability, and it allows the assessment of the peripheral properties of the neuromuscular system without the involvement of the central nervous system. Thus, M wave has been extensively used as a measure of sarcolemmal excitability. For example, M wave has been used to determine the effects of fatiguing contractions on action potential neuromuscular propagation, showing decreased M wave conduction velocity, and increased amplitude, and duration [9], which have been linked to fatigability-induced regulation of Na+/K+ processes [10]. Moreover, M wave has been used to determine the success of stimulating the motoneuron pool when eliciting other types of responses, e.g., the Hoffman reflex [11], as well as when estimating the number of motor units (MUs) in patient populations [5], [7]. Finally, the maximal M wave might be used as a normalizing factor for surface EMG activity to account for peripheral influences on the signal amplitude [8].
The M wave represents the sum of motor unit action potentials (MUAPs), generated by MUs that are synchronously elicited by the electrical stimulus. However, in surface EMG recordings, the M wave is composed of MUAPs arising from MUs with different innervation zone locations and different conduction velocities [12], [13]. This results in a time dispersion of individual MUAPs in an M wave and can also lead to misinterpretations of the results due to partial M wave amplitude cancelation, where the positive phases of some MUAPs partially overlap with the negative phases of the other MUAPs. Thus, the M wave analysis at the MU level provides an important new investigation path that may increase our understanding of the neuromuscular mechanisms affecting the M wave [14].
Typically, MUAPs are obtained by single MU recordings with intramuscular electrodes (iEMG) that are highly selective [15], but limit the number of MUs that can be recorded to typically two or three MUs per person [16], [17], [18]. Such a low number of MUs might not be representative of the whole MU pool in a muscle and the analyses of individual MU firings in the M waves recorded by intramuscular EMG are lacking, at least to our knowledge.
With significant advances in surface high-density EMG (HD-EMG) [19], [20] and in EMG decomposition algorithms [21], [22], [23], identification of firings of a large number of MUs in isometric voluntary contractions became possible [24], [25]. In these conditions, MUs fire largely asynchronously, allowing for estimation of the HD-EMG mixing process and its mathematical inversion [21], [22], [23]. This results in so called MU filters [32] that support the separation of contributions from different MUs. But in M waves MUs firings are highly synchronized hindering the estimation of the HD-EMG mixing process directly from M wave recordings. Yet, in isometric conditions the HD-EMG mixing process and corresponding MU filters can be estimated from voluntary contractions and, afterwards, applied to the HD-EMG recordings of M wave. The aim of this study was to demonstrate the feasibility of this process and assess the accuracy of the M waves decomposition into contributions of individual MUs in simulated conditions with known MU firing patterns, as well as in experimental HD-EMG recordings of the soleus muscle.

A. Synthetic HD-EMG
Synthetic HD-EMG signals were generated by a multilayer cylindrical volume conductor model [26], comprising 40 mm thick bone, 30 mm thick muscle, 4 mm thick fat and a 1 mm thick skin layer. Ten fusiform muscles were simulated, with an elliptical cross-sectional area (depth and width radius set to 15 and 30 mm, respectively), an average fiber length of 130 mm and spread of the fiber endings of 5 mm. The innervation zone was positioned at the longitudinal center of the muscle fibers. In each muscle, 165656 fibers were first uniformly distributed across the muscle cross-section and 200 MUs with sizes ranging from 24 to 2408 fibers were distributed randomly and uniformly in the muscle tissue (Fig. 4). The MU size and recruitment threshold followed the Henneman's size principle [27], with many small and progressively fewer big MUs. MU conduction velocity was normally distributed with the mean value of 4.0 ± 0.3 m/s, whereat all the fibers belonging to the same MU shared the same conduction velocity.
For each muscle, 10 s long voluntary contractions at 10%, 30%, 50%, 70% and 90% of maximal voluntary contraction (MVC) were simulated, resulting in 105, 155, 178, 193 and 200 active MUs, respectively. Their firing patterns were generated as per the model of [28], with the parameters adapted to the biceps brachii muscle. The MU recruitment range was set to 0 -80% MVC. The MUs increased their firing rate linearly with a simulated contraction level from 8 Hz at recruitment to 35 Hz at the MVC level. For each MU, the coefficient of variation for the interspike interval was set to 20%.
A voluntary contraction was followed by 100 simulated M waves. In each simulated M wave, firing latencies of all the MUs were normally distributed with the mean set to 10 ms and Standard Deviation (SD LAT ) set to 1.3 ms, 0.7 ms, 0.3 ms and 0.1 ms, respectively (Fig. 1). To support unambiguous discrimination between the firing patterns of different MUs in the simulated M waves, we then removed one third of the firings of every MU, for each MU in randomly selected M waves. Of note, with SD LAT set to 0.1 ms and HD-EMG sampled at 2048 Hz, the Standard Deviation of all the simulated MU latencies per M wave was within 5% of the intersample interval, meaning that the simulated MUs shared a high portion of the firing times. Thus, had no simulated MU firings been removed from the simulated M wave, we would not have been able to discriminate between simulated firing patterns of different MUs in the M waves and determine which MU was identified by the M wave decomposition. Consequently, we would not have been able to reliably quantify the accuracy of the MU identification.
The HD-EMG recordings were simulated with an array of 10×9 electrodes. The interelectrode distance was set to 5 mm and the electrode diameter to 0.5 mm. For each simulated HD-EMG channel, we first convolved the channel specific MUAP with the simulated MU firing pattern. Afterwards, the resulting MUAP trains from all the MUs were summed up, yielding the simulated HD-EMG channel. Colored Gaussian noise (bandwidth 20-500 Hz) was added to the simulated HD-EMG signals with a signal-to-noise ratio (SNR) of 25 dB. The signals were generated at a sampling rate of 4096 Hz and, afterwards, downsampled to 2048 Hz. Thus, the simulated MU firings occurred also between the HD-EMG samples, replicating the conditions in experimental HD-EMG.

B. Experimental Recordings
Twelve healthy males (age: 33.6 ± 5.8 years) volunteered to take part in the study. The study procedures were in accordance with the Declaration of Helsinki and were approved by the Medical Ethics Committee of the Republic of Slovenia (n 0120-84/2020/4). The procedures were explained to the participants, who gave written informed consent prior to their participation.
The experiment started with the preparation of the skin by shaving and abrading it with abrasion paste. Subsequently, a semi-disposable adhesive electrode array (GR08MM1305, OT Bioelettronica, Torino, Italy, 5 × 13 electrodes with 8 mm interelectrode distance), was placed over the soleus. A reference electrode (5×3 cm, T3545, OT Bioelettronica, Torino, Italy) was placed on the tuberositas tibiae just under the patellar ligament, and the signal was grounded by a water-soaked strap (WS2, OT Bioelettronica, Torino, Italy) fixed around the ankle. The HD-EMG signals were amplified using a 16-bit amplifier (Quattrocento, OT Bioelettronica, Torino, Italy) and sampled with a frequency of 5120 Hz (OTBiolab+, OT Bioelettronica, Torino, Italy). Additionally, a pair of EMG electrodes (biEMG; interelectrode distance of 25 mm; Covidien 24mm, Walpole, MA, USA) were placed on the soleus, laterally to the HD-EMG electrode array, for live monitoring of muscle activity. The reference electrode (50×100 mm, 00734, Compex, Guildford, Surrey, United Kingdom) was placed over the patella. The biEMG signal was sampled at 4000 Hz using the PowerLab and LabChart data acquisition toolbox (both ADInstruments, Bella Vista, New South Wales, Australia). Both signals (HD-EMG and biEMG) were band-pass filtered (10 -500 Hz).
The participants were seated in a custom-made ankle dynamometer with their hips, knees and ankles flexed at 90°. Contractions were performed with the right leg, the foot of which was strapped over the metatarsal bone on a wooden foot plate connected to a force sensor. The lateral malleolus was aligned with the dynamometer axis of rotation, and a rigid blockade attached to the dynamometer was pressed on the thigh just above the knee joint, to restrict plantar flexion movement. The participants were instructed to maintain a steady and relaxed posture, with the forearms resting on the dynamometer support ( Fig. 2). Minimal extraneous movement was encouraged during the active phase of the procedures.
Each session started with a warm-up, consisting of seven 5-second-long submaximal isometric plantar flexions with intensities ranging from 50 to 100% of the perceived maximal force, followed by 30 seconds of rest. After an additional 60 seconds of rest, two MVC trials were performed (with 120 seconds rest in between). The participants were instructed to increase their effort to the maximal over 1-2 seconds and to maintain maximal effort for approximately 3 seconds. Loud verbal encouragement was provided to facilitate the validity of a MVC measurement. Afterwards, the participants performed 20 seconds long contractions at 10, 20, and 30% MVC, and 15 seconds long contractions at 40, 50, and 70% MVC. A rest period of 120 seconds was provided between submaximal contractions. The entire protocol of voluntary isometric contractions was repeated at the end of the experiment.
Stimulations were performed either at rest (REST), or whilst performing isometric plantar flexions at 10% (C10) or 20% (C20) of MVC. Single rectangular electrical impulses (1 ms long [29]) were delivered to the tibial nerve via a custom-built current driven electrical stimulator (NeoStim1, EMF-Furlan, Ljubljana, Slovenia). The optimal stimulation position in the popliteal fossa was first assessed with a stimulation pen, which was then substituted with a self-adhesive cathode (Covidien 24mm, Walpole, Massachusetts, USA) fixed additionally with an elastic band. Anode (50×90 mm: MyoTrode PLUS, Globus, Italy) was placed on the patella. Afterwards, stimulation intensities started at 10 mA and were increased by 2 mA until it was clear that the long-latency responses (H-reflex) had peaked, i.e., the response amplitude was lower compared to the observed peak amplitude for four consecutive stimuli. After that, the intensity was increased in 5 mA increments until no further increase in peak-to-peak M wave amplitude was observed despite additional increments in electrical current. At each stimulation intensity, four responses were evoked at a frequency of 0.2 Hz. At least 52 M waves were elicited for each level of background muscle activity (REST, C10 and C20). 80 ms long epochs of biEMG signals were plotted after each electrical stimuli, allowing realtime visual feedback and quality control of the biEMG response, but were not further processed in this study.

C. Data Analysis
HD-EMG signals were modeled by convolutive data model introduced by Holobar and Zazula [21]: where h ij is L samples long MUAP of the j-th MU, as detected by the i-th HD-EMG electrode. MU spike train was defined as with δ denoting the unit-sample pulse and τ j (k) denoting the time of the k-th MU firing.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
The Convolution Kernel Compensation (CKC) technique [21] was used to identify the MU spike trains t j (n) from the HD-EMG of voluntary contractions. Each recording of a specific contraction level was decomposed independently. The Pulse-to-Noise Ratio (PNR) introduced and validated in [30] was used to assess the accuracy of each identified MU. Only the MUs with PNR ≥ 28 dB were kept for manual inspection by an expert. During the manual inspection, the segmentation of MU spike trains into MU firings as provided by the CKC technique was assessed and manually corrected, using the DEMUSE tool software (University of Maribor, Slovenia, version 6.0). The MUs with irregular firing patterns were discarded (firing rate < 8 Hz or coefficient of variation of interspike interval > 0.4).
In the synthetic HD-EMG, the identified MU firings were compared to the simulated ones and the precision and sensitivity of MU firing identification were calculated [31]: where TP, FP and FN denote the number of true positive, false positive and false negative MU firings, respectively, with a tolerance of MU firing match set to 0.5 ms.
In the experimental conditions the HD-EMG signals from different voluntary contraction levels were concatenated, and the MU filter was calculated for each identified MU at every contraction level, using the procedure described in [32]: is a vector of HD-EMG channels, extended by factor F = 15 [21], C −1 Y is the matrix inverse of the correlation matrix of Y, and n p denotes the firing time of a given MU in a voluntary contraction.
Afterwards the MU filters were applied to the concatenated HD-EMG signals Y(n), identifying the MU spike trains t(n) in all the voluntary contractions [31]: The spike trains were inspected manually and segmented into MU firing patterns. As the same MU could be identified from different voluntary contraction levels, we compared the rates of agreement in MU firings (tolerance of MU firing match set to 0.5 ms), and all the MUs that matched in more than 30% of firings were treated as duplicates of the same MU. Only one of them (the one with the highest PNR value) was kept.
Afterwards, the MU filters were recalculated as in (4), taking the longest possible temporal support of the MU spike train where the individual MU exhibited PNR ≥ 30 dB. The remaining part of the MU spike train was not used for MU filter recalculation. Finally, using the (5), the MU filters were applied to the extended HD-EMG signals simulated or recorded during the M waves, yielding the MU spike trains in the elicited contractions. The spike trains were segmented into firing patterns using the segmentation threshold computed from voluntary contractions by the CKC technique. In C10 and C20 conditions, only spikes appearing at least 5 ms and at most 25 ms after the stimuli were considered as potential MU firings in M waves. All the remaining MU firings were considered to represent voluntary muscle activation and were ignored. Afterwards, the segmentation results were inspected and edited by an expert, who evaluated the presence of MU crosstalk in the identified MU spike train t(n) (i.e., the level of contributions of other MUs to a given MU spike train, as suggested in [30]) during both voluntary and elicited contractions, and corrected the results of the automatic spike train segmentation.
To test the robustness of the MU filter transfer and identify potential changes of MUAP shapes in experimental conditions, we also applied MU filters to the HD-EMG recordings obtained from the voluntary contractions after the M wave protocol, identifying the spike trains of the MUs that were active before the M wave protocol. Again, the spike trains were segmented into MU firing patterns using the segmentation threshold calculated by the CKC method [21]. When the same MU was identified in more than one voluntary contraction after the M wave protocol, the spike train with the highest PNR value was considered. Afterwards, we pairwise compared the PNR values of MU spike trains before and after the M wave protocol.
In the experimental recordings, we analyzed the number of MUs identified per induced M wave, the latencies of MU firings with respect to the time of stimuli, the Standard Deviations of the MU firing latencies calculated across the M waves for each MU, and the Standard Deviations of the MU firing latencies calculated across all the MUs active in each M wave.
Because the CKC method identifies the MU spike trains with the spikes positioned close to the MUAP peaks in the HD-EMG [21], all the spikes of the given MU get delayed for the same amount of time (for a few ms). As different MUs exhibit different MUAP shapes, these delays differ among the MUs. To compensate for this methodologically induced dispersion of MU firing latencies, we estimated MUAPs by spike triggered averaging of HD-EMG signals from concatenated voluntary contractions, using the identified MU firings as triggers. We then used custom-made MATLAB script to plot the multichannel MUAPs of each MU along with the relative position of the MU spikes used as triggers and visually estimated the delay of identified MU spikes with respect to the earliest detected start of a MUAP in the HD-EMG (i.e., the earliest deflection of MUAP waveform from baseline). We subtracted this delay from the MU latencies in the M wave recordings. Afterwards, we applied the 2D Laplacian filter (LP) to the HD-EMG signals to remove the contributions of distant MUs to individual HD-EMG channels. Noteworthy, unlike in intramuscular signals, not all the detected MUs get identified by HD-EMG decomposition techniques. For example, the identified MUs account for 20-50% of energy from decomposition of voluntary contraction, with the remaining energy representing unidentifiable small and/or distant MUs [25]. In order to quantify the representativeness of MU identification in M waves in both monopolar and LP-filtered HD-EMG signals, we calculated the energy accounted for (EAF) by M wave decomposition [25]: where E is mathematical expectation and y i (n) is the i-th HD-EMG channel. m ij (n) is the MUAP train of the j-th MU in the i-th HD-EMG channel, obtained by convolving estimated MUAPs with the identified MU spike train, taking the subsample arithmetic into account to align the MUAPs in m ij (n) with the MUAPs in y i (n).
Finally, to further test the sensitivity and specificity of the decomposition in experimental conditions, we introduced perturbations to the HD-EMG signals and MU filters, respectively. Both perturbations have been applied for each MU individually, i.e., in each test only one MU and its filter were affected at a time, whereas the tests were sequentially conducted across all the identified MUs. The first and relatively small perturbation aimed to test the sensitivity of identified spikes to small changes of MUAP shapes. For this reason, we added a zero mean random shift of MU firing times n p in (4) with standard deviation set to 0.5 ms. This resulted in limited changes of MUAPs as estimated by spike triggered averaging but in a relatively large changes of MU spike heights in identified MU spike trains (see Results section). The second perturbation subtracted the experimentally estimated MUAP train from HD-EMG signals, as in (6). Afterwards, we repeated the MU spike train estimation in (5) and compared the mean height of MU spikes in the identified spike train before and after the perturbation. For comparison reasons, we also quantified the MUAP changes due to perturbations, as estimated by spike triggered averaging of HD-EMG. Note that this is directly related to p Y T (n p ) factor in (4). In this analysis, the decrease of spike heights after the perturbation served as indication that only one MU (the perturbated one) contributed to the identified MU spikes, thus the identified spikes were not crosstalk from other MUs.
Statistical analysis was performed in MATLAB (version R2020a). A Lilliefors test was used to check for normality of the data distribution. The test rejected the normal distribution of the reported results. In addition, the sets of MUs identified differed across different experimental conditions. Therefore, Kruskal-Wallis test was used for unpaired comparisons, and a Friedman test was used for paired comparison. Tukey's honestly significant difference was used when significant differences were detected. For all the statistical comparisons the level of significance was set at p < 0.05. As can be seen in Fig. 4, the identified MU firings in M waves were mainly those belonging to MUs closer to the surface, of larger sizes and greater conduction velocities.

A. Synthetic HD-EMG
When analyzing the precision and sensitivity of MU firing identification in the M waves, MUs were clustered based on recruitment threshold (0-30, 30-60, 60-80% MVC; Fig. 5). In   4. Comparison of positions, sizes, and conduction velocities of MUs whose firings were (red) and were not (blue) identified from synthetic M waves with SD LAT set to 0.7 ms. Bottom right panel depicts the muscle cross-section and MU territories (circles) as accumulated across all ten simulated muscles. In total, 164 (8.2%) out of 2000 simulated MUs were identified from all ten muscles. agreement with the previously presented results on voluntary contractions [31], the low-threshold and small MUs were difficult to identify accurately compared to higher-threshold MUs across all synchronization levels (Fig. 5). Indeed, the MUs with recruitment threshold between 60 and 80% MVC demonstrated higher sensitivities (p < 0.01) and precisions (p<0.03) than MUs with recruitment thresholds < 30% MVC (Fig. 5), regardless of the MU synchronization level.
In monopolar HD-EMG, the identified MUs accounted for 15 ± 5%, 20 ± 6%, 25 ± 6% and 28 ± 5% of M wave energy for SD LAT of 0.1 ms, 0.3 ms, 0.7 ms and 1.3 ms, respectively. When M waves were filtered by Laplacian filter, these figures Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. increased to 36 ± 7%, 44 ± 11%, 49 ± 10% and 52 ± 9%, which was significantly greater compared to monopolar (Friedman's test: p < 10 -5 ). The remaining energy was contributed by small and distant MUs (Fig. 4). In voluntary contractions, the identified MUs accounted for 37 ± 6% and 54 ± 6% of energy of monopolar and LP-filtered HD-EMG signals. These results are comparable to the ones reported in other voluntary muscle contractions [25].

B. Experimental Recordings
After concatenating the results of the voluntary contraction levels before the M wave protocol and after eliminating MU duplicates, we identified 42.6 ± 11.2 (range 24 -64) unique MUs, and, thus, MU filters per participant. The PNR values of the spike trains used for MU filter calculation were 33.7 ± 3.2 dB. When applied to the M waves, MU filters yielded the spike trains with the PNR values of 31.5 ± 8.5 dB, whereas in voluntary contractions after the M wave protocol, the identified MU spike trains exhibited the PNR of 32.6 ± 3.6 dB, which is above the recommended value of 30 dB [30], and corresponds to a > 90% accuracy of MU firing identification. All the MUs identified before the M wave protocol were also identified after the M wave protocol. An example of the identified MU firing patterns in representative isometric voluntary contractions before (top panel) and during the M wave protocol (bottom panel) is depicted in Fig. 6.
Examples of the identified MU firings in three consecutive M waves elicited in a single subject at REST (top row), C10 (central row) and C20 (bottom row) conditions are depicted in Fig. 7. Although each depicted M wave was decomposed individually by applying the MU filters learned from the voluntary contractions, the MU firing patterns in different M waves exhibit relatively high consistency (see Fig. 8 for visual comparison of two superimposed M waves).
For the purpose of statistical analyses, we considered the average of the first four (minimal stimulation intensity) and the last four (maximal stimulation intensity) elicited M waves in all the participants, respectively. The minimal and maximal stimulation intensities were chosen for this analysis to standardize the electrical input used to trigger the responses. Physiologically, the responses to the minimal and maximal intensities indicate the motor threshold and maximal muscle response, respectively. The maximal root-mean-square (RMS) values of the elicited M waves as calculated across the LPfiltered HD-EMG channels did not differ across REST, C10 and C20 for M waves elicited with minimal (p > 0.40, Fig. 9(a) and maximal stimulation intensities (p > 0.96, Fig. 9(b)). Similarly, the corresponding number of identified MU firings in the M waves elicited by minimal (Fig 9(c)) and maximal (Fig  9(d)) stimulation intensities did not differ across conditions (p > 0.05).  In monopolar HD-EMG signals, the identified MUs accounted for 11 ± 9%, 10 ± 6% and 8 ± 5% of the M wave energy in REST, C10 and C20, respectively. When M waves were filtered by LP filter, the energy accounted for increased to 30 ± 14%, 25 ± 11% and 18 ± 6%, respectively (p < 10 -5 ). The reasons for lower values in C10 and C20 conditions are likely due to the employed spike segmentation procedure, where any MU firing appearing outside the interval from 5 to 25 ms after the stimuli were considered to represent the voluntary muscle activity and, thus, ignored in our analysis. In voluntary contractions, the identified MUs accounted for 30 ± 9% and 32 ± 10% of energy for monopolar and LP-filtered HD-EMG signals.
In order to indicate better the relation between the global M wave metrics and the individual MU behavior, we do not report the latencies from the stimuli to the start of the M wave, as they were not significantly different across the REST, C10 and C20 conditions. Instead, following the previously employed methodology [33], the M wave latencies were measured from the stimuli to the maximal positive peak of each M wave, whereat the maximums were calculated across all the LP-filtered HD-EMG channels. The latencies of the peak of the LP-filtered M waves were smaller in C10 than in REST and C20 condition (p < 0.02, Fig. 10(a)). Concomitantly, the MU firing latencies were the smallest in C10 condition (p < 0.01, Fig 10(b)). The Standard Deviations of the MU firing latencies across all the firings of individual MUs were the smallest in REST condition (Fig. 10(c), p < 0.02), whereas the Standard Deviations of the MU firing latencies across all the MU firings per M wave were not significantly different in REST, C10 and C20 conditions (p ≥ 0.07, Fig. 10(d)).
Finally, the comparison of spike heights before and after perturbations are depicted in Fig. 11. All the identified MUs Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. ) and across all the MU firings per M wave (panel (d)) in REST, C10 and C20 conditions, respectively. In panel (a), for the elicited M wave, the LP-filtered HD-EMG channel with the maximal P2P amplitude of M wave was determined, respectively, and the latencies were calculated from the stimuli to the maximum positive peak. The horizontal black lines denote the statistically significant differences as revealed by the Kruskal-Wallis test, followed by the Tukey's honestly significant difference procedure (p<0.05).
during voluntary contraction and 97.7% of MUs during M waves decreased their spikes after MU filter perturbation, even though the MUAP shapes as estimated by spike triggered averaging of HD-EMG after the relatively small perturbation of MU firing times were often highly similar to the MUAP shapes before the perturbation (Fig. 11). After MUAP subtraction, all the MUs during voluntary contraction and 94.9% of MUs during M waves decreased their spikes. In both cases, the spike heights decreased significantly (Friedman test, p<10 -20 ). These results suggest that MU filters of the identified MU are specific to this individual MU and that the identified spikes are not the crosstalk from other MUs active in the M wave.

IV. DISCUSSION
In this study, we aimed to assess the feasibility and accuracy of the M wave decomposition into contributions of individual MUs in simulated conditions with known MU firing patterns, as well as in experimental HD-EMG recordings of the soleus muscle. The results suggest that the MU filters can be estimated by the CKC method [21] from the HD-EMG signals of voluntary contractions and applied to the HD-EMG recordings of the elicited M waves, identifying the firing times of individual MUs. As demonstrated by our results, the MU firings are not completely synchronized in percutaneous nerve stimulation (Fig. 7). These Fig. 11. Comparison of MUAP shapes (panels (a)-(c)) and peak-topeak (P2P) amplitudes (panel (f)) before (blue) and after the perturbation of MU filters (red) and MUAP subtraction (yellow) in voluntary contractions, respectively. As demonstrated in panels (c) and (f), subtraction of MUAPs was not perfect, due to errors in estimation of MUAPs by spike triggered averaging and errors in their alignment to HD-EMG signals. Namely, both spike-triggered averaging and MUAP alignment in HD-EMG signals are sensitive to MUAP superposition and to the activity of MUs that were not identified by CKC decomposition. As demonstrated in panels (d) and (e), the identified spike trains decrease significantly after MU filter perturbation and MUAP subtraction, respectively. differences in MU firing latencies are likely caused by differences in nerve and muscle fiber conduction velocities and in location and dispersion of innervation zones among MUs. When combined with the MU anatomical differences (e.g., location and size of MU territory and fibers length), these differences are obviously large enough to support identification of at least a portion of MU firings in M waves, recorded by HD-EMG.
In synthetic HD-EMG, the level of synchronization of simulated MUs was similar, or even higher, compared to that observed in the experimentally collected data (SD LAT ranged from 0.1 to 1.3 ms in synthetic and from 1.2 to 2.3 ms in experimentally collected data). MUs with a voluntary recruitment threshold of 30% MVC or lower (Fig. 5) were difficult to identify, especially at the two highest synchronization levels (SD LAT of 0.1 and 0.3 ms). At these recruitment thresholds, MUs are small, resulting in small MUAP amplitudes in the HD-EMG signals, and are thus more difficult to identify in both the evoked and voluntary contractions [31]. For the MUs with recruitment threshold > 30% MVC, the precision of MU firing identification was always higher than 90% (Fig. 5), demonstrating the high specificity of MU firing identification. However, the sensitivity was systematically lower than precision, ranging from 78.2% to 92.6% (Fig. 5). This indicates that some MU firings might be missed in the M wave decomposition, especially at high synchronization levels. Noteworthy, the largest number of MUs were identified from the M waves with the lowest simulated MU synchronization level (Fig. 5).
In experimental conditions the validation of decomposition accuracy is a challenge, as the true firing times of MUs are unknown. In voluntary contractions, an approach with both surface and intramuscular EMG recorded simultaneously and decomposed independently has been proposed for addressing this issue [34]. In this approach, MUs identified via both surface and intramuscular EMG are analyzed for agreement between firing times, and high agreement rates are considered to be a solid proof of decomposition accuracy. However, the number of MUs that are identified jointly from intramuscular and surface EMG has been demonstrated as relatively low [25]. Therefore, intramuscular EMG signals were not recorded in the present study and the validation of MU identification accuracy was done by means of indirect measures only.
First, we established the stability of the MU filters by controlling the changes of the MUAP shapes from two isometric voluntary sessions performed before and after the M wave protocol. Although the PNR values of the spike trains decreased from preto post-session, they remained above the recommended values of 30 dB [30], still guaranteeing the reliable identification of MU firings.
Secondly, the MU firing patterns detected during different M waves appeared to be highly consistent (Figs. 7 and 8). Indeed, when calculated across the M waves, SD LAT was 1.2 ± 1.2 ms, 1.5 ± 1.3 ms and 1.7 ± 1.6 ms at REST, C10 and C20, respectively (Fig. 10(c)). The SD LAT increased significantly with the voluntary contraction level (Fig. 10(c)), indicating an increase in the excitation of α motoneurons caused by voluntary contraction. This observation was further confirmed by the decrease of M wave peak latency in C10 condition compared to REST ( Fig. 10(a)) and by simultaneous decrease of MU latencies from REST to C10 condition ( Fig. 10(b)).
When measured across the identified MUs per M wave, the SD LAT were 2.1 ± 1.7 ms, 1.9 ± 1.6 ms and 2.3 ± 2.1 ms at REST, C10 and C20, respectively ( Fig. 10(d)). Whilst we did not find any published data on MU firings in M waves, the observed Standard Deviations of latencies are comparable to, or even smaller, than MU firings identified in the H-reflex with intramuscular EMG (individual MU mean 2.4 ms; across MU mean 3.4 ms; [35]).
After a relatively small perturbation of their MU filter, identified MU spikes significantly decreased their height, even though, when estimated by spike triggered averaging, the MUAP shapes after the perturbation were similar to the MUAP shapes before the perturbation (Fig. 11). The spike heights decreased significantly also after MUAP subtraction. These results complement the ones on synthetic HD-EMG, suggesting that also in experimental conditions MU filters were specific enough to reliably detect the spikes of individual MU. For at least 94.9% of tested MUs, identified spikes originated from perturbated MUs and not from the crosstalk of other MUs. Nevertheless, it also showed that the problem of MUAP similarities, and thus MU spike train crosstalk, increases in elicited contractions when compared to voluntary contractions (Fig. 11(d) and (e)). This is expected as both the level of MU synchronization and the number of activated MUs increase in elicited contractions. Noteworthy, the described permutation and MUAP subtraction methodology are signal and MU specific and allow us to identify and remove the problematic MUs from the M wave analysis.
Although decomposition identified high threshold MU firings at low stimulation intensities in some cases (Fig. 6 -bottom left), this pattern was not consistent across all the experimental conditions and all the participants. It has been previously shown indirectly (via twitch analysis; [36], [37]) and directly (via spike trigger averaging of intramuscular EMG signals; [38]), that electrical stimulation over the motor point of the muscle might result in reversal of the MU recruitment order. Indeed, electrical stimulation of a muscle is believed to first activate axons with large diameters, innervating high threshold MUs [39]. However, the recruitment order in response to direct stimulation of the motor nerve is less clear as revealed by the diverse results of other studies that used indirect measures to assess the MU recruitment order in nerve stimulations [36], [37], [40]. Our results suggest the MUs of the soleus muscle are not recruited in an orderly fashion in response to percutaneous nerve stimulations. However, further investigations are required to characterize all the factors that influence the MU recruitment order in response to nerve stimulation. This study aimed at demonstrating the feasibility of MU firing identification in M waves. Thus, the exact quantification of MU recruitment order is outside the scope of this study.
As the compound M wave increased progressively up to the maximal amplitude, the number of identified MUs increased significantly ( Fig. 9(c) and (d)). For LP filtered M waves in REST condition, the energy accounted for was 30 ± 14% and was comparable to 36 ± 7% of M wave energy accounted for when decomposing synthetic HD-EMG signals with SD LAT of 0.1 ms and lower than 52 ± 9% of M wave energy accounted for in synthetic HD-EMG signals with SD LAT of 1.3 ms. Although relatively low at the first glance, these numbers are comparable to the ones reported for voluntary contractions [25].
We offer the following explanations for these observations: First, the methodological procedure relies on transferring the MU filter learned from the voluntary contraction to the elicited contraction. MUs that are not reliably detected during the voluntary contractions cannot be tracked during the M waves. Therefore, the proposed method is limited to identification of MUs near the recording site (Fig. 4) and is likely not able to detect the MUs with the highest recruitment thresholds that are active during high contraction levels only (>70% of MVC in our experimental signals). Second, low threshold MUs tend to have smaller MUAPs in the HD-EMG, and are also difficult to detect in the presence of larger MUs [31] and high MU synchronization (Fig. 5). Third, the MUs that share similar MUAP shapes during voluntary contractions and get merged in estimated MU spike trains by the CKC method [31] are recognized (e.g., by low PNR value or/and irregular firing pattern) and removed in the process of voluntary MU firing detection. Thus, also these MUs cannot be tracked in M waves. Noteworthy, the number of active and identified MUs differs considerably among different muscles [31] and, therefore, also between synthetic and experimental conditions in our study. Fourth, MUAP estimation by spike triggered averaging and, consequently, MUAP alignment between y i (n) and m ij (n) in (4), have limited accuracy (Fig. 11) that depends on noise, background MU activity and the number of identified MU firings, among others. Collectively, these limitations are reflected in energies accounted for by HD-EMG decomposition and in Figs. 7 and 8, where the amplitudes of the sum of decomposed MUAPs are much lower than the amplitudes of the compound M wave. Fifth, the proposed methodology detects the MU firings close to the MUAP peak and not at the very beginning of the MUAP. This can be compensated partially by detecting the start of a MUAP manually, as revealed by the spike triggered averaging of HD-EMG, but this does not compensate for the potential delays due to MUAP propagation from the neuromuscular junction to the closest HD-EMG electrode. As the exact locations of the innervation zone and conduction velocities differ in different MUs, so may the delays of identified MU firings. This might have a significant impact on the estimated MU firing latencies ( Fig. 10(d)). Sixth, the experimental results of our study are presently limited to the soleus muscle, and thus replication of proposed analysis is needed in other muscles. Finally, in experimental conditions, no ground truth about the MU firing patterns is available. Therefore, the level of conclusiveness of provided evidence remains limited by the tests performed in our study. Hopefully, in the future, MU firing identification from much more selective intramuscular EMG will become available, supporting additional analysis of the proposed HD-EMG based approach.

V. CONCLUSION
We demonstrated that the identification of MU firings in electrically elicited M waves recorded by HD-EMG is feasible. Although the limitations presented herein suggest that the decomposition of electrically elicited M waves into the contribution of single MUs reflects only a portion of the whole muscle response to electrical stimulation of the nerve, several tens of MUs can be identified. Therefore, this methodology offers new possibilities of in vivo investigation of the human motor system. To the best of our knowledge, this is the first demonstration that the noninvasive MU firing identification in the M waves is possible.