1 Introduction

Central Neuropathic Pain (CNP) is a chronic pain arising from any kind of injury to Central Nervous System (CNS) [1]. However, the most common reasons behind CNP are Spinal Cord Injury (SCI) [2], Stroke [3] and Parkinson disease [4]. Many studies have shown that presence of CNP produces changes in the resting state brain activity [5, 6] , such as increased resting state electroencephalography ( EEG) power and slower down of dominant frequency of α brainwaves [6,7,8,9]. Although, there are several pharmacological drugs, with different side effects, used for treatment of CNP [10], medication of CNP still a great struggle. There is no broadly consensus on most appropriate drug for CNP, and the treatment is often a matter of trials [11]. Additionally, pharmacological drugs usually aims at tackling with symptoms not with cause of CNP [12], and symptoms respond inadequately to drugs [13]. The side effects and limited efficacy of existing pharmacotherapy of CNP have increased the need to non-pharmacological treatments based on neuromodulation and neurostimulation such as meditation, transcranial Direct Current Stimulation (tDCS) [14] and neurofeedback [15].

Neurofeedback Training can be alternative to pharmacological drugs for treatment of CNP [7, 16]. By definition, neurofeedback is an approach to training a brain’s participant to reinforce a desired brainwave or to inhibit a brainwave thought to be associated with the pain occurrence [17]. Neurofeedback usually carried out employing EEG electrodes attached to the scalp to provide the participant with almost instantaneous assessment of the brain functioning in a visual or an auditory form [18]. The ability to self-regulate brain activity [7] is the root cause made it is possible to apply neurofeedback for treatment various types of chronic pain, including CNP [8, 16]. Neurofeedback training applied for treatment of CNP aims at rewarding alpha (α) and inhibiting theta (Ɵ) and beta (β) brainwaves [16]. The efficacy of neurofeedback training, as a pain management therapy, is basically evaluated by change in pain intensity, scores of behavioral questions, and changes in EEG waves before and after treatment [19, 20]. Machine learning (ML) based approach may be applied to develop automatic method in order to evaluate the efficacy of neurofeedback training applied for treatment of CNP, as the presence of CNP has signatures of EEG [21, 22]. The most challenging task in this regard is to extract discriminative features characterizing CNP in raw EEG[23]. However, performing motor imagery (MI) tasks during EEG recording may result in obtaining more discriminative features distinguishing CNP in EEG, as the amplitude and energy of EEG of patients with CNP significantly change during MI of a limb[21]. The more discriminative features associated with CNP in training data result in more precise and accurate brain computer interface(BCI) for detecting CNP [23].

At present, research studies designed for EEG investigation based on ML are growing actively in research community. Most active application of ML in EEG based studies is the classification of EEG patterns associated with different disorders such as epilepsy[24], and mental activities such as MI [25]. Acharya has investigated the possibility of applying ML for diagnosis of epileptic EEG and showed that classification accuracies of applied classifiers were ranged from 88.1 to 98.1% [26]. Maksimenko applied ML to classify EEG associated with MI, and the classification accuracy reached up to 95% [25]. Moreover, ML has been used to analyze post-treatment EEG based on pre-treatment EEG in order to assess the efficiency of clozapine therapy in adult patients suffering from chronic schizophrenia [27]. Accordingly, we carried out this study to investigate the possibility of applying ML algorithms for automating the process of analyzing post-neurofeedback EEG. The principal goal in this study was to develop effective BCI-classifier in order to automatically evaluate the efficacy of neurofeedback training applied for treatment of CNP. Therefore, a variety of ML algorithms were fed with training data extracted from pre-neurofeedback EEG in order to find out the optimal BCI-classifier. With the aim of obtaining more discriminative features associated with CNP, participants in this study were instructed to perform MI of limbs during EEG recording. For further optimization of BCI-classifier, classifiers training carried out with various EEG channels combinations.

This study does not assume that the CNP will disappear after neurofeedback, but there would be a significant reduction in the pain to the extent that EEG signatures associated with CNP would have no significant influence in post-neurofeedback EEG [28]. Based on this, our hypothesis is that EEG data recorded after neurofeedback should be classified in no-pain or able-bodied group when provided as an input to the developed BCI-classifier.

2 Methods

2.1 Participants

The study is divided in two phases. In the first phase which called as Before Neurofeedback(BNF), experiments were performed in three groups of participants named as able-bodied (AB), paraplegic patients with (PWP) and with no (PNP) neuropathic pain. In all, thirty subjects participated in the BNF phase of this study (10 in each group). The age of participants was between 18 and 55 years (AB group: 3F, 7M, age 39.6 ± 10.2; PNP group: 2F, 8M, age 43.8 ± 9.1; PWP group: 3F, 7M, age 46.2 ± 9.4). In the second phase which called as After Neurofeedback(ANF), experiments were performed on five PWP participants (recruited from the first phase) who have completed neurofeedback training provided for the management of pain. Informed consent was obtained from all participants, and the ethical approval to the trials was granted from the University of Strathclyde Ethical Committee and from National Health Service Ethical Committee for Greater Glasgow and Clyde.

2.2 BNF data acquisition

Neuroscan EEG machine with 61 electrodes (according to standard 10–10 locations) was used to record task-related EEG at sampling rate of 1000 Hz. Participants were trained to perform three types of motor imaginary tasks which are imagination of the right hand waving (RH), left hand waving (LH), and tapping with both feet (F) for 3 sec excluding 1 sec resting and 1 sec preparation periods. Therefore, the length of each epoch is 5 sec long. In total, EEG data was recorded with 60 trials of each motor imaginary task (180 trials in total) carried out by each subject.

2.3 Neurofeedback training

Neurofeedback training was provided to five PWP participants with the aim of reinforcing α band power, and inhibiting Ɵ and β. The alpha band power was reinforced because it was reported in some previous neurofeedback studies that increasing α band has promising results on chronic pain including CNP [29,30,31]. Additionally, the group of PWP participants in our study had the dominant α frequency on average 1 Hz lower than that of PNP participants. Therefore PWP participants were trained to increase the energy of higher α band (9–12 Hz) which does not include lowest α frequencies at 8Hz. The β and Ɵ power was inhibited because these frequencies bands are believed to be positively correlated with pain [21, 30].

Neurofeedback training was conducted firstly with an audio neurofeedback during which patients sat in front of a computer screen listening to music. The level of music was louder when α power was under the threshold, set at 110% of the baseline value, and quieter when α power was above the threshold. Following the audio neurofeedback, visual neurofeedback was provided to the patients in order to decrease β and Ɵ power to a value under the threshold set at 90% of the baseline value. In order to visualize α, β and Ɵ power changes during neurofeedback, three bars were presented on the computer screen. The changes of bars size represent a change of relative EEG power in three selected frequency bands as shown in Fig. 1. The bars color turned green when the EEG power of the chosen frequency band was in the desired range and turned red otherwise. The bar in the middle represented α band and turned green when the power was above the threshold set at 110% of the baseline value. Two sidebars represented β and Ɵ frequency bands and turned green when the power was below the threshold set at 90% of the baseline value.

Fig. 1
figure 1

The graphical user interface used for online evaluation of neurofeedback training. Horizontal black dotted lines show threshold levels. Middle bar shows power of dominant frequency band (alpha), which turned green when the power was above the threshold (reinforced). Side bars represent theta, and beta frequency bands that were suppressed

As the five PWP participants, who completed 20 or more neurofeedback sessions, reported to have a sensation of pleasant warmth, the ANF EEG power in three targeted frequency bands (α, β and Ɵ) was compared to BNF EEG power. Throughout the comparison, it was noticed that there was a significant change in ANF EEG power as the β and Ɵ were suppressed and α power increased after neurofeedback. The change of EEG power after neurofeedback in two patients, (PWP1, PWP2), as represented in Fig. 2 was more noticeable.

Fig. 2
figure 2

The EEG power after and before neurofeedback training in patient PWP1 and PWP2. Solid black line refers to BNF (before neurofeedback training) EEG power, and dotted black line refers to EEG power after neurofeedback (ANF) in two patients( PWP1, PWP2)

Following neurofeedback training of five PWP participants, EEG data was obtained from them in order to test the accuracy of BCI-classifier developed, in first phase of this study, based on the pre-neurofeedback EEG. ANF data was obtained from the five participants in the same way and with the same EEG machine used in the first phase. The participants were instructed to perform MI of limbs during EEG recording in this phase as well.

2.4 Data preprocessing

EEG data was down-sampled to 250 Hz. An infinite impulse response high pass filter with 12db cutoff frequency was set to 1 Hz and a notch filter was applied between 48 and 52 Hz to eliminate line noise at 50 Hz. Before performing the analysis, EEG data were re-referenced to the average reference because the evoked activity, or signal of interest, has a more or less fixed time-delay to the stimulus, while the ongoing EEG activity behaves as additive noise. The averaging procedure will enhance the signal-to-noise ratio.

2.5 Features extraction

In this study time–frequency decomposition was performed in a frequency range 3–55 Hz using a sinusoidal wavelet with minimum 3 wavelet cycles per data window. Following this, ERD/ ERS was averaged for the theta (4–8 Hz), alpha (8–12 Hz), and beta (16–24 Hz) frequency bands in four time windows (0.4–0.8s, 0.8–1.2s, 1.2–1.6s, 1.6–2.4s). These frequency bands and time windows were selected on the basis of results reported in our previous study [32].

Data normalization considered to be one of the main processing steps for learning tasks, and it can be defined as scaling technique that helps finding new range from an existing one range. In this study, Min-Max normalization technique was applied to rescale data to [0, 1] range.

2.6 Classification

The classification was performed on EEG data recorded in the first phase of this study (BNF). Decision tree(DT), Naïve Bayes(NB), Support vector machine(SV), K-nearest neighbor(KNN), and Artificial Neural Network(ANN) algorithms are employed in this study as classifiers in order to find the most appropriate algorithm to build the predictive model for studying the efficacy of NF training. The classification process was carried out with four different combinations of electrodes which were called as C1, C2 and C3 that included 61, 27 and 9 electrodes respectively as shown in Fig. 3. The combination C1 was chosen because studies reported global effect of pain while C2 was chosen on the basis of definition of NP reported lesion to the somatosensory cortex. The combination C3 were chosen on the basis of somatotopoical representation of upper and lower limbs and cortical shift of paralyzed limbs reported in the literature. In this study, the efficiency of each classifier was evaluated by taking 70% and 30% data for training and testing respectively.

Fig. 3
figure 3

The locations of selected combinations of EEG electrodes on the scalp

2.6.1 DT classifier

DT algorithm builds a classification tree with nodes, branches and leaves (leaf defined as a single class which doesn’t split any more) based on selected predictors set. The DT classifier in this study was grown using fitctree function that splits categorical predictors applying the exact search algorithm. The number of discriminative predictors and their levels decide the number of the nodes that determines the size of tree. The numbers of nodes were 21, 101, and 43 in case of classification of MI data of RH for three group combinations AB vs PWP, AB vs PNP, and PNP vs PWP respectively. However, in case of data of MI of F, the numbers of nodes were 77, 69, 67 for all three groups’ combinations AB vs PWP, AB vs PNP, and PNP vs PWP respectively. Similarly, the numbers of nodes were 33, 25 and 7 in case of classification of MI data of LH. The reason for choosing DT classifier for this study is its effectiveness in implementation and it has been used for analyzing EEG data in many other studies as well [33].

2.6.2 NB classifier

The Normal distribution was employed to estimate the distribution of the data since it is the easiest to work with because it is only required to estimate the mean and the standard deviation from training data. Based on this, NB classifier was selected for this study as it is powerful and fast learning algorithm. Moreover, NB is reliable approach to analyzing EEG [34].

2.6.3 K-NN classifier

K-NN classifier does not use any model to fit but is only depending on memory, and a new instance is classified based on the closest training samples present in the feature space. When a test data is entered, it is assigned to the class that is most common amongst its k nearest neighbors [26].The optimal value of k (number of nearest neighbors) in this study was 3; and this value found through experiments. KNN classifier was chosen for this study as it was reported that it gained high accuracy in EEG classification based studies [35].

2.6.4 SVM classifier

SVM classifier in this study was used for binary classification with linear Kernal Function and supports sequential minimal optimization. SVM classifier is used on a large scale in EEG based studies for being powerful classifier based on pattern recognition [36, 37]. Additionally, SVM for binary classification is considered to be more effective for EEG classification based on pattern recognition approach [38].

2.6.5 ANN classifier

The ANN classifier applied in this study as a classifier was selected with the following parameters. The number of neurons in the hidden layer was chosen according to experiments that revealed that the optimal number of neurons in our case of EEG classification is 20 neurons. The number of neurons in the output layer was one, because the ANN classifier was selected to perform binary classification. The network training function was trainscg that updates weight and bias values according to the scaled conjugate gradient method. The ANN trained on back propagation algorithm. ANN classifier was selected for this study based on literature review showed that ANN classifier gained high accuracy in classification of EEG patterns [37]. Moreover, ANN was applied successfully for classification of EEG patterns related to MI [25].

3 Results

Grayscale images in Fig. 4 show classification accuracy of between each two groups calculated for all three types of MI tasks. Each image shows how the classification accuracy of MI patterns associated with CNP changes with the type of participants group, classifier and electrodes combination. Image (A) represents the classification accuracy of RH-MI associated with CNP, and it shows that the highest classification accuracy of RH-MI was gained with AB vs PWP groups, with electrodes combination C1 and with all classifiers especially with SVM which obtained accuracy of 99±0.49%, and the lowest accuracy was obtained with PWP vs PNP group and C3. All classifiers obtained relatively low accuracy with PWP vs PNP except SVM that obtained accuracy of 77 ± 1%. Image (B) illustrates the classification accuracy of F-MI associated with CNP, showing that highest accuracies in classification of F-MI were obtained with AB vs PWP and AB vs PNP groups, with C1 and with SVM classifier. SVM gained accuracy of 95 ± 1.1% and 97 ± 0.6% with two groups respectively; and the lowest accuracy was obtained with C3 and PWP vs PNP with all classifiers; and SVM obtained accuracy of 88 ± 7%.

Fig. 4
figure 4

The classification accuracy of five classifiers per three combinations of electrode (C1, C2, C3) calculated between each two group (AB vs PWP, AB vs PNP and PWP vs PNP) during MI of RH (subfigure (A)), F (subfigure (B)), LH (subfigure (C)). The classifiers’ names are Decision Tree (DT), Naïve Bayes (NB), Support Vector Machine (SVM) and Artificial Neural Network (ANN). Black colour in all three subfigures showed high accuracy while light grey color represents low value of accuracy

Image (C) demonstrates the classification accuracy of LH-MI associated with CNP, showing that AB vs PWP, C1 and with SVM gained the highest accuracy of 98±0.36%; and the lowest accuracy of 91 ± 1% was gained with PWP vs PNP, C3 and with SVM.

Figure 5 showed the sensitivity and specificity of SVM based BCI-classifier. Sensitivity refers to percentage of positive cases (AB or PNP) predicted correctly as specificity refers to percentage of negative cases (PWP) predicted correctly with combination 1 (C1). It is obviously noticed that specificity is always higher than sensitivity and that implies that the ability of SVM based classifier to detect PWP cases is higher.

Fig. 5
figure 5

Average of specificity and sensitivity for data classification between each two group (AB vs PNP; subfigure (a), AB vs PWP; subfigure (b), PNP vs PWP; subfigure (c),) calculated with SVM classifier with electrodes combination (C1). Black bars represent sensitivity and white bars represent specificity

Figure 6shows the prediction accuracy of SVM-classifier developed in the first phase of this study. SVM-classifier in this case was fed with ANF-EEG data obtained from five PWP who received neurofeedback training applied for treatment of CNP. With the analysis of ANF-EEG power, it was clear that there was change in EEG power especially in α, Ɵ and β frequency bands as shown in Fig. 2. Therefore, we assumed that the success of developed SVM-classifier lies in its ability to detect these changes in EEG power and to classify ANF-EEG in no-pain, or able-bodied group. It can be noticed from subfigures (a) to (c) in Fig. 6 that pain patients after receiving neurofeedback training are mainly classified into healthy(AB) or no-pain groups showing the good performance of neurofeedback training. Further, it can be noticed that accuracy of classification of ANF-EEG into no-pain or healthy group is the highest when analysis is performed with electrode combination 1 (C1) i.e. with 61 electrodes placed on the scalp.

Fig. 6
figure 6

The prediction accuracy (Mean ± STD) of SVM based classifier applied for detecting positivepredictive values in EEG data obtained from PWP after receiving neurofeedback training; a The positivepredictive values are (AB) for AB vs PWP; b positive predictive values (PNP) for PWP vs PNP subsets;c positive predictive values (AB) for AB vs PNP subsets

4 Discussion

This study was designed to investigate the possibility of using ML to develop automatic method for evaluating the efficacy of neurofeedback training applied for treatment of CNP. Performance of classifiers applied to classify MI-EEG patterns associated with CNP therefore was explored in order to select the classifier gained highest accuracy to be used for achieving the purpose of this study.

The presence of CNP is positively associated with increasing in the power of β, Ɵ bands [22] and a shift of dominant α frequency towards lower range [28]. Therefore, neurofeedback was applied in this study to increase α band power, and suppress Ɵ and β power [16]. The efficacy evaluation of neurofeedback applied for treatment of CNP is traditionally achieved by experts based on the visual interpretation. The visual comparison of pre- and post-neurofeedback as shown in Fig. 2 is not quite reliable method as it is subject to human error, as well as, it is time consuming and costly method. Therefore, we carried out this study to investigate the possibility of developing automatic method, based on ML, to evaluate the efficacy of neurofeedback applied for treatment of CNP.

ML techniques seem to have a possibility to evaluate the effectiveness of neurofeedback based treatments since neurofeedback training produces changes in the power of a subject’s brain waves [16]. The most challenging task in this regard is to obtain features distinguishing the presence of CNP in EEG record in order to be used as a training data for the classifiers. However, MI of a limb can produce more discriminative features associated with the presence of CNP as it was shown in the previous study [21] that CNP is associated with significant increased event-related desynchronization (RED) in the theta, alpha, and beta bands during MI of a limb. Therefore, the participants, in the first and second phase of this study, were instructed to perform three types of MI in order to find out most discriminative features which may result in developing more accurate classifier.

The possibility of applying ML for evaluating the efficacy of neurofeedback therapy was investigated in this study by developing BCI classifier based on the features extracted from pre-neurofeedback EEG records. The performance of the developed classifier was tested by features extracted from post-neurofeedback EEG records. Post-neurofeedback EEG records were visually interpreted by experts to make sure that neurofeedback training resulted in significant reduction of CNP.

From a practical point of view, classification accuracy of EEG data greatly depends on the discriminative relevant features extracted from EEG recording representing complex and multiple electrical activities of the brain, and employed classifier also significantly influences the classification accuracy. Furthermore, some research studies have shown that irrelevant EEG channels may add redundant information that result in reducing classification accuracy of EEG [25, 39]. Therefore, in this study, twelve features were selected based on previous study [32] which has shown that presence of CNP produces changes of brainwaves in three bands of frequency; they are the theta (4–8 Hz), alpha (8–12 Hz), and beta (16–24 Hz), and these changes were time-locked in four periods; they are 0.4–0.8 s, 0.8–1.2 s, 1.2–1.6 s and 1.6–2.4 s. Furthermore, in an attempt to detect irrelevant channels to be avoided, we carried out process of classification with many different combinations of channels (electrodes) and as a result, we found that performance of classifiers were higher with three combinations (C1, C2 and C3) as illustrated in Fig. 3. The all combinations were placed on the scalp according to standard 10–10 location where C1 included 61 electrodes which were chosen based on studies reported global effect of pain while C2 was chosen on the basis of definition of neuropathic pain reported lesion to the somatosensory cortex. The combination C3 were chosen on the basis of somatotopoical representation of upper and lower limbs and cortical shift of paralyzed limbs reported in the literature.

DT, SVM, KNN, NB and ANN based classifiers were applied to classify extracted features in this study. The classifiers were selected on the basis of many different studies in which these classifiers were successfully applied for EEG classification [25, 40,41,42,43,43]. Moreover, these classifiers are very powerful and fast learning algorithms [23]. In order to find the most appropriate classifier for classifying MI-EEG patterns related to CNP in this study, performances of different classifiers were compared with each other.

Classification results as demonstrated in Fig. 4 shows how the classification accuracy changes with the classifier, electrodes combination and type of participants groups. It can be noticed that for all electrodes combinations (C1 to C3) (Fig. 3) and for all three types of MI tasks, the highest accuracy is mainly obtained with SVM classifier. The reason behind high accuracy gained by SVM can be attributed to mathematical model of SVM algorithm which analyzes data depending on constructing separating hyper-planes which maximize the margin between the input data classes. It can also be noticed that among all electrode combinations the combination 1 (C1) showed highest accuracy followed by C2 and C3 respectively because size of training data changes with number of channels and it was proven that number of EEG channels significantly influences classification accuracy [25]. The EEG channels are more, the classification accuracy is higher. Furthermore, comparison of sensitivity and specificity of SVM as illustrated in Fig. 5 shows that SVM classifier’s ability to correctly predict negative cases (PWP) is higher than its ability to predict positive cases (AB, PNP) because extracted features from PWP are more discriminative than features of AB or PNP. These outcomes therefore emphasize the results of previous studies [32, 44, 45] which have shown that the presence of CNP itself improves the performance of MI based classifier.

Since SVM based classifier in this study gained the highest accuracy with all combinations of electrodes, with all groups (AB, PWP and PNP) and with all classes (RH, F and LH), it was the selected classifier for evaluating efficacy of neurofeedback training applied for treatment of CNP.

Results of SVM based evaluation of CNP neurofeedback therapy as shown in Fig. 6 demonstrate that SVM classifier may effectively predict outcomes of neurofeedback applied for treatment of CNP. Performance of SVM based classifier significantly changes with electrodes combinations, participants groups, and type of pre-treatment MI-EEG patterns used for developing SVM based classifier. Electrodes combinations significantly influence the performance of SVM classifier since electrodes number in a combination determines the size of training data; and from a practical view of point, fairly increasing training data size improved performance of SVM classifier, as well as, electrodes locations on the scalp decide captured ERD response for MI[45].

Participants groups and type of MI-EEG patterns associated with CNP (BNF) also noticeably influenced performance of SVM based classifier as illustrated in Fig. 6 showing that the highest percentage of positive predicted values were obtained with RH followed by LH and F classes respectively, and with AB versus PWP followed by PWP vs PNP and AB versus PNP subsets respectively. The influence of MI-EEG subset over the classifier accuracy can be attributed to different ERD response of different subsets. Consequently, it was found that features extracted from PWP subsets are more discriminative. Moreover, the largest difference in classification accuracy among three groups (PWP, PNP and AB) was noticed with RH class. Therefore, BNF-EEG obtained from AB vs PWP of all three classes especially RH class can be used to develop effective SVM based classifier for automatically evaluating the efficacy of neurofeedback training applied for treatment of CNP.

5 Conclusion

This research study was the first study designed based on ML for developing automatic method for evaluating the efficacy of neurofeedback training applied for treatment of CNP. The results of this study support our hypothesis which states that ANF- EEG data (EEG obtained from PWP participants after receiving neurofeedback) should be classified in PNP or AB group when provided as an input to the developed BCI-classifier. The classification accuracy of ML algorithms applied in this study, changes with participants groups (training data), type of MI task, classifier, and electrodes combinations. However, highest classification accuracy was achieved with SVM classifier, (AB vs PWP) group, MI-EEG of RH, and 61 electrodes combinations. Study findings show that ML can be used effectively to validate neurofeedback training protocols. We therefore believe that ML approach can be applied to automatically evaluate the performance of neurofeedback training protocols designed for treatment of CNP. Additionally, we believe that this approach can be effectively applied in developing MI based BCIs dedicated for helping paraplegic patients to live their lives independently. However, we suggest increasing the diversity of BNF-EEG data as well as the size and diversity of ANF-EEG data to get BCI-classifier with more reliable performance.