Altered brain functional network dynamics in obsessive-compulsive disorder.

Obsessive-compulsive disorder (OCD) is a debilitating and disabling neuropsychiatric disorder, whose neurobiological basis remains unclear. Although traditional static resting-state magnetic resonance imaging (rfMRI) studies have found aberrant functional connectivity (FC) in OCD, alterations in whole-brain FC and topological properties in the context of brain dynamics remain relatively unexplored. The rfMRI data of 29 patients with OCD and 40 healthy controls were analyzed using group independent component analysis to obtain independent components (ICs) and a sliding-window approach to generate dynamic functional connectivity (dFC) matrices. dFC patterns were clustered into three reoccurring states, and state transition metrics were obtained. Then, graph-theory methods were applied to dFC matrices to calculate the variability of network topological organization. The occurrence of a state (State 1) with the highest modularity index and lowest mean FC between networks was increased significantly in OCD, and the fractional time in brain State 1 was positively correlated with anxiety level in patients. State 1 was characterized by having positive connections within default mode (DMN) and salience networks (SAN), and negative coupling between the two networks. Additionally, ICs belonging to DMN and SAN showed lower temporal variability of nodal degree centrality and efficiency in patients, which was related to longer illness duration and higher current obsession ratings. Our results provide evidence of clinically relevant aberrant dynamic brain activity in OCD. Increased functional segregation among networks and impaired functional flexibility in connections among brain regions in DMN and SAN may play important roles in the neuropathology of OCD.

INTRODUCTION

Obsessive–compulsive disorder (OCD) is a debilitating and disabling neuropsychiatric disorder characterized by intrusive thoughts and repetitive or ritualistic behaviors (Goodman, Grice, Lapidus, & Coffey, 2014; Li et al., 2020). It has a lifetime prevalence of 2–3% in the general population (Abramowitz, Taylor, & McKay, 2009; Li et al., 2011). As an emerging subspecialty of radiology, psychoradiology has led to our better understanding of the complex brain alterations in patients with OCD (Huang et al., 2019), where magnetic resonance imaging (MRI) studies have provided evidence of structural and functional alterations in fronto‐striatal circuitry, notably abnormal gray matter volume in the anterior cingulate cortex, prefrontal cortex, striatum, and thalamus (Milad & Rauch, 2012), and dysconnectivity within and between frontal and subcortical regions (Endrass & Ullsperger, 2014; Li et al., 2014; Milad & Rauch, 2012).

Recent investigations of large‐scale brain networks using resting‐state functional MRI (rfMRI) found extensive alterations in brain functional connectivity (FC) in OCD (Gursel, Avram, Sorg, Brandl, & Koch, 2018). These findings generally support the “triple‐network” model in which abnormal FC patterns within and between the default mode network (DMN), executive control network (ECN), and salience network (SAN) represent core pathophysiological mechanisms of OCD (Fan et al., 2017; Posner et al., 2017; Stern, Fitzgerald, Welsh, Abelson, & Taylor, 2012). Graph‐theoretical analysis of the brain connectome has demonstrated decreased global metrics such as small‐worldness, as well as altered regional metrics such as clustering coefficient and nodal degree in frontal and subcortical regions (Armstrong et al., 2016; Beucke et al., 2013; Zhang et al., 2011), improving understanding of the topological organization of functional brain networks in OCD.

Conventional rfMRI FC analysis treats brain connectivity patterns as if they are stationary during scan sessions. However, the brain is a highly dynamic system with non‐stationary neural activities and rapidly‐changing neural interactions (Fu et al., 2018), and this information is lost in static FC analytic approaches. Although the static analysis approach has provided valuable information concerning brain functional organization in OCD, the dynamic functional connectivity (dFC) approach can provide a complementary understanding of the dynamic changes in functional activity in large scale networks and their alteration in neuropsychiatric illnesses (Hutchison et al., 2013). Neural dynamics is especially important during resting periods where the mental activity is not directly modulated by task demands (Fu et al., 2018). dFC fluctuations are related to changes in cognitive and affective processes (Thompson et al., 2013), and altered functional dynamics of the brain have been associated with changes in cognitive and emotional functions whose alterations characterize psychiatric disorders (Anticevic, Schleifer, & Youngsun, 2015). In some cases, dFC features have provided a higher predictive accuracy for diagnostic and prognostic purposes compared with static FC metrics (Rashid et al., 2016).

Several approaches have been developed to examine the dynamic properties of FC. The sliding‐window method is most widely performed (Preti, Bolton, & Van De Ville, 2017), in which, instead of generating an average correlation between brain regions across the entire period of data acquisition, a separate FC matrix for each time window is calculated. Then, the variability of FC over time windows is examined to index dFC. These changes in brain connectivity states are systematic, alternating between discrete states that can be identified statistically using cluster analysis. The potential biomarker utility and clinical relevance of dFC patterns have been demonstrated in several psychiatric disorders, including schizophrenia and bipolar disorder (Du et al., 2017; Zhang et al., 2018), autism (Yao et al., 2016), and major depression (Wu et al., 2019). The graph theory based approach applied to dFC has further clarified the whole brain and nodal patterns of dynamic brain connectivity alteration in neuropsychiatric disorders (Yu et al., 2015). Measures of FC between pairs of brain regions are necessarily constrained by the topology of the whole brain, which has not been taken into account in typical FC analysis focused on a sequential comparison of pairs of brain regions (Shine & Poldrack, 2018).

