Robust estimation of cortical similarity networks from brain MRI

Abstract

Structural similarity is a growing focus for magnetic resonance imaging (MRI) of connectomes. Here we propose Morphometric INverse Divergence (MIND), a new method to estimate within-subject similarity between cortical areas based on the divergence between their multivariate distributions of multiple MRI features. Compared to the prior approach of morphometric similarity networks (MSNs) on n > 11,000 scans spanning three human datasets and one macaque dataset, MIND networks were more reliable, more consistent with cortical cytoarchitectonics and symmetry and more correlated with tract-tracing measures of axonal connectivity. MIND networks derived from human T1-weighted MRI were more sensitive to age-related changes than MSNs or networks derived by tractography of diffusion-weighted MRI. Gene co-expression between cortical areas was more strongly coupled to MIND networks than to MSNs or tractography. MIND network phenotypes were also more heritable, especially edges between structurally differentiated areas. MIND network analysis provides a biologically validated lens for cortical connectomics using readily available MRI data.

Fig. 1. Estimation of MIND.

As input, we used the mesh reconstructions of the cortical surface generated from T1w MRI scans by FreeSurfer’s recon-all command. This surface can be described by a set of vertices (163,842 vertices per hemisphere for the fsaverage template). Each vertex was characterized by five structural MRI features: CT, SA, Vol, MC and SD. To estimate the similarity between cortical areas, we standardized each MRI feature across all vertices and then aggregated all the MRI metrics for all vertices within each cortical area (defined by a prior parcellation template) to form a regional multivariate distribution. We then compiled a pairwise distance matrix using a k-nearest neighbor density algorithm to estimate the symmetrized KL divergence, also known as Jeffrey’s divergence, between each pair of regional multivariate distributions. Finally, we transformed the KL divergence KL(a,b) for regions a and b to estimate the inter-areal MIND similarity, bounded between 0 and 1, with higher values indicating greater similarity. Illustrative distributions for regions a and b are shown as scatter plot matrices, with diagonal panels showing the marginal univariate distribution for five structural features and the off-diagonals showing each pairwise bivariate relationship. Bottom row: visualization of a group mean MIND similarity matrix and cortical surface maps of two elementary MIND network phenotypes—that is, edges between cortical nodes (the top 2% are shown here) and weighted nodal degree, calculated as the average edge weight for each of 318 cortical nodes defined by the DK parcellation.

Fig. 2. Cortical similarity connectomes: MIND networks and MSNs compared.

a,b, Illustrative MIND network and MSN from the same randomly sampled participant in the ABCD cohort. LH, left hemisphere; RH, right hemisphere; MS, morphometric similarity. c, Cortical surface maps of group mean weighted degree for the MIND networks and MSNs. d, Scatter plot representing the positive correlation between edge weights of the group mean MIND networks and MSNs. e, The distributions of pairwise correlations of network edges between subjects for MIND networks and MSNs, for all pairs of 10,367 subjects. f, The correlation between MIND networks and MSNs constructed using 1–5 additional random features of Gaussian noise (for n = 150 random subjects). The solid line represents mean values, with shading representing empirical 95% confidence interval (CI). g, The fraction of total inter-hemispheric connections represented at different network densities for both group mean MIND networks and MSNs. h, Parcellation consistency of MIND network phenotypes at nodal level (weighted degree) and at edge level. The left plot shows the correlation between weighted degree estimated by each of the possible pairs of three parcellation templates: DK, DK-318 and HCP. To calculate between-parcellation correlations, each vertex was assigned the weighted degree of the region within which it was located, for each parcellation, and the correlation was calculated between the resulting vectors of vertex-wise values. The right plot shows the correlation between 2,278 network edges calculated using the 68-region DK parcellation or by using the finer-grained DK-318 parcellation to estimate 50,403 edges and coarse graining (DK-318 interp.) to match the number of edges in the original DK network. i, The fraction of edges between two regional nodes of the same cytoarchitectonic class over a range of network densities. In g and i, shading represents the 95% CI estimated by population bootstrapping, and the solid line represents the mean over all bootstrapped results. In all panels, except as noted in h, the DK-318 parcellation was used to define 318 cortical regions of approximately equal volume.

