Motor Unit-Driven Identification of Pathological Tremor in Electroencephalograms

Background: Traditional studies on the neural mechanisms of tremor use coherence analysis to investigate the relationship between cortical and muscle activity, measured by electroencephalograms (EEG) and electromyograms (EMG). This methodology is limited by the need of relatively long signal recordings, and it is sensitive to EEG artifacts. Here, we analytically derive and experimentally validate a new method for automatic extraction of the tremor-related EEG component in pathological tremor patients that aims to overcome these limitations. Methods: We exploit the coupling between the tremor-related cortical activity and motor unit population firings to build a linear minimum mean square error estimator of the tremor component in EEG. We estimated the motor unit population activity by decomposing surface EMG signals into constituent motor unit spike trains, which we summed up into a cumulative spike train (CST). We used this CST to initialize our tremor-related EEG component estimate, which we optimized using a novel approach proposed here. Results: Tests on simulated signals demonstrate that our new method is robust to both noise and motor unit firing variability, and that it performs well across a wide range of spectral characteristics of the tremor. Results on 9 essential (ET) and 9 Parkinson's disease (PD) patients show a ~2-fold increase in amplitude of the coherence between the estimated EEG component and the CST, compared to the classical EEG-EMG coherence analysis. Conclusions: We have developed a novel method that allows for more precise and robust estimation of the tremor-related EEG component. This method does not require artifact removal, provides reliable results in relatively short datasets, and tracks changes in the tremor-related cortical activity over time.

To the best of our knowledge, all of the existing studies assessed corticomuscular coupling by computing the coherence between cortical activity, recorded with EEG or magnetoencephalograms (MEG), and an estimate of muscle activity derived from the surface EMG (2, 3, 6-13, 15, 16). Since the coherence function reveals a linear relationship between two signals at a given frequency (17), coherence at the tremor frequency is assumed to indicate tremor-related EEG activity. Although robust to noise at different frequencies, coherence only provides an indirect measure of corticomuscular coupling, and does not enable tracking changes in tremor properties over a short time scale. Furthermore, it requires off-line processing of relatively long EEG and EMG recordings, which need to be cleaned of artifacts beforehand. This limits the comparison of the tremor-related cortical activity across conditions and diseases.
Besides the coherence function, the cortical tremor component could also be potentially identified using blind source separation (BSS) algorithms. For example, Delorme et al. (18) identified, using independent component analysis (ICA) techniques, a number of components in the EEG activity of healthy subjects performing a working memory task. However, no group has demonstrated the feasibility of separating the tremor component from other brain activity. Notably, all the ICA/BSS algorithms proposed so far build on general assumptions about the EEG properties such as independence of the identified components, and would not exploit the specific characteristics of tremor, such as the relationship between sensorimotor cortical activity and muscle activity. As a result, they would suffer from large inter-trial and inter-subject variability in convergence toward the specific (tremor-related) EEG component.
With the exception of Gallego et al. (15), where we used cumulative motor unit spike trains (CST) to characterize the neural drive to the muscle, the authors of all of the studies mentioned above used the rectified surface EMG as an estimator of muscle activity. However, recent studies have shown that the CST of several (e.g., ≥5) motoneurons that innervate a muscle provide a more accurate representation of the synaptic inputs to a motoneuron pool than the EMG envelope (19)(20)(21). In addition, rectification of the surface EMG may or may not enhance the detection of synaptic inputs to the pool depending on the muscle contraction level (22). Indeed, as shown in Farina et al. (22) rectification is preferable over the raw EMG only at low contraction levels. Therefore, methods based on the CST rather than on the traditional surface EMG analysis are likely to identify tremor-related cortical activity more reliably.
In this study, we present and validate a novel method to identify tremor-related cortical activity. This method builds on the assumption that the tremor-related EEG component is stochastically phase-locked to the motor unit firings in a tremulous muscle. This implies close to linear relationship between the cortical activity and the firings of the pool of motor units that form a muscle (15,21,23). This relationship has been experimentally demonstrated before by the existence of EEG-EMG coherence at the tremor frequency (2, 3, 6-13, 15, 16, 24).
In our method, motor unit spike trains are identified from non-invasively recorded multichannel surface EMG recordings using the Convolution Kernel Compensation (CKC) technique (25)(26)(27). These firings are then used to construct the phaselocked estimator of the tremor-related activity in EEG. By using simulated data, we show that our method tracks reliably the tremor activity for a wide range of physiologically realistic conditions. We then apply this method to recordings from nine essential tremor (ET) and nine Parkinson's disease (PD) patients, and show that the method outperforms the traditional coherence approach when detecting tremor-related cortical activity.

