Graph theory analysis of cortical thickness networks in adolescents with d‐transposition of the great arteries

Abstract Objective Adolescents with d‐transposition of the great arteries (d‐TGA) who had the arterial switch operation in infancy have been found to have structural brain differences compared to healthy controls. We used cortical thickness measurements obtained from structural brain MRI to determine group differences in global brain organization using a graph theoretical approach. Methods Ninety‐two d‐TGA subjects and 49 controls were scanned using one of two identical 1.5‐Tesla MRI systems. Mean cortical thickness was obtained from 34 regions per hemisphere using Freesurfer. A linear model was used for each brain region to adjust for subject age, sex, and scanning location. Structural connectivity for each group was inferred based on the presence of high inter‐regional correlations of the linear model residuals, and binary connectivity matrices were created by thresholding over a range of correlation values for each group. Graph theory analysis was performed using packages in R. Permutation tests were performed to determine significance of between‐group differences in global network measures. Results Within‐group connectivity patterns were qualitatively different between groups. At lower network densities, controls had significantly more long‐range connections. The location and number of hub regions differed between groups: controls had a greater number of hubs at most network densities. The control network had a significant rightward asymmetry compared to the d‐TGA group at all network densities. Conclusions Using graph theory analysis of cortical thickness correlations, we found differences in brain structural network organization among d‐TGA adolescents compared to controls. These may be related to the white matter and gray matter differences previously found in this cohort, and in turn may be related to the cognitive deficits this cohort presents.


| INTRODUCTION
Congenital heart disease (CHD) is the most commonly occurring congenital anomaly (Tennant, Pearce, Bythell, & Rankin, 2010). Due to improvements in medical and surgical care, a steadily increasing proportion of those born with CHD are surviving into adolescence and adulthood. Research has increasingly focused on the behavioral and neuropsychological deficits present throughout development in this cohort (Bang et al., 2013;Bellinger et al., 2011;Cassidy, White, DeMaso, Newburger, & Bellinger, 2015;Heinrichs et al., 2014;Marino et al., 2012). Recently, evidence of differences in brain structure of adolescents with CHD has accumulated, and some of these differences have been associated with cognition (Rivkin et al., 2013;Rollins et al., 2014;M von Rhein et al., 2014). D-transposition of the great arteries (d-TGA) is a form of CHD that is corrected by the arterial switch operation using cardiopulmonary bypass in early infancy. Brain abnormalities can be seen on MRI preand postoperatively, as well as in utero (Clouchoux et al., 2013;Licht et al., 2009;Limperopoulos et al., 2010;. Although white matter is most often affected, reduced cortical gray matter volume in the frontal and parietal lobes is also present several months after surgery Watanabe et al., 2009). Our group, using region-of-interest analyses, has shown that these measurable differences in brain structure have not returned to normal by adolescence in d-TGA patients (Rivkin et al., 2013;Rollins et al., 2014;Watson et al., 2016). While this regional approach to studying brain anatomical differences is helpful, it does not take into account the global organization of the brain, that is, its network structure. In a separate analysis of a subset of the same d-TGA adolescents, we have shown that differences exist in global brain organization based on white matter connectivity (Panigrahy et al., 2015).
To the best of our knowledge, analysis of gray matter connectivity has not been employed to discern related organizational networks in CHD patients of any age and may contribute to a more complete picture of the structural brain differences in this cohort. Here, we use a graph theoretical approach to analyze brain networks based on cortical thickness correlations to compare brain structure in a group of adolescents born with d-TGA corrected surgically in early infancy with that of typically developing control adolescents.

| Subjects
Adolescents in the d-TGA group were recruited from the Boston Circulatory Arrest Study, as previously described (Bellinger et al., 1995;Newburger et al., 1993). In brief, d-TGA subjects underwent the arterial switch procedure before 3 months of age between April 1988 and February 1992 at Boston Children's Hospital (BCH) (Bellinger et al., 1995;Newburger et al., 1993). Exclusion criteria included: known risk factors for brain disorders (e.g., history of closed head injury with loss of consciousness), any contraindication to acquisition of MRI data (e.g., metal implants), Trisomy 21, and adolescents with forms of CHD other than d-TGA requiring surgical correction.
The criteria used to recruit healthy control subjects were adapted from those of the NIH MRI study of normal brain development (Almli, Rivkin, McKinstry, & Brain Development Cooperative Group, 2007;Evans & Brain Development Cooperative Group, 2006). This study was approved by the Institutional Review Board and adhered to institutional guidelines and the Declaration of Helsinki. Parents provided informed consent and adolescents provided assent.

