On the Prediction of Motor Unit Filter Changes in Blind Source Separation of High-Density Surface Electromyograms During Dynamic Muscle Contractions

We study the changes of Motor Unit (MU) filters in MU identification from high-density surface electromyograms recorded during isokinetic dynamic contractions of biceps brachii muscle. We demonstrate that these changes can be predicted for limited changes of the joint angle by linearly extrapolating previously recorded changes, allowing for the linear prediction-correction paradigm of MU filter updating. We then demonstrate the efficiency of this paradigm by implementing MU filter updating by the Kalman filter and integrating it into the previously published Convolution Kernel Compensation (CKC) MU identification method. When compared to the original CKC method and the previously published cyclostationary CKC method devoted to MU identification in repeated dynamic contractions, the Kalman based MU filter prediction yielded a superior precision of MU firing tracking in dynamic contractions. In the case of relatively fast biceps brachii contractions with full elbow flexion and extension in 2s, the Kalman based MU filter prediction tracked 21.3 ± 1.8 MUs with an average sensitivity of 95.6 ± 7.0% and precision of 96.5 ± 3.5%. In the same conditions, the original CKC method identified 7.1 ± 2.0 MUs with an average sensitivity of 62.7 ± 20.1% and precision of 98.1 ± 3.9%, whereas cyclostationary CKC tracked 18.9 ± 2.0 MUs with an average sensitivity of 91.8 ± 12.2% and precision of 94.7 ± 5.1%.


I. INTRODUCTION
Identification of neural codes that govern human movement is gaining relevance in our aging society. The number of studies clarifying the potential benefits and applications behind such identification is increasing rapidly with a particular focus on population coding [1]- [3]. Namely, the rapid development of systems for acquisition and processing of multichannel intramuscular [4], [5] or surface electromyographic (EMG) recordings [6] enabled the identification of the codes of several tens of simultaneously active motor nerves [7]- [9]. In this measuring setup, the skeletal muscles act as natural amplifiers of neural codes. Indeed, each motor nerve innervates The associate editor coordinating the review of this manuscript and approving it for publication was Ghulam Muhammad . several tens to hundreds of muscle fibers, forming with them the basic functional unit, the so-called Motor Unit (MU) [10], [11]. In healthy conditions, each electrical pulse conducted by the motor nerve triggers the electrical response in all the innervated muscle fibers, the so-called Motor Unit Action Potential (MUAP). Due to the relatively large number of innervated muscle fibers, MUAPs have an amplitude much higher than neuron action potentials, and can thus be, detected reliably at relatively large distances from their origin. Indeed, the amplitude of MUAPs is high enough to be detected by electrodes on the surface of the skin, above the investigated muscle [12].
The Central Nervous System (CNS) controls the exerted muscle force by activating/deactivating the MUs and by modulating their firing rates. MUAPs of different MUs are 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/ superimposed linearly in time and space, forming highly complex EMG patterns [13]. Namely, the EMG amplitude depends on the number of activated MUs and their firing rates, but also on the distance of the MUs from the EMG uptake electrodes [13]. The process of neural code estimation from EMG separates the impacts of MU location in the muscle tissue from their excitation patterns and, therefore, eases the interpretation of muscle excitation [7]- [9]. However, neural code estimation from the EMGs is not a trivial task, and computer-aided methods have been designed to estimate and invert the MUAP mixing process in EMG mathematically [7], [8], [14]- [16]. These techniques typically require multichannel EMG recordings and identify the MU filter, which, when applied to a surface EMG, yields firing patterns of individual MU. Several tens of simultaneously active MUs can be identified in this way. However, these methods have been largely limited to stationary EMG recordings in which the shape of the MUAPs does not change significantly [7], [8], [14]- [16]. In dynamic muscle contractions, as well as in the case of severe muscle fatigue, the MUAPs change in time, and this increases the complexity of the EMG decomposition to contributions of individual MUs [11]. Different strategies have been proposed to tackle this challenge [14], [17]. They all include some sort of MU tracking. For example, dynamic muscle contractions can be repeated in time, resulting in cyclostationary MUAP shapes. In [17], this property was exploited to update the MU filter estimates across different contraction repetitions, leading to successful tracking of the MU firing patterns in repeated dynamic conditions.
In this study, we extend the EMG mixing model tracking by the ability of MU filter prediction. We first demonstrate the ability of MU filter prediction by linear extrapolation of previously recorded MUAP changes. Afterwards, we propose the Kalman filter based MU filter updating in dynamic EMG conditions. We demonstrate the efficiency of this proposed approach on the dynamic contractions of the biceps brachii muscle with simulated and, thus, known MU firings that were convolved with experimentally recorded MUAPs. We apply the proposed decomposition scheme to simulated EMG signals with different levels of noise and different speeds of muscle contractions. We then test the proposed Kalman based MU filter updating scheme systemically against the internal parameters selection. Finally, we compare the decomposition method proposed herein with the previously published Convolution Kernel Compensation algorithm [8], [18], designed for isometric muscle contractions, and with the recently proposed cyclostationary CKC algorithm [17], designed for repeated dynamic contractions. The simulated conditions also allow for MU filter calculation out of known MUAP shapes in each time moment. Therefore, we also evaluate this MU filter updating out of known MUAP shapes, demonstrating that its performance is inferior to the proposed Kalman-based MU filter updating.

