Large-scale dysfunctional white matter and grey matter networks in patients with social anxiety disorder.

Dysfunction of large-scale brain networks has been implicated in social anxiety disorder (SAD); most work has focused on grey matter (GM) functional connectivity (FC) abnormalities, whereas white matter (WM) FC alterations remain unclear. Here, using a K-means clustering algorithm, we obtained 8 GM and 10 WM functional networks from a cohort dataset (48 SAD patients and 48 healthy controls). By calculating and comparing FC matrices between SAD group and healthy controls, we demonstrated disrupted connections between the limbic and dorsal prefrontal, lateral temporal, and sensorimotor networks, and between the visual and sensorimotor networks. Furthermore, there were negative correlations between HAMD scores and limbic-dorsal prefrontal and limbic-sensorimotor networks, and between illness duration and sensorimotor-visual networks. These findings reflect the critical role of limbic network, with its extensive connections to other networks, and the neurobiology of disordered cognition processing and emotional regulation in SAD.

Introduction

Social anxiety disorder (SAD) is the most common anxiety disorder, and causes a major social and economic burden (den Boer, 2000; Stein and Stein, 2008); 12-month and lifetime prevalence is 7 and 12%, respectively (Kessler et al., 2005; Ruscio et al., 2008). SAD patients are characterized by persistent and extreme fear, anxiety, and avoidance in social settings (Zhang et al., 2020). Patients with SAD can be debilitated, with significant functional impairments especially in relationships and social experiences (Mizzi et al., 2021). However, the neural substrates of SAD are not well-defined. With the rapid development of psychoradiology which is an emerging subspecialty of radiology, an abundance of magnetic resonance imaging (MRI) studies have provided key information for brain funtional changes in SAD (Gong, 2020, Li et al., 2021, Lui et al., 2016, Zhang et al., 2022).

Functional magnetic resonance imaging (fMRI) studies have helped us understand SAD in terms of impaired communication between large-scale functional networks. Especially implicated is the limbic network: the limbic system is involved in the ‘emotional’ processing of input from sensory systems and relaying it to the telencephalon, notably the prefrontal cortex and cingulate gyrus (Catani et al., 2013; Morgane et al., 2005). Deficits in these pathways might be the driving force behind pathological anxiety (Cannistraro and Rauch, 2003). Disturbances in functional connectivity (FC) between the limbic and prefrontal networks are a core feature of SAD (Mizzi et al., 2021), and abnormal FC has also been found between the limbic and temporal networks, and between the sensorimotor and visual networks (Liao et al., 2010a; Mizzi et al., 2021).

As well as known abnormalities in grey matter (GM), a contributing factor in the pathophysiology of SAD is disturbance in the white matter (WM) fiber tracts which interconnect the GM, notably the uncinate fasciculus (connecting prefrontal cortex and limbic regions) and superior longitudinal fasciculus (passing from the frontal cortex to the parietotemporal cortex) (Jenkins et al., 2016, Sun et al., 2015). There is growing evidence for neural activation and functional organization in WM (Ding et al., 2018; Gawryluk et al., 2014; Gore et al., 2019; Li et al., 2019; Peer et al., 2017). Specific WM activation has been found during visual, tactile, and motor task fMRI (Fabri et al., 2011; Gawryluk et al., 2011) and in resting-state fMRI (Ding et al., 2016; Marussich et al., 2017; Mezer et al., 2009), and the synchronization of neural activity in WM tracts and cortical GM regions suggests a common functional encoding (Ding et al., 2018). Furthermore, large-scale interacting WM functional networks correspond to WM tracts identified by diffusion tensor imaging (Peer et al., 2017). Interactions between WM and GM activation are therefore a powerful new tool to investigate physiology and pathophysiology, which has been used in schizophrenia (Jiang et al., 2019a), major depressive disorder (Zhao et al., 2020), and benign epilepsy (Jiang et al., 2019b).

This study therefore aimed to investigate FC changes in GM and WM functional networks and their interaction in SAD, and to define the relationship between these and demographic and clinical variables. According to previous reports of network abnormalities, we hypothesized that SAD patients would demonstrate widespread disruption within or between the GM and WM networks, particularly in the limbic network and prefrontal network, and that these would be associated with the clinical characteristics of SAD.

Results

Participants

