Non-invasive Decoding of the Motoneurons: A Guided Source Separation Method Based on Convolution Kernel Compensation With Clustered Initial Points

Introduction

A connectivity is generated between the motor-related areas of the brain during movement (Sporns et al., 2004; Rubino et al., 2006; Kim et al., 2017). A top-down structure in motor control is used particularly during an upper limb movement (with many degrees of freedom), where the brain and the muscles are functionally combined so that the muscle receives and electrically amplifies the resultant neural commands from the motor system (Haggard, 2008).

Over the last two decades, the non-invasive techniques for assessment of cortical and muscular activity have demonstrated significant progress and linear (e.g., coherence analysis) or non-linear (e.g., mutual information) methods have been used to analyze the relations between electromyographic (EMG), and electroencephalographic (EEG) signals (Chen et al., 2008; Meng et al., 2008). Due to the non-linear transfer functions of motor neurons, non-linear methods are more suitable for such a dependence analysis (Hashimoto et al., 2010; Ioannides and Mitsis, 2010) and for analysis of neural transfer functions between the central nervous system and pool of active motor units (MUs) (Negro and Farina, 2011; Gallego et al., 2015a,b).

Despite this progress in understanding of neural codes, the cortico-muscular coupling studies are still largely dependent on interferential EMG signal or its rectification for the assessment of motor neuron pool behavior (Gross et al., 2001; Schelter et al., 2009; Artoni et al., 2017). This assessment is not straightforward since the amplitude of surface EMG demonstrates several anatomical properties of the muscles, significantly interfering with the neural commands (Farina et al., 2010; Farina and Holobar, 2015, 2016). Indeed, amplitude or frequency of EMG are considerably affected by many factors, such as muscle anatomy, low-pass filtering of the subcutaneous tissues and MU Action Potential (MUAP) cancelation (Merletti and Farina, 2016). As result, they are not precisely related to the ongoing motoneuron activities and provide only a crude estimate of the central neural commands (so called neural drive) to the skeletal muscles (Farina et al., 2004, 2014a; Merletti and Farina, 2016). Therefore, the assessment of cortico-muscular connectivity and neuromuscular coupling via EEG–rectified EMG relationship is non-trivial and should be used with precaution.

Direct analysis of neural codes by decomposing the EMG (Webster et al., 2016), also known as neural decoding (Farina et al., 2014b; Karimimehr et al., 2017), is an alternative to EMG amplitude estimation (Gallego et al., 2017; Úbeda et al., 2017). It represents a paradigm shift, because it enables direct assessment of the neural drive to muscles (Farina and Holobar, 2015, 2016; Karimimehr et al., 2017). In fact, MU identification can be thought as a spike sorting algorithm (typical for computational neuroscience) applied to the outer layer of the human motor system (Balasubramanian and Obeid, 2011; Pani et al., 2016). The results of this MU identification can not only be used in EMG-EEG coupling analysis, but also in variety of research areas such as clinical neurophysiology for diagnosing neuromuscular disorders (Wheeler et al., 2006; Povalej BrŽan et al., 2017), sports and behavioral science (Merletti and Parker, 2004), movement science (Winter, 2009), robot-assisted rehabilitation (Savc et al., 2018), brain machine interface (Werner et al., 2016), and prosthesis control (Yoshida et al., 2010; Farina et al., 2014a).

In a typical experimental setup, EMG signals are detected using either conventional invasive intramuscular electrodes or non-invasive surface electrodes (Merletti and Parker, 2004). Both intramuscular (iEMG) and surface EMG (sEMG) signals include MUAP trains, superimposed into interferential signal patterns. Although iEMG provides some advantages, such as recording from deep muscles, it also has problems, such as discomfort and high selectivity. sEMG is thus a good alternative, particularly in sport sciences and studies of children (Ghaderi and Marateb, 2017). On the other hand, substantial MUAP superposition occurs in sEMG signals (Farina and Holobar, 2015, 2016). Also, surface MUAP shapes from different MUs are rather similar due to the volume conductor effects (Chen and Zhou, 2016). Thus, surface EMG decomposition is considered a very difficult task (Zhou and Rymer, 2004).

Variety of sEMG decomposition methods have been introduced in the literature. Generally, they either use shape-based algorithms, also called template matching (Xu et al., 2001; Gazzoni et al., 2004; Garcia et al., 2005; De Luca et al., 2006; Ren et al., 2006; Kleine et al., 2007; Nawab et al., 2008; Winslow et al., 2009; Siqueira Júnior and Soares, 2015) or the blind source separation (BSS) algorithm (Holobar and Zazula, 2004, 2007; Glaser et al., 2013; Ning et al., 2013, 2015; Chen and Zhou, 2016; Negro et al., 2016; Savc et al., 2018).

Although complementary fusion often results in more reliable findings (Durrant-Whyte, 1988), it is not yet sufficiently clear how to combine template-based and BSS sEMG decompositions, especially, as the aforementioned two methodological approaches differ substantially in the data model assumptions. For example, template-based algorithms assume a few recorded EMG channels and relatively low (or progressively increasing) muscle excitation levels in order to (progressively) identify the MUAP templates. They then use MUAP peel-off approach and combinatorial methods supported with artificial intelligence algorithms to identify each individual MU firing (Nawab et al., 2010). BSS approaches, on the other hand, model the mixing process of MUAPs in EMG, invert it and apply the inverse mixing procedure directly to the recorded EMG signals. For this reason, they do not rely on the initialization of MUAP templates, but do require relatively large number of recorded EMG channels to identify all the MU firings at once. In other words, BSS approaches estimate the MU filter directly in the space of MU spike trains and, once MU filter is estimated, apply it very efficiently (in terms of computational costs) to the recorded EMG signals, yielding all the MU firings at once (Holobar and Farina, 2014). Combination of these two fundamentally different sEMG decomposition approaches is non-trivial as it imposes complex methodological steps in which many parameters must be fine-tuned.

In this study, we focus on a hybrid sEMG decomposition approach that combines the advantages of both template-based and BSS decomposition approaches. Briefly, we use the single-pass density-based clustering algorithm to identify possible cluster representatives in different sEMG recording channels. Unlike traditional BSS-based decomposition algorithms, in which the initial search positions are randomly selected and checked on a trial-and-error basis, we use cluster representative samples as initial points. Similar approach was proposed in Ning et al. (2015), where k-means method was used to cluster MUAPs and, therefore, improve the initial estimate of MU filter in classical CKC approach. On the other hand, Chen et al. (2016, 2018a) and (Chen and Zhou, 2016) proposed the progressive FastICA framework along with the advanced peel-off and valley-seeking process that efficiently enhances the initialization of the MU filters and speeds-up the decomposition.

We then use the Kalman filter and a modified Gradient CKC (gCKC) algorithm (Holobar and Zazula, 2008) to further increase the rate of convergence of MU filter identification. Moreover, we introduce a new non-linear soft-thresholding algorithm to reduce the False Positives (FP) and False Negatives (FN) in MU spike post-processing. Thus, the proposed hybrid algorithm can be considered a “Guided Blind Source Separation” method.

The rest of the paper is organized as follows: in the next section, information about the signals and formulation of methods used in this study is presented. Section Results provides the results of the proposed method. The discussion is provided in section Discussion, along with conclusions.

Materials and Methods

Simulated Signals

A planar volume conductor model was used for generating synthetic sEMG signals (Farina and Merletti, 2001). Muscle, fat and skin tissues were used in the non-homogeneous and anisotropic volume conductor. It included an inimitable, and semi-infinite muscle layer with cross-section of 30 mm (transversal) × 15 mm (depth). The average fiber length, isotropic subcutaneous and skin layer thickness were 130 mm, 4 mm and 1 mm, respectively. Each MU had a random number of fibers uniformly distributed between 24 and 2,048 with the circular territories of 20 fibers/mm². The conduction velocities were normally distributed (4.0 ± 0.3 m/s). In the initial recruitment, each MU discharged at 8 pulses per second (pps) (Fuglevand et al., 1993). Its discharge rate increased linearly with excitation (0.3 pps per % of muscle excitation). A grid of 9 × 10 (9 columns, 10 rows) monopolar surface electrodes with 5-mm inter-electrode distances in both directions was simulated. Fifteen sEMG signals with length of 16 s were generated and sampled at 4,096 Hz. Muscle excitation was set to 10, 30, and 50% Maximum Voluntary Contraction (MVC), yielding 262, 388, and 446 active motor units. Noteworthy, regardless the decomposition algorithm used, we can identify only superficial motor units, whereas small and distant motor units contribute to the physiological noise. In addition to these simulated physiological noise, colored Gaussian zero-mean noise with the signal-to-noise ratio (SNR) of 10, 20, and 30 dB and the bandwidth of 20–500 Hz was added to the raw surface EMG signals. Overall, 45 simulated signals were analyzed in this study. The simulated signals are available online: https://doi.org/10.6084/m9.figshare.5808291.

The Proposed Algorithm

