Motor Unit Identification From High-Density Surface Electromyograms in Repeated Dynamic Muscle Contractions

本文引入了MUAP变化的MU特异性测量，该测量可以直接在已识别的MU spikes中使用，并且不依赖于肌肉几何形状或关节角度的任何外部测量。PNR同样适用于动态下的分解效果评估。所引入的方法对MUAP变化和肌肉收缩速度表现出相当大的稳健性。

由于EMG混合过程的复杂性，通常由卷积多输入多输出系统建模，个体MU发放模式的识别在很大程度上仅限于等长肌肉收缩，其中肌肉和周围组织的几何形状不会发生显著变化。在这种条件下，由电极获取的MU动作电位（MUAP）是相对稳定的，并且主要的EMG非平稳性来自MU的募集和放电速率调制。在动态收缩中，MUAP 形状改变的两个因素：

• 肌肉缩短会影响动作电位在肌肉连接处终止的潜伏期，从而导致MUAP成分不传导；
• 肌肉纤维和摄取电极之间的距离随着肌肉缩短而变化。这两个因素都会影响检测到的MUAP形状。

---

2024/3/27

• 该算法的步骤是怎样的？

  该算法首先是选取一段小区间，这段小区间上的 $\Delta{H}$ 可以视为噪声，然后使用传统的 CKC 算法可以分解得到不同的 MU。关键在于怎样去追踪这个 MU。算法分为两步：

  • 先寻找 PNR 值接近的区间，这样的区间和初始区间有着相似的混合矩阵（这样的区间可不可以直接用 gCKC 算法分解出来？）。
  • 然后需要寻找别的区间，这涉及到 MUAP 的缓慢变化，因此需要设置一个加权函数 f，当权重在右边时，右边的相关向量占主要成分，所以会以右边的 MUAP 波形为基础去寻找下一个区间

---

Due to the complexity of the EMG mixing process, which is typically modelled by convolutive multiple-input-multiple-output system [12]–[15], identification of individual MU firing pattern has been largely limited to isometric muscle contractions, in which the geometry of the muscle and surrounding tissue does not change considerably. In such conditions, the MU action potentials (MUAPs), acquired by the uptake electrodes are relatively stationary and the main EMG nonstationarity comes from the MU recruitment and firing rates modulation.

上述表明，之前的 MU 分解都限制在等长收缩（肌肉的几何尺寸不变，周围的组织变化很小），在这种情况下，MUAP 是相对静态的，主要的非静态来自于 MU 的募集和发放频率的调制。如果连发放频率和募集都是恒定的，那么就说收缩是完全静态的。

hdEMG的动态建模

在等长非疲劳性肌肉收缩中，可以假设混合矩阵 $\mathbf{H}(n)$ 是静止的。在这种情况下，经典的 hdEMG 分解技术，如 CKC，可以应用于整个 hdEMG 信号。在缓慢的动态肌肉收缩和长时间的等长收缩中，由于几何变形或肌肉疲劳，MUAP 形状逐渐变化，因此，hdEMG 需要被划分为 $\mathbf{H}(n)$ 没有显著变化的足够短的时期。因此，经典的 hdEMG 分解技术需要独立地应用于每个信号历元。在典型的设置中，信号区间重叠 50%，以支持从不同的 hdEMG 区间识别的 MU 尖峰序列的不同片段能够有效合并。在中等和快速的动态收缩中，$\mathbf{H}(n)$ 变得非平稳，但在某种意义上仍然遵循局部性原理，即它逐渐变化，并且这些变化可以用连续函数来建模。

注意本文对 $H(n)$ 进行了重新建模：

$$\begin{gathered} \text{y(n)} =\bar{\mathbf{H}}\bar{\mathbf{t}}\left(n\right)+\Delta\mathbf{H}(n)\bar{\mathbf{t}}\left(n\right)+\boldsymbol{\omega}\left(n\right) \\ =\bar{\mathbf{H}}\bar{\mathbf{t}}(n)+\mathbf{z}(n)+\omega(n) \end{gathered}$$

