Evaluating tDCS Intervention Effectiveness via Functional Connectivity Network on Resting-State EEG Data in Major Depressive Disorder

Major depressive disorder (MDD) is a mental illness marked by an enduring state of melancholy, guilt, worthlessness, and hopelessness, which may increase the risk of suicidal behavior[1]. MDD is the most common type of unipolar affective disorder, characterized by abnormal brain activity[2]. According to data released by the World Health Organisation (WHO), over $350$ million individuals worldwide suffer from MDD and the number of new cases has increased by almost $18\%$ over the past ten years[3].

tDCS involves the application of low-intensity electrical currents to specific brain regions via scalp electrodes, resulting in the modulation of neuronal excitability and synaptic activity[4]. Unlike other neuromodulation techniques, such as selective serotonin reuptake inhibitors[5], tDCS is portable, well-tolerated, and relatively inexpensive, making it an attractive option for both research and clinical applications.

tDCS is considered for treating MDD because of evidence suggesting problems in neural circuitry (abnormalities in the connections between different brain regions) and neurotransmitter dysregulation (imbalance in the levels or activity of neurotransmitters) contributing to depression. Neuroimaging studies have consistently identified alterations in brain structure and function, particularly within corticolimbic circuits involved in emotional regulation and mood processing, in individuals with MDD[6]. Dysfunction within these circuits, including the dorsolateral prefrontal cortex (DLPFC)-limbic system pathway, has been implicated in developing and maintaining depressive symptoms[7, 8].

The antidepressant effects of tDCS in MDD have been investigated in numerous clinical trials and meta-analyses, with accumulating evidence supporting its efficacy and tolerability as a standalone treatment or adjunctive therapy[9]. Additionally, studies have explored the potential of tDCS to induce sustained antidepressant effects through repeated or maintenance treatments, as well as its combination with other therapeutic modalities to enhance treatment outcomes[10].

As research into transcranial Direct Current Stimulation (tDCS) for Major Depressive Disorder (MDD) continues to expand, several critical challenges remain. These include optimizing stimulation parameters, elucidating underlying neurobiological mechanisms, and identifying predictors of treatment response. Addressing these challenges is essential for advancing our understanding of tDCS as a therapeutic intervention for MDD and enhancing its clinical efficacy.

This paper aims to contribute to our understanding of the alterations in neural circuitry induced by tDCS and its potential role in the management of depression.

Functional connectivity analysis and network neuroscience are valuable tools in neuroimaging for understanding brain function. tDCS has gained interest for its potential in treating MDD. Though its exact mechanisms are complex, research has highlighted key neurobiological processes involved in tDCS-induced antidepressant effects. Functional connectivity examines statistical associations between brain regions, revealing networks underlying cognition and emotion. Network neuroscience characterizes brain networks mathematically, assessing properties such as average path length and clustering. In depression-related research, these approaches have shown altered connectivity and network organization, where tDCS has been found to modulate these networks.

Nitsche et al. (2008) proposed that anodal stimulation, by promoting neuronal depolarization, enhances synaptic strength and induces long-term potentiation (LTP) in targeted brain regions, and cathodal stimulation, by decreasing neuronal depolarization, induces long-term depression (LTD) and reduces cortical excitability [4]. This modulation of neural activity within mood-relevant circuits is thought to restore dysfunctional neural circuits associated with MDD and alleviate depressive symptoms.

Animal studies have been instrumental in elucidating the neuroplasticity-related changes induced by tDCS. For example, Monai et al. (2016) demonstrated that tDCS leads to enduring alterations in synaptic plasticity and neurotransmitter levels, particularly within the glutamatergic and gamma-aminobutyric acid (GABA) systems [11]. These findings suggest that tDCS may exert its therapeutic effects by modulating neuronal excitability and synaptic strength in specific brain regions implicated in mood regulation.

Research on functional connectivity in MDD has revealed several consistent findings. Mulders et al. (2015) identified increased connectivity within the anterior default mode network, between the salience network and the anterior default mode network, and decreased connectivity between the posterior default mode network and the central executive network [12]. Fattahi et al. (2021) further explored these changes, finding significant differences in connectivity between various network pairs in MDD patients with suicidal thoughts [13]. Zhi et al. (2018) identified the disrupted topological organization of dynamic functional network connectivity in MDD, with patients spending more time in a weakly connected state associated with self-focused thinking [14].

The brain is an extremely complex multivariate dynamical network. Network neuroscience employs mathematical tools to model these intricate network systems using the concept of nodes (also called vertices) and edges[15]. In the context of EEG connectivity, nodes represent electrodes, while edges represent the connectivity between electrode pairs. Thus, brain networks and brain functional connectivity are depicted in either matrix or graph form.