The structure of the proposed algorithm is depicted in Figure 1. Briefly, the signal was passed to the first-order band-pass Butterworth filter with the cut-off frequencies of 20 and 500 Hz in the forward and reverse direction. Furthermore, a whitening method was performed on sEMG signals (Thomas et al., 2006). Then, the whitened signal was used for initial point estimation as well as for the modified gCKC. Moreover, the convergence rate of the gCKC algorithm was improved by the Kalman filter. The detailed information about each aforementioned step is provided in the sequel.

FIGURE 1

Figure 1. The schema of the proposed decomposition algorithm. The sEMG signal is first filtered using band-pass and whitening filters. The template-based clustering algorithm is then used to identify the initial points for the modified gCKC algorithm. Such a clustering algorithm includes the segmentation and high-resolution alignment of the up-sampled signal. Then, a modified density-based clustering OPTICS algorithm is used to automatically locate cluster representatives (i.e., templates) in different recording channels. Such information is combined for different channels and the peak samples of the decimated templates are used to initialize the modified g-CKC algorithm. Finally, Kalman filtering and optimized peak finding is implemented to increase the efficiency of the algorithm and also to reduce the decomposition errors.

Whitening

Whitening, referred to as “convolutive sphering” (Thomas et al., 2006) was used as the first step of sEMG decomposition. Assume the following signal model,

$\begin{array}{ll}\text{X}=\text{W}\text{Z}& (1)\end{array}$

where Z is the N × M extended EMG signal matrix in which each row is a delayed repetition of one of EMG channels as proposed in Holobar and Zazula (2004) and Thomas et al. (2006) and each column corresponds to a time sample. In this study, the number of delayed repetitions of each signal is dependent to the sampling frequency and number of channels and in this work it was fixed to 30. The number of Whitening matrix W could be obtained, provided that the covariance matrix of X at time lag zero is equal to the identity matrix (Belouchrani et al., 1997):

$\begin{array}{ll}W=\text{U}{\text{D}}^{-\frac{1}{2}}{\text{U}}^T& (2)\end{array}$

where D is a diagonal matrix obtained by the eigenvalue decomposition of the covariance matrix of Z and U is the modal matrix:

$\begin{array}{ll}E\left\{{\text{Z}\text{Z}}^T\right\}={\text{UDU}}^T& (3)\end{array}$

with E{ZZ^T} denoting the covariance matrix of Z.

Initial Point Estimation

Similar to the shape-based sEMG decomposition approach in De Luca et al. (2006), the time samples of the sEMG signals were up-sampled to 10 KHz to increase the time resolution (Figure 1). The resampling was based on band-limited approximation of the original signal, produced by inserting zeros into the discrete Fourier transform of the waveform (Crochiere, 1979). We also included steps from Stearns and Hush (2011) to reduce the end effects caused by the zero-insertion.

Next, the signal detection threshold, required by our segmentation, was calculated as follows: the absolute values of the EMG samples were calculated and sorted (vector s). The square root values of cumulative sum of s were then calculated (ss). The maximum sample index i, where ss(i) multiplied by K exceeded the s(i) was found and s(i) was used as the detection threshold. In our study, K was empirically set to 4.0. The segmentation was then performed on each recording channel independently, using the estimated detection threshold on 2 ms intervals. Two-millisecond long signal segments, whose peak value was more than the detection threshold, were entitled as active segments (AcSs) (McGill et al., 2005).

Then, a high-resolution peak alignment method was used to align the detected AcSs on the highest peak (McGill and Dorfman, 1984). The time lags of the aligned AcSs were used as features for clustering. The density-based clustering was performed by using Ordering points to identify the clustering structure (OPTICS) algorithm (Ankerst et al., 1999; Daszykowski et al., 2002). This algorithm defines the clusters as areas of higher density. Among different density-based clustering algorithms, OPTICS was shown to be suitable for iEMG decomposition in which clusters have different dispersion (Marateb et al., 2011b). It is a single-pass clustering algorithm and unlike K-means or any other aggregative clustering method, does not require multiple runs on different predefined number of clusters (Daszykowski et al., 2002; Larose, 2005). The original OPTICS algorithm requires two parameters, namely the minimal number of points in a cluster MinPts and the Euclidean distance ε. A point p is called a “core point” if there are at least MinPts other points in its ε-neighborhood (N_ϵ(p)). A point p is directly density-reachable from another point q, if q is a core point and p is in the ε-neighborhood of q i.e. pϵN_ϵ(q). Point p is density-reachable from another point q if there is a chain of points p₁,…,p_n where p_i+1 is directly density reachable from p_i such that p₁ = q and p_n = p.

By setting the parameter ε to the maximum distance between any points in the dataset, each point in the dataset will be the core point. Thus, the modified OPTICS algorithm only requires the parameter MinPts. In our case, MinPts was set to 40 for 10 s long signal epochs. OPTICS detects varying-density clusters (i.e., clusters with different dispersions) in the dataset, such that most similar points (AcSs) become neighbors. In the other words, similar points are grouped together into hierarchical clusters. In this procedure, the first AcS is selected as the current point. Its reachability distance (RD) is set to INF. Its nearest neighbor in respect to RD is then set as the current point. Its RD is calculated as the smallest density reachable distance from the previous AcS (Marateb et al., 2011b). This process is repeated until all the points are processed in this way. Using this method, the entire feature space is transferred to a two-dimensional plot in which x-axis corresponds to the ordered points and y-axis depicts the RD of the ordered points. Each RD valley in such plot is a possible cluster, with the point corresponding to the minimum RD in a valley as the corresponding cluster representative. In our study, the possible valleys were automatically extracted using the adaptive method proposed by Marateb et al. (2011b, 2012).

Modified CKC

Denote by X(n, m) the (n,m)-th element of the signal matrix X (the m-th sample of channel n) and by X(m) the m-th column of matrix X. The rows in matrix X are EMG channels and their delayed repetitions (Holobar and Zazula, 2004; Thomas et al., 2006). In the CKC algorithm (Holobar and Zazula, 2004, 2007), the spike train ${\check{t}}_j$ of the j-th MU is estimated as

$\begin{array}{ll}{\check{t}}_j=R_{X{\check{t}}_j}\prime R_{XX}^{-1}X& (4)\end{array}$

where ${\text{R}}_{\text{X}\text{X}}^{-1}$ is the inverse of autocorrelation matrix of X, ${\text{R}}_{\text{X}{\check{t}}_j}$ is the cross-correlation vector of X and ${\check{t}}_j$ and the operator ′ denotes the matrix transpose. The cross-correlation vector ${\text{R}}_{\text{X}{\check{t}}_j}$ is unknown, but can be iteratively estimated by gCKC (Holobar and Zazula, 2007, 2008):

$\begin{array}{ll}{\text{R}}_{\text{X}{\check{t}}_{j,g}}={\text{R}}_{\text{X}{\check{t}}_{j,g-1}}+\eta (g){(f({\check{t}}_{j,g}){\text{X}}^{\prime })}^{\prime }& (5)\end{array}$

where ${\check{t}}_{j,g}$ is the j-th MU spike train, estimated in the g-th iteration, η is the step-size, f (x) = α(x) log (1 + x²) and α(x) is the attenuation coefficient introduced herein to improve the FP and FN errors in gCKC. The initial value of α is set to one but is then determined for each MU firing in each gCKC iteration. Namely, we observed that sometimes gCKC iterations erroneously amplify signal artifacts. Thus, after each iteration, we search for a spike detection threshold that yields the maximum Pulse-to-Noise Ratio (PNR) (Holobar et al., 2014).

$\begin{array}{ll}PNR\left(j,k\right)=10\times log\left(\frac{E({\check{t}}_{j,g}){|}_{{\check{t}}_{j,g}\ge thr}}{E({\check{t}}_{j,g}){|}_{{\check{t}}_{j,g}<thr}}\right)& (6)\end{array}$

where thr is the pulse detection threshold and PNR(j,k) is the Pulse-to-noise-ratio of the j-th MU in the k-th iteration. This is in fact a one-dimensional optimization problem and we efficiently solved it using a greedy search algorithm (Feo and Resende, 1995). The fitness function is the PNR while the optimization variable is the threshold thr. When this threshold is estimated, its 95% CI (Confidence Interval) is estimated and spikes within thr ± 95% CI of thr are detected (marginal spikes). The cluster representative (CR) of the spike samples surpassing the upper 95% CI threshold is formed and then if the corresponding sEMG AcS of a marginal spike has a high correlation with AcS of a CR, parameter α is set to one, otherwise it is set to 0.9. In this way, template-based methods are also used in the proposed modification of the previously introduced gCKC iterations (in addition to initial point estimation) (Holobar and Zazula, 2008). This operation is, in principle, equivalent to non-linear soft thresholding (Krzakala et al., 2016).

Kalman Filter

In our method, Equation (4) is substituted by expectation maximization Kalman filter (EM-Kalman filter), proposed by Bensaid et al. (2010). For this purpose, the gCKC is defined as the following state-space model joint to the EM-Kalman filter:

