Adaptive HD-sEMG decomposition: towards robust real-time decoding of neural drive

Objective. Neural interfacing via decomposition of high-density surface electromyography (HD-sEMG) should be robust to signal non-stationarities incurred by changes in joint pose and contraction intensity. Approach. We present an adaptive real-time motor unit decoding algorithm and test it on HD-sEMG collected from the extensor carpi radialis brevis during isometric contractions over a range of wrist angles and contraction intensities. The performance of the algorithm was verified using high-confidence benchmark decompositions derived from concurrently recorded intramuscular electromyography. Main results. In trials where contraction conditions between the initialization and testing data differed, the adaptive decoding algorithm maintained significantly higher decoding accuracies when compared to static decoding methods. Significance. Using “gold standard” verification techniques, we demonstrate the limitations of filter re-use decoding methods and show the necessity of parameter adaptation to achieve robust neural decoding.


Introduction
Force generation in human skeletal muscles is governed by the activity of constituent motor units (MUs).Each MU is comprised of a single alpha motor neuron and the set of muscle fibers that it innervates, where a single axonal action potential initiates a tension-generating contractile twitch in the innervated fibers.The discharge pattern of a MU population thus encode the neural drive underlying gross muscular contraction [1,2].Historically, the precise activation times of individual MUs were only attainable via manual or semi-automatic spike sorting of electromyography (EMG) signals measured from indwelling electrodes [1,[3][4][5][6].More recently, convolutive blind source separation techniques have been developed to automatically extract MU spike trains from high-density surface electromyography (HD-sEMG) [7][8][9].Such methods yield detailed neural information in a non-invasive manner and are capable of extracting far more MUs compared to the spike sorting of intramuscular EMG (iEMG) [10].For these reasons, HD-sEMG decomposition has garnered considerable interest in studies on neurophysiology, motor control and neuromuscular disorders [10][11][12][13].In particular, MU decomposition offers practical advantages over established modes of human-machine interfacing (HMI) due to the access to higher neural information without the need of invasive procedures [14][15][16].
Traditionally, decomposition yields a set of separation vectors (MU filters) that distill HD-sEMG into underlying source activities.This process relies on repeated execution of iterative numerical methods over observations spanning substantial periods of time, typically 10 s or more [8,17].Hence, such batch decomposition algorithms are unsuitable for realtime deployment.Instead, reapplication of batchdecomposed MU filters to real-time measurements has been a commonly adopted approach [16,18,19].However, these techniques assume surface MU action potentials (sMUAPs) to remain consistent.In reality, factors such as fatigue, contraction intensity, and joint position alter the expression of sMUAPs on the skin surface [20][21][22][23].To tackle this challenge, decoding algorithms that adapt to new data have been developed [17,24].However, these methods have been tailored to specific conditions and are yet to be evaluated against the gold standard reference of iEMG-decomposed spike trains.
Here we propose a real-time MU decoding algorithm that updates the MU filter and signal preprocessing transforms as new action potentials of the observed MU emerge.The algorithm was evaluated on HD-sEMG recordings pertaining to isometric wrist extension contractions that vary across contraction intensities and joint angles.The accuracy of the algorithm was verified using reference spike trains manually decomposed from concurrently recorded fine-wire iEMG.

HD-sEMG decomposition
The decomposition techniques employed in this work are based on a convolutive mixture model for EMG generation: where z i (k) is the value of HD-sEMG channel i at time instant k. τ j (k − l) is the pulse train of MU j while a ij (l) encodes its respective action potential.L is therefore the maximum duration of impulse responses that is considered in the model and ε i (k) is additive noise inclusive of the activities of unextractable MUs.The algorithm for batch decomposition is described in detail in [8] though, here, a brief overview will be given for completeness.Firstly, the HD-sEMG signals are extended by appending their timedelayed versions as additional observations.This conditions the data for the FastICA algorithm [25], which normally decomposes instantaneous mixtures, to handle convolutive mixtures [26].Further preprocessing of the observations includes zero-phase component analysis (ZCA) sphering, which aids in the convergence of FastICA [26].The batch algorithm then extracts underlying source activities in a sequential manner, thereby estimating the firing intervals of MUs responsible for the generation of the observed HD-sEMG.Each source signal is extracted as: where z(k) is the extended observation vector and is the sphering matrix, calculated as the inverse square root of the covariance matrix of extended obsrvations, Finally, b is the spatiotemporal filter that extracts the MU contribution.The processes involved for simultaneously estimating b and ŝ(k) in a blind manner can be broken down into the Extraction and Refinement step.The Extraction step employs FastICA which iterates a fixed-point algorithm with an objective function optimizing the sparsity of ŝ(k).Orthogonalization of MU filters is performed at every iteration to ensure convergence to new sources.In the Refinement step, the MU filters and spike trains are iteratively updated to optimize the silhouette measure (SIL), a value which measures the accuracy of the separation [8,23].
As per [7], each iteration first involves peak detection on the estimated source signal.Spike classification is then performed where the kmeans++ algorithm is used to distinguish peaks as either spikes or noise, with cluster centroids c hi and c lo , respectively.Finally the MU filter is re-calculated as the cross-correlation between the sphered, extended observations and the current estimated spike train: where zΨ represents members in a set of extended observations corresponding to time instants of estimated spikes, Ψ = {z (k) : τ (k) = 1}.The Refinement step is thus repeated until the SIL value of the re-estimated source ceases to increase and sources with a final SIL value above a minimum acceptance threshold are deemed as viable MU pulse trains.

Online decomposition 2.2.1. Static decoding
So far, the most prominent approach to estimating MU activities in real-time is through the reuse of the MU filter, b, and pre-process transforms, Σ −1/2 zz and µ z, as presented by Barsakcioglu and Farina [27].These are initially obtained from training data and then continuously reapplied to new windows of extended data in the same manner as equation (2).Detected peaks in the estimated source signal are further sorted as either spikes or noise peaks.Rather than using the kmeans++ algorithm, this is simply determined by a threshold set at the midpoint between the spike and noise centroids, c hi and c lo , also retained from batch decomposition of the training data.To accommodate for deviations between the conditions of the training data and the new, unseen data, this decision boundary may be altered by a relaxation factor, 0 ⩽ α ⩽ 1: for the adaptive decoding algorithm can therefore be written as: where v is now the MU filter and Σ −1 zz is the inverse of the observation covariance matrix.
With each new data window, temporary transforms are first derived from the updated statistics of extended observations: where zwin (k) is the kth sample in the new window of extended data and 0 ⩽ λ ⩽ 1 is the forgetting factor that controls the influence of new data.An initial estimation of the source signal is then calculated: where Σ * −1 zz is the inverse of the temporary covariance matrix obtained from equation (7).
Peak detection is then conducted on the estimated source signal of the tracked MU, ŝwin (k).To identify potential new spike instances for learning, the sparseness property of MU firing patterns can be leveraged.With a-priori knowledge regarding the short timespan that the data window corresponds to, any strong responses to the MU filter, meaning potential spikes, will appear as outliers in the distribution of rectified peak amplitudes.Hence, rectified peak amplitudes with z-scores above a certain threshold will have their corresponding extended observation vectors added to set Ψ mem .Functionally, Ψ mem is implemented as a first-in-first-out storage buffer of constant size, initialized from extended observations corresponding to spike events in the training data.As candidate spikes are detected from new data windows and new observation vectors are added to Ψ mem , past observations are discarded.With each update of Ψ mem , the MU filter is recalculated using equation ( 9): Equation ( 9), similar to equation (3), updates the MU filter via cross-correlation between new extended observations and the |Ψ mem | most recently estimated spikes.
The spike and noise centroids, which are used for online spike detection outside of the adaptation algorithm, are subsequently updated.The spike centroid is recalculated via equation (10) which corresponds to the squared average amplitude of peaks extracted from the observations stored in Ψ mem .The noise centroid is then updated as a λ-weighted merging of the past c lo and the average of noise peak amplitudes detected from the new source signal window: where η win is the set of observations corresponding to noise peaks detected inŝ win (i.e.peak observations not added to Ψ mem ).Algorithm 1 summarizes this entire process for real-time updating of decomposition parameters.Prior to implementation, there are three static parameters that need to be defined: the threshold z-score value for accepting new observations into Ψ mem , the rate of forgetting, λ, and the cardinality of Ψ mem .As in past studies, such parameters were selected empirically [17,24].For the results obtained in this work, the corresponding values used were 3.3, 0.985 and 110, respectively, based on initial testing.Add new spike observations to Ψmem while discarding an equal number of oldest observations.8.
Build η win from observations corresponding to noise peaks extracted from new data.9.
Accept updated inverse covariance matrix: Accept updated statistics:

Experimental setup
Five able-bodied subjects were recruited for the experiment, four male, one female, ages 29-34, all right-handed.The study was approved by the local ethical board of Aalto University (approval number D/505/03.04/2022).Prior to the experiments, all subjects gave their written informed consent in accordance with the Declaration of Helsinki.Subjects were seated for the duration of the experiment with their dominant upper-limb placed in a specialized tabletop rig designed to constrain the wrist joint at various angles of extension (figure 1).Forces generated by isometric contractions pertaining to wrist extension were measured with a load cell (TAS606, HT Sensor Technology, China) at a sampling rate of 100 Hz.Prior to the insertion of fine-wire electrodes, the subject's maximum voluntary contraction (MVC) forces were measured at wrist joint angles corresponding to 0%, 12.5% and 25% of their maximal extension, with 0% relating to a neutral wrist position.MVC was calculated as the averaged maximal force from three MVC contractions of 1.5 ms long with short breaks in between each contraction to prevent fatigue.
Stainless steel/silver (SS/Ag) wires with polytetrafluoroethylene insulation (Spes Medica s.r.l, Italy) were used as intramuscular electrodes.The wires had a diameter of 0.11 mm with the final 3-5 mm of the recording tips stripped of the insulating material.Three insertion points were targeted, centered at the bulk of the extensor carpi radialis brevis (ECRB) and aligned down the muscle axis at approximately 4 mm intervals.Location of the ECRB was guided by [28] and palpation during wrist extension and radial deviation movements.The fine-wires were inserted as pairs (bipolar configuration) using 25G cannulae to a depth targeting MUs proximal to the skin surface.Signal inspection was conducted after the insertion of each electrode pair.If the signal was invalid (short-circuited, excessive noise, low selectivity or no viable units detected) and could not be remedied by light manipulation of the fine-wires, the wires were removed and another insertion of new electrodes was made slightly lateral to the original insertion point.The maximum overall number of insertion attempts was bounded to five for the sake of subject comfort.The experiment only proceeded so long as at least one valid iEMG channel was attained.The bipolar iEMG signals were preamplified by an adapter (ADx5JN, OT Bioelettronica, Italy) with a gain of 5, and acquired by a bioamplifier (Quattrocento, OT Bioelettronica, Italy) with a fixed gain of 150 at 10 240 Hz with 10-4400 Hz hardware bandpass filtering.Subsequent processing of iEMG signals included high-pass filtering with a 250 Hz cutoff to lower baseline noise and to produce sharper action potentials [5,29].
Placement of the overlaying HD-sEMG matrix was conducted approximately 8 minutes after the final fine-wire insertion.This allowed for sufficient coagulation, minimizing the leakage of blood or plasma to the surface recording site.A 64-channel rectangular electrode matrix (GR08MM1305, OT Bioelettronica, Italy) with 8 mm inter-electrode distance was placed on top of the ECRB, centered above the iEMG insertion sites (figure 1).Two reference electrodes (Neuroline 720, Ambu A/S, Denmark), one for the pre-amplifier and one for the bioamplifier, were placed at the medial epicondyle and olecranon process.The HD-sEMG signals were buffered by a pre-amplifier (AD64F, OT Bioelettronica, Italy) prior to being acquired by the same benchtop amplifier used for iEMG at 150 gain, 10 240 Hz with 10-4400 Hz hardware bandpass filtering.Pre-processing of the HD-sEMG signals for automatic decomposition included downsampling to 2048 Hz and bandpass filtering with 10-900 Hz cut-offs.
Prior to the commencement of recordings, subjects were asked to perform slow dynamic wrist extensions, up to 25% of maximum range of movement, to allow the settling-in of fine-wire electrodes and HD-sEMG matrix.The recording and cueing of contractions were facilitated by a custom Matlab R2021b (MathWorks Inc.USA) framework.All subject cues, along with the real-time force feedback, were displayed on a computer screen.

Experimental protocol
Isometric wrist extension contractions with trapezoidal force profiles (5 s ramp, 20 s plateau) were recorded at different joint angles and different force levels.To ensure iEMG decomposability, which relies on low to moderate signal complexity [4,5], contraction intensities were kept at low levels.For Figure 2. iMUAPs and sMUAPs extracted from different contraction conditions.Up to 3 unique iMUAP shapes are utilized for manual matching of MUs extracted from different contraction conditions.While each individual iMUAP profile is still susceptible to change across the angle and force conditions, causing potential matching ambiguities, the presence of multiple time-locked and distinct profiles facilitate matching of MUs that may share similar iMUAP profiles in one particular channel.Variation in the sMUAP profiles across contraction conditions is also observed, resulting in sub-optimal extraction of source activities when using static decoding algorithms.(a) iMUAPs and sMUAPs of MU B1 extracted from all 5 contraction conditions that it was detected in.From darkest to lightest plot lines, the displayed MUAPs correspond to angle/force combinations of 0%/5%, 0%/10%, 25%/5%, 25%/10% and 25%/15%, respectively.(b) iMUAPs and sMUAPs of MU B2 obtained from the same contraction conditions as displayed in (a).(c) iMUAPs and sMUAPs of MU C1 extracted from 3 different contraction conditions.From darkest to lightest plot lines, the displayed MUAPs correspond to angle/force combinations of 0%/5%, 12.5%/7.5% and 25%/10%, respectively.(d) iMUAPs and sMUAPs of MU C2 obtained from the same contraction conditions as displayed in (c).
subjects A and B, contractions were recorded at force levels of 5%, 10% and 15% MVC at wrist joint angles of 0% and 25% maximal extension.For subjects C, D and E, contractions of 5%, 7.5% and 10% MVC were recorded at 0%, 12.5% and 25% maximal wrist extension.Recordings progressed from 0% to 25% extension while the order of force levels recorded was randomized.Three repetitions were recorded for each contraction condition.

Obtaining iEMG decomposition benchmarks 2.5.1. Extraction of MU activity concurrent in iEMG and sEMG
To identify MUs present in both surface and intramuscular signals, a two-stage semi-automatic technique was employed.For each repetition, a set of MUs and their respective spike trains are first extracted from HD-sEMG via the batch decomposition method described in section 2.1.The resultant spike intervals were then used to trigger action potentials in the iEMG signals.Here, MUs whose activities are present in both the concurrently recorded HD-sEMG and iEMG signals will trigger distinct intramuscular motor unit action potentials (iMUAPs).
Typically, these are mono and polyphasic waveforms with peaks well above the baseline noise [4].On the other hand, MUs that were only extractable via HD-sEMG decomposition will trigger flat iMUAPs (peak-to-peak amplitudes <2 µV).In this way, units present in both surface and intramuscular recordings are identified.The spike-trains of such MUs were then imported to EMGlab [29], a Matlab-based spike annotation software, for manual correction by an experienced operator such that a high-confidence benchmark is obtained.

Tracking MUs across contraction conditions
MUs were matched by the same experienced operator through visual comparison of their multi-channel iMUAPs.As each iEMG channel consisted of a bipolar measurement, activity from a single source manifests as action potentials that vary greatly in profile across channels but are, albeit, time-locked.Thus, a single MU may be characterized by up to three distinct action potentials triggered by the same spike instances.Examples are shown in figure 2 where such iMUAP profiles may be used to manually match MUs across contraction conditions.

Pseudo-online testing
In the pseudo-online tests, multiple trials were conducted to gauge the robustness of the proposed adaptive MU decoding algorithm across different contraction conditions.In each trial, the decoding algorithm was initialized from one repetition and then applied to extract MU activity in another.Here, data was fed in windows of 200 ms and in time increments of 100 ms, thereby simulating real-time deployment.For comparative purposes, the static decoding technique (section 2.2.1) was also tested using different spike threshold relaxation values, from α = 0 to α = 0.5 in increments of 0.1.
Since the decoding algorithms were to be compared in scenarios where the conditions of the training data differed from those of the test data, only MUs with high-confidence iEMG-decomposed benchmarks (obtained by methods described in sections 2.5.1 and 2.5.2) in at least two force levels for at least two angle conditions were selected for this analysis.Table 1 lists the MUs selected for this testing along with the contraction conditions in which they were detected in.For each eligible MU, all pair-wise combinations of training and testing repetitions were analyzed.
For each trial, the estimated spike train was compared to the iEMG-decomposed spike train using the Rate-of-Agreement (RoA) metric: where O is the number of spikes that were only detected by the online decomposition algorithm and I is the number of firing instances exclusive to the iEMG decomposition while C is the number of spikes that were identified in both estimations of MU activity.
In addition, two metrics that are analogous to False Negative Rate (FNR) and False Discovery Rate (FDR), when considering the iEMG-decomposed spike train as ground truth, were calculated for each trial:

Statistical analysis
To detect statistically significant differences between the decoder performances, repeated-measures analysis of variance (RM-ANOVA) was conducted on the RoA values obtained from pseudo-online testing.The normality of the results was verified with Shapiro-Wilks's testing while the assumption of sphericity was tested with Mauchly's test.In cases where the sphericity assumption was not satisfied, Greenhouse-Geisser correction was applied to the RM-ANOVA.
If the choice of decoder was found to have a significant effect, post-hoc pair-wise comparisons (Tukey-Kramer) were conducted.In all analyses, significance levels of 0.05 were used.
In addition to analyzing the full set of results, three auxiliary analyses, 'Intra-condition' , 'Interangle' and 'Inter-force' , were conducted on different subsets of the data.In the Intra-condition analysis, only the trials where the training and testing data had identical angle/force conditions were considered.In the Inter-angle analysis, trials where the training and testing data were recorded from identical force conditions, but had different angle conditions, were considered.Finally, the Inter-force analysis focused on trials where the training and test recordings consisted of contractions with identical angle but different force conditions.

Results
Eighteen MUs (table 1) were found to satisfy the inclusion criterion for this study.This yielded 6073 In the Inter-angle analysis (N= 1192), only trials with differing joint angle conditions, but identical force levels, between the training and test repetitions are included.For the Inter-force analysis (N = 1394), the included trials feature identical joint angle levels between the training and test repetitions but differ in force level.RM-ANOVA was conducted on the RoA results with the decoding algorithm selected as the independent variable.Statistically significant differences between the online decoders was found in all the analyses.Pairwise comparisons between decoders with statically significant differences are indicated by ' * ' in the grids above the bar-charts.The proposed adaptive algorithm is shown to achieve superior robustness as high RoA is maintained when tested contraction conditions differ from those used for decoder initialization.(Middle and bottom) FNR and FDR for each the decoder is also shown.Adaptive decoding is shown to negate the trade-off in increasing α, which lowers the occurrence of false negatives at the expense of higher false positive rates.trials in total, of which, 831, 1192, and 1394 were included in the Intra-condition, Inter-angle, and Inter-force analyses, respectively.Figure 3 shows the results from pseudo-online testing of the static and adaptive decoders in terms of RoA with iEMGreferenced benchmark decompositions.Statistically significant differences between decoder performances were detected in all analyses (F(  2.00, 438.80) = 78.94,p < 0.001 for Global, Intracondition, Inter-angle, and Inter-force, respectively).In post-hoc comparisons, the proposed adaptive decoding algorithm significantly outperformed static decoding for all tested α values (0-0.5) in the Global, Inter-angle, and Inter-force analyses (p < 0.001 for all comparisons).Overall, α = 0.3 gave the best static decoder performance with an RoA of 77.1% ± 25.2% in the Global analysis.Still, this was exceeded by the adaptive decoder by 6.7% ± 0.2%.Similarly, in the Inter-angle and Inter-force analyses, the best static decoding performances were 75.7% ± 25.0% (α = 0.3) and 80.8% ± 22.4% (α = 0.2) RoA, respectively.The adaptive decoder also outperformed these by 8.0% ± 0.4% and 5.1% ± 0.4%, respectively.In the Intra-condition analysis, static decoding is shown to still perform well with α = 0.1 and 0.2 yielding the highest average RoAs of 94.1% ± 7.5% and 93.9% ± 7.5%, respectively.Adaptive decoding marginally underperformed these by −0.7% ± 0.2% and −0.6% ± 0.2%, respectively.
Figure 3 also shows decoding performances in terms of FNR and FDR.Here, the effect of increasing α is clearly shown.By relaxing the spike amplitude threshold, fewer spikes are missed (lower FNR) but in turn, more noise peaks are misclassified as spikes (higher FDR).In contrast, the adaptive decoding algorithm maintains low rates of either misclassification types.Compared to static decoding with α = 0.3, which yielded 14.2% ± 13.7% FNR and 16.6% ± 21.8% FDR in the Global analysis, adaptive decoding achieved lower misclassifications by −0.5% ± 0.2% and −8.1% ± 0.2%, respectively.The adaptive decoding algorithm therefore resolves this trade-off between FNR and FDR. Figure 4 shows how this is achieved by comparing the source activities extracted  Cells close to the heatmap diagonals represent results used for the 'Intra-condition' analysis as the decoders are tested on data pertaining to the same contraction conditions used for initialization.In such cases, static decoding with no relaxation of the spike detection threshold (α = 0, left) maintains high RoA with iEMG-referenced benchmarks.However, when tracking MU activity in contraction conditions that differ from the training data, RoA decreases.In some cases, this can be remedied by increasing α (middle), though the majority of decoding accuracies remain poor.The proposed adaptive decoder (right) therefore offers the best robustness with the majority of trials yielding RoA above 90%.However, the algorithm may not always compensate for some large changes in contraction conditions as shown by the few trials with low RoA (<50%) results.For instance, the algorithm can fail to converge to an appropriate filter when presented with contractions of wrist angle/force level combinations of 25%/10% when using decoding parameters initialized from contractions of 0%/5%.via online decoding with static and adapting parameters.By updating the MU filtering parameters as new data is received, a clear separation of spike peaks and noise peaks in the extracted source is maintained.
The RoA values of all trials regarding a single MU are shown in figure 5. Results obtained via static decoding with α = 0 and 0.3, which yielded the best overall static decoding performance, and the proposed adaptive decoding algorithm are shown.The adaptive algorithm yielded similar decoding accuracies as static decoders in the trials that fall under the Intra-condition analysis (cells proximal to the diagonals of the heatmaps).However, in the majority of trials where the training condition does not match the test condition, adaptive decoding achieved a much higher RoA with the iEMG-decomposed benchmarks.Still, there remain cases where adaptation is unable to compensate for the large changes to the sMUAPs.

Discussions
We have proposed an adaptive algorithm for decoding MU activity from HD-sEMG that continuously updates its internal parameters in real-time, as new measurements are acquired.Using experimental data, we demonstrate the performance of our proposed algorithm in tracking MU activities across isometric contraction conditions that vary in joint angle and intensity.In comparison to the static, nonadaptive decoding algorithm, adaptive decoding was shown to be more robust to such changes.This was verified against benchmark spike trains manually decomposed from iEMG signals that were recorded concurrently with the HD-sEMG.In terms of RoA between decoder estimations and the iEMGreferenced benchmarks, the adaptive decoder significantly outperformed static decoders across all tested spike threshold relaxation values (figure 3).Even when the test trials only differed from the training trial by one factor (Inter-angle or Inter-force), adaptive decoding was shown to be beneficial.Nonetheless, static decoding was still effective in estimating MU activity from contractions that are similar to the training data (Intra-condition).
As contraction intensity and joint angle change, so do the sMUAP profiles (figure 2).This renders sphering transforms and MU filters derived from different contraction conditions to be sub-optimal for accurate source estimation [30,31], resulting in missed spikes.In [23], local batch optimization of MU filters allowed for the accurate decomposition of MU activity during dynamic contractions.However, this required prior knowledge regarding the periodicity of the dynamic contractions.Here, we demonstrate that relaxation of spike acceptance thresholds can help compensate for changes to sMUAPs but this also causes an increase in false spike identifications, as evidenced by the inverse relationship between the FNR and FDR results obtained using static decoding (figure 3).Conversely, the adaptive decoder maintains low rates of both false negative and false positive errors.This is achieved by adaptation of pre-process transforms and the MU filter as new spikes are estimated, which helps maintain a distinct separation between spike and noise peaks (figure 4).
Previous studies on adaptive, real-time decoding algorithms include [17] and [24], both of which are based on the convolutional kernel compensation algorithm [7].In [17], an adaptive decoding algorithm was tested on a set of isometric contractions, ranging from 5% to 20% MVC, that were recorded from the tibialis anterior of eight subjects.Compared to spike trains extracted via batch decomposition, the real-time algorithm achieved an average RoA of 83%.In this work, we have achieved comparable accuracies (figure 3) in the more challenging scenario of decoding across contraction conditions.In [24], a separate algorithm was applied to decode simulated and experimental dynamic contractions.While dynamic contractions better represent the user input of HMI applications, only pulse-to-noise ratio was used to gauge decoder accuracy.In our study, the proposed algorithm has been directly verified using benchmark spike trains decomposed from intramuscular signals which remain the gold standard in the field [32,33].
Beyond algorithmic adaptation, in [34], the online decomposition algorithm was extended by incorporating a self-administered enhancement process utilizing the FitzHugh-Nagumo resonance model.This extension aimed to enhance MU source signals.When decomposing synthetic HD-sEMG signals where MU recruitment, sMUAP amplitude, and additive noise were varied, significant improvement over the baseline decomposition algorithm was achieved (88.70% ± 4.17% vs. 92.43%± 2.79%).Hence, employing physiologically-inspired models for further signal processing, in conjunction with adaptive decomposition, may yield even greater decomposition robustness.

Limitations and future work
Currently, our proposed algorithm has only been tested in a pseudo-online manner as verification of decoding accuracy against manually decomposed benchmarks necessitates offline procedures.Nonetheless, the algorithm is appropriate for realtime deployment.In this study, data windows were advanced in time increments of 100 ms while the average execution time of parameter adaptation, along with spike estimation, was 57.1 ± 14 ms on a desktop computer (Intel Xeon W-2133, 3.6 GHz, 36 GB RAM, Microsoft Windows 10, 64 bit).The main computational cost to the algorithm lies in the computation of the inverse covariance matrix.Despite this, prior testing has shown that omission of this step is detrimental to the overall effectiveness of the adaptive decoding algorithm.One way to reduce computational demand is to reduce the extension factor.In this work, we have used an extension factor of 16 to align with established works [8,35].However, past investigations suggest that lower extension factors can be employed and still retain user-intention estimation performances in HMI applications [36].While this work focuses on verifying the accuracy of the proposed adaptive decoding algorithm against iEMGreferenced benchmarks, deployment in real-time interfacing applications is left for future works.
Manual and semi-automatic decomposition of iEMG by experienced operators has been extensively verified from decades of research [4,37].For this reason, comparison with the decomposition of concurrently recorded iEMG, sometimes referred to as a variant of the two-source method [38][39][40], remains the most convincing means of experimental validation with regards to the decomposition of surface signals [32,38,41,42].However, low signal complexities are required to ensure the decomposability of iEMG [4,5].Hence, force levels were kept to 15% MVC and below in this study.Indeed, past works reliant on iEMG for verification have been similarly constrained to low level contractions [17,38,41,42].While we have directly verified that our proposed adaptive algorithm can compensate for sMUAP changes at lower force levels, further verification at higher force levels will have to depend on indirect measures of decomposition accuracy or simulations [24,31,34].For instance, the application of static filter re-use in contractions of up to 70% MVC have been previously investigated via such methods [31].
The issue of contraction intensity is closely linked to the nature of action potential superposition.The batch decomposition algorithm [8], which our work extends upon, is based on blind source separation which inverts the mixing process modelled by equation ( 1) and accounts for action potentials from different sources overlapping in time and space.The calculation of spatiotemporal filters capable of extracting MU activities despite such overlaps is aided by: 1) the high number of concurrent observations afforded by the high-density electrode grid, 2) statistical power afforded by longer recordings.The former aspect aids in discrimination of sources with spatially overlapped action potentials [43] while the latter permits computation of a sphering transformation which decorrelates the extended observations.This spatiotemporal decorrelation localizes sMUAPs in space and time, aiding in their separation.However, at higher contraction intensities, signal complexity from sMUAP superposition can increase significantly due to heightened MU recruitment and rate coding.Moreover, distinguishing sMUAPs between unique sources becomes more difficult due to the low-pass filtering effects of tissue volume conduction [43].This limitation is common to many blind source separation-based decomposition algorithms [44] though it may be overcome with signal acquisition methods that offer greater spatial selectivity [43].
The adaptive decoding algorithm may not always compensate for large, sudden changes in sMUAP profiles.As show in figure 5, when the difference between the training and test contraction conditions is significant, the adaptive algorithm may fail to converge to the correct filtering parameters.Here, the decoding algorithm is presented with an abrupt change from one isometric contraction to another, whereas, in practice, such changes occur in a continuous manner.Hence, future works will also focus on the application of adaptive decoding over experimental data pertaining to dynamic contractions.

