Identification of Motor Unit Firings in H-Reflex of Soleus Muscle Recorded by High-Density Surface Electromyography

We developed and tested the methodology that supports the identification of individual motor unit (MU) firings from the Hoffman (or H) reflex recorded by surface high-density EMG (HD-EMG). Synthetic HD-EMG signals were constructed from simulated 10% to 90% of maximum voluntary contraction (MVC), followed by 100 simulated H-reflexes. In each H-reflex the MU firings were normally distributed with mean latency of 20 ms and standard deviations (SDLAT) ranging from 0.1 to 1.3 ms. Experimental H-reflexes were recorded from the soleus muscle of 12 men (33.6 ± 5.8 years) using HD-EMG array of ${5}\times {13}$ surface electrodes. Participants performed 15 to 20 s long voluntary plantarflexions with contraction levels ranging from 10% to 70% MVC. Afterwards, at least 60 H-reflexes were electrically elicited at three levels of background muscle activity: rest, 10% and 20% MVC. HD-EMGs of voluntary contractions were decomposed using the Convolution Kernel Compensation method to estimate the MU filters. When applied to HD-EMG signals with synthetic H reflexes, MU filters demonstrated high MU identification accuracy, especially for $\text {SD}_{\text {LAT}}>{0.3}$ ms. When applied to experimental H-reflex recordings, the MU filters identified 14.1 ± 12.1, 18.2 ± 12.1 and 20.8 ± 8.7 firings per H-reflex, with individual MU firing latencies of 35.9 ± 3.3, 35.1 ± 3.0 and 34.6 ± 3.3 ms for rest, 10% and 20% MVC background muscle activity, respectively. Standard deviation of MU latencies across experimental H-reflexes were 1.0 ± 0.8, 1.3 ± 1.1 and 1.5 ± 1.2 ms, in agreement with intramuscular EMG studies.


I. INTRODUCTION
T HE Hoffmann (or H) reflex is one of the most studied reflexes in humans and has been extensively used in physiological [1] and pathophysiological investigations [2]. The H-reflex is evoked by low-intensity electrical stimulation of the peripheral mixed nerve by depolarising mostly Ia afferent fibers from muscle spindles primary endings, which have monosynaptic excitatory projections to the heteronymous motoneuron (innervating the muscle from which the afferents emanate) and is, therefore, considered the electrical equivalent of the monosynaptic stretch reflex. It was widely believed that the magnitude of the H-reflex can be used to measure the motoneuron pool excitability [3]. However, since the H-reflex magnitude is sensitive to subtle changes in pre-and postsynaptic events and to the concomitant motoneuron activity, it cannot merely estimate the excitability of the motoneuron pool [4]. Nevertheless, its sensitivity allows studies of neural conductions across proximal segments and different facilitatory and inhibitory somatosensory responses, like the presynaptic, recurrent and reciprocal inhibition, which play an important role in motor control [4]. Thus, the H-reflex can be used as a probe to explore complex spinal neuronal pathways [4].
The H-reflex can be easily recorded using bipolar surface EMG; however, the H reflex is not a fully monosynaptic reflex and events from di-, oligo-, and poly-synaptic pathways can alter the magnitude of the reflex. For example, the stimulation of the mixed nerve might co-activate Ib afferents (from the Golgi tendon organ) and Renshaw cells. Due to a poor time resolution of the H-reflex, it has been suggested that single motor unit (MU) recordings are needed to correctly interpret the contribution of facilitatory and inhibitory pathways [5], [6]. Indeed, the effects of di-, oligo-and poly-synaptic contributions have been shown using single MU recordings, for example during vibratory stimuli [7] and active maintenance of standing [8].
Although single MU behavior can be recorded from surface electrodes [9], the behavior of individual MU is usually studied with intramuscular EMG (iEMG). Intramuscular electrodes are highly selective, allowing the analysis of individual MU action potentials (MUAPs) [10]. But this selectivity also represents a limitation as iEMG often reflects the activity of only a small number of MUs close to the recording site [11]. In fact, iEMG reflex studies are mainly based on analysis of 1-2 MUs per person [12], [13], [14], [15]. In addition, decomposition and interpretation of iEMG at higher forces is difficult [15], therefore the majority of iEMG studies are conducted at forces <10% of maximum voluntary contractions (MVC). This means that the iEMG reflex studies will likely be limited to a smaller proportion of the motor pool and to the MUs with lower recruitment thresholds.
In the last 20 years, there have been significant improvements in surface High-Density EMG (HD-EMG) technology [16], [17]. Along with the development of acquisition systems, decomposition algorithms also emerged supporting the identification of individual MU firing times from HD-EMG recorded in voluntary contractions [18], [19], [20], [21]. The study of Yavuz and colleagues [15] offered promising results on induced reflexes, demonstrating that the MU behavior after the electrically evoked reflex activity can be successfully identified. This provides an accurate quantification of reflex responses for a large number of individual MUs, but does not address the decomposition of the H-reflexes into contributions of individual MUs.
Indeed, the study of evoked compound muscle action potentials is challenging from the perspective of HD-EMG decomposition algorithms since the MU firings and, thus, MUAPs, are highly synchronized, hindering the performance of decomposition algorithms. In the evoked potentials, the level of MU firing synchronization largely exceeds those observed in voluntary contractions and likely also the ones seen in severe pathological tremor [22].
The H-reflex remains one of the most important techniques to assess spinal mechanisms in humans. Being able to decompose the HD-EMG of H-reflex and assess the contributions of a large number of individual MUs offers possibilities to investigate the complex interplay of di-, oligo-and poly-synaptic actions in human motoneurons. Therefore, the aim of this study was to propose the procedure that utilizes the information from decomposition of HD-EMG signals recorded during voluntary contractions (so called MU filters) and applies it to the recordings of elicited contractions to decompose the H-reflexes into contributions of individual MUs. Furthermore, we aimed to assess the reliability and accuracy of this proposed approach in simulated cases with known MU firing patterns as well as in experimental HD-EMG recordings of soleus muscle. We selected the soleus muscle since many of H-reflex studies have been conducted on this muscle [4].

