Brain structural covariance networks in obsessive-compulsive disorder: a graph analysis from the ENIGMA Consortium.

Brain structural covariance networks reflect covariation in morphology of different brain areas and are thought to reflect common trajectories in brain development and maturation. Large-scale investigation of structural covariance networks in obsessive-compulsive disorder (OCD) may provide clues to the pathophysiology of this neurodevelopmental disorder. Using T1-weighted MRI scans acquired from 1616 individuals with OCD and 1463 healthy controls across 37 datasets participating in the ENIGMA-OCD Working Group, we calculated intra-individual brain structural covariance networks (using the bilaterally-averaged values of 33 cortical surface areas, 33 cortical thickness values, and six subcortical volumes), in which edge weights were proportional to the similarity between two brain morphological features in terms of deviation from healthy controls (i.e. z-score transformed). Global networks were characterized using measures of network segregation (clustering and modularity), network integration (global efficiency), and their balance (small-worldness), and their community membership was assessed. Hub profiling of regional networks was undertaken using measures of betweenness, closeness, and eigenvector centrality. Individually calculated network measures were integrated across the 37 datasets using a meta-analytical approach. These network measures were summated across the network density range of K = 0.10-0.25 per participant, and were integrated across the 37 datasets using a meta-analytical approach. Compared with healthy controls, at a global level, the structural covariance networks of OCD showed lower clustering (P < 0.0001), lower modularity (P < 0.0001), and lower small-worldness (P = 0.017). Detection of community membership emphasized lower network segregation in OCD compared to healthy controls. At the regional level, there were lower (rank-transformed) centrality values in OCD for volume of caudate nucleus and thalamus, and surface area of paracentral cortex, indicative of altered distribution of brain hubs. Centrality of cingulate and orbito-frontal as well as other brain areas was associated with OCD illness duration, suggesting greater involvement of these brain areas with illness chronicity. In summary, the findings of this study, the largest brain structural covariance study of OCD to date, point to a less segregated organization of structural covariance networks in OCD, and reorganization of brain hubs. The segregation findings suggest a possible signature of altered brain morphometry in OCD, while the hub findings point to OCD-related alterations in trajectories of brain development and maturation, particularly in cingulate and orbitofrontal regions.

Introduction

Three decades of neuroimaging research support the view that structural brain abnormalities in obsessive-compulsive disorder (OCD) do not merely involve alterations in discrete brain regions, but rather are best characterized in terms of altered networks of brain structures (Boedhoe ). More specifically, brain-based models of OCD have emphasized the role of the cortico-striato-thalamo-cortical loops and have also suggested the involvement of fronto-limbic, fronto-parietal and cerebellar regions (Menzies ; Milad and Rauch, 2012; de Wit ; Piras ; van den Heuvel ; Boedhoe , 2018; Fouche ). Most studies of brain networks in OCD have used resting state functional MRI (rs-fMRI) (Soriano-Mas and Harrison, 2017; Gürsel ), with alterations evident in intra-network connections of fronto-limbic and fronto-striatal networks (Anticevic ; Gottlich ; Posner ; Armstrong ; de Vries ; Takagi ). Furthermore, a meta-analysis of rs-fMRI studies comparing OCD to healthy controls found decreased intra-network connectivity of the fronto-parietal and salience networks, as well as reduced inter-network connectivity between the salience, fronto-parietal and default-mode networks (Gürsel ).

Brain structural covariance networks reflect intra-individual (Yun ; Seidlitz ) or inter-individual (Alexander-Bloch ; Kaczkurkin ; Wannan ) covariation in morphology of different brain areas, which may in turn point to common trajectories in brain development and maturation (Yun , 2016; Hunt ). Such networks may focus on a range of morphological features including regional brain volume (Spreng ), cortical thickness (Solé-Casals ), cortical surface area (Sharda ), and cortical white-grey contrast (Makowski ), as well as the paired or conjoint patterns between different brain regions (Seidlitz ; Hoagey ) Brain structural covariance has been estimated using Pearson’s correlation coefficient (Seidlitz ; Solé-Casals ; Wannan ), partial least squares (Hoagey ; Spreng ), non-negative matrix factorization (Kaczkurkin ), and inverse exponential of the difference between z-score transformed brain morphological values (Wee ; Yun , 2016), among others. Structural covariance networks are more similar to patterns of functional connectivity than the architecture of white matter connections, suggesting that areas that co-vary in morphological characteristics also belong to the same functional network (Zielinski ; Soriano-Mas ). Such networks are thought to be shaped by genetic and environmental influences from early childhood (Richmond ) and may continue to be reshaped during the lifespan (Alexander-Bloch ; Aboud ; Qi ) by a range of trophic influences (Ferrer ; Draganski ; Mechelli ).

Inter-individual brain structural covariance networks have been explored in a few studies of OCD and healthy controls. For example, Pujol found a negative association between relative volume reduction for OCD (compared to healthy controls) in the medial prefrontal-insulo-opercular cortical regions and relative volume enlargement of ventral striatum, suggesting that abnormal brain morphology in OCD might be distributed in coordinated fashion across diverse brain regions. In addition, a recent mega-analysis found higher covariance between volumes of left putamen and left frontal operculum, and higher covariance between volumes of right amygdala and ventromedial prefrontal cortex in OCD compared to healthy controls (Subira ). Further, local cortical gyrification (associated with cortical maturation)-based structural covariance network demonstrated lower covariance among mainly ventral brain regions in OCD compared to healthy controls (Reess ). However, few studies have explored intra-individual brain structural covariance networks in OCD; consequently our understanding of the factors that influence changes in global and regional network characteristics within individuals with OCD is limited.