Alterations in whole‐brain FC and topological properties in the context of brain dynamics remain relatively unexplored in OCD. Gursel et al. using the sliding‐window approach found peak dysconnectivity (i.e., peak reduction in correlation between two networks across time windows) between bilateral fronto‐parietal networks (FPN) and between the left FPN and the SAN in patients with OCD (Gursel et al., 2020). This study focused on the peak FC impairment of specific networks (FPN, DMN, and SAN) rather than the whole brain, and other key indices of dFC were not obtained. Liu et al. used spatial group independent component analysis (GICA) and the sliding‐window approach to assess whole‐brain dFC, and found an altered number of transitions among different dynamic states in OCD patients that correlated with obsessive–compulsive symptom severity (Liu et al., 2020). This study used a low model order GICA and did not explore alterations of dynamic topological properties. Another study found that dynamic topological characteristics of a repetitive transcranial magnetic stimulation (rTMS) target node (dorsolateral prefrontal cortex) at baseline were positively associated with the reductions in distress after rTMS in patients with OCD (Douw et al., 2020). However, their study focused on the target node of rTMS rather than the whole brain network and did not explore differences in related dynamic global and nodal topological metrics relative to healthy individuals. Investigation of whole brain resting‐state dynamic properties may allow the identification of altered functional brain dynamics and their topological metrics in patients with OCD.

In the present study, we used rfMRI and a sliding‐window approach to compare whole brain dynamics in patients with OCD and healthy controls. Analyses of dFC states and their dynamic global and nodal topological organization were used to evaluate brain dynamics. We hypothesized that patients with OCD would show abnormal temporal properties of dFC states and altered variability in network topological indices in brain regions of DMN and SAN, and that altered functional brain dynamic metrics would be related to the OCD symptom severity.

MATERIALS AND METHODS

Participants

This study was approved by the local research ethics committee of West China Hospital of Sichuan University and informed written consent was obtained from all participants. Diagnosis of OCD was determined by experienced clinical psychiatrists using the Structured Clinical Interview for DSM‐IV (SCID) (Patient Edition). We recruited 29 patients with OCD, 15 drug‐naive and 14 who had previously received psychiatric medications (four clomipramine, four paroxetine, three fluoxetine, and three sertraline) but had been medication free for more than 2 weeks. Treatment was not controlled for this study; prior treatment discontinuation was based on independent physician advice or patient decisions to terminate treatment. The predominant obsessions/compulsions of patients with OCD were determined according to the five clinical dimensions defined by Mataix‐Cols, Rauch, Manzo, Jenike, and Baer (1999)): 21 aggressive/checking, 5 contamination/cleaning, 2 symmetry/ordering, 1 sexual/religious, and none with hoarding symptoms. Clinical symptoms of patients with OCD were assessed using the Yale‐Brown obsessive–compulsive scale (Y‐BOCS), 14‐item Hamilton anxiety rating scale (HARS), and 17‐item Hamilton depression rating scale (HDRS) (Table 1).

TABLE 1

Demographic and clinical characteristics of study participants

	OCD (n = 29)	HCs (n = 40)	χ2/t value	p‐value
Gender (number)	19M, 10F	25M, 15F	0.066	.797
Age (years)	27.8 ± 9.4 (18–52)	27.9 ± 9.2 (18–52)	−0.032	.975
Education (years)	13.9 ± 2.9 (8–19)	12.7 ± 3.5 (5–19)	1.422	.160
Duration (years)	6.2 ± 5.5 (1–23)	–	–	–
Total Y‐BOCS score	23.0 ± 5.4 (16–33)	–	–	–
Obsessive subscale score	17.2 ± 4.6 (10–28)	–	–	–
Compulsive subscale score	5.9 ± 5.8 (0–16)	–	–	–
HARS	7.7 ± 3.0 (3–19)	–	–	–
HDRS	9.7 ± 2.8 (5–17)	–	–	–

Note: Values were given as mean ± SD (range). p value of gender was obtained by chi‐square test and p values of age and education years were obtained by two‐sample t test.

Abbreviations: F, female; HARS, Hamilton anxiety rating scale; HCs, healthy controls; HDRS, Hamilton depression rating scale; M, male; OCD, obsessive–compulsive disorder; Y‐BOCS, Yale‐Brown obsessive–compulsive scale.

Forty healthy controls (HCs) were recruited from the local region through poster advertisements and screened using the SCID‐I/Non‐patient Edition to identify any history of psychiatric illness. HCs had no known history of psychiatric illness among their first‐degree relatives. Exclusion criteria for all participants included: (a) age younger than 18 years or older than 60 years; (b) any history of OCD, neurodevelopmental, affective or psychotic illness; (c) significant systemic or neurologic illness or neurosurgery; (d) history of substance abuse or dependence disorders; and (e) pregnancy and other MR contraindications. MR images were inspected by two experienced neuroradiologists to confirm the absence of gross brain abnormalities. The OCD and HC groups did not differ significantly with respect to age (p = .975), gender (p = .797), or years of education (p = .160). All participants were right‐handed (Table 1). Some of the rfMRI data used in the present study have been used in a recently published paper focusing on altered effective anterior cingulate connectivity in OCD (Long et al., 2020).

MRI acquisition