II. METHODS
In this section we first define the HD-EMG data model used (Section II.A). We then describe the generation of synthetic HD-EMG with known MU firing patterns that support precise quantification of MU identification accuracy in both voluntary and elicited contractions (Section II.B). In Section II.C, we describe experimental protocol followed to record the soleus muscle. Finally, in Section II.D, we describe the proposed MU identification procedure and the measures used to assess its accuracy in both synthetic and experimental HD-EMG.

A. HD-EMG Model
HD-EMG signals were modeled by convolutive multipleinput-multiple-output system [23]: where, for the j -th MU, h i j is L samples long MUAP as detected by the i -th electrode and t j (n) denotes the MU spike train: with δ denoting the unit-sample pulse and the k-th firing of the MU appearing at time τ j (k). In voluntary contractions, the firing times of different MUs are largely asynchronous, but in elicited contractions, like H-reflexes, the external stimulus triggers the synchronous activation of several MUs.

B. Synthetic HD-EMG
To test the efficiency and accuracy of MU identification in controlled H-reflex conditions, we used the multilayer cylindrical volume conductor model [24] and generated synthetic HD-EMG signals with known MU firing patterns. The volume conductor that supports the simulation of fusiform muscles comprised 40 mm thick bone, 30 mm thick muscle, 4 mm thick fat and 1 mm thick skin layer [19]. Ten fusiform muscles with 200 MUs each were simulated [25]. Average fiber length was 130 mm with the fiber ending spread set to 5 mm [25]. In each simulated muscle, MU territories were randomly distributed in the elliptical muscle cross-sectional area of 1413 mm 2 (depth of 30 mm, width of 60 mm). The MU size (from 24 to 2408 fibers per MU [26]) and recruitment threshold (from 0 to 80% of MVC [26]) were distributed according to the Henneman's size principle [27], with many small MUs and exponentially fewer large MUs. Innervation zone was positioned at the center of the fiber length. Single fiber action potentials of the fibers belonging to the same MU shared the conduction velocity and were summed up to form the MUAP. Across simulated MUs, conduction velocity was normally distributed with the mean of 4.0 ± 0.3 m/s [25].
Two different contraction types were simulated in each muscle, namely 10 s long voluntary contraction at 10% (105 active MUs), 30% (155 MUs), 50% (178 MUs), 70% (193 MUs) and 90% of MVC (200 active MUs), respectively, followed by 100 elicited H-reflexes. During the voluntary contraction, the MU firing patterns were generated by the model proposed in [28], with the simulator parameters adapted to biceps brachii muscle. When recruited, MUs started firing at 8 Hz and linearly increased their firing rates with muscle excitation level to the maximum of 35 Hz [28]. In each simulated MU, the interspike variability was normally distributed with the coefficient of variation set to 20%.
In the phenomenological model of H-reflexes, 178 smallest MUs were simulated as active (orderly recruitment, equivalent to simulated voluntary contraction level of 50% of MVC). This level was selected to study both specificity and sensitivity of MU identification. In each H-reflex, the MU firings were normally distributed with mean latency of 20 ms and standard deviation (SD LAT ) set to 1.3 ms, 0.7 ms, 0.3 ms and 0.1 ms. Although such high synchronizations are not expected in experimental condition (and were also not observed in our study), they were used to test the ability of the proposed decomposition procedure to cope with very high levels of MU synchronization and to further reveal its limitations.
However, such a high synchronization rates pose a challenge to the estimation of MU firing identification accuracy. Namely, when all the firings of several MUs are highly synchronized, firing pattern of one simulated MU becomes almost identical to the firing patterns of the other MUs, which makes it difficult to determine which MU was truly identified by HD-EMG decomposition. To cope with this challenge and to generate unique MU firing pattens for each MU also in 100 simulated H-reflexes, we removed one-third of randomly selected MU firings in each simulated H-reflex. This resulted in 118 out of 200 simulated MUs being active in each H-reflex (∼60 % of all the MUs), comparable to the 50% of MUs observed active in the H-reflex of soleus activated by maximal stimulation of Ia afferents [29]. The generated MU firing patterns are exemplified in Fig. 1.
Two-dimensional array of 10 × 9 electrodes was simulated, with the interelectrode distance of 5 mm. At each electrode, the interferential HD-EMG signal was generated by convolving the MU firing patterns with corresponding MUAPs and summing up the contributions of different MUs. Colored (20-500 Hz) Gaussian noise with the signal to noise ratio of 25 dB was added to each simulated HD-EMG signal. All the HD-EMG signals were generated at sampling frequency of 4096 Hz and downsampled to 2048 Hz to simulate the occurrence of MU firings between the two HD-EMG samples.