Three patients were excluded because of excessive head movement (translation > 2.5 mm or rotation >2.5⁰). Thus, we finally included 48 SAD patients (male = 31; female = 17; mean age = 27.7) and 48 HC subjects (male = 31; female = 17; mean age = 24.3). The demographic and clinical characteristics are summarized in Table 1. Among the 48 patients, 30 were medication-naïve, 12 had taken selective serotonin reuptake inhibitors (SSRIs) with a median time since medication of 40 days (range 10 days - 3.5 years), and for the remaining 6 patients’ medical information was unavailable. There were no significant differences between two groups in age and gender (p > 0.05).

Table 1

Sample characteristics

Characteristic	Mean	SAD (N = 48)	Range	Mean	HC (N = 48)	Range	P
		SD			SD
Age (years)	27.7	7.8	18–50	24.3	4.1	18–50	0.959†
Gender (Male/Female)	31/17	–	–	31/17	–	–	1.000‡
Disease duration (years)	8	6.3	1–30	–	–	–	–
LSAST	67.8	26.2	23–139	–	–	–	–
LSASF	33.5	12.6	14–68	–	–	–	–
LSASA	34.3	14.8	4–71	–	–	–	–
HAMD∗	13.0	7.5	0–32	–	–	–	–
HAMA∗	14.9	7.7	1–36	–	–	–	–

†p-value by two-sample t-test. ‡p-value by two-tailed chi-square test. ∗Data for 7 subjects were missing. Abbreviations: HAMA, Hamilton Anxiety Rating Scale; HAMD, Hamilton Depression Rating Scale; HC, healthy controls; LSAST, LSASF and LSASA, total score and fear and avoidance factor of Liebowitz Social Anxiety Scale (LSAS); N, number; SAD, social anxiety disorder; SD, standard deviation

Identification of 8 grey matter and 10 white matter functional networks

The clustering results identified 8 GM and 10 WM stable functional networks (the largest number of networks with Dice coefficient > 0.9) (Figures S1 and S2). The 8 GM networks were: GM1, default mode network (DMN) including the medial frontal, posterior cingulate, precuneus, and lateral parietal cortices; GM2, mesencephalic cerebellum network including the mesencephalic cerebellar cortex; GM3, sensorimotor network (SMN) including the precentral, postcentral, and paracentral cortices; GM4, limbic network including the amygdala, hippocampus, temporal pole, and orbital frontal cortices (OFC); GM5, dorsal attention network (DAN), including the superior parietal, posterior inferior temporal, and middle frontal cortices; GM6, frontal-parietal network (FPN) including the inferior parietal, middle inferior temporal, and superior and middle prefrontal cortices; GM7, visual network (VIS) including the occipital cortex; and GM8, dorsal prefrontal network including the dorsal anterior cingulate/medial superior frontal, frontoinsular, and anterior parietal cortices (Figure S1 and Table 2).

Table 2

GM and WM functional networks

Number	Network name	Main region
GM1	Default mode network	Medial frontal, posterior cingulate, precuneus, and lateral parietal cortices
GM2	Mesencephalic cerebellum network	Mesencephalic cerebellar cortex
GM3	Sensorimotor network	Precentral, postcentral, and paracentral cortices
GM4	Limbic network	Amygdala, hippocampus, temporal pole, and the orbital part of frontal cortices
GM5	Dorsal attention network	Superior parietal, posterior inferior temporal, and middle frontal cortices
GM6	Frontal-parietal network	Inferior parietal, middle and inferior temporal, and superior and middle prefrontal cortices
GM7	Visual network	Occipital cortex
GM8	Dorsal prefrontal network	Dorsal anterior cingulate/medial superior frontal, frontoinsular, and temporoparietal cortices
WM1	Corona radiata network	Internal capsule, centrum semiovale
WM2	Medial and orbital prefrontal network	Medial and orbital prefrontal white matter
WM3	Deep bundle network	Internal and external capsule, cerebral peduncle
WM4	Anterior cingulum network	Anterior and subcallosal cingulum white matter, rostrum and genu of corpus callosum
WM5	Sensorimotor network	Precentral, postcentral, and paracentral white matter
WM6	Lateral temporal network	Temporal white matter (inferior longitudinal fasciculus)
WM7	Middle cingulum network	Retrosplenial cingulum white matter, splenium of corpus callosum
WM8	Visual network	Occipital white matter
WM9	Cerebellar network	Cerebellar and brainstem white matter
WM10	Posterior cingulum network	Splenium of corpus callosum, parietal white matter