Past studies have highlighted the application of graph theory in understanding brain networks. Chung (2021) and Vecchio et al. (2017) both explored the scale-free and small-world properties of these networks, with the latter also discussing the potential for graph theory to aid in understanding brain disconnection and monitoring treatment impact [16, 17]. Additional studies have underscored the significance of graph theory in unravelling the architecture, development, and evolution of brain networks[18, 19].

Studies of resting-state EEG data in MDD using functional connectivity and graph theoretical approaches suggest that depressed patients exhibit altered brain functional connectivity patterns compared to healthy individuals. These studies highlight the potential of multilayer brain functional network frameworks for analyzing abnormal brain interaction patterns in depression[20]. Furthermore, the effect of tDCS in modulating resting-state functional connectivity, particularly in regions distal to the stimulation site, has been emphasized[21]. Research utilizing graph theory to analyze brain networks in MDD has revealed several key findings. Hasanzadeh et al. (2020) and Hasanzadeh et al. (2017) both found that MDD patients exhibit a more randomized structure in their brain networks, with higher node degree and strength [22, 23]. Sun et al. (2019) further confirmed these findings by identifying deficiencies in right hemisphere function and a randomized network structure in MDD [24]. Ye et al. (2015) added to this finding by demonstrating that MDD patients have higher local efficiency and modularity in their brain networks, as well as altered nodal centralities in specific brain regions, implying that emotional and cognitive functions are significantly affected by MDD [25].

These insights support the hypothesis that statistically significant differences exist in binarized network structures in MDD. These differences are quantified using network neuroscience measures such as small-worldness, hubness, and asymmetry. This approach provides a promising avenue for further exploration and understanding of MDD pathology.

Individuals diagnosed with MDD were administered tDCS to manage symptoms. Resting-state EEG data were recorded before and after the intervention.

1. 1.

   Pre-processing: Pre-processing of the EEG data involved several steps. First, all EEG data in both groups were average re-referenced to minimize the impact of common noise sources[27]. Following this preprocessing, band-pass filtering was applied with a frequency range of $0.01-45$ Hz to isolate the EEG signals within the desired frequency band and to remove unwanted noise and artefacts. Then, the EEG data were cropped and epoched into segments with a duration of $4$ seconds each, with an overlapping window of $2$ seconds between consecutive epochs[28].

2. 2.

   Power against frequency plot: We plotted power against frequency for both pre-intervention and post-intervention groups, within the frequency range between $0.01-45Hz$.

3. 3.

   Division of frequency bands and parcellation scheme The division of frequency bands was as follows: theta ranges from $4-8Hz$, alpha ranges from $8-13Hz$ and beta ranges from $13-30Hz$. We included Fp1, AF3, F3 and F7 in the left frontal region; Fp2, AF4, F4 and F8 in the right Frontal region; T7 in the left temporal region; T8 in the right temporal region; FC5, FC1, C3, CP1 and CP5 in the left central region; FC2, FC6, C4, CP2 and CP6 in the right central region; P7, P3, PO3 and O1 in the left parietal-occipital region; P4, P8, PO4 and O2 in the right parietal-occipital region.

4. 4.

   Phase-based connectivity: In phase-based connectivity, we used frequency domain features. Any signal may be analyzed in the frequency domain using the Fourier transform, which provides a frequency-domain representation of EEG data[29]. However, the Fourier transform has some limitations[30]. Changes in frequency content in the signal over time are difficult to interpret using the Fourier transform, and it works well primarily for stationary signals. Since brain signals are dynamic and non-stationary, time-frequency decomposition is required to retain the advantages of both the time and frequency domains, though with some sacrifice in temporal and frequency precision[31]. We used the Morlet wavelet (Gabor wavelet) for time-frequency analysis.

5. 5.

   Phase lag index (PLI): When two signals are perfectly in phase, they are synchronized. PLI works by focusing on the consistency of the phase differences between signals, rather than the absolute phase difference. With sensitivity towards the directionality of the phase difference between two signals, it measures the consistency of the phase lead/lag relationship between the signals, rather than just the magnitude of the difference. The fundamental idea is to disregard phase locking centred around 0 phase difference to exclude volume conduction effects (at the risk of ignoring true instantaneous interactions). This also applies to phase locking at $\pi,2\pi$ and so on, i.e., repeating at every $\pi$, also given as 0 mod $\pi$. The PLI ranges between $0$ and $1$, with $0$ indicating no coupling due to volume conduction and $1$ indicating true, lagged interaction. We calculated it using multitaper which employs multiple windows to estimate the spectrum and thereby reduces the variance of the measurement. Mathematically it is expressed as follows[32]:

   $$PLI=\left|n^{-1}\sum_{t=1}^{n}sgn\left(Im\left(e^{i\left(\phi^{j}-\phi^{k}\right)_{t}}\right)\right)\right|$$

   (1)

   The vectors are not averaged with the PLI; instead, the sign of the imaginary part of the cross-spectral density is averaged. The PLI will be large if all phase angle variations are on one side of the imaginary axis. Conversely, the PLI will be zero if half of the phase angle discrepancies are positive and half are negative relative to the imaginary axis.