Conclusions
In conclusion, we have developed an adaptive MU decoding algorithm that adapts to new data in real-time.Using high-confidence in-vivo-referenced benchmarks, the proposed algorithm was demonstrated to be more accurate in decoding MU activities across varying states of isometric contractions.This work therefore paves the way towards robust, realtime non-invasive neural interfacing.

Figure 1 .
Figure 1.Experimental setup.(a) Three fine-wire electrode pairs inserted into a subject's ECRB.(b) A 64-channel high-density surface electrode grid placed above the iEMG insertion sites shown in (a).(c) Subject with iEMG and HD-sEMG electrodes attached to their dominant arm which has been placed inside the force measurement rig.Task cues are shown on the computer screen in front of the subject.

Figure 3 .
Figure 3. (Top) Performance of static and adaptive MU decoding methods (SD and AD, respectively) in terms of RoA with iEMG-referenced benchmark decompositions.Error bars indicate standard deviation.The Global analysis encompasses all conducted trials (N= 6073).The Intra-condition analysis (N= 831) considers only trials in which contraction conditions of the training repetition used for decoder initialization matched those of the test repetition.In the Inter-angle analysis (N= 1192), only trials with differing joint angle conditions, but identical force levels, between the training and test repetitions are included.For the Inter-force analysis (N = 1394), the included trials feature identical joint angle levels between the training and test repetitions but differ in force level.RM-ANOVA was conducted on the RoA results with the decoding algorithm selected as the independent variable.Statistically significant differences between the online decoders was found in all the analyses.Pairwise comparisons between decoders with statically significant differences are indicated by ' * ' in the grids above the bar-charts.The proposed adaptive algorithm is shown to achieve superior robustness as high RoA is maintained when tested contraction conditions differ from those used for decoder initialization.(Middle and bottom) FNR and FDR for each the decoder is also shown.Adaptive decoding is shown to negate the trade-off in increasing α, which lowers the occurrence of false negatives at the expense of higher false positive rates.

