Abstract
Pathological tremor is a common but highly complex movement disorder, affecting ~5% of population older than 65 years. Different methodologies have been proposed for its quantification. Nevertheless, the discrimination between Parkinson’s disease tremor and essential tremor remains a daunting clinical challenge, greatly impacting patient treatment and basic research. Here, we propose and compare several movement-based and electromyography-based tremor quantification metrics. For the latter, we identified individual motor unit discharge patterns from high-density surface electromyograms and characterized the neural drive to a single muscle and how it relates to other affected muscles in 27 Parkinson’s disease and 27 essential tremor patients. We also computed several metrics from the literature. The most discriminative metrics were the symmetry of the neural drive to muscles, motor unit synchronization, and the mean log power of the tremor harmonics in movement recordings. Noteworthily, the first two most discriminative metrics were proposed in this study. We then used decision tree modelling to find the most discriminative combinations of individual metrics, which increased the accuracy of tremor type discrimination to 94%. In summary, the proposed neural drive-based metrics were the most accurate at discriminating and characterizing the two most common pathological tremor types.
1. Introduction
Pathological tremor is the most common movement disorder and is strongly increasing in incidence and prevalence with ageing. About ~4% of the elderly (>65 years) suffer from essential tremor (ET), whereas ~1% of the population of age >50 years suffer from Parkinson’s disease (PD) [1].
The central origin of ET and PD tremor is widely accepted but not fully understood. ET is believed to originate at the cerebellothalamocortical pathways [2, 3], although the involvement of the inferior olive has also been suggested due to its rhythmic properties [4, 5]. In PD patients, the loss of dopaminergic nigrostriatal neurons is believed to cause abnormal oscillations in the loop linking the cortex, basal ganglia, and thalamus [6]. Despite their potentially different origins, ET and PD frequently demonstrate similar symptoms and are often difficult to discriminate, especially at early stages. Up to 20% of patients with ET develop the signs typical for PD, whereas 30 to 50% of patients diagnosed with ET do not have ET [2, 7, 8].
Current methods for clinical tremor diagnosis and quantification are based on the mechanical demonstration of tremor (the tremulous movement), mainly with the help of movement disorder scales, such as the Fahn-Tolosa-Marin [9] and the UPDRSIII scale [10]. Some of the key clinical features for diagnosing ET and PD are that ET tends to increase in postural and kinetic tasks, whereas PD tremor is frequently present at rest. However, 20–30% of the ET patients also have clinically noticeable rest tremor [11], an observation that highlights the complexity of diagnosing these two tremor types. As to the tremors themselves, the mechanisms behind the observed differences (and similarities) between ET and PD are not well understood. In addition, tremor diagnosis is mainly based on the subjective clinical assessment of the tremulous movement (i.e., without quantifying its characteristics), whereas the pathophysiology of the underlying tremor is not fully understood, mainly due to the highly complex interactions between the central and peripheral structures of the nervous system [12].
Different computer-aided methods for tremor diagnosis have already been proposed [13–18]. Among them, the difference in the energy of the higher tremor harmonics in movement recordings, typically performed with inertial sensors, best discriminates ET and PD patients [17, 18]. These harmonics are usually computed after using relatively simple frequency filtering techniques which separate the higher tremor frequencies from the relatively low frequencies of voluntary movement [14, 15, 17, 19]. However, the tremulous movement is nonlinearly related to the underlying muscle activation and results from complicated interactions among the many muscles that control the same joint. These interactions between muscles are also yet not fully understood. For example, it has been recently demonstrated that not all the muscles are equally affected by ET [20] and that the interaction between antagonist muscles is different across conditions and patients, and intrinsically nonstationary [21].
In this study, we applied a sophisticated blind source separation algorithm [22–24] to the surface electromyograms (EMG) of PD and ET patients to estimate the neural drive that reflects the net input from the nervous system [25] to the wrist flexors and extensors. Our rationale was that by using a more detailed measurement of nervous system function than the classically used rectified EMG or resultant movement, we could develop better metrics for tremor characterization and discrimination. We systematically analyzed the properties of the neural drive in these two groups of patients, proposed novel pathological tremor characterization metrics, and analyzed their power to discriminate between ET and PD tremor. Finally, we compared the discriminative power of the proposed methodologies to that of previously proposed methods and demonstrated that our new metrics largely outperform them.
The manuscript is organized as follows. Section 2 justifies multichannel EMG recordings and the estimation of the neural drive to muscle using blind source separation of these signals in pathological tremor. Section 3 describes the experimental protocol and our novel methodology for quantifying and classifying pathological tremor types. Section 4 systematically analyzes the statistical power of the described metrics in discrimination of ET and PD patients. Section 5 discusses the main results and limitations and concludes the manuscript.
2. Neural Drive to Skeletal Muscles and Tremor Demonstration in Surface Electromyograms
The central nervous system causes movement by modulating the firing patterns of the motoneurons that innervate task-relevant muscles. The net output from all the motoneurons is referred to as neural drive to the muscle [26] and is the ultimate control signal mediating movement generation [27]. In heathy muscles, each motoneuron firing causes a motor unit action potential (MUAP) that makes the innervated fibers contract. MUAPs can be measured by surface EMG [28]. Subcutaneous tissue progressively filters electrical potentials with increasing distance from the muscle fibers, inducing spatial variability of MUAPs across the skin surface [28].
In pathological tremor, the neural drive to muscle has classically been estimated from bipolar surface EMG recordings [14, 15, 17]. In this procedure, the EMG signals are first rectified and their power spectrum is calculated to estimate the power at the main tremor frequency and its higher harmonics. While appropriate for coarse tremor assessment and for many pioneering tremor studies (e.g., [29]), this procedure of computing the power spectrum of the rectified EMG did not provide the most accurate insight into the complex and nonstationary neural drive underlying tremulous muscle dynamics in our study (see Results in Section 4). There are a number of reasons that support this statement as well as our results. First, not all the motor units are active during pathological tremor, and the muscle activation pattern, estimated from bipolar EMG recordings, depends significantly on the electrode position (Figure 1) and on motor unit distribution within the muscle tissue [27]. Second, in dynamic conditions as tremor, muscle fibers move relative to the skin, causing unpredictable changes in the detected MUAPs [25, 27, 28]. Third, the spatial distribution of the motor units varies significantly among subjects, and there is no universally optimal position for bipolar EMG recordings [28]. Fourth, the algebraic summation of biphasic MUAP trains causes EMG amplitude cancellation, which leads to an unpredictable underestimation of the effective neural drive [30–32]. For these reasons, the neural drive to muscle is better characterized using two-dimensional high-density (2D hdEMG) arrays of surface EMG electrodes, with several tens of surface electrodes placed over a substantial portion of the skin above the investigated muscle [25, 27]. This, however, generates large quantities of recorded data.
(a)
(b)
Different methodologies have been proposed to cope with this challenge [28]. Given that the neural drive to muscle is the net activity of all the innervating motoneurons, the best estimate of the neural drive is derived from the motoneuron (and, consequently, motor unit) spike trains themselves [25]. For this purpose, the highly variable MUAPs that mask the changes in the neural drive need to be separated from the motor unit spike trains. This is accomplished by using advanced signal decomposition algorithms [23, 24, 32–34]. However, not all of these algorithms are applicable to hdEMG data in pathological tremor; in pathological tremor, motor unit firings tend to be significantly more synchronized than in healthy conditions [21, 23], causing bursts of EMG activity (Figure 2) that are very difficult to decompose into the contributions of individual motor units. The Convolution Kernel Compensation (CKC) technique has been recently demonstrated to successfully solve this problem [23], providing a solid ground for fully automatic and highly accurate assessment of the tremulous neural drive to muscles.
3. Methodology
3.1. Experimental Protocol
Fifty-four subjects (27 ET and 27 PD patients), all suffering from mild to severe tremor, were recruited at Hospital Universitario 12 de Octubre, Madrid, Spain. The ET patients were 12 females and 15 males with mean age of years (range 43–80 years). The PD patients were 10 females and 17 males with mean age of years (range 44–80 years). Tremor severity in the most affected limb ranged from 10 to 56 (mean ± SD: ) according to the Fahn-Tolosa-Marin scale in ET patients and from 5 to 51 (mean ± SD: ) according to the UPDRS scale in PD patients. The mean disease duration was years (range 2–65 years) for the ET patients and years (range 1–13 years) for the PD patients. Seventeen ET patients were taking antitremor drugs, which in all cases were withheld for at least 12 h before the recordings. Sixteen PD patients were taking antitremor drugs and continued their medications during the recordings. This is the group of patients that was included into the induction of the decision model and its tenfold cross-validation.
Based on the method we devised to objectively identify tremor, 27 ET and 27 PD patients exhibited tremor in both rest and postural tasks studied, and 27 ET and 22 PD patients exhibited bilateral tremor. Our tremor identification method consisted in computing the power spectrum of the signal of interest (the summed motor unit spike trains, the rectified surface EMG, or wrist movement) and summing the peaks at the basic tremor frequency and its first two higher harmonics. If this sum was ≥50% of the total power of the signal, we considered that recording to exhibit tremor.
To further test the performance of our decision tree model, we recruited six additional ET and PD patients, which were blindly diagnosed by the proposed model and three expert clinicians. Importantly, both diagnoses were completely independent, and no data from these six patients were included into decision model induction. The ET patients were 1 female and 2 males, with mean age years (range 47–80 years); their mean disease duration was years (range 10–27 years). The PD patients were 2 females and 1 male, with mean age years (range 62–80 years); their mean disease duration was years (range 5–14 years). Tremor severity in the most affected limb ranged from 10 to 34 (mean ± SD: ) according to the Fahn-Tolosa-Marin scale for the ET patients and from 15 to 23 () according to the UPDRS scale for the PD patients. These patients formed the test group.
All patients received a detailed explanation of the study and gave written informed consent prior to participation. The study was conducted in accordance with the Declaration of Helsinki and approved by the local ethics committee. The experimental protocol comprised measuring the inertial and EMG data from both forearms during standard tremor triggering tasks (see, e.g., [2, 6, 11, 14, 16–18, 21]). Each patient performed three repetitions of the following 30 s long standard tasks:(i)Rest task (RE): the patient was sitting relaxed, with both arms supported in elbow and wrist region. The patient’s wrists were relaxed with no noticeable voluntary activity throughout the entire duration of the task.(ii)Arms outstretched against gravity (AO): after 5 seconds of rest, the patient outstretched his/her arms and maintained them parallel to the ground until the end of the task. The fingers were separated and the hands were held at the same height as the shoulders.(iii)Arms outstretched against gravity with weight (WE): it is the same as the AO task with additional weight load of 1 kg applied to both hands. The patient had the fingers outstretched.(iv)Arms supported + postural tremor elicited (PO): it is the same as RE task, but with both wrists held extended against gravity.Wrist movement was assessed as the difference between the measured acceleration values in the axis parallel to the forearm [19]. Accelerometers (Tech MCS, Technaid S.L., Madrid, Spain) were fixed with surgical tape over the hand dorsum and the distal third of the forearm (on the dorsal side, close to the wrist) on each hand. Data were sampled at 100 Hz by a 12-bit A/D converter and low-pass filter (cut-off frequency of 20 Hz).
Surface hdEMG signals were recorded from the wrist flexors and extensors with 5 × 13 electrode arrays (LISiN–OT Bioelettronica, Torino, Italy, 8 mm interelectrode distance). The electrode arrays were centred over flexor carpi radialis and extensor digitorum communis, respectively. Before electrode placement, the skin was lightly abraded using abrasive paste (Meditec–Every, Parma, Italy) and cleaned. In order to improve the electrode-skin contact, each of the electrode cavities in the array was filled 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 2048 Hz and 12-bit resolution (LISiN–OT Bioelettronica, Torino, Italy). EMG and movement recordings were resampled to 2048 Hz and synchronized offline, using a common clock signal.
The hdEMG signals were decomposed into motor unit spike trains using the CKC algorithm [24]. The pulse-to-noise ratio (PNR) metric [23] was used to assess motor unit identification accuracy and all the motor units with PNR < 30 dB (corresponding to accuracy < ~90% [24]) were discarded.
3.2. Assessment of the Tremor in the Neural Drive to Muscles
As mentioned above, tremor results from the pathological activation of several synergistic muscles, and oftentimes, of antagonists as well. Due to the central origin of tremor, the descending drives to nearby muscles will possibly exhibit a similar tremor component [35], or at least its delayed reflection via afferent feedback loops [20, 21]. However, the characteristics of the neural drive to different muscles have been rarely studied in pathological tremor [23].
In this study, the neural drive to muscle was assessed by summing up individual motor unit spike trains into cumulative spike train (CST) (see, e.g., [20, 21, 35]). Afterwards, the CSTs of the different muscles were low-pass filtered (cut-off frequency of 20 Hz). In order to emphasize the rhythmic motor unit activity related to tremor and minimize the asynchronous motor unit activity caused by the concurrent voluntary contraction, we calculated the auto- and cross-correlation functions among all pairs of filtered CSTs. Representative results from one ET patient and one PD patient are depicted in Figures 3 and 4. Note that, in Figure 3(a), the neural drive to the extensors of the right wrist exhibited both tremulous and voluntary motor unit activity. The proposed auto- and cross-correlation functions filtered out the voluntary activity and emphasized the tremulous neural drive (Figure 3(b)). Moreover, usage of the auto- and cross-correlation functions of the neural drive also circumvented the ambiguity of our hdEMG decomposition method when estimating the amplitude of the neural drive. This well-known ambiguity is common to all blind source separation algorithms [22, 24] and implies that the amplitude of the neural drive can only be determined up to an unknown scalar factor (Figures 3, 4, and 6). This prevents the assessment of the absolute tremor power in the neural drive; however its relative power (with respect to the voluntary drive) can be easily assessed. The cross-correlation functions assess the interactions between the neural drives to different muscles and are, thus, insensitive to this amplitude ambiguity [25].
(a)
(b)
(a)
(b)
We observed different neural drive patterns for the two groups, with the ET patients exhibiting a larger degree of shared neural drive across muscles than the PD patients. To quantify this phenomenon, we considered the following neural drive-based metrics for ET and PD discrimination:(N1)Tremor amplitude in the neural drive: average difference between the maxima and minima of the autocorrelation function of the neural drive to each individual muscle, calculated over correlation lags from −1 to +1 s.(N2)Tremor frequency in neural drive: dominant frequency in the power spectrum of the neural drive to each individual muscle.(N3)Dominant hand tremor symmetry: absolute values of cross-correlation coefficients between the neural drive to the wrist extensors and the neural drive to the wrist flexors in the dominant tremor hand, averaged over correlation lags from −1 to +1 s (Figure 3(b), panel N3).(N4)Interlimb tremor symmetry in extensors: absolute values of the cross-correlation coefficients between the neural drive to the extensors of the left and the right wrist, averaged over correlation lags from −1 to +1 s (Figure 3(b), panel N4).(N5)Interlimb tremor symmetry in flexors: absolute values of the cross-correlation coefficients between the neural drive to the flexors of the left and the right wrist, averaged over correlation lags from −1 to +1 s (Figure 3(b), panel N5).(N6)Global tremor symmetry: absolute values of the cross-correlation coefficients between the neural drives correlation lags, averaged over all the muscle pairs and over correlation lags from −1 to +1 s (all panels in Figures 3(b) and 4(b)).(N7)Normalized global tremor symmetry: the global tremor symmetry for the postural task (AO, WE, or PO), divided by the global tremor symmetry for the RE task.
3.3. Assessment of the Spatial Distribution of the Neural Drive to a Muscle
CST captures the temporal variability of neural drive to a muscle, while averaging out the differences among identified motor units [36]. Thus, the CST ignores how the tremor input is distributed across the motoneurons innervating the muscle. The distribution of the tremor input across motor units may be studied by quantifying the pairwise motor unit synchronization. During voluntary contractions, motor units discharge asynchronously, supporting smooth modulation of force [28]. In contrast, during pathological tremor, the level of motor unit synchronization increases significantly [20, 21, 23]. This can be quantified by assessing the distribution of backward and forward motor unit recurrence times in pairs of simultaneously active motor units [37]. The forward (backward) recurrence time is defined as the distance from the current motor unit discharge to the first next (closest previous) discharge of another motor unit (Figure 5(a)). In normal conditions, the recurrence times are uniformly distributed over the observed time support, indicating low motor unit synchronization. In pathological tremor, the motor unit recurrence times are grouped together (Figure 5(b)), reflecting highly synchronous motor unit activity.
(a)
(b)
(a)
(b)
We analyzed the distributions of motor unit recurrence times by calculating their histogram. We computed the 99% confidence limit assuming a uniform distribution and then identified the peaks exceeding this limit [37]. We then defined the pairwise motor unit synchronization index as the ratio between the peak area and the sum of all the histogram bins (Figure 5(b)). Finally, we used this index to define the following metrics of spatial tremor variability:(N8)Coefficient of variation for the pairwise motor unit synchronization: the standard deviation of the distribution of pairwise motor unit synchronization indexes divided by their mean value, calculated over all the pairs of identified motor units per muscle and per task.(N9)Normalized coefficient of variation for the pairwise motor unit synchronization: the N8 metric in the AO task, divided by the N8 metric in RE task. The normalized coefficient of variability was averaged over all four investigated muscles for each patient.
3.4. Assessment of Tremor in Wrist Movements
Tremor characterization with inertial sensor recordings is very appealing due to the simplicity, robustness, and cost-effectiveness of these sensors [38]. Inertial sensor data have been used in many tremor diagnosis studies, demonstrating promising results in terms of ET and PD discrimination [13, 16, 18]. However, movement data reveal very limited information on the underlying tremor mechanisms. It then follows that the ability to discriminate between complicated ET and PD cases with metrics extracted from movement data is likely to be limited. To verify this hypothesis, we computed, separately for the dominant and nondominant tremor hands in each recorded task, the following metrics [18]: the basic tremor frequency and its power, the number of spectral peaks at the higher harmonics of the basic tremor frequency, and their mean power. We also computed the asymmetry of the autocorrelation function [13], along with the statistical signal characterization (SSC) metric proposed in [16].
We also defined the following new movement-related metrics, calculated using the base 10 logarithm of the harmonic peak powers in the data from the dominant and nondominant tremor hand, respectively:(M1)Mean log power of all harmonics: the mean of the logarithmic power of all the tremor harmonic peaks in the power spectrum of the movement recordings during a single task.(M2)Average of mean log power of all the harmonics: the average of the mean of the logarithmic power of all the harmonic peaks across different tasks.In order to make the differences between ET and PD patients most evident, and to further improve the patient classification performance, the values of all the movement metrics in the postural tasks (AO, WE, and PO) of individual patient were normalized by their corresponding values during the RE task. In this normalization, we subtracted the metric in the RE task from the corresponding values in the AO, WE, and PO tasks, respectively. We chose the repetition with the highest tremor amplitude for each of the AO, WE, and PO tasks, respectively, and the repetition with the lowest tremor amplitude for the RE task. This yielded one value for each type of task (AO, WE, and PO) for each patient.
3.5. Assessing the Discriminative Power of the Tremor Metrics
The presented metrics were combined with metrics describing the patient’s demographics (age, gender) and tremor etiology (tremor onset, tremor duration, tremor severity, results on UPDRS/Fahn-Tolosa scale, family history, medications, and predominance). In total, 194 tremor metrics were analyzed. The discriminatory power of each extracted metric was evaluated using the receiver operating characteristics (ROC) curve. In particular, we used the following metrics: area under the ROC curve (AUC) with its corresponding 95% confidence interval (95% CI); specificity (Sp); sensitivity (Se); positive predictive value (PPV); and negative predictive value (NPV). The ROC curve describes the performance of different models based on varying threshold values, and thus allows selecting the optimal operating point [39]. In our study, the optimal operating point was the best trade-off between sensitivity and specificity. Sensitivity denotes the ability of a prediction model to correctly classify all positive cases and specificity is the ability to correctly classify negative cases. In our study, ET patients were indicated as positive cases and PD patients as negative cases. PPV and NPV values provide additional insight into the robustness of the generated prediction models. PPV or precision is the proportion of true positives (ETs) among all the patients classified as positive. Analogously, NPV or false discovery rate is the proportion of true negatives (PDs) among all the patients classified as negatives. Additionally, a Mann–Whitney test was used to compare the extracted metrics between PD and ET patients. The basic level of statistical significance was set at .
Finally, the aforementioned movement and neural drive metrics were tested against their discriminative power in decision support systems. We were interested in methods that provide a symbolic representation of the extracted knowledge and that implement a simple data mining model. Therefore, we gave preference to binary logistic regression and decision tree modeling algorithms instead of black-box approaches (e.g., neural networks, support vector machines). However, as in most of the statistical modeling methods, the drawback of logistic regression is that all the observations that contain missing values are ignored, resulting in a reduced data set. This can substantially weaken the predictive power of these models. This limitation is mitigated in decision tree algorithms, where surrogate splitting rules enable the use of other predictor variables to perform a split for observations with the missing values. We thus chose four different decision tree models, namely, CART [39], CHAID [40], C5.0 [41], and QUEST [42]. These four different algorithms differ basically in the splitting criteria. The CART (Classification and Regression Tree) algorithm uses the Gini index as splitting criterion [39]. The CHAID (Chi-Square Automatic Interaction Detector) algorithm uses chi-square tests of independence for evaluating an association between categorical variables [40]. The C5.0 algorithm, which is an improved version of the original Quinlan’s C4.5 algorithm [41], defines each split based on the metric that provides the maximum information gain. In the QUEST (Quick, Unbiased, Efficient Statistical Tree) algorithm, the splits are determined by quadratic discriminant analysis using the selected input features on groups formed by the target categories [42].
The decision tree models were evaluated using tenfold cross-validation (see details in [43]). In this method, the original dataset is randomly partitioned into ten subsets. Then, in each of ten runs, a single subset is retained as testing set, while the remaining nine subsets are used as training data for computing the prediction model. Using this method, all the observations are used for both training and testing; each observation is used once for testing. The results from all ten models are then averaged to produce a single estimation of model performance.
4. Results
In total, ~20,500 motor units were identified with an accuracy above 90% according to the PNR metric [21]. The number of identified motor units did not depend significantly on the type of tremor triggering task (Kruskal-Wallis test, ). A relatively large number of motor units was also identified in the rest condition in ET patients and during the postural tasks in PD patients. Figure 6 depicts the average rest and postural tremor amplitude in the 27 ET and 27 PD patients as derived from the neural drive, the spatially averaged rectified EMG and accelerometer data. Despite the statistically significant differences in tremor amplitude between ET and PD patients, none of these differences accurately discriminated ET and PD patients. This observation demonstrates the complexity of ET and PD tremor [2].
First we assessed the discriminative power of each of the metrics described above. Due to the large number of tested metrics, we only present the results for the tasks and metrics with highest discriminatory power (Tables 1 and 2). Tremor amplitude (N1) and tremor frequency (N2) in the neural drive did not significantly differ across ET and PD patients (Table 1, Figure 6). However, tremor symmetry in the dominant tremor hand (N3) was significantly different between ET and PD patients, especially during the AO and WE tasks. These differences further increased when we normalized the dominant tremor hand symmetry for the AO task with its corresponding value for the RE task (Table 1). The neural drive metric that best discriminated across tremor types was the normalized global tremor symmetry (N7) in the AO task (AUC = 0.94, Se = 85.2, Sp = 96.2), followed by the same tremor symmetry metric during the WE task (AUC = 0.90, Se = 85.2, Sp = 100.0). The second most discriminative metric was the normalized coefficient of variation of the pairwise motor unit synchronization (N9) during the AO task (Table 1). Notably, neural drive metric N7 outperformed all movement-related metrics at tremor discrimination (Table 2), as expected, because it captures more subtle aspects of the neural mechanisms of each condition. Figure 7 shows the N6, N7, N8, and N9 metrics for the AO and RE tasks. Note that the normalization of metrics during postural tremor triggering tasks with their values during the RE task significantly increased the discriminative power.
(a)
(b)
(c)
To test the added value of analyzing the neural drive via cumulative spike trains rather than EMG data, we compared the discriminative power of the normalized global tremor symmetry (N7) metric in the AO task as assessed from (1) the neural drive (low-pass filtered cumulative MU spike train); (2) a single rectified EMG channel (the channel with the maximal EMG energy per muscle); and (3) the spatial average of all the rectified EMG channels in the investigated muscle. We employed the same exact methodology in all three cases. The symmetry index computed from the channel with maximum EMG energy misclassified 8 ET and 10 PD patients (Se = 70.4, Sp = 63.0, PPV = 65.5, NPV = 68.0). The symmetry index obtained from all rectified EMG channels performed slightly better, but still misclassified 7 ET and 7 PD patients (Se = 74.1, Sp = 74.1, PPV = 74.1, NPV = 74.1). In clear contrast, the symmetry index calculated from the neural drive only misclassified 4 ET and 2 PD patients. These results support our statement that characterizing the neural drive to muscle based on the cumulative motor unit spike train outperforms other tested measures, likely because it provides the most precise characterization of the effective output of the innervating motoneurons.
Among the movement metrics, the highest AUC (AUC = 0.83, Se = 88.9, Sp = 70.4) was obtained with the M2 metric (mean log power of all the harmonics) of the nondominant tremor hand, averaged over the AO, WE, and PO tasks and normalized by the RE task (Table 2). The metric M1 during the AO (for the nondominant tremor hand) and WE tasks (for the dominant tremor hand) and normalized on RE task also achieved good performance (Table 2). Overall, the discriminative power of the movement metrics was slightly lower in the dominant than in nondominant tremor hand.
After testing the discriminative power of individual metrics, we compared the performance of the different decision tree algorithms. Among the tested decision tree algorithms (Section 3.5), the best predictive power was obtained using the CHAID algorithm (Figure 8). The overall accuracy of the induced model on the whole dataset was 94.4%. All the PD patients were correctly classified and only 3 ET patients were misclassified (AUC = 0.99 (0.98, 1.00), Se = 88.5, Sp = 100.00, PPV = 100.0, NPV = 90.0). The tenfold cross-validated risk estimate for the final tree was calculated as the average of the risks for all of the ten trees (Risk = 0.189, Std. Error = 0.054).
Only three metrics were automatically selected by the final CHAID model:(i)The normalized global tremor symmetry for the AO task normalized by the same metric for the RE task (N7)(ii)The normalized coefficient of variation for the motor unit synchronization index (N9) during the AO task in dominant tremor hand(iii)The mean log power of all the movement harmonics, averaged over the AO, PO, WE, and RE tasks, for the nondominant tremor hand (M2).Notably, these three metrics were proposed in this study, which highlights the potential impact of these novel methods. The most discriminative metric, the neural drive symmetry for the AO task normalized by the RE task (N7), was placed in the root of the decision tree. The cut-off point = 1.555 clearly divided 25 PD patients, with parameter value lower or equal to , from 23 ET patients (Figures 7 and 8). Thus, PD patients exhibited higher neural drive symmetry in the RE task than in the AO task, whereas the opposite was true for the ET patients (note that this result could not have been drawn from the rectified surface EMG). The other two metrics were used in the second level of the decision tree, thereby improving classification accuracy to accurately diagnose 2 additional PD patients. The ET patients exhibited lower variability in motor unit synchronization for the RE task than the PD patients (Figure 7). The interpretation of M2 splitting rules is more difficult, mainly due to the complex relationship between muscle activation dynamics and movement discussed above.
To further test the performance of our decision tree system for tremor diagnosis, we applied it on the data from the six additional patients (test group) that had not been included in the training datasets. Remarkably, all these six additional patients received the model-based diagnosis that matched the clinical diagnosis. This result highlights the potential of our model to generalize to new (unprocessed) patients.
5. Discussion
In this study we tested 194 tremor metrics (many of them proposed herein) with the goal of developing a computer-based system to diagnose ET and PD tremor. Despite this large number of metrics, the most successful data mining model relied only on three metrics we proposed: the global tremor symmetry in the AO task normalized to the RE task (N7); the coefficient of variation for the pairwise motor unit synchronization (N9) in the AO task normalized to the RE task; and the mean log power of all the harmonics averaged over the AO, PO, WE, and RE tasks (M2) in the nondominant tremor hand. The same three tremor metrics were also proposed by logistic regression (results not shown). Critically, the two metrics that best discriminated tremor types were based on the analysis of motor unit activity. This indicates that accurate computer-aided tremor diagnosis benefits from examining fine aspects of neural control that are not detectable with raw surface EMG recordings.
There are several possible explanations why the neural drive metrics (especially tremor symmetry) are more informative than the wrist movement metrics. First of all, the tremulous neural drive may be subtle and obscured by substantial simultaneous voluntary and/or involuntary cocontraction of several antagonistic muscles, which may stabilize the controlled joint and, thus, mask the tremor component of movement. Second, due to the intrinsic complexity of the tremor dynamics [21], the tremor may shift from one muscle to another or may affect the whole group of muscles during the entire trial. When those muscles control the same joint, the changes of their net mechanical effect can be relatively small or relatively complex and, thus, difficult to interpret and quantify from movement measurements. Third, the interactions among synergistic and antagonistic muscles likely affect the harmonics of the wrist movement, but this effect is highly complex and difficult to interpret from wrist movements. Assessments of the neural drive to individual muscles provide more precise insight into the underlying neural mechanisms of tremor, which likely explains why it is the basis for the most discriminative metrics.
In pathological tremor, the motor unit firing patterns that ultimately cause the tremulous movement depend on the strength of the descending common tremor drive to the motoneuron pool; the concomitant voluntary drive, if present; afferent feedback from receptors within the same muscle and from other muscles; the intrinsic motoneuron properties; and the level of synaptic noise [20, 21]. By summing motor unit spike trains together into CST, the intrinsic noise components of many motoneurons average out, revealing their time-varying common input [21, 36, 44, 45]. Autocorrelation function of neural drive further emphasizes the rhythmic (synchronous) motor neuron activities and suppresses the asynchronous ones. When extended to several muscles, this analysis of neural drive may provide also an insight into the tremor drive projections to neighboring muscles [35] or even compartments within the same muscle [46].
Critically, this detailed assessment of the neural drive is unlikely to come from classical EMG measures, such as the EMG power spectrum or the rectified EMG, because, contrary to CSTs, these measures reflect both the tremor-related input to the muscle and the undesired effects of muscle anatomy and geometry on the EMG [25, 28]. The latter aspects are not relevant for the neural control of movement and thus constitute “biological noise” from a data processing point of view. This “noise” is mostly due to the negative influence of the MUAP shapes (see, e.g., [25, 30–32, 47]), caused by the heterogeneity of subcutaneous tissues and the (close to) random motor unit distribution within the muscle tissue [27]. In fact, the EMG power spectrum can be directly expressed as the sum of the MUAP spectra, weighted by motor unit discharge rates [25]. The CSTs, analyzed in this study, are not sensitive to MUAP shapes (nor to muscle geometry and anatomy) and thus constitute much more appropriate means for characterizing the tremulous drive than the rectified EMG. In this regard, here we clearly demonstrated the negative effect that a single surface electrode placement has on the assessment of pathological tremor characteristics (Figure 1). We also showed that metrics based on the rectified EMG have less discriminative power than metrics based on the summed motor unit spike trains. This was true even when we minimized the impact of electrode location by averaging the rectified EMG over all 60 EMG channels in the electrode matrix whereat the electrode matrix covered most of the skin above the muscle of interest (see Section 4).
Possibly due to the aforementioned very high sensitivity of our tremor identification metrics, we found that all the tested ET and PD patients had tremor in the tested rest and postural conditions (Figures 6 and 7). In several cases, the tremor was subclinical and almost invisible to the naked eye. Thus, further studies of this observation are required. However, all the three selected metrics, namely, the global symmetry (N6), the coefficient of variation for pairwise motor unit synchronization metrics (N8), and the mean log power of all the movement harmonics (M2), increased significantly their discriminative power when we normalized their values in the postural tasks by their value in the rest task (Tables 1 and 2, Figure 7). Critically, this indicates that, while tremor properties (and thus the values of our metrics) are similar across PD and ET patients during an individual task (Figure 7), they have very different trends when switching from rest to posture. In our study, both the global tremor symmetry and the coefficient of variation for the pairwise motor unit synchronization increased in the ET patients, whereas the increase was much smaller (or there even was a decrease) in the PD patients. Similar trends were observed in the wrist movement metrics.
An interesting observation from all of the experiments is that none of the demographics/clinical parameters (age, gender, tremor onset, tremor duration, severity of tremor, results on UPDRS/Fahn-Tolosa scale, family history, medications, and predominance) was included in the models built by stepwise algorithms. This suggests that the proposed tremor metrics complement clinical parameters in discriminating ET and PD, thereby demonstrating their potential to aid diagnosis in complicated cases.
We also want to note that one of the three ET patients that were misclassified by the CHAID decision tree model was rediagnosed as PD more than one year after he was originally diagnosed by his neurologist. Another misclassified ET patient was rediagnosed as dystonic tremor with ET. The third misclassified ET patient experienced very mild tremor in all of the investigated tasks. These observations further highlight the potential of the proposed tremor metrics and the proposed decision model to aid clinical diagnosis of tremors. Also, importantly, these objective measures open the door to deriving metrics of patient response to different treatments and to drawing objective comparisons among interventions. Even though clinical diagnosis will remain the gold standard for tremor diagnosis, our data suggest that neural drive-based metrics may identify patient-specific peculiarities that are otherwise hard to detect, especially during the early stage of the disease when misdiagnosis is most common [7, 8]. Furthermore, due to their objectiveness and precision, these metrics may help guide patient follow-up and treatment.
There are two practical limitations that make the clinical use of neural drive metrics less attractive than movement metrics. First, the acquisition of hdEMG signals is more cumbersome and time consuming than movement recordings, which can be easily performed using inertial measurement units. Second, several muscles need to be recorded to compute the proposed metrics. These limitations may be mitigated in the near future by recent developments of unobtrusive and wireless surface EMG sensors, such as epidermal electronics [48]. We further acknowledge that, in this study, we only measured each patient in a single session. Thus, we were not able to assess the test-retest reliability. For the same reason, we could not evaluate the long-term tracking of tremor characteristics. Finally, we focused on the discrimination of ET and PD patients because these are the most prevalent types of tremor [1, 49]. The analyses of other types of tremor such as dystonic tremor are left for future work.
In summary, here we have introduced novel metrics to quantify and discriminate PD and ET tremor. In particular, we proposed metrics that characterize the temporal (CSTs), spatial (the normalized coefficient of variation for the pairwise motor unit synchronization), and cross-muscle tremor variability (neural drive symmetry). We showed that these metrics provide highly accurate insight into the characteristics of ET and PD tremor, and that they outperform metrics based on the rectified EMG and movement data. By demonstrating the discriminative power of the neural drive-based metrics on 30 ET and 30 PD patients, we confirmed our hypothesis that neural drive assessment provides novel insight into these pathologies.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was supported by the Commission of the European Union (Grant agreement no. ICT-2011.5.1-287739 “NeuroTREMOR: A novel concept for support to diagnosis and remote management of tremor”), by the Slovenian Research Agency (research core funding no. P2-0041 and the Project J2-7357 “Exact quantification of muscle control strategies and co-activation patterns in robot-assisted rehabilitation of hemiparetic patients”), and by Spanish government, MINECO (Research Grant DPI2015-72638-EXP). J. A. Gallego received funding from the Commission of the European Union (FP7-PEOPLE-2013-IOF-627384).