6. 6.

   Small world network index (SWI): The “small-world” phrase was derived from Stanley Milgram’s 1967 study which demonstrated that people are connected through an average of six degrees of separation in social networks, revealing unexpectedly short paths between individuals[33]. Expanding upon this concept, Watts and Strogatz introduced the small-world network [34]. Before small-world networks, researchers typically modelled real-world networks using either random or regular networks. The concept of small-world networks revolutionized this understanding by showing that many real-world networks exhibit a small-world property, characterized by a high clustering coefficient like regular networks and a short average path length like random networks. Bassett and Bullmore implemented this model for brain anatomical and functional networks[35] demonstrating its capacity to support both segregated (specialized) and distributed (integrated) information processing. The SWI of a network is calculated as follows:

   $$\text{SWI}=\frac{C}{L_{\text{avg}}}$$

   (2)

   where $C$ is the clustering coefficient of a network and $L_{\text{avg}}$ is the average shortest path length of the network.

   Permutation comparison of random network SWI, pre-intervention SWI and post-intervention SWI: We compared the SWI metric of the random network with that of the pre-intervention and post-intervention groups using the following algorithm:

   1. i.

      The algorithm iterates through a range of 200 thresholds ($\Theta$) from $0-1.5$ for each frequency band ($\beta_{f}$), where $\beta_{f}\in\{\theta,\alpha,\beta\}$. The upper bound of $1.5$ is based on experimental observations indicating that values exceeding this range result in completely disconnected random networks.

   2. ii.

      At each threshold ($\Theta$), we compute the threshold value ($\Theta_{\text{thresh}}$) for both pre-intervention and post-intervention groups ($\Gamma$), where $\Gamma\in\{\text{pre},\text{post}\}$ and for each frequency band $\beta_{f}\in\{\theta,\alpha,\beta\}$. We define $\Phi$ as the functional connectivity matrix, which depends on both the group and frequency band i.e., $\Phi(\Gamma,\beta_{f})$. Using $\Phi$ and $\Theta$, we calculate $\Theta_{\text{thresh}}$:

      $$\Theta_{\text{thresh}}=\text{median}(\Phi)+\Theta\cdot\sigma(\Phi)$$

      (3)

      where, $\Theta_{\text{thresh}}$ ranges for each connectivity matrix ($\Phi$) from $\text{median}(\Phi)$ to $\text{median}(\Phi)+1.5\cdot\sigma(\Phi)$.

      The different $\Theta_{\text{thresh}}$ values for both treatment groups in $\Gamma$ and in each $\beta_{f}$ facilitate the calculation of the binary matrix ($\mathcal{B}$). This matrix is defined as a function of $\Gamma$, $\Theta_{\text{thresh}}$, and $\beta_{f}$:

      $$\mathcal{B}(\Gamma,\beta_{f},\Theta_{\text{thresh}})=\Phi(\Gamma,\beta_{f})>\Theta_{\text{thresh}}$$

      (4)

      Furthermore, we compute the number of connections ($N_{\text{conn}}$) and the SWI in $\mathcal{B}(\Gamma,\beta_{f},\Theta_{\text{thresh}})$.

   3. iii.

      Firstly, for each $\beta_{f}$ and $\Theta_{\text{thresh}}$, we generate $100$ Erdős–Rényi random models ($\mathcal{R}$). The number of connections in these randomly created networks is set to the average $N_{\text{conn}}$ of $\mathcal{B}$ for both the groups in $\Gamma\text{ i.e., }\{\text{pre},\text{post}\}$. We then compute the SWI ($\text{SWI}_{\text{random}}(\mathcal{R}_{\text{i}})$), clustering coefficient ($\text{C}_{\text{random}}(\mathcal{R}_{\text{i}})$) and average path length ($\text{L}_{\text{avg, random}}(\mathcal{R}_{\text{i}})$) for all $100$ random networks where $i\in\{0,1,2,\ldots,99\}$.

      Secondly, we randomly select two distinct networks $\mathcal{R}_{j}$ and $\mathcal{R}_{k}$ from the $100$ generated networks. The selected networks must exhibit some clustering ($C\neq 0$) and not have all nodes connected ($L_{\text{avg}}\neq 0$). We then calculate $\text{SWI}_{\text{permutations}}(x)$ as follows:

      $$\text{SWI}_{\text{permutations}}(x)=\frac{\text{SWI}_{\text{random}}(\mathcal{R}j)}{\text{SWI}_{\text{random}}(\mathcal{R}_{k})}$$

      (5)

      This process is repeated $1000$ times ($x\in[0,999]$), generating a small-world network distribution of random networks centred around $1.0$. These $1000$ permutations represent random networks with equivalent topological properties but shuffled connections, creating a null distribution of small-world network values, $\text{SWI}_{\text{permutations}}(x)$. We compare the observed small-world network value of the original networks for both groups against this null distribution. The observed real small-world network value is computed as:

      $$\text{SWI}_{\text{real}}(\Gamma,\beta_{f})=\frac{\text{SWI}(\Gamma,\beta_{f})}{(\frac{\text{$\mu$}(C_{\text{random}}(\mathcal{R}))}{\text{$\mu$}(L_{\text{avg, random}}(\mathcal{R}))})}$$

      (6)

      where $\text{SWI}(\Gamma,\beta_{f})$ is calculated on $\Phi(\Gamma,\beta_{f})$. This comparison enables us to assess the significance of the network’s small-world properties relative to random networks. Our threshold algorithm adapts a statistical approach for selecting a threshold based on normalized small-worldness[36]. We conducted multiple experiments with this modified algorithm to optimize threshold selection by accounting for small-worldness in both pre-intervention and post-intervention groups.

      Finally, we compute the $p-value$ for both groups by comparing $\text{SWI}_{\text{real}}$ against the $\text{SWI}_{\text{permutations}}$ distribution.

   4. iv.

      Optimal threshold selection: Z-scores are calculated for each $\text{SWI}_{\text{real}}$ at each threshold for both groups to quantify the deviation of observed small-world network metrics from the null distribution. Ideally, we should select a threshold for which the small-world network metrics (SWI) are high for both groups and the difference between them is maximized. However, we have chosen a threshold for which the normalized small-world network metric exhibits the maximum difference between the groups (Figure 5).

