Multi-scale detection of hierarchical community architecture in structural and functional brain networks

TL;DR

This paper introduces a multi-scale community detection method for brain networks that captures hierarchical structures across different topological levels, improving understanding of structural and functional brain organization.

Contribution

The authors develop a multi-scale extension of modularity maximization for hierarchical community detection in brain graphs, applicable to synthetic and real neuroimaging data.

Findings

Structural brain networks exhibit more topological scales than functional networks.

The method identifies conserved community organization across hierarchical levels.

Multimodal analysis reveals scales where structural and functional communities align or differ.

Abstract

Community detection algorithms have been widely used to study the organization of complex systems like the brain. A principal appeal of these techniques is their ability to identify a partition of brain regions (or nodes) into communities, where nodes within a community are densely interconnected. In their simplest application, community detection algorithms are agnostic to the presence of community hierarchies, but a common characteristic of many neural systems is a nested hierarchy. To address this limitation, we exercise a multi-scale extension of a community detection technique known as modularity maximization, and we apply the tool to both synthetic graphs and graphs derived from human structural and functional imaging data. Our multi-scale community detection algorithm links a graph to copies of itself across neighboring topological scales, thereby becoming sensitive to conserved community organization across neighboring levels of the hierarchy. We demonstrate that this method allows for a better characterization of topological inhomogeneities of the graph's hierarchy by providing a local (node) measure of community stability and inter-scale reliability across topological scales. We compare the brain's structural and functional network architectures and demonstrate that structural graphs display a wider range of topological scales than functional graphs. Finally, we build a multimodal multiplex graph that combines structural and functional connectivity in a single model, and we identify the topological scales where resting state functional connectivity and underlying structural connectivity show similar versus unique hierarchical community architecture. Together, our results showcase the advantages of the multi-scale community detection algorithm in studying hierarchical community structure in brain graphs, and they illustrate its utility in modeling multimodal neuroimaging data.

Equations

$$\begin{gathered} Q=\sum _{ij}\left[\left(A_{ij}-\gamma P_{ij}\right)\right]\delta \left(g_i,g_j\right) \\ , \end{gathered}$$

$$\begin{gathered} Q=\frac{1}{2\mu }\sum _{ijxy}\left\{\left(A_{ij}-{\gamma }_x{P}_{ij}\right){\delta }_{xy}+{\delta }_{ij}{\tau }_{jxy}\right\}\delta \left(g_{ix},g_{jy}\right) \\ , \end{gathered}$$

$$\begin{gathered} Q=\frac{1}{2\mu }\sum _{ijlr}\left\{\left(A_{ijl}-{\gamma }_l{P}_{ijl}\right){\delta }_{lr}+{\delta }_{ij}{\omega }_{jlr}\right\}\delta \left(g_{il},g_{jr}\right) \\ , \end{gathered}$$

$$\begin{gathered} Q=\frac{1}{2\mu }\sum _{ijlrxy}\left\{\left(A_{ijl}-{\gamma }_l{P}_{ijl}\right){\delta }_{(lx,ry)}+{\delta }_{(ix,jy)}{\omega }_{jlrx}+{\delta }_{(il,jr)}{\tau }_{jxy}\right\}\delta \left(g_{ilx},g_{jry}\right) \\ , \end{gathered}$$

Full text

††journal: Neuroimage

1 Introduction

Hierarchical organization is a common motif in information processing systems (Bassett et al., 2010). The local embedding of similar processing units within groups that are then iteratively combined into larger and larger subsystems (Simon, 1991) provides a unique solution to the problem of balancing information segregation (within a group at a single scale) and integration (between groups across multiple scales) (Park and Friston, 2013). Such an organization is observed in very large-scale computer circuits and computing architectures (Ozaktas, 1992; Chen, 2016), cellular communication systems (Akyildiz et al., 2005), and social messaging systems (Moody and White, 2003). Across these various real-world information processing systems, hierarchical organization can additionally offer robustness to damage (Zhang et al., 2007; Helbing et al., 2006), and a complex and diverse repertoire of system functions (Hilgetag and Hütt, 2014; Valverde et al., 2015) that promote optimal and efficient information processing (Kinouchi and Copelli, 2006; Beggs, 2008) and transmission.

While these previous examples are all man-made systems, hierarchical organization is also present in natural information processing systems. A quintessential example is the brain – whether dissected in a nematode worm such as C. elegans, or non-invasively measured in a healthy adult human (Lohse et al., 2014). Importantly, hierarchies in these systems can occur in both time (Chaudhuri et al., 2014; Siebenhuhner et al., 2013) and space (Bassett et al., 2010), and can exist in the clustering of gene expression (Barabasi and Oltvai, 2004; Conaco et al., 2012; Arcila et al., 2014; Henzler et al., 2013) or the groupings of neuronal cell types in lamina and columns (Sümbül et al., 2014). Arguably one of the most complex types of architecture in the brain is hierarchical network architecture (Bassett and Siebenhuhner, 2013; Betzel and Bassett, 2016). Here, brain regions serve as nodes and structural or functional connections serve as edges in the network. Both structural and functional networks in the brain are critical conduits for information flow, processing, transmission, and cognitive computations more generally (Sporns, 2010). Importantly, hierarchical network structures can give rise to critical dynamics (Werner, 2009), where the behavioral repertoire of the neural system can be maximized with very few degrees of freedom (Chialvo et al., 2008). Yet, despite its fundamental importance, our understanding of the hierarchical organization in the brain remains limited, in part due to the fundamental nature of complex networks: they defy visual interpretation, and instead require computational algorithms to characterize.

Algorithmic methods to identify hierarchical network architecture must overcome the challenge of identifying embedded processing units within local groups. Particularly useful candidates include community detection methods, which currently dominate the study of brain networks (Rubinov and Sporns, 2010). Community detection techniques can take on many forms (Porter et al., 2009a; Fortunato, 2010), but perhaps the most common in the context of neuroimaging data is modularity maximization (Newman and Girvan, 2004). In this approach, nodes are partitioned into communities such that nodes within a community are more likely to connect to one another than expected in a random network null model (Newman, 2010). Importantly, the size of communities identified can be tuned by a structural resolution parameter, which titrates the relative difference between the real intra-community density and that expected in the null model (Reichardt and Bornholdt, 2006; Porter et al., 2009b). Therefore, sweeping across a range of resolution parameters offers glimpses into the hierarchical organization of the graph (Fenn et al., 2009, 2012); however, since the communities are identified independently at each point along the sweep, a secondary algorithm is required to track or link communities across topological scales, for example based on the similarity between communities in neighboring slices.

To address this limitation, we use a multi-scale community detection algorithm recently developed in applied mathematics (Mucha et al., 2010) to retrieve the underlying hierarchical organization of both artificial graphs and graphs representing human brain connectivity. We find that multi-scale community detection carefully preserves local information about the stability of sub-communities in the graph, enabling a thorough description of its hierarchical levels. Perhaps even more interestingly, we can uncover communities that remain stable across topological scales, and we can characterize their longevity and frequency. This approach offers unique advantages – such as sensitivity to community longevity – that extend more common methods that sweep across global topological scales (Fenn et al., 2009; Pons and Latapy, 2005) with independent estimates. Finally, we present methods for statistical assessment of the identified hierarchical communities, and we further offer an approach for the estimation of a consensus partition across the hierarchy.

We exercise and apply this multi-scale approach to better understand the putative hierarchical community organization of patterns of white matter pathways (SC) estimated from diffusion spectrum imaging, and of functional connections (FC) estimated from resting state functional magnetic resonance imaging. Across 60 healthy adult individuals, we show that SC is topologically heterogeneous, displaying a varying number of stable communities across brain regions. In contrast, we show that the hierarchical organization of FC is flatter, displaying a smaller number of stable communities across scales. Building on these observations, we probe the spatial embedding of communities in each modality separately, and then compare and contrast the modalities with one another. Our work offers a roadmap for the use of multi-scale community detection in revealing hierarchical network structure in structural and functional brain graphs, and in assessing their relationships to one another. In future research, this technique could be combined with multilayer approaches to better understand multimodal hierarchical architectures in health and disease.

2 Methods