C. Experimental Recordings
Twelve healthy recreationally trained males (age: 33.6 ± 5.8 years, mass: 79 ± 4.8 kg, height: 181± 4.9 cm) participated in this study. They were involved in moderate aerobic and anaerobic exercises twice a week. Exclusion criteria were acute injuries in the upper or lower extremities, locomotor dysfunctions, cardiovascular or neurological conditions. The protocol was approved by the Medical Ethics Committee of the Republic of Slovenia [n 0120-84/2020/4]. All procedures were performed according to the Declaration of Helsinki. Each participant provided written informed consent prior to the experiment.
Prior to applying the HD-EMG electrodes, the participant's skin was prepared: the hairs were shaved from the calf of the right leg; the skin was lightly abraded using an abrasive paste and cleansed. EMG activities of soleus muscle were recorded using semi-disposable adhesive matrix of 5 × 13 electrodes (GR08MM1305, OT Bioelettronica, Italy, interelectrode distance of 8 mm). The electrode array was mounted on the muscle belly covering the central portion of the soleus muscle, with the long side of the array aligned with the longitudinal axis of the muscle. This electrode position has been suggested as the position where the H reflex with the largest amplitude can be observed on the soleus muscle [30]. A water-soaked strap reference electrode (WS2, OT Bioelettronica, Italy) was fixated around the ankle. An additional ground electrode (5 × 3 cm, T3545, OT Bioelettronica, Italy) was mounted on the tuberositas tibiae just under the patelar ligament.
An additional bipolar EMG (biEMG) recording was used to monitor activity of soleus and to visually inspect the quality of the H-reflexes. Recording electrodes (Covidien 24mm, Walpole, USA) were mounted on soleus in a standard bipolar configuration at an interelectrode distance of 25 mm [31]. The electrodes were placed laterally to the HD-EMG array. The reference electrode (50×100 mm, 00734, Compex, Guildford, UK) was placed over the patella.
HD-EMG signals were collected with the Quattrocento acquisition system at 5120 Hz and 16-bit resolution and recorded with the OTBioLab+ software (both OT Bioelettronica, Italy). The biEMG signals were collected at 4000 Hz using the PowerLab toolbox and LabChart 8 software (both ADInstruments, Australia). Signals were band-pass filtered (10 -500 Hz) for both systems.
The participants sat in the ankle dynamometer (Wise Technologies, Ljubljana, Slovenia) equipped with a force sensor with their hips, knees and ankles flexed at 90 • ( Fig. 2A). The lateral malleolus of the tibia was aligned with the machine's axis of rotation and the plantarflexion movement was restricted using the machine fixation system, which uses blockade pressed on the thigh, just above the knee joint. The foot was additionally fixed using straps over the metatarsal bones. The participants were instructed to maintain a relaxed position with the forearm placed on the support of the dynamometer. To minimize the variability of the signal, they were instructed to maintain the head in the same position without clenching the teeth or squeezing their hands.
In a warm-up phase, participants contracted their plantar flexors seven times for 5 s (30 s between contractions) progressively increasing subjectively defined effort from 50% of MVC at the first trial to a quasi-maximal level at the last trial. After 60 s rest, they were verbally encouraged to contract their plantar flexors maintaining the MVC level for 3s. The MVC recording was repeated after 120 s rest. After another 120 s of rest, participants were instructed to increase and maintain torque for 20 s at 10%, 20%, 30% and for 15 s at 40%, 50% and 70% of MVC (120 s rest between contractions). The target torque intensity was visually displayed as a beam on the monitor placed in front of the participants. A second beam of a different color was modulated by the participants via exerting plantar flexion torque. The same procedure was repeated after the H-reflex protocol.
H-reflexes were elicited in the right soleus muscle by a custom built current driven high voltage electrical stimulator (TMG-S1, EMF-Furlan, Slovenia) delivering single rectangular electrical impulses (1 ms [32]) to the tibial nerve. A rather short 5 s interstimulus interval was used to minimize the duration of the sustained contraction avoiding fatigue and minimizing the inhibitory effect of post activation depression [33]. The anode (50 × 90 mm, MyoTrode PLUS, Globus, Italy) was placed over the patella. The optimal stimulation point of the tibial nerve in the popliteal fossa was assessed by a stimulation pen. The pen was then substituted for a self-adhesive electrode which served as the cathode (Covidien 24mm, Walpole, USA). Since the amplitude of the H-reflex is very sensitive to small changes in electrode and participant's position, we stimulated the tibial nerve with a stimulation intensity that caused a response with H-reflex on the ascending part of the HM curve and the M wave visible. The reproducibility of the M wave was used to monitor the stability of the stimulation [2]. A stable M-wave suggests that a constant number of motor nerve fibers and thus Ia afferents were excited by the same stimuli [34]. H-reflexes were evoked at three levels of background muscle activity: during rest (REST) and plantar flexion at 10% (C10) and 20% of MVC (C20), respectively. These conditions took place in a random order for each participant. At each contraction level, at least 60 H-reflexes were induced at a constant stimulation intensity. The biEMG was used to plot 80 ms long intervals of the signal after each electrical stimuli allowing real-time visual inspection of responses. Here, particular attention was paid to the amplitude of the small M wave preceding the H-reflex [35]. The recordings with an M wave amplitude exceeding the target by ±2 SDs were discarded. Apart from this, biEMG were not further analyzed in this study.