Fig. 3. Structural similarity from MRI compared to axonal connectivity from tract tracing in the macaque brain.

a, Correlation between structural similarity edge weights, in MIND networks or MSNs derived from macaque MRI, and axonal connectivity edge weights derived from tract tracing in five connectomes: the {29 × 29}, {29 × 91}, {40 × 40} and {40 × 91} versions of the Markov parcellation, with the number of target and source regions, respectively, indicated in each case^,, and the whole-cortex connectome based on the separate RM parcellation^,,. The five connectomes contained n = 536, n = 1,615, n = 978, n = 2,229 and n = 3,267 edges, respectively. Shading indicates 95% confidence interval (CI). Asterisks indicate significantly increased correlation with tract-tracing data for MIND networks compared to MSNs, determined by bootstrapping network edges and performing a two-sided test on the difference in tract-tracing correlations: *P = 0.0018 and **P < 0.001, uncorrected. b, Correlation between tract-tracing {40 × 40} weights and MIND network or MSN edge weights over a range of tract-tracing network densities. Shading represents 95% CI. c, Scatter plot of the correlations between tract-tracing weights and MRI similarities (MIND or MS) for the set of edges connecting each regional node to the rest of the connectomes; thus, each point represents the correspondence between tract-tracing weights and structural similarity for each region in the {40 × 40} connectome (averaged for afferent and efferent connections; see Methods for details). The dashed line y = x highlights that similarities estimated by MIND were generally more strongly correlated with tract-tracing weights (above the line of identity) than morphometric similarities. d, Radar plots of the stability of the correlation between axonal connectivity, again from the {40 × 40} connectome, and structural similarity from MSNs or MIND networks, estimated over all possible input feature sets with one or two missing features. Missing features are noted at each radial position, with the radius from the center indicating correlation with tract-tracing weights. Best-case correlations for each type of structural similarity network estimated using all six MRI features are shown as dashed lines: MY, myelination (T1/T2 ratio).

Fig. 4. Predicting age from structural similarity and DWI tractography human brain networks.

a, Pairwise correlation between the edges of age-specific MIND networks, computed by averaging over subjects grouped by age in years. All age-specific group-level networks were highly correlated (r > 0.94) but nonetheless demonstrated a clear age-dependent progression—that is, age-specific group networks became less similar when compared across larger age gaps. b, Comparison of the performance of models trained to predict age using nodal weighted degree or edge weights of either MIND networks or MSNs in the HCP-D cohort (ages 8–21 years). The dataset was split into 10 training and test sets (90:10 ratio, all test sets non-overlapping); each point indicates the performance on one test set. The y axis indicates the partial Spearman correlation between predicted and true age, corrected for sex, Euler index (a proxy for scan quality) and a global connectivity coefficient (the sum over all matrix elements) to mitigate confound effects. Further training and evaluation details are provided in Methods. Significantly differential performance between MSNs and MIND networks (**P < 0.01, FDR corrected from paired two-sided t-tests) was observed for both edge-based (P = 0.004) and degree-based (P = 0.004) models. c, An analogous plot to b, comparing the performance of models trained on MIND networks, MSNs or DWI tractography connectomes in the HCP-YA cohort (ages 21–35 years). Networks in b used the DK-318 parcellation, whereas networks in c were based on the HCP 360-region parcellation to match the publicly available DWI tractography connectomes (**P < 0.01 and *P < 0.05, FDR corrected from paired two-sided t-tests). Exact P values for edge-based models were as follows: MSN versus DTI (P = 0.11), MSN versus MIND (P = 0.12) and DTI versus MIND (P = 0.01). Exact P values for degree-based models were as follows: MSN versus DTI (P = 0.25), MSN versus MIND (P = 0.008) and DTI versus MIND (P = 0.002). For b and c, box plots indicate data quartiles, and whiskers indicate the full data range, excluding outliers. NS, not significant.

