The battery of neuropsychological tests administered to subjects in V0 and V21, included the MMSE [29] that is a measure of global cognition, the FCSRT [30], a memory test that controls attention and cognitive processing to better investigate episodic memory performance, as well as other tests focused on attention and executive functions [31], such as the Trail-Making Test A and B (TMT A and B; [32]), the Digit span forward and backward, and the Frontal Assessment Battery (FAB; [33]). Additional instruments included the 15-item version of the McNair Frequency of Forgetting Questionnaire [34], the Geriatric Depression Scale (GDS; [35]), and the State-Trait Anxiety Inventory (STAI A and B; [36]) for assessing, respectively, memory complaints, mood, and anxiety disturbances.

At V0, and upon completion of the neuropsychological tests, EEG data resting state was acquired for 150 seconds: 75 seconds with eyes closed, and 75 seconds with eyes open looking at a fixation cross on the screen. A similar EEG acquisition was performed at V21, the post-training visit, for the subjects who completed the full training protocol. For both sessions (V0 and V21), the same EEG amplifier with the same recording parameters were used: a 20-electrodes cap with a traditional 10–20 electrodes layout and a NIC® EEG amplifier system (Neuroelectrics, Barcelona, Spain). The left earlobe electrode was used as the EEG reference [37]. The impedance of all electrodes was kept below 10 k$\mathrm{\Omega }$. The sampling rate of EEG signals with the NIC system was 500 Hz, and no filter was applied during the acquisition phase. For EEG data analysis at pre- and post-training visits, we considered only the five electrodes on fronto-central and parieto-central sites (Cz, Pz, Fz, C3, and C4), in accordance with the low-density EEG system used for neuromodulation. More details are given in the next section.

The training started up to seven days after the pre-training visit (V0). It was performed in three blocks of 10 minutes per biomarker (30 minutes per session) during 20 visits over three months.

In practical terms, when the subject arrived, he/she was invited to sit in front of a screen. For the first session, the investigator explained to the subject how the EEG is performed and how the feedback will appear on the screen. Then, the technician places the EEG system (Neuroelectrics, Barcelona, Spain) on the scalp and two electromyography (EMG) sensors on the trapezius muscles of each subject (from either the real or the sham feedback group). The feedback information was provided visually to each subject as a picture of a tree on her/his computer screen (see Fig. 2 ). There were three different tree models: for each session, a different tree was assigned randomly to each biomarker. Subjects were only instructed to concentrate in order to make the tree grow. All subjects were asked to make the tree grow the same way; however, the feedback presented to subjects of Group A was based on their brain dynamics, while the feedback presented to subjects of Group B was based on their muscular activity. All precautions were taken to prevent the experimenters and subjects from identifying the task type (neurofeedback or shame/control). The task design was performed with the MATLAB Psychtoolbox (version 3, https://www.psychtoolbox.net/) [38], and the EEG synchronized and the stimuli presentation on the screen was performed with the System Level Simulations (SLS) library (https://fr.mathworks.com/matlabcentral/fileexchange/134292-wireless-system-level-simulations-sls-event-logger-) [39].

Fig. 2.

Brain-Computer Interface training. Left: Picture of two subjects during a session of the training phase. Right: The image of the feedback task is displayed on the screen.

At the beginning of each visit, one minute of the resting state EEG was recorded. These data were used to compute the threshold values of the subject’s biomarkers. During training, each biomarker was computed, in real-time, from the recorded signals. If the actual computed value of the biomarker was higher than the threshold, the tree grew up; otherwise, the tree size remained constant. The same protocol was followed for subjects submitted to sham feedback.

For neuromodulation, we built an EEG-neurofeedback BCI based on a small set of EEG electrodes. Indeed, as mentioned earlier in the introduction, our aim was to study the feasibility of EEG-neurofeedback in a constrained context with a limited number of electrodes for a low-cost and portable BCI. We considered five fronto-central and parieto-central electrodes: Fz, Cz, Pz, C3, C4. These electrodes were chosen based on experiments that we had carried out on an EEG database of 22 subjects with Subjective Cognitive Impairment (SCI) and 28 mild AD patients. EEG data were recorded during resting state eyes-closed condition, using a Deltamed digital EEG acquisition system (http://medicalsystems.ma/) with 30 scalp electrodes positioned over the whole head according to the 10–20 international system [19]. To discriminate between SCI and AD groups, we considered different combinations of a small number of electrodes among the 30 available, on which we computed an average Phase Synchrony measure [40] at different frequency bands. Then, such average synchrony value was used as input of a Linear Discriminant Analysis classifier to evaluate the discrimination performance in terms of specificity (SCI correctly classified) and sensitivity (AD correctly classified). The best performances (specificity = 81.8% and sensitivity = 71.4%) were obtained using EEG electrodes of fronto-central and parieto-central sites, namely Cz, Pz, Fz, C3, and C4.

Thirty-seven subjects were initially included and randomly assigned into two groups, called Group A and Group B for preserving the double-blind approach. Six subjects were desisted for personal and scheduling reasons at the beginning of the protocol, thus explaining the unbalanced number of subjects in the groups: 17 subjects in Group A and 14 subjects in Group B.

Neuropsychological and EEG data analyses were performed on the data collected in pre- and post-training visits from the 31 subjects who completed the study. As summarized in Table 1, a significant difference in terms of gender ratio was observed between the two groups, with a more significant proportion of women in Group A (p = 0.05, Fisher’s test), while no significant differences were found between the groups for age and education level. Note that categorizations of the education level into “high” and “low” or “intermediate” are as indicated in [31]. In short, more than 80% of our participants had an education level score $\ge$7, which corresponds to the label “high education” (i.e., “higher than up-secondary”, according to the International Standard Classification of Education, 2011).

Due to the high number of artifacts on eyes-open EEG recordings, the EEG analysis was performed only on the data collected from the eyes-closed condition [41], and considering only Fz, Cz, Pz, C3 and C4 electrodes, as explained in Section 2.2.3.

The pre-processing of the resting-state EEG signals was performed on MATLAB version 2016b (The MathWorks, Natick, MA, USA), using the Brainstorm software [42, 43] in combination with custom scripts developed internally within the MATLAB environment (MathWorks). The sampling rate of the data during the acquisition was 500 Hz, and no downsampling was applied in the pre-processing steps.

Data cleaning includes space separation to remove eye blink and saccade artifacts. Then, EEG signals were notch-filtered at 50 Hz to remove the power line noise and band-pass-filtered in the four frequency bands corresponding to theta [4–8 Hz], alpha [8–13 Hz], beta [13–25 Hz], and gamma [35–45 Hz] bands. The filtered data, which had a total duration of 75 seconds per acquisition and per subject, were segmented into 2-second epochs. All the epochs with a kurtosis below one or above five, or a peak-to-peak amplitude exceeding 100 µV, were considered rejected from the analysis. The remaining epochs were visually checked, and the bad segments were identified and rejected. The final number of epochs without artifacts per acquisition ranged from 32 to 37 for each subject. For the analysis, we kept the same amount of data from each subject, thus retaining 32 epochs of 2 seconds per subject.

Finally, values of the three biomarkers of interest were extracted from the pre-processed epochs at pre- and post-training visits, and considering only the five electrodes (Fz, Cz, Pz, C3, C4): (i) the PAF that measures the frequency with the highest magnitude within the alpha rhythm sub-band spectrum after estimating Welch’s Power Spectral Density (PSD), on 1 s-window with 25% overlapping; (ii) Gamma-band synchronization that relies on the computation of Phase Synchrony [40] between all pairs of EEG signals in gamma band; (iii) the TBR calculated by dividing average theta power by average beta power after a Welch’s PSD estimation on 1 s-window with 25% overlapping. More precisely, the PAF and TBR biomarkers were computed on each epoch from the signal captured by a single electrode; hence, five features were extracted per epoch from the PAF and five from the TBR. By contrast, Gamma-band synchronization requires two electrodes; 10 pairwise combinations of five electrodes exist, and ten features were computed from that biomarker. Therefore, for each subject, 20 features were extracted per epoch.

To assess the effectiveness of our EEG-neurofeedback framework, we adopted two statistical approaches for data analysis based on classical statistical tests and machine learning methods [44], as explained in the following sections.

Since most neuropsychological and EEG data were not found to be normally distributed, non-parametric tests were performed in the analysis (for further details, see section 3). All statistical analyses were conducted using R version 3.6.1 (R Development Core Team, 2019, https://www.r-project.org/), and plots were generated with the ggplot2 package (https://cran.r-project.org/web/packages/ggplot2/index.html) [45]. The level of statistical significance was set at p 0.05 for all tests. To avoid type-I error, correction for multiple comparisons was completed by using the Bonferroni method [46].

To evaluate the effectiveness of EEG-neurofeedback training, we conducted both inter-group and intra-group data analyses. Therefore, four classifiers were designed for four different purposes: (i) discriminating subjects of Group A from subjects of Group B, using the EEG epochs recorded during V0 (inter-group analysis at V0); (ii) discriminating subjects of Group A from subjects of Group B, using the EEG epochs recorded during V21 (inter-group analysis at V21); (iii) discriminating epochs recorded from subjects of Group A during V0 from epochs recorded from the same subjects during V21 (intra-group analysis, V0 versus V21 within Group A); (iv) discriminating epochs recorded from subjects of Group B during V0 from epochs recorded from the same subjects during V21 (intra-group analysis, V0 versus V21 within Group B).

As above-mentioned, each subject (among the 31) has 32 epochs at session V0 and 32 epochs at session V21, all contributing to the data set available for classifier design. On each epoch, 20 features were extracted (five PAF values, five TBR values, and 10 Gamma-band synchrony values). To compute such features and train the classifiers, we used the SIGMABox (https://github.com/tmedani/SIGMAbox) [47], a homemade MATLAB toolbox for EEG signal processing and classification.

As all features may not be equally relevant for the classification, a subsequent feature selection procedure was applied to discard irrelevant features. The procedure has two steps: Orthogonal Forward Regression (OFR) for ranking the features in order of decreasing relevance, followed by the random probe method to find the threshold rank for rejecting irrelevant features [48]. The random probe method consists of generating random variables (called probes) by shuffling the components of the vectors of the candidate variables (corresponding to the real features), ranking the candidate variables and the probes in a single ranked list, and estimating the acceptable risk value that a random variable can explain the output process more reliably than one of the selected real features.

For inter-group analyses (classifiers (i) and (ii)), the number of epochs available for the design of the classifier was 32 epochs per subject, hence 992 examples were available for classifier design. For intra-group analyses, the number of epochs available for the design was 64 epochs per subject, as the epochs to be classified were recorded during V0 and V21, hence 1088 and 896 examples were available for classifiers (iii) and (iv), respectively. Table 2 summarizes the number of subjects and epochs.

The subjects of Group A and Group B were randomly assigned to the training set (70%) and the test set (30%). The four classifiers were designed to perform epoch-wise classification from the features; the classification of subjects was performed according to the following rule: a subject was assigned to the class to which the majority of his/her epochs (at least 17) had been assigned by the epoch-wise classifier.

We used Support Vector Machines (SVM) with linear soft-margin classifiers [44]. The estimation of the regularization constant of the SVM was performed by leave-one-subject-out cross-validation. Due to the small size of the test set, the estimation of the performance of the classifier may vary with the partition of the examples into training set and test set. In order to decrease the variance of that estimation, 100 different partitions of the complete data set were performed, ascertaining that the number of inclusions Ni of each subject i in the test set was such that 5 $\le$ Ni $\le$ 8. The whole design procedure was performed for each partition, thereby designing 100 classifiers. After completion of the iterations, the number of correct classifications for each subject when she/he was present in the test set, was computed according to the following rule: if the number of correct classifications of subject i was larger than the integer part of Ni/2, the subject was considered correctly classified. For instance, a subject with Ni = 5 was considered correctly classified if she/he was correctly classified at least 3 times; a subject with Ni = 8 was considered correctly classified if she/he was correctly classified at least 5 times.

To investigate whether the two groups (A and B) were different regarding their EEG activity at V0 (before the 20 training sessions), Mann-Whitney non-parametric tests were applied on the 20 features extracted from the three biomarkers of interest. For each subject and each of the 20 features, a value was computed for each of the 32 epochs; then the mean of these 32 values was used as the biomarker associated to each subject. As shown in Supplementary Table 4, irrespective of the biomarker or the electrode, there were no significant differences between groups at V0. This confirms that Group A and Group B were not significantly different regarding the three biomarkers of interest at the beginning of the study.

Inter-group analysis was also performed at V21 to determine whether the 20 training sessions had a significantly different effect on the EEG dynamics of the two groups (A vs B). Mann-Whitney non-parametric tests were applied on features (same procedure as presented in the previous section). As shown in Supplementary Table 5, we observed a significant difference between groups regarding the Gamma-band synchronization for three pairs of electrodes: Pz/C4, Pz/Fz, and C4/Fz, with greater values of gamma-synchrony for Group A. A significant difference between groups was also found for this biomarker when considering the mean over all electrodes (p 0.001). At V21, the groups were therefore significantly different with regard to the Gamma-band synchronization, suggesting that subjects of Group A and Group B did not modulate this neural biomarker similarly post-training.

To investigate changes in brain dynamics between pre- and post-training sessions (V0 vs V21), we compared, for each group, the feature values computed at both V0 and V21. Thus, Mann-Whitney non-parametric tests were applied separately for each group (A and B), and each biomarker. As shown in Supplementary Table 6, there were significant changes regarding the features computed for the Gamma-band synchronization for both Group A and Group B. However, for Group B (sham feedback), the difference was only observed in one pair of electrodes (Pz/Fz, p 0.001), whereas, for Group A (real neurofeedback), a more global effect was corroborated by a significant increase of Gamma-band synchronization for the mean of all channels (p = 0.006). These results are consistent with the observation reported in the previous section, thus confirming that subjects of Group A and Group B did not modulate the Gamma-band synchronization similarly post-training.

We trained a SVM classifier to discriminate subjects of Group A from those of Group B based on their EEG data at V0. Using the three biomarkers, 20 features were extracted and the most relevant ones were selected with the OFR algorithm and the random probe, as explained in Section 2.3.4. A 20% risk of selecting a candidate variable although it is irrelevant turned out to be a reasonable choice, resulting in selecting seven features (among the 20). Table 3 shows that Gamma-band synchronization biomarker is predominant in the selected variables.

These selected seven features were then used by the classifier for inter-group classification, to determine which epochs belonged to Group A or Group B. Table 4 reports the confusion matrix resulting from the cross-test procedure described in Section 2.3.4.

The accuracy of the classifier estimated by cross-test was 54.8%. It is good practice to compare the results of a classifier to the results that would have been obtained by a baseline classifier (called “zero-classifier”), which assigns all examples to the most populated class (Group A in our case), irrespective of the features. If a classifier does not significantly outperform this baseline, it indicates a limitation in classification ability. In this instance, the accuracy of our classifier matches that of the “zero-classifier”, implying an inability to distinguish between the two groups based on the three selected biomarkers. This strongly suggests that the EEG biomarkers computed from both groups (A and B) were similar at V0 before training.

Here, the purpose of the classifier was to discriminate subjects of Group A from those of Group B based on their EEG data at V21. Interestingly, at V21, among the 20 features extracted from the three biomarkers, seven relevant features were also selected with the OFR algorithm and the random probe method, all belonged to the Gamma-band synchronization. Such selected features were then used as inputs to the SVM classifier.

As presented in the confusion matrix resulting from cross-test (Table 5), the SVM classifier was able to discriminate subjects from Group A and Group B with a global accuracy of 71.0%. As the obtained accuracy is greater than 54.8% (the “zero-classifier”, see section 3.3.1), this result strongly suggests that the selected Gamma-band synchronization features were able to distinguish, post-training, subjects of Group A (with real neurofeedback) from those of Group B (sham feedback).

To validate the feature selection procedure performed by the random probe method, three classifiers (i.e., one classifier per biomarker) were designed to perform the same inter-group analysis at V21, without the application of feature selection.

As presented in Table 6, the classifier whose feature inputs were extracted from the Gamma-band synchronization had an accuracy of 71%, while the accuracy of the classifiers trained with features extracted from either the TBR or the PAF was 61.3% (estimated by cross-test).

These results confirm that Gamma-band synchronization biomarker is more pertinent for inter-group classification at V21 than the other two biomarkers. Interestingly, the accuracy of the classifier whose features were computed from the Gamma-band synchronization alone (see Table 6) considering all the 10 features, was equal to the accuracy of the classifier that considered only the seven selected features among the 10 Gamma-band synchronization features. Again, this finding confirms the result of the feature selection by the random probe method, namely that the selected Gamma-band synchronization features was particularly relevant to discriminate subjects of the two groups based on their respective epochs.

Finally, we aimed to determine whether there was a significant difference between V0 and V21, within each group separately. In other words, we wondered whether subjects of Group A and Group B, separately, were able to self-modulate their brain dynamics on the basis of the three biomarkers of interest.

Table 7,8 show that the classifiers that used features extracted from the Gamma-band synchronization and the PAF were able to discriminate between epochs recorded at V21 and those recorded at V0, only for subjects of Group A (accuracies of 67.6% and 61.8%, for the Gamma-band synchronization and the PAF, as shown in Table 7, versus accuracies of 35.7% and 32.1% regarding the same biomarkers in Group B, as shown in Table 8). By contrast, the classifiers using features extracted from the TBR were unable to discriminate V0 from V21 in Group A (Table 7) and Group B (Table 8).

It is worth noting that these results were obtained in a double-blind setting: the information that subjects who received real neural feedback belonged to Group A, and subjects who received sham feedback belonged to Group B, was disclosed after completion of data collection and analysis.