Figure 4 .
Figure 4.Estimated source signals and spike trains of MU E1 using static (top row) and adaptive (bottom row) decoders.Only extraction of the first 4 s out of the full 12 s recordings are shown.Estimated spikes that agree with the corresponding iEMG-referenced decomposition are indicated by green circles.Spike estimates that disagree are indicated by red x's while missed spikes are indicated by purple crosses.For static decoders, spike estimations using α = 0 and 0.3 are shown.Spike amplitude thresholds are plotted as black dotted lines and dash-dot lines for the relaxed threshold (α = 0.3).The plots on the left side show the application of decoding parameters initialized from a wrist angle/force level combination of 0%/7% on a contraction of 25%/10%.The right side plots show results with swapped initialization and test data.With static decoding, the extracted signals are noisy and result in numerous misclassifications.The amount of missed spikes can be reduced by relaxing the spike amplitude threshold but this results in higher occurrences of misidentified spikes.With adaptive decoding, continuous updating of decomposition parameters maintain a clear separation of spike and noise peaks which result in higher decoding accuracies.

Figure 5 .
Figure 5. RoA results from tracking MU E1 using static and adaptive decoders.All possible pairwise combinations of training and testing repetition are shown.Cells close to the heatmap diagonals represent results used for the 'Intra-condition' analysis as the decoders are tested on data pertaining to the same contraction conditions used for initialization.In such cases, static decoding with no relaxation of the spike detection threshold (α = 0, left) maintains high RoA with iEMG-referenced benchmarks.However, when tracking MU activity in contraction conditions that differ from the training data, RoA decreases.In some cases, this can be remedied by increasing α (middle), though the majority of decoding accuracies remain poor.The proposed adaptive decoder (right) therefore offers the best robustness with the majority of trials yielding RoA above 90%.However, the algorithm may not always compensate for some large changes in contraction conditions as shown by the few trials with low RoA (<50%) results.For instance, the algorithm can fail to converge to an appropriate filter when presented with contractions of wrist angle/force level combinations of 25%/10% when using decoding parameters initialized from contractions of 0%/5%.

Table 1 .
Catalog of MUs and trial conditions used for this study.18 MUs were found to be identifiable, both in HD-sEMG and iEMG, in at least two levels of joint angle and force.Hence, they satisfied the inclusion criterion for this study.MUs sharing the same letter in their designation were extracted from the same subject and set of recordings.
o = MU concurrently detected in HD-sEMG and iEMG, x = MU not concurrently detected, -= trial not recorded for subject.