Fig. 5. Structural similarity and transcriptional co-expression networks.

a, Gene co-expression networks (upper triangles) compared to MIND networks and MSNs (lower triangles). LH, left hemisphere. b,c, Correlation between the gene similarity network and the group MSN or MIND network, at the level of edges (b) and weighted nodal degrees (c). Shading represents 95% confidence interval (CI) of best-fit line. P values (uncorrected) were based on a two-sided spin test that generated a distribution of null correlations from random spatial network rotations (Methods). d, Stability of the correlation between structural and transcriptomic networks constructed from all subsets of the six postmortem brain gene expression datasets available. Also included are results on univariate MIND networks derived from cortical thickness and consensus DTI connectivity from the HCP-YA dataset. For each number of donors included, all combinations of transcriptional networks were constructed (without gene filtering), and the mean edge correlation was calculated for each network type. There were $\begin{array}{c}\left(\begin{array}{c}6\hfill \\ n\hfill \end{array}\right)\hfill \end{array}$ possible networks created for n = 1, 2…6 included donors. Shading indicates the minimum and maximum value of the association observed for each number of included donors. e, Cortical brain maps of MIND weighted node degree beside a weighted gene expression map derived from PLS analysis of the covariation between degree, or ‘hubness’, of MIND nodes and gene expression. PLS1 explained a significant amount of covariance (62%, P_spin = 0.01, two-sided spin test) between these two modalities. f, Cell type enrichment of the weighted, ranked gene list from PLS analysis of covariation between MIND degree and gene expression, using the median loading rank within one of seven sets of genes, each characteristic of a canonical class of cells in the central nervous system: excitatory neurons (Neuro-Ex, P = 0.049), inhibitory neurons (Neuro-In, P = 0.17), endothelial cells (Endo, P = 0.81), astrocytes (Astro, P = 0.48), microglia (Micro, P = 0.02), oligodendrocytes (Oligo, P = 0.42) and oligodendroglial precursor cells (OPC, P = 0.30). The zero position on the x axis represents the median position of all 15,633 genes (position 7,816), with negative ranks indicating genes that have expression positively correlated with MIND node degree—that is, overexpressed at highly connected MIND network hubs. P values were FDR corrected after a two-sided permutation test controlling for both spatial autocorrelation in the brain MRI data and correlation structure in gene expression (*P < 0.05; see Methods for details).

Fig. 6. Estimating heritability, h², of five regional MRI metrics and structural similarity network phenotypes derived from them.

a, Twin-based heritability ($h_{twin}^2$) of regional MRI metrics (SA, CT, Vol, MC and SD) and of weighted nodal degree for MIND networks and MSNs; each point represents one of 318 cortical areas. b, SNP-based heritability ($h_{SNP}^2$) for weighted degree of MIND networks and MSNs (n = 318 regions). Box plots in a and b indicate data quartiles, and whiskers indicate the full data range, excluding outliers. c, Scatter plot of twin-based versus SNP-based h² estimates for weighted degree of MIND networks and MSNs; each point represents a regional node in the cortical network. P values (uncorrected) were based on two-sided spin tests. Shading represents 95% confidence interval (CI) of best-fit lines. d, Cortical map of the regional $h_{twin}^2$ for MIND network degree. e,f, The strongest (e) and weakest (f) 1% of MIND edges and their corresponding $h_{twin}^2$ estimates. g, Scatter plot of $h_{twin}^2$ versus MIND network edge weights, with fitted line indicating significant negative correlation; each point is an edge in the network. The indicated P value was based on a two-sided spin test.