MULTICHANNEL EMG-DRIVEN IDENTIFICATION OF TREMOR-RELATED ACTIVITY
We first assume that EEG signals are linear mixtures of brain oscillations (rhythms) and noise (Figure 1). Then, using this data model, we build a linear minimum mean square error (LMMSE) estimator of the cortical tremor activity that exploits the large synchronization of motor unit firings in pathological tremor.

Data Model
Assume the mixing model depicted in Figure 1, and denote the M EEG channels by y (n) = y 1 (n) , y 2 (n) . . . y M (n) T , where the n-th sample of the i-th channel appears in the i-th row of y (n). The model inputs, s j (n), represent the brain rhythms in the EEG. For example, in a normal condition, these inputs would reflect the alpha, beta or gamma rhythms. As in other BSS/ICA EEG studies, artifacts, such as blinking artifacts, would also be considered a model input. In the case of pathological tremor, one or more of these inputs reflects tremor activity. The mixing model in Figure 1 can be expressed in a matrix form as.
contains blocks of L samples from J sources and the noise vector ω (n) = ω 1 (n) , ω 2 (n) . . . ω M (n) T is modeled as a zero-mean represent the non-tremoric brain rhythms in the EEG. The impulse responses a mj project the j-th source s j (n) to the m-th EEG channel y m (n), whereas ω m (n) denotes the additive noise in the m-th EEG channel.
ergodic Gaussian process, spatially and temporarily independent of the activity in s(n). The mixing matrix contains stationary impulse responses a mj = a mj (0) · · · a mj (L − 1) that convolve the j-th input source and add the result to the m-th EEG channel.
The model described in Equation (1) is typically underdetermined, with more inputs than measurements. However, low energy inputs can always be modeled as physiological noise and, thus, included in ω (n). We can also extend y (n) by adding F-1 delayed repetitions of each EEG channel: y (n) = y 1 (n), y 1 (n − 1) . . . y 1 (n − F + 1), y 2 (n) . . .
This increases the number of outputs in model (1) and, more importantly, compensates potentially different time delays of same source in different EEG channels (see EMG-Driven Decomposition of EEG section for details). This is important because the same rhythm may be present at different EEG electrodes at close but different time lags (this assumption is further confirmed by the results in Results section). The extended input vector s and mixing matrix A now change to: For the purpose of mathematical derivations in Appendix, we will further represent the EEG recordings y(n) as analytic signals. Note that this does not alter the assumed mixing model, and can be readily fulfilled by applying Hilbert transform to the EEG signals.