The ENIGMA-OCD Working Group has collaborated on developing a large database of structural brain imaging in OCD and healthy controls, providing a unique opportunity to undertake such an exploration. Here we constructed intra-individual structural covariance networks from region of interest-based brain morphological features using 37 datasets worldwide (n = 1616 for OCD; n = 1463 for healthy controls), and investigated network topology using a graph theory approach. The current study aimed to capture the intra-individual distribution of brain morphological changes (Wee ; Yun , 2016) in OCD across 33 cortical surface areas, 33 cortical thickness values, and six subcortical volumes (Kremen ; Amlien ; Sussman ; Vijayakumar ; Krongold ; Schmaal ). Thus edge weights of the intra-individual structural covariance networks were estimated in proportion to the similarity between two brain morphological features in terms of deviation from healthy controls (i.e. z-score transformed). Networks were characterized at the global level using measures of network segregation (clustering coefficient and modularity), network integration (global efficiency), and their balance (small-worldness), as well as at the regional level using betweenness, closeness, and eigenvector centralities (Lancichinetti and Fortunato, 2009; Rubinov and Sporns, 2010; Cao ; Palaniyappan ; Vriend ). For preservation of the network edge weights-related information in the derived graph metrics, the global and regional graph metrics were summed across the network density range of K = 0.10–0.25 (Uehara ).

Previous neuroimaging studies of global network metrics have reported more (Zhang , 2014), less (Shin ; Armstrong ; Jung ; Reess ), or similar levels (Reess ) of segregated organization of white matter-based structural connectivity networks, resting state functional connectivity networks, or local gyrification index-based structural covariance networks in individuals with OCD, compared to healthy controls. These inconsistent findings raise the need for larger-scale meta-analysis. Therefore, the current study aimed to assess the level of global network segregation, as determined by the global clustering coefficient, using the largest dataset of structural covariance networks in OCD to date.

Materials and methods

Samples

This study included 37 datasets from 26 international research institutes participating in the OCD Working Group of the ENIGMA (Enhancing NeuroImaging and Genetics through Meta-Analysis) Consortium used in the meta-analytic between-group comparisons of OCD and healthy controls in terms of the subcortical volumes (Boedhoe ), cortical surface area and cortical thickness (Boedhoe ), in addition to the cortical and subcortical asymmetry (Kong ). Each dataset included demographic and neuroimaging data from OCD and healthy controls, as well as OCD clinical data (Table 1 and Supplementary material). The diagnosis of psychiatric disorders including OCD and other comorbid disorders (if any) was made using a structured or semi-structured interview; the Structured Clinical Interview for DSM-IV [SCID-I (First ); n = 23 datasets], the Mini-International Neuropsychiatric Interview [MINI (Sheehan ); n = 6 datasets], the Anxiety Disorder Interview Schedule [ADIS (Silverman ; Grisham ); n = 2 datasets], or the Schedule for Affective Disorders and Schizophrenia for School-Aged Children: Present and Lifetime Version [K-SADS-PL (Kaufman and Schweder, 2003); n = 7 datasets] (Table 1 and Supplementary material). Comorbid lifetime depressive disorder was present in 256 individuals with OCD, and comorbid lifetime anxiety disorder was present in 267 (Table 1 and Supplementary material). At the time of MRI acquisition, 721 individuals with OCD were on psychotropic medication and 881 were not. Age of illness onset of OCD was 18.8 ± 9.1 years, and illness duration was 10.8 ± 10.1 years (n = 1415). Severity of obsessive-compulsive symptoms was assessed with the Yale-Brown Obsessive-Compulsive Scale (Y-BOCS; for patients aged ≥18) or Children’s Y-BOCS (CY-BOCS; for patients aged <18); the mean score of 24.2 ± 6.8 (n = 1581) indicated a moderate to severe range of symptoms in the study population. All local institutional review boards permitted the use of extracted numerical measures for meta-analysis.

Table 1

Demographic and clinical information

Study	Study PI	Study site	MRI field strength, T	Total, n		Age, mean (SD)		Sex, male / female		Comorbid lifetime depression (OCD) n (%)	Comorbid lifetime anxiety (OCD) n (%)	Y-BOCS total (OCD), mean (SD)	Medicated OCD, n (%)	Illness duration (OCD), mean (SD)
				HC	OCD	HC	OCD	HC	OCD