A. HDEMG EXPERIMENTS AND SIMULATIONS
In order to facilitate the objective evaluation of decomposition performance, we simulated the dynamic HDEMG signals by using previously recorded MUAPs of biceps brachii during slow dynamic contractions [17]. In short, five healthy males (age of 34.4 ± 5.4 years, height of 1.77 ± 0.05 m and weight of 77.2 ± 5.5 kg) participated in the experiment. Before participating, the subjects were provided with a detailed explanation of the study and gave their written informed consent. The study was conducted in accordance with the Declaration of Helsinki, and was approved by the local Ethics Committee.
The subjects had their dominant arm put into a custom made brace with a rotatable lever and performed slow 80 s long isokinetic dynamic contractions in the sagittal plane, starting at elbow fully extended and ending at elbow fully flexed. Using a pulley system, 1 kg weight was fixed to the lever, thus providing a constant torque of 3 Nm to the subject's arm. An array of 13×5 electrodes (OT Bioelettronica, Torino, Italy, interelectrode distance of 8 mm) was centered over the short head of the dominant biceps brachii muscle. Surface EMG signals were amplified (EMG-USB2 amplifier, OT Bioelettronica, Torino, Italy), filtered between 20 and 700 Hz, and sampled at 2048 Hz and 12-bit resolution.
MUAP shapes at different muscle shortening levels were identified offline. For this purpose, the consecutive 10 s long epochs of surface EMG signals were decomposed by the CKC method [17]. The epoch overlapping was set to 50 % and overlapped epoch sections were used to pair the identified MUs firing patterns on different epochs and to identify the entire MU firing pattern. Afterwards, MUAPs of 252 identified MUs were tracked through the dynamic contractions by Spike Triggered Averaging (STA) of 36 consecutive 10 s long EMG epochs that overlapped for 7.8 seconds [17]. This resulted in 36 MUAP shapes corresponding to 36 muscle shortening levels (level 1: muscle fully extended, level 36: muscle fully shortened).
The MU firing patterns were simulated by the model proposed in [19]. The model parameters were adapted to the properties of the biceps brachii muscle. In total, 200 MUs were assumed to be active in the detection volume of the surface electrodes. When recruited, MUs began firing at 8 Hz and increased their firing frequency linearly to 35 Hz at the maximum muscle excitation. The MU interspike interval variability was normally distributed with the coefficient of variation set to 20%. The MU recruitment order was distributed exponentially, with many lowthreshold MUs and progressively fewer high-threshold MUs. The last MU was recruited at 80% of maximum muscle excitation. Stationary 10% and 30% muscle excitation were simulated, resulting in 105 and 155 active MUs per contraction.
Different muscles were simulated by selecting the set and the order of MUs randomly from experimentally estimated MUAPs of 252 MUs. As explained in [17], the experimentally estimated MUAP library contains only MUs that were identifiable from surface EMGs. In order to compensate for this bias and also simulate the contributions from many unidentifiable MUs (so-called physiological noise), we followed the procedure described in [17]. In particular, we scaled the MUAPs by a factor β, where x is a random integer uniformly distributed on the interval [1,NoMUs], and NoMUs stands for the number of simulated MUs. A representative example of estimated MUAP shapes is depicted in Fig. 1.

B. HDEMG MODEL AND MU FILTERS
In dynamic muscle contractions, multi-channel surface EMG can be modeled by [17], [18]: where with δ denoting the unit-sample pulse and the k-th MUAP of the j-th MU appearing at time τ j (k). The mixing matrix H(n) = [H 1 (n) . . . H N (n)] comprises all the L samples long MUAPs as detected by the surface electrodes [17], [18], where for j = 1 . . . N with h ij (n, l) defining the l-th sample of the j-th MU's MUAP as detected by the i-th uptake electrode at the n-th EMG sample.