EMG-Driven Decomposition of EEG
By temporarily neglecting the impact of noise ω (n), the LMMSE estimator of the input s j (n) is given by (25) where superscript H denotes conjugate transpose, c s j y =E s j (n)y(n) is the cross-correlation vector between s j (n) and y(n), and c s j s =E s j (n)s(n) is the cross-correlation vector between s j (n) and s(n). Matrices C s and C y denote the correlation matrices of s(n) and y(n), respectively. This LMMSE estimator is Bayesian optimal in the minimum square error sense, also in the presence of noise, but requires a priori knowledge of the cross-correlation vector c s j y . In experimental conditions c s j y is not known and needs to be estimated from the measurements.
To obtain an estimate of c s j y , we developed a method that is based on the assumption that, in an affected muscle, the motor unit firings are phase-locked to the tremor-related cortical activity. This assumption follows from the observations of significant coherence (linear relationship) between the EEG and the rectified EMG (1-10, 13, 14, 16) or, more directly, between the EEG and the motor unit CST (15). Remarkably, the number of motor units needed to represent population behavior during tremor is not very large, which in practice will favor the implementation of the proposed method. For example, in Negro and Farina (21) and Gallego et al. (15), the coherence between the CST and the EEG was nearly maximal with only five motor units and did not increase significantly when more motor units were added to the CST.
To obtain an estimate of c s j y , we first define the firing pattern of the r-th tremor-affected motor unit: where δ(τ ) denotes the unit-sample impulse, d r is the common time lag (in samples) of the pulses in t r (n) due to the transmission from motor cortex, d rk is the intrinsic motor unit firing time variability (in samples), and is frequently defined as the k-th realization of Gaussian random variable d rk ∼ N(0, σ d r ), f r is the motor unit firing frequency, f s is the sampling frequency and K r is the number of firings in the observed time interval. Note that, in order to simplify analytical derivations, we ignored the doublets in the motor unit firing pattern (8). This has small impact on the results presented herein, as doublets can always be modeled as two different instances of the k-th motor unit firing with different d rk values. Furthermore, as shown in Results section, when respecting its physiologically induced range, d rk has almost negligible impact on the proposed tremor estimation.
The cross-correlation between the EEG component s j (n) and y (n) can be estimated as (see derivation in the Appendix) t r (η) y(η) (9) where N is the number of signal samples and factor α is defined in the Appendix). Due to the amplitude ambiguity of source components extracted by BSS (28), the scalar factor α can be neglected and the EEG-component that is phased-locked to the firings of the J motor units expressed aŝ By knowingŝ j (n) cross-correlation vectorĉ s j y in (9) may be recalculated asĉ whereŜ j f = F ŝ j (n) is Fourier transform ofŝ j (n) and g (x) = x * |x| denotes the element-wise product of element x with its absolute value. Equations (10,11) are then iteratively recalculated until the convergence is reached. In our study, iterations were stopped when the second norm ofŝ j (n) changed for <0.1%. These iterations emphasize the spectral peaks in the extracted EEG component and suppress the wideband frequency components with low energies.
Note that (9) is only true when none of the other oscillatory inputs in s(n) is a higher harmonic of the input s j (n). In the opposite case, (10) would identify both the tremor-related EEG component and its higher harmonics as one joint input (see explanations in the Appendix).
The presented method still needs a good approximation of the motor unit firing times,t r (n), in order to accurately estimateĉ s j y . This can be obtained from high-density surface EMG recordings using the CKC decomposition technique (25,27,(29)(30)(31), which has already been demonstrated to be highly robust to high levels of motor unit synchronization (27).

SIMULATIONS AND EXPERIMENTAL RECORDINGS
The presented method was validated on a set of synthetic signals and on experimental recordings from 18 tremor-affected patients.

Simulations
First, we tested the proposed method in a simple model that generated EEG-like oscillations as mixtures of sinusoidal sources with time-varying amplitudes. Our goal was to test the method's ability to accurately reconstruct such EEG-like sources from their convolutive mixtures, and to study its sensitivity to its three main parameters: extension factor F in (4), motor unit firing variability d rk and signal-to-noise ratio (SNR). The simulated signals comprised 10 (J = 10) mutually orthogonal sinusoids s j (n) and their first higher harmonics as input signals: where a(n) is an amplitude modulator generated by filtering white noise with a second order low-pass Butterworth filter with cut-off frequency of 1 Hz. The amplitude B was set equal to 1, whereas the amplitude of the first harmonic, H 1 , was varied across simulations and was set to 0, 0.2, 0.4, 0.6, 0.8, and 1. This way we simulated the experimentally observed ratios between the basic tremor frequency and its first harmonic (9,13,32). The frequency f j of the oscillatory inputs was set to 5+j/2 Hz with j = 1, 2 . . . 10, and the phase φ j was randomly selected from the interval [0, 2π]. The sampling frequency was set to 1024 Hz and each simulation lasted 30 s. We assumed that the first oscillatory input, s 1 (n), represented the tremor-related component we wanted to detect. Next, we simulated the spike trains of ten motor units, t r (n) with r = 1, 2 . . . 10, by finding the local maxima of the first generated oscillatory source, s 1 (n). We imposed a corticospinal delay d r = 10 ms in Equation (8) to simulate the physiological delays (due the transmission from motor cortex to the output of the motoneuron pool) between the tremor-related EEG component s 1 (n) and the simulated motor unit firing patterns t r (n). Finally, the firing variability d rk of each individual motor unit t r (n) was modeled as Gaussian random variable d rk ∼ N(0, σ d r ) (33, 34) (see Appendix). The standard deviation σ d r was set to 0, 10, and 20% of the average inter-spike interval in t r (n), respectively (35)(36)(37)(38). Ten Monte Carlo simulation runs were performed for each value of H 1 and σ d r , resulting in a total of 180 (10 × 6 × 3) simulation runs.
We computed 15 synthetic channels y(n) in Equation (1), with the mixing matrix A having dimensions 15 × 50. The length of the impulse responses a mj (l) in (3) was set to L = 5 samples. To simulate different delays in the representation of the oscillatory inputs s j (n) across the 15 channels, each element a mj had one randomly selected element set to a non-zero random value from a normal distribution N(0, 1), whereas the remaining four elements were set to zero. In each simulation run, five different realizations of zero-mean random Gaussian noise were added to the measurements, with the SNR ranging from 0 to 20 dB in steps of 5 dB. This resulted in 900 (180 × 5) simulated sets of signals. Note that our simple generative model does not realistically represent actual EEG signals; our goal was to use the simulation results to identify the range of parameter values that was adequate for the experimental data analysis. All the simulations were performed in Matlab version 8.6.0.267246 (R2015b).

Experimental Recordings
We recorded data from nine ET patients (four females and five males; age, mean ± SD: 70 ± 6 years, range 61-79 years) and nine PD patients (three females and six males; age, mean ± SD: 64 ± 14 years, range 44-88 years) at Hospital 12 de Octubre, Madrid, Spain. In the ET patients, tremor severity ranged from mild (two patients) to severe (three patients), with a mean score of 36 ± 12 (mean ± SD; range 20-51) according to the Fahn-Tolosa-Marin scale (39). In the PD patients, tremor severity also ranged from mild (five patients) to severe (two patients), with a mean score of 12 ± 6 (mean ± SD; range 5-23) according to the UPDRSIII scale (40). All the participants included in the study gave their written informed consent after full explanation of the procedure. The study, which was conducted in accordance with the principles of the Helsinki declaration of 1975, was approved by the ethical standards committee on human experimentation at the University Hospital "12 de Octubre" (Madrid).
Right wrist kinematics were recorded with inertial measurement units (IMUs) comprising three-axis accelerometers, magnetometers and gyroscopes (Tech MCS, Technaid S.L., Madrid, Spain). These sensors were fixed with surgical tape over the hand dorsum and the distal third of the forearm (on the dorsal side, close to the wrist), respectively, with one of their axes aligned to that of the wrist. Data were sampled at 100 Hz by a 12-bit A/D converter and low pass filtered (<20 Hz). Wrist kinematics was assessed as the difference between the measured accelerations in the axis parallel to the wrist (41).
Surface EMG signals were recorded from the right wrist flexors and extensors with 13 × 5 electrode grids (LISiN-OT Bioelettronica, Torino, Italy, 8 mm interelectrode distance). The electrode grids were centered over flexor carpi radialis and extensor digitorum communis, respectively. Before the placement of the electrode grid, the skin was lightly abraded using abrasive paste (Meditec-Every, Parma, Italy) and cleansed afterward. Electrical conductivity was ensured by filling each of the electrodes in the grids with conductive gel (Meditec-Every, Parma, Italy). A soaked bracelet placed around one of the wrists was used as reference. The surface EMG signals were amplified as bipolar recordings along the direction of the fibers, band-pass filtered (3 dB bandwidth, 10-750 Hz), and sampled at 2,048 Hz by 12-bit A/D converter (LISiN-OT Bioelettronica, Torino, Italy). We synchronized the EEG, EMG, and IMU recordings using a common clock signal, which was fed into each acquisition systems. The rising edge of the first and last clock signal pulses were identified using a purposely-developed Matlab script. Data were then resampled to 2,048 Hz using a routine that incorporated an anti-aliasing filter.

Data Analysis
In the experimental recordings, individual motor unit firing patterns were identified from the multichannel surface EMG using the CKC algorithm (25,27,29). The pulse-to-noise ratio metric (PNR) (42) was used to assess the accuracy of firing estimation for each identified motor unit. Only reliably identified motor units (PNR > 30 dB; accuracy of motor unit firing estimation >90%) were used for further analysis (42), whereas all the remaining motor units were discarded.
We estimated the tremor EEG component with Equation (7), using the firings of all the identified (experimental recordings) or simulated (simulations) motor units to estimate the crosscorrelationĉ s j y between the oscillatory components and the EEG signals, as defined in Equation (10). We tested different extension factors, from F = 1 to F = 15. Due to the amplitude ambiguity (see Appendix), the estimateŝ 1 (n) was further normalized to yield a unit norm. Finally, Equations (10, 11) were iteratively applied until the convergence criterion was reached (the second norm ofŝ j (n) changed for <0.1%).
The delay between the motor unit CST and the estimated tremor EEG componentŝ 1 (n) was estimated as the argument of the cross spectrum at the basic tremor frequency f b (17,43).
In the simulations, we further assessed the accuracy of our method by computing the normalized mean square error (NMSE; Equation 13) and the cross-correlation coefficient (CC) between the simulated and estimated tremor inputs,ŝ 1 (n) and s 1 (n), after both signals were aligned in time: In the experimental data, we compared the coherence between the extracted tremor component and the CST in each muscle to the coherence between the CST and the Laplacian-filtered EEG (15) (without any artifact rejection applied). In all these cases the CST was smoothed by convolving it with a 25 ms long Gaussian window. The 99% confidence limit of the coherence function was obtained as (17): 1 − (0.01) 1/(L−1) (14) where L is the number of disjoint 1-s segments used in the spectral estimation. Finally, we computed the relative power H1/(H1+B) of the first tremor harmonic with respect to the basic tremor frequency in the estimated tremor EEG component,ŝ 1 (n).
Due to their non-Gaussian distribution (Lilliefors test, p > 0.05), the non-parametric Kruskal-Wallis test was used to compare the differences between the ET and PD patient groups, whereas the Wilcoxon signed rank test was used for paired comparisons. The significance level was set to p < 0.05 and p < 0.01, respectively (see Results section for details).   A,B) show the NMSE and the cross-correlation (CC) between the estimated and simulated tremor source, and the estimated H1/(H1+B) ratio as a function of parameters SNR and σ dr for simulated H1/(H1+B) ratios of 0.0 and 0.286, respectively. Results are averaged over ten simulation runs and depicted as a mean (colored surfaces) and mean + SD (black lines). Note that SD was negligibly small. Figure 2A shows a representative example of detection of the tremor component in the simulated signals, demonstrating that the proposed method accurately extracts the tremor-related component from the synthetic signals. Figures 2B,C show the correlation between the estimated and simulated sources, the NMSE and the estimated H1/(H1+B) ratio as a function of extension factor F, for two different SNR and σ d r values. Note the high correlation and small NMSE between the simulated and the estimated tremor-related component, and the accuracy with which the ratio H1/(H1+B) was detected when extension factor was set to F = 3 or higher. For a SNR of 20 dB, extension factors from F = 2 to F = 5 were optimal, whereas for lower SNRs, extension factors from F = 4 to F = 9 were optimal ( Figure 2C). In both cases the estimated H1/(H1+B) ratio was largely independent from the extension factor in the interval F = 3 to F = 10. All metrics degraded slightly when the model looked too many samples back in the past (F > 10). Based on these results and on the coherence analysis of experimental data (Figure 6), we selected an extension factor F = 8 for further analyses. Note that Figure 6 indicates that our results would hold across a broad range of values of F, from F = 5 to F = 10. Figures 3A,B summarize the NMSE, CC and H1/(H1+B) ratio as functions of the SNR and motor unit firing time variability σ d r at two different simulated H1/(H1+B) ratios. In both cases, the NMSE decreased and the CC increased as the SNR increased, whereas they did not change significantly with σ d j , which suggests that the simulated intrinsic variability in motor unit firing did not affect the source estimation. The estimated H1/(H1+B) ratio did not change significantly with the SNR or σ d r and was always very close to the simulated H1/(H1+B) ratio. Namely, when averaged over all simulated SNRs, σ d r and H1/(H1+B) ratios, the difference between the simulated and estimated H1/(H1+B) ratio was 0.01 ± 0.05.