7. 7.

   Hubness: We also implemented the concept of ’hub’ or ’hubness’ from graph theory in this section[37]. A hub is a significant node in a network with numerous incoming and outgoing branches or connections to other nodes. In our case, a hub within an undirected network is characterized as a node with extensive connections to other nodes, with many nodes being indirectly linked through the hub. Consequently, if a hub node is removed, the network’s functionality is significantly disturbed. Hubness implies a connectivity measurement indicating a node’s potential to become a hub in a network. It is calculated for the difference matrix between the pre-intervention binarized group $\mathcal{B}(\Gamma_{\text{pre}},\beta_{f},\Theta_{\text{thresh}})$ and the post-intervention binarized group $\mathcal{B}(\Gamma_{\text{post}},\beta_{f},\Theta_{\text{thresh}})$. This involves counting the number of edges for each node, establishing a threshold to segregate channels with higher hubness values and identifying these channels as cortical regions of interest.

8. 8.

   Asymmetry and the number of links: We created a difference matrix by subtracting the pre-binarized matrix from the post-binarized matrix. This approach allows us to focus on the activity that increases post-intervention compared to pre-intervention in specific regions[38]. The asymmetry ratio between the left and the right brain hemispheres is calculated by comparing the number of edges connecting to any left hemispheric nodes with those connecting to any right hemispheric nodes in the binary matrix obtained from the optimal threshold for each frequency band. We analyze the number of links after binarization to characterize the connectivity differences between both groups.

In this study, we evaluated the efficiency of tDCS for MDD by analyzing resting-state EEG data and assessing functional connectivity using network neuroscience across different frequency bands. A novel aspect of our research was modifying the binarizing thresholding algorithms to provide a detailed understanding of network topology changes following tDCS treatment. Despite these advancements, limitations such as a small sample size, data variability, and reliance on thresholding methods must be addressed. Future research with larger sample sizes and longitudinal designs is necessary to confirm these findings and explore the long-term effects of tDCS on brain connectivity in MDD.

Our findings contribute to the literature on the neural mechanisms of MDD and the potential for tDCS as a therapeutic intervention. We observed significant changes in brain network randomization in the beta band post-tDCS. Additionally, differences in PSD in the alpha band were noted between pre-intervention and post-intervention groups. Moreover, an increased number of connectivity links in the theta and alpha bands were found pre-intervention which were more scattered but became more localised in the DLFPC region in all the bands. Furthermore, left hemispheric asymmetry increased following tDCS. These insights could inform the development of personalized treatment strategies for MDD and other psychiatric conditions. Overall, while tDCS shows promise as a non-invasive treatment for MDD, further research is needed to validate and refine its clinical application.

Funding: This research received no external funding.
Conflict of interest: The authors declare no conflict of interest.

The authors thank all the twelve participants who facilitated the analysis of the tDCS intervention and the technicians for the extensive data collection procedure. We acknowledge using Elicit: Notebook¹¹1https://elicit.com/notebook/ to assist in generating relevant research papers.