C. MU SPIKE IDENTIFICATION
In the CKC method, the inversion of the EMG mixing process (2) is embedded into the so-called MU filter f j that, when applied to EMG, yields the estimate of the MU spike train of the j-th MU [8], [18]: where E() stands for mathematical expectation and C y = E(y(n)y T (n)) represents the correlation matrix of extended EMG measurements [8]. The CKC method estimates the MU filter and the MU spike train iteratively directly out of the EMG signals recorded during isometric contractions:f With thet j (n) denoting the delta train (3), the crosscorrelation vector E(t j (n)y(n)) represents the STA of EMG measurements, and contains F samples long epochs of the jth MU's MUAP, as detected in all the EMG channels. In this study, we tested extension factor F = 20, as this value is close to the ones recommended [8], [9].

D. LINEARITY OF MUAP CHANGES
In dynamic muscle contractions, the EMG-detected MUAP shapes are changing with the muscle geometry, and so does the optimal MU filter f j . In moderate dynamic contractions, these changes are gradual and continuous, and if we could track or even predict them, even when only on F samples long MUAP epochs used in MU filters (6), we could forecast the changes of the optimal MU filter, and, therefore, improve the estimation of MU spike trains in dynamic conditions. In order to verify the hypothesis of piecewise linearity of MUAP changes, we used the linear extrapolation of the experimentally measured MUAP shape changes on biceps brachii shortening levels c, c + 1 and c + 2 to predict the MUAP changes on the shortening levels linearly from c + 3 onward. We then used the predicted MUAP changes at shortening levels c + 3, c + 4, c + 5. . . to form the MU filter f j and track the MU firing patterns across the simulated 20 s long muscle shortening (from muscle fully extended to muscle fully shortened). In each shortening level, we measured the sensitivity, precision and Pulse-to-Noise Ratio (PNR) [20] of VOLUME 9, 2021 MU spike trains reconstructed by Eq. (5): where TP stands for the number of true positive firings, FP for the number of false positive firings and FN for the number of false negative firings. Tolerance with respect to the reference value was set to 0.5 ms, meaning that identified MU firing was classified as TP when it was identified within ±0.5 ms from its reference firing.
The PNR metric was defined as [20]: where E(x|ˆt j (n) r ) and E(x|ˆt j (n)<r ) denote the mean across all time moments in which the jth MU is estimated to have or not to have fired, respectively [20]. This estimation is achieved by applying the threshold-based segmentation to the estimated MU spike traint j (n), with r denoting the threshold value [8].
The method was tested on 4 simulations of 10% muscle excitation, each having 105 different active MUs. In this way, we assessed both the assumption on the linearity of MUAP changes and its impact on MU tracking in dynamic conditions. The results of these tests are listed in Subsection III-A.

E. MUAP TRACKING BY KALMAN FILTER
Encouraged by the results of the piecewise linearity of MUAP changes in Subsection III-A, we designed a Kalman filter that predicts and tracks the changes of the MF × 1 MU filter f j across different muscle shortening levels: where f K ,j (k) is the state of the Kalman filter at the k-th MU firing, u j (k) is the control input, w(k) is the state noise with covariance matrix Q, and v(k) is measurement noise with covariance matrix R.
The state transition matrix A and measurement matrix B were both set to identity matrix I of size MF × MF, whereas control input u j (k) was set to The noise co-variance matrices were set empirically to Q = 0.01 · I and R = I. The MU filter f STA,j (n) = E(t j (n)y(n)) as estimated by STA of the EMG on the interval [n − , n] [8], [18] was used as an observation z K ,j in the Kalman filter Eq. (11), and the estimated state of the Kalman filter was used in place of f j (n) in Eq. (5) to identify the MU spike train on the [n − − , n + ] interval. The precision (7), sensitivity (8) and PNR (9) of identified MU spikes were calculated in each [n − − , n + ] interval.
In order to increase the number of identified MUs and test the ability of MU filter tracking in dynamic contractions thoroughly, the f STA,j (n) used for initialization of the Kalman filter state at the beginning of each decomposition run was estimated from the CKC-based decomposition of simulated 5 second long isometric contractions at 30% excitation at muscle fully extended. The control input was initialized to zero (u j (0) = 0). Different values of the parameter were tested in different test runs, with excitation level set to 30% and set to 205 samples (corresponding to the time length of 0.1 s), 410 samples (0.2 s), 1024 samples (0.5 s) and 2048 samples (1 s), respectively. The parameter was fixed to 512 samples, corresponding to 0.25 s. Three different values of parameter α (α = 0, 0.5 and 1) and two different speeds of muscle contractions were tested, with two and ten repeated muscle shortenings and extensions in 20 s long simulated EMG signals. In each simulation, three different levels of additive colored noise (bandwidth of 20-500 Hz) were simulated, with Signal-to-Noise Ratio (SNR) set to ∞ dB, 20 dB and 10 dB, respectively. Finally, ten different muscles were simulated, selecting the experimentally estimated MUAPs for each simulated muscle randomly. This resulted in 3 (α) ×2 (speed of muscle shortening) ×3 (SNR) ×10 (simulated muscles) = 180 simulated sets of EMG signals and 180 × 4 (selection of parameter) = 720 EMG decomposition runs of the Kalman-based CKC algorithm.