2.1 Participants

Sixty participants (28 male, 32 female) were recruited locally from the Pittsburgh, Pennsylvania area as well as the U.S. Army Research Laboratory in Aberdeen, Maryland. Participants were neurologically healthy adults with no history of head trauma, neurological or psychological pathology. Participant ages ranged from 18 to 45 years old (mean age, 26.5 years). Informed consent, approved by the Institutional Review Board at Carnegie Mellon University and in compliance with the Declaration of Helsinki, was obtained in writing for all participants. Pittsburgh participants were financially compensated for their time.

2.2 MRI acquisition

All 60 participants were scanned at the Scientific Imaging and Brain Research Center at Carnegie Mellon University on a Siemens Verio 3T magnet fitted with a 32-channel head coil. An MPRAGE sequence was used to acquire a high-resolution (1 mm3 isotropic voxels, 176 slices) T1-weighted brain image for all participants. DSI data was acquired following fMRI sequences using a 50 min, 257-direction, twice-refocused spin-echo EPI sequence with multiple $q$ values ($TR=11,400ms$, $TE=128ms$, voxel size 2.4 mm3, field of view 231 $\times$ 231 mm, $b$-max 5000 s/mm2, 51 slices). Resting state fMRI (rsfMRI) data consisting of 210 T2*-weighted volumes were collected for each participant with a BOLD contrast with echo planar imaging (EPI) sequence (TR 2000 ms, TE 29 ms, voxel size 3.5 mm3, field of view 224 $\times$ 224 mm, flip angle $79$ degrees).

Head motion is a major source of artifact in resting state fMRI data (rsfMRI). Although recently developed motion correction algorithms are far more effective than typical procedures (Satterthwaite et al., 2013; Power et al., 2014; Pruim et al., 2015; Ciric et al., 2016), head motion was additionally minimized during image acquisition with a custom foam padding setup designed to minimize the variance of head motion along pitch and yaw directions. The setup also included a chin restraint that held the participant’s head to the receiving coil itself. Preliminary inspection of EPI images at the imaging center showed that the setup minimized resting head motion to 1 mm maximum deviation for most subjects. Only 2 subjects were excluded from the final analysis because they moved more than 2 voxels multiple times throughout the imaging session.

2.3 Diffusion MRI reconstruction