| MRI acquisition
Subjects were scanned on identical GE Twinspeed 1.5 Tesla (T) systems (General Electric, Milwaukee, WI, USA) with a quadrature head coil at either BCH or Beth Israel Deaconess Medical Center (BIDMC). The volumetric series for each subject was acquired using a Spoiled Proton Gradient Recalled (SPGR) sequence with parameters: TR/TE = 35 ms/6 ms, flip angle = 45 degrees, acquisition matrix = 256 × 256, FOV = 220 mm, slice thickness = 1.5 mm, with resultant voxel size = 0.86 × 0.86 × 1.5 mm 3 . The images were inspected by a radiologist to assure data quality and detect structural abnormalities (e.g., tumors, stroke, etc.).

| Cortical thickness calculation
Images were processed using Freesurfer v5.0 (A.A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital). The technical details are described elsewhere Fischl & Dale, 2000;. Briefly, MRI images are first partitioned into white matter, gray matter, and cerebrospinal fluid.
The outer pial surface of the brain is calculated, as is the surface comprising the white matter/gray matter junction. Cortical thickness is obtained by taking the distance between these two surfaces at every data point. Finally, the cortical surface is parcellated into distinct units based on gyral and sulcal anatomy (Fischl et al., 2004). The Desikan-Killiany atlas, which contains 34 regions per hemisphere, was used for the parcellation (Desikan et al., 2006). Mean cortical thickness was obtained for all regions for each subject.