III. RESULTS
The results of linear prediction of MUAP changes are reported in Subsection III-A, whereas the Subsection III-B reports the results of MU tracking by the Kalman filter.

A. LINEAR MUAP PREDICTION
We have selected three different starting combinations of adjacent shortening levels for linear extrapolation of MUAP changes, starting from shortening levels 1-3, 13-15 and 24-26. As the results in Figs. 2, 3 and 4 show, we were able to predict MUAP shapes across different numbers of additional shortening levels, thus confirming the piecewise linearity of MUAP shape changes. The results of testing across different simulated muscles are presented in Table 1. On average, the linear prediction of the MUAP changes supported reliable MU firing tracking for additional 7.3 ± 3.0, 7.7 ± 3.6 and 7.4 ± 2.6 shortening levels, when starting from shortening levels 1-3, 13-15 and 24-26, respectively.

B. MU IDENTIFICATION BY KALMAN FILTER
The results of MU tracking by the Kalman MU filter were compared to the a) MU filter constructed out of known MUAP shapes (grand truth used in the EMG simulations, denoted as Original MUAP or Orig. MUAP in the sequel), b) MU filter f STA,j (n) = E(t j (n)y(n)), as estimated by the STA of the EMG on the interval [n − , n], and c) to the cyclostationary CKC (csCKC) MU identification proposed in [17]. The representative examples of estimated MU spike trains are depicted in Fig. 5.  The performance of the algorithm was first evaluated against the selection of and α parameters. The best performance of the Kalman based CKC algorithm was observed with set to 2048 samples and α set to 0. Therefore, we report the results further for these parameter selections.
We then measured the numbers of MUs that were identified with overall sensitivity of at least 0.8 as calculated on the whole length of simulated EMG signals, their PNR, sensitivity and precision for different speeds of muscle contraction. The numbers of identified MUs were accumulated across different simulated muscles. The values of PNR, sensitivity and precision metrics were accumulated over different [n − − , n + ] intervals, and over identified MUs. The Shapiro-Wilk test rejected the distribution normality of the results. Therefore, a Friedman non-parametric test with Bonferroni correction was used for pairwise comparison. Statistical significance was set at p < 0.05. The results are presented in Figs. 6 -8 and Table 2.