1	Beucke	Berlin, GER	1.5	54	57	32 (11)	33 (11)	23 / 31	31 / 26	11 (19)	6 (11)	20 (7)	23 (40)	16.2 (11)
2	Cheng	Kunming, CHN	1.5	28	16	32 (8)	32 (12)	8 / 20	5 / 11	4 (25)	6 (38)	31 (7)	10 (63)	4.2 (5.2)
3	van den Heuvel	Amsterdam, NLD	1.5	35	37	31 (8)	35 (9)	12 / 23	11 / 26	11 (30)	6 (16)	23 (6)	0 (0)	20 (11.8)
4	Hoexter	San Paulo, BRA	1.5	9	38	28 (6)	31 (9)	5 / 4	17 / 21	20 (53)	24 (63)	28 (6)	8 (21)	17.4 (10.3)
5	Kwon	Seoul, KOR_01	1.5	103	45	24 (4)	25 (5)	57 / 46	34 / 11	0 (0)	0 (0)	20 (6)	11 (24)	7.3 (5.2)
6	Kwon	Seoul, KOR_02	1.5	45	34	25 (5)	29 (7)	29 / 16	19 / 15	1 (3)	0 (0)	24 (6)	0 (0)	9.9 (7.1)
7	Mataix-Cols	Stockholm, SWE	1.5	28	34	36 (11)	39 (11)	9 / 19	15 / 19	9 (26)	9 (26)	25 (8)	14 (41)	20.5 (14.9)
8	Menchon	Barcelona, ESP	1.5	55	95	32 (10)	35 (9)	22 / 33	47 / 48	15 (16)	19 (20)	25 (6)	91 (96)	13.9 (9.9)
9	Morgado	Braga, Portugal	1.5	51	58	28 (6)	27 (8)	19 / 32	27 / 31	–	–	26 (6)	58 (100)	–
10	Nakamae	Kyoto, JPN	1.5	48	81	30 (8)	32 (9)	25 / 23	37 / 44	18 (22)	8 (10)	25 (6)	39 (48)	6.7 (6.8)
11	Reddy	India	1.5	20	29	26 (6)	28 (7)	14 / 6	16 / 13	–	–	25 (9)	0 (0)	5.7 (5.3)
12	Benedetti	Milan, ITA	3	23	22	29 (11)	35 (11)	19 / 4	13 / 9	0 (0)	0 (0)	31 (6)	13 (59)	18.7 (12)
13	Cheng	Kunming, CHN	3	72	40	26 (4)	33 (11)	20 / 52	21 / 19	13 (33)	37 (93)	28 (6)	25 (63)	5.4 (5.8)
14	Denys	Amsterdam, NLD	3	15	14	38 (12)	34 (11)	6 / 9	1 / 13	4 (29)	1 (7)	27 (6)	9 (64)	14.9 (13.5)
15	van den Heuvel	Amsterdam, NLD	3	30	32	39 (11)	39 (11)	12 / 18	16 / 16	17 (53)	13 (41)	21 (6)	0 (0)	25.9 (12.9)
16	Koch	Munchen, GER	3	71	75	30 (9)	31 (10)	28 / 43	28 / 47	0 (0)	0 (0)	21 (6)	45 (60)	14.3 (10.6)
17	Kwon	Seoul, KOR	3	89	90	26 (7)	27 (7)	54 / 35	56 / 34	2 (2)	1 (1)	27 (7)	2 (2)	7.7 (6.7)
18	Nakamae	Kyoto, JPN	3	39	34	30 (7)	33 (10)	19 / 20	12 / 22	7 (21)	3 (9)	22 (7)	0 (0)	8.1 (6.1)
19	Nakao	Fukuoka, JPN	3	31	66	39 (13)	37 (10)	11 / 20	30 / 36	22 (33)	0 (0)	23 (6)	59 (89)	12.2 (9.3)
20	Nurmi	Los Angeles, USA	3	22	45	31 (12)	34 (11)	14 / 8	22 / 23	9 (20)	16 (36)	25 (4)	12 (27)	23 (10.8)
21	Reddy	India	3	139	201	26 (5)	30 (7)	86 / 53	107 / 94	31 (15)	15 (7)	26 (6)	82 (41)	7.3 (5.4)
22	Simpson	New York, USA	3	31	30	28 (8)	30 (8)	17 / 14	17 / 13	10 (33)	7 (23)	26 (4)	0 (0)	15.1 (8.7)
23	Spalletta	Rome, ITA	3	95	71	38 (11)	36 (11)	54 / 41	45 / 26	8 (11)	8 (11)	23 (9)	65 (92)	16.6 (11.4)
24	Stein	Cape Town, ZAF	3	25	21	31 (11)	31 (11)	10 / 15	11 / 10	0 (0)	0 (0)	23 (4)	9 (43)	17.9 (11.3)
25	Tolin	Conneticut, USA	3	32	27	48 (12)	32 (12)	7 / 25	18 / 9	11 (41)	12 (44)	23 (5)	21 (78)	–
26	Walitza	Zurich, CHE	3	15	13	33 (9)	31 (7)	4 / 11	7 / 6	6 (46)	7 (54)	18 (10)	6 (46)	12.8 (10)
27	Wang	Shanghai, CHN	3	35	47	26 (8)	30 (9)	18 / 17	23 / 24	0 (0)	0 (0)	25 (5)	0 (0)	6.5 (5.5)
28	Lazaro	Barcelona, ESP	1.5	29	29	15 (2)	14 (2)	14 / 15	18 / 11	0 (0)	5 (17)	22 (6)	15 (52)	2.1 (1.8)
29	Arnold	Ontario, CAN	3	11	34	12 (2)	13 (2)	6 / 5	20 / 14	7 (21)	10 (29)	21 (8)	21 (62)	4.2 (2.6)
30	Gruner	Conneticut, USA	3	17	10	14 (2)	15 (2)	8 / 9	9 / 1	2 (20)	6 (60)	27 (5)	6 (60)	–
31	Hoexter	San Paulo, BRA	3	26	27	12 (2)	13 (2)	15 / 11	16 / 11	6 (22)	20 (74)	27 (5)	12 (44)	5.5 (2.4)
32	Huyser	Amsterdam, NLD	3	20	20	14 (3)	14 (2)	8 / 12	6 / 14	7 (35)	9 (45)	26 (5)	0 (0)	3.1 (2.6)
33	Lazaro	Barcelona, ESP	3	43	53	15 (2)	15 (2)	23 / 20	30 / 23	3 (6)	14 (26)	19 (7)	42 (79)	2.5 (2.1)
34	Nurmi	Los Angeles, USA	3	36	53	13 (2)	13 (3)	18 / 18	29 / 24	1 (2)	2 (4)	24 (4)	7 (13)	–
35	Reddy	India	3	10	14	14 (3)	14 (2)	5 / 5	8 / 6	1 (7)	3 (21)	22 (7)	12 (86)	1.5 (0.9)
36	Soreni	Ontario, CAN	3	20	18	11 (3)	13 (2)	10 / 10	7 / 11	0 (0)	0 (0)	23 (4)	0 (0)	–
37	Walitza	Zurich, CHE	3	11	6	16 (2)	16 (1)	6 / 5	5 / 1	0 (0)	0 (0)	18 (10)	4 (67)	5 (2.4)

A more detailed version of this table is provided in the Supplementary material. A dash indicates data were not available.

BRA = Brazil; CAN = Canada; CHE = Switzerland; CHN = China; ESP = Spain; GER = Germany; HC = healthy control; ITA = Italy; KOR_01/02 = South Korea site 1/2; NLD = the Netherlands; PI = principal investigator; SWE = Sweden; Y-BOCS = Yale–Brown Obsessive Compulsive Scale; ZAF = South Africa.

Image acquisition and processing