$\begin{array}{l}{\check{x}}_{k|k-1}={\text{F}}_k{x}_{k-1|k-1}+w_k\\ =\left[\begin{array}{cc}{\text{I}}_N& \eta (g)(\text{X}{f}^{\prime }(t_{j,g,k-1}))\\ 0& 1\end{array}\right]\left[\begin{array}{c}{\text{R}}_{\text{X}{t}_{j,g,k-1}}\\ 1\end{array}\right]\\ +w_k{\check{\text{t}}}_{j,g,k}(m)={\check{\text{x}}}_{k|k-1}\text{H}(m)& (7)\end{array}$

where ${\check{\text{x}}}_{k|k}$and k are estimated state and current Kalman filter iteration, respectively, g is related to the g-th iteration of the gCKC algorithm and $w_k=C_w{e}_k=C_w[e_k(1),e_k(2),\dots ,e_k(N),0]\prime$ is the normally distributed noise vector of size (N+1) × 1 with variance ${\sigma }_k^2$. C_w is a diagonal (N+1) × (N+1) covariance matrix of noise and $e_k=V_k\prime {\check{x}}_{k|k-1}$. According to the Bensaid et al. (2010), V_k is covariance matrix of ${\check{\text{x}}}_{k|k-1}$ and due to the whitening process, V_k could be considered identity matrix. Thus, ${\text{w}}_k={\text{C}}_{\text{w}}{\check{\text{x}}}_{k|k-1}$.

${\text{F}}_k=\left[\begin{array}{cc}{\text{I}}_N& \eta (g)(\text{X}{f}^{\prime }({\check{t}}_{j,g,k}))\\ 0& 1\end{array}\right]$ is the (N+1)×(N+1) state (or system) matrix with I_N denoting the identity matrix and N being the number of channels. ${\check{\text{x}}}_{k|k}$ and H(m) are defined by Equations (8) and (9):

$\begin{array}{ll}{\check{\text{x}}}_{k|k}=\left[\begin{array}{c}{\text{R}}_{\text{X}{\check{t}}_{j,g,k}}\\ 1\end{array}\right]& (8)\end{array}$

$\begin{array}{ll}\text{H}(m)=\left[\begin{array}{c}\text{X}(m)\\ 0\end{array}\right]& (9)\end{array}$

In this regard, a priori error covariance P and ${\sigma }_k^2$ could be defined by Equations (10) and (11), respectively.

$\begin{array}{ll}{P}_{k|k-1}=F_k{P}_{k-1|k-1}{F}_k\prime +C_w{C}_w\prime & (10)\end{array}$

$\begin{array}{ll}{\sigma }_k^2=\frac{1}{M-1}{\displaystyle \sum _{m=1}^M\left[H\prime (m)\left({\check{x}}_{k|k}{\check{x}}_{k|k}^{\prime }+P_{k|k}\right)H(m)-{\check{t}}_{j,g,k}\left(m\right)\right]}& (11)\end{array}$

where H = [H(1), H(2), H(3)…H(M)] and its size is (N + 1) × M. Kalman gain (K) is defined with Equation (12).

$\begin{array}{ll}{K}_k(m)=P_{k|k-1}H(m){\left(H\prime \left(m\right)P_{k|k-1}H\left(m\right)+{\sigma }_k^2\right)}^{-1}& (12)\end{array}$

We then update gCKC algorithm and error covariance as:

$\begin{array}{ll}{\check{\text{x}}}_{k|k}={\check{\text{x}}}_{k|k-1}+{\text{K}}_k({\check{t}}_{j,g,k}-{\check{t}}_{j,g,k-1})& (13)\end{array}$

$\begin{array}{ll}{P}_{k|k}=P_{k|k-1}-K_{k}H\prime P_{k|k-1}& (14)\end{array}$

When considering the state-space system matrix, h_m is the output matrix in the state-space model and K_k = [K_k(1), K_k(2), …, K_k(M)] of size (N+1) × M is the Kalman gain in the k-th iteration. ${\check{t}}_{j,g,k}$ in the first gCKC state is generated with the fixed point algorithm and one iteration of gCKC (Holobar and Zazula, 2007) and P_{k−1|k−−1} is initialized as I_N+1 i.e., the identity matrix. The estimation of the observation noise power ${\sigma }_k^2$is achieved by maximizing the Log-Likelihood function $log(P({\check{t}}_j|\check{\text{x}},{\sigma }_k^2))$ relative to ${\sigma }_k^2$ (Bensaid et al., 2010).

In a nutshell, we added additional noise factor in gCKC formulation and estimated it during each iteration since, practically, colored noise exists instead of white noise and gCKC algorithm, which works by correlation procedure, enhances colored noise (Li, 1990). This approach is used in other fields for finding noise parameters (Schwartz et al., 2015, 2016; Ge and Kerrigan, 2017) and also neuro spike decoding (Xue et al., 2017). In this light, EM-Kalman filter could prevent enhancing noise in ${\check{\text{x}}}_{k|k}$ in each iteration of gCKC and decrease the FP errors. This hypothesis will be, at least partially, confirmed by the results presented in section Results as the newly presented method detects substantially more MUs than CKC method.

Motor Unit Identification

The non-linear soft thresholding method mentioned in section Modified CKC, was applied to the reconstructed MU spike trains to identify MU firing times. Moreover, extracted MUs were monitored in terms of regularity of the firing times and their distinguishability from the background noise. MUs with the PNR <20 dB were excluded from further validations (Holobar et al., 2014). The following firing time parameters were extracted and used for MU exclusion: the number of inconsistent firings times n_I(number of pulses that are fired more than once in every 40 ms), the detection probability p_d(the probability of MU inter-spike interval having normal distribution as assessed by the method proposed by McGill Kevin, 1984) and the mean discharge rate (MDR). A MU was excluded when n_I was higher than 50, or p_d was lower than 50% or MDR was higher than 35 Hz. The underlying assumptions were the following: during constant force isometric contractions, MU firing rates lie between 5 and 25 Hz. Moreover, the inter-spike intervals have a unimodal distribution which is approximately Gaussian (Clamann, 1969; Andreassen and Rosenfalck, 1980). The pseudocode of the proposed HDsEMG decomposition algorithm is available as the Supplementary Material.

Performance Assessment

The identified MU firings were compared with the simulated firings. When at least 30% of the firings were time-locked to within ± 0.5 ms, a MUs was considered to be identified by our algorithm (Marateb et al., 2011a). For each identified MU, the accuracy of our new decomposition algorithm was assessed using the parameters of signal detection theory (e.g., TP (True Positive; correct firings), FN (False Negative; missed firings), and FP (False Positive; erroneous firings)). TP was the number of firings matching the simulated firings within ± 0.5 ms. FP was the number of firings not matching any simulated firing to within ± 0.5 ms. FN was the number of firings of the simulated MU that did not match any firings of the identified MU. Then, the performance indices Sensitivity (Se), Precision (Pr), and the Rate of agreement (RoA) were calculated (Holobar et al., 2010, 2014) (Table 1).

TABLE 1

Table 1. Decomposition validation parameters.

Moreover, the signal-to-interference ratio (SIR) was reported as the overall quality of the sEMG decomposition (Table 1). In fact, SIR estimates the percentage of the variance of the energy of the single-differential sEMG signals explained by the decomposition (Holobar et al., 2010). Also, the overall decomposability of a MU in the entire recording signals (cDI) was measured by the sum of norm of individual decomposability indices (DIs) (Holobar et al., 2010) (Table 1), calculated over all the EMG channels, normalized by the number of channels. The program was run on an Intel Core i7-8700 3.2 GHz CPU with 32 GB of RAM.

Statistical Analysis

Continuous variables were reported as mean ± std. The level of statistical significance was set to P = 0.05. The normality of the variables was assessed using the Shapiro-Wilk test. We used generalized estimating equation (GEE) method (Hardin, 2005) to model factors associated with the repeated responses (RoA and the running time of the algorithm). The paired t-test was used to identify if there is a significant difference between the estimated MDR and CoV in the decomposed signals, compared with those of the simulated MU firings. The Pearson correlation coefficient (r) was used to assess the association between two normal variables. The statistical analysis was performed using SPSS version 16 (SPSS for Windows, Released 2007, Chicago, IL, USA, SPSS Inc.).

Results

The average number of MUs, identified per contraction and SNR and also the accuracy of the decompositions of the assessed MUs are listed in Table 2. The accuracy of the clustering algorithm is also shown in this table as the average number of identified clusters. The cluster's representatives as well as their 20 nearby spikes in the RD plots, were compared with the simulated MU firings. When at least 75% of the firings in a cluster was time locked within 0.5 ms with firings of a simulated MU firings, the cluster was marked as correctly identified. Therefore, the reported number of identified clusters shows how many initial points were related to correct MU clusters. For comparison, the number of identified MUs and their decomposition accuracy of the gCKC (Holobar and Zazula, 2007) are also shown. Overall, our algorithm identified 16.41 ± 4.18 MUs with RoA of 84.00 ± 0.06%. The average number of modified gCKC iterations required to identify individual MU was 54.11 ± 3.06.