Abbreviations: GM, grey matter; WM, white matter. See also Figures S1 and S2.

The 10 WM networks were: WM1, corona radiata network including the internal capsule, and centrum semiovale; WM2, medial and orbital prefrontal network including the medial, and orbital prefrontal WM; WM3, deep bundle network including the internal and external capsule, and cerebral peduncle; WM4, anterior cingulum network (ACN) including the anterior and subcallosal cingulum WM, rostrum and genu of corpus callosum; WM5, sensorimotor network (SMN) including the precentral, postcentral, and paracentral WM; WM6, lateral temporal network including the temporal WM (inferior longitudinal fasciculus); WM7, middle cingulum network including the retrosplenial cingulum WM, and splenium of corpus callosum; WM8, visual network (VIS) including the occipital WM; WM9, cerebellar network including the cerebellar and brainstem WM; and WM10, posterior cingulum network including the splenium of corpus callosum, and parietal WM (Figure S2 and Table 2).

Functional connectivity differences between the groups

Between-group FC differences in the three kinds of networks are shown in Figure 1 and Table 3. GM network FC changes: in SAD patients compared with HC, GM3 and GM8 showed significantly decreased FC with GM4 (p < 0.05, Bonferroni correction); no GM networks showed increased FC in SAD. WM network FC changes: in SAD patients compared with HC, WM5 showed significantly decreased FC with WM8 (p < 0.05, Bonferroni correction); no WM networks showed increased FC in SAD. GM-WM network FC changes: in SAD patients compared with HC, there was significantly decreased FC between GM4 and WM6 and between GM7 and WM5 (p < 0.05, Bonferroni correction); no GM-WM networks showed increased FC in SAD. Subgroup analysis identified no significant FC changes between medicated patients and medication-naïve patients after Bonferroni correction (p > 0.05, Figure S3).

Figure 1

Decreases in functional connectivity (FC) between grey matter (GM) networks (yellow links), between white matter (WM) networks (blue links), and between grey and white matter networks (green links) in patients with SAD compared to healthy controls

Table 3

Decreases in functional connectivity between networks in SAD Group compared with HC

Functional connectivity between brain networks	FC values∗	HC	t	P
	SAD
Functional Connectivity between GM networks
GM4-GM3	0.61 ± 0.15	0.71 ± 0.20	−3.44	< 0.001
GM4-GM8	0.51 ± 0.18	0.65 ± 0.19	−4.03	< 0.001
Functional connectivity between WM networks
WM5-WM8	0.60 ± 0.20	0.73 ± 0.13	−3.66	< 0.001
Functional connectivity between GM and WM networks
GM4-WM6	0.72 ± 0.14	0.83 ± 0.09	−4.37	< 0.001
GM7-WM5	0.61 ± 0.20	0.73 ± 0.12	−3.59	< 0.001

∗ Data are presented as means ± SD. Dashes (−) indicate FC between the two identified networks (the networks are defined in Table 2). The table lists Student t-values and associated p-values; all p-values survived Bonferroni correction. Abbreviations: FC, functional connectivity; HC, healthy controls; GM, grey matter; SAD, social anxiety disorder; WM, white matter.

Correlation between altered FC and clinical variables

In the SAD group, HAMD score was negatively correlated with GM3-GM4 FC (r = −0.324, p = 0.039) (Figure 2A) and GM4-GM8 FC (r = −0.443, p = 0.004) (Figure 2B); and illness duration was negatively correlated with GM7-WM5 FC (r = −0.314, p = 0.036) (Figure 2C) and WM5-WM8 FC (r = −0.362, p = 0.014) (Figure 2D). No significant associations were found between FC and age or LSAS scores (p > 0.05, Table S1).

Figure 2

Relationships between the functional connectivity (FC) values and illness duration and clinical symptoms

(A) FC between sensorimotor network (GM3) and the limbic network (GM4) was negatively associated with HAMD scores.

(B) FC between dorsal prefrontal network (GM8) and the limbic network (GM4) was negatively associated with HAMD scores.