Structural T1-weighted brain MRI scans were acquired and processed at each study site. For acquisition parameters of each site see Supplementary Table 1. All parcellations were performed with fully automated segmentation software FreeSurfer version 5.3. (Fischl, 2012), following standardized ENIGMA protocols (http://enigma.usc.edu/protocols/imaging-protocols/). To ensure quality control, we visually inspected the segmentations of 68 (34 left and 34 right) cortical grey matter regions and seven subcortical regions based on the Desikan-Killiany atlas (Desikan ) and statistically evaluated the data for outliers (Boedhoe , 2018). We excluded the volume values of bilateral entorhinal cortices and the nucleus accumbens because of segmentation issues (as calculation of intra-individual brain structural covariance networks requires every region of interest to be adequately measured in each participant; inclusion of regions of interest with relatively poorer quality segmentations would effectively decrease sample size).

Intra-individual cortical-subcortical structural covariance networks

As illustrated at ‘step 1’ in Fig. 1, bilaterally-averaged values (where brain regions were poorly segmented in one hemisphere, the value from the contralateral hemisphere was used as a proxy) of 33 cortical surface area regions of interest, 33 cortical thickness regions of interest, and six subcortical volume regions of interest, were corrected for age, sex, and individual brain size (Vuoksimaa ) per dataset (n = 37). The resulting residuals were then z-score transformed using mean and SD values of each region of interest calculated from healthy controls (to derive the degrees of brain morphological variations per region of interest relative to the ‘average healthy controls’ values). Finally, a measure of joint variation (which is not the same as the classical statistical definition of covariance) between the 72 morphometric features (33 cortical surface area values, 33 cortical thickness values, and six subcortical values) represented the edge-weights (distributed between 0 and 1) of the network and was calculated using the following formula (Yun , 2016):

Figure 1

Schematic description of the study procedures: construction of intra-individual brain structural covariance networks. HC = healthy controls; L = left; M = mean; R = right; ROI = region of interest; SD = standard deviation.

[Intra-individual brain structural covariance (joint variation) between the ith (for i = 1 to 72) and j-th (for j = 1 to 72) regions of interest in the k-th (for k = 1 to ‘total number of participants per dataset’) participant] = 1/exp{[(z-transformed value of i-th region of interest in k-th participant) − (z-transformed value of j-th region of interest in k-th participant)]2} (1)

Graph theory approach: single subject level

Global network characteristics

Intra-individual structural covariance networks were thresholded (using ‘threshold_proportional.m’ function in network density range of K = 0.05–0.30; with interval of 0.01) and binarized (using the ‘weight-conversion.m’ function; e.g. when we applied a density threshold of K = 0.10, the edge weights in the network were sorted into numerical order and a cut-off was applied to retain only the strongest 10% of edges with edge weights converted to ‘1’ and edges weights for other remaining edges becomes ‘0’ (Fig 2, steps 2A and 3A). From these thresholded and binarized networks, four global metrics were determined: (i) global clustering (a tendency for brain regions to segregate into locally interconnected triplets of neighbouring nodes); (ii) global modularity (a measure of the segregation of the network into communities where nodes are more strongly connected with each other than nodes outside the community because of similar morphological characteristics; this measure is operationalized as the most frequently occurring value over 500 runs of estimation using ‘modularity_und.m’) (Newman, 2006; Reichardt and Bornholdt, 2006); (iii) global efficiency (how well on average each node is connected to all others based on the minimum number of steps nodes are separated from each other); and (iv) small-worldness (a measure of balance between the degree of segregation versus integration in brain network) using the Brain Connectivity Toolbox (Rubinov and Sporns, 2010) in MATLAB R2017a (Weinberg ; Das ; Zaremba ).

Figure 2

Schematic description of the study procedures. (A) Calculation of graph theory metrics from the intra-individual brain structural covariance networks at single-subject level and (B) meta-analytic integration of graph theory metrics for 37 datasets. HC = healthy controls; ROI = region of interest.

Among the diverse network density levels of K = 0.05–0.30 (with density interval of 0.01), only in the narrower network density levels of K = 0.10–0.25, three criteria of (Uehara ) (i) network connectedness (> 80% of nodes remain connected to other nodes within the network); (ii) modular organization (modularity > 0.3); and (iii) small-world organization (small-worldness > 1) were satisfied for >95% of the intra-individual structural covariance networks comprising each dataset (n = 37). Therefore, these network density levels of K = 0.10–0.25 (density interval = 0.01) were selected for the between-group comparison of global network characteristics, community membership detection, and hub profiling using the regional network characteristics (Fig. 2, step 3A). Estimation of the global network characteristics was done using Brain Connectivity Toolbox (https://www.nitrc.org/projects/bct/) in MATLAB R2017a.

Detection of community membership

In addition, we assessed community membership (Fortunato, 2010) for each structural covariance (joint variation) network. For thresholded (K = 0.10–0.25) and binarized intra-individual structural covariance (joint variation) networks, detection of communities [i.e. densely connected subgroups of nodes in a network (Power )] was conducted using the InfoMap algorithm (Rosvall and Bergstrom, 2007; Fortunato, 2010; Power ; Kawamoto and Rosvall, 2015). First, a participant-level co-classification matrix (Dwyer ) that represented the fraction of network density level, in which each pair of nodes was clustered into the same community according to the InfoMap algorithm (Rosvall and Bergstrom, 2007; Kawamoto and Rosvall, 2015), was generated. Second, the InfoMap algorithm was applied to this co-classification matrix to generate a participant-level consensus of community membership (Fornito ). All procedures other than the InfoMap-based community estimation were done using MATLAB R2017a software (https://kr.mathworks.com).

Hub profiling and regional network characteristics

Principal brain regions that could be essential indicators of brain morphological changes within the network were assessed using hub profiling, which provided three local network measures: (i) betweenness centrality (the frequency of a node being located in the shortest path for each pair of two other nodes in a network); (ii) closeness centrality (the ease with which one node can reach all other nodes within a network); and (iii) eigenvector centrality (a self-referential measure of centrality that reflects the presence of connectedness of one node to other nodes with high eigenvector centrality) (Rubinov and Sporns, 2010) (Fig 2, step 2B). As distribution of these local network measures does not follow normal distribution in a scale-free network, prior to the between-group comparison and meta-analysis, these regional centrality metrics were rank-transformed using the ‘tiedrank.m’ function of MATLAB R2017a and were averaged in the network density range of K = 0.10–0.25 to be re-ranked at participant-level; participant-level hubs were selected as top-10 ranked nodes in two or three centralities. All of the procedures described above were conducted using the Brain Connectivity Toolbox (Rubinov and Sporns, 2010) and MATLAB R2017a software (https://kr.mathworks.com).

Meta-analysis of graph metrics

Meta-analysis of between-group differences in global network characteristics across the whole dataset (n = 37; Fig. 2, step 3A) was performed using a random-effects meta-analytic model (Hedges and Vevea, 1998; Kambeitz ) incorporating the bias-corrected standardized mean difference (SMD = Hedges’ g) between OCD and healthy controls for each of the four global network characteristics (summated over the network density range of K = 0.10–0.25) that satisfied network connectedness, modular organization, and small-world organization; see ‘Graph theory approach: single subject level’ section). Summary effect sizes were calculated with restricted maximum-likelihood estimator (REML) (Raudenbush, 2009; Viechtbauer, 2010). Estimates for heterogeneity were assessed with the I2 value (Raudenbush, 2009). For all analyses, a significance level of P < 0.01 was used, i.e. P < 0.05/5 number of global network characteristics (= 4) plus local network-related measure of the Dice coefficient (= 1; see section below) (Kambeitz ). All statistical analyses were conducted using the R package ‘metafor’ version 2.0.0 (Viechtbauer, 2010).

Community membership

First, summation of network-transformed community profiles for each individual provided dataset-level co-classification matrices (in which higher edge weights indicated that two nodes were clustered in the same community across a large proportion of participants in dataset) for OCD and for healthy controls (Fornito ). Second, consensus of community membership at dataset level (for OCD and healthy controls separately) was estimated by applying the InfoMap algorithm to the weighted and thresholded (at density level of K = 0.10) version of the dataset-level co-classification matrices. Third, dataset-level consensus community profiles of OCD and healthy controls were binarized, multiplied by the square root of participants number per dataset, and summed to generate the meta-analytic co-classification matrices of OCD or healthy controls (n = 37). Finally, a weighted and thresholded (at density level of K = 0.10) version of these meta-analytic co-classification matrices underwent InfoMap-based community detection, to determine the meta-analytic consensus community profile for brain structural covariance networks of OCD and healthy controls. All procedures other than the InfoMap-based community estimation were performed using MATLAB R2017a.

Hub profiling and regional network characteristics

In the current study, hub profiling was done to find the principal brain regions that could be essential indicators of intra-individual distribution of brain morphological changes (= deviation from healthy controls) based on the three local metrics of betweenness, closeness, and eigenvector centralities (Fig. 2). Three rank-transformed centralities (betweenness, closeness, and eigenvector) were rank-transformed at the participant level, and were averaged in the network density range of K = 0.10–0.25. The top 10 ranked nodes (i.e. 10 nodes illustrated in Fig. 4 and Supplementary Fig. 1) for two or three centralities as calculated from the summation of participant-level centrality values within each dataset (n = 37) were classified as dataset-level hubs for OCD or healthy controls. Finally, meta-analytic hub scores for all network nodes (= 33 cortical surface area values + 33 cortical thickness values + six subcortical volumes) were calculated by summing the values of [(presence (= 1) or absence (= 0) of network nodes in the hub profile of each dataset) × (square root of participants number per dataset)] across the whole dataset (n = 37) for OCD and healthy controls separately; top-10 ranked nodes for this meta-analytic hub score were defined as meta-analytic hubs for OCD or healthy controls, respectively.

Figure 4

Meta-analysis of community membership and hubs. (A) Healthy bontrols (HC); and (B) OCD. Spheres represent nodes [= bilaterally-averaged values of 33 cortical surface areas (CSAs), 33 cortical thickness (CT), and six subcortical volumes (vol)] comprising the intra-individual structural covariance network. Larger spheres represent hubs, and differential colours were used to denote the spheres (or network nodes) segregated as different modules.

Between-group comparison of rank-transformed centrality values at the dataset level (n = 37) was performed using the Wilcoxon rank sum test. Nodes that showed statistically significant differences between OCD and healthy controls (P < 0.05) were recoded into MNI coordinates using brainGraph (https://cran.r-project.org/web/packages/brainGraph), and underwent coordinate-based meta-analysis, i.e. activation likelihood estimation (ALE), using gingerALE version 2.3.6. (Eickhoff ). In this ALE-based meta-analysis, nodes that showed significant effect sizes [cluster-level corrected threshold of P < 0.05 (family-wise error, FWE); cluster-forming threshold at voxel level of P < 0.001] for between-group differences in two or three centralities were considered valid (Fig. 2, step 3B).

Lastly, to explore the difference in hubs in terms of their topographical location between OCD and healthy controls, we also calculated the Dice similarity coefficient (Dice, 1945), a measure of the degree of overlap between each participant-level hub profile versus the reference (= hub profile of healthy controls per dataset). For meta-analysis, the bias-corrected SMD (Hedges’ g) of Dice similarity coefficient (i) between the healthy controls and OCD (37 dataset) as well as (ii) between unmedicated OCD and medicated OCD (12 dataset in which >10 participants existed for all of the two subgroups) were calculated and entered into a random-effects meta-analytic model (Schmidt ; Kambeitz ). Summary effect sizes were calculated with REML (Raudenbush, 2009; Viechtbauer, 2010), and estimates for the amount of heterogeneity were assessed by way of the I2 value (= the percentage of total variability across dataset that is due to heterogeneity than by chance) (Higgins ). For all analyses, a significance level of P < 0.05 (two-tailed) was used (Kambeitz ) and all statistical analyses were conducted using the R package ‘metafor’ version 2.0.0 (Viechtbauer, 2010).

Influence of comorbid lifetime depressive or anxiety disorders in patients with OCD

Thirty-five (of 37) datasets provided information about comorbid lifetime depressive and anxiety disorders in OCD individuals; meta-analysis of global network characteristics and Dice coefficients was conducted to assess between-group differences in (i) OCD with and without comorbid lifetime depressive disorder (n = 10 datasets, in which n > 10 for both OCD subgroups); and (ii) OCD with and without comorbid lifetime anxiety disorders (n = 7 datasets, in which n > 10 for both OCD subgroups).

Influence of medication

Twenty-seven (of 37) datasets provided information about medication status (= presence or absence of psychotropic medication prescribed at the time of MRI data acquisition) of OCD individuals; meta-analytic integration for the between-group comparison of regional network characteristics (= centralities) between medicated OCD versus unmedicated OCD was undertaken for these datasets. Furthermore, meta-analytic integration of between-group differences for global network metrics and Dice coefficients were conducted using results retrieved from 12 datasets (in which n > 10 for both medicated and unmedicated subjects).

Influence of OCD illness duration

Fisher’s z-transformed correlation coefficients between the OCD illness duration and four global network metrics were calculated per dataset (n = 32 datasets). Each of these correlation coefficients per dataset and per global network characteristics were meta-analytically integrated using the same pipeline as for the global network characteristics. Likewise, Spearman correlation coefficients between the OCD illness duration and rank-transformed (betweenness, closeness, or eigenvector) centrality measures were also calculated per dataset. Meta-analysis of the dataset-level nodes that showed significant correlation with OCD illness duration (P < 0.05) was performed using using gingerALE version 2.3.6 [P < 0.05 (cluster-level FWE)] (Eickhoff ).

Data availability

De-identified data are available from the corresponding author upon reasonable request.

Results

Patients with OCD versus healthy controls

Demographic and clinical characteristics

A total of 37 datasets worldwide (n = 1616 for OCD; n = 1463 for healthy controls) were included in this study. Demographic and clinical characteristics for each dataset are described in Table 1 and Supplementary material. Between-group (OCD versus healthy controls) statistical tests for age (using the independent t-test) and sex ratio (using the chi-squared test) did not show statistically significant differences between OCD and healthy controls (P > 0.05) for 31 (83.8%) and 34 datasets (91.9%), respectively. On the other hand, years of education (information available for 27 datasets) were fewer in OCD compared to healthy controls (P < 0.05) in 10 (27.0%) datasets.

Meta-analysis of global network characteristics for the intra-individual brain structural covariance networks (Table 2 and Fig. 3A–D) showed lowered global clustering and modularity in OCD compared to healthy controls (all P’s < 0.01). Global efficiency and small-worldness did not differ significantly between OCD and healthy controls (all P’s > 0.01). When the sample was divided into two groups (adults and adolescents), and analyses run in each, these findings continued to hold true (Table 2). Additional meta-analyses using years of education as a moderator did not show any significant influence of this variable (all P’s > 0.05) on either the global network metrics of global clustering (Q = 1.456, df = 2, P = 0.483), modularity (Q = 0.819, df = 2, P = 0.664), global efficiency (Q = 0.673, df = 2, P = 0.714), and small-worldness (Q = 0.139, df = 2, P = 0.933), or on the Dice similarity coefficient (Q = 1.447, df = 2, P = 0.485).

Table 2

Meta-analysis of global network characteristics and Dice similarity coefficients

	logSMD	k	z	P-value	95% CI	I 2 (%)	Q	P
OCD versus HC
Global clustering coefficient (total)	0.77	37	−6.94	<0.001	0.72 to 0.83	0.01	44.8	0.149
Adults (≥18 years)	0.79	27	−5.89	<0.001	0.73 to 0.85	<0.001	26.8	0.418
Adolescents (<18 years)	0.66	10	−3.16	0.002	0.50 to 0.85	45.7	16.5	0.058
Modularity (total)	0.82	37	−5.21	<0.001	0.77 to 0.89	0.01	43.1	0.194
Adults (≥18 years)	0.84	27	−4.28	<0.001	0.78 to 0.91	0.01	22.3	0.670
Adolescents (<18 years)	0.68	10	−2.63	0.009	0.51 to 0.91	54.0	19.1	0.025
Small-worldness (total)	0.92	37	−2.39	0.017	0.85 to 0.98	0.001	26.2	0.886
Adults (≥18 years)	0.93	27	−1.82	0.069	0.86 to 1.01	<0.001	18.1	0.872
Adolescents (<18 years)	0.84	10	−1.82	0.068	0.70 to 1.01	<0.001	7.2	0.621
Global efficiency (total)	0.98	37	−0.54	0.586	0.91 to 1.05	0.02	38.5	0.358
Adults (≥18 years)	0.97	27	−0.68	0.494	0.89 to 1.06	10.7	32.6	0.174
Adolescents (<18 years)	1.05	10	0.50	0.621	0.87 to 1.26	<0.001	5.3	0.809
Dice similarity coefficient (total)	0.48	37	−14.36	<0.001	0.43 to 0.53	39.35	58.3	0.011
Adults (≥18 years)	0.49	27	−11.97	<0.001	0.44 to 0.55	45.7	49.6	0.004
Adolescents (<18 years)	0.41	10	−9.18	<0.001	0.34 to 0.50	<0.001	2.9	0.969
OCD patients with versus without lifetime comorbid depressive disorder
Global clustering coefficient	0.89	10	−1.13	0.257	0.73 to 1.09	11.6	10.1	0.344
Modularity	0.90	10	−0.91	0.365	0.72 to 1.13	26.5	11.9	0.217
Small-worldness	0.96	10	−0.44	0.659	0.80 to 1.15	<0.001	6.6	0.678
Global efficiency	1.00	10	−0.04	0.966	0.83 to 1.20	<0.001	5.8	0.764
Dice similarity coefficient	1.04	10	0.32	0.751	0.84 to 1.28	21.4	13.2	0.155
OCD patients with versus without lifetime comorbid anxiety disorder
Global clustering coefficient	0.99	7	−0.09	0.929	0.79 to 1.25	<0.001	5.1	0.531
Modularity	0.96	7	−0.39	0.695	0.76 to 1.20	<0.001	1.8	0.934
Small-worldness	1.00	7	0.01	0.993	0.79 to 1.27	3.1	6.6	0.357
Global efficiency	1.03	7	0.24	0.814	0.82 to 1.29	<0.001	6.2	0.403
Dice similarity coefficient	1.15	7	0.96	0.338	0.87 to 1.52	31.2	8.3	0.215
Medicated OCD versus unmedicated OCD
Global clustering coefficient	0.95	12	−0.63	0.531	0.82 to 1.11	0.00	11.03	0.441
Modularity	0.94	12	−0.83	0.408	0.8 to 1.09	0.00	8.73	0.647
Small-worldness	0.99	12	−0.08	0.934	0.83 to 1.18	15.25	12.97	0.295
Global efficiency	0.85	12	−1.66	0.097	0.7 to 1.03	28.05	13.32	0.273
Dice similarity coefficient	1.06	12	0.72	0.474	0.91 to 1.24	2.16	6.90	0.807
	Correlation coefficient	k	z	P-value	95% CI	I2	Q	P
Illness duration in OCD
Global clustering coefficient	−0.03	32	−0.85	0.393	−0.11 to 0.04	40.13	53.37	0.008
Modularity	−0.05	32	−1.32	0.188	−0.12 to 0.02	34.70	48.57	0.023
Small-worldness	−0.02	32	−0.67	0.584	−0.10 to 0.05	34.00	46.28	0.038
Global efficiency	−0.02	32	−0.59	0.558	−0.07 to 0.04	0.00	20.64	0.921

CI = 95% confidence interval; I2 = total heterogeneity/total variability; k = number of studies included in given meta-analysis; log SMD = log-transformed standardized mean difference; P = P-value of heterogeneity test; P-value = P-value of random effect model (REML); Q = heterogeneity score; z = z-score.

Figure 3

Forest plots of the meta-analysis of global graph metrics comparying the OCD and healthy control groups. (A) Global clustering, (B) small-worldness, (C) modularity, (D) global efficiency, and (E) dice similarity coefficient. HC = healthy controls; ROI = region of interest.

Community membership analysis detected that the healthy controls network had six modules (or subgroups within the network), while the OCD network had three modules, indicative of less global network segregation. The six community modules of the healthy controls network (Fig. 4 and Supplementary Fig. 1) were module 1 [the principal (31 nodes); including six hubs of cortical surface area for medial orbitofrontal, caudal middle frontal, and parahippocampal cortices, as well as cortical thickness for posterior cingulate, pars triangularis, and insula], module 2 [cingulate-parietal-inferior frontal (13 nodes)], module 3 [subcortical (six nodes); including two hubs named pallidal and caudate volumes], module 4 [frontal pole-occipital (six nodes); including cortical thickness for cuneus as hub], module 5 [paracentral-temporal (six nodes); including a hub of paracentral cortical thickness], and module 6 [insula-perisylvian (five nodes)]. As smaller communities with less than four nodes (<5% of total nodes) were excluded, six nodes comprising module 2 for healthy controls [cortical surface area of caudal-rostral anterior cingulate and lateral orbitofrontal cortices, in addition to cortical thickness of paracentral, superior parietal, and supramarginal cortices] were not classified in these communities.

In contrast, community membership of individualized structural covariance networks for OCD (Fig. 4 and Supplementary Fig. 1) showed just three modules: module 1 [in which eight OCD hubs for cortical surface area of superior temporal sulcus (module 1 in healthy controls), posterior cingulate (module 2 in healthy controls), rostral middle frontal-insular-superior temporal (module 6 in healthy controls), and pericalcarine cortices, as well as cortical thickness of caudal anterior cingulate-frontal pole (module 1 in healthy controls) included], module 2 [comprising cortical thickness of inferior parietal lobule-precuneus (module 2 in healthy controls) in addition to cuneus-lingual-pericalcarine gyri (module 4 in healthy controls)], and module 3 (includes a hub named hippocampal volume).

Regional network characteristics

Of the 10 hubs for the OCD network (Fig. 4 and Supplementary Fig. 1), only one node, i.e. cortical thickness of postcentral cortex [member of the paracentral-temporal module in healthy controls; fifth community (red square) in Supplementary Fig. 1], was found among the 10 healthy controls hubs. Meta-analysis of Dice similarity coefficients showed lower Dice similarity coefficient in OCD compared to healthy controls (Table 2 and Fig. 3E), indicating that the nodes classified as hubs differed between OCD and healthy controls. In terms of the centralities, compared to healthy controls, rank-transformed centrality of caudate nucleus volume was lower in OCD (healthy controls hub; Fig. 5A and Supplementary Fig. 3).

Figure 5

Meta-analysis of regional network characteristics (= rank-transformed betweenness, closeness, and eigenvector centralities). (A) Comparing OCD and healthy controls (HC); (B) comparing medicated OCD with unmedicated OCD; and (C) estimating the degrees of relationship with illness duration for OCD. CSA = cortical surface areas; CT = cortical thickness.

No significant differences in global network characteristics or Dice similarity coefficients were found between OCD with comorbid lifetime depression versus OCD without lifetime depression, nor between OCD with comorbid lifetime anxiety disorder versus OCD without lifetime anxiety disorder (Table 2).

Influence of medication for OCD

No significant differences in global network characteristics or Dice similarity coefficients were found between medicated and unmedicated OCD (Table 2). The structural covariance networks of healthy controls, medicated OCD, and unmedicated OCD demonstrated five, three, and two modules (or subgroups within the network), respectively (Supplementary Fig. 2).

OCD illness duration did not show significant correlations with global network characteristics (Table 2). However, OCD illness duration showed significant positive relationships with centrality (Fig. 5C and Supplementary Fig. 5) of cortical thickness for caudal anterior cingulate (OCD hub), cortical surface area for posterior cingulate (OCD hub), and cortical surface area of lateral orbitofrontal cortex (non-hub). Furthermore, OCD illness duration showed significant negative correlations with centrality of the cortical surface area for parahippocampal cortex (healthy control hub), cortical thickness for the frontal pole, cortical surface area for superior temporal and pericalcarine cortices (OCD hubs), cortical thickness for inferior parietal lobule, and cortical surface areas for inferior temporal and cingulate isthmus cortices (non-hubs).

Discussion

The current meta-analysis of 37 datasets from 26 sites worldwide is the largest investigation of structural covariance networks in OCD to date. Two main findings emerged. First, we observed lower clustering, modularity, and small-worldness of OCD brain structural covariance networks compared with healthy controls, with community membership analysis confirming a less segregated organization of the global structural covariance network of OCD patients. Second, hub profiling demonstrated reduced centralities in subcortical volumes of caudate nucleus and thalamus as well as cortical surface area of paracentral cortex in OCD. Alterations in hub organization were associated with both medication status and illness duration. These novel findings are important; the first suggests a possible signature of altered brain morphometry in OCD compared to healthy controls, and the second provides evidence for OCD-related alterations in trajectories of brain development and maturation.

Lower clustering, modularity and small-worldness, but normal global efficiency, are indicative of lower global segregation, but spared global integration of OCD networks. In particular, lower modularity might be related to over-connectedness of certain nodes and diminished ability of the network to adapt flexibly (Guye ). This finding is consistent with previous observations of abnormal brain network segregation in functional networks in OCD (Zhang ). Small-worldness relates to an optimal network organization that combines regional specialization and efficient global (Watts and Strogatz, 1998; Latora and Marchiori, 2001; Lefort-Besnard ). Thus, despite intact global efficiency, decreased levels of small-worldness and modularity in OCD point to a disrupted hierarchical network architecture.

Global network findings were not impacted by medication status or illness duration. This contrasts with previous research, which although based on functional MRI data, suggested that abnormal global network characteristics may depend on psychotropic treatment (Shin ). Although it is theoretically possible that the effects of psychotropic medication on OCD brain morphology differ in the acute versus chronic stage of pharmacotherapy so that there the net result over time is one of no change, there is little evidence to support this idea. In our view, a more plausible conclusion is that the lower global network segregation found here may represent a possible signature of altered brain morphometry in OCD. Further research is needed to confirm this.

The study also found reduced centralities of caudate nucleus and thalamic volumes in OCD compared to healthy controls. This is in line with our previous multicentre mega-analysis, which showed increased thalamic volume in OCD compared to healthy controls, even though only in the paediatric patients (Boedhoe ). Likewise, caudate nucleus and thalamus showed marked expansion in OCD and in their unaffected siblings compared to healthy controls, suggesting genetic contributions to altered brain morphology (Shaw ). Of note, meta-analytic integration of task-related functional MRI studies demonstrated OCD-specific differences in functional activation of the caudate nucleus. Similarly, nodal efficiency of the caudate nucleus was reduced in OCD in a white matter-based structural connectivity network (Zhong ), in line with a resting state functional connectivity profile that showed increased intra-subcortical modular connections for caudate nucleus and thalamus in OCD (Vaghi ).

Our data emphasize that alterations in hubs in OCD are associated with illness duration. This is consistent with previous work suggesting brain-related changes during the development of OCD (van den Heuvel ). In particular, we found that centralities of brain regions including the cortical thickness of caudal anterior cingulate as well as the cortical surface areas for posterior cingulate and lateral orbitofrontal cortices, were associated with longer illness duration in OCD. As an interface between sensorimotor, limbic and executive networks, the caudal anterior cingulate plays a major role in attentional control (Margulies ) and self-referential sensorimotor processing (Jung ; Mao ), the posterior cingulate cortex and connected default mode network supports internally-directed cognition, participates in the control of arousal state, and interacts with other brain regions for attentional modulation and conscious awareness (Leech and Sharp, 2014). Orbitofrontal regions have also previously been emphasized in OCD. The hub findings reported here point to OCD-related alterations in trajectories of brain development and maturation, particularly in cingulate and orbitofrontal regions. However, these hypotheses will require confirmation in longitudinal studies.

This study has some limitations that deserve emphasis. First, the current study analysed datasets that were extracted from brain MRI data collected from 26 international research institutions using diverse acquisition parameters (Boedhoe , 2018), which may have introduced systematic biases. Nevertheless, our meta-analytical approach took into account differential site effects. Second, although all brain segmentation results underwent quality check procedures prior to extraction of numerical values, we were unable to implement motion correction of structural images, and it is theoretically possible that estimates of group differences are inflated by uncorrected motion. Nevertheless, there is no reason to suspect increased motion in either group. Third, in the calculation of intra-individual structural covariance networks, the current study applied the bilaterally-averaged values of 33 cortical surface area regions of interest, 33 cortical thickness regions of interest, and six subcortical volume regions of interest and therefore did not explore the homologous connectivity between the brain regions. However, we would like to emphasize that patterns of brain cortical-subcortical morphological asymmetry in adult OCD are not significantly different from healthy controls (Kong ). Fourth, the current study did not explore the possible effect of other clinical features such as the severity of depressive or anxiety symptoms, and IQ score, on the brain morphological features, because of the lack of sufficient information. Fifth, this was a cross-sectional study and any conclusions regarding developmental trajectories are necessarily tentative.

Taken together, this study showed that the structural covariance networks of individuals with OCD are less segregated and show a reorganization of brain hubs, compared to healthy controls. These findings support the hypothesis that OCD brain abnormalities are best described at the network level and involve alterations in the hierarchical structure of the brain. The segregation findings here are important insofar as they suggest a possible signature of altered brain morphometry OCD, while the hub findings are useful in emphasizing the importance of OCD-related alterations in trajectories of brain development and maturation, particularly in cingulate and orbitofrontal regions.

Funding

This research was funded by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2017R1D1A1B03028464) and the Basic Research Laboratory Program through the National Research Foundation of Korea (NRF) (Grant no. 2018R1A4A1025891).

ENIGMA is supported, in part, by an NIH grant (U54 EB020403) for big data analytics.

Competing interests

P.D.A (Alberta Innovates Translational Health Chair in Child and Youth Mental Health). D.P.H. is now an employee of Genentech, Inc. working on projects unrelated to this publication. D.M-C. receives royalties for contributing articles to UpToDate, Wolters Kluwer Health and fees from Elsevier for editorial tasks (all unrelated to the submitted work). H.B.S. (Biohaven research support for clinical trial; Royalties from UpToDate, Inc and Cambridge University Press). N.S. (Lundbeck-IIT). P.M.T. has received a research grant from Biogen, Inc. unrelated to the topic of this paper. No further conflict of interest is reported.

Click here for additional data file.