IV. DISCUSSION
We demonstrated that the blind source separation methods, such as CKC, utilize portions (epochs) of MUAPs to form the MU filter, which, when multiplied by extended EMG channels, yields the spike train (3) of a given MU. The MUAP shapes on these epochs, and, thus, the MU filter, change during the dynamic contractions, due to the changes in muscle geometry. In this study, we demonstrated that in biceps brachii muscles, these changes are piecewise linear, at least to the extent that allows for prediction of these changes by linear extrapolation. We were able to track the MU spikes across approx. 20 % of elbow joint angle range by linearly VOLUME 9, 2021  extrapolating the MUAP changes of three consecutive MUAP shapes, estimated from approx. 9 % of elbow joint angle range.
We then exploited these findings to design and test the prediction of MU filter changes by a Kalman filter. For this purpose, we first used the CKC method [8], [18] to estimate the ME spike train on a portion of the EMG signal and used the Kalman filter to predict the MU filter changes and correct them with MUAPs shape estimations from dynamic contractions. The MUAP shapes used for correction by the Kalman filter were estimated by using the newly identified MU firings in the STA of EMGs, thus tracking the MUAP changes across the whole investigated range of elbow angles. This MUAP update phase depends on several internal parameters, and we performed an extensive set of simulations to identify their optimal values in different dynamic conditions. First of all, we investigated the length of the EMG epoch used in STA. On the one hand, we wanted this epoch to be as long as possible in order to increase the quality of MUAP estimates. On the other hand, we wanted to keep it short enough to track fast MUAP changes during the dynamic contractions. Thus, the length of the EMG epoch in STA needs to be set as a compromise between these two factors. In our study, we investigated relatively slow muscle contractions (full elbow flexion in 5 s) as well as relatively fast muscle contractions (full elbow flexion in 1 s). Furthermore, we tested three different noise levels. For both speeds and all three noise levels tested, the EMG epoch length of 2048 samples was optimal. These results were expected for slow contractions, but are a bit surprising for fast muscle contractions. They indicate that the increased variability of MUAP shape as estimated by STA on shorter EMG epochs is more detrimental for MU tracking than the gradual MUAP changes during the dynamic contractions of the biceps brachii.
We then tested whether we could speed up the tracking of MUAP changes by including the control input in the Kalman filter described by Eq. (10), and setting it to a portion of the difference between the last two consecutive MUAP shape predictions. The relative contribution of this control input was tuned by parameter α, setting its values to 0, 0.25 and 0.5. While all three tested values allowed for successful MU tracking across the whole movement range, the highest number of accurately tracked MUs was achieved by setting the parameter α to 0, that is, by ignoring the control input in Eq. (10).
We then compared the performance of the Kalman based MU filter prediction to the STA of the same EMG epoch (the very input into the Kalman filter), to the ground truth MUAP shapes used in EMG simulations, and to the previously introduced cyclostationary CKC method that exploits the repeated contractions, to gather the MU spikes and estimate the MU filer across all the repetitions of the same joint angle [17]. The results demonstrated the superior performance of the Kalman based MU filter update proposed herein. The Kalman based CKC method yielded a significantly higher number of accurately tracked MUs (Table 1, Figs. 6 -8). Also, the sensitivity and precision of MU firing identification were among the highest. The PNR value of MU spike trains, identified by the Kalman based CKC method, was often significantly lower than in the case of the cyclostationary CKC method, but still well above (for SNR > 10 dB), or at least very close (SNR = 10 dB), to the recommended value of 30 dB [20]. The other two methods (the STA based MU filter update and original MUAP based MU filter update) yielded a significantly lower number of accurately tracked MUs at lower sensitivity. Precision was still high for those two methods but very likely at the cost of lower sensitivity. These results demonstrate the advantages of dynamic MUAP tracking in comparison to MUAP estimation by STA. On the other hand, the MU filter estimated from original MUAP shapes appeared to be specifically adapted to the current MUAP shape and was not able to identify the firings from nearby HDEMG epochs.
Our study has several limitations. First, all the MUs studied herein are from the biceps brachii muscle only and likely superficial, as they were identified by decomposing a surface EMG. Therefore, the conclusions of this study cannot be generalized to other muscles.
Second, we investigated relatively low contraction levels. It remains to be tested whether the presented results are also representative for high contraction levels, and, therefore, high threshold MUs. The same consideration also applies to MUAPs in different pathologies, as, in this study, we limited our investigations to healthy young males.
Also important, we used relatively short 5 seconds long isometric contractions to identify the MUs and then tracked the changes of their MUAPs by Kalman filter. This setup is always possible in experimental isokinetic conditions but requires an isometric calibration phase. The impact of the length of isometric EMG signals in the calibration phase was not assessed in our study.
Finally, we investigated the MUAP shape tracking on relatively short 20 samples long MUAP epochs, used in the CKC method. The CKC method centered these MUAP epochs automatically between the minimum and maximum MUAP peaks (Fig. 1). Thus, we cannot generalize the claims of linear MUAP changes to the entire MUAP shapes. For the same reasons, we cannot claim that the Kalman based MUAP prediction is also efficient in tracking of other physiological variables, for example, in MU conduction velocity estimation.

V. CONCLUSION
We demonstrated that, when observed on relatively short epochs, the MUAP changes in dynamic contractions of biceps brachii can be predicted by linear extrapolation. This allows for the exploitation of linear prediction-correction schemes in MU firing tracking across dynamic muscle shortenings. The results presented herein are promising and demonstrate the increased performance of the state-of-the-art MU tracking and neural codes investigation in dynamic conditions. The latter is currently lacking in many application fields, ranging from Neurophysiology, Neurology, Sport Sciences to Ergonomics and even the Gaming industry.