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

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–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–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–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–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 real-time deployment. Instead, reapplication of batch-decomposed 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–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.

2.1. HD-sEMG decomposition

The decomposition techniques employed in this work are based on a convolutive mixture model for EMG generation:

Equation (1)

where $z_(k)$ is the value of HD-sEMG channel i at time instant k. $\tau_(k-l)$ is the pulse train of MU j while $a_(l)$ encodes its respective action potential. L is therefore the maximum duration of impulse responses that is considered in the model and $\varepsilon_(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 time-delayed 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:

Equation (2)

where $\tilde}(k)$ is the extended observation vector and $\boldsymbol_}} = \mathbb \left[ \tilde}(k) \right]$ is the vector of subtractive means for centering $\tilde}(k)$. $\boldsymbol_}\tilde}}^$ is the sphering matrix, calculated as the inverse square root of the covariance matrix of extended obsrvations, $\boldsymbol_}\tilde}} = \mathbb \left[ \left(\tilde} \left(k \right) - \boldsymbol_}} \right) \left(\tilde} \left(k \right) - \boldsymbol_}} \right)^ \right]$. Finally, b is the spatiotemporal filter that extracts the MU contribution.

The processes involved for simultaneously estimating b and $\hat(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 $\hat(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 chi and clo , respectively. Finally the MU filter is re-calculated as the cross-correlation between the sphered, extended observations and the current estimated spike train:

Equation (3)

where $\tilde}_}$ represents members in a set of extended observations corresponding to time instants of estimated spikes, $ = \}\left( k \right): \hat (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.

2.2. Online decomposition2.2.1. Static decoding

So far, the most prominent approach to estimating MU activities in real-time is through the re-use of the MU filter, b, and pre-process transforms, $\boldsymbol_}\tilde}}^$ and $\boldsymbol_}}$, 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_\mathrm$ and $c_\mathrm$, 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 \unicode \alpha \unicode 1$:

Equation (4)2.2.2. Adaptive MU decoding

While relaxation of the decision boundary may decrease the likelihood of missed spikes, it can also lead to an increase in false positives. To address this, we propose adaptation of the MU filter and preprocess transforms as potential spike events are detected. Here, adaptation of decoding parameters occurs in parallel to the online decoding process. The updating of the MU filter is performed in a manner similar to the Refinement step of the batch algorithm. Since this does not rely on FastICA, re-computation of ZCA sphering transforms is unnecessary. Rather, only the inverse of the sample covariance matrix is needed. The equivalent of source extraction (equation (2)) for the adaptive decoding algorithm can therefore be written as:

Equation (5)

where v is now the MU filter and $\boldsymbol_}\tilde}}^$ 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:

Equation (6)Equation (7)

where $\tilde}_} \left( k \right)$ is the $k}$ sample in the new window of extended data and $0 \unicode \lambda \unicode 1$ is the forgetting factor that controls the influence of new data. An initial estimation of the source signal is then calculated:

Equation (8)

where $\boldsymbol_}\tilde}}^$ 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, $\hat_\mathrm(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 time-span 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 $_}$. Functionally, $_}$ 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 $_}$, past observations are discarded. With each update of $_}$, the MU filter is recalculated using equation (9):

Equation (9)

Equation (9), similar to equation (3), updates the MU filter via cross-correlation between new extended observations and the $|_}|$ 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 $_}$. The noise centroid is then updated as a λ-weighted merging of the past $c_\mathrm$ and the average of noise peak amplitudes detected from the new source signal window:

Equation (10)Equation (11)

where $\mathrm_}$ is the set of observations corresponding to noise peaks detected in $\hat_}$ (i.e. peak observations not added to $_}$).

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 $_\mathrm$, the rate of forgetting, λ, and the cardinality of $_\mathrm$. 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.

Algorithm 1. Adaptation of online MU decoding parameters1. Calculate $\boldsymbol^*$ with equation (6).2. Calculate $\boldsymbol_}\tilde}}^$ with equation (7).3. Calculate $\hat_\mathrm(k)$ with equation (8).4. Perform peak detection on $\hat_\mathrm(k)$.5. Sort peak as spikes or noise based on z-scores of their rectified amplitudes.6. if Spike peaks detected then 7.  Add new spike observations to $_\mathrm$ while discarding an equal number of oldest observations.8.  Build $\mathrm_}$ from observations corresponding to noise peaks extracted from new data.9.  Calculate the new MU filter, v, with equation (9).10.  Calculate new spike and noise centroids with equations (10) and (11).11.  Accept updated inverse covariance matrix: $\boldsymbol_}\tilde}}^ \leftarrow \boldsymbol_}\tilde}}^$.12.  Accept updated statistics: $ \boldsymbol_}} \leftarrow \boldsymbol_}}^*$, $\boldsymbol_}\tilde}} \leftarrow \boldsymbol_}\tilde}}^$ 13. end if 2.3. 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.

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.

Download figure:

Standard image High-resolution image

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 pre-amplified 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 cut-off 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.

2.4. 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 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.

2.5. Obtaining iEMG decomposition benchmarks2.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.

2.5.2. 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.

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).

Download figure:

Standard image High-resolution image 2.6. 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.

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.

Angle (%)Force (%)MUA1A2B1B2B3C1C2C3C4C5D1D2E1E2E3E4E5E605oooooooooooxooooox7.5-----ooooooxoooooo10ooooooooxooxoooxxo15oxxxx-------------12.55-----oooooxxoooooo7.5-----ooooooooooooo10-----ooooooooooooo15------------------255ooo
View original article
JOURNAL OF NEURAL ENGINEERING

留言 (0)

沒有登入
gif
format_indent_decrease