TABLE 2

Table 2. Decomposition accuracy.

The firing characteristics of the MUs, identified by the decomposition program and simulated MUs were very similar (Table 3), regardless the tested level of excitation and SNR.

TABLE 3

Table 3. The firing statistics of the decomposed MUs.

The sensitivity, RoA and precision of the proposed algorithm vs. PNR for different excitation levels and SNR values are depicted in Figure 2. Also, the distribution of the precision of the proposed decomposition algorithm compared with that of the gCKC for 30 dB SNR at 30% excitation and 20 dB SNR at 50% excitation is shown in Figure 3.

FIGURE 2

Figure 2. The sensitivity, precision, and Rate of Agreement (RoA) of the proposed algorithm vs. PNR (in dB). Representative plot is provided for each SNR level (10, 20, and 30 dB) at each simulated level of muscle excitation (10, 30, and 50% MVC).

FIGURE 3

Figure 3. The histogram of the precision of the proposed decomposition algorithm (black), compared to the one form gCKC (blue) for 30 dB SNR at 30% excitation level (left), and 20 dB SNR at 50% excitation level (right).

Figure 4 shows an example of the MUAP trains identified during seven seconds of a contraction at 10% MVC and 30 dB SNR. In this case, 14 MUs were identified with the average RoA of 89.00 ± 0.89 (%). Figure 5 shows the spatial distribution of Single-differential sEMG MUAP of MU 14 form Figure 4 as estimated by the spike-trigger averaging of the HDsEMG signals over the estimated MU firing times.

FIGURE 4

Figure 4. MU spike trains, identified from the simulated sEMG signal with 30 dB SNR and 10% excitation level (red), and the simulated firings (black). Each vertical line indicates one MU firing.

FIGURE 5

Figure 5. Single-differential sEMG MUAP waveforms of MUs from different sEMG channels. The corresponding MU was identified with accuracy >90%.

Discussion

Neural decoding is used to identify how the electrical activity of neurons generates responses in the brain (Jacobs et al., 2009). In the case of motor system, this can be performed either at the level of motor nerves or motor units. Although, the methods are very similar, the former is called “spike sorting” in computational neuroscience while the latter is known as “EMG decomposition,” or “decoding of the neural drive to the muscles” (Rey et al., 2015; Webster et al., 2016; Karimimehr et al., 2017). Such a decomposition algorithm could be used in variety of applications, including prosthesis control (Farina et al., 2014a) or robot-assisted neurorehabilitation (Savc et al., 2018). The structure of the spike sorting algorithms is, in principle, similar to that of iEMG decomposition methods, where the recording electrodes are close to the MU fibers. However, this is not the case for sEMG decomposition, where electrodes are located over the skin at a relatively large distance from the muscle fibers.

sEMG decomposition is considered a very difficult task when low selectivity of traditional sEMG electrodes, low-pass filtering of MUAPs, their overlapping, and shape similarities are considered (Farina et al., 2004). In this study, we proposed a new algorithm that combines the advantages of two different approaches to sEMG decomposition, namely the template-based and the BSS algorithms. Classification procedures are critical steps of the template-based algorithms. sEMG MUAPs are rather similar in shape and their classification is challenging (De Luca et al., 2006; Nawab et al., 2010). BSS algorithms, on the other hand, combines MU activities in the space of MU spikes in order to identify all the MU firings (Holobar and Farina, 2014). Many BSS algorithms cannot guarantee the identification of the same set of sources in multiple runs, because they rely (at least partially) on stochastic algorithms, for example on random initialization of vector ${\text{R}}_{\text{X}{\check{t}}_j}$in Equation (5).

Our algorithm, is fully deterministic since it identifies the initial points using a new clustering algorithm (Figure 1) rather than random initialization. Other CKC-based algorithms, use the time samples of the peaks of the sEMG signal as the possible starting point of the gradient-based optimization. However, sEMG contains highly overlapped MUAPs, and selected peaks of the signal usually correspond to the overlapped MU activity. This overlapped MU activity is then separated in several iterations of gCKC algorithm (Equation 5). On the other hand, if the initial ${\text{R}}_{\text{X}{\check{t}}_j}$ estimate is not related to any MU firings, but rather to the movement artifact or line interference, the CKC algorithm cannot amplify the MU spikes and the noise is amplified instead. Thus, being fully deterministic in our sense means that we know where a MU firing exists and we start the CKC gradient optimization from the initial point with a much higher MU identification potential than gCKC algorithm.

Similar aggregative clustering methods have been proposed in the past for sEMG decomposition (Ning et al., 2015). However, aggregative clustering methods, such as K-means or Fuzzy C-means do not have satisfactory “rerun” stability due to poor reproducibility based on random initialization (Jayaram and Klawonn, 2013). Furthermore, they may get stuck in local minimum owing to a two-stage iterative algorithm used to minimize the sum of point-to-center distances over all clusters. Finally, they have high computational complexity due to the need of running the algorithm several times with the pre-defined number of clusters (Duda et al., 2001). Therefore, we used a density-based OPTICS algorithm, which is a single-pass algorithm designed particularly for high-dimensional data (Ankerst et al., 1999; Daszykowski et al., 2002). This algorithm can automatically identify the optimal number of clusters when a proper automatic valley detection algorithm is used (Marateb et al., 2011b). Moreover, it adaptively identifies clusters of different sizes and variable dispersions and densities (Loh and Park, 2014).

In sEMG, there are many active MUs. This causes the overlap between MUAPs and thus the number of isolated segments reduces significantly, especially at higher contraction levels. Although each recording sEMG channel provides a different perspective of the MUAPs, the performance of the initial point detection generally drops at higher contraction levels (Table 2; The number of identified clusters variable). This is one of the limitations of the proposed algorithm. Using different spatial filters could improve the performance of this step of the algorithm which will be investigated in our future work.

Moreover, several cluster representatives identified by the OPTICS algorithm may be related to the same MU. This is why the average number of clusters identified by the OPTICS algorithm is higher than the number of identified MUs by the entire algorithm (Table 2). Since the MUAP shapes are very similar over the skin and largely submerged into the physiological noise, it is not possible to merge them unless their firing times are identified by the CKC gradient loops. On the other hand, due to the high overlapping between different MUAPs and the similarity between the shapes of different MUAPs in the sEMG signal, it is not optimal to use traditional classifiers (e.g., the minimum-distance classifier) based on the clusters obtained by the OPTICS algorithm. Thus, the combination of the OPTICS and gCKC algorithm is required, as proposed in our manuscript.

In our study, results of clustering were used as initial points in a modified g-CKC algorithm (Figure 1). In fact, our algorithm is a Guided-Source Separation (GSS) method with the following modifications introduced into the gCKC method:

(1) An adaptive soft-thresholding method was designed based on stochastic optimization, correlation analysis and estimation theory on the spike times as to reduce the FP (by suppressing noise spikes) and FN (by amplifying valid marginal spikes) error rates;

(2) The Kalman filtering was used to improve the convergence rate of the algorithm and to reduce the noise enhancement in gCKC iterations.

The performance of our algorithm was compared with that of the previously published gCKC algorithm (Table 3) since gCKC was shown to be accurate with different muscle architectures (Holobar et al., 2010; Marateb et al., 2011a; Webster et al., 2016) and different pathological conditions including essential tremor (Gallego et al., 2015a), Parkinson's disease (Holobar et al., 2012), diabetes (Watanabe et al., 2013), after stroke (Li et al., 2015), targeted muscle reinnervation (Farina et al., 2014c), and cleft lip surgery (Radeke et al., 2014). The previously introduced PNR, an efficient objective signal-based metric of gCKC decomposition accuracy (Holobar et al., 2014), was also reported for our decomposition results (Table 2, Figure 2). Although the sensitivity and RoA of the proposed algorithm is less than that of the gCKC algorithm, we identified more MUs at comparable Precision (Table 2; The average precision of 0.96 ± 0.05 compared with 0.95 ± 0.05 for gCKC). In fact, as mentioned in section Kalman filter, the Kalman filter decreases the FP errors and such errors have a direct effect on the Precision (Table 1). Therefore, FP error reduction goal was achieved by Kalman filter and the adaptive soft thresholding method. In our future work, we intend to focus on the reduction of FN errors. Moreover, we tested our algorithm with and without the Kalman filter. With the Kalman filter, the number of iterations was reduced in average down to 48% compared to the case without the Kalman filter. Moreover, when the Kalman filter was added, two more MUs were correctly identified, on average, whereas the running time of our algorithm was significantly reduced (P < 0.001; GEE).

