A Habenula Neural Biomarker Simultaneously Tracks Weekly and Daily Symptom Variations during Deep Brain Stimulation Therapy for Depression

2.1 Participant

This study included a male TRD patient aged 36-40 years old (see Note S1 for detailed patient medical information) participating in a clinical trial of LHb DBS treatment starting in October 2021. The patient provided informed consent for participation in the clinical trial. This study received approval from the Ethics Committee of Zhejiang University School of Medicine Second Affiliated Hospital (protocol number 20210218). It was registered at www.clinicaltrials.gov (NCT05716555), where detailed information regarding the inclusion and exclusion criteria can be accessed.

At the beginning of the clinical trial, two independent psychiatrists evaluated the patient’s psychotic symptoms using the 17-item Hamilton Depression Rating Scale (HAMD), the Montgomery Asberg Depression Scale (MADRS), and the Hamilton Anxiety Rating Scale (HAMA) as baseline assessments. In addition to the psychiatric assessments, the patient underwent a comprehensive physical examination, various mental scale assessments, and a magnetic resonance imaging (MRI) examination. We carefully ensured that other psychiatric diagnoses outlined in the Diagnostic and Statistical Manual of Mental Disorders-Fifth Edition (DSM-5) were excluded.

2.2 Surgical procedure

A standard DBS implantation procedure was employed. Bilateral quadripolar electrodes (1200-40, SceneRay, Suzhou, China) were surgically implanted in the LHb under local anesthesia (Figure 1B). The DBS electrodes had a diameter of 1.27 mm and a lead length of 400 mm. Each electrode’s four contacts measured 1.5 mm in length with a spacing of 0.5 mm. The LHb targeting was guided by preoperative MRI sequences. After confirming the absence of stimulation side effects through intraoperative testing, an implantable pulse generator (SR1101, SceneRay) was placed under general anesthesia. The DBS device was also capable of recording and wireless transmitting LFP signals (Figure 1C).

2.3 DBS treatment process and symptom evaluations

During the bilateral DBS treatment process, we made multiple adjustments to the stimulation parameters to achieve the best therapeutic effect. We divided the treatment process into six stages based on the alterations of stimulation parameters (Figure 1D): 1) the “Preop” stage, the time before the DBS electrode implantation; 2) the “Off-1” stage, patient recovery with DBS turned off; 3) the “1 Hz” stage, activation of 1 Hz stimulation; 4) the “Off-2” stage, DBS turned off because of unnoticed power off; 5) the “20 Hz” stage, re-activation of 20 Hz stimulation; 6) the “130 Hz” stage, activation of 130 Hz stimulation. More details can be found in Note S2. The entire duration of DBS treatment spanned 41 weeks (starting from DBS electrode implantation).

The efficacy of DBS treatment was evaluated from two perspectives: clinician evaluation and self evaluation. For clinician evaluation, a psychologist blinded to the current stimulation parameters and their adjustments evaluated the patient’s depression and anxiety symptoms on a weekly basis using standardized rating scales (HAMD, MADRS, HAMA). Response is defined as a 50% or greater improvement on the HAMD score from the pre-treatment baseline. Remission is defined as achieving a HAMD score of 7 or less. The psychologist also evaluated the patient’s emotional blunting and cognitive functioning during the treatment (see Note S3). For self evaluation, the patient used the Visual Analogue Scale (VAS) for depression (VAS-D) and anxiety (VAS-A) to self-report the symptom severity. Self-reported VASA and VAS-D had been used to assess the rapid effects of antidepressants [48]. To facilitate daily data collection, we established an online questionnaire system where the patient could conveniently complete the self-reports via the smartphone or computer.

2.4 LFP signal recording, signal processing, and feature extraction