(C) FC between visual network (GM7) and the sensorimotor network (WM5) was negatively associated with disease duration.

(D) FC between visual network (WM8) and the sensorimotor network (WM5) was negatively associated with disease duration. The blue line represents the linear trend line, and shaded areas 95% confidence intervals. Abbreviations: GM, grey matter; WM, white matter. See also Table S1.

Discussion

To the best of our knowledge, this study is the first to reveal in patients with SAD the intrinsic functional networks of both WM and GM and their interactions. We identified 8 GM and 10 WM intrinsic functional networks, and found several FC abnormalities in SAD patients compared with healthy controls: the GM of the limbic network (GM4) showed decreased FC with the sensorimotor network (GM3), the dorsal prefrontal network (GM8), and the lateral temporal network (WM6); the WM of the sensorimotor network (WM5) showed decreased FC with the GM and WM of the visual network (GM7 and WM8); FC between the GM3-GM4 networks and between the GM4-GM8 networks were negatively correlated with HAMD scores; and FC between the GM7-WM5 networks and between the WM5-WM8 networks were negatively correlated with illness duration. We discuss the possible physiological significance of these abnormalities briefly in turn.

Decreased connectivity between the limbic (GM4) and dorsal prefrontal networks (GM8)

The limbic network plays a vital role in processing emotional stimuli in social fear and anxiety (Kim and Yoon, 2018; Liao et al., 2011), and the prefrontal regions are associated with attention regulation of threat distractors and interpretation of potentially threat-related stimuli (Bishop, 2007). The limbic and prefrontal networks are collectively the ‘fear circuitry’, which plays a critical role in SAD (Brühl et al., 2014), in whom the limbic system is hypersensitive to social stimuli (e.g., facial emotion tasks), and under-regulated by the prefrontal network (Brühl et al., 2014; Etkin and Wager, 2007). Altered FC between the two networks is the most consistent finding in SAD, especially between the amygdala and prefrontal areas (Mizzi et al., 2021). This limbic-dorsal prefrontal network disruption may explain this under-regulation (Beaton et al., 2008; Bishop, 2007; Blackford et al., 2014; Schwartz et al., 2012). Moreover, the FC between the limbic and dorsal prefrontal networks reflects treatment response in SAD: decreased FC between amygdala and prefrontal regions was correlated with less symptom reduction following 12-week cognitive behavioral therapy (Klumpp et al., 2014); and oxytocin, a neuropeptide that modulates anxiety, stress and social behaviors, was found to normalize the decreased amygdala-frontal connectivity (Dodhia et al., 2014). In addition, reduced FC between GM4 and GM8 was related to higher HAMD scores, consistent with idea that that limbic-dorsal prefrontal networks could play a role in the failure of emotional regulation which is typical of SAD patients.

Decreased connectivity between the limbic (GM4) and lateral temporal networks (WM6)

Cognitive models propose that biases in information processing may contribute to the maintenance of SAD (Heimberg, 1995; Mellings and Alden, 2000). In previous work socially anxious individuals demonstrated biases in their perceptions of negative valence from faces (Coles et al., 2008; Gutiérrez-García and Calvo, 2017) and exhibited hypervigilance-avoidance in response to negative expressions (Claudino et al., 2019). The temporal lobe helps visual processing and establishing object recognition, and its neurocytes are known to have responding to face images (Rolls, 2000). The limbic network receives perceptual messages from the temporal network and participates in appraising the emotional significance of stimuli and guiding social decisions and social behavior (Baron-Cohen et al., 1994). The temporal network, together with limbic network (amygdala and orbitofrontal gyrus), is predominantly involved in the processing of negative expressions (Iidaka et al., 2001). Previous resting-state fMRI studies have reported decreased FC between the amygdala and temporal gyrus in anxiety disorder patients (Hahn et al., 2011; Jung et al., 2018; Liao et al., 2010b). The temporal network identified in our study mainly involves the inferior longitudinal fasciculus (ILF) extending to the uncinate fasciculus, which interconnects the limbic regions such as the amygdala and OFC (Kiernan, 2012; Latini, 2015; Olson et al., 2007). Furthermore, diffusion tensor imaging studies in SAD have shown reduced fractional anisotropy in ILF (Tükel et al., 2017) and uncinate fasciculus (Baur et al., 2013), which indicates disrupted WM microstructure. We speculate that the reduced FC between the GM of limbic network and WM of lateral temporal network may be related to the disrupted WM microstructure of the lateral temporal network, and may be interpreted as the neural underpinnings of deficit of negative expression processing in SAD patients.