The results of our statistical analysis showed that RoA was significantly associated with the SNR parameter i.e., the signal quality (P < 0.001; GEE) but not with the muscle excitation level i.e., signal complexity (P > 0.05; GEE). RoA was significantly and highly correlated with cDI at 10 dB and 20 dB SNR (P < 0.001; r > 0.7), but not at 30 dB SNR. This suggests that the proposed algorithm would work even if MUAPs were relatively similar in shapes when the quality of the signal was high. Moreover, there was no significant differences between the MDR and CoV values estimated from the decomposed sEMG signals and the corresponding simulated MU firings (P > 0.05; paired t-test). Thus, the decoded neural information can be considered reliable. Although measuring the MU identification accuracy indirectly, the PNR parameter was highly correlated with all the performance indices, namely with RoA, Sensitivity, and Precision (Figure 2; P < 0.001, r > 0.8). Thus, MUs with high PNR values can be considered as accurately identified. This is in agreement with the results in Holobar et al. (2014); Martinez-Valdes et al. (2016), and Watanabe et al. (2016).

In our study, simulated data was used for validation. Although, the used volume conductor model has been widely discussed in the literature (e.g., in McGill, 2004; Zazula and Holobar, 2005; Holobar et al., 2009; Zalewska, 2009), it could not completely resemble the experimental sEMG signals (De Luca and Nawab, 2011). Thus, further tests on experimental signals from different skeletal muscles in different experimental conditions (i.e., different contraction levels, force profiles, muscle geometries etc.) are required before the reported results can be generalized. This is a common challenge to all the reported HDsEMG decomposition algorithms as the results from one experimental setup cannot easily be generalized to all experimental conditions. In this study, we introduced the novel methodological steps and tested their efficiency against previously introduced CKC method that has been extensively validated over the past decade and applied to many (but not all) experimental conditions. Further generalization of the results reported herein exceeds the scope of this study and is left for the future work.

Indeed, the validation of the decomposition on experimental sEMG signals is, in principle, controversial. Using cross-checking, also known as two-source method, between the decomposition of concurrently recorded iEMG and sEMG signals, could only assess the accuracy of common MUs. Partitioning the sEMG channels into two groups and cross-checking them against each other, on the other hand, is also problematic as the information they convey is highly correlated (Webster et al., 2016). However, the objective measure PNR could be reported and used as the MU reliability index by the decomposition program (Holobar et al., 2014). Moreover, a novel validation approach proposed by Chen et al. could serve as a supplement to the conventional two-source methods (Chen et al., 2018b).

Finally, the proposed GSS algorithm has relatively low sensitivity to MUs, although it identifies more MUs than gCKC (Table 2). It identifies only a small portion of active MUs. This is a common limitation to all EMG decomposition approached introduced by now. The proposed algorithm is also limited to off-line analysis and cannot easily be converted to the online neural decoding algorithms (Glaser et al., 2007; Karimimehr et al., 2017). Its online implementation requires improvements in the rate of convergence of gradient-based optimization (Equation 5). For this purpose, the condition number of the Hessian matrix could be analyzed and proper diagonal scaling could be used to reduce the number of iterations (Beck, 2014), but this additional step has not been thoroughly investigated yet. Moreover, detection of firing-time inconsistency could be added to the OPTICS clustering algorithm in order to increase the number of MUs identified by the entire algorithm. We will address these methodological improvements in our future work.

In conclusion, we proposed a new framework for sEMG decomposition. The proposed algorithm is promising and a new offline tool in clinical neurophysiology.

Author Contributions

MRM, HM, SK, and AH participated in conceptualization, data curation, and validation. MRM, HM, SK, MAM, JK, and AH participated in formal analysis, investigation, and methodology. HM, MAM, and AH participated in project administration and resources. MRM, HM participated in software, and visualization. HM and AH supervised the project. MAM and AH acquired funding. MRM, HM, JK, and AH contributed to writing the original draft and SK, MAM revised the manuscript. All authors read and approved the final version of the manuscript and agreed for all aspects of the work.

Funding

This work was supported by the Spanish Ministry of Economy and Competitiveness- Spain project DPI2017-83989-R, by the University of Isfahan and by the Slovenian Research Agency (projects J2-7357 and L7-9421 and Programme funding P2-0041).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer PZ declared a past co-authorship with one of the authors, AH, to the handling Editor.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2019.00014/full#supplementary-material

References

Andreassen, S., and Rosenfalck, A. (1980). Regulation of the firing pattern of single motor units. J. Neurol. Neurosurg. Psychiatr. 43, 897–906. doi: 10.1136/jnnp.43.10.897

CrossRef Full Text | Google Scholar

Ankerst, M., Breunig, M. M., Kriegel, H.-P., and Sander, J. (1999). “OPTICS: ordering points to identify the clustering structure”, in ACM Sigmod Record: ACM, eds A. Deshpandeh, Z. Ives, and V. Raman (New York, NY: Association for Computing Machinery), 49–60. doi: 10.1145/304181.304187

CrossRef Full Text | Google Scholar

Artoni, F., Fanciullacci, C., Bertolucci, F., Panarese, A., Makeig, S., Micera, S., et al. (2017). Unidirectional brain to muscle connectivity reveals motor cortex control of leg muscles during stereotyped walking. Neuroimage 159, 403–416. doi: 10.1016/j.neuroimage.2017.07.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Balasubramanian, K., and Obeid, I. (2011). Fuzzy logic-based spike sorting system. J. Neurosci. Methods 198, 125–134. doi: 10.1016/j.jneumeth.2011.03.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Beck, A. (2014). Introduction to Nonlinear Optimization : Theory, Algorithms, and Applications with MATLAB. Philadelphia, PA: Society for Industrial and Applied Mathematics; Mathematical Optimization Society. doi: 10.1137/1.9781611973655

CrossRef Full Text | Google Scholar

Belouchrani, A., Abed-Meraim, K., Cardoso, J. F., and Moulines, E. (1997). A blind source separation technique using second-order statistics. IEEE Trans. Signal Process. 45, 434–444. doi: 10.1109/78.554307

CrossRef Full Text | Google Scholar

Bensaid, S., Schutz, A., and Slock, D. T. M. (2010). “Single microphone blind audio source separation using EM-Kalman filter and short+long term AR modeling,” in Latent Variable Analysis and Signal Separation, LVA/ICA 2010. Lecture Notes in Computer Science, Vol. 6365, eds V. Vigneron, V. Zarzoso, E. Moreau, R. Gribonval, and E. Vincent (Berlin; Heidelberg: Springer). doi: 10.1007/978-3-642-15995-4_14

CrossRef Full Text

Chen, C. C., Hsieh, J. C., Wu, Y. Z., Lee, P. L., Chen, S. S., Niddam, D. M., et al. (2008). Mutual-information-based approach for neural connectivity during self-paced finger lifting task. Hum. Brain Mapp. 29, 265–280. doi: 10.1002/hbm.20386

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, M., Holobar, A., Zhang, X., and Zhou, P. (2016). Progressive fastICA peel-off and convolution kernel compensation demonstrate high agreement for high density surface EMG decomposition. Neural Plast. 2016:3489540. doi: 10.1155/2016/3489540

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, M., Zhang, X., Chen, X., and Zhou, P. (2018a). Automatic implementation of progressive fastICA peel-off for high density surface EMG decomposition. IEEE Trans. Neural Syst. Rehabil. Eng. 26, 144–152. doi: 10.1109/TNSRE.2017.2759664

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, M., Zhang, X., and Zhou, P. (2018b). A novel validation approach for high-density surface EMG decomposition in motor neuron disease. IEEE Trans. Neural Syst. Rehabil. Eng. 26, 1161–1168. doi: 10.1109/TNSRE.2018.2836859

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, M., and Zhou, P. (2016). A novel framework based on fastICA for high density surface EMG decomposition. IEEE Trans. Neural Syst. Rehabil. Eng. 24, 117–127. doi: 10.1109/TNSRE.2015.2412038

PubMed Abstract | CrossRef Full Text | Google Scholar