Simulated Data Analysis
The delay between the estimated and simulated sources,ŝ j (n) and s j (n) was largely independent of the SNR, σ d r , and the simulated H1/(H1+B) ratios, averaging 0.4 ± 1.4 ms across all combinations of parameters. When a 10 ms corticospinal delay was imposed between the motor unit firings and their cortical drive s j (n), the estimated delay averaged 11.0 ± 1.6 ms. This implies a 1.0 ± 1.6 ms difference with the simulated 10 ms delay. Despite this estimate being quite accurate, we want to note that the current simulations do not generate signals as complex as the recorded EEGs. Nor did the simulations incorporate the delays due to propagation of the motor unit action potentials along the muscle fibers or due to EMG decomposition with CKC (25)(26)(27). All these factors contribute to the unknown global delay between 0 and ∼15 ms and cannot be easily estimated in experimental conditions. Thus, we do not expect the experimental estimates of corticospinal delay to be as accurate as those obtained based on the model.   Figure 4 shows an example of EMG decomposition in a representative ET patient, along with the smoothed CST. Table 1 summarizes the number of motor units that were detected from the surface EMG with PNR > 30 dB and then used for the identification of the tremor EEG component. On average, 7.7 ± 5.2 and 8.6 ± 6.3 motor units per contraction were identified for the ET and PD patients, respectively. The average number of firings per motor unit was 160 ± 100 for the ET and 218 ± 130 for the PD patients ( Table 1). In each contraction, all the accurately identified motor units per muscle were used to estimate the tremor EEG component (see Appendix). Since we recorded EMGs from the wrist extensors and flexors of the right arm, two estimates of tremor EEG components were extracted per each task repetition. Figure 5 shows a representative example of the estimated tremor EEG component. Figure 5A depicts the estimated tremor EEG component and how it relates to the smoothed CST, both in the time (left plots) and the frequency domain (right plot). The estimated EEG component exhibits clear tremorrelated activity with peaks both at the basic frequency and the first harmonic of that observed in the pooled motor unit firings and the wrist kinematics. Figure 5B shows time-frequency domain contour plots of the extracted EEG component, the smoothed CST and the wrist kinematics, as reference. In the presented case, the tremor-related EEG activity preceded the pooled motor unit activity and the observed wrist movements, which manifested simultaneously, but this was not always the case.