Other decreased connectivity

We also found decreased FC between the limbic network (GM4) and the GM of SMN (GM3), and between the WM of SMN (WM5) and the VIS networks (GM7 and WM8). Reciprocal connections between limbic and perceptual networks are critical pathways associated with anxiety; deficits in either pathway might drive pathological anxiety (Cannistraro and Rauch, 2003). Indeed, the primary motor gyrus may relate to the preparation of motor responses to affective material in the context of the emotional experience (Hardee et al., 2017). Experimentally, stimulating the somatosensory cortex by transcranial magnetic stimulation can affect emotion recognition in social face recognition tasks (Pourtois et al., 2004). Furthermore, the somatosensory network, together with the limbic area, regulates the mode of fear responses to threat-related information (Kropf et al., 2019), and altered FC between the amygdala and sensorimotor processing areas may relate to sensitization to emotion-related cues (Sandman et al., 2020). Resting-state MRI has revealed reduced connectivity between the amygdala and SMN, which indicated a decreased ability to regulate emotion (Pagliaccio et al., 2015). The negative association between the limbic-SMN FC and HAMD scores provides further evidence for dysregulation of emotion in SAD patients during the resting state.

The SMN and VIS represent the perceptual networks and play a critical role in lower-level cognitive processing (Liu et al., 2015). They are together responsible for information communication with the external environment (Pang et al., 2021). Aberrant activity in VIS and SMN has been found both in the resting state and responding to social threat in SAD patients. For example, an inconsistent functional activity pattern, i.e. hyperactivity in VIS and hypoactivity in SMN, has been identified in SAD patient during the resting state (Qiu et al., 2015). In addition, during a social anxiety imagery task, anxiety-related neural responses in VIS and SMN were enhanced in SAD patients after treatment (Kilts et al., 2006). Reduced connectivity between the SMN and VIS areas in our study may reflect aberrant perceptual reactivity responding to social intercourse, consistent with previous studies (Liao et al., 2010a; Liu et al., 2015).

Furthermore, the exploratory correlation analysis found a negative association between illness duration and the FC in SMN-VIS networks, such that longer-duration patients show more reduction in FC. The SMN and VIS play a critical role in motional arousal at the perception level (Kohn et al., 2014; Nummenmaa et al., 2012; Sabatinelli et al., 2013). Previous studies in anxiety disorder found that participants with sensory hypersensitivity showed elevated trait anxiety and state anxiety (Engel-Yeger and Dunn, 2011), and this sensory sensitivity was correlated with duration of current episode (Serafini et al., 2017). We suspect that as the disease progresses, progressive FC damage might be related to gradually impaired visual and sensory function in SAD. However, considering the paucity of data points for longer durations, these results should be regarded with caution. Future work using larger samples and longitudinal designs are needed to investigate the relationships between progressive FC alterations and disease state.

Functional networks correspond to structural white-matter tracts

We clustered several functional networks showing overlap with structural tracts such as corona radiata network, cingulum network, and lateral temporal network (inferior longitudinal fasciculus). We defined several white matter networks according to the spatial anatomy and the correspondence between our white matter networks and the known resting-state grey matter networks (Yeo et al., 2011), such as sensorimotor network, visual network. A previous study measured the similarity between 20 anatomical DTI tracts based on the John Hopkins University (JHU) white-matter tractography atlas and 12 WM functional networks by calculating the percentage of voxels in the functional network identified as belonging to that tract; WM functional networks showed a high overlap with structural tracts, although several functional networks extended across several tracts (Peer et al., 2017). Broadly consistent with this, we clustered some white matter functional networks that overlapped with structural tracts, whereas some white matter networks showed clustering patterns similar to resting-state grey matter networks.

In conclusion, this study explored the interactions between large-scale brain networks in SAD patients, revealing whole brain functional dysconnectivity between limbic network and dorsal prefrontal network, lateral temporal network, and SMN, which were correlated with the severity of clinical symptoms. In addition, the decreased FC between the SMN and VIS correlated with illness duration, suggesting a progressive disturbance of perceptual reactivity to social stimuli in SAD. Furthermore, although the physiological origins of the WM fMRI signal remain unclear, our study uncovered abnormalities of GM-WM connectivity which help extend our understanding of the neurophysiologic basis of SAD.