Clamann, H. P. (1969). Statistical analysis of motor unit firing patterns in a human skeletal muscle. Biophys. J. 9, 1233–1251. doi: 10.1016/S0006-3495(69)86448-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Crochiere, R. E. (1979). “A general program to perform sampling rate conversion of data by rational ratios,” in Programs for Digital Signal Processing (Digital Signal Processing Committee of the IEEE Acoustics, Speech, and Signal Processing Society, ed C. J. Weinstein (New York, NY: IEEE Press, Programs), 8.2-1–8.2-7.

Daszykowski, M., Walczak, B., and Massart, D. L. (2002). Looking for natural patterns in analytical data. 2. Tracing local density with OPTICS. J. Chem. Inf. Comput. Sci. 42, 500–507. doi: 10.1021/ci010384s

PubMed Abstract | CrossRef Full Text | Google Scholar

De Luca, C. J., Adam, A., Wotiz, R., Gilmore, L. D., and Nawab, S. H. (2006). Decomposition of surface EMG signals. J. Neurophysiol. 96, 1646–1657. doi: 10.1152/jn.00009.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

De Luca, C. J., and Nawab, S. H. (2011). Reply to farina and enoka: the reconstruct-and-test approach is the most appropriate validation for surface EMG signal decomposition to date. J. Neurophysiol. 105, 983–984. doi: 10.1152/jn.01060.2010

CrossRef Full Text | Google Scholar

Duda, R. O., Hart, P. E., and Stork, D. G. (2001). Pattern Classification. New York, NY: Wiley.

PubMed Abstract | Google Scholar

Durrant-Whyte, H. F. (1988). Sensor models and multisensor integration. Int. J. Rob. Res. 7, 97–113. doi: 10.1177/027836498800700608

CrossRef Full Text | Google Scholar

Farina, D., and Holobar, A. (2015). Human- machine interfacing by decoding the surface electromyogram [Life Sciences]. IEEE Signal Process. Mag. 32, 115–120. doi: 10.1109/MSP.2014.2359242

CrossRef Full Text | Google Scholar

Farina, D., and Holobar, A. (2016). Characterization of human motor units from surface EMG decomposition. Proc. IEEE 104, 353–373. doi: 10.1109/JPROC.2015.2498665

CrossRef Full Text | Google Scholar

Farina, D., Holobar, A., Merletti, R., and Enoka, R. M. (2010). Decoding the neural drive to muscles from the surface electromyogram. Clin. Neurophysiol. 121, 1616–1623. doi: 10.1016/j.clinph.2009.10.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Farina, D., Jiang, N., Rehbaum, H., Holobar, A., Graimann, B., Dietl, H., et al. (2014a). The extraction of neural information from the surface EMG for the control of upper-limb prostheses: emerging avenues and challenges. IEEE Trans. Neural Syst. Rehabil. Eng. 22, 797–809. doi: 10.1109/TNSRE.2014.2305111

PubMed Abstract | CrossRef Full Text | Google Scholar

Farina, D., and Merletti, R. (2001). A novel approach for precise simulation of the EMG signal detected by surface electrodes. IEEE Trans. Biomed. Eng. 48, 637–646. doi: 10.1109/10.923782

PubMed Abstract | CrossRef Full Text | Google Scholar

Farina, D., Merletti, R., and Enoka, R. M. (2004). The extraction of neural strategies from the surface EMG. J. Appl. Physiol. 96, 1486–1495. doi: 10.1152/japplphysiol.01070.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Farina, D., Merletti, R., and Enoka, R. M. (2014b). The extraction of neural strategies from the surface EMG: an update. J. Appl. Physiol. 117, 1215–1230. doi: 10.1152/japplphysiol.00162.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Farina, D., Rehbaum, H., Holobar, A., Vujaklija, I., Jiang, N., Hofer, C., et al. (2014c). Noninvasive, accurate assessment of the behavior of representative populations of motor units in targeted reinnervated muscles. IEEE Trans. Neural Syst. Rehabil. Eng. 22, 810–819. doi: 10.1109/TNSRE.2014.2306000

PubMed Abstract | CrossRef Full Text | Google Scholar

Feo, T. A., and Resende, M. G. (1995). Greedy randomized adaptive search procedures. J. Glob. Optim. 6, 109–133. doi: 10.1007/BF01096763

CrossRef Full Text | Google Scholar

Fuglevand, J., Winter, D. A., and Patla, A. E. (1993). Models of recruitment and rate coding organization in motor-unit pools. J. Neurophysiol. 70, 2470–2488. doi: 10.1152/jn.1993.70.6.2470

PubMed Abstract | CrossRef Full Text | Google Scholar

Gallego, J. A., Dideriksen, J. L., Holobar, A., Ibanez, J., Glaser, V., Romero, J. P., et al. (2015a). The phase difference between neural drives to antagonist muscles in essential tremor is associated with the relative strength of supraspinal and afferent input. J. Neurosci. 35, 8925–8937. doi: 10.1523/JNEUROSCI.0106-15.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Gallego, J. A., Dideriksen, J. L., Holobar, A., Ibanez, J., Pons, J. L., Louis, E. D., et al. (2015b). Influence of common synaptic input to motor neurons on the neural drive to muscle in essential tremor. J. Neurophysiol. 113, 182–191. doi: 10.1152/jn.00531.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Gallego, J. A., Dideriksen, J. L., Holobar, A., Rocon, E., Pons, J. L., and Farina, D. (2017). “Neural control of muscles in tremor patients,” in Converging Clinical and Engineering Research on Neurorehabilitation II: Proceedings of the 3rd International Conference on NeuroRehabilitation (ICNR2016), October 18-21, 2016, Segovia, Spain, eds. J. Ibáñez, J. González-Vargas, J. M. Azorín, M. Akay and J. L. Pons (Cham: Springer International Publishing), 129–134. doi: 10.1007/978-3-319-46669-9_24

CrossRef Full Text | Google Scholar

Garcia, G. A., Okuno, R., and Azakawa, K. (2005). A decomposition algorithm for surface electrode-array electromyogram. IEEE Eng. Med. Biol. Magaz. 24, 63–72. doi: 10.1109/MEMB.2005.1463398

PubMed Abstract | CrossRef Full Text | Google Scholar

Gazzoni, M., Farina, D., and Merletti, R. (2004). A new method for the extraction and classification of single motor unit action potentials from surface EMG signals. J. Neurosci. Methods 136, 165–177. doi: 10.1016/j.jneumeth.2004.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Ge, M., and Kerrigan, E. C. (2017). Noise covariance identification for non-linear systems using expectation maximization and moving horizon estimation. Automatica 77, 336–343. doi: 10.1016/j.automatica.2016.11.011

CrossRef Full Text | Google Scholar

Ghaderi, P., and Marateb, H. R. (2017). Muscle activity map reconstruction from high density surface EMG signals with missing channels using image inpainting and surface reconstruction methods. IEEE Trans. Biomed. Eng. 64, 1513–1523. doi: 10.1109/TBME.2016.2603463

PubMed Abstract | CrossRef Full Text | Google Scholar

Glaser, V., Holobar, A., and Zazula, D. (2007). “An approach to the real-time surface electromyogram decomposition,” in 11th Mediterranean Conference on Medical and Biomedical Engineering and Computing 2007, eds T. Jarm, P. Kramar, and A. Zupanic (Heidelberg: Springer), 105–108. doi: 10.1007/978-3-540-73044-6_27

CrossRef Full Text | Google Scholar

Glaser, V., Holobar, A., and Zazula, D. (2013). Real-Time motor unit identification from high-density surface EMG. IEEE Trans. Neural Syst. Rehabil. Eng. 21, 949–958. doi: 10.1109/TNSRE.2013.2247631

PubMed Abstract | CrossRef Full Text | Google Scholar

Gross, J., Kujala, J., Hamalainen, M., Timmermann, L., Schnitzler, A., and Salmelin, R. (2001). Dynamic imaging of coherent sources: studying neural interactions in the human brain. Proc. Natl. Acad. Sci. U.S.A. 98, 694–699. doi: 10.1073/pnas.98.2.694

PubMed Abstract | CrossRef Full Text | Google Scholar

Haggard, P. (2008). Human volition: towards a neuroscience of will. Nat. Rev. Neurosci. 9, 934–946. doi: 10.1038/nrn2497

PubMed Abstract | CrossRef Full Text | Google Scholar

Hardin, J. W. (2005). “Generalized Estimating Equations (GEE),” in Encyclopedia of Statistics in Behavioral Science, eds B. S. Everitt and D. C. Howell (Hoboken, NJ: John Wiley & Sons Inc.). doi: 10.1002/0470013192.bsa250

CrossRef Full Text

Hashimoto, Y., Ushiba, J., Kimura, A., Liu, M., and Tomita, Y. (2010). Correlation between EEG-EMG coherence during isometric contraction and its imaginary execution. Acta Neurobiol. Exp. 70, 76–85.

PubMed Abstract | Google Scholar

Holobar, A., and Farina, D. (2014). Blind source identification from the multichannel surface electromyogram. Physiol Meas. 35, R143–R165. doi: 10.1088/0967-3334/35/7/R143

PubMed Abstract | CrossRef Full Text | Google Scholar

Holobar, A., Farina, D., Gazzoni, M., Merletti, R., and Zazula, D. (2009). Estimating motor unit discharge patterns from high-density surface electromyogram. Clin. Neurophysiol. 120, 551–562. doi: 10.1016/j.clinph.2008.10.160

PubMed Abstract | CrossRef Full Text | Google Scholar

Holobar, A., Glaser, V., Gallego, J. A., Dideriksen, J. L., and Farina, D. (2012). Non-invasive characterization of motor unit behaviour in pathological tremor. J. Neural Eng. 9:056011. doi: 10.1088/1741-2560/9/5/056011

PubMed Abstract | CrossRef Full Text | Google Scholar

Holobar, A., Minetto, M. A., Botter, A., Negro, F., and Farina, D. (2010). Experimental analysis of accuracy in the identification of motor unit spike trains from high-density surface EMG. IEEE Trans. Neural Syst. Rehabil. Eng. 18, 221–229. doi: 10.1109/TNSRE.2010.2041593

PubMed Abstract | CrossRef Full Text | Google Scholar

Holobar, A., Minetto, M. A., and Farina, D. (2014). Accurate identification of motor unit discharge patterns from high-density surface EMG and validation with a novel signal-based performance metric. J. Neural Eng. 11:016008. doi: 10.1088/1741-2560/11/1/016008

PubMed Abstract | CrossRef Full Text | Google Scholar

Holobar, A., and Zazula, D. (2004). Multichannel Blind source separation using convolution kernel compensation. IEEE Trans. Signal Process. 55, 4487–4496. doi: 10.1109/TSP.2007.896108

CrossRef Full Text | Google Scholar

Holobar, A., and Zazula, D. (2007). “Gradient convolution kernel compensation applied to surface electromyograms,” in Independent Component Analysis and Signal Separation: 7th International Conference, ICA 2007, London, UK, September 9-12, 2007. Proceedings, eds. M. E. Davies,. C. J. James, S. A. Abdallah, and M. D. Plumbley (Heidelberg: Springer), 617–624. doi: 10.1007/978-3-540-74494-8_77

CrossRef Full Text | Google Scholar

Holobar, A., and Zazula, D. (2008). “On the selection of the cost function for gradient-baseddecomposition of surface electromyograms,” in 2008 30th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (IEEE: Vancouver), 4668–4671.

PubMed Abstract | Google Scholar

Ioannides, A. A., and Mitsis, G. D. (2010). Do we need to consider non-linear information flow in corticomuscular interaction? Clin. Neurophysiol. 121, 272–273. doi: 10.1016/j.clinph.2009.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobs, A. L., Fridman, G., Douglas, R. M., Alam, N. M., Latham, P. E., Prusky, G. T., et al. (2009). Ruling out and ruling in neural codes. Proc. Natl Acad. Sci. U.S.A. 106, 5936–5941. doi: 10.1073/pnas.0900573106

PubMed Abstract | CrossRef Full Text | Google Scholar

Jayaram, B., and Klawonn, F. (2013). “Can fuzzy clustering avoid local minima and undesired partitions?,” in Computational Intelligence in Intelligent Data Analysis, eds. C. Moewes and A. Nürnberger (Heidelberg: Springer), 31–44. doi: 10.1007/978-3-642-32378-2_3

CrossRef Full Text | Google Scholar

Karimimehr, S., Marateb, H. R., Muceli, S., Mansourian, M., Mañanas, M. A., and Farina, D. (2017). A real-time method for decoding the neural drive to muscles using single-channel intra-muscular EMG recordings. Int. J. Neural Syst. 27:1750025. doi: 10.1142/S0129065717500253

PubMed Abstract | CrossRef Full Text | Google Scholar

Kevin, C. M. (1984). A Method for Quantitating the Clinical Electromyography. Ph.D., Stanford University.

Kim, B., Kim, L., Kim, Y.-H., and Yoo, S. K. (2017). Cross-association analysis of EEG and EMG signals according to movement intention state. Cogn. Syst. Res. 44, 1–9. doi: 10.1016/j.cogsys.2017.02.001

CrossRef Full Text | Google Scholar

Kleine, B. U., Van Dijk, J. P., Lapatki, B. G., Zwarts, M. J., and Stegeman, D. F. (2007). Using two-dimensional spatial information in decomposition of surface EMG signals. J. Electromyogr. Kinesiol. 17, 535–548. doi: 10.1016/j.jelekin.2006.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Krzakala, F., Ricci-Tersenghi, F., Zdeborova, L., Zecchina, R., Tramel, E. W., and Cugliandolo, L. F. (2016). Statistical Physics, Optimization, Inference, and Message-Passing Algorithms : École de Physique des Houches : Special Issue, 30 September-11 October 2013. Oxford: Oxford University Press.

Larose, D. T. (2005). “Hierarchical and k-means clustering,” in Discovering Knowledge in Data, ed D. T. Larose (Hoboken, NJ: John Wiley & Sons Inc.), 147–162. doi: 10.1002/0471687545.ch8

CrossRef Full Text

Li, W. (1990). Mutual information functions versus correlation functions. J. Stat. Phys. 60, 823–837. doi: 10.1007/BF01025996

CrossRef Full Text | Google Scholar

Li, X., Holobar, A., Gazzoni, M., Merletti, R., Rymer, W. Z., and Zhou, P. (2015). Examination of poststroke alteration in motor unit firing behavior using high-density surface EMG decomposition. IEEE Trans. Biomed. Eng. 62, 1242–1252. doi: 10.1109/TBME.2014.2368514

PubMed Abstract | CrossRef Full Text | Google Scholar

Loh, W.-K., and Park, Y.-H. (2014). “A survey on density-based clustering algorithms,” in Ubiquitous Information Technologies and Applications: CUTE 2013, eds Y-S. Jeong, Y-H. Park, C-H. Hsu and J. J. Park (Heidelberg: Springer), 775–780. doi: 10.1007/978-3-642-41671-2_98

CrossRef Full Text | Google Scholar

Marateb, H. R., McGill, K. C., Holobar, A., Lateva, Z. C., Mansourian, M., and Merletti, R. (2011a). Accuracy assessment of CKC high-density surface EMG decomposition in biceps femoris muscle. J. Neural Eng. 8:066002. doi: 10.1088/1741-2560/8/6/066002

PubMed Abstract | CrossRef Full Text | Google Scholar

Marateb, H. R., Muceli, S., McGill, K. C., Merletti, R., and Farina, D. (2011b). Robust decomposition of single-channel intramuscular EMG signals at low force levels. J. Neural Eng. 8:066015. doi: 10.1088/1741-2560/8/6/066015

PubMed Abstract | CrossRef Full Text | Google Scholar

Marateb, H. R., Rojas-Martínez, M., Mansourian, M., Merletti, R., and Villanueva, M. A. M. (2012). Outlier detection in high-density surface electromyographic signals. Med. Biol. Eng. Comput. 50, 79–89. doi: 10.1007/s11517-011-0790-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Martinez-Valdes, E., Laine, C. M., Falla, D., Mayer, F., and Farina, D. (2016). High-density surface electromyography provides reliable estimates of motor unit behavior. Clin. Neurophysiol. 127, 2534–2541. doi: 10.1016/j.clinph.2015.10.065

PubMed Abstract | CrossRef Full Text | Google Scholar

McGill, K. C. (2004). Surface electromyogram signal modelling. Med. Biol. Eng. Comput. 42, 446–454. doi: 10.1007/BF02350985

PubMed Abstract | CrossRef Full Text | Google Scholar

McGill, K. C., and Dorfman, L. J. (1984). High-resolution alignment of sampled waveforms. IEEE Trans. Biomed. Eng. 462–468. doi: 10.1109/TBME.1984.325413

PubMed Abstract | CrossRef Full Text | Google Scholar

McGill, K. C., Lateva, Z. C., and Marateb, H. R. (2005). EMGLAB: An interactive EMG decomposition program. J. Neurosci. Methods 149, 121–133. doi: 10.1016/j.jneumeth.2005.05.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Meng, F., Tong, K. Y., Chan, S. T., Wong, W. W., Lui, K. H., Tang, K. W., et al. (2008). Study on connectivity between coherent central rhythm and electromyographic activities. J. Neural Eng. 5, 324–332. doi: 10.1088/1741-2560/5/3/005

PubMed Abstract | CrossRef Full Text | Google Scholar

Merletti, R., and Farina, D. (2016). Surface Electromyography : Physiology, Engineering and Applications. Hoboken, NJ: IEEE Press; Wiley. doi: 10.1002/9781119082934

CrossRef Full Text | Google Scholar

Merletti, R., and Parker, P. A. (2004). Electromyography: Physiology, Engineering, and Non-invasive Applications. Hoboken, NJ: John Wiley & Sons. doi: 10.1002/0471678384

CrossRef Full Text | Google Scholar

Nawab, S. H., Chang, S.-S., and De Luca, C. J. (2010). High-yield decomposition of surface EMG signals. Clin. Neurophysiol. 121, 1602–1615. doi: 10.1016/j.clinph.2009.11.092

PubMed Abstract | CrossRef Full Text | Google Scholar

Nawab, S. H., Wotiz, R. P., and De Luca, C. J. (2008). Decomposition of indwelling EMG signals. J. Appl. Physiol. 105, 700–710. doi: 10.1152/japplphysiol.00170.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Negro, F., and Farina, D. (2011). Linear transmission of cortical oscillations to the neural drive to muscles is mediated by common projections to populations of motoneurons in humans. J. Physiol. 589, 629–637. doi: 10.1113/jphysiol.2010.202473

PubMed Abstract | CrossRef Full Text | Google Scholar

Negro, F., Muceli, S., Castronovo, A. M., Holobar, A., and Farina, D. (2016). Multi-channel intramuscular and surface EMG decomposition by convolutive blind source separation. J. Neural Eng. 13:026027. doi: 10.1088/1741-2560/13/2/026027

PubMed Abstract | CrossRef Full Text | Google Scholar

Ning, Y., Zhu, S., Zhu, X., and Zhang, Y. (2013). “A hybrid multi-channel surface EMG decomposition approach by combining CKC and FCM,” in: 2013 6th International IEEE/EMBS Conference on Neural Engineering (NER). (San Diego, CA), 335–338. doi: 10.1109/NER.2013.6695940

CrossRef Full Text | Google Scholar

Ning, Y., Zhu, X., Zhu, S., and Zhang, Y. (2015). Surface EMG decomposition based on K-means clustering and convolution kernel compensation. IEEE J. Biomed. Health Inform. 19, 471–477. doi: 10.1109/JBHI.2014.2328497

PubMed Abstract | CrossRef Full Text | Google Scholar

Pani, D., Barabino, G., Citi, L., Meloni, P., Raspopovic, S., Micera, S., et al. (2016). Real-time neural signals decoding onto off-the-shelf DSP processors for neuroprosthetic applications. IEEE Trans. Neural Syst. Rehabil. Eng. 24, 993–1002. doi: 10.1109/TNSRE.2016.2527696

PubMed Abstract | CrossRef Full Text | Google Scholar

Povalej BrŽan, P., Gallego, J., Romero, J., Glaser, V., Rocon, E., Benito-León, J., et al. (2017). New perspectives for computer-aided discrimination of parkinson's disease and essential tremor. Complexity 2017:4327175. doi: 10.1155/2017/4327175

CrossRef Full Text | Google Scholar

Radeke, J., Van Dijk, J. P., Holobar, A., and Lapatki, B. G. (2014). Electrophysiological method to examine muscle fiber architecture in the upper lip in cleft-lip patients. J. Orofac. Orthop. 75, 51–61. doi: 10.1007/s00056-013-0193-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, X., Hu, X., Wang, Z., and Yan, Z. (2006). MUAP extraction and classification based on wavelet transform and ICA for EMG decomposition. Med. Biol. Eng. Comput. 44, 371–382. doi: 10.1007/s11517-006-0051-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Rey, H. G., Pedreira, C., and Quian Quiroga, R. (2015). Past, present and future of spike sorting techniques. Brain Res. Bull. 119, 106–117. doi: 10.1016/j.brainresbull.2015.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Rubino, D., Robbins, K. A., and Hatsopoulos, N. G. (2006). Propagating waves mediate information transfer in the motor cortex. Nat. Neurosci. 9, 1549–1557. doi: 10.1038/nn1802

PubMed Abstract | CrossRef Full Text | Google Scholar

Savc, M., Glaser, V., Kranjec, J., Cikajlo, I., Matjačić, Z., and Holobar, A. (2018). Comparison of convolutive kernel compensation and non-negative matrix factorization of surface electromyograms. IEEE Trans. Neural Syst. Rehabil. Eng. 26, 1935–1944. doi: 10.1109/TNSRE.2018.2869426

PubMed Abstract | CrossRef Full Text | Google Scholar

Schelter, B., Timmer, J., and Eichler, M. (2009). Assessing the strength of directed influences among neural signals using renormalized partial directed coherence. J. Neurosci. Methods 179, 121–130. doi: 10.1016/j.jneumeth.2009.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwartz, B., Gannot, S., and Habets, E. A. (2015). Online speech dereverberation using Kalman filter and EM algorithm. IEEE/ACM Trans. Speech Lang. Process. 23, 394–406. doi: 10.1109/TASLP.2014.2372342

CrossRef Full Text | Google Scholar

Schwartz, O., Gannot, S., and Habets, E. A. (2016). An expectation-maximization algorithm for multimicrophone speech dereverberation and noise reduction with coherence matrix estimation. IEEE/ACM Trans. Speech Lang. Process. 24, 1495–1510. doi: 10.1109/TASLP.2016.2553457

CrossRef Full Text | Google Scholar

Siqueira Júnior, A. L. D., and Soares, A. B. (2015). A novel method for EMG decomposition based on matched filters. Res. Biomed. Eng. 31, 44–55. doi: 10.1590/2446-4740.0643

CrossRef Full Text | Google Scholar

Sporns, O., Chialvo, D. R., Kaiser, M., and Hilgetag, C. C. (2004). Organization, development and function of complex brain networks. Trends Cogn. Sci. 8, 418–425. doi: 10.1016/j.tics.2004.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Stearns, S. D., and Hush, D. R. (2011). Digital Signal Processing With Examples in MATLAB®. Boca Raton, FL: CRC Press.

Google Scholar

Thomas, J., Deville, Y., and Hosseini, S. (2006). Time-domain fast fixed-point algorithms for convolutive ICA. IEEE Signal Process. Lett. 13, 228–231. doi: 10.1109/LSP.2005.863638

CrossRef Full Text | Google Scholar

Úbeda, A., Del Vecchio, A., Sartori, M., Yavuz, U. S., Negro, F., Felici, F., et al. (2017). “Corticospinal coherence during frequency-modulated isometric ankle dorsiflexion,” in Converging Clinical and Engineering Research on Neurorehabilitation II: Proceedings of the 3rd International Conference on NeuroRehabilitation (ICNR2016), October 18-21, 2016, Segovia, Spain, eds J. Ibáñez, J. González-Vargas, J. M. Azorín, M. Akay, and J. L. Pons (Cham: Springer International Publishing), 135–140.

Google Scholar

Watanabe, K., Gazzoni, M., Holobar, A., Miyamoto, T., Fukuda, K., Merletti, R., et al. (2013). Motor unit firing pattern of vastus lateralis muscle in type 2 diabetes mellitus patients. Muscle Nerve 48, 806–813. doi: 10.1002/mus.23828

PubMed Abstract | CrossRef Full Text | Google Scholar

Watanabe, K., Holobar, A., Kouzaki, M., Ogawa, M., Akima, H., and Moritani, T. (2016). Age-related changes in motor unit firing pattern of vastus lateralis muscle during low-moderate contraction. Age 38:48. doi: 10.1007/s11357-016-9915-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Webster, J. G., Marateb, H. R., and McGill, K. C. (2016). “Electromyographic (EMG) decomposition,” in Wiley Encyclopedia of Electrical and Electronics Engineering, ed J. G. Webster (Hoboken, NJ: John Wiley & Sons Inc.). doi: 10.1002/047134608X.W8296

CrossRef Full Text

Werner, T., Vianello, E., Bichler, O., Garbin, D., Cattaert, D., Yvert, B., et al. (2016). Spiking neural networks based on OxRAM synapses for real-time unsupervised spike sorting. Front. Neurosci. 10:474. doi: 10.3389/fnins.2016.00474

PubMed Abstract | CrossRef Full Text | Google Scholar

Wheeler, K. R., Chang, M. H., and Knuth, K. H. (2006). Gesture-based control and EMG decomposition. IEEE Trans. Syst. Man Cyber. 36, 503–514. doi: 10.1109/TSMCC.2006.875418

CrossRef Full Text | Google Scholar

Winslow, J., Dididze, M., and Thomas, C. K. (2009). Automatic classification of motor unit potentials in surface EMG recorded from thenar muscles paralyzed by spinal cord injury. J. Neurosci. Methods 185, 165–177. doi: 10.1016/j.jneumeth.2009.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Winter, D. A. (2009). Biomechanics and Motor Control of Human Movement. Hoboken, NJ: John Wiley & Sons. doi: 10.1002/9780470549148

CrossRef Full Text | Google Scholar

Xu, Z., Xiao, S., and Chi, Z. (2001). ART2 neural network for surface EMG decomposition. Neural Comput. Appl. 10, 29–38. doi: 10.1007/s005210170015

CrossRef Full Text | Google Scholar

Xue, M., Wu, H., Zeng, Y., and Li, Y. (2017). “Unsupervised neuron spike decoding for macaque's finger position via EM algorithm,” in Information, Cybernetics and Computational Social Systems (ICCSS), 2017 4th International Conference (IEEE: Dalian), 599–604. doi: 10.1109/ICCSS.2017.8091485

CrossRef Full Text | Google Scholar

Yoshida, K., Farina, D., Akay, M., and Jensen, W. (2010). Multichannel intraneural and intramuscular techniques for multiunit recording and use in active prostheses. Proc. IEEE 98, 432–449. doi: 10.1109/JPROC.2009.2038613

CrossRef Full Text | Google Scholar

Zalewska, E. (2009). Insight into the motor unit activation and structure properties gained from EMG signal analysis. Clin. Neurophysiol. 120, 449–450. doi: 10.1016/j.clinph.2008.12.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Zazula, D., and Holobar, A. (2005). An approach to surface EMG decomposition based on higher-order cumulants. Comput. Methods Programs Biomed. 80(Suppl. 1), S51–60. doi: 10.1016/S0169-2607(05)80006-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, P., and Rymer, W. Z. (2004). MUAP number estimates in surface EMG: template-matching methods and their performance boundaries. Ann. Biomed. Eng. 32, 1007–1015. doi: 10.1023/B:ABME.0000032463.26331.b3

PubMed Abstract | CrossRef Full Text | Google Scholar