Heart Rate Variability Biofeedback Alleviates Subthreshold Depression by Reconstructing the Brain–Heart Axis via Habenular Network Functional Connectivity Modulation

Figures

• 

  FIGURE 1. Workflow of the study. (a) Study procedure: Including participant recruitment, HRV‐BF training, and rs‐BOLD fMRI scanning. (b) Study assessments: Comprising scale assessment, cardiovascular measurement, and bilateral habenula‐based functional connectivity analysis. The red markers indicate the anatomical location of the bilateral habenula in axial, coronal, and sagittal planes. FC, functional connectivity; HAMD‐17, 17‐item Hamilton Depression Rating Scale; HR, heart rate; HRV, heart rate variability; PHQ‐9, Patient Health Questionnaire‐9; PSQI, pittsburgh sleep quality index.

• 

  FIGURE 2. Group differences in HR, BP, and HRV parameters between individuals with StD and HC. (a–c) SBP, DBP, and HR; (d–e) lnSDNN and lnRMSSD; (f, g) HF power and LF power; (h) HRV amplitude. Variables including SBP, HR, HF power, and LF power followed a normal distribution and were analyzed using independent samples t‐tests. Non‐normally distributed variables were analyzed using the Mann–Whitney U test. Data are presented as violin plots with embedded boxplots (median and interquartile range) to illustrate the distribution of data and central tendency. * p < 0.05, **p < 0.01, ***p < 0.001. StD group (n = 34), HC group (n = 32). DBP, diastolic blood pressure; HF, high frequency; HR, heart rate; HRV amp., heart rate variability amplitude; LF, low frequency; lnRMSSD, natural logarithm of the root mean square of successive differences; lnSDNN, natural logarithm of the standard deviation of NN intervals; SBP, systolic blood pressure.

• 

  FIGURE 3. Habenula‐ROIs functional connectivity analysis (a, e, i, m, q) Axial, coronal, and sagittal views showing the anatomical locations of five ROIs: NAc, VeP, AMY, GPi, and SN. ROIs are marked in red. (b, f, j, n, r) Partial correlation analyses (Spearman) between Hb‐ROI zFC values and HAMD‐17 scores, controlling for sex and age. (c, g, k, o, s) Partial correlation analyses (Spearman) between Hb‐ROI zFC values and PHQ‐9 scores, controlling for sex and age. (d, h, l, p, t) ROC curves evaluating the discriminative ability of Hb‐ROI zFC values in distinguishing StD from HC. Hb‐NAc (AUC = 0.876, 95% CI, 0.699–1.000, p < 0.001), Hb‐VeP (AUC = 0.876, 95% CI, 0.725–0.980, p < 0.001), Hb‐AMY (AUC = 0.791, 95% CI, 0.556–0.974, p = 0.080), Hb‐GPi (AUC = 0.778, 95% CI, 0.582–0.941, p = 0.002) and Hb‐SN (AUC = 0.706, 95% CI, 0.484–0.892, p = 0.056). Total participants: 26 (StD group = 9, HC group = 17). AMY, amygdala; GPi, globus pallidus internus; NAc, nucleus accumbens; SN, substantia nigra; VeP, ventral pallidum.

• 

  FIGURE 4. Whole‐brain voxel‐wise functional connectivity analysis based on bilateral habenula (Hb). Left panel: Brain regions showing significantly increased functional connectivity with bilateral Hb in healthy controls (HC) compared to individuals with subthreshold depression (StD) (HC > StD). Right panel: Brain regions showing significantly increased functional connectivity with bilateral Hb after HRV‐BF intervention compared to pre‐intervention (HRV‐BF post > HRV‐BF pre). Statistical threshold was set at voxel‐wise p < 0.001 (uncorrected), with a cluster‐level p < 0.05 (FWE corrected). Color intensity (from red to yellow) indicates the strength of functional connectivity with bilateral Hb. StD group (n = 8, female), HC group (n = 11, female).

• 

  FIGURE 5. Effects of heart rate variability biofeedback on the autonomic nervous system. (a–c) Changes in cardiovascular physiological parameters: SBP (a), DBP (b), HR (c). (d, e) Changes in HRV time‐domain parameters: LnSDNN (d) and lnRMSSD (e). (f, g) Changes in HRV frequency‐domain parameters: HF power (f) and LF power (g). (h) Changes in HRV amplitude. (i) Changes in respiration rate. Paired‐samples t‐test was used to assess pre‐post differences, n = 8. * p < 0.05, **p < 0.01, ***p < 0.001.

• 

  FIGURE 6. Effects of HRVBF on depressive symptoms, sleep quality, and bilateral habenular‐ROIs functional connectivity. (a–c) changes in scale assessments, including HAMD‐17 (a), PHQ‐9 (b), and PSQI (c). Paired raincloud plots summarize the intervention effect after HRVBF. Each plot integrates a half‐violin plot (showing kernel density estimation), a boxplot (indicating the median and interquartile range), individual paired data points with connecting lines (depicting within‐subject changes), and central mean ± SD bars (highlighting overall group‐level trends). n = 8. (d–h) show comparisons of zFC values between the bilateral Hb and ROIs before and after HRV‐BF training, including Hb‐NAc, Hb‐VeP, Hb‐AMY, Hb‐SN, and Hb‐GPi. Data are displayed as boxplots representing the median and interquartile range. n = 8. Normality of the pre‐post difference (delta values) was tested for each parameter (a–h). Variable dHAMD‐17 did not meet the normality assumption and was analyzed using the Wilcoxon signed‐rank test. The remaining variables were analyzed using paired‐samples t‐tests. * p < 0.05, ** p < 0.01, *** p < 0.001. (i) Functional connectivity changes (delta zFC value) of Hb with key regions…