Experimental Recordings
We investigated how the number of EEG samples in the time domain (extension factor F) influenced the accuracy of the estimation of the tremor-related cortical activity. To this end, we computed the coherence between the tremor EEG component and the smoothed CST for increasing values of F (from F = 1 to F = 15). As shown in Figure 6A, the coherence first increased, but it saturated around extension factor F = 8, in agreement with the simulations with lower SNR ratios. Note, however, that the increase in coherence as more EEG samples were included was not significant after F = 5. Figure 6B demonstrates that the proposed method significantly outperforms the classical coherence between Laplacian-filtered EEG and spatially averaged rectified EMG, low pass filtered at 15 Hz. Figure 7 shows an example of how the proposed method performs compared to the traditional approach of computing the coherence between an estimate of muscle activity (in this case the smoothed CST) and the spatially filtered EEG signals.
We calculated the coherence function for the entire recordings (1 s windows with 50% overlapping), and found no significant coherence values in any EEG channel. In contrast, the coherence between the CST and the EEG component extracted using the proposed method (extension factor: F = 8) was significant. These findings generalized to the all the tested ET and PD patients ( Figure 6B). Our proposed method thus outperforms classic coherence approaches.
In 52% of cases studied (28 of 54), the tremor-related EEG component preceded the CST by 11.0 ± 6.4 ms, whereas in the remaining 48% of cases, the CST preceded the extracted EEG component by 11.0 ± 5.9 ms. All the delays were clustered on the interval between −30 ms and + 30 ms, and we observed no significant difference between the delays in PD and ET patients (P > 0.05, Kruskal-Wallis test). These latency values are in agreement with previous studies (7)(8)(9)(10)13), notwithstanding the limitation listed in the Discussion section.
Finally, a significant difference was observed in the H1/(B+H1) ratio of the extracted tremor EEG between PD and ET patients (Figure 8), also in agreement with previous studies (9,13). Coherence between the CST of the right wrist extensor and the estimated tremor EEG component as a function of the extension factor F. For reference, the average coherence with the EEG signals over C1 is also depicted. The results are averaged over all the trials from the PD and ET patients and reported as mean (bars) and SD (whiskers). (B) Coherence between the Laplacian-filtered EEG and the CST in comparison with the coherence between the extracted tremor EEG component and the CST (our proposed method, F = 8): Two results for CST-EEG coherence are depicted in each graph, one for coherence between CST and spatially averaged EEG ("average EEG") and one for coherence between CST and EEG channel at C1 position ("C1"). Similarly, two results for coherence between CST and extracted tremor component are depicted, one for tremor component, extracted from EEG by using the CST-based estimation described in the Appendix ("tremor CST") and one for tremor component, extracted from EEG by replacing the CST in Appendix by the spatially averaged rectified EMG, low pass filtered at 15 Hz ("tremor rect. EMG"). Results are plotted separately for each investigated muscle. Superscripts *p < 0.05 and **p < 0.01 denote statistically significant difference as assessed by Wilcoxon signed rank test.