Limitations of the study

This study has some limitations. First, the relatively small sample size: larger studies are needed to validate the results. Second, the heterogeneous pharmacological profiles (12 SAD patients on different medication): although the subgroup analysis found no FC difference between medicated patients and medication-naïve patients, arguing against significant effects of SSRIs on FCs, the unequal sample sizes (30 patients were medication-naïve, whereas 12 patients were treated) may affect statistical power (Rusticus and Lovato, 2014). Longitudinal studies will be needed to clarify whether FC changes before and after treatment in SAD. Third, the precise physiological mechanism underling the BOLD signal in white matter is still incompletely understood (Gawryluk et al., 2014). Fourth, the precise relationship between white matter functional networks and white matter tracts is unclear, and future studies are needed to define this.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the corresponding author, Qiyong Gong (qiyonggong@hmrrc.org.cn).

Materials availability

The study did not generate new unique reagents.

Experimental model and subject details

Participants

Fifty-one patients (male = 33; female = 18; mean age = 27.8) with social anxiety disorder (SAD) who met the Diagnostic and Statistical Manual of Mental Disorders, 4th ed. (DSM-IV) criteria were recruited (First et al., 1997). Symptoms of fear and avoidance were assessed using the self-administered Liebowitz Social Anxiety Scale (LSAS). Symptoms of anxiety and depression were assessed using Hamilton Anxiety Rating Scale (HAMA) and Hamilton Depression Rating Scale (HAMD). No patients had psychiatric comorbidities. Forty-eight healthy controls (HC) (male = 31; female = 17; mean age = 24.3) were recruited by poster advertisements from the local area, and the SCID non-patient edition was used to confirm the lifetime absence of psychiatric illness in both HC subjects and their first-degree relatives. Exclusion criteria for both groups, assessed by two experienced psychiatrists in a detailed clinical interview, were: (1) Age <18 or >60 years; history of (2) head injury, (3) alcohol or substance abuse, (4) neurologic illness, (5) serious medical or surgical illness; and (6) current pregnancy and (7) claustrophobia or other contraindications to magnetic resonance imaging (MRI) examination. All participants were right-handed.

This study was approved by the ethics committee of West China Hospital, and written informed consent was obtained from all participants.

Method details

MRI acquisition

All scans were performed at West China Hospital in a 3.0 T magnetic resonance scanner (Siemens Trio, Erlangen, Germany) with an 8-channel head coil. The participants lay supine. Foam pads were used to stabilize the head, and soft earplugs to reduce noise. Resting-state fMRI images were obtained using the following gradient-echo planar imaging sequence: repetition time (TR) 2000 ms; echo time (TE) 30 ms; flip angle 90°; field of view (FOV) 240 × 240 mm2; data matrix 64 × 64; voxel size 3.75 × 3.75 × 5 mm3; slice thickness 5 mm; 205 volumes. High-resolution T1-weighted images were acquired using a spoiled gradient recalled sequence: TR 1900 ms; TE 2.26 ms; TI 900 ms; flip angle 9°; FOV 256 × 256 mm2; data matrix 256 × 256; voxel size 1 × 1 × 1 mm3; 176 sagittal slices with thickness 1 mm. During the scan, subjects were instructed to keep their eyes closed, not fall asleep or think of anything in particular, and keep as still as possible.

MRI preprocessing

