In this study, the auditory oddball task EEG data were collected from healthy subjects (HC) and patients with DD and CPP. First, using the PLI coupling to construct the brain functional connectivity network, the connectivity differences between groups of different brain regions in different frequency bands are compared and analyzed using graph theory and CNN. Then, the CNN algorithm is used to classify diseases based on EEG feature extraction in different frequency bands, providing a new method for objective evaluation of CPP. The processing flow of relevant data is shown in Figure 1.

In this study, twenty-two patients with CPP (13 females, 46 ± 9.62 years), twenty-two patients with DD (11 females, 46 ± 12.33 years), and twenty-three HCs (13 females, 49 ± 10.26 years) were recruited from the Shanghai Changhai Hospital, Shanghai Hongkou Mental Health Center and University of Shanghai, respectively, (Shanghai, China). All participants gave their written informed consent after the experimental procedure had been carefully explained. The research has been approved by the local research ethics committee.

Inclusion and exclusion criteria of subjects: chronic primary pain concerning the diagnostic criteria proposed by ICD-11: pain in one or more anatomical regions that (1) persists or recurs for longer than 3 months; (2) associated with significant emotional distress or functional disability (interference with activities of daily life and participation in social roles); (3) cannot be better accounted for by another chronic pain condition (American Psychiatric Association, 2013). Depressive disorder diagnosis was established according to the Diagnostic and Statistical Manual of Mental Disorders, 5th edition (DSM-V) criteria, as assessed by the structured clinical interview for DSM-V and Hamilton Depression Scale (HAMD). The two categories of patients also need to meet the following inclusion criteria: (1) Age 18–65 years old; (2) Junior high school degree or above; (3) Not taking medicine at the initial diagnosis or more than 6 months after drug withdrawal; (4) Volunteer to participate in this study and sign the informed consent form. Exclusion criteria: (1) Serious cognitive impairment or hearing disability, unable to cooperate in the completion of project-related assessment tests; (2) Serious physical diseases; (3) Severe psychiatric symptoms; (4) Abuse of psychoactive substances; (5) Suffering from diseases that can cause secondary chronic pain.

All subjects were instructed to fill in SDS for self-evaluation of depression level (Table 1). Two groups of patients were assessed by the physician with HAMD. The subjective and multidimensional experience of pain in CPP patients was quantitatively measured using the Short-Form McGill Pain Questionnaire (Dworkin et al., 2009). It comprises three subscales: a pain rating index (PRI) describing the qualities of pain, a 10 cm visual analog scale (VAS) describing the intensity of averaged daily pain during the past 2 weeks, and a present pain intensity (PPI) index describing the intensity of current pain.

Auditory oddball task was applied. The stimulus materials for the auditory oddball task were composed of 30 target stimuli (1,200 Hz, 75 dB, 50 ms) and 200 standard stimuli (1,000 Hz, 75 dB, 50 ms). The stimuli were arranged in a pseudo-random order with an interval of 1,000–1,500 ms between each stimulus. Before the experiment, there was a training part consisting of 10 tones, including two target stimuli. Subjects were asked to identify the tones with a low probability of occurrence, count them in silence, and then report to the researchers. The experimental program was written by E-Prime 1.0 software of Psychology Software Tools in the United States. The presentation mode of the stimulus sequence is shown in Figure 2.

The experimental data acquisition equipment is 32 channels high-density EEG acquisition equipment made by Brain Productions Company in Germany. The recording software is based on the Vision Recorder system developed by the above company. During the whole experiment, the electrode was within the range below 30 k Ω, and the sampling rate was 1,000 Hz. The EEG cap was set with a 0.1–100 Hz band-pass filter in the default setting of the headset. To reduce the calculation cost during data processing, the data was down-sampled to 250 Hz. The impedance of recording electrodes was kept at less than 5 k Ω. The collected data were preprocessed using the EEGLAB toolbox (v2019.1) of MATLAB (R2020a) software: select 23 electrodes were placed on the scalp based on the international 10–20 electrode placement system and take Cz as the reference electrode. See Figure 1 for the selected electrodes and functional area division (Delorme and Makeig, 2004). Take the target stimulus as 0 times, select −0.2 s∼0.8 s data as a trial, and remove the components of eye and muscle electricity through independent component analysis.

As Figure 3 shows: (1) First, EEG data is collected and preprocessed; (2) To establish the functional connection, it is necessary to obtain the time series of brain activities in different regions; (3) According to these sequences, aggregation measurement is used to calculate the correlation of these sequences, so that the brain network is represented as a correlation matrix; (4) The consistency threshold method is applied to the matrix (30% of the connections are reserved) to retain the most significant part of the features so that the connection strength in different brain regions can be visually expressed; (5) Finally, the graph theory features are calculated and analyzed.

Before calculating PLI, the analytic signal based on the Hilbert transform is used to determine the instantaneous phase, to calculate the phase synchronization. The calculation formula of PLI is (Stam et al., 2007):

Where < > represents the expected value, ∅_x(f)−∅_y(f) represents the phase synchronization of x and y channels/brain region signals at f frequency. The value range of PLI is [0,1], where 1 represents complete phase synchronization and 0 represents no phase synchronization.

Graph theory is a method that can be applied to brain networks. It can describe the topological structure of complex networks and the changes in different network metrics in networks. The diagram is composed of a group of nodes (electrode array) and their connections (edges). For the quantification of graph topology, there are many measures. In this study, the following indicators were selected to conduct in-depth research on the constructed functional brain network.

This paper uses a compact CNN that can be used for various EEG signal classification tasks, including event-related potentials (ERP) (Prabhu et al., 2020). In addition, this article improves on the CNN and reduces the first AvgPool2d layer kernel size and stride by half, enabling CNN to classify the PLI matrix (22 × 22). Its architecture is shown in Table 2. The input of this model is raw EEG data, including channel number, and time sample. To limit the number of trainable parameters, the architecture adopts depth and separable convolution. The initial combination of 2D convolution and depth convolution allows each temporal filter to learn spatial filters (Schirrmeister et al., 2017). At the same time, the number of spatial filters learned from each feature map is controlled by the depth parameter. After each convolution, batch normalization is performed to achieve model stability. In addition, the dropout layer is used to significantly reduce overfitting (Mukhtar et al., 2021). The final multi-category classification layer uses the SoftMax function.

When using CNN for training, this study uses 10-fold cross-validation processing method, in which one random fold is used as the test set, the other twofolds are used as the validation set, and the rest as the training set. It also ensures that three types of subjects are evenly distributed in each fold and that the data of each subject only exists in the same fold. At the same time, each dataset (training set, validation set, and test set) is disordered before the training stage (Figure 4).

The network is trained using the backpropagation algorithm (Hung and Adeli, 1993) with a batch size (the number of training samples in iteration) of 32. An optimization algorithm, namely adaptive moment estimation (Adam) (Mukhtar et al., 2021) is adopted in this work to update the parameter of the proposed network structure. It was observed it enables the network to converge at a faster rate thereby improving the efficiency of the training process.

The following equation is used to update the 1st-moment estimate (Kingma and Ba, 2014). All operations on vectors are element by element.

Where m, t, c, θ_t and ∂⁡c∂⁡θ are defined as the 1st-moment vector, timestep, cost function, resulting parameters (weights), and gradient, respectively. The parameters β₁ and β₂ represent the exponential decay rates which are chosen to be 0.9 and 0.999, respectively, (Kingma and Ba, 2014).

The following equation is used to update the 2nd-moment estimate (Kingma and Ba, 2014).

The following equations are used to compute the 1st and 2nd-moment estimates, respectively, (Kingma and Ba, 2014).

In addition, the following equation is used to update the weights of links connecting the layers (Kingma and Ba, 2014).

Where α and ϵ denote the learning rate and numerical value and are set at 1×10⁻⁴ and10⁻⁸, respectively, in this work. The variable λ, known as the regularization parameter, is also one of the essential parameters during training to prevent data overfitting. It is tuned to 0.2 in this work.

To avoid overfitting and improve the generalization, the dropout (Srivastava et al., 2014) technique is applied to the fully-connected layers 6 and 12. During training for each mini-batch, some of the neurons from these layers are selected randomly and dropped. This forces the model to learn from a subset of input features and not the entire input features. The value is set to 0.5, i.e., the probability of a neuron being retained during training is 50% and the probability of a neuron being rejected is 50%.

In this paper, the negative log-likelihood loss function is used to deal with multi-classification problems (de Boer et al., 2005). The input is the logarithmic probability value. For the batch dataD(x,y) containing N samples, x is the output of the neural network and is normalized and logarithm zed. y is the category label corresponding to the sample, and each sample may be one of C categories.

l_n is the loss corresponding to the nth sample, the value range is [0, C-1]:

w is used for sample imbalance between multiple categories:

The default value of reduction is mean, and the corresponding l(x,y) is:

F1 Score is an indicator used to measure the accuracy of the two-class model in statistics. It can be seen as a weighted average of the model’s accuracy and recall. Its value range is [0, 1]. The larger the value is, the better the model is.

In multi-classification problems, Macro-F1 is usually used (Supratak et al., 2017).

Where n is the number of categories and i is the number of categories.

SPSS software was used for statistical analysis. One-way analysis of variance analysis (ANOVA) was used for the comparative analysis of brain networks of three groups, where group category was used as a factor, and the average value of each channel after the superposition of the brain network connection matrix was used as a dependent variable. The Tukey method was used as a post hoc analysis for testing differences between groups.

One-way ANOVA found that the brain networks of the three groups are statistically different except for the Beta band (Table 3). Further post hoc analysis based on the Tukey method (Table 4) shows that in the Gamma band there was difference in network characteristics between CPP and HC subjects but not between CPP and DD or DD and HC. In Delta and Alpha bands, a significant difference was present between CPP and DD, DD and HC, but not between CPP and HC. In the Theta frequency band only, the P-values between groups were all less than 0.05, indicating that different groups possess independent brain network characteristics in this frequency band.

To further explore the specific differences in brain networks among the three groups, we carried out a functional connectivity matrix analysis of CPP, DD, and HC and compared the differences after subtraction and the consistency threshold method. See Figure 5 for Theta band results (see Supplementary material for results in other bands). The results show that in the CPP group, the Theta frequency band connection between parieto-occipital electrodes (P3, O2) and prefrontal electrodes (FP1, FP2) is significantly enhanced (the connection strength of 4–5); in the HC group, the connections between electrode Cp5 and O2, electrode FC2 and F4 are slightly enhanced (the connection strength of 3–4); in DD group, the connection strength is the weakest (the connection strength of 1–3). Further difference comparison shows that the CPP group has significantly increased connectivity between the prefrontal area and the posterior parietal area compared with the DD group, while the difference between CPP and HC groups is mainly reflected in the stronger connectivity of CP5-O2 and F4-FC2; Compared to HC group, the connectivity strength of brain network in DD subjects is inhibited, mainly in the form of weakened connectivity between the parieto-occipital area and the prefrontal area.

We further analyzed and compared the global feature revealed by E_glob, as shown in Figure 6. The E_glob of DD subjects in Delta and Alpha bands is higher than the other two groups of subjects, while it is inhibited in the Theta band. This is also verified by the classification results of the PLI matrix using CNN (see Section “3.3. Classification accuracy of different frequency bands”). The difference between CPP and HC subjects is mainly reflected in the significant enhancement of the Theta band.

At the node level, DC, BC, E_loc, and CC of the three groups in the Theta band is also different, as shown in Figure 7. The DC intensity of CPP subjects is higher in the right frontal and left parieto-occipital region than in the other two groups. This result is consistent with the aforementioned trend of brain network connection intensity in the Theta band. In the distribution of BC, the intensity in the central area of CPP subjects is lower than that of the other two groups, and the intensity of BC in the prefrontal and temporal areas of CPP and DD is lower than that of HC subjects. Similarly, the E_loc intensity in the central area of CPP subjects is also lower than that of the other two groups, and the overall E_loc intensity of DD subjects is lower than that of the other two groups. In the distribution of CC, the overall distribution of the three groups is similar, but the left frontal lobe intensity is higher in the CPP and DD groups.

CNN was first used to classify the EEG of three groups of subjects. Further, this paper used the improved CNN model to classify the PLI matrix of three groups of subjects to verify the above brain network analysis results. After 1,000 training epochs, the loss rate, and verification loss rate of the two classifications began to converge to 0 and 0.8, respectively, while the training accuracy and verification accuracy converged to 1 and 0.3, respectively, which proved that the network structure was stable (Figure 8). Figure 9 shows that both sets of classifications have achieved ideal accuracy in six frequency bands: in the classification of EEG, the Gamma band has the highest accuracy of 85.01%, followed by the Theta band with 81.37%; in the classification of PLI, the highest frequency band of Theta is 79.64%, followed by Delta and Gamma with 77.81 and 77.47%. The results of Macro-F1 in Figure 10 show that the classification model of EEG is better than that of PLI. Except for the Beta frequency band, the Macro-F1 evaluation of other frequency bands are above 0.65, and the PLI classification model is also above 0.6 except for the Beta frequency band.