DISCUSSION
In this study, we derived and validated a new method for the extraction of the tremor-related EEG activity in the case of pathological tremor. The method builds on the physiological coupling between the tremor-related cortical activity and the neural drive to the muscle (the output of the motoneurons that innervate a muscle). In particular, our method combines the motor unit spike trains identified in the decomposition of high-density surface EMG recordings to build an estimator of the tremor-related EEG component. We applied it to EEG recordings to demonstrate its feasibility, but it could also be used for analyzing magnetoencephalographic (MEG) data.
The proposed method was tested on simulated data and on recordings from 9 PD and 9 ET patients. In the simulations, our method detected the simulated tremor component with great accuracy, as indicated by the low NMSE and high crosscorrelation values. The small difference between the simulated and estimated H1/(B+H1) ratio (Figures 2, 3, global average error of 0.006 ± 0.053 for simulated H1/B, ranging from 0 to 1) further demonstrates the fidelity of the estimated tremor component. Our method also yielded very accurate estimates of the delay between the motor unit population activity and the simulated EEG (the average error was 1.0 ± 1.6 ms for a simulated delay of 10 ms).
In the experimental data, the extracted tremor EEG component exhibited clear similarities with the recorded kinematics and motor unit population activity, both in the time and frequency domains. The ground truth about the estimated EEG tremor component is not available in experimental conditions. However, we believe our method performed well because the estimated EEG component exhibited significantly larger coherence with the identified population of motor units than the spatially filtered EEG signals, which is the standard approach (1-10, 13, 16). This observation indicates that the proposed method is likely to help studying the neural mechanisms of tremor. Indeed, our method always identified tremor-related activity in the EEG, while in many of the investigated cases (34 of 54) we did not find significant coherence between the spatially filtered EEG signals and the identified population of motor units. Note that this observation is in agreement with reports that several tens of second long recordings are needed to obtain robust results in standard coherence analysis (2,(6)(7)(8)(9)(10)13), whereas our datasets were only 30 s long.
Movement artifacts are an important potential confound when studying corticospinal coupling using coherence techniques. We performed two complementary analyses to discard the presence of movement artifacts. First, we tested whether the tremor components were present across many spatially filtered EEG channels, as it would be the case if they resulted from movement artifacts. We calculated the coherence between the CST and each spatially filtered EEG channel. As reported above, in 34 out of 54 cases we did not find significant coherence between the spatially filtered EEG signals and the identified population of motor units. In the remaining 20 cases, significant coherence at the tremor frequency or at its higher harmonics was observed on one or two EEG channels only. As a second control, we examined whether the EEG-CST delays depended on the basic tremor frequency. Finding a significant association between these two parameters would indicate a potential mechanical coupling. Our results ruled out this possibility: the EEG-CST delays lied within the −30 to 30 ms interval. These values did not overlap with the range of delays potentially indicating an artifact (from ∼45 to 100 ms; interval defined by the maximum and minimum tremor frequencies, 11 and 5 Hz, respectively). Therefore, our control analyses indicate that the identified tremor component is unlikely to originate from movement artifacts. Note that these extensive tests are necessary every time the presented methodology is used as, similar to the classical coherence analysis, the movement artifacts could completely mask any tremor-related activity in cortex.
The only parameter that needs to be chosen in our method is the extension factor F in Equation (4). In simulated conditions, the optimal extension factor was dependent on the SNR (of the input sources in the simulated EEG signals) with the optimal values between F = 2 and F = 5 for SNR of 20 dB, whereas larger values, between F = 4 and F = 9, proved to be optimal in the case of lower SNRs (Figure 2). We decided to choose extension factor F = 8 for subsequent analyses. This choice was confirmed during the analysis of the experimental signals, because the coherence between the estimated tremor-related EEG component and the smoothed CST reached the plateau region at F = 8, whereas the overall increase in coherence was not significant after F = 5 ( Figure 6A). The observation that in all the cases studied F = 1 yielded significantly worse results (lower coherence), indicates that the extension of the convolutive model (1) helps in coping with either the existence of different delays in the representation of brain rhythms across different EEG channels, or with the convolutive nature of EEG mixtures. Since the computational complexity of the proposed method increases with the square of the extension factor F, it is to our advantage that the preferred value of F is relatively small.
Regarding the neurophysiological results of this study, we found that the relative power of the first tremor harmonic compared to the basic tremor frequency is greater in PD than ET patients, regardless of the investigated muscle (Figure 8). This is in agreement with other studies using EEG-EMG coherence (9,13). The observed delay between the estimated tremor EEG component and the pooled motor unit firings also agrees with previously reported values. Several studies in ET and PD patients reported a bidirectional interaction between the primary sensorimotor area of cortex and the affected muscles, with an efferent and afferent delay between 10 and 30 ms (7-10, 13). In our dataset, the EEG activity preceded the motor unit firings in half of the cases, and in the other half followed it. This is likely due to the fact that the primary motor and sensory cortices are next to each other and the limited spatial resolution of the EEG makes their activities hard to disentangle. We want to emphasize that these results were obtained using significantly shorter datasets (30 s vs. the typically ≥60 s long signals employed in other studies), and avoiding the need of manually discarding epochs with artifacts.
Results in Figure 6B suggest that the strength of our method derives from the direct use of the CST in the identification of the tremor-related cortical activity. One of the reasons for this is that the CST provide a more accurate representation of the common synaptic input to the muscles than rectified EMG as it eliminates the influence of frequency components introduced by the motor unit action potentials (22,44). In the case of tremor, the CST has most of its power at the frequency of the tremor and its harmonics (15,45). Thus, our approach averages out artifacts and other non-physiological factors (42).
Our corticospinal latency results are consistent with previous studies (7)(8)(9)(10)13). However, they must be interpreted with caution. The convolutive mixing models used to represent the EMG and EEG recordings, which are critical for accurate source separation (25,27,42), may introduce a temporal uncertainty to the reconstructed spike trains and tremor-related EEG components. We estimate this uncertainty to be about ±5 ms for each reconstructed source. Moreover, the propagation of the motor unit action potentials along the muscle fibers from the innervation zone to the uptake electrodes may introduce additional few ms delay. This could potentially further decrease the accuracy of EEG-CST delay estimation. In the current study, we used arrays of several tens of surface electrodes, whereas many previous studies were based on bipolar EMG recordings. The propagation of the motor unit action potentials may differ substantially across these two setups, and may also be muscle specific. Thus, the delays estimated in our study cannot easily be compared to the ones in other studies.
The availability of our method to automatically assess the accuracy with which each motor unit spike train is identified is also of critical importance because this accuracy is then reflected in the extracted tremor-related EEG component (see Appendix). Our group demonstrated in Holobar et al. (42) that motor units with PNR > 30 dB exhibit accuracy > 90% in identification of their firing patterns and in this study, we carefully utilized this knowledge to increase the accuracy of EEG component identification. In the future, we will investigate the minimal number of EEG channels required for accurate detection of tremor-related EEG activity, since it is likely that not all the EEG channels included in this study contribute significantly to tremor identification.
The proposed method is also computationally efficient. The most time consuming step is its first stage (surface EMG decomposition), which typically requires a few minutes of processing time on regular PC for 30 s long measurements. EEG decomposition in Equations (10, 11) is performed quickly.
The method does require multichannel EMG recordings from a muscle, increasing the experimental costs. However, multichannel EMG acquisition demonstrated significant progress in the recent years and became an important source of information in neurophysiology, neurology, sport sciences, prosthetics and ergonomics, to name just a few major scientific fields. Thus, it is likely that the price of multichannel acquisition systems will decrease in the near future.
We limited our study to the EEG decomposition of pathological tremor. The latter is a specific neurological disorder that is characterized by clear spectral peaks in acquired EEG, EMG, and inertial data. It is currently unclear to what extent the presented methodology is applicable to investigations of other types of pathological tremor (e.g., dystonic or cerebellar tremor) or to other disorders, such as multiple sclerosis, stroke and traumatic brain injuries and overactive thyroid, especially as tremor frequently accompanies these disorders. All these questions need to be systematically addressed in separate studies.
In conclusion, we have presented a novel method for estimating tremor-related cortical activity. This method uses pooled motor unit firings to directly extract the tremor component from cortical recordings. Based on the presented results, we believe that our method constitutes a significant step forward in the current state-of-the-art as: (a) it is the first method that directly extracts the tremor component from EEG recordings; (b) it successfully tracks time changes in the tremorrelated cortical activity and has a potential for online tremor detection.

AUTHOR CONTRIBUTIONS
AH, JG, ER, JR, JB-L, JP, and VG participated in conceptualization and data acquisition. AH, JK, and VG participated in formal analysis and methodology. AH, ER, JB-L and JP participated in project administration and resources. All authors contributed to writing the original draft and revised the manuscript. All authors read and approved the final version of the manuscript and agreed for all aspects of the work.