Preprocessing used SPM8 (http://www.fil.ion.ucl.ac.uk/spm) and DPARSF (http://rfmri.org/DPARSF) (Chao-Gan and Yu-Feng, 2010) in MATLAB (MathWorks, Natick, MA). High-resolution T1-weighted images were reoriented, skull-stripped and segmented into GM, WM, and cerebrospinal fluid (CSF). Resting-state fMRI data was preprocessed in the following steps: removal of the first 10 time-points; slice-time correction; head motion correction (subjects with head motion >2.5 mm and >2.5⁰ degrees were excluded, and motion “spikes” with framewise displacement >1 mm were removed as separate regressors); removal of linear trend and band-pass filtering (0.01–0.15 Hz); multiple regression to remove nuisance variance (e.g., 24-parameter motion correction and mean CSF signals), excluding the WM and global brain signals to retain the signals of interest; spatial normalization to Montreal Neurological Institute (MNI) space and resampling to 3 × 3 × 3mm3; spatial smoothing, separately in the WM and GM, with 4 mm full-width half-maximum (FWHM).

Functional network clustering

Network clustering was performed using code by Peer et al. (2017). Briefly, the segmentation results of each subject’s anatomical images were used to create group-level WM and GM masks. For each subject, we divided each voxel into GM, WM, or CSF based on the maximum probability from the three segmentation images, resulting in individual WM, GM, and CSF masks. The individual masks were averaged to obtain the percentage of subjects for each voxel in which it was classified as WM or GM. Voxels identified as WM in > 60% of subjects were used to create the WM mask. Voxels identified as GM in > 20% of subjects were used to create the GM mask which aims to avoid missing cortical regions with high variability, but exclude any voxels included in the white matter mask. Thus, the mask contains almost all GM voxels, even in cortical regions with high variability. Harvard-Oxford Atlas was used to removed subcortical areas from the WM mask to classify the deep brain structures correctly. This finally yielded two group-level masks: a WM mask with 18,649 voxels and a GM mask with 42,859 voxels.

Next, white matter and grey matter functional networks were identified by a clustering method. First, in order to reduce the computational complexity of the clustering, voxels in the masks were subsampled using an interchanging grid method. The key element of this method is extracting any second voxel along the image rows and columns, and shifting it by 1 between slices to avoid missing entire columns of data (Peer et al., 2017), resulting in 4,729 nodes in the WM mask and 4,637 nodes in the GM mask. Next, Pearson’s correlation analysis was performed between each WM voxel and each subsampled node to generate a group-level WM correlation matrix (18,649 × 4,729); an analogous process yielded a 42,859 × 4,637 GM correlation matrix. Then the group-level correlation matrices were randomly divided into 4-folds (18,649 × 1182 matrix per fold in WM and 42,859 × 1159 per fold in GM). K-means clustering (distance metric-correlation, 10 replicates) was employed to classify voxels with similar connectivity patterns, with numbers of clusters (K) ranging from 2-22. Note that clustering was performed on each fold of each number of clusters.

Finally, Dice coefficient was computed to assess the stability of the clustering results. For each clustering number (from 2 to 22), if clustering results in all 4-folds are roughly similar, clustering solutions can be considered stable. Therefore, to evaluate the similarity between clustering solutions of 4-folds, an adjacency matrix was computed for each fold, and these adjacency matrices were compared using Dice coefficient, obtaining an average Dice coefficient for each clustering number. We give preference to the maximum number of clusters with good stability (Dice coefficient > 0.9), because the largest number of clusters would offer the most detailed level of description of functional networks (Jiang et al., 2019b; Peer et al., 2017). Functional networks constructed by above methods have been shown to be highly symmetrical and valid (Peer et al., 2017).

Functional network construction

For each participant, each functional network’s mean time series (obtained by averaging the fMRI time series across all voxels in each network) were extracted from the clustering results. Next, Pearson’s correlation coefficients were computed, network-by-network. Fisher’s z transformation was applied to improve the normality of correlation coefficients. Finally, three FC matrices were constructed for each subject: between WM and WM networks, between GM and GM networks, and between WM and GM networks.

Quantification and statistical analysis

To detect significant FC within each group, we performed one-sample t-tests with Bonferroni correction (p < 0.05) for multiple comparisons. To detect significant differences of FC between SAD patients and HC, we used two-sample t-tests with Bonferroni correction (p < 0.05) for multiple comparisons. We then extracted the FC values with group differences, and in SAD group performed Pearson correlation analyses against age, illness duration, and symptom severity scores. To examine the potential effects of medication, we conducted a subgroup analysis between medication-naïve and medicated patients. To compare demographic data between the two groups we used two-sample t-tests for age and chi-square tests for gender, implemented in SPSS software.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
MRI data and code	This paper	Available from the corresponding author upon reasonable request
Software and algorithms
MATLAB	MathWorks	https://www.mathworks.com/
SPM8	FIL	http://www.fil.ion.ucl.ac.uk/spm
DPARSF	RFMRI.ORG	http://rfmri.org/DPARSF