D. Data Analysis
Both the simulated and experimentally recorded HD-EMG during voluntary contractions were decomposed in monopolar mode, using previously described and extensively tested Convolution Kernel Compensation (CKC) technique [16], [17], [18], [19], [20], [21], [23]. HD-EMG signals from each contraction level were decomposed independently (Pseudocode 1). The quality of automatic decomposition was assessed by previously validated Pulse-to-Noise Ratio (PNR), a metric that efficiently assesses MU identification accuracy [19]: where E is mathematical expectation,t j (n) t j (n)=1 denotes the samples of MU spike train at the times of MU firings, whereast j (n) t j (n)=0 denotes the remaining samples of the MU spike train (the so called baseline noise). Following the recommendations in [19], all the MUs with PNR < 28 dB were discarded. The firing patterns of the remaining MUs were manually examined and edited by an expert to improve the quality of MU identification [36]. In this process, MUs with abnormal average firing rates (<8 Hz) or highly irregular firing patterns (coefficient of variation of interspike interval >0.4) were discarded. For each identified MU its MU filter was estimated by using Eq. (4) in Pseudocode 1. The MU filter represents the weights of the linear spatio-temporal HD-EMG channel combinations (Eqs. (4) and (5) in Pseudocode 1) yielding the estimation of individual MU spike train [37], [38]. Noteworthy, in isometric conditions, a MU filter can be estimated from HD-EMG recordings of one contraction and applied to HD-EMG recordings of another contraction of the same muscle, yielding the spike train of corresponding MU in the second contraction [37], [38]. In our study, the MU filters were estimated from HD-EMG recordings of voluntary is a vector of concatenated HD-EMG channels, extended by factor F = 15 or similar [23], C −1 Y is the inverse of the correlation matrix of Y, and n p denotes the firing time of a given MU in a voluntary contraction. Apply MU filter to the concatenated HD-EMG signals to estimate the MU spike train t (n) on the entire temporal support: contractions and, afterwards, applied to recorded H-reflexes (Pseudocode 1).
To prevent the identification of the same MU from different voluntary contraction levels, the MU firing patterns were mutually compared (Pseudocode 1, step 6). The duplicates of MUs with the smallest PNRs were discarded. For each identified MU, only the reliably estimated portion of the spike train, that is the largest time interval on which the spike train yielded PNR ≥ 30 dB, was used to recalculate the MU filter (Pseudocode 1, step 5), further improving its quality [38].
Afterwards, the MU filters were applied to HD-EMG recordings of H-reflexes, identifying the MU spike trains in elicited contractions (Fig. 3). MU spike train samples were segmented into MU firings or baseline noise (indicating no MU firing) using the previously introduced threshold-based spike segmentation of the CKC method [23]. Finally, the segmentation of identified spike trains was manually inspected by an expert, who visually evaluated and mutually compared the presence of MU crosstalk in the spike trains identified from voluntary and elicited contractions, and, when required, manually corrected the results of automatic spike train segmentation, as in previous studies [18], [19], [23], [39].
To test the robustness of MU filters in the experimental conditions, we applied the MU filters, estimated from voluntary contractions before the H-reflex protocol to the HD-EMG signals recorded during the voluntary contractions after the H-reflex protocol. We used threshold-based MU spike segmentation built into the CKC method [23] to identify the MU firings after the H-reflex protocol and calculated two separate PNRs for each tracked MU, one for spike train identified before and one for spike train identified after the H-reflex protocol.
In simulated conditions, the true MU firings are known. Therefore, we used precision and sensitivity, with the tolerance of firing identification set to ± 0.5 ms [38] to directly assess the accuracy of MU firing identification in H-reflexes: where TP, FP and FN stand for the number of true positive, false positive and false negative MU firings, respectively.
We also analyzed the number of identified MUs and their voluntary recruitment threshold (known from simulations). Noteworthy, we do not report the accuracy of MU firing identification during the voluntary contractions, as this has already been reported in previous studies [18], [19], [39]. In experimental conditions, the ground truth MU firings are not known. Thus, we used indirect measures of MU identification accuracy. First, we analyzed the number and the voluntary recruitment threshold of the MUs identified as active in H-reflexes. Second, we analyzed the MU firing latency (as measured from the corresponding stimuli) and two different standard deviations (SD) of MU latencies, one calculated across elicited H-reflexes per MU and one calculated across all the identified MU firings in each H-reflex. Third, we analyzed the peak-to-peak (P2P) amplitude of elicited H-reflexes and compared it to the number of identified MU firings. Fourth, we compared the latencies of the first positive peak of H-reflex (measured from the time of corresponding stimulus) to the MU firing latencies across different conditions (REST, C10 and C20). Although the H-reflex latency is usually measured from the stimuli to the start of the H-reflex, we decided to report the latency of the peak as it better reflects the trends noticed in MU latencies (see Fig. 8), and is consistent with the approach used previously to infer on the MU activity in elicited contractions [40].
It has been previously demonstrated on voluntary contractions [39], [41] that the HD-EMG decomposition identifies firing patterns of MUs that are relatively close to the uptake electrode, whereas more distant MUs constitute the physiological noise and cannot be identified. Due to the proposed concept of MU filter transfer from voluntary to elicited muscle contractions, this limitation also applies to the study of H-reflexes presented herein. To compensate for this and to better demonstrate the yield of H-reflex decomposition, the H-reflexes and MUAPs in HD-EMG signals depicted in the figures were filtered by 2D Laplacian (LP) spatial filter which removes the contributions of distant MUs and emphasizes the contributions of MUs close to the uptake electrodes. As the fifth indirect measure, we calculated the energy accounted for by HD-EMG decomposition (separately for voluntary and elicited contractions, and for monopolar and LP-filtered HD-EMG signals): where y i (n) stands for the i -th HD-EMG channel and m i j (n) denotes the estimated MUAP train of the j -th MU in the i -th HD-EMG channel. To calculate the m i j (n), we first estimated the MUAP of the j -th MU by spike-triggered averaging the i-th HD-EMG channel using the MU firings identified from voluntary contractions as triggers. We then convolved the MUAP with the identified firing pattern of the j -th MU. Another important limitation to consider was the methodological dispersion of the identified firings of different MUs. Namely, CKC method builds on convolutive mixing model and optimizes the decomposition results by identifying the MU spike train with the spikes positioned close to the MUAP peaks and not to the very start of the MUAP in HD-EMG [23]. In this optimization, all the spikes of the given MU spike train get delayed for the same amount of time (for a few ms). As different MUs exhibit different MUAP shapes [42], these methodological delays of spike trains differ among the identified MUs. To minimize this, we manually assessed the positions of the detected spikes for each MU with respect to the start of the earliest detected MUAP across all the monopolar HD-EMG channels and subtracted the estimated delay from the identified MU firing latencies in H-reflexes. As in the case of EAF calculation (7), MUAPs were estimated by spike triggered averaging of HD-EMG signals during voluntary contractions.
The results are reported as mean ± SD. Statistical analysis was performed in MATLAB (version R2020a). Lilliefors test was used to test for the normal distribution of results. When the normal distribution was rejected, Friedman's test was used for paired and Kruskal-Wallis test for unpaired comparison. When normal distribution was not rejected, ANOVA was used. Tukey's honestly significant difference procedure was applied when significant differences were detected. In all the statistical tests, the results were considered statistically significant when p ≤ 0.05. Fig. 1. 178 MUs with the recruitment threshold ≤50% of MVC were active in simulated H-reflexes. Their MU spike trains as identified by MU filters are exemplified in Fig. 3.

