Study participants

Eligible PD patients and HC were right-handed individuals between ages 40 and 80 years. Participants in the PD group were diagnosed as idiopathic PD according to the UK Brain Bank Criteria [45]. PD participants at regular follow-up at the Movement Disorders Institute at the Chaim Sheba Medical Center, Ramat Gan, Israel were included. Some were taking part in a repetitive transcranial magnetic stimulation (rTMS) study with brain monitoring; however, the current analysis used only data obtained before rTMS treatment. HC were recruited through advertisement in Sheba Medical Center. Exclusion criteria were: 1) an unstable medical condition, or 2) a history of epilepsy, dementia, cerebrovascular disease, previous head injury, brain tumor, or any craniosurgical intervention; 3) active depressive or psychotic symptoms, current drug abuse or alcoholism, antipsychotic treatment. Eligible participants had to score below 14 on the Beck Depression Inventory (BDI) [46] and above 25 on the mini-mental state exam (MMSE) [47] or above 23 on the Montreal Cognitive Assessment (MoCA) [48]. The study was approved by the institutional Helsinki committee and all participants gave their written informed consent.

PD participants on stable anti-Parkinsonian therapy for at least one month before enrollment were subdivided into 2 groups. The first, ESPD, that formed the neuromarker identification group, consisting of patients who had been diagnosed for 2 years or less, were on Hoehn and Yahr stage 1 or 2 [49] and could be on medical treatment for PD except for Levodopa. A second group, established PD, included patients who had been diagnosed for over 2 years and could be treated also with Levodopa. Data from this group were used for further exploration of the neuromarker.

Study procedures

All study procedures were performed at the Movement Disorders Institute. For the current study, each participant visited the site on 2 or 3 occasions within a period of 2 weeks: once for clinical evaluations and twice for EEG recording (in case of 2 visits, clinical evaluation was performed on the first visit). In the current analysis, only data from the first EEG recording visit were used.

Clinical assessments

All PD patients were evaluated according to the Hoehn and Yahr staging scale and the motor examination (part III) of the Unified Parkinson’s Disease Rating Scale (UPDRS) [50] in the off-medication condition. Additional clinical evaluations included the MMSE or the MoCA, and the BDI. For cognitive screening, PD patients were administered the MMSE or the MoCA, while HC completed only the MMSE. For analysis, we used a conversion from MoCA to MMSE following Roalf and colleagues [51,52].

Stimuli and experimental paradigm

EEG recording was conducted in the off-medication condition, while participants performed the AOB and VGNG cognitive tasks (participants in the rTMS study completed 3 additional tasks not reported here). The AOB task included a series of 400 tones, 80% frequent (2000 Hz), 10% target (1000 Hz) and 10% novel complex sounds, which varied across trials (S1A Fig). Sounds were presented at a rate of 1 every 1.5 sec at 70 dB SPL. Each task contained two blocks of 200 trials separated by a 1 minute break and took approximately 15 minutes to complete. Subjects were required to respond to the target tone by pressing a key. During the VGNG task, letters were presented for 150 ms, either target or non-target (white on black background, height of 1.64 degrees, flanked by white vertical lines which remained constant on the screen throughout the block). The next letter was presented after a blank period of 850–2450 ms. There were 400 trials in total. Participants were instructed to press the button with their right index finger as quickly and accurately as possible in response to the frequent target letters (B, C, D, E, F, G) presented in the center of the screen (the Go condition; 80% of all trials) and not to respond to the rare non-target letter ‘X’ (the No-Go condition; 20% of all trials) (S1B Fig). For both tasks, there was a short practice period before EEG recording started. Only epochs with correct responses were analyzed.

EEG recordings

Participants were seated in a quiet room, at a distance of 70 cm from a 19-inch LCD monitor. They were instructed to avoid eye and body movements as much as possible and maintain their gaze at the center of the screen during task performance. Clinical evaluation and EEG recording took approximately 3 hours to complete. High-density EEG was recorded using the 64 Ag-AgCl electrodes BrainCap TMS with Multitrodes system (BrainProducts, Gilching, Germany). The reference electrode was FCz, and data were sampled at 250 Hz.

EEG processing

The EEG signals were recorded and cleaned as previously described [42]. Briefly, the EEG was referenced to average mastoid electrodes, and signals were band-pass filtered into overlapping physiological frequency bands of delta (0.5–4 Hz), theta (3–8 Hz), alpha (7–13 Hz) and beta (12–30 Hz). Then the recordings were cut into epochs based on stimulus onset and response, and averaged across trials of the same condition, producing ERPs. For the VGNG data, epoch segments (200 ms pre-stimulus to 800 ms post-stimulus) were averaged separately for the “Go” and the “No-Go” stimuli. In addition, motor-related activity was assessed by analyzing response-locked ERPs. Epoch segments (400 ms pre-response to 500 ms post-response) were averaged for the response-locked Go condition. For the AOB data, epoch segments (200 ms pre-stimulus to 1200 ms post-stimulus) were averaged separately for the “Frequent”, “Novel” and “Target” sounds. These ERPs were used as input data for the BNA analysis.

Brain Network Analytics (BNA)

The BNA methods have been previously described [42]. Briefly, BNA produces a set of event-related spatiotemporal parcels (hereafter referred to as ‘STEPs’) of the ERP activity. Each STEP contains information about the magnitude, temporal and spatial features of the ERP including amplitude, latency, spatial location (left-right and posterior-anterior, each on a scale of -1 to 1) and a topographic similarity score (S2 Fig). After STEPs are generated for each participant in each frequency band of the ERP, clustering is applied to generate group STEPs that represent spatiotemporal events common to at least 70% of the subjects in the group. A pool of healthy subjects divided into age bins (120 subjects per bin) was used to generate each STEP. For individual participants, STEPs corresponding to group STEPs in a matching process are selected and their attributes are calculated (i.e., BNA scores). The topographic similarity scores represent the similarity of the participant’s STEP to the normative group’s STEP in terms of shape and dynamics around the peak, after spatially and temporally aligning the peaks. STEP attributes served as input data for the ML algorithm as discussed below. In addition to the BNA scores, the ERP variability (ERPv) was calculated. This measure considers the variability of the ERP waveform (at a frequency band of 0.5–30 Hz) across multiple trials in a single electrode and is based on the standard error (SE). Low ERPv reflects good repeatability across multiple trials, while high ERPv reflects the variability of the neural response to the same stimulus. The ERPv is averaged across electrodes in each of nine scalp regions as detailed in S2 Table.

Machine-learning based classification

We used 199 engineered features extracted by the BNA analysis from the HC and ESPD groups as input data for the ML model training to identify a neuromarker associated with ESPD. To lessen over-fitting, we used the False Positive Rate (FPR) feature selection method to reduce the number of input features for training a logistic regression (LR) classifier. The LR classifier predicted the probability of being healthy (i.e., [1- probability for ESPD]). The classification model was evaluated with a 10-fold cross-validation stratified by group (in each iteration, data of 3 HC and 2 ESPD were left out). In each iteration, the FPR feature selection method was applied, followed by the LR classifier based on the selected features. Hence, each cross-validation iteration resulted in a model, which coefficients were used to calculate the scores for the left-out subjects (i.e., cross-validation scores). The scores from the cross-validation iterations were also used for the ROC analysis. Finally, we ran feature selection and fitted a model on the full data set of the HC and ESPD groups, to get the final model. The coefficients of the final model were used to calculate the linear term for the established PD group (see below). Fig 1A shows the weight of each feature in each iteration (i.e., fold). The weight of each feature in the final model was used to determine its importance. In the logistic function: (1) b_i is the coefficient for BNA feature i; x_i is an observation value for BNA feature i. The linear term, marked in grey in Eq (1), is a weighted sum of the input features.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. The BNA neuromarker characteristics.

(A) Weights of each feature in each iteration (fold) of the cross-validation process. The consistency of the selected features in each fold can be observed. The names of the 5 most important features of the final model are indicated at the top. (B) BNA neuromarker scores for participants in the ESPD and the HC groups. The scores were obtained from the cross-validation process and represent [1- probability for ESPD]. Whiskers extend to the highest and lowest observations, boxes represent 25–75 interquartile ranges, and horizontal bars represent medians. Mann-Whitney U test (normal approximation): U = 120, Z = -3.38, p-value = 0.0007. (C) Receiver-operating characteristic (ROC) curve for the BNA neuromarker. The cross-validation scores were used for the ROC analysis. topo, topographic similarity; PD, Parkinson’s disease; HC, healthy controls.

https://doi.org/10.1371/journal.pone.0261947.g001

As an exploratory analysis, the relevance of the BNA score to PD symptoms was assessed by correlating the outcome of the linear term of the final LR model with the motor UPDRS (mUPDRS) scores from the established PD group. We took the linear term values rather than the probability values for checking the relevancy of the model to PD symptoms, since we used a linear correlation analysis and considering the non-linear nature of the logit function. It is important to note that the data of the established PD group were not part of the training set.

Statistical analysis

Statistical analyses were carried out using JMP version 12 software (SAS Institute Inc., Cary, NC, 1989–2019). The majority of the demographic, clinical and task performance measures were not normally distributed and therefore they were compared between groups by the non-parametric Wilcoxon rank-sum test or a Chi-square test. P-values <0.05 were considered statistically significant. Performance measures were adjusted for multiple comparisons with Bonferroni correction. Response time was measured using the median response time of correct responses. Correlations were assessed using Pearson’s correlation coefficient.

Cognitive task performance

Behavioral results of the VGNG and AOB tasks in the ESPD and HC groups, are shown in Table 2. One ESPD patient performed the AOB below chance level and was therefore not included in the AOB performance analysis. There were no significant differences in the AOB task performance between groups. In the VGNG task, accuracy was significantly lower in the ESPD group than in the HC group.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Results of the cognitive tasks performance.

https://doi.org/10.1371/journal.pone.0261947.t002

Discriminant ability of the BNA neuromarker

The output of the BNA system was analyzed by the ML algorithm (see Methods and Fig 1A), which identified a BNA neuromarker composed of 15 features. Box plots with the neuromarker scores from the cross-validation process for the ESPD and HC groups are shown in Fig 1B. The cross-validation scores are more representative in terms of model validity, in the case a validation set is lacking.

The discriminant ability of the BNA neuromarker was assessed using a receiver-operating characteristic (ROC) analysis (Fig 1C). Again, the cross-validation scores were used for the ROC analysis. The area under the ROC curve was 0.79 (95% CI: 0.65, 0.93). At the cutoff point that maximizes sensitivity while specificity is still above chance level, sensitivity was 74% and specificity was 73%. The ROC curve data in S3 Table shows various sensitivity-specificity values for different cutoffs.

While the control’s gender was balanced, the ESPD group was predominantly male (84%). This could potentially lead to a false neuromarker identification due to gender differences in the groups rather than disease-related difference. To check this possibility, we first examined the difference between the neuromarker scores of males and females in the HC group, which is gender balanced. A Wilcoxon rank-sum test showed no statistically significant effect (Z = -0.83, p-value = 0.41). Second, we ran 5 iterations of the ROC analysis with balanced gender across groups, i.e., taking each time data of randomly chosen three females (with no repeat) from the HC group and all the 15 HC males. The average area under the curve was 0.77 (range 0.76–0.78). At the cutoff point that maximizes sensitivity, while specificity is maximal above chance level for this sensitivity value, sensitivity was 0.74 (in all iterations) and average specificity was 0.69 (range 0.67–0.72). These results are very similar to those obtained with all HC females and thus suggest that gender differences are not the basis of the neuromarker discrimination between patients and HC.

Components of the BNA neuromarker

The weights of the 15 features of the BNA neuromarker are shown in Fig 2A. The most important features, as determined by their weights, included the P50 amplitude from the No-Go condition and the N100 latency from the Novel condition, both of which are associated with early sensory processing; the P200 amplitude from the Frequent condition and the P200 topographic similarity from the Novel condition, associated with filtering of information; and P-200 topographic similarity from the response-locked Go condition, activity which precedes the motor response and is suggestive of reflecting motor-related processes. These features were ERPs in the alpha frequency band. Other features of the BNA neuromarker were ERPv scores of frontal and central locations from the response-locked and stimulus-locked VGNG ERPs. Two other components with relatively modest contributions to the BNA neuromarker were related to cognitive control and novelty processing (Novel N2 P-A in alpha band) [53] and Go processing including response preparation and execution [54] (Go P300 P-A in delta band; a “target” in the sense that it requires a response) [55]. See S4 Table for a summary of physiological explanations of the features.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Features of the BNA neuromarker.

(A) Weights of the BNA neuromarker features. The box indicates the five most important features. Weights are absolute value model weights normalized by their sum. (B) Box plots of the BNA scores of the five most important neuromarker features in participants with ESPD and HC. See Fig 1 for boxplot description. Points outside the whiskers are outliers, exceeding 1.5 times the interquartile range below the 1st quartile or above the 3rd quartile. Values are provided in S1 Table. AOB, auditory Oddball; VGNG, visual Go/No-Go; Topo, topographic similarity; ERPv, event related potential variability; μV, micro volts; ms, milliseconds.

https://doi.org/10.1371/journal.pone.0261947.g002

The differences between the ESPD and the HC groups, in the neuromarker features with the highest importance, are shown in Fig 2B. Compared with HC, ESPD patients had higher No-Go P50 amplitude, slower Novel N100 latency, lower Frequent P200 amplitude, and lower Novel P200 and response-locked Go P-200 topographic similarity.

Further exploration of the neuromarker

Twenty established PD patients went through the same procedures as the ESPD patients (see S5 Table for demographic and relevant clinical characteristics). The established PD patients’ scores from the linear term of the LR model (see Methods) correlated significantly with their mUPDRS scores (r = 0.70, p = 0.0008, 95% CI: 0.36, 0.88; Fig 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Correlation between the BNA neuromarker linear scores and the mUPDRS scores.

Correlation between the BNA neuromarker linear scores and the mUPDRS scores in an independent established PD group (N = 19; one subject was assessed only in the on-state and was not included). mUPDRS, the motor examination (part III) of the Unified Parkinson’s Disease Rating Scale.

https://doi.org/10.1371/journal.pone.0261947.g003

Lastly, we calculated for each early-stage patient the levodopa equivalent daily dose (LEDD), a measure of the total daily medication a patient receives, which enables comparison across patients with different drug regimens [56]. LEDD was then correlated with the neuromarker score. As can be seen in Fig 4, there was no correlation between the neuromarker and the LEDD (r = 0.075, p = 0.760, 95% CI: -0.39, 0.51).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. BNA Neuromarker scores as a function of LEDD.

LEDD, levodopa equivalent daily dose.

https://doi.org/10.1371/journal.pone.0261947.g004