High‐resolution three‐dimensional structural MRI and rfMRI data were collected using a 3.0 T MRI system (EXCITE, General Electric) with an eight channel phased array head coil. High‐resolution T1‐weighted scans were acquired using a spoiled gradient recall sequence with the following parameters: repetition time/echo time (TR/TE) = 8.5 ms/3.4 ms, flip angle = 12°, slice thickness = 1.0 mm, field of view = 240 × 240 mm2 with an acquisition matrix comprising 256 readings of 128 phase encoding steps that produced 156 contiguous coronal slices. Structural MRI images were automatically interpolated in‐plane to 512 × 512, which yielded an in‐plane resolution of 0.47 × 0.47 mm2. The rfMRI data were obtained via a gradient‐echo echo‐planar imaging (EPI) sequence (TR/TE = 2000 ms/30 ms, flip angle = 90°, slice thickness = 5 mm with no gap, field of view = 240 × 240 mm2, matrix size = 64 × 64, voxel size = 3.75 × 3.75 × 5 mm3. Each brain volume comprised 30 axial slices to cover the whole brain, and 200 volumes preceded by five dummy volumes (without data generated) collected for shimming to eliminate field inhomogeneities, leading to a total scanning time of 410 s. During scanning, subjects were instructed to keep their head still and relax with their eyes closed without falling asleep or systematic thought. Earplugs and foam padding were used to reduce noise and head motion.

Data preprocessing

Preprocessing was performed using the toolbox for Data Processing & Analysis of Brain Imaging (rfmri.org/DPABI). The first 10 volumes of rfMRI data were removed to reduce equilibration effects, leaving a total of 190 volumes for statistical analysis. The remaining functional images underwent slice‐timing correction and were realigned to reduce displacement between volumes. Individual structural images were co‐registered to the mean of functional images after realignment. The transformed structural images were then segmented into gray matter, white matter, and cerebrospinal fluid. The diffeomorphic anatomical registration through exponentiated lie algebra tool was used to compute transformations from individual native space to the Montreal Neurological Institute (MNI) template and to resample functional images to 3 × 3 × 3 mm3 resolution. Then, the normalized data were spatially smoothed using a 6 mm full‐width at half‐maximum Gaussian kernel. Participants with a maximum head displacement more than 1.5 mm, maximum rotation greater than 1.5°, or mean framewise displacement (FD) (Jenkinson) larger than 0.2 mm were excluded from the analysis. No participant was excluded according to these head motion criteria. FD did not differ between the two groups (0.043 ± 0.031 for patients, 0.053 ± 0.026 for controls, t = 1.48, p = .14).

GICA analysis and identification of independent components

The flowchart for dFC state analysis and dynamic topological analysis is shown in Figure 1. To identify intrinsic connectivity networks, we performed spatial GICA using the GIFT toolbox (mialab.mrn.org/software/gift). GICA is a data‐driven approach to decompose rfMRI data into functionally homogeneous regions (Kiviniemi et al., 2009), enabling a whole‐brain analysis without resorting to atlas‐defined regions of interest that may merge or imprecisely delineate functionally distinct areas (Allen, Erhardt, Wei, Eichele, & Calhoun, 2012). ICA‐based studies have identified independent components (ICs) that resemble established functional cortical and subcortical networks (Damoiseaux et al., 2006; Zuo et al., 2010). Before ICA, principal component analysis for subject‐specific data reduction was performed, in which 150 principal components were identified using standard economy‐size decomposition. In data reduction for group comparisons, the concatenated subject‐reduced data were decomposed into 100 ICs.

FIGURE 1

Flowchart of dynamic functional connectivity (FC) state analysis and dynamic topological analysis. The steps included: (a) 100 independent components (ICs) were obtained by group independent component analysis (GICA), then 51 of 100 ICs were characterized as meaningful and assigned to eight brain networks; (b) a sliding‐window approach was used to segment the whole scan series into consecutive windows, and the FC covariance matrices of the 51 ICs were computed in each window; (c) k‐means clustering was used to identify three dFC states, and three state transition metrics for these states were calculated for each subject; and (d) in graph theory analysis, each FC matrix in each window was binarized using a series of sparsity thresholds (from 0.10 to 0.37, with an interval of 0.01), and then the coefficient of variation (CV) of the area under the curve (AUC) of both global and nodal graph metrics was calculated across all windows for all study participants

Reliability and stability of the Informax GICA algorithm in ICASSO (Himberg, Hyvarinen, & Esposito, 2004) was ensured by repeating the algorithm 20 times, and using the most central run to reconstruct subject‐specific time courses and spatial maps of each IC using the GICA back reconstruction algorithm (Calhoun, Adali, Pearlson, & Pekar, 2001). The group ICs of the 20 runs were clustered to estimate their reliability, and ICs with a quality index (I q) more than 0.7 were selected (Ma et al., 2011). Using a one sample t test across all subjects for each IC, we obtained the t‐map of each IC, which was thresholded with t > mean (μ) + 4SD (σ) (Allen et al., 2011). The μ in the formula represents the mean value of all voxels in the t‐map for each IC, and σ representing the SD of the values for each IC. Labels and spatial maps of each IC are presented in Figure S1, and the peak coordinates of ICs are provided in Table S1.

We characterized 51 of 100 ICs as meaningful based on the following criteria (Allen et al., 2014): peak coordinates of spatial maps located in gray matter, with minimal spatial overlap with white matter, vessels, ventricles, or susceptibility artifacts, and time courses characterized by a high dynamic range. To provide a more informative description of the organization and primary spatial features of the 51 ICs, we then sorted the selected 51 ICs into different functional networks to classify the ICs. The ICs can be sorted either spatially or temporally. We chose spatially sorting in the present study by comparing the IC's spatial image with a template. This step was implemented by using the “sort components” tool of the GIFT toolbox and the Stanford functional ROI template (findlab.stanford.edu/functional_ROIs.html) (Shirer, Ryali, Rykhlevskaia, Menon, & Greicius, 2012) as the reference to sort the 51 ICs into different functional networks. There are four types of sorting criteria, including multiple regression, correlation, kurtosis, and maximum voxel. Among them, the multiple regression criterion is most usually used, which can select several template regions as a whole.

The Stanford template includes 90 ROIs that belong to eight predefined networks. We performed multiple regression eight times, each time selecting all ROIs belonging to one network of the template as a whole, and spatial map intensity of all 51 ICs was selected as regressors for each network of the template. Then we obtained the coefficient of determination (R 2) of each IC's spatial map intensity for the spatial distribution of every network of the template (Table S1). R 2 reflects the proportion of variance in the dependent variable that was predictable from the independent variable. Therefore, a higher R 2 indicates a higher similarity between an IC and a specific network of the template. Each IC was assigned to the network to which it had the highest R 2. We also visually confirmed that IC linkages to brain networks by establishing that the peak coordinates of ICs were located within the template for this network. This procedure guided assigning each of the 51 ICs into one of the eight functional brain networks for the subsequent analyses: auditory network (AN), DMN, ECN, language network (LAN), SAN, subcortical network (SC), sensorimotor network (SMN), and visual network (VN) (Figure 2).

FIGURE 2

Fifty‐one independent components (ICs) were divided into eight functional networks for all participants in the present study. The networks included the auditory network (AN), default mode network (DMN), executive control network (ECN), language network (LAN), salience network (SAN), subcortical network (SC), sensorimotor network (SMN), and visual network (VN)

Computation of dFC

dFC was examined using the temporal dFNC toolbox in GIFT. Before computation for windowed matrices, additional postprocessing steps were performed for time courses of all ICs to regress out the influence of noise sources including: (a) detrending (regressing linear, quadratic, and cubic trends); (b) despiking using 3D DESPIKE; (c) low‐pass filtering using a high‐frequency cut‐off of 0.15 Hz; and (d) reducing influence of six parameters of head movement using regression.

A sliding‐window approach was used to explore the time‐varying changes of FC within the 51 ICs identified during rfMRI scans. We chose a 22‐TR window (44s) because previous studies have suggested that windows of 30 to 60s can successfully capture patterns of resting‐state fluctuations of dFC (Preti et al., 2017). We used a Gaussian (σ = 3 TRs) function to create a tapered window, slid step‐wise by one TR along the scan image time series, and then computed the 51 × 51 pairwise FC matrix by Pearson's correlations in every sliding window. To promote sparsity in estimations, the L1 norm penalty was implemented in the LASSO framework with 50 repetitions (Friedman, Hastie, & Tibshirani, 2008). Then the correlation values of the 51 × 51 pairwise functional matrices were converted to z‐values with Fisher's z‐transformation to improve normality and comparability. We removed variance from the windowed FC correlations related to age, sex, education years, and mean FD for each participant.

dFC states analysis

To assess reoccurring dFC patterns, k‐means clustering was performed on the 51 × 51 FC matrices of all sliding windows for all participants. The k value was varied from two to 10 to identify the optimal value. To reduce the potential bias of the initial random selection of cluster centroids, the k‐means clustering algorithm was iterated 50 times for each k value. In this procedure, we used the Manhattan distance (Aggarwal, Hinneburg, & Keim, 2001) to measure the similarity of FC matrices between time windows. Then we performed an analysis for cluster number validity (using the best run of each 50 iterations under different k values and Dunn's index) (Havens, Bezdek, Keller, & Popescu, 2008), which determined the optimal k value to be three in the present study (Figure S2). The best run across the 50 iterations of k = 3 was kept for further analysis. We used the three cluster centroids of all study participants to represent three reoccurring FC states (Figure 3a). For visualization of dFC patterns in the two groups, respectively, we calculated the group‐specific centroid of each state by averaging subject‐specific centroids of each group (Figure 4a).

FIGURE 3

Cluster centroids and characteristics of each dynamic functional connectivity (FC) state under 22‐TR window size for all participants. (a) Cluster centroids for each state. (b) Modular distribution for each state. (c) The radar map and line graph of the mean FC strength within and between networks for all three states

FIGURE 4

Centroids of dynamic functional connectivity (FC) states and connections with the top 5% in FC strength in patients with obsessive–compulsive disorder (OCD) and healthy controls (HC) under the window size of 22‐TR

As the interpretation of the meaning of dFC states is essential in dynamic studies, we described two aspects of the three dFC states: the global integration level of the 51 ICs and the connection characteristics of the eight functional networks. First, we calculated the modularity index Q (Yue et al., 2017) of 51 ICs in each state to describe the integration level of each state. Modularity is defined as the ability of a graph to be subdivided into modules that are maximally connected within modules and sparsely connected between modules (Newman & Girvan, 2004), which thus measures both integration within modules and separation between modules. The Brain Connectivity Toolbox (www.brainconnectivity-toolbox.net) was used to calculate the modularity index Q using a normal Louvain community detection algorithm. A larger Q value indicates that the 51 functional ICs in the brain are more likely to aggregate into specific modules (Wu et al., 2019). Further, to describe the connection characteristics of each functional network formed by disparate ICs, we calculated the mean FC value within and between the eight networks investigated in the present study for each state using the cluster centroid.

To quantify temporal properties of the FC states, we assessed three state transition metrics (Kim et al., 2017): (a) fractional time, which represents the proportion of time spent in each state as measured by percentage; (b) mean dwell time, which represents the average duration of time intervals spent in each state before switching into other states; and (c) number of transitions, representing the number of switches between states throughout the scan acquisition. Nonparametric permutation tests (10,000 iterations) were used to assess differences of the fractional time, mean dwell time in each state and number of transitions between groups. Because of the wide age range of participants in our study, we also tested whether there were age‐by‐diagnosis interactions by using a general linear modal on the state transition metrics between groups. To test for the possible effect of prior psychopharmacological treatment, we compared the differences of three state transition metrics between drug‐naive and previously medication‐treated patients. Significance was set at p < .05 and the false discovery rate (FDR) correction was used for all group comparisons.

To evaluate the consistency and validity of the k‐means clustering at different window sizes, we used a 30‐TR window and repeated the above dFC states analysis. We calculated Pearson's correlation coefficients between the cluster centroids under the two different window sizes, and the states with the highest correlation coefficient were defined as the same state as in the primary analysis (Wu et al., 2019).

Dynamic topological metrics analysis

To explore the dynamic topological organization of the resting‐state functional brain network, we applied a graph theory approach to examining topological metrics across all sliding windows using GRETNA software (www.nitrc.org/projects/gretna). There are many ways to define nodes for the analysis of a topological network, something that has often been done using a brain atlas. Alternatively, ICA is a data driven approach that reveals ICs in the data (represented by spatial maps of regions with associated time courses), which avoids assumptions about brain maps and region borders defined primarily on anatomy rather than function. The inter‐relationship of identified ICs can be examined with graph theory approaches (Smith et al., 2011). Using this strategy, in the present study, the 51 ICs were defined as functionally independent nodes and the FC between pairs of ICs as edges.

To achieve this aim, we first binarized the FC matrices of all windows with a series of sparsity thresholds. In graph theory, sparsity is defined as the ratio of the number of existing edges divided by the maximum possible number of edges in a network. This approach normalized networks to have the same number of edges by applying a subject‐specific FC threshold and minimized possible discrepancies in estimation of FC strength (He et al., 2009). Instead of selecting a single threshold, we thresholded each FC matrix repeatedly over a wide range of sparsity levels, and computed the area under the curve (AUC) for each network metric across sparsity thresholds. This provided an overall scalar value for the topological characterization of brain networks independent of any specific sparsity threshold. The selection of sparsity thresholds was based on the following criteria: (a) the average degree (the degree of a node is the number of edges connected to the node) over all nodes under a specific sparsity threshold was larger than 2 × log(N) (He, Chen, & Evans, 2007), where N represents the number of nodes (N = 51 in the present study); and (b) the small‐worldness scalar sigma (σ) of the thresholded network was larger than 1.1 for all subjects (Zhang et al., 2011) (Figure S3). The generated thresholds ranged from 0.10 to 0.37 (with an interval of 0.01) for the subsequent analyses.

At each sparsity threshold, we calculated both global and regional network metrics in each 51 × 51 pairwise FC matrix of all sliding windows for all subjects (Zhang, Wang, Wu, et al., 2011). Global measures included: (a) small‐world global metrics of clustering coefficient (C), normalized clustering coefficient (γ), characteristic path length (L), normalized characteristic path length (λ), and small‐worldness (σ); and (b) measures of global (E ) and local network efficiency (E ). The regional properties were nodal degree centrality and nodal efficiency. The interpretation of these topological network measures is presented in Table S2. For the dynamic networks, we calculated the coefficient of variation (CV = SD/mean) of AUC of network parameters across all sliding windows to assess the variability of topological metrics over time. The CV is useful for comparing the variability of features that are measured on different units or that exhibit different means (Pauly & Smaga, 2020).

The nonparametric permutation approach (10,000 iterations) was used to test for dynamic topological property differences in the AUC of each metric (small‐world metrics, network efficiency, and nodal degree, centrality, and efficiency) between groups. To test for statistical significance of group differences in the temporal variability of global network properties, we set the significance at p < .05 with FDR correction for these multiple comparisons. To test for significance of group differences in the larger number of nodal properties, a false positive adjustment approach was used to control the false‐positive rate to one per analysis (Fornito, Yoon, Zalesky, Bullmore, & Carter, 2011; Lynall et al., 2010). We also tested whether there were age‐by‐diagnosis interactions in dynamic topological metrics. To test for the possible effect of prior psychopharmacological treatment, we compared the differences in global and nodal dynamic topological metrics between drug‐naive and previously medication‐treated patients.

Correlational analyses

Because of the nonnormality of dynamic measures, we performed Spearman's partial correlation analyses to examine relations of altered network properties (state transition metrics and CV of topological metrics) with illness duration and OCD symptom severity, treating head motion (FD Jenkinson) and demographics (age, sex, and education years) as covariates. Nominal significance was set at p < .05 for these exploratory analyses conducted for descriptive and heuristic purposes.

RESULTS

General temporal properties of dFC states

Using k‐means clustering, we identified three dFC states in analyses of all study participants: State 1 occurred in a moderate frequency (31%), State 2 occurred at a low frequency (14%), and State 3 occurred most frequently (55%). The correlation matrix of centroids in each state is shown in Figure 3a.

The three states showed different characteristics with respect to their modularity. The modularity analysis (Figure 3b) showed that State 1 had the highest Q (0.421), and ICs in State 1 aggregated primarily into two functional modules, a module mainly involving ICs in AN, ECN, and SAN, and another module mainly involving ICs in DMN, LAN, SC, SMN, and VN. State 2 had the lowest Q (0.148). ICs in State 2 also tended to aggregate into two modules, one mainly involving ICs in AN, SAN, SMN, and VN, and another mainly involving ICs in DMN, ECN, LAN, and SC. State 3 had a moderate Q (0.246). ICs in State 3 aggregated into four modules, module 1 mainly involving ICs in ECN and SC, module 2 mainly having ICs in AN and SAN, module 3 mainly involving ICs in DMN, SMN, and VN, and module 4 mainly involving ICs in DMN and LAN.

The three states also show different mean FC within and between eight functional networks, and the mean FC within networks in each state was higher than the mean FC between networks (Figure 3c and Table S3). State 1 and State 3 showed similar mean within‐network FC, while mean within‐network FC was the highest in State 2. The mean FC between a network and all other networks was weakest in state 1, the strongest in State 2, and moderate in State 3. In State 1 using a one sample t test, there were significant but modest negative between‐network FC effects between DMN and other networks such as SAN (predominantly, mean FC = −0.029, 95% confidence intervals [CI] = −0.030 ~ −0.027, t = −36.302, p < .001) and ECN (mean FC = −0.017, 95% CI = −0.020 ~ −0.015, t = −15.988, p < .001).

Group differences in temporal properties of dFC states

Patients with OCD and controls showed similar cluster centroids. However, there were patient‐control differences in key dFC features. In patients with OCD, the total occurrences of States 1, 2, and 3 were 41, 11, and 48%, respectively. In HC, the total occurrences of the three States were 24, 16, and 60%, respectively, which differed significantly from findings in patients (χ2 = 387.11, p < .001) (Figure 4a). In State 1, the top 5% FC with the strongest connections were mainly located within DMN and SAN with positive coupling, and between DMN and SAN (predominantly) or between DMN and other task‐positive networks such as ECN, SAN, or LAN with negative coupling. In State 2, the top 5% of FC metrics were scattered between many networks with positive coupling. In State 3, the top 5% FC were distributed within and between many networks with positive coupling (Figure 4b).

For all analyses, p values are reported with correction for multiple comparisons, except where uncorrected p values were presented as noted for specific exploratory studies of clinical or demographic correlations with brain features. For state transition metrics, patients with OCD showed both greater fractional time (p = .034, 95% CI = 0.049 ~ 0.338, Cohen's d = 0.611, Power = 0.694) and mean dwell time (p = .017, 95% CI = 3.706 ~ 28.197, Cohen's d = 0.752, Power = 0.860) in State 1 than controls; in States 2 and 3, there were no differences between OCD patients and controls (Figure 5a,b). We did not find a significant difference in the number of state transitions between the two groups (Figure 5c). Age‐by‐diagnosis interactions on all altered state transition metrics were not significant (Table S4). Subgroup analysis showed no differences between drug‐naive and previously medication‐treated patients in altered state transition metrics (Table S5). In patients with OCD, fractional time of State 1 was positively correlated with HARS scores (r = .487, uncorrected p = .017, 95% CI = 0.101 ~ 0.760, Power = 0.825) (Figure 5d). Other state transition metrics were not significantly correlated with clinical scales or duration of illness.

FIGURE 5

State transition vectors in patients and controls, and the correlation between state transition vectors and clinical symptom ratings, in patients with obsessive–compulsive disorder. * indicated p‐value < 0.05

In validation analysis with the window size set to 30‐TR with other parameters unchanged, three dFC states for all the subjects were revealed (Figure S4A). State 2 under 30‐TR window size and State 1 under 22‐TR window size (r = .962, p < .001), state 1 under 30‐TR window size and State 2 under 22‐TR window size (r = .975, p < .001, and State 3 under both window sizes (r = .976, p < .001) provided similar characterization of dFC states (Table S6). The statistical comparison of state transition metrics showed the same between‐group differences under both window sizes (Figure S4b–d) and the fractional time of State 2 under 30‐TR window size (corresponding to State 1 under 22‐TR window size) was positively correlated with HARS scores (r = .492, uncorrected p = .012, 95% CI = 0.036 ~ 0.739, Power = 0.835) (Figure S4e).

Dynamic topological metrics

Comparison of dynamic topological metrics between patients with OCD and controls failed to identify significant group differences in CV of AUC of small‐world metrics (for σ, p = .922; for γ, p = .700; for λ, p = .367) or network efficiency (for E , p = .780; for E , p = .219). For temporal variability of nodal metrics, patients showed decreased variability of degree centrality in IC28 (p = .004, 95% CI = − 0.100 ~ − 0.026, Cohen's d = 0.659, Power = 0.759), IC 36 (p = .006, 95% CI = − 0.117 ~ − 0.041, Cohen's d = 0.778, Power = 0.882), and IC 57 (p = .019, 95% CI = − 0.086 ~ − 0.021, Cohen's d = 0.656, Power = 0.755), as well as decreased variability of efficiency in IC 28 (p = .004, 95% CI = − 0.055 ~ − 0.014, Cohen's d = 0.843, Power = 0.926) and IC 36 (p = .005, 95% CI = − 0.058 ~ − 0.020, Cohen's d = 0.795, Power = 0.895) compared with controls (Figure 6a). IC 28 (peak MNI coordinate: 1, −70, −40) mainly located in precuneus, IC 36 (−3, −78, 33) located in precuneus/cuneus, and IC 57 (33, 46, 25) located in the right middle frontal gyrus (Figure 6c). The age‐by‐diagnosis interactions of altered dynamic topological metrics were not significant (Table S4). Subgroup analysis showed no differences between drug‐naive and previously medication‐treated patients for any altered dynamic topological metrics (Table S5).

FIGURE 6

Group comparisons of dynamic topological metrics (a) and the correlations between dynamic topological metrics and illness duration and clinical symptom ratings in OCD patients (b and c). * indicated p‐value < .05

Correlation analysis showed that CV of AUC of nodal degree centrality and efficiency in IC 36 was negatively correlated with the duration of illness (r = −.468, uncorrected p = .018, 95% CI = − 0.669 ~ − 0.239, Power = 0.785; r = −.457, uncorrected p = .022, 95% CI = − 0.584 ~ − 0.232, Power = 0.760, respectively), and CV of AUC of nodal degree centrality in IC 57 was negatively correlated with obsessive subscale score (r = −.398, uncorrected p = .049, 95% CI = − 0.805 ~ −0.032, Power = 0.615) (Figure 6b). Other dynamic topological metrics showed no correlation with clinical scales or illness duration for patients with OCD.

DISCUSSION

In this study, we focused on temporal properties of whole brain dFC states in OCD. We found three discrete intrinsic connectivity configurations: a most frequent and moderately interconnected state (State 3), a least frequent and most strongly interconnected state (State 2), and a moderately frequent state with the lowest mean FC between networks (State 1). There were three major findings associated with OCD in our study: (a) Patients with OCD spent more fractional and dwell time in State 1, a state characterized by the highest modularity index with positive connections within DMN and SAN and negative coupling between DMN and SAN; (b) OCD patients demonstrated decreased temporal variability of nodal degree centrality and efficiency in the precuneus and middle frontal gyrus; and 3) Clinical relevance of these findings was reflected in a positive correlation between fractional time of state 1 and anxiety level, a negative correlation between nodal dynamic properties of precuneus and illness duration, and a negative correlation between dynamic nodal degree of the middle frontal gyrus and obsession symptom severity in patients with OCD. These findings indicate that the altered dynamic properties of dFC states and the variability of nodal topological metrics mainly involving DMN and SAN may be critical in the pathophysiology of OCD.

dFC states

In the three dFC states, we observed the lowest overall mean FC between brain networks in State 1 compared with the other two states. Compared with controls, patients with OCD spent more time in State 1, characterized by positive connections within DMN and SAN and negative coupling predominantly between DMN and SAN. These observations indicate an abnormal functional relationship between DMN and SAN as a core feature underpinning the neurophysiology of OCD. The DMN is activated when a person is “mind‐wandering” or with internally focused self‐reflection, and deactivated when externally focused on cognitive tasks and activities (Li et al., 2020; Raichle & Snyder, 2007). The SAN plays an important role in switching between processing modes, triggering the activation of executive regions, and disengagement from DMN regions when an external event requiring attention is detected (Menon, 2011; Sweeney et al., 1996).

Previous static rfMRI studies in patients with OCD found decreased FC between DMN and SAN and other task‐positive networks (Chen et al., 2018; Cyr et al., 2020), and a meta‐analysis of static FC studies found hypoconnectivity among FPN, SAN, and DMN (Gursel et al., 2018). The present findings suggest that a higher occurrence of a state with the lowest FC between networks might be the neural substrate for the decreased static FC between DMN and SAN in OCD. Reduced DMN‐SAN connectivity may contribute to decreased attention to the outside world (Posner et al., 2017) and the reduced behavioral flexibility and poor insight in patients with OCD (Fan et al., 2017). In our study, the ICs in State 1 tended to be subdivided into two functional modules with a high modularity index, indicating nodes that are maximally connected within the module but minimally connected with other modules. In State 1, DMN and SAN belonged to different modules, paralleling the low mean FC between DMN and SAN. Therefore, the increased fractional and mean dwell time, lowest FC between networks especially between DMN and SAN, and highest modularity in State 1 in patients indicates increased periods of excessive functional segregation in OCD and potentially a reduced ability to flexibly shift out of State 1.

In OCD patients, the fractional time of State 1 was positively correlated with anxiety ratings, which are a characteristic of OCD (Swartz et al., 2014; Weidt et al., 2016). Thus, from a clinical perspective, the inability to shift from State 1 might reflect a parallel inability to shift thinking about obsessive concerns that leads to increased anxiety (Hirschtritt, Bloch, & Mathews, 2017). Previous evidence also has shown abnormal functional interaction between DMN, central executive network (CEN), and SAN, and altered FC between SAN and CEN to be significantly correlated with anxiety level in OCD patients (Fan, Zhong, Gan, et al., 2017). Our findings on dFC states demonstrated that the functional network of patients with OCD has lost normal dynamic rhythms, and that this alteration with increased dwell time in a brain state with isolated network activity was related to clinical manifestations of the illness. These alterations may be considered as a potential neural marker for OCD and provide a novel target for intervention.

In contrast to findings in State 1, we did not find significant between‐group differences in the occurrence of States 2 and 3, or number of transitions among those states. These findings differ somewhat from findings of Liu et al. (2020)), who found an altered number of transitions. The patients in their study were all first‐episode and treatment‐naive, while half of patients in our study were previously medicated. In addition to patient cohort differences, methodological differences in the analysis stream may have contributed to these differences. Firstly, Liu et al. using a low model order GICA so that the brain rfMRI data were decomposed into 28 ICs, while we used a higher model order in our GICA that resulted in more refined delineations. Secondly, FC matrices of all sliding windows for all participants were clustered into four states in the Liu et al. study, altering state transition metrics.

Temporal variability of topological metrics

Global topological metrics of dFC were not robustly altered in patients with OCD, suggesting an intact temporal variability of the overall organization of functional brain networks. However, the temporal variability of nodal topological metrics, including degree centrality and efficiency of the precuneus in DMN and degree centrality of the right middle frontal gyrus in SAN, was decreased in patients compared with controls. In general, dFC in resting‐state analysis reflects the fluctuating exchange of information among brain regions at rest which occurs in an organized sequential way as shown in our analyses. Thus, decreased temporal variability of nodal properties of these regions in our study suggests reduced efficiency or flexibility in their communication with other brain regions. This pattern parallels the clinical expression of being locked into an obsessive preoccupation or behavioral pattern with high anxiety, with reduced flexibility in shifting mental sets and limited impact of external experiences to shift patients back to a more adaptive mode of function.

Interestingly, an electroencephalogram study of healthy participants found that regions in SMN, SAN, and DMN exhibit stronger mean FC strength and lower temporal variability than other brain areas (Shou et al., 2020). Regions that tend to maintain a high level of FC are likely to be less easily modulated based on contextual demands (Allen et al., 2014), perhaps contributing to their vulnerability in ways that increase the cognitive and behavioral rigidity that characterizes OCD. Previous rfMRI studies have found significant activation of the left precuneus (Menzies et al., 2008) as well as greater activity (Bohon, Weinbach, & Lock, 2020) and regional homogeneity in the middle frontal gyrus in patients with OCD (Yang et al., 2019), areas where we observed reduced dynamic flexibility. These observations are consistent with a relation between increased activity and connectivity with decreased temporal variability of these core regions in DMN and SAN, and that together these alterations may contribute to lower behavioral flexibility in OCD.

Furthermore, we found negative correlations between dynamic variability of regions in DMN and SAN and both illness duration and obsession symptom ratings in patients with OCD. The clinical relevance of the correlations can be interpreted in light of the function of these brain regions. The precuneus resides in the medial posterior parietal cortex and is a major constituent of the DMN involved in the generation of conscious self‐perception (Raichle et al., 2001). Patients with OCD have increased self‐referential thinking, often rethinking recent events and repetitively imagining future events in an attempt to increase control and decrease uncertainty (Pinto, Steinglass, Greene, Weber, & Simpson, 2014). The relation of this alteration with illness duration offers preliminary evidence for progressive brain abnormalities in this network in OCD.

The right middle frontal gyrus is part of the SAN, and it links the ventral and dorsal attention networks by acting as a “circuit‐breaker” responsible for the flexible modulation of endogenous and exogenous attention (Corbetta, Patel, & Shulman, 2008). Patients with OCD are less able to use available external information to counteract their negative thoughts and fears, which might be related to an imbalance or reduced flexibility in brain systems that subserve internally and externally focused attentional states (Stern & Taylor, 2014). Therefore, altered dynamics of the precuneus and right middle frontal gyrus may contribute to the clinical persistence and severity of symptoms in patients with OCD. Dynamic topological analysis revealed decreased temporal variability of brain regions located within DMN and SAN, and correlational analyses revealed that these were related to illness duration and symptom severity in OCD, suggesting these brain regions could be treated as biomarkers for illness progression and may have implications for precision medicine in patients with OCD.

Limitations

Certain limitations of our study need to be considered. First, the sample size was not large, perhaps providing insufficient statistical power to detect more modest alterations, and the generalization of the present results to other centers needs to be interpreted with caution. While replication is required, our validation analysis using different window sizes suggested that our findings were not a result of the specific window size choice used in our primary analyses. Second, 14 of 29 patients in the present study previously received psychopharmacological treatment. Although we found no differences between drug‐naive and previously treated patients in all dynamic metrics, the potential influences of prior medication on brain function cannot be completely excluded given limited statistical power to identify such differences. Longitudinal studies with controlled treatment are needed to establish medication effects on brain functional dynamic metrics investigated in the present study. Third, dFC analyses are particularly sensitive to head motion (Laumann et al., 2017). While we applied stringent inclusion criteria and found no group differences in head motion parameters, and motion parameters were regressed out in statistical analysis, head motion remains a potential source of artifact. Fourth, the scanning length of rfMRI in our study was not long. A longer scan with higher temporal resolution may provide a more refined assessment of dynamic FC features in OCD (Fu et al., 2018). Finally, although the sliding‐window method has been widely used, controversy remains regarding the optimal method to best capture dynamic fluctuations (Hindriks et al., 2016). Evaluating dFC with more novel methods and metrics might contribute to future advances in this field.

CONCLUSION

In summary, we found three dFC states in rfMRI scans, and patients with OCD exhibited a higher occurrence of State 1 than controls, characterized by the lowest mean FC between DMN and SAN and highest modularity. Additionally, we observed decreased temporal variability of topological organization, and potentially reduced functional flexibility, in brain regions within the DMN and SAN. This alteration, related to illness duration and symptom severity, may represent an important aspect of the pathophysiology of OCD. Our investigations of the temporal dynamics of brain FC states and topological networks facilitate understanding of the neural bases of the clinical presentation of OCD.

CONFLICT OF INTEREST

Dr. Sweeney has consulted to VeraSci. Other authors declare no biomedical financial interests or potential conflicts of interest.

AUTHOR CONTRIBUTIONS

Fei Li, Yanchun Yang, and Qiyong Gong conceptualized the project. Fei Li, Lekai Luo, and Qiyong Gong designed the study and drafted the manuscript. Fei Li, Lekai Luo, Qian Li, Wanfang You, Yuxia Wang, Bin Li, Wanjie Tang, and John A. Sweeney contributed to literature search, data collection, and analysis, as well as data interpretation. Fei Li, Lekai Luo, John A. Sweeney, and Qiyong Gong critically revised the manuscript. All authors approved the final version of the manuscript.

Table S1 Peak activation of 51 independent components in the present study.

Table S2. Mathematical definitions and interpretations of topological network measures used in the present study.

Table S3. The value of mean functional connectivity (FC) within and between networks in three states for all participants.

Table S4. Age by diagnosis (patients vs. controls) interactions using a general linear modal univariate analysis on altered state transition vectors and dynamic topological metrics.

Table S5. Comparison of state transition vectors and dynamic topological metrics between drug‐naive and previously medication‐treated patients with obsessive–compulsive disorder under the 22‐TR window size.

Table S6. The correspondence between pairs of centroids in three dynamic functional connectivity states under the two different window sizes (22TR and 30TR).

Click here for additional data file.

Figure S1 The spatial map of the selected 51 independent components (ICs) were sorted into eight functional networks including auditory network (AN), default mode network (DMN), subcortical network (SC), salience network (SAN), visual network (VN), executive control network (ECN), and sensorimotor network (SMN).

Click here for additional data file.

Figure S2 k‐means clustering analysis assessing the reoccurring dFC patterns. The optimal cluster number was determined to be three, in which Dunn's index was highest. While the five‐cluster solution had a similar Dunn value, the lower order three‐cluster was preferred for primary analyses because of its greater simplicity.

Click here for additional data file.

Figure S3 The scalar of small‐worldness sigma (σ) as a function of the sparsity threshold. Sigma value decreased with the increment of sparsity. In this study, the sparsity threshold was defined as 0.10 ≤ S ≤ 0.37, in which the averaged degree of each threshold over all nodes was larger than 2 × log(51) and the sigma value of each thresholded network was larger than 1.1 for all participants.

Click here for additional data file.

Figure S4 Results of the validation analysis under the window size of 30‐TR. (A), cluster centroids for each state in all participants. State 1 accounted for 15%, state 2 accounted for 35%, and state 3 accounted for 50%. (B and C), patients with obsessive–compulsive disorder (OCD) showed increased fractional time (p = 0.047, 95% confidence intervals [CI] = 0.023 ~ 0.307, Cohen's d = 0.495, Power = 0.643) and increased mean dwell time (p = 0.036, 95% CI = 1.177 ~ 31.014, Cohen's d = 0.511, Power = 0.666) than controls in state 2 (*, indicating P < 0.05), which was equivalent to State 1 in the primary 22‐TR analysis. (D), the overall number of transitions did not show a significant group difference. (E), for patients, the fractional time of State 2 was positively correlated with score of Hamilton Anxiety Rating Scale (HARS), and other state transition metrics were not correlated with symptom severity or illness duration.

Click here for additional data file.