Example of simulated MU firings during voluntary contractions and H-reflexes is depicted in
For statistical analysis of the decomposition performance, MUs were clustered into five different groups, with clustering based on identification of MUs from simulated voluntary contractions at a given contraction level (10%, 30%, 50%, 70% and 90% of MVC), and collapsed across all ten simulated muscles. Several MUs with recruitment threshold <50% of MVC were identified also at 70% and 90% voluntary contractions, whereas none of the MUs with recruitment threshold >50% of MVC was identified as active during simulated H-reflexes (Fig. 4, top row). Except for MUs with recruitment threshold <10% MVC, precision and sensitivity of MU firing identification in H-reflexes were found to be very high for each of the simulated MU synchronization levels (Fig. 4, central  and bottom rows).

B. Experimental Recordings
We identified 42.6 ± 11.2 (range 24 -64) unique MUs per participant from HD-EMG signals recorded before the H-reflex protocol. By using their MU filters, we were able to identify the firing patterns of all these MUs also in the voluntary contractions after the H-reflex protocol (Fig. 5). The MU spike trains before and after the H-reflex protocol had PNR of 33.7 ± 3.2 dB and 32.6 ± 3.6 dB, respectively. Although both sets of PNRs were significantly different (Friedman's test: p<0.01), the PNR values were above the recommended  HD-EMG signals and MU firing patterns during the 30% MVC voluntary contraction (grey line) before (left) and after the H-reflex recording (right). Different colored bars denote the firing times of different MUs. For clarity reasons only three HD-EMG channels are depicted. We used the CKC method to identify MU firing patterns before the H-reflex protocol and estimate the MU filters. The MU filters of all the MUs identified before the H-reflex protocol were applied to the HD-EMG signals recorded after the H-reflex protocol.
value of 30 dB [19], suggesting > 90% accuracy of MU firing identification. Fig. 6 depicts the MU firing patterns during five consecutive H-reflexes in REST, C10 and C20, respectively. Although consecutive H-reflexes differ in the P2P amplitudes and in the  In Panels A and B the Lilliefors test rejected the normal distribution. Therefore, the Kruskal-Wallis test was used. In panels C and D the normal data distribution was not rejected. Therefore, we used repeated measures ANOVA. around 1 ms, 1.3 ms and 1.5 ms, in REST, C10 and C20, respectively (Fig. 8.D).
The P2P amplitude of the maximal H-reflex, calculated across all the LP-filtered HD-EMG channels for each H-reflex, was significantly larger in C20 than in REST and C10 (p < 10 −8 ; Fig. 7.A). In agreement with P2P amplitudes of H-reflexes, the number of identified MU firings was significantly higher in C20 than in other two conditions (p < 10 −6 ;  Fig 7.B). There were moderate (∼0.5) but significant (p<0.01) mean cross-correlations among the amplitudes of H-reflexes on different LP-filtered HD-EMG channels demonstrating that the elicited H-reflexes exhibited both temporal (within the same HD-EMG channel) and spatial (across HD-EMG channels) variability (Fig. 7.C). Conversely, the maximum crosscorrelation (pooled across LP-filtered HD-EMG channels) between the P2P H-reflex amplitude and the number of identified MU firings per H-reflex reached the values around 0.92 (Fig. 7.D), indicating that the number of identified MU firings per H reflex was significantly correlated with the H-reflex P2P amplitude.
Not all the MUs contributing to the elicited H-reflexes were identified. This was further demonstrated by EAF, where 30 ± 21 % and 30 ± 16 % of H-reflex energy was explained by identified MUs in monopolar and LP-filtered HD-EMG signals recorded during REST condition, respectively. In C10 condition, the EAF increased to 39 ± 15% and 39 ± 15%, respectively, whereas in C20 EAF was 37 ± 22 % and 41 ± 18 % in monopolar and LP-filtered HD-EMG signals, respectively. For comparison, in voluntary contractions, the identified MUs accounted for 30 ± 9% and 32 ± 10% of energy of monopolar and LP-filtered HD-EMG signals, respectively. This is in agreement with other studies, where small and/or distant MUs were typically difficult to identify [39], [41].
Employing the previously described methodology [40], the latencies of the first positive peak of the maximal H-reflex (maximum calculated across LP filtered HD-EMG channels per elicited H-reflex) decreased significantly with the voluntary contraction level (p ≤ 0.045; Fig. 8.A). Similarly, the mean MU latencies (the mean per MU as calculated across all the identified firings in H-reflexes) demonstrated a significant decrease with voluntary contraction level (p ≤ 0.05; Fig. 8.B). Note that latencies from the stimuli to the start of the H-reflex did not show significant differences in REST, C10, C20 (results not shown). Conversely, increases were observed in standard deviations of individual MU latencies as calculated across different H-reflexes (p < 0.035; Fig. 8.C) and across the firings of different MUs identified per individual H-reflex (p < 10 −5 ; Fig. 8.D).

IV. DISCUSSION
In this study we showed that the MU filters estimated by CKC method [23] from HD-EMG recordings of voluntary contractions when applied to HD-EMG recordings of elicited H-reflexes identify the firing times of individual MUs. We investigated the feasibility and accuracy of the H-reflex decomposition in simulated fusiform muscle with known MU firing patterns as well as in experimental HD-EMG recordings of soleus muscle, which is a pennate muscle. Noteworthy, in fusiform muscles MUs have parallel muscle fibers and more similar MUAP shapes than in pennate muscles. Therefore, HD-EMG signals from fusiform muscles are usually more difficult to decompose than the signals from the pennate muscles [37], [38]. By reporting on both muscle geometries, we provide a wide demonstration of the applicability of the proposed H-reflex decomposition. Nevertheless, the provided results on synthetic and experimental HD-EMG signals should be compared with caution.
Results on synthetic HD-EMG demonstrated high accuracy of MU firing identification also at very high MU synchronization levels. Indeed, the highest level of simulated MU synchronization was about one order higher than the one observed experimentally (SD LAT ranged from 0.1 to 1.3 ms in synthetic and from 2.7 to 3.3 ms in experimental conditions). These levels were selected to provide a better insight into the limitations of the proposed methodology. At the highest two synchronization levels (SD LAT of 0.1 and 0.3 ms) MUs with recruitment thresholds <10% MVC were not successfully identified (Fig. 4). The highest synchronization level (SD LAT of 0.1 ms) was also problematic for MUs with recruitment threshold <30% of MVC, as only two MUs were identified. These MUs have relatively small size and small MUAP amplitudes in HD-EMG signals. As such they are difficult to identify at 50% of MVC, in both voluntary and elicited contractions [39], [41]. This is further confirmed by statistical results presented in Fig. 4, where the highest number of MUs was identified by MU filters that were extracted from voluntary contractions comparable to the levels of muscle excitation during the elicited H-reflexes. At higher voluntary contraction levels, we were still able to extract the MU filters of lower threshold MUs, but their number decreased progressively with the decrease in their recruitment threshold (Fig. 4, top panels). This agrees with the results of MU filter transfer efficiency in voluntary contractions [38]. Of note, precision of MU identification in H-reflexes was relatively high. Sensitivities were consistently lower than precisions but were above 90% for MUs with recruitment threshold >30% MVC and simulated SD LAT > 0.3 ms. This indicates that a few MU firings might be missed in the H-reflex decomposition, but the false alarm rates (defined as 1-Precision) are relatively low (<1 % for MUs with recruitment threshold >30% MVC).
In experimental conditions assessment of H-reflex decomposition accuracy is a challenge as there is no ground truth about the MU firing times. The two source technique has been proposed in the past [16]. This technique simultaneously records both intramuscular and surface EMG and decomposes them independently. The rate of agreement between the firing times of MUs identified from both surface and intramuscular EMG is then considered to be a proof of accuracy. But it has been shown in previous studies that the number of MUs that are jointly identified from intramuscular and surface EMG is relatively low [19], [39]. In our study, we did not record the intramuscular EMG signals. Therefore, we assessed MU identification accuracy by different indirect measures.
First, we verified the stability of MU filter transfer from voluntary session before the H-reflex protocol to voluntary session after the H-reflex protocol. In this way we controlled for potential changes in MUAP shapes that could happen between both sessions. The PNR values of spike trains decreased from pre to post session but stayed above the recommended threshold of 30 dB [19], indicating reliable identification of MU firings (Fig. 5).
Next, we compared the P2P amplitudes of recorded H-reflexes (Fig. 7.A) to the number of identified MUs per H-reflex (Fig. 7.B) and showed that they are significantly correlated (Fig. 7.D). In fact, the correlation coefficient between the P2P H-reflex amplitude and the number of identified MUs per H-reflex was considerably higher that the correlation coefficient between the P2P H-reflex amplitudes on different HD-EMG channels (Figs. 7.C and 7.D).
Afterwards, we studied the consistency of MU firing patterns in elicited H-reflexes. When calculated across H-reflexes the standard deviations of individual MU latencies were 1.0 ± 0.8 ms, 1.3 ± 1.1 ms and 1.5 ± 1.2 ms at REST, C10 and C20 ( Fig. 8.C). This is comparable to 2.4 ± 1.4 ms measured by intramuscular recordings of electrically elicited H-reflexes in soleus [43]. When measured across identified MUs per H-reflex, the standard deviations of MU latencies were 2.7 ± 1.0 ms, 3.3 ± 1.7 ms and 3.3 ± 1.5 ms at REST, C10 and C20 (Fig. 8.D). This is comparable to 3.4 ms (range from 1.5 to 5.4 ms) obtained with intramuscular recordings by Burke and colleagues [43].
The mean MU firing latencies decreased with voluntary contraction level (Fig. 8.B), in agreement with the decreased latencies of H-reflex peak (Fig. 8.A). These changes can be explained by an increase in the composite excitatory postsynaptic potentials (EPSPs) of α motoneurons caused by voluntary contraction. When the EPSPs increase, the motoneurons get closer to the activation threshold and require lower level of additional electrical excitation to fire [2]. This is further confirmed by an increase of both peak-to-peak H-reflex amplitudes and the number of identified MU firings per H-reflex with voluntary contraction level (Figs. 7.A and 7.B). At the same time, the standard deviations of MU firing latencies increased with the level of voluntary contraction (Figs. 8.C and 8.D). This can again be explained by the facilitation of α motoneurons caused by the summation of EPSPs induced by the voluntary contraction and the electrically elicited stimuli.
Our methodology has several limitations. First, the presented procedure builds on transferring the MU filter from voluntary to elicited contractions. MUs that are not detected during the voluntary contractions cannot be tracked during the H-reflexes. Both voluntary contractions and H-reflex likely share the recruitment order and, therefore, the pool of activated MUs and this has been exploited by the methodology proposed herein. It remains to be verified to what extent the proposed methodology is also applicable to other evoked potentials, like M waves or motor evoked potentials (MEPs) and to different stimulation parameters, such as intensity and frequency of stimulation and stimuli duration.
Second, the effect of the MU crosstalk in different experimental conditions of H-reflex elicitations is yet to be investigated and understood. This could have significant impact on spike segmentation strategy. In this study, we relied on the segmentation strategy developed for voluntary contractions, which might be suboptimal for elicited contractions.
Third, HD-EMG decomposition identifies the MUs that are sufficiently large and relatively close to the measuring electrodes. MUs that are further away cannot be identified and, therefore, constitute physiological noise. On the other hand, the large high threshold MUs, that are either not active or not identified during the voluntary contractions will not be tracked in H-reflexes. It is therefore expected that identified MUs will never fully explain the energy of H-reflexes. The use of a spatial filter (e. g. Laplacian filter used in this study) decreases but does not fully resolve this problem.
Finally, H-reflex decomposition is likely facilitated by the dispersion of MU firing latencies caused by the differences in conduction velocities and/or length of afferent and efferent axons in the investigated muscle. In this respect, the proposed method might yield better results in distal than in proximal muscles. In this study, we limited our H-reflex decomposition to soleus muscle, which has been frequently studied in the past. H-reflex decompositions in other muscles are yet to be investigated.
In conclusion, we demonstrated that the identification of MU firings in electrically elicited H-reflexes is feasible. Although there are several limitations yet to be addressed, the methodology presented herein is opening new ways of in vivo human motor system investigations.