这里的 $\Delta\mathbf{H}\left(n\right)\mathbf{\bar{t}}\left(n\right)$ 可以视作模型误差，$\bar{\mathbf{H}}$ 是 $\mathbf{H}(n)$ 在区间内的平均。在重复非疲劳实验中，$\Delta\mathbf{H}$ 和 $\mathbf{H}$ 是重复的，但是 $\mathbf{z}(n)$ 不是重复的，因为还取决于 $\bar{\mathbf{t}(n)}$ 。但是结果表明 $\mathbf{z}(n)$ 的相关矩阵能部分反应周期重复性。

MUST的估计

总结一下估计方法：分为区间内和区间外，区间内是寻找 PNR 值接近的区间（本质是寻找 $\mathbf{H}$ 接近的区间，$\mathbf{H}$ 代表的是 MUAP，$\mathbf{H}$ 接近说明 MUAP 接近，也就是说 MUAP 近似没有改变），然后重新求解参数。也就是说，只要把区间选的足够小，那么就可以把 hdEMG 视为静态的。

区间内估计

1. 第 $j$ 个 MU 的 spike train 的估计等式为：

$$\hat{t}_j\left(n\right)=\mathbf{c}_{\mathrm{t_j}\mathbf{y}}^{\mathrm{T}}\mathbf{C_y}^{-1}\mathbf{y}(n)$$

1. 在区间[a1, a2]上，通过梯度下降法得到 $\mathbf{c}_{t_j\boldsymbol{y}}$ ，进而估计 $\hat{t_j}\left(n_a\right)$
2. 重求 $\mathbf{c}_{t_j\boldsymbol{y}}$ 为 $\begin{aligned}\mathbf{c}_{t_j\boldsymbol{y},n_a}&=\sum_{n_a\in[a_1,a_2]}t_j\left(n_a\right)\boldsymbol{y}\left(n_a\right)\end{aligned}$ ，并重新估计 $\hat{t}_{j,n_a}\left(n\right)$
3. 寻找 PNR 接近的区间，在全部区间上重新求步骤 3 中的参数（PNR 会随着 $\Delta\mathbf{H}(n_a)$ 和 $\Delta\mathbf{H}(n_b)$ 之间差异的增大而减小，因此 PNR 只是一个代替指标，用来寻找 $\mathbf{H}(n_a)$ 和 $\mathbf{H}(n_b)$ 接近的区间。等式中的 PNR 值取决于所有 hdEMG 通道中第 j 个 MU 的 MUAP 的形状。因此，局部 PNR 允许对具有相对恒定的 $\mathbf{c}_{t_j\boldsymbol{y}}$ 和 $\mathbf{C}_y$ 的时间间隔进行MU特定的选择。

Δa was fixed to 1/5 of the contraction ramp (see the next section), whereas the number of decomposition runs in step 12 was fixed to 100

区间外估计

通过一个加权函数来使得区间中不同位置的 spike 被放大，最后把所有识别到的所有 spikes 都加在一块。移动区间 [a1, a2]，不断寻找。

实验

仿真数据集：

仿真数据集就是先逐段分解出 MUAP，通过重叠来匹配MUAP，从而得到一个 MUAP 库，这个库中的 MUAP 有一个动态的变化，使用这个库中的 MUAP 来合成 repeated dynamic EMG（注意分解的范式和合成的范式是不一样的）。

• 志愿者进行了80秒长的缓慢等速收缩，从肘部完全伸展开始，到肘部完全弯曲结束。肌肉收缩的范式包括完全伸展-屈曲斜坡。带通滤波范围：20-700 Hz
• 记录的信号被划分为10s长的历元，具有50%的重叠。通过 CKC 分解技术独立地分解每个历元，并成对比较所识别的MU尖峰序列，以识别整个80s记录信号上的MU发射模式。他们的 MUAP 是通过36个连续10秒长的 EMG 周期的尖峰触发平均值来估计的，这些周期重叠7.8秒。（注意，计算 MUST 和 MUAP 时的步长时不一样的，但是都是以10 s为一个 epoch）
• 仿真了 60 s的肌肉兴奋。在生成了 spike 后（MU firing rate 峰值为 35，意味着最小的 spike interval 为 30ms 左右），与 MUAP 库中的动态波形卷积得到仿真的 dynamic EMG。斜坡的数量设置为2、3、6和12，分别对应在15、10、5和2.5秒内实现肌肉完全伸展/屈曲。为了符合小的，远的 MU 常常淹没在生理噪声中这种情况，所有的 MUAP 被乘以了一个系数。