| Network construction
All statistics were performed in R v3.2 (R Core Team, 2017), using functions in the packages igraph and brainGraph (https://cran.r-project. org/web/packages/brainGraph) (Csardi & Nepusz, 2006;Kolaczyk & Csárdi, 2014). First, a general linear model was specified for each brain region, with mean cortical thickness as the outcome variable and age, sex, and scanner location (BCH or BIDMC) as covariates. Next, Pearson correlation coefficients between the model residuals for all pairs of regions were calculated, creating an adjacency matrix of size 68 × 68 for each group.
The adjacency matrix of each group was binarized by thresholding and removing any correlations lower than the threshold. Negative correlations were not considered, as these likely do not represent real anatomic connections in the brain (Gong et al., 2012). To ensure equal network sizes for both groups, the thresholds were chosen to result in a specific density (i.e., ratio of the number of edges present in the network to total possible number of edges); a range of densities from 0.05 to 0.40 (step size: 0.01) was investigated. Since correlations were generally larger in the control group, the correlation threshold for each specified density was correspondingly larger (i.e., an equal threshold for each group would have resulted in a higher density in the control network). The networks created from these matrices were un-directed, un-weighted, and simple (i.e., no loops).

| Network metrics
Vertex-(i.e., region-) and graph-level metrics were calculated for both groups at each density. For visualization and group analysis purposes, a density of 22% was chosen, as this was the lowest density for which at least 95% of vertices were connected for both groups. This density is within the range used in several other studies (Bernhardt et al., 2011;He et al., 2008;Khundrakpam et al., 2013;Nie, Li, & Shen, 2013).

| Vertex importance
Vertex degree (the number of connections of a vertex), betweenness centrality (the number of shortest paths a vertex lies on), and nodal efficiency were used as measures of vertex importance. A vertex was considered to be a hub if its betweenness centrality was at least one standard deviation greater than the mean across all vertices for that density (Bernhardt et al., 2011;Tijms et al., 2013;Wang et al., 2013). Since regions classified as hubs may change across densities, we report those regions which were classified as hubs in at least half (i.e., 18/36) of the densities investigated. We also investigated rich-club organization for each group (Baker et al., 2015;Colizza, Flammini, Serrano, & Vespignani, 2006;van den Heuvel & Sporns, 2011). A rich-club is a group of vertices with high degree that are significantly more likely to be connected to each other compared to a set of equivalent random graphs. For both groups, we calculated the rich-club coefficient, ϕ, for each degree (from 1 to the maximum degree present in the network). We normalized ϕ (denoted ϕ norm ) by dividing by the average over a set of 1,000 equivalent random graphs (for each group and density). The random graphs were generated by randomly rewiring edges in the group-specific graphs for 10,000 iterations, keeping constant the graph's density and degree sequence (Maslov & Sneppen, 2002). Rich-club organization is considered present if ϕ norm > 1 for a range of degree thresholds. To determine a degree boundary, we used the "rich-core" algorithm of Ma and Mondragón (2015), which sorts the vertices from highest to lowest degree; the boundary is calculated as the local maximum of the function of degree rank and number of connections to higher-degree vertices (Ma & Mondragón, 2015). For simplicity, we used the maximum boundary value across subject groups.

| Network segregation and integration
Network segregation was assessed with three metrics. Modularity measures the strength of a given network partition. Higher modularity indicates that vertices belonging to the same network module (or community) are more connected to each other than they are to vertices of a different module. The Louvain algorithm was used to partition the networks into communities and compute the modularity (Blondel, Guillaume, Lambiotte, & Lefebvre, 2008). Degree assortativity is a related metric that measures the strength with which vertices of similar degree connect to one another; higher assortativity indicates that high-degree vertices are more likely to connect to other high-degree vertices compared to low-degree vertices. We also introduce lobe assortativity, which measures the number of inter-lobar connections relative to intra-lobar connections. This is equivalent to calculating the modularity of the network if it were a priori partitioned into the major lobes of the brain (i.e., frontal, parietal, temporal, occipital, insula, and cingulate). Higher values of lobe assortativity are present in networks with relatively fewer inter-lobar connections.
Small-worldness represents a balance between segregation and integration (Watts & Strogatz, 1998). Small world parameters clustering coefficient (C; the tendency of a vertex's neighbors to be connected to one another) and characteristic path length (L; the average of shortest path lengths between all vertices) were calculated, along with the average of each parameter from all random graphs (denoted C rand and L rand , respectively) for each group (Humphries & Gurney, 2008;Watts & Strogatz, 1998). The small-world index (σ) is the ratio of the normalized C to the normalized L (calculated as γ = C/C rand and λ = L/L rand , respectively), and a network is considered to possess the "small world" property if σ > 1.
Since the networks in this study are generated from correlations, they will tend to have a higher-than-expected level of clustering (Hosseini & Kesler, 2013;Zalesky, Fornito, & Bullmore, 2012). As a result, the random graphs generated by a simple rewiring procedure may not be entirely appropriate, as they will have very low clustering by design (Newman, 2010). Thus, as an alternative we generated (for each density and each group) 100 random networks while controlling for global clustering using a Markov Chain process (Bansal, Khandelwal, & Meyers, 2009). We then calculated an alternate small-world index, ω: (Telesford, Joyce, Hayasaka, Burdette, & Laurienti, 2011) Here, C latt is the mean clustering coefficient of a set of equivalent lattices; i.e., the graphs generated from the Markov Chain process in this case with L rand as previously described. A network is considered a small-world network if −0.5 ≤ ω ≤ 0.5; networks with ω closer to −1 are more similar to a lattice, and networks with ω closer to 1 are more similar to a random network. Networks with ω = 0 are considered to have a balance between global clustering and characteristic path length.

| Network closeness
Edge distances were calculated as the Euclidean distance in MNI coordinates (in mm) between centroids of pairs of connected regions (Alexander-Bloch, Vértes, et al., 2013;Bassett et al., 2008;He et al., 2007). Vertex distances were calculated as the mean distance of all edges connecting a given vertex to all other vertices (Alexander-Bloch, Vértes, et al., 2013). Similarly, characteristic path length (L) serves as a measure of the global closeness of a network.

| Asymmetry and hemispheric efficiency
A measure of asymmetry, the asymmetry index, was calculated as the difference in the number of left and right hemisphere intrahemispheric connections, divided by the average number of intrahemispheric connections of both hemispheres. An asymmetry index <0 indicates that the network has more intrahemispheric connections in the right compared to the left hemisphere. Additionally, we separated the group networks into isolated left and right hemisphere networks (Iturria-Medina et al., 2011;Li et al., 2015). We then calculated global efficiency (the inverse of the edges' shortest path lengths, averaged over all edges) and local efficiency (the efficiency of a subnetwork comprising a vertex's neighbors, averaged across all vertices) of each hemisphere separately.

| Network analysis
Between-group difference in the set of inter-regional correlations was determined by a two-sample t-test. Permutation testing was performed to assess group differences in global network measures (i.e., number of hubs, modularity, assortativity, clustering coefficient, characteristic path length, edge asymmetry, global and local efficiency, and vulnerability) and vertex-level measures (degree, betweenness centrality, and nodal efficiency). Each subject was randomly assigned to one of two groups (of the same size as the d-TGA and control groups), and then we followed the procedure for network construc-

| Subjects
Demographic and medical characteristics of the 92 d-TGA subjects and 49 control subjects have been described previously and is shown in Table 1; of note, the d-TGA subjects were significantly older and more likely to be male than control subjects (Watson et al., 2016).

| Cortical thickness covariance
The set of inter-regional correlations of cortical thicknesses was significantly greater in the control group compared to the d-TGA group

| Network hubs
The distribution of hubs differed both quantitatively and qualitatively between groups. Regions determined to be hubs in at least half of all densities are shown for both groups in Figure 1, overlaid onto the brain.
There were 9 hubs in the control group and 4 hubs in the d-TGA group (Table 2, top and bottom, respectively). Hubs in both groups tended to be in the right hemisphere (8/9 control; 3/4 d-TGA). Only the right superior frontal gyrus and right fusiform gyrus were common to both groups.
The d-TGA group rarely demonstrated more hubs than the control group across densities; however, the between-group difference in AUC for the number of hubs was not significant (p AUC = .08). Furthermore, there were no significant between-group differences in any measures of vertex importance (degree, betweenness centrality, and nodal efficiency).

| Rich-club organization
Rich-club organization was present in both groups for a range of degree thresholds. The "rich-core" degree boundary was k = 19 for the control group (containing 26 vertices, or 38% of the total network) and k = 16 for the d-TGA group (containing 27 vertices, or 40% of the total network).
When thresholding by the maximum across groups (i.e., k = 19), as in the distribution of hub regions, rich-club regions were predominantly righthemispheric in the control group (Figure 2, top). However, the distribution of rich-club regions in the d-TGA group was relatively symmetric across hemispheres (Figure 2, bottom). Based on permutation analysis, vertex degree was significantly greater in the control group compared to the d-TGA group in right pars triangularis and right superior temporal gyrus (p = .02 and p = .002 uncorrected, respectively). Qualitatively, the d-TGA group has more fronto-frontal and frontoparietal connections than the control group, and the control group has more fronto-and parieto-temporal connections. This is reflected in the lobe assortativity, which was greater in the d-TGA group at every network density tested, though the differences were not statistically significant ( Figure 4). Additionally, both degree assortativity and modularity were consistently higher in the d-TGA group, although the difference was not statistically significant across all densities (p AUC = .08 and p AUC = .28 for assortativity and modularity, respectively). Several other global network measures are also plotted against density in Figure 4.  suggesting these networks possess the small-world property and are closer to a lattice than a random network (Telesford et al., 2011).

| Network closeness
Characteristic path length was significantly lower in the control group across densities, indicating that regions are topologically closer to one another compared to the d-TGA group (p AUC = .026). Edge distance was significantly different between groups at several network densities (Figure 7, top), but not across all densities (p AUC = .14). In all cases, the median edge distance was higher in the control group than in the d-TGA group. Figure 7 (bottom) shows a histogram and density plot of edge distances for both groups at the network density of 0.22; while the between-group difference was not statistically significant at this density (p FDR = .11), the control group tended to have more longdistance connections than the d-TGA group.

| Asymmetry and hemispheric efficiency
Asymmetry was negative and lower in the control group at all network densities (Figure 4), indicating a greater number of intrahemispheric connections in the right hemisphere than in the left. This betweengroup difference was statistically significant at multiple network densities, including a density of 0.22 (p = .04). Across all densities, the between-group difference did not reach statistical significance (p AUC > .05). Additionally, a rightward asymmetry is evident in the control network's rich club, whereas the d-TGA network has rich-club regions distributed more evenly across hemispheres (Figure 2).
In the individual hemispheres, global efficiency was significantly higher in the control group than the d-TGA group for the right hemisphere at multiple densities (statistically significant across all densities, p AUC = .02) but was lower for the left hemisphere (statistically significant across all densities, p AUC = .03) (Figure 8, top). In the control group, global efficiency was higher in the right than the left hemisphere at every density; in the d-TGA group, global efficiency was lower in the right hemisphere or nearly equal to the left hemisphere. Local efficiency was higher in the left hemisphere for the d-TGA group, but lower in the right hemisphere (except for two densities), compared to the control group (Figure 8,bottom). The between-group difference in AUC of local efficiency was not significant for either hemisphere (p AUC = .52 and p AUC = .37 for the left and right hemispheres, respectively). In the control group, local efficiency was higher in the right hemisphere compared to the left at every density; in the d-TGA group, there was no consistent pattern present.

| Network robustness
At a density of 0.22, global vulnerability was significantly higher in the d-TGA group compared to the controls (group difference = −0.049; p = .027); a similar relationship was seen at multiple densities ( Figure 4). For the random failure and targeted attack analyses, at a density of 0.22, there were no significant differences between groups ( Figure 9). Both groups' networks were resilient against random failure of vertices and edges, and were similarly resilient to targeted attacks of vertices and edges.

| DISCUSSION
Using a graph theoretical approach to analyze cortical thickness networks, we found both qualitative and quantitative differences in brain organization between a group of adolescents with d-TGA and healthy controls. The d-TGA network tended to be more segregated, with higher modularity and assortativity. The control network had more long-range connections and was significantly more asymmetric toward the right hemisphere, both in terms of number of connections and global and local efficiency. The control network F I G U R E 2 Rich-club regions in the control and d-TGA groups at a density of 0.22. Rich-club regions and their connections are displayed for the control (top) and d-TGA (bottom) groups. The left column depicts a sagittal view of the left hemisphere, the center column depicts an axial view, and the right column depicts a sagittal view of the right hemisphere. Colors of individual vertices are based on membership in brain lobes -red: frontal; green: parietal; blue: temporal; orange: occipital; yellow: insula also had more hub regions, which, in turn, demonstrated more longrange connections distributed across both hemispheres and all the major lobes of the brain. Finally, global vulnerability was higher in the d-TGA network for a range of network densities, indicating a lower resilience to vertex "attacks." However, networks in both groups did display small-world and rich-club organization in addition to having similar modular structure. In summary, these results suggest that compromised brain networks in d-TGA adolescents possess less efficient and less integrated information processing as compared to control adolescents.
Congenital heart disease affects brain structure and development as early as the 3 rd trimester in utero causing features of developmental immaturity (Clouchoux et al., 2013;Dimitropoulos et al., 2013;Licht et al., 2009;Limperopoulos et al., 2010;Miller et al., 2007). Both preand postoperatively, the most common location of brain injury involves white matter, sometimes manifested as periventricular leukomalacia (PVL) (Beca et al., 2013;Gaynor, 2004). White matter damage early in development can adversely affect neuronal number and organization in gray matter (Inder et al., 1999;Leviton & Gressens, 2007;Volpe, 2009  Graph theory analysis of gray matter networks has been successful in differentiating other patient groups from healthy controls. An analysis of adults with temporal lobe epilepsy showed similar global organization but altered hub distribution relative to healthy controls (Bernhardt et al., 2011). Similarly, in adolescents with scoliosis, both the patient and control networks showed small-world organization, but an altered hub distribution; interestingly, there was an asymmetry in hub location matching the patient group's thoracic curves, suggesting a functional relevance of hub regions (Wang et al., 2013). In young adults with a history of childhood maltreatment, Teicher, Anderson, Ohashi, and Polcari (2014) found that regions involved in emotional processing possessed lower centrality in the proband group as compared to controls; the groups also had significantly different rich-club organization (Teicher et al., 2014). Taken (Khundrakpam et al., 2013). In the control group, there was strong right hemispheric occipito-and parietofrontal connectivity between hubs for most densities; these regions are broadly involved in attention and visuospatial functions. Altered asymmetry in cortical GM gyrification, particularly in frontal and temporal regions, has been found in HLHS fetuses (Clouchoux et al., 2013). The mechanisms contributing to such asymmetries remain unknown; however, the relatively reduced blood flow in the right internal carotid artery as compared to the left, may render cerebral tissue on that side more vulnerable to injury during periods of significant reduction in oxygen delivery in hypoxic or ishemic conditions (Bogren, Buonocore, & Gu, 1994). Finally, edge distances tended to be greater in the control group; long-distance connections may provide "shortcuts," and support improved integration among brain regions underlying different cognitive functions (Kaiser, 2011).
We, and others have reported that patients with d-TGA and other CHD have a high incidence of attention deficit hyperactivity disorder (ADHD), as well as deficits in executive function, attention, and visuospatial functioning (Bellinger et al., 2003(Bellinger et al., , 2011Razzaghi, Oster, & Reefhuis, 2015;von Rhein et al., 2015). The right hemisphere's importance in attention has been established for several decades, particularly in the study of patients with spatial neglect (Brain, 1941 Gainotti, Messerli, & Tissot, 1972;Mesulam, 1981).
Unilateral neglect is more commonly the result of right-hemisphere lesions and tends to be more severe compared to neglect consequent to left hemisphere lesions. In addition to study of stroke patients, abnormalities in right hemispheric structure and function have been found in children and adults with ADHD (Almeida et al., 2010;Epstein, Conners, Erhardt, March, & Swanson, 1997;Makris et al., 2007;Valera, Faraone, Murray, & Seidman, 2007;Vance et al., 2007). Network efficiency is representative of parallel information processing and fault tolerance (Latora & Marchiori, 2001

| CONCLUSIONS
Structural brain networks in adolescents with d-TGA, as measured by inter-regional cortical thickness correlations, differ in several aspects from a group of control subjects. Globally, both groups possessed small-world architecture, but network segregation tended to be higher in the d-TGA group. The control network was more asymmetric, containing more connections in the right hemisphere, and had a more efficient right hemisphere than the d-TGA group. Additionally, the d-TGA group had fewer long-range connections. Locally, the Red asterisks indicate a significant (p < .05) group difference based on permutation testing (N = 1,000); blue asterisks indicate a trend (p < .1). LH and RH indicate left and right hemispheres, respectively F I G U R E 9 Robustness analyses for the control and d-TGA groups at a density of 0.22. Relative maximal connected component sizes plotted as a function of: upper left, the percent of total edges removed in a random failure analysis; upper right, the percent of total vertices removed in a random failure analysis; lower left, the percent of total edges removed in a targeted attack analysis; and lower right, the percent of total vertices removed in a targeted attack analysis d-TGA group had fewer hub regions which were connected to closer brain regions than in the control group. Taken together, these differences in structural connectivity based on cortical thickness indicate differences in brain organization that likely relate not only to gray matter connectivity but also to underlying white matter differences we have seen in these patients. Further, these network differences constitute candidate measures to test for associations with the cognitive differences already identified between typically developing adolescents and those born with and treated for d-TGA.

ACKNOWLEDGMENTS
This work was supported by funding from the National Institutes of Health, National Heart, Lung, and Blood Institute (R01 HL77681 and HL41786), The Children's Heart Foundation (Chicago), the Farb Family Fund, and the Kostin Family Innovation Fund.

CONFLICT OF INTEREST
There have been no identified potential conflicts of interest.