DSI Studio (http://dsi-studio.labsolver.org) was used to process all DSI images using a $q$-space diffeomorphic reconstruction method (Yeh and Tseng, 2011). A nonlinear spatial normalization approach (Ashburner et al., 1999) was implemented through 16 iterations to obtain the spatial mapping function of quantitative anisotropy (QA) values from individual subject diffusion space to the FMRIB 1 mm fractional anisotropy (FA) atlas template. QA is an orientation distribution function (ODF) based index that is scaled with spin density information that permits the removal of isotropic diffusion components from the ODF to filter false peaks, facilitating the resolution of fiber tracts using deterministic fiber tracking algorithms. For a detailed description and comparison of QA with standard FA techniques, see Yeh et al. (2013). The ODFs were reconstructed to a spatial resolution of 2 mm3 with a diffusion sampling length ratio of 1.25. Whole-brain ODF maps of all 60 subjects were averaged to generate a template image of the average tractography space.

2.4 Fiber tractography and analysis

Fiber tractography was performed using an ODF-streamline version of the FACT algorithm (Yeh et al., 2013) in DSI Studio (September 23, 2013 and August 29, 2014 builds). All fiber tractography was initiated from seed positions with random locations within the whole-brain seed mask with random initial fiber orientations. Using a step size of 1 mm, the directional estimates of fiber progression within each voxel were weighted by $80\%$ of the incoming fiber direction and $20\%$ of the previous moving direction. A streamline was terminated when the QA index fell below $0.05$ or had a turning angle greater than $75$ degrees. We performed a region-based tractography to isolate streamlines between pairs of regional masks. All cortical masks were selected from an upsampled version of the original Automated Anatomical Labeling Atlas (AAL) (Tzourio-Mazoyer et al., 2002; Desikan et al., 2006) containing 90 cortical and subcortical regions of interest but not containing cerebellar structures or the brainstem. This resampled version contains 600 regions and is created via a series of upsampling steps in which any given region is bisected perpendicular to its principal spatial axis in order to create 2 equally sized sub-regions (Hermundstad et al., 2013, 2014). The final atlas contained regions of an average size of 268 voxels (with a standard deviation of 35 voxels). Diffusion-based tractography has been shown to exhibit a strong medial bias (Croxson et al., 2005) due to partial volume effects and poor resolution of complex fiber crossings (Jones and Cercignani, 2010). To counter the bias away from more lateral cortical regions, tractography was generated for each cortical surface mask separately.

2.5 Resting state fMRI preprocessing and analyses

SPM8 (Wellcome Department of Imaging Neuroscience, London) was used to preprocess all rsfMRI collected from 53 of the 60 participants with DSI data. To estimate the normalization transformation for each EPI image, the mean EPI image was first selected as a source image and weighted by its mean across all volumes. Then, an MNI-space EPI template supplied with SPM was selected as the target image for normalization. The source image smoothing kernel was set to a FWHM of 4 mm, and all other estimation options were kept at the SPM8 defaults to generate a transformation matrix that was applied to each volume of the individual source images for further analyses. The time-series was up-sampled to a 1Hz TR using a cubic-spline interpolation. Regions from the AAL600 atlas were used as seed points for the functional connectivity analysis (Hermundstad et al., 2013, 2014). A series of custom MATLAB functions were used to do the following: (1) extract the voxel time series of activity for each region, (2) remove estimated noise from the time series by selecting the first five principle components from the white matter and CSF masks.

2.6 Data preprocessing

Both the DSI and BOLD data were used to construct $(N\times N,N=600~{}regions)$ structural and functional networks. We then studied the hierarchical community structure of these graphs using multi-scale community detection.

2.6.1 Functional network construction

Following prior work (Bassett et al., 2011b; Mantzaris et al., 2013; Bassett et al., 2013b, 2014), we estimated the dynamic functional connectivity between all region pairs using a wavelet coherence (Grinsted et al., 2004). We choose the wavelet decomposition based on its denoising properties (Zhang et al., 2016) and its utility in estimating statistical similarities between long memory time series such as those observed in resting state fMRI data (Achard et al., 2008). We observe two distinct bands of high coherence: $0.24-0.17Hz$ and $0.16-0.08Hz$. We focus on the $0.16-0.08Hz$ band due to known sensitivity to underlying neural activity (Hutchison et al., 2013). Coherence amplitudes were averaged over all frequencies and time points within the selected band to construct the average band-passed coherence for each pair of regions resulting in a single $N\times N$ functional connectivity (FC) adjacency matrix per subject (Chai et al., 2017).

2.6.2 Structural network construction

The individual subject’s structural connectivity (SC) matrix represents the fiber count between all region pairs. Commonly the fiber count values are normalized by the region size, so that the values of the SC matrix reflect the density of the white matter streamlines constructed between two regions (Hagmann et al., 2008; Bassett et al., 2011a; Gu et al., 2015; Betzel et al., 2016). However, we circumvented this issue by using the AAL600 atlas (Hermundstad et al., 2013, 2014), which was purposefully designed to contain similarly-sized regions. Due to the heavy tailed nature of the edge weight distribution (see Appendix Fig. 1) we applied a log transform to the edge weights $\textstyle(log(SC+1)/max(log(SC+1)))$ (see Fig. 3A) to increase the discriminibility of low edge weights and well as to increase the comparability to the edge weight distribution of functional connections.

2.7 Community detection

Common community detection algorithms can be used to partition a graph into clusters, where nodes tend to be more tightly connected to other nodes in their same cluster than to nodes in other clusters (Porter et al., 2009a; Fortunato, 2010). In the context of network data (or other relational data that can be represented as a network), we adopt common parlance and refer to these clusters as communities. A common heuristic for the identification of community structure in network data is the optimization of a quality function, which measures the relative density of the intra- versus inter-community edges (Newman and Girvan, 2004; Reichardt and Bornholdt, 2004). One particularly popular quality function is the modularity quality function (Newman, 2006), which can be defined as:

[TABLE]

where for a graph of $N$ nodes, $A$ is the $N\times N$ weighted adjacency matrix, the $ij^{th}$ element of the adjacency matrix indicates the weight of the connection between node $i$ and node $j$, the Kronecker delta $\delta\left(g_{i},g_{j}\right)=1$ if the community assignment of node $i$ and node $j$ ($g_{i},g_{j}$) are identical ($g_{i}=g_{j}$) and zero otherwise, $\gamma$ is the structural resolution parameter, and $P_{ij}$ is the expected weight of the $ij^{th}$ edge between node $i$ and node $j$ under a specified null model. Importantly, by maximizing this quality function, one can identify a partition of nodes into communities; however, identifying the optimal partition is NP-hard, and therefore the problem is usually solved with clever heuristics such as the Louvain-like locally greedy algorithm (Blondel et al., 2008). To account for the near degeneracy of the modularity landscape (Good et al., 2010), the algorithm is used to optimize the modularity quality function multiple times, and results are only reported that remain consistent over those optimizations (Bassett et al., 2013a).

The Newman-Girvan null model (Girvan and Newman, 2002) is the most commonly used null model in modularity maximization. It can be defined as: $P_{ij}=\frac{k_{i}k_{j}}{2m}$, where $k_{i}=\sum_{j}A_{ij}$ is the strength of node $i$ and $m=\frac{1}{2}\sum_{ij}A_{ij}$. In short, this null sets the expectations of an edge based on the strength of its nodes. The choice of structural resolution parameter $\gamma=1$ is common (Lancichinetti and Fortunato, 2011; Berry et al., 2011; Traag et al., 2011), however it only represents the community organization at a single topological scale (Bassett et al., 2013a; Lohse et al., 2014). Because graphs often display hierarchical organization, it is frequently useful to explore community structure in a graph over a range of values for $\gamma$ (Fenn et al., 2009, 2012; Bassett et al., 2013a; Lohse et al., 2014). When a graph has a particularly salient topological scale at which community structure is strongest, this parameter sweep can be used to identify the “optimal” structural resolution parameter value at which this community structure can be identified (Bassett et al., 2013a). However, for a graph that has hierarchical structure in which multiple topological scales are equally salient, this approach can fail to identify a single “optimal” structural resolution parameter value.

2.8 Hierarchical community detection

To study hierarchical community structure in graphs, we suggest a method based on optimizing the modularity quality function across all neighboring topological scales simultaneously. We achieve this by creating a multilayer network (Mucha et al., 2010), linking duplicates of the graph at each $\gamma$ value to the graphs at neighboring $\gamma$ values. For a schematic representation of the proposed multilayer graph, see Fig. 1. Formally, we define the multi-scale modularity quality function as:

[TABLE]

where the $ij^{th}$ element of the adjacency matrix $\mathbf{A}$ indicates the weight of the connection between node $i$ and node $j$, the Kronecker delta $\delta\left(g_{ix},g_{jy}\right)=1$ if the community assignments of node $i$ from scale $x$ and node $j$ from scale $y$ ($g_{ix},g_{jy}$) are identical ($g_{ix}=g_{jy}$) and zero otherwise, $\gamma_{x}$ is the structural resolution parameter at layer $x$, $P_{ij}$ is the expected weight of the edge between node $i$ and node $j$, $\tau_{jxy}$ is the topological scale coupling parameter which indicates the strength of the links between neighboring topological scales (as represented by layers), the total edge weight in the network is $\mu=\frac{1}{2}\sum_{jy}K_{jy}$, the strength of node $j$ in layer $y$ is $K_{jy}=k_{jy}+c_{jy}$, the intra-layer strength of node $j$ in layer $y$ is $k_{jy}=\sum_{i}A_{ij}$, and the inter-layer strength of node $j$ in layer $y$ is $c_{jy}=\sum_{x}\tau_{jxy}$.

Following prior work (Wymbs et al., 2012; Bassett et al., 2013a, 2015a; Papadopoulos et al., 2016), here we choose the expected values of the edge $P_{ij}$ uniformly for all edges as the average strength of all nodes: $P_{ij}$ equals some constant. This is sometimes referred to as the geographic null model (Bassett et al., 2015a; Papadopoulos et al., 2016). We choose the constant pragmatically, and separately for the structural matrices versus functional matrices. In the structural matrices, we noted that the $min(A_{ij})$ is relatively consistent across the subjects in our sample, since the smallest value of streamline count is 1 streamline, while $mean(A_{ij})$ is quite different across subjects in our sample. By contrast, in functional matrices, we noted that the $mean(A_{ij})$ was relatively consistent across subjects in our sample, while the $min(A_{ij})$ was not. In order to maintain the greatest sensitivity to structure that is conserved across subjects in the sample, we therefore chose $P_{ij}=min(A_{ij})$ for structural matrices and $P_{ij}=mean(A_{ij})$ for functional matrices. Note that this means that the exact value of $\gamma$ used across structural and functional matrices is not directly comparable, while its relative value is.

In this manuscript, we only examine the multi-scale community structure of the FC and SC graphs at a low value of the topological scale coupling parameter $(\tau=0.5)$ where the community organization of the neighboring topological scales exhibit relatively small dependencies on one another. While not the focus of this paper, it is important to note that the hierarchical community detection framework we describe and exercise can also be extended to other multi-layer and temporal graphs. In the latter case, our framework can be used to link communities across different temporal scales; see Appendix.A for details.

2.9 Statistics of multi-scale sommunities: Stability and consensus communities

The multi-scale community detection algorithm identifies many communities that could span several topological scales, each here represented as a layer in the multilayer network. In other similar multilayer contexts, it is crucial to be able to assess the stability of the identified communities across scales (Fenn et al., 2009), under the assumption that stable communities are of particular interest (Lambiotte et al., 2014). In our multi-scale framework, we measure the stability of individual nodes’ allegiance to their communities across scales: for node i the stability of its allegiance to community ’X’ is calculated as the number of $\gamma$ values where node i belongs to community ’X’ divided by the total number of slices (i.e., all structural resolution parameter values examined). Higher values of stability indicate that a node belongs to a single community across a greater number of layers, indicating its participation in a wider range of topological scales in the hierarchy.

Importantly, a node’s stability can be calculated as a function of the value of the structural resolution parameter, $\gamma$. For example, suppose node i is assigned to community ’X’ at $\gamma=1$ and community ’Y’ at $\gamma=2$. The stability of node i at the point $\gamma=1$ is then equal to the number of $\gamma$ values where node i belongs to community ’X’ divided by the total number of slices; by contrast, the stability of the same node i at the point $\gamma=2$ is equal to the number of $\gamma$ values where node i belongs to community ’Y’ divided by the total number of slices. Thus, in fact we can calculate a stability matrix that encodes the stability of each node at each value of the structural resolution parameter (e.g., Appendix Fig. 1A). In this matrix, highly diverse patterns of stability are indicative of topological heterogeneity in the graph. By constrast, less diverse patterns of stability are indicative of topological homogeneity in the graph.

Because we employ a heuristic to maximize the modularity quality function (Blondel et al., 2008), the identified partition of the multilayer network into multi-scale communities can change at each iteration (Good et al., 2010). To establish a robust, representative partition across these iterations, we perform the following steps: (i) we maximize the modularity quality function many times $(n=100)$ to adequately sample the modularity landscape, (ii) for every pair of brain regions, we calculate the average probability of two nodes appearing in the same community (which we refer to as the intra-layer community allegiance) from the multi-scale partitions for all layers, (iii) for every brain region, we calculate the average probability of it appearing in the same community across two neighboring layers (which we refer to as the inter-layer community allegiance) from the multi-scale partitions for all neighboring layers, (iv) we identify the nodes (and layers) with reliable inter-and intra-layer community allegiance by comparing average community allegiance values with that of a null model. The average community allegiance of the null model was generated from randomizing community labels from step (i). Then, (v) we create a consensus multi-layer graph where the values of the intra-layer and intra-layer edges correspond to the average community allegiance of the edges that were found to be significantly different from the null model. All the non-significant edges were removed from the multi-layer graph. Finally, (vi) the consensus partition is identified from the multi-layer consensus graph using the multi-layer community detection algorithm with parameter values $\gamma=1$, and $\omega=1$.

3 Results

3.1 Hierarchical community organization of synthetic graphs

To illustrate the method, we begin with a synthetic hierarchical graph that is constructed so as to contain clear community structure across a range of topological scales (Fig. 2A). Specifically, the graph displays identifiable community structure across four topological scales, with nested clusters of 3, 9, and 27 nodes. Heterogeneity is introduced by adding gradients in the values of edge weights such that not all clusters of a given size have the same average weight. To uncover the hierarchical community structure in this synthetic graph, we first applied an existing approach: a maximization of a single-layer modularity quality function (Newman, 2006) with the Newman-Girvan null model (Girvan and Newman, 2002) using a Louvain-like locally greedy algorithm (Blondel et al., 2008). Across different values of the structural resolution parameter ($\gamma$), we observe that the communities identified appear to change frequently, with different communities being present at different values of $\gamma$ (Fig. 2B). Prior work suggests that a reasonable method to choose the “optimal” value for the structural resolution parameter is to identify a range of $\gamma$ over which the community structure does not change appreciably (Fenn et al., 2009; Bassett et al., 2013a). Applying that approach to these data, one might identify $9<\gamma<15$ as a range of $\gamma$ values over which the community structure is relatively stable (Fig. 2B). Yet, the community structure present in this range of $\gamma$ values alone reveals little about the planted hierarchical organization and the topological inhomogeneities across nodes, as seen in Fig. 2A.

To overcome this limitation, we apply a multi-scale community detection method using the hierarchical algorithm described in the Methods section. We observe that the community structure displays branching across layers (or values of $\gamma$; Fig. 2C), meaning that a community in one layer can branch into two or more subcommunities in the next layer. Tracking the changes in a node’s community allegiance as a result of branching into subcommunities gives us information about the local topology of the synthetic graph. These results suggest that the multi-scale community detection technique can accurately uncover planted hierarchical communities in synthetic graphs. To gain further intuition regarding the performance of the method, we also apply the approach to three other synthetic graphs with differing architectures (Appendix Fig. 1). Again, we observe that the community structure displays branching across layers and that the number of branches is indicative of the number of local hierarchical scales in the synthetic graph. Moreover, across all synthetic graphs, we can observe that communities with higher-valued edge weights appear at higher values of $\gamma$ than communities with lower-valued edge weights. Together, these examples highlight the utility of the multi-scale community detection method for revealing hierarchically organized communities.

3.2 Hierarchical community organization of white matter structure in the human brain

For a given subject, the structural connectivity (SC) matrix generated from the fiber count between 600 brain parcels is sparse, with an average streamline count of $3.57$ and a standard deviation of $37.83$ (see Methods). Intuitively, this sparsity indicates that relatively few brain regions share direct fiber connections with one another. Moreover, the distribution of edge weights is heavy tailed (see Fig. 1), ranging from only a few streamlines per node pair to several thousand per node pair. To better visualize the architecture of the SC matrix, we apply a log transform $\textstyle(log(SC+1)/max(log(SC+1)))$ (see Fig. 3A). Then, we use the multi-scale community detection method to determine whether the SC matrix displays hierarchical community structure, and if it does, to characterize that structure both qualitatively and quantitatively.

In single subjects, we observe that the SC graphs display hierarchical organization where communities branch into smaller sub-communities across a range of topological scales (Fig. 3B). These characteristics of single-subject SC graphs are recapitulated at the group level. By performing consensus clustering, we can estimate a hierarchical decomposition that is most characteristic of all subjects within the group. We observe a similar hierarchical community structure, indicating a high degree of similarity (and low variance) across subjects (Fig. 3C). In the Appendix Figs. 1 and 2, we show that these group-level communities tend to be composed of spatially proximal, or connected, regions indicating that topological clustering is nontrivially related to spatial location. An interesting exception to this general trend is the existence of a few spatially distributed communities located in the fronto-striatal circuitry that bridge frontal cortex and the striatum.

The hierarchical community organization in the SC graphs can be described by characteristic curves of community number and size as a function of resolution. Specifically, as the structural resolution parameter value increases, we observe a rapid increase in the number of non-singleton communities (maximum average number of SC communities = 143.20 $\pm$ 7.3 std) and an analogous drop in the average community size (Fig. 3D). As the value of $\gamma$ increases farther, the number of communities branching into singletons increases, and thus the number of non-singleton communities decreases. Together, these trends indicate that communities begin to branch at low $\gamma$ values, and continue to branch as $\gamma$ increases (full width at half max = 42.49 $\pm$ 2.9, Fig. 3D). Far from haphazard, this global branching process is highly structured, with a large number of nodes maintaining their allegiance to hierarchical communities over a long range of $\gamma$ before branching (see Fig. 3B).

3.3 Hierarchical community organization of functional connections in the human brain

The functional connectivity (FC) values were calculated based on the average wavelet coherence between regional time series, and unlike the SC matrix values, they exhibit a normal distribution with an average value of $\approx 0.49$ and relatively small standard deviation $(\approx 0.04)$ (Fig. 1. C$\&$D). Intuitively, this narrow range of edge weights means that FC graphs will display hierarchical organization over a smaller range of $\gamma$ values. To better sample the FC hierarchical structure, we therefore assessed community organization over a smaller range of $\gamma$ values, with $\gamma\in[0.95,1.7]$ in the FC matrices as opposed to the $\gamma in[0.0133,1]$ used for the SC matrices. Importantly, these ranges were chosen pragmatically so as to map out the entire curve from partitions containing a single community (lowest value in the $\gamma$ range), to partitions containing only singletons (highest value in the $\gamma$ range).

In single subjects, we observe that the FC graphs display hierarchical organization where communities branch into smaller sub-communities across a range of topological scales (Fig. 4B). Although multi-scale communities can be detected reliably in single subjects, the group-level consensus reveals less robust hierarchical organization (Fig. 4C). Indeed, at the group level, communities extend over smaller $\gamma$ ranges before branching (note the speckled nature of the community allegiance matrix). This observation indicates that there is relatively low inter-subject similarity of the hierarchical communities observed in FC graphs. In fact, a large number of brain regions (especially subcortical regions) fail to display any significant intra-layer community allegiance over most topological scales (see Appendix Figures 3 and 4).

The hierarchical community organization in the FC graphs can be described by characteristic curves of community number and size as a function of resolution. Consistent with the trends observed in the SC graphs, as the structural resolution parameter value increases, we again observe a rapid increase in the number of non-singleton communities (maximum average number of FC communities = 94.44 $\pm$ 20.66 std) and an analogous drop in the average community size (Fig. 4D). As the value of $\gamma$ increases farther, the number of communities branching into singletons increases, and thus the number of non-singleton communities decreases. Again, the global branching process is highly structured, with a large number of nodes maintaining their allegiance to hierarchical communities over a long range of $\gamma$ before branching. Nevertheless, in comparison to the SC graphs, the FC graphs display this hierarchical community organization over a smaller range of $\gamma$ values (full width at half max = 28.86 $\pm$ 7.1, Fig. 3D).

3.4 Heterogeneity in the hierarchical community organization of functional and structural brain networks

Next we aim to explicitly characterize similarities and differences in the hierarchical community organization of functional and structural brain networks. We begin by focusing on the notion of community longevity or stability across topological scales, and we estimate the average number of stable communities that each node belongs to across layers. Such a computation depends on first choosing a mathematical definition of what constitutes a “stable” community. Pragmatically, we choose a parametric definition in which a stable community is defined as a community that exists across more than $x$% of the $\gamma$ range studied. We refer to the $x$% as a stability threshold. Intuitively, graphs with pronounced multi-scale hierarchical organization display many stable communities per node across a wide range of stability thresholds, while graphs with weak multi-scale hierarchical organization display very few stable communities per node. In addition to the number of stable communities, it is also of interest to quantify the variance of this number across nodes in the network. Consequently, graphs with a large variance of these values across nodes are characterized by greater topological heterogeneity than graphs with a smaller variance of these values across nodes.

Applying these analyses and statistics to graphs extracted from imaging data, we observe that the brain’s structural and functional connectivity graphs are indeed hierarchical, with the vast majority of nodes displaying stable allegiance to communities over more than one topological scale (Fig. 5). First considering only SC graphs, we observe that at low stability thresholds, each brain region is allied to approximately $8$ communities across topological scales, while at higher stability thresholds, a brain region may only be allied to $1$ community. Across brain regions, we also observe high variance; at low stability thresholds, the number of communities to which a region allies ranges from approximately $5$ to approximately $12$, indicating a high degree of heterogeneity in the SC graphs. Finally, we observe a marked similarity between the community stability curves at the subject level and at the group level (Fig. 5A$\&$B), providing further evidence of inter-subject similarity of hierarchical community structure in SC graphs.

Next, considering the FC graphs, we observe that brain regions exist in a smaller number of stable communities across layers (Fig. 5C$\&$D). The comparison between the FC and SC subject-level results demonstrate that although the average curves appear similar in shape, the average number of stable communities is approximately 1.5 times higher in SC than FC graphs. These observations suggest that FC graphs display a flatter hierarchical community organization, characterized by a smaller number of topological scales. Across brain regions, we again observed relatively high variance; at low stability thresholds, the number of communities to which a region allies ranges from approximately $3$ to approximately $8$, indicating a high degree of heterogeneity in the FC graphs. The community stability curves also display a marked difference at the subject and group levels, again providing evidence of inter-subject variability of hierarchical community structure in FC graphs.

3.5 Reliable detection of a region’s consistent allegiance to communities across topological scales

In the previous sections, we first observed the hierarchical nature of community structure in structural and functional brain graphs, and then we determined the number of topological scales that characterize each node’s community allegiance profiles. In this section, we seek to better understand the fine-scale features of the multi-scale network model that allow for reliable estimation of these topological scales in individual brain regions. Importantly, the multi-scale network model is not simply an agglomeration of weighted adjacency matrices. Rather, it explicitly stitches these matrices together with inter-layer connections that link a brain region in one layer to itself in the preceding and following layers. These inter-layer links allow for the quantitative assessment of communities across layers in a statistically principled manner.

Explicit inter-layer links within the multi-scale model motivate an effective description of inter-layer consistency in a node’s allegiance to communities. In particular, over the large number of optimizations of the modularity quality function that must be performed to adequately sample the underlying landscape, it is of interest to quantify how frequently a node remains in the same community across two adjacent layers (here representative of topological scales). We define a reliable inter-layer association as occurring when a node remains in the same community across two adjacent layers for a greater number of optimizations than expected by chance (see Methods). We observe that a large number of reliable inter-layer associations can be identified in structural brain graphs at both the subject and group levels. This is true particularly for lower values of $\gamma$, indicating the presence of a $\gamma$ range over which hierarchical community assignments can be reliably detected (Fig. 6A). While a number of cortical and most of the subcortical structures including dienchephalone and limbic system display reliable inter-layer association across a small range of $\gamma$ increments at the group-level, several bilateral clusters in the frontal, parietal, and temporal cortex display a notably higher range (Fig. 6B). In FC graphs, we observe reliable inter-layer associations over a much smaller range of $\gamma$ at the subject level (Fig. 6 A). At the group level, we not only observe that reliable inter-layer associations occur over a small range of $\gamma$, but also that there are very few reliable associations at all (Fig. 6C). Together, these findings underscore both the flatter hierarchical nature of FC graphs and the greater inter-subject variability in comparison to SC graphs. Yet, the identified regions from the SC(FC) graphs with reliable inter-layer associations is a reflection of the fact that networks of structures with similar and/or strong structural (functional) connectivity are present across subjects.

3.6 Homogeneity versus heterogeneity of topological scales in structural and functional brain graphs

The previous results indicate that we can reliably detect hierarchical community structure in structural and functional brain graphs, and that the two types of graphs display differing degrees of topological heterogeneity. To better understand this heterogeneity, particularly across nodes in the network (or regions in the brain), we examine the stability of communities more closely. Specifically, for node i, we measure the stability of its allegiance to community ‘X’ by calculating the fraction of layers in which node i belongs to community ‘X’ (see Methods). This formulation allows us to define a stability matrix by replacing the community label with the stability of the node’s allegiance to the community (Appendix Fig. 1A-B). This matrix quantifies the stability of a node at a given structural resolution parameter value. In this matrix, highly diverse patterns of community allegiance stability are indicative of topological heterogeneity in the graph. By contrast, less diverse patterns of community allegiance stability are indicative of topological homogeneity in the graph.

To quantity homogeneity versus heterogeneity, we decomposed the stability matrix using a principle component analysis, such that each component indicated a coherent pattern of communities across scales (Appendix Fig. 1C-D). Intuitively, graphs with greater topological heterogeneity require a larger number of principle components to explain a given amount of variance in the community stability matrix compared to more homogeneous graphs. In both SC and FC graphs, we observed that a handful of components explained most of the variance in the community stability matrix (Appendix Fig. 1E-F). In SC graphs, eight principle components explain more than $95\%$ of the variance in the stability matrix. The first component is marked by the stability profile of the singleton nodes, and the second component highlights the lower half of the $\gamma$ range where most of the larger communities reside. In FC graphs, only five principle components explain more than $95\%$ of the variance in the stability matrix. The group-level analysis shows that on average significantly ($t$-test, $p<0.001$) smaller number of principle components ($4.86\pm 0.78$) explain subjects FC stability matrices compared to that of the SC graphs ($9.49\pm 1.26$). These results provide converging evidence that SC graphs display a greater topological heterogeneity while FC graphs display a greater topological homogeneity in hierarchical community structure.

3.7 Comparison between the hierarchical community organization of the brain’s structural and functional connectivity

In previous sections, we demonstrate that the hierarchical organization of the SC and FC graphs differ both in terms of the presence and stability of communities across topological scales. Yet it is also important to note that the identified communities across modalities do share some similarities, perhaps supporting the notion that structure provides the scaffold for emergent functional dynamics. To better understand the similarities between the hierarchical community organization of SC and FC graphs, we begin by summarizing the community structure at each $\gamma$ value in each modality as an $N\times N$ allegiance matrix, where the $ij^{th}$ element indicates the fraction of times that node $i$ and node $j$ are placed in the same community over all optimizations of the multilayer modularity quality function.

Next, we calculate the Pearson correlation coefficient between the allegiance matrix of SC at a given $\gamma$ and the allegiance matrix of FC at a given $\gamma$, for all possible $\gamma$ pairs (Fig. 7A). This approach enables us to capture the degree to which densely connected communities of brain regions in FC similarly echo their underlying SC, thereby providing insight into the structural drivers of global dynamics. We observe that SC and FC communities show high similarity at medium topological scales (Fig. 7B), suggesting that it is not simply the case that densely structurally connected regions are also functionally connected. Instead, these results suggest that medium-sized bundles that link the densely connected (and commonly local) brain regions allow global functional synchronization between relatively large ensembles (Fig. 7C).

3.8 Explicit multimodal investigation using a multiplex, multi-scale graph

While the comparisons thus far between structural and functional brain graphs have been illuminating, it is natural to ask whether there is a more principled and model-based approach to comparing the two modalities within the multilayer framework. Indeed, the multilayer framework does allow additional graphs to be interconnected along distinct dimensions. Thus, it is possible to construct a graph where one dimension hard-codes topological scale (as done throughout the earlier sections of this paper), and a second dimension that hard-codes imaging modality (e.g., structural connectivity and functional connectivity).

Here we construct exactly this multiplex graph to more formally study the multi-scale nature of both the structural and functional connectivity matrices within the same model (Appendix Fig. 1A). In addition to inter-scale links $\tau$, this model also contained inter-modality links $\kappa$ that link a node in one scale and one modality to itself in the same scale in a different modality. We optimize the modularity quality function in this multiplex case to identify the hierarchical community structure of the SC and FC graphs. Importantly, as $\kappa$ is tuned down, community structure is allowed to be independent in the SC and FC graphs (Appendix Fig. 1B). In contrast, when $\kappa$ is tuned up, community structure is forced to be consistent across the two types of graphs (Appendix Fig. 1F). In other words, by employing higher values of $\kappa$, we are able to extract community structure that is most representative of the graphs in both imaging modalities.

Interestingly, we observe that this cross-modality community structure appears more similar to the community structure of the SC graphs when they were studied independently, than to the community structure of the FC graphs when they were studied independently. This phenotype can occur in community detection when the community structure in one graph is stronger than the community structure in the other graph. To investigate and more thoroughly quantify this observation, we study the similarity between the allegiance matrices of the multiplex SC-FC graph, the FC graph alone, and the SC graph alone, as a function of the topological scale ($\gamma$ value) at which they were constructed (Appendix Fig. 2). We observe that the hierarchical structure of the multiplex SC-FC graph is similar to that of the FC graph alone only in a narrow range of topological scales, consistently across $\kappa$ values (Appendix Fig. 2A-B). In contrast, we observe that the hierarchical structure of the multiplex SC-FC graph is similar to that of the SC graph alone across a wide range of topological scales, and consistently across $\kappa$ values (Appendix Fig. 2C-D). These results suggest that the joint optimization is more heavily influenced by the hierarchical community structure in the SC graph than it is by that of the FC graph.

Importantly, these results are reported over a single subject, and thus it is critical to ask to what degree these insights hold over the entire participant cohort. To address this question, we perform the same set of analyses, but instead of using the single-subject allegiance matrices, we use the group-level allegiance matrices. In general we observe consistent results at the group scale. Specifically, the hierarchical community structure of the muliplex SC-FC graphs at high $\kappa$ values are consistently reminiscent of the SC structure (highest observed correlation approximately $r=0.85$), along a range of topological scales (Fig. 8A). And they are reminiscent of the FC structure to a weaker degree (highest observed correlation approximately $r=0.41$), along a much narrower range of topological scales (Fig. 8B). Notably, the similarity between either SC or FC and the multiplex SC-FC graphs is higher than between SC and FC alone (highest observed correlation approximately $r=0.25$). The FC and SC communities at coarse topological scales show the highest similarity (Fig. 8 C-F), however they diverge at higher topological scales as the SC communities branch into smaller local communities (Appendix Fig. 1), whereas the FC communities at higher topological scales are more spatially distributed. These results confirm at the group level that the joint optimization is more heavily influenced by the hierarchical community structure in the SC graph than it is by that of the FC graph. In addition, regions that display comparable community allegiance between FC and SC graphs such as the subcortical nodes and some clusters within the medial frontal and medial occipital cortices are also identifiable in the multiplex SC-FC graph’s multi-scale communities as they maintain their community allegiance across $\gamma$ increments. Overall the observed similarity between the multiplex SC-FC graph’s communities and the FC and SC communities provides converging evidence that both modalities share major organizational features. Nevertheless the multiplex SC-FC graph’s communities are comparable to the original FC and SC communities at different hierarchical scales, which highlights the differences in the hierarchical community organization of the SC and FC graphs.

4 Discussion

The human brain is a complex system that can be fruitfully represented as a graph or network in which brain regions correspond to network nodes and structural or functional connections between regions correspond to network edges (Bullmore and Sporns, 2009). Recent observations have pointed to the fact that both structural and functional brain networks may have community structure (Sporns and Betzel, 2016): the presence of densely interconnected groups of regions that might support specific cognitive functions (Power et al., 2011; Yeo et al., 2011; Meunier et al., 2009). Moreover, evidence suggests that these communities exist over multiple topological scales (Bassett et al., 2010), with larger communities potentially being composed of smaller communities (Betzel and Bassett, 2016). Yet a comprehensive characterization of this putative hierarchical community structure in structural and functional brain graphs has remained difficult largely due to inadequacies in existing analytical paradigms and computational tools. Here we address these limitations by exercising a multi-scale community detection algorithm (Mucha et al., 2010), and applying it to both structural brain networks estimated from diffusion imaging and functional brain networks estimated from resting state fMRI. Using novel statistics including community stability and inter-scale reliability, we show that structural brain graphs display a wider range of topological scales than functional graphs. We also illustrate the utility of this method in examining multimodal graphs that combine both structural and functional connectivity information. Our work illustrates the practical utility of multi-scale community detection in revealing hierarchical community structure in brain graphs, and opens the door for future investigations of this structure in both health and disease.

4.1 Detecting multi-scale community structure

Characterization of multivariate dependencies across spatio-temporal scales is critical for a fundamental understanding of observable dynamics across systems as diverse as the climate system (Steinhaeuser et al., 2012) and the human brain (Betzel and Bassett, 2016). The multi-scale community detection algorithm that we exercise here reveals the hierarchical community organization of a graph by assuming dependence between neighboring topological scales (Mucha et al., 2010). A marked advantage of this approach compared to conventional single-scale algorithms is that it provides a statistically principled answer to the question: “Is a community at one scale the same as or different from a community at another scale.” Perhaps even more importantly, the approach provides an estimate of the stability of local topological structure, and therefore a pragmatic means of identifying model parameter values that maximize the consistency of locally stable communities across several topological scales. These local estimates of community stability (unlike the global measures of community stability that have been previously defined (Delvenne et al., 2010; Pons and Latapy, 2005; Arenas et al., 2008; Ronhovde and Nussinov, 2009; Karrer et al., 2008)) are robust to topological heterogeneities (Danon et al., 2006) in the form of communities of different sizes with different average edge weights. When applied to human brain networks, we find that the local community stability estimates allow characterization of communities that are stable across a range of topological scales. Taken together, our study offers not only a methodological approach to studying hierarchical community structure in graphs, but also a set of statistical methods to characterize the observed structure and to compare it across different classes of graphs, either treated separately or combined into a multiplex model.

4.2 Multi-scale community structure in the human brain’s white matter architecture

Structural brain graphs estimated from diffusion imaging data tend to be sparse, and the edge weight distributions tend to be heavy-tailed (Lohse et al., 2014; Bassett et al., 2011a; Hagmann et al., 2008). These characteristics can occur when a complex topology is embedded into a 3-dimensional space (Bassett et al., 2010), in such a way as to enhance the efficiency of information transmission (Bullmore et al., 2009) while decreasing the cost of the wiring (Bullmore and Sporns, 2012). Interestingly, prior work has also offered initial evidence that the complexity of structural connectivity is in part due to the fact that it is organized in a hierarchically modular fashion (Lohse et al., 2014), which is thought to support its information processing capabilities (Simon, 1991). Here we use a principled mathematical modeling approach to more exactly identify hierarchical community structure in structural brain graphs. Our results demonstrate that structural connectivity is characterized by heterogeneous multi-scale communities, and by nodes that form stable hierarchical communities across a range of topological scales. Multi-scale communities appear to be largely consistent across different subjects and also tend to be spatially localized. The observation of both regional heterogeneity and diverse topological scales indicates that the application of single-scale community detection techniques is likely to produce an overly-simplified picture of the brain’s organization.

While the majority of multi-scale communities were spatially localized, communities in basal ganglia-thalamo-cortical circuitry were not. Over several topological scales, subcortical structures including basal ganglia (mainly putamen, palladum, and caudate) and anterior thalamus as well as several frontal neocortical areas were identified within the same community. While frontal and subcortical structures are spatially distributed, it is commonly known that much of the cortex (including both allocortex and isocortex) projects to the striatum (Swanson, 2000), although not all projections are entirely reciprocal. These consistent projections can manifest as structural motifs that are accessible to community detection algorithms. Other complementary algorithms based on tools from algebraic topology, including the notions of persistent homology (Giusti et al., 2016), have demonstrated that basal ganglia-thalamo-cortical connections are among the few most common cycles identifiable in structural brain graphs (Sizemore et al., 2016). These cycles and motifs are known to play key roles in rhythmic gain control, and in the gating and integration of information across the brain (Womelsdorf et al., 2014; Rajan et al., 2016). For example, basal ganglia influence cortical states and behavior via dopaminergic inputs to thalamus, thereby enabling the integration of information characteristic of reinforcement learning (Schultz et al., 1997; Schultz, 1998). Indeed, basal ganglia input is modulatory and also serves to gate higher-order relay signals that are propagated through cortico-thalamo-cortical loops (Sherman and Guillery, 2006). Together, the unique role that the basal ganglia-thalamo-cortical pathways serve in driving brain states is made possible through the unique structural fingerprint of subcortical regions.

4.3 Multi-scale community structure in the human brain’s resting state functional connectivity

The pattern of phase-locking between regional BOLD time-series over a period of several minutes demonstrates that many regions display high functional connectivity with one another (Yaesoubi et al., 2015). While the resulting graph is relatively dense and homogeneous, it nevertheless displays some amount of hierarchical community structure. In single subjects, the hierarchical consensus analysis revealed a small range where multi-scale community structure is reliably identified. Interestingly, the group-level consensus analysis showed that the majority of brain regions failed to produce reliable inter-layer links across individuals, indicating the high degree of inter-subject variance in hierarchical community organization. These results are not entirely unexpected in light of the mounting evidence for both inter-session and inter-subject variability in resting state functional connectivity (Gonzalez-Castillo et al., 2014), particularly that located in heteromodal association areas (Mueller et al., 2013; Finn et al., 2015). Speculatively, it is possible that some of this inter-subject variability is due to the fact that these heteromodal association areas are more susceptible to and likely influenced by environmental factors, a fact highlighted by research on the postnatal period where they display protracted development during a time period of high plasticity (Mueller et al., 2013; Brun et al., 2009; Zilles and Amunts, 2013). Thus, the observed inter-subject variability in the multi-scale community structure of functional brain graphs could provide important fodder for a fundamental understanding of the principles of brain wiring, evolution, and ontogenetic development (Zilles and Amunts, 2013).

4.4 A comparison of hierarchical community structure in functional and structural graphs

It has long been observed that resting state functional connectivity shows statistically similar patterns to those observed in underlying structural connectivity (Hagmann et al., 2008; Vincent et al., 2007; Sporns, 2013). Yet, we have little understanding of how exactly the anatomical connections gives rise to observed functional interactions (Hermundstad et al., 2013, 2014). Evidence suggests that the relationship between structure and function is likely quite indirect, with the broader network adjacent to direct structural paths being critical to healthy dynamic couplings between brain regions (Goñi et al., 2014; Becker et al., 2015). Our work supports this notion by demonstrating that the highest similarity in hierarchical community structure between the two modalities is found between the fine topological scales of the functional graph and the relatively coarse topological scales of the structural graph, which takes into account broader anatomical network organization. The differential scales of function and structure that map onto one another can in part be explaind by the observation that structural graphs on average display $\approx 1.5$ times more topological scales than functional graphs. Together these results suggest that (i) functional connectivity dynamics are not strictly bound to or constrained within direct anatomical projections, but instead extend to spatially distributed circuits, and (ii) structural connections as estimated by white matter tractography display hierarchical community structure that can support long range functional coupling (Werner, 2009; Valverde et al., 2015).

It is important to note that although we found similarities between the hierarchical community structure of functional and structural graphs at different topological scales, the overall magnitude of the similarity was relatively small. One fundamental property of brain connectivity that might explain this relative independence of structural and functional graphs is temporal dynamics (Mattar et al., 2016; Hutchison et al., 2013; Bassett et al., 2011b; Preti et al., 2016). Indeed, while we have here studied a static functional graph that represents patterns of co-activation over several minutes, in reality the brain displays time-varying patterns of functional connectivity (Calhoun et al., 2014) that can track changes in cognitive processes (Braun et al., 2015; Mattar et al., 2015; Chai et al., 2016) and behavior (Vatansever et al., 2015; Gerraty et al., 2016; Bassett et al., 2015b). The nature of these dynamics suggest that the answer to the question “how does structure constrain function?” depends nontrivially on the time scale of the function (or functional connectivity pattern) in question. Indeed, recent work has demonstrated that synchronization patterns in hierarchical modular structures such as the brain may appear at different topological scales depending on the time scale of their interaction (Arenas et al., 2006; Villegas et al., 2014; Betzel et al., 2013; Betzel and Bassett, 2016; Honey et al., 2007; Zhou et al., 2006). Therefore, comparing the topology of the multi-scale functional and structural graphs can prove useful for understanding the properties of anatomical projections that are critical for the emergence of multi-scale functional patterns.

4.5 Explicit multimodal models of the brain’s hierarchical community structure

While the method that we propose and exercise is applicable to brain graphs constructed from a single imaging modality, it is also flexible and generalizable to questions that require fusion of brain graphs constructed from two or more modalities. We illustrate the method’s utility in this class of problems by exploring the joint optimization of modularity across a multiplex network composed of both functional and structural layers (Kivelä et al., 2014). This construction enables us to highlight the characteristics of community structure that are echoed across the two imaging modalities: diffusion imaging and resting state functional MRI. Our results reveal that at high inter-modal coupling strengths, the community organization of the multi-modal graph merges across modalities to form a hybrid structure. Similarity analysis revealed that the multi-modal graph shows highest similarity to functional graphs at relatively fine topological scales and to structural graphs at relatively broad topological scales. The narrow range of topological scales at which this hybrid structure was identified highlights the fact that community structure in functional dynamics extends beyond direct white matter projections. This work complements prior efforts to bridge functional and structural connectivity patterns using tools and techniques that span the purely qualitative and the exquisitely quantitative: these approaches include direct superposition, fiber tracking from functional parcels, and regression analysis relating functional, anatomical, and behavioral data (Rykhlevskaia et al., 2008). One recent study identified communities separately for each modality, and then subsequently maximized a cross-modularity function to identify a community partition shared by structure and function at medium-to-fine topological scales (Diez et al., 2015). Our work extends these findings by assessing hierarchical community organization in a multiplex network composed of both functional and structural layers. Future work could seek a better understanding of how joint optimization of graphs with varying topologies across modalities can lead to the robust detection of shared features.

4.6 Methodological considerations

Several methodological considerations are pertinent to this work. First, although the multi-scale community detection algorithm identifies the entire hierarchical community structure simultaneously, the resolution at which we study the topology of the graph is limited by the number of layers in the multi-scale graphs. This means that for very large graphs, high resolution multi-scale community detection can be computationally intensive. This property can be especially limiting in multi-modal graphs where two graphs show widely different hierarchical organizations, thus making it time-intensive to accurately sample the modularity landscape. Second, one of the greatest advantages of using conventional community detection algorithms is that they provide a single partition of nodes into communities. Using this partition we can extract additional information regarding the role that individual nodes play within the graph by calculating summary statistics including the participation coefficient and within-module degree $z$-score (Guimera and Amaral, 2005). However, there are currently no equivalent summary statistics for multi-scale community organization, and therefore the utility of this approach could be significantly enhanced by the parallel development of such statistics. Third, the modularity quality functions that we studied in this work incorporate a uniform null model (Bassett et al., 2013a; Wymbs et al., 2012; Bassett et al., 2015a; Papadopoulos et al., 2016) instead of the more traditional Newman-Girvan null model (Newman, 2010). We chose this null to preserve the one-to-one relationship between the stability of the local communities and the absolute edge values. We observed that the Newman-Girvan null model changes this relationship, thereby altering the relative estimates of local stability. Future work should explore different null models for multi-scale and multi-modal communities (Paul and Chen, 2016), with the goal of determining their relative utility in extracting hierarchical community structure from brain graphs. Finally, it is worth noting that due to computational limitations, we only examined low values of the inter-layer (scale) dependence. Future work could extend our observations by assessing hierarchical community structure apparent at different values of the inter-layer weight, or when linking layers either ordinally or categorically (Mucha et al., 2010).

5 Conclusion

In this work, we examined a multi-scale community detection algorithm and its advantages for uncovering the hierarchical organization of synthetic and real world graphs. By assuming dependence between the adjacent topological scales, the multi-scale algorithm links the communities persisting over several scales, thereby effectively uncovering hierarchical community organization in graphs. We demonstrated the statistical robustness of this hierarchical organization by defining notions of community stability and inter-scale reliability. After exercising the method on synthetic graphs, we next examined and compared the hierarchical community organization of structural brain graphs and functional brain graphs estimated from diffusion imaging and resting state functional MRI, respectively. Compared to the functional graphs, the structural graphs displayed a higher degree of topological heterogeneity with a more pronounced hierarchical organization as evidenced by a higher average number of stable communities across topological scales. With the exception of basal ganglia-thalamo-cortical circuitry, the structural communities across topological scales tended to be spatially localized, where nodes within a community were located in close physical proximity to one another. Interestingly, functional communities displayed weak similarity to structural communities at coarse topological scales, and this dissimilarity became more pronounced at finer topological scales as spatially distributed functional communities emerged. These statistical differences between the spatially distributed functional communities and spatially localized structural communities were also apparent in an explicit multi-modal extention of our method, which performs a joint optimization of modularity across a multiplex network composed of both functional and structural layers. Taken together, these results illustrate the practical utility of multi-scale community detection in revealing hierarchical community structure in single-modality and multi-modality brain graphs.

6 Acknowledgments

This work was supported by the Army Research Laboratory through contract number W911NF-10-2-0022. DSB would also like to acknowledge support from the John D. and Catherine T. MacArthur Foundation, the Alfred P. Sloan Foundation, the Army Research Office through contract number W911NF-14-1-0679, the National Institute of Health (2-R01-DC-009209-11, 1R01HD086888-01, R01-MH107235, R01-MH107703, R01MH109520, 1R01NS099348 and R21-M MH-106799), the Office of Naval Research, and the National Science Foundation (BCS-1441502, CAREER PHY-1554488, BCS-1631550, and CNS-1626008).The content is solely the responsibility of the authors and does not necessarily represent the official views of any of the funding agencies.

Appendix A Spectral community retection of rynamic (multi-slice) networks

The single layer modularity quality function has been generalized to multi-slice networks to identify communities in multiplex or time-dependent networks. Formally, the multi-slice modularity quality function can be defined as

[TABLE]

where the adjacency matrix of layer $l$ has components $A_{ijl}$, the element $P_{ijl}$ gives the components of the corresponding layer-$l$ matrix for the null model, $\gamma_{l}$ is the structural resolution parameter of layer $l$, $g_{il}$ gives the community assignment of node $i$ in layer $l$, $g_{jr}$ gives the community assignment of node $j$ in layer $r$, $\omega_{jlr}$ gives the connection strength (i.e., an inter-layer coupling parameter) from node $j$ in layer $r$ to node $j$ in layer $l$, the total edge weight in the network is $\mu=\frac{1}{2}\sum_{jr}K_{jr}$, the strength (i.e., weighted degree) of node $j$ in layer $l$ is $K_{jl}=k_{jl}+C_{jl}$, the intra-layer strength of node $j$ in layer $l$ is $k_{jl}=\sum_{i}A_{ijl}$, and the inter-layer strength of node $j$ in layer $l$ is $k_{jl}=\sum_{r}\omega_{jlr}$.

Here we extend the multi-scale framework to multi-slice networks. Formally the multi-slice mulit-scale modularity quality function we study can be defined as

[TABLE]

where the adjacency matrix of layer $l$ has components $A_{ijl}$, the element $P_{ijl}$ gives the components of the corresponding layer-$l$ matrix for the null model, $\gamma_{l}$ is the structural resolution parameter of layer $l$, $g_{il}$ gives the community assignment of node $i$ in layer $l$, $g_{jr}$ gives the community assignment of node $j$ in layer $r$, $\omega_{jlrx}$ gives the connection strength from node $j$ in layer $r$ to node $j$ in layer $l$ at scale layer $x$, $\tau_{jxy}$ gives the topological scale coupling parameter which indicates the strength of the links between neighboring topological scales (as represented by layers), from node $j$ in scale layer $x$ to node $j$ in scale layer $y$. the total edge weight in the network is $\mu=\frac{1}{2}\sum_{jrx}K_{jrx}$, the strength (i.e., weighted degree) of node $j$ in layer $l$ and scale layer $y$ is $K_{jly}=k_{jly}+C_{jly}+T_{jly}$, the intra-layer strength of node $j$ in layer $l$ and scale layer $y$ is $k_{jly}=\sum_{i}A_{ijl}$, and the inter-layer strength of node $j$ in layer $l$ and scale layer $y$ is $C_{jly}=\sum_{r}\omega_{jlry}$, and the inter-scale strength of node $j$ in layer $l$ and scale layer $y$ is $T_{jly}=\sum_{x}\tau_{jlxy}$.

Finally we provide a synthetic example in Fig. 1 to show how the multi-scale community detection algorithm links communities across temporal scales and to uncover the relationships between them.

Appendix B Hierarchical community organization of synthetic graphs

Here we provide synthetic examples of graphs where each node can be identified locally within a community at three different topological scales. Next we create variations in the hierarchical organization of the graphs by systematically introducing edge strength inhomogeneities. As seen in Fig. 1, multi-scale communities and the relative stability of communities across scales clearly uncovers the planted relationship (as well as the inhomogeneity profile) across the nodes.

Appendix C Principal components analysis of the SC and FC stability matrices

Here we use principal components analysis (PCA) to assess the stability profiles of nodes across $\gamma$ increments and measure the topological heterogeneity of the graphs. We used the number of components that account for more than $95\%$ of the variance in the stability matrices as a proxy for topological heterogeneity. Low numbers indicate that most nodes display similar stability profiles and therefore the graph is relatively topologically homogeneous. Conversely, higher numbers indicate that most nodes display diverse stability profiles and therefore the graph is relatively topologically heterogeneous (Fig. 1).

Appendix D Hierarchical community organization of a multiplex SC and FC graph

The FC and SC communities share similar community organization, and joint-optimization of FC and SC graphs (i.e. SC-FC multiplex graphs) can in theory be used to evaluate these similarities. That said, the community organization of the SC-FC multiplex graph is highly dependent on the inter-modality coupling parameter, $\kappa$. In Fig. 1, we demonstrate that at smaller $\kappa$ values the FC-SC graph yields two separate community structures for FC and SC components of the graph; however, for higher $\kappa$ values they both share the same hybrid hierarchical community structure. The direct comparison between the community allegiance matrices of the FC, SC, and SC-FC graphs provided in Fig. 2 shows the effect of increasing inter-modality coupling parameter on the SC-FC community structure.

Appendix E Sorting nodes based on multi-scale community allegiance

A note on visualization. Sorting the nodes of the adjacency matrices based on their community allegiance allows us to visualize communities of densely connected nodes. Nevertheless for the multi-scale communities, the order of the nodes can change depending on the topological scale. In an effort to by-pass this limitation and enhance the visualization of these communities we sort the nodes based on the similarity of their community assignments across scales. Specifically, we perform optimal leaf ordering (optimalleaforder.m) for hierarchical clustering (linkage.m) using the distances (pdist.m) calculated between the community assignments of each pair of nodes. All multi-scale community plots in this manuscript were generated using this node sorting algorithm.

Appendix F Multi-scale group consensus communities in structural and functional brain graphs

Here we provide the group consensus multi-scale community results for the structural (Fig. 1) and functional graphs (Fig. 3). One salient feature of the group consensus SC multi-scale community is that the communities are overwhelmingly made up of neighboring brain structures across the entire range of topological scales. To highlight the spatial proximity of the communities of the structural connectivity graphs, we identified (bwconncomp.m) and only presented the communities with more than one cluster in Fig. 2 while removing all the other communities with only one cluster of brain regions. Next, we tested the statistical significance of these observations via permutation test ($N=1000$) across $\gamma$ increments Fig. 5. Our results demonstrate that the number of SC communities with more than one cluster is significantly ($p<0.01$, Bonferroni corrected for multiple comparisons) smaller than that of the null distribution (generated by changing the assignment of nodes to communities uniformly at random) for all the increments of $\gamma$ (except $\gamma=3$). Unlike structural graphs, functional graphs fail to yield group level consensus results for a large number brain regions, including several subcortical and cortical structures (Fig. 3). We highlighted these structures separately in Fig. 4.

Appendix G Distribution of edge weights in the structural and functional brain graphs

The distribution of the edges in the structural and functional connectivity matrices are notably different. While the distribution of SC edges are extremely heavy-tailed, the FC edges are close to a Gaussian distribution (Fig. 1).