实验数据

• 5个人，4男1女
• 两个64通道电极，13×5，5列电极平行于肌纤维

在每次试验中，受试者进行重复的膝盖屈伸运动，将1公斤的重量固定在主腿的脚踝上。通过Fastrak运动跟踪器（Polheus，USA）测量膝关节角度，参考膝关节角度以稳定的速率从100度变化到160度，并在10秒内（中速）和20秒内（缓慢收缩）恢复。这里180度对应于膝盖的完全伸展。踝关节以100度角固定。每个试验包括八次（中速）或四次（慢速）连续重复所述膝盖运动，在整个试验过程中连续记录肌电图，并连续跟踪膝盖角度。至少 5 分钟间隔用来休息。

结果

• All the identified MUs extensively changed their MUAPs with muscle shortening.

• MUAP 的改变程度不仅特定于肌肉和收缩，还特定于 MU，这说明不同条件下 MUAP 的改变程度不同：

  the extent of MUAP changes in dynamic contractions is not only muscle and contraction, but also MU specific.

  但是， 有的 MU 即使 MUAP 改变了， 也不影响它的 MUST 的识别，能够在整个周期上识别出来，如后图所示：

  the impact of the MUAP changes on identified MU spikes is MU specific.

• 分解对MU是特异性的，有的 MU 可以在整个收缩周期被分解出来，但有的 MU 只能在某一段强度下分解出来。

  

• 脉冲噪声比 PNR，先前引入用于测量等长肌肉收缩中MU识别的准确性，也可应用于动态条件。PNR反应了 spike train 中 spike 的高度，spike 的高度越高，那么 PNR 值就越大。。sensitivity 和 precision 与 PNR 成正相关，false alarm 和 PNR 值成负相关。PNR > 30 dB 时，灵敏度和准确度 > 90%，误报率 < 1%，这个在等长收缩时的特性在动态收缩中依然保留了下来，尽管明确指出 PNR 值在肌肉收缩、MUAP 改变时会受到影响。因此，PNR ≥ 30 dB 的标准应使用在局部 PNR 中，而在灵敏度 > 90%的 MU 中经常观察到相当低的全局PNR值（图6）。

• 发力斜坡的长度并不能决定分解出的 MU 个数、PNR 值，以及灵敏度和准确度，说明这种算法对肌肉的收缩速度是鲁棒的。但是兴奋程度越高识别到的 MU 越多。

• 与所提出的循环平稳CKC相比，经典的 CKC 方法识别出更多的 MUs ，但灵敏度明显较低(图5)。这与图 3 中识别出的 spike 序列的例子一致。

• 实验范式设置参考：MU 分解对曲线重复次数的敏感性是特定于肌肉和收缩水平的。一般来说，运动速度越高，所需的运动重复次数就越大。值得注意的是，运动重复不需要具有与伪码1的步骤6中实现的MU跟踪相同的持续时间，该MU跟踪在每个运动重复中自动确定用于MU尖峰识别的最佳时间支持。

• 与传统 CKC 的区别：传统的 CKC 是作用于整个信号，通常不能完整地识别出所有的 spikes。CKC 算法的循环次数一般设置为100。而 csCKC 算法识别出的 spike 更多。

  

• PNR 这个标准可以用于任何复杂的肌肉收缩。

还有哪些工作可以做

• 没有试图解释同心和偏心期MU发射率的差异，也没有试图解释股外侧肌和股直肌中MU发射模式之间的差异。这些解释需要大量的重复肌肉收缩，超出了本研究的范围。
• 未来需要进一步的研究来评估所述方法在不同骨骼肌中识别MU的适当性，并分析其对肌肉收缩重复次数的敏感性。
• Validation of dynamic surface EMG decomposition is currently an open research problem.