After activating the 1 Hz stimulation, we collected daily LFP signals (30 minutes per day) concurrent with daily self-reported VAS-D and VAS-A (details in section 2.3). LFP signals were recorded at a sampling rate of 1000 Hz. Notably, stimulation was deactivated during the signal acquisition process. We reconstructed the electrode positions using MRI and selected two contacts in the left hemisphere for bipolar recording of a single LFP channel. The patient was instructed to attempt daily LFP recording and VAS-D/VAS-A reporting. Throughout the entire 41-week (287-day) treatment, the patient was able to activate LFP recording and report VAS-D and VAS-A on 122 days distributed across 26 weeks. Therefore, the subsequent analyses focused on the LFP signals, VAS-D, and VAS-A scores recorded from these 122 days, and the HAMD, MADRS, and HAMA scores recorded for the 26 weeks.

Custom MATLAB scripts (MathWorks Inc., Natick, MA, USA) were used to preprocess the LFP signals. The LFP signals were first band-pass filtered from 1 to 30 Hz using a Butterworth filter of order 12 to avoid the noise observed in higher frequency bands. Then, we divided the daily 30-minute LFP signals into 10-second epochs with a 50% overlap. Next, we used a standard procedure (details in Note S4) to remove bad epochs from daily LFP signals (example temporal traces and spectrum of preprocessed LFP epochs were shown in Figure 1C).

For each remaining LFP epoch, we computed its spectral domain (SD) and temporal domain (TD) features. SD features included PSD of the four bands (δ (1-4 Hz), θ (4-8 Hz), α (8-12 Hz), and β (12-30 Hz)) and phase-amplitude coupling (PAC) for six specific pairs of coupling. TD features included fourteen temporal domain features used in previous study [49], e.g., Hjorth mobility, singular value decomposition (SVD) Fisher information, Hurst exponent, etc. These features capture the temporal properties of LFP from probabilistic distribution and information theory perspectives and have been widely used in brain signal analyses [50,51]. As a result, we obtained 24 features, comprising 10 SD features and 14 TD features for each LFP epoch. Details of these 24 features are included in Table S1 and Note S4. Finally, we averaged each feature across LFP epochs within the same day and obtained a single averaged 24-dimensional LFP feature vector. Our subsequent analyses were based on the daily LFP features as computed above.

2.5 Identification of neural biomarker

We first conducted Spearman’s rank correlation analyses between LFP features and symptoms. For each day, we correlated each daily LFP feature with the daily VAS-D and VAS-A self-reports. For each week, we computed the average of LFP features across the days that belonged to this week, resulting in weekly LFP features; we then correlated each weekly LFP feature with the weekly clinical evaluation scales HAMD, HAMA, and MADRS. Bonferroni correction was used to adjust for multiple comparisons.

Next, we used a data-driven method to identify an LHb neural biomarker of depression symptoms, where we built a machine learning model to use LFP features to classify high and low depression symptom states.

First, we defined the high and low depression symptom states of the patient by k-means clustering the weekly depression scales HAMD and MADRS similar to prior work [17]. Among the total 26 weeks (122 days) of LFP data, 7 weeks (29 days) of LFP data belonged to the low depression symptom state (labeled 0), 4 weeks (22 days) of LFP data belonged to the high depression symptom state (labeled 1). The remaining 15 weeks (71 days) were unlabeled and used as test data for subsequent biomarker tracking evaluation (see next section).

Second, based on the labeled data, we built a machine learning model to use the LFP features to classify high and low depression symptom states. We constructed six machine learning models: logistic regression (LR), multilayer perceptron (MLP), adaptive boosting (AdaBoost), support vector machine (SVM), random forest (RF), and linear discriminant analysis (LDA). We trained and tested these models using 5-fold cross-validations that were repeated 200 times, where we computed the averaged cross-validated classification accuracy, specificity, sensitivity, F1 score, and Receiver Operating Characteristic (ROC) Area Under the Curve (AUC) score as the performance metrics. The model with the highest accuracy was selected for further analysis.

Third, the chosen model was retrained with all labeled data, leading to a “neural biomarker model”. This model takes the LFP feature as input and outputs the decision variable as the neural biomarker value (e.g., in the LR model, the decision variable was computed from the decision probability via the inverse sigmoid function). This allows us to compute a neural biomarker value for any given LFP feature. Higher neural biomarker values indicate more severe depression symptoms.

In essence, our identified neural biomarker aggregates spectral and temporal domain features from the LHb LFP signal to classify high and low depression states during DBS treatment.

2.6 Evaluation of the neural biomarker

We evaluated the identified neural biomarker in terms of (1) tracking the temporal dynamics of weekly symptom variations; (2) tracking the temporal dynamics of daily symptom variations; (3) reflecting changes in symptom variations induced by DBS frequency alterations.

First, we investigated tracking the temporal dynamics of weekly depression and anxiety symptom scales that were not used in neural biomarker identification. We took the daily LFP features as inputs to the neural biomarker model and computed the output daily neural biomarkers. We then averaged the daily neural biomarkers belonging to the same week to compute the weekly neural biomarkers. We next correlated the weekly neural biomarker values with the weekly HAMD, MADRS, and HAMA scores, respectively, using Spearman’s rank correlation analysis with explained variance (EV) as an estimation. We further analyzed the temporal dynamics in the neural biomarker and symptoms, using the dynamic time warping (DTW) distance [52] to measure the temporal tracking ability of the neural biomarker. Both the neural biomarker values and the symptom scales were normalized to a range of 0 to 1. We used a size three Sakoe–Chiba warping window in the DTW analysis following prior work [17,53]. A smaller DTW distance represents better temporal tracking. To determine the significance of the computed DTW distance, we randomly shuffled the temporal sequence of the neural biomarker 10,000 times and used the corresponding shuffled DTW distances as the null hypothesis distribution for computing the P value.

Second, we investigated tracking the temporal dynamics of the daily VAS-D and VAS-A self-reports, which were also not used in neural biomarker identification. Similar to the weekly case, daily LFP features were used to generate daily neural biomarkers, which were then correlated with daily VAS-D and VAS-A reports. DTW was again used to assess the temporal tracking of daily depression symptom variations. Third, we qualitatively compared trends in weekly neural biomarkers and depression ratings across three DBS frequency alterations (1 Hz to Off-2, Off-2 to 20 Hz, 20 Hz to 130 Hz). We used the two-sided Wilcoxon rank-sum test to check whether there was a significant difference between the two stages before and after alteration. We also averaged the neural biomarker values and depression ratings across five time periods for each case: 1) from the beginning of this stage to two weeks before the alteration week; 2) during the week before the alteration week; 3) during the alteration week; 4) during one week after the alteration week; 5) averaged from two weeks after the alteration week to the end of this stage.

3.1 LHb DBS improved the patient’s clinical symptoms, which were significantly correlated with LHb LFP features

We first examined the TRD patient’s symptom changes throughout the LHb DBS treatment process. At the beginning of treatment, the patient’s baseline HAMD score was 20, MADRS score was 25, and HAMA score was 16. In terms of the weekly clinical ratings (Figure 2A), the patient responded at week 14 (HAMD score dropped to 10; MADRS score dropped to 19; HAMA score dropped to 8) and achieved remission by the end of the 41-week treatment (HAMD score was 7; MADRS score was 9; HAMA score was 6). The daily self-reports followed a similar decreasing trend (Figure 2B). Such a consistent trend was confirmed by the strong positive correlation between the daily self-reports and weekly clinical ratings (Spearman’s ρ > 0.5, P < 0.05 for all pair-wise correlations; see Table S2 and Figure S1 for details). Besides alleviating the symptoms based on the weekly clinical ratings and daily self-reports, we also found improvement in emotional blunting and cognitive functioning (Note S3 and Tables S3-S5).

• Download figure
• Open in new tab

Figure 2.

Changes of weekly and daily symptom scores during the LHb DBS treatment and their correlations with LFP temporal and spectral domain features. (A) Changes of weekly clinical ratings during the treatment. The vertical dashed lines represent different treatment stages indicated by the x-axis labels. (B) Changes of daily self-reports during the treatment. The vertical dashed lines represent different treatment stages indicated by the x-axis labels. (C) Heatmap of the correlation coefficients. Each cell shows the correlation coefficient (CC) value between one LFP feature (y-axis) and one symptom score (x-axis), and cells marked with * indicate the coefficients that are significantly different from zero (Bonferroni corrected p<0.05). (D) Positive and negative correlation examples with the weekly HAMD score and daily VAS-D score as shown by the yellow boxes in (C).

During the LHb DBS treatment process, we recorded daily LHb LFP signals. Therefore, we investigated how the LFP features correlated with the patient’s symptom changes. We found that many of the temporal and spectral domain LHb LFP features were significantly correlated with the weekly clinical ratings and daily self-reports (Figure 2C). For example, Hurst exponent exhibited the strong correlations with both weekly HAMD scores (Figure 2D, Spearman’s ρ = −0.85, Bonferroni corrected P = 4.7 × 10^!“) and daily VAS-D scores (Figure 2D, Spearman’s ρ = −0.76, Bonferroni corrected P = 3.8 × 10^!##). These results show that LHb LFP temporal and spectral domain features were strongly correlated with weekly and daily depression symptom scores, indicating that it is feasible to identify an LHb neural biomarker of depression symptoms from the LFP temporal and spectral domain features.

3.2 Accurate classification of high and low depression symptom severity states led to the identification of an LHb neural biomarker

We next used the LFP temporal and spectral features to identify a neural biomarker that can classify high and low depression symptom severity states. We started by defining a state of high symptom severity and a state of low symptom severity via clustering the weekly depression scales HAMD and MADRS (Figure 3A). The high symptom state (7 weeks) had an average HAMD score of 12.8 and an average MADRS score of 22.5, while the low symptom state (4 weeks) had an average HAMD score of 6.3 and an average MADRS score of 9.0. We then used LFP temporal domain and spectral domain features from these 11 weeks to classify the high and low symptom severity states via six machine learning models in cross-validation. Among these six models, the LR model performed better than other more complicated models (Figure 3B and table S6). Specifically, for the LR model, the cross-validated classification accuracy was 0.973 ± 0.002 (Mean ± SEM), the specificity was 0.961 ± 0.003, the sensitivity was 0.988 ± 0.002, the F1-score was 0.970 ± 0.002, and the AUC score was 0.974 ± 0.001, which were all significantly higher than other models (Wilcoxon rank sum test, Bonferroni corrected P < 0.05 for all comparisons), suggesting that the LR model was best suited for classifying the collected data. We thus selected the LR model for further analysis. Then, we retrained the LR model using all labeled data, resulting in the neural biomarker model. The neural biomarker model takes the LFP features as input and outputs the model decision value as the identified neural biomarker, with higher values indicating worse depression symptoms.

• Download figure
• Open in new tab

Figure 3.

Accurate classification of high and low depression symptom severity states by the identified neural biomarker. (A) Clustering of depression symptom severity states. We clustered the HAMD and MADRS scores to obtain two distinct symptom states: a high depression symptom severity state (shaded in red) and a low depression symptom severity state (shaded in blue). Each point represents data of a week. The average HAMD and MADRS scores of the two clusters are also indicated in the figure. (B) Classification performance of different classifiers. The bar represents mean and the whiskers represent the 95% confidence interval. Six classification models (different colors) were compared in terms of five performance matrices (different x-axis groups): accuracy, specificity, sensitivity, F1-Score and AUC. The best model was indicated by a yellow star for each metric. Classification model names and abbreviations: logistic regression (LR), multilayer perceptron (MLP), adaptive boosting (AdaBoost), support vector machine (SVM), random forest (RF), and linear discriminant analysis (LDA). (C) Classification accuracy of individual features and all LFP features by the logistic regression model. The individual SD features are in green, and their indices are ordered based on classification accuracy. The individual TD features are in blue, and their indices are ordered by classification accuracy. The combination of all LFP features is in brown. The best individual feature was indicated by a yellow star for each domain. (D) The absolute coefficients of each feature in the logistic regression model trained with all LFP features. The individual feature indices are the same as in (C).

We further investigated how different features contributed to identifying the neural biomarker by performing separate classifications for each feature (Figure 3C). The best temporal domain feature was the Hurst exponent, with a cross-validated classification accuracy of 0.932. While the best spectral domain feature was the PSD of the β band, with an average accuracy of 0.803. Overall, temporal domain features (average accuracy: 0.781) outperformed spectral domain features (average accuracy: 0.705). Notably, combining all the LFP features yielded superior performance compared to individual features. Consistently, by investigating the logistic regression coefficients of the neural biomarker model (Figure 3D), we found that the features with better classification accuracy also had larger coefficients in the neural biomarker model. These results show that the temporal domain and spectral domain features had supplementary information that both contributed to the identification of the neural biomarker, with temporal features having a stronger influence.

3.3 The identified neural biomarker simultaneously tracked the temporal dynamics of weekly and daily depression symptom variations during LHb DBS treatment

After identifying the neural biomarker of depression symptoms, we evaluated its ability to track the temporal dynamics of slow (weekly) and fast (daily) depression symptom variations during the LHb DBS treatment. For slow weekly variations, we used the identified neural biomarker model to compute weekly neural biomarker values (see Methods section 2.6 for details). We used the weekly neural biomarker values to predict the associated weekly clinical ratings, where we strictly excluded the weekly data that were used to identify the neural biomarker (i.e., the prediction was based on new unseen data not used in training the neural biomarker model). We found that the weekly neural biomarker values significantly predicted the HAMD scores (Figure 4A, EV=0.74, P = 1.1 × 10^!$). Further, considering the temporal dynamics in detail by using the DTW distance analysis (see Methods section 2.6), we found that the weekly neural biomarker significantly tracked the temporal dynamics of weekly HAMD score variations (random shuffle P = 0.00001). Consistently, the weekly neural biomarker values significantly predicted the MADRS scores (Figure 4B, EV=0.34, P = 0.039), and tracked the temporal dynamics in the DTW distance analysis with marginally significant statistics (random shuffle P = 0.1079). Conversely, the weekly neural biomarker values did not predict the HAMA scores (Figure 4C, EV=0.03, P = 5.2 × 10^!%) or tracked the temporal dynamics (random shuffle P = 0.5865).

• Download figure
• Open in new tab

Figure 4.

Accurate tracking the temporal dynamics of weekly depression symptom variations by the identified neural biomarker. (A) Left: correlation between the identified neural biomarker values and the weekly HAMD scores unseen in neural biomarker identification. EV: Explained variance. Middle: neural biomarker tracking of the weekly HAMD score dynamics over time. Right: the DTW distance analysis result for evaluating the significance of tracking in the middle panel. Smaller DTW distance represents better tracking. Note that we normalized both the neural biomarker values and the symptom scales to a range of 0 to 1 using min-max normalization for better visualization. (B) same as (A) but for the weekly MADRS scores. (C) same as (A) but for the weekly HAMA scores.

For fast daily variations, we used the identified neural biomarker model to compute daily neural biomarker values and used the daily neural biomarker values to predict the associated daily self-reports (again, data not used in training the neural biomarker model). We found that the daily neural biomarker values significantly predicted the VAS-D scores (Figure 5A, EV=0.63, P = 1.3 × 10^!#&) and showed significant tracking of VAS-D dynamics (random shuffle P = 0.0001). By contrast, while the daily neural biomarker values predicted the VAS-A scores (Figure 5B, EV=0.51, P = 5.3 × 10^!%“) but the daily neural biomarker did not track VAS-A dynamics (random shuffle P = 1.00).

• Download figure
• Open in new tab

Figure 5.

Accurate tracking the temporal dynamics of daily depression symptom variations by the identified neural biomarker. (A) Left: correlation between the identified neural biomarker values and the daily VAS-D scores unseen in neural biomarker identification. EV: Explained variance. Middle: neural biomarker tracking of the daily VAS-D score dynamics over time. Right: the DTW distance analysis result for evaluating the significance of tracking in the middle panel. Smaller DTW distance represents better tracking. (B) same as (A) but for the daily VAS-A scores.

In summary, the results show that the identified neural biomarker significantly tracked the temporal dynamics of both weekly and daily variations in depression symptoms during the LHb DBS treatment and specifically tracked depression symptoms rather than anxiety symptoms.

3.4 The identified neural biomarker reflected changes of depression symptoms in response to DBS parameter alterations

A useful neural biomarker for DBS also needs to reflect the effect of different DBS parameters. We thus finally evaluated if the identified neural biomarker could reflect changes in depression symptoms in response to DBS parameter alterations. We applied three different DBS frequencies during the treatment: 1 Hz, 20 Hz, and 130 Hz. For the DBS alteration from 1 Hz to stimulation off (Off-2, the DBS device shut down due to unnoticed power off), there was a trend of increasing for the neural biomarker, HAMD, and MADRS while the statistical tests were not significant due the limited sample size (Figure 6A, 1Hz v.s. stimulation off, normalized mean±s.e.m., neural biomarker: 0.912 ± 0.055 v.s. 0.956 ± 0.026, P=0.35; HAMD: 0.562 ± 0.062 v.s. 0.688 ± 0.036, P=0.16; MADRS: 0.656 ± 0.031 v.s. 0.734 ± 0.053, P=0.35), which indicated a rebound trend of depression symptoms due to the disruption of DBS treatment. For the DBS alteration from stimulation off (Off-2) to 20 Hz stimulation (Figure 6B), the neural biomarker, HAMD, and MADRS consistently decreased (Figure 6B, stimulation off v.s. 20 Hz, neural biomarker: 0.905 ± 0.095 v.s. 0.545 ± 0.059, P=0.04; HAMD: 0.875 ± 0.125 v.s. 0.375 ± 0.072, P=0.05; MADRS: 0.875 ± 0.000 v.s. 0.606 ± 0.061, P=0.08). Specifically, the neural biomarker, HAMD and MADRS scores all decreased at the week of DBS frequency alteration, further decreased one week after the alteration and continued to decrease with more obvious changes after week two. For the DBS alteration from 20 Hz to 130 Hz stimulation (Figure 6C), the neural biomarker, HAMD, and MADRS also consistently decreased (Figure 6C, 20 Hz v.s. 130 Hz, neural biomarker: 0.524 ± 0.060 v.s. 0.088 ± 0.029, P=0.001; HAMD: 0.323 ± 0.054 v.s. 0.042 ± 0.026, P=0.004; MADRS: 0.573 ± 0.056 v.s. 0.062 ± 0.016, P=0.001). More specifically, the neural biomarker, HAMD and MADRS scores already showed a trend of decreasing before the DBS frequency alteration, and the alleviated symptoms stayed relatively stable during the alteration week, at week one after the alteration, and the same stable trend continued after week two. The results suggested that the 1 Hz DBS did not induce an obvious change in depression symptoms, while the 20 Hz and 130 Hz DBS had more meaningful effects. The results further demonstrated that the identified neural biomarker indeed reflected the different change patterns in depression symptoms when the DBS frequencies were altered.

• Download figure
• Open in new tab

Figure 6.

The consistent trend between the neural biomarker and the depression symptom changes when DBS frequencies were altered. (A) DBS frequency was altered from 1 Hz to Off-2. Top: comparison of the values of the biomarker, HAMD scores and MADRS scores between the 1Hz stage and Off-2 stage. Two-sided Wilcoxon rank-sum test was used for significance test. Row two: changes in the identified neural biomarker time-locked to the DBS frequency alteration week (vertical dashed line). Row three: changes in HAMD time-locked to the DBS frequency alteration week. Bottom: changes in MADRS time-locked to the DBS frequency alteration week. (B) same as (A) but for DBS frequency alteration from Off-2 to 20 Hz. (C) same as (A) but for DBS frequency alteration from 20 Hz to 130 Hz.

4.1 A data-driven LHb neural biomarker for tracking slow and fast depression symptom variations during DBS treatment

A mechanism-driven neural biomarker for depression is currently lacking mainly because the neural circuitry underlying depression has not been clearly delineated [54]. Therefore, current neural biomarkers of depression symptoms for tracking DBS effects have largely used data-driven machine learning methods to map LFP features to depression symptom ratings [17,31,32]. The usefulness of data-driven neural biomarkers critically depends on the data used to identify the neural biomarker. For example, a recent work [31] focused on a cingulate neural biomarker that was trained with and accordingly predicted longer-term (on the time scale of weeks) clinical ratings. On the other hand, another recent work [17] only trained and tested a neural biomarker with shorter-term (on the time scale of minutes) self-reports. Our work is unique in that while we identified our LHb neural biomarker based on weekly clinical ratings, we demonstrated that the neural biomarker predicted not only weekly clinical ratings (data not used in identification) but also daily self-reports (data again not used in identification). The results suggested that our LHb neural biomarker could track the temporal dynamics of both slow and fast depression symptom variations, which was useful for developing new aDBS strategies that are robust across different time scales. It is worth noting that our LHb neural biomarker specifically tracked the temporal dynamics of weekly and daily depression symptom scores but not the anxiety symptom scores. It suggests that despite the overlapping of depression-related and anxiety-related brain networks [55], LHb neural activity is mainly related to depression, which is supported by prior animal studies [21,22].

Both population-level and personalized neural biomarkers of depression symptoms have been identified for tracking DBS effects. Population-level neural biomarkers are derived from data collected from several patients and have the benefits of being directly applicable to a new patient and robust interpretability of the neural biomarker’s biophysical mechanism across patients [31]. By contrast, personalized neural biomarkers are derived from data collected from an individual patient, which is more powerful in capturing the unique characteristics of depression symptoms in each patient, especially given the large inter-individual variability in depression-related brain networks [56]. With the emerging capability of recording more data within a single patient using mobile devices, personalized neural biomarker models can be more accurate in tracking the temporal dynamics of depression symptom variations. With such trends, a personalized neural biomarker has been identified and used for realizing aDBS targeting VC/VS [17]. Our study identified a personalized neural biomarker that achieved accurate classification and tracking of depression symptom variations during the DBS treatment targeting LHb, confirming the usefulness of personalization. Nevertheless, population-level and personalized neural biomarkers can complement each other. For example, one can leverage a large amount of population data to train an interpretable population-level neural biomarker model, followed by fine-tuning with personalized data to further improve its accuracy, which is an important future research direction.

4.2 Importance of fusing LFP temporal domain and spectral domain features to identify the neural biomarker

Previous studies have exclusively used LFP spectral domain features to identify the neural biomarker for depression [17,31,32]. By contrast, we found that both the LFP temporal domain and spectral domain features contributed to our neural biomarker and that the temporal domain features contributed relatively more than the spectral domain features. The main contributing temporal domain feature was the Hurst exponent. The Hurst exponent measures “long-term memory” in temporal dynamics and is associated with the autocorrelation of a time series [57]. Our results suggest that the LHb LFP signal’s auto-correlations might change during DBS. The main contributing (top 2) spectral domain feature was the PSD of the β band. The LHb β band oscillation was also found to correlate with depression symptoms before DBS treatment in a previous study [27]. These findings highlight the significance of the β band oscillations in LHb as related to depression. Both the Hurst exponent and the β band oscillation features reflect the abnormal synchronization of LHb neural ensembles underlying depression and might be related to the abnormal burst spiking phenomena of LHb neurons found in rodents exhibiting depression-like behaviors [21,22]. The mechanism underlying the contributing LFP temporal domain and spectral domain features requires further investigation by future research.

4.3 DBS frequency deferentially modulates the depression symptom and neural biomarker

Prior studies have shown that DBS frequency can significantly influence the treatment efficacy for TRD, e.g., high-frequency DBS has generally yielded better treatment outcomes than low-frequency DBS [23,58,59]. Consistently, we discovered that a very low DBS frequency of 1 Hz was not effective in alleviating the depression symptoms in our patient, but higher frequencies of 20 Hz and 130 Hz were more effective. Beyond the depression symptom ratings, we additionally found that the neural biomarker was also consistently modulated by the different DBS frequencies. DBS frequency might influence the release of neurotransmitters in depression-targeted pathways [60], thus modulating the neural biomarker and the depression symptoms. However, similar to other DBS targets, the optimal DBS frequency at LHb is still also an open question that requires further research.