participants
In this study, 140 depression patients were recruited with permission and informed consent. All subjects had no other neurological disorders (stroke, brain tumour, severe concussion, migraine, etc.), cardiovascular (myocardial infarction, arrhythmia, etc.) diseases, and no obvious impairment in vision and hearing. The experiment was conducted under the supervision of professional psychologists. No identifiable data of the participants were recorded during the whole study. Based on the HAMD [2] and consultation with a senior psychiatrist, the severity of depression of participants was diagnosed. On the psychologist's advice, 58 subjects with scores were labelled as mild depression, while 82 subjects with scores above were marked as severe depression. More details of the dataset are shown in Table 1.
Table 1
Subject information of our data
|
Severe depression
|
Mild depression
|
Number
|
82
|
58
|
Gender
|
54 females, 28 males
|
39 females, 19 males
|
Age(years)
|
32.48\(\pm\)17.71
|
36.65\(\pm\)19.97
|
fNIRS data acquisition
In the experiments, fNIRS produced by Hitachi is used to obtain data. It emits frequency-conversion near-infrared rays (wavelength: \(659\pm 20\)nm, \(830\pm 20\)nm) to penetrate the scalp, skull, and cerebrospinal fluid to irradiate the cerebral cortex. The device has 52 pairs of near-infrared emitters and receivers, forming 52 channels, including 18 in the prefrontal lobe and 17 in the left and right frontal lobes. During the verbal fluency task (VFT), neural activity leads to increased CBF, causing changes in brain oxyhemoglobin (HbO) and deoxygenated haemoglobin (HbR) [6]. According to Lambert-Beer [18], the fNIRS device converts the measured change in near-infrared rays into the change of the HbO and HbR. During the VFT, subjects wore the near-infrared device on their heads and performed three tests within 60 seconds, with 30 seconds of silence before and after the tasks. The task lasted 120 seconds, with a sampling frequency of 10 Hz, and each subject collected \(120\times 10\times 52\times 2\) sample data. Figure 1 shows the experimental data collection process.
Methods
To predict the severity of depression, we developed a prediction model for mild depression and severe depression. The model involves data preprocessing, feature extraction, feature selection, feature fusion, and classification prediction, as illustrated in Fig. 2.
Data pre-processing
We derive the change in concentration of HbO and HbR calculated by the Beer-Lambert [19] by fNIRS device, and these hemodynamic parameters are considered sensitive indicators for studying psychiatric disorders [20]. To remove physiological oscillations (such as respiration(0.2–0.4 Hz) and heartbeat(0.5-2.0 Hz)), we apply a third-order Chebyshev type II filter [21] with a cut-off frequency of 0.01 Hz and a stop frequency of 0.2 Hz. However, this filtering cannot eliminate jitter noise and motion-related artefacts during signal acquisition. Empirical Mode Decomposition (EMD) and Ensemble Empirical Mode Decomposition (EEMD) methods have been used in some studies to eliminate motion artefacts, and baseline drifts in physiological signals [22, 23]. In this paper, a complete ensemble empirical mode decomposition with adaptive noise-wavelet threshold combined denoising method (CEEMDAN-WPT) is utilized to eliminate noise in fNIRS data and achieve good results. The process of the CEEMDAN-WPT algorithm is outlined in Fig. 3, and the specific implementation steps are as follows.
1. Suppose
\({E}_{i}(\bullet )\) is the
\(i\)-th IMF component operation operation calculated by the EMD algorithm [n]. Adding white noise that conforms to the standard distribution to the original signal yields
\(I\) initial signal
\({x}^{i}\left(t\right)\):
$$\begin{array}{c}{x}^{i}\left(t\right)=x\left(t\right)+{\epsilon }_{i}{\omega }^{i}\left(t\right) \left(1\right)\end{array}$$
where \(x\left(t\right)\) is the original signal, \({\epsilon }_{i}\) is the white noise amplification factor, and \({\omega }^{i}\left(t\right)\) is the white noise that conforms to the standard Gaussian distribution.
2. EMD decomposition is used for all initial signals to obtain the first modal component:
$$\begin{array}{c}\stackrel{-}{IM{F}_{1}\left(t\right)}=\frac{1}{I}\sum _{i=1}^{I}IM{F}_{1}^{i}\left(t\right) \left(2\right)\end{array}$$
Then calculate the remaining components of the first order:
$$\begin{array}{c}{r}_{1}\left(t\right)=x\left(t\right)-\stackrel{-}{IM{F}_{1}\left(t\right)} \left(3\right)\end{array}$$
3. Add \({E}_{1}\left({\omega }^{i}\right(n\left)\right)\) to \({r}_{1}\left(t\right)\) to get the new signal \({r}_{1}\left(t\right)+{\epsilon }_{1}{E}_{1}\left({\omega }^{i}\right(n\left)\right)\), perform EMD decomposition again to calculate the second residual component:
$$\begin{array}{c}\stackrel{-}{IM{F}_{2}\left(t\right)}=\frac{1}{I}\sum _{i=1}^{I}{E}_{1}\left({r}_{1}\left(t\right)\right)+{\epsilon }_{1}{E}_{1}\left({\omega }^{i}\left(n\right)\right) \left(4\right)\end{array}$$
Then calculate the residual components of the second order:
$$\begin{array}{c}{r}_{2}\left(t\right)={r}_{1}\left(t\right)-\stackrel{-}{IM{F}_{2}\left(t\right)} \left(5\right)\end{array}$$
4. Repeat the above steps, and when the number of poles of the obtained residual signal is less than two, the algorithm ends decomposition. At this point, the original signal \(x\left(t\right)\) is decomposed into:
$$\begin{array}{c}x\left(t\right)=\sum _{k=1}^{K}\stackrel{-}{IM{F}_{k}}+R\left(t\right) \left(6\right)\end{array}$$
5. Carry out wavelet threshold denoising for the noisy high-frequency component \(IM{F}_{k}\):
$$\begin{array}{c}\left\{\begin{array}{c}y=W\left(IM{F}_{k}\right)\\ \tilde{y}=\tilde{D}\left(y,\lambda \right)\\ \stackrel{-}{IM{F}_{k}}=\stackrel{-}{W}\left(\tilde{y}\right)\end{array}\right. \left(7\right)\end{array}$$
In the formula, \(W\) and \(\stackrel{-}{W}\) are respectively the wavelet packet change and its inverse operation. \(\tilde{D}\) and \(\lambda\) are the threshold function and threshold magnitude, respectively, and \(\stackrel{-}{IM{F}_{k}}\) is the high-frequency component of denoising processing.
6. The denoising high-frequency component \(\stackrel{-}{IM{F}_{k}}\) and the low-frequency component \(IMF\) are reconstructed to obtain the denoising signal \({x}^{{\prime }}\left(t\right)\).
Feature Extraction
a. Time Feature
Considering that the actual task time is only 60 seconds during the VFT, we extract the temporal features (TF), including time domain, frequency domain features, and information entropy during this task time. Time domain features consist of maximum, minimum, mean, rectification average, skewness, and peak factors. Mean squared frequency is the selected frequency domain feature. Power spectral entropy and singular spectral entropy are the information entropy features, calculated as follows:
$$\begin{array}{c}{H}_{f}=-\sum _{k=0}^{k-1}({S}_{k}\sum _{k=1}^{N}{S}_{k})\text{log}\left({S}_{k}\sum _{k=1}^{N}{S}_{k}\right) \left(8\right)\end{array}$$
In the formula, \({H}_{f}\)represents the power spectral entropy, and \({S}_{k}\) represents the energy division of the signal in the frequency domain.
$$\begin{array}{c}{H}_{t}=-\sum _{i=1}^{l}\left(\frac{{\lambda }_{i}}{\sum _{i=1}^{l}{\lambda }_{i}}\right)\text{log}\left(\frac{{\lambda }_{i}}{\sum _{i=1}^{l}{\lambda }_{i}}\right) \left(9\right)\end{array}$$
Where, \({H}_{t}\) is the singular spectral entropy and \({\lambda }_{i}\) is the singular value spectrum obtained by the singular value decomposition of the signal.
b. Correlation Feature
In our study, the correlation among channels was extracted as the correlation features (CF), which was obtained by calculating the Pearson correlation coefficient [11] of the change of HbO and HbR concentration between 18 channels in the prefrontal lobe of each subject. Its calculation formula is as follows:
$$\begin{array}{c}Cor\left(X,Y\right)=\frac{\sum XY-\frac{\sum X\bullet \sum Y}{N}}{\sqrt{\left(\sum {X}^{2}-\frac{{\left(\sum X\right)}^{2}}{N}\right)\left(\sum {Y}^{2}-\frac{{\left(\sum Y\right)}^{2}}{N}\right)}} \left(10\right)\end{array}$$
where \(X\) and \(Y\) represent the change of hemoglobin concentration in different channels of the subject, respectively.
Therefore, a total of \(18\times 2\times 9=324\) TF and \({C}_{18}^{2}\times 2=306\) CF can be extracted from the HbO and HbR concentration change curves of 18 channels in the prefrontal lobe of each subject.
Feature Selection
In the classification task, when the amount of sample data is small and the number of features is large, it will cause a disaster of dimensionality and reduce the accuracy. The recursive feature elimination with cross-validation (RFECV) method [24] is utilized to select the optimum features from TF and CF. The RFECV method can obtain the importance of each feature through RFE [25] and find the optimal number of features based on cross-validation accuracy. Compared with dimension reduction methods such as PCA [26], the RFECV method can exclude redundant feature information and determine the most impactful features for achieving optimal classification performance. In our prediction model, we used the REFCV method twice. First, the optimal 47 TF and 54 CF were screened separately. Second, after fusing the optimal CF and TF selected in the first step, the best 35 features for discriminating mild and severe depression were screened again.
Classification models
Given the limited number of depression data samples, we utilized machine learning to construct the classification prediction. Four classical algorithms are utilized to build categorical predictive models, including two supervised learning methods (logistic regression (LR) [27] and support vector machines (SVM) [28]), an ensemble learning algorithm (random forest (RF) [29]), and artificial neural networks (Multi-layer perceptron (MLP) [30]). The experiment separated the datasets into two parts: the 60% training set and the 40% testing set. In order to find the best prediction model for each mode, the grid search [31] based on hyperparameter turning [32] was selected to determine the model parameter with the highest accuracy. During model training, validation testing was done using five-fold cross-validation to avoid over-fitting.
Model metrics
The signal-to-noise ratio (SNR) and root mean square error (RMSE) are introduced to evaluate the noise reduction effect of the CEEMDAN-WPT algorithm. For the classification models, we calculated the area under the ROC curve (AUC) of different models in contrast to their performance differences. Several other indicators included specificity, sensitivity, accuracy, and F1 score.