The modular organization of human anatomical brain networks: Accounting for the cost of wiring

Brain networks are expected to be modular. However, existing techniques for estimating a network’s modules make it difficult to assess the influence of organizational principles such as wiring cost reduction on the detected modules. Here we present a modification of an existing module detection algorithm that allowed us to focus on connections that are unexpected under a cost-reduction wiring rule and to identify modules from among these connections. We applied this technique to anatomical brain networks and showed that the modules we detected differ from those detected using the standard technique. We demonstrated that these novel modules are spatially distributed, exhibit unique functional fingerprints, and overlap considerably with rich clubs, giving rise to an alternative and complementary interpretation of the functional roles of specific brain regions. Finally, we demonstrated that, using the modified module detection approach, we can detect modules in a developmental dataset that track normative patterns of maturation. Collectively, these findings support the hypothesis that brain networks are composed of modules and provide additional insight into the function of those modules.


Supplementary Materials and Methods
In the main text we presented an analysis of the modular organization of human structural brain networks.Specifically, we proposed a modification of the well-known modularity quality function (1), wherein we replace the null connectivity model with one that depends upon the brain's spatial embedding.In this supplement we present a number of additional analyses that demonstrate the robustness of our results.

Including versus excluding subcortical regions
In the main text we focused on a network composed of N = 1014 brain regions.These regions were based on a subdivision of the so-called Desikan-Killany atlas (2).That atlas consists of 68 cortical regions and 14 subcortical regions.The subdivision was constructed by dividing the 68 cortical regions into 1000 approximately equal volume regions (3) while not sub-dividing subcortical regions at all.As a consequence, the 1000 cortical regions tend to be of smaller volume than the 14 subcortical regions.In general, regions with larger volume will tend to have higher degree which can bias the topology of the resulting network.To avoid these biases a number of recent studies using the same parcellation have opted to exclude all sub-cortical regions from analysis (4)(5)(6)(7).Here, we demonstrate that the optimal consensus modules we Figure S1: Comparison of communities with and without subcortical regions.To test whether the inclusion of subcortical regions influenced detected modules, we constructed a group-representative network comprised of N = 1000 cortical regions and performed modularity maximization using the SPTL model on this network.(A) We found that the optimal cortical partition was detected at γ = 2.8.Comparing the optimal cortical and subcortical + cortical partitions, we observed that they were highly similar to one another (z rand = 260.26;maximum value of null distribution z rand = 3.40, 10000 permutations).(B) Subcortical + cortical parcellation association matrix ordered by optimal subcortical consensus partition.(C) Cortical parcellation association matrix ordered by optimal subcortical consensus partition.describe in the main text are relatively robust to our decision of whether to include or exclude subcortical regions in the network.
In parallel to our analysis presented in the main text, we generated and analyzed a grouprepresentative connectivity matrix comprised of only the N = 1000 cortical regions.To differentiate this matrix from the one in the main text, we refer to the network analyzed in the main text as SUBCORTICAL + CORTICAL and the cortical only network as CORTICAL.Our analysis of the CORTICAL network was performed in precisely the same manner as our analysis of the SUBCORTICAL + CORTICAL network -e.g.we selected the optimal resolution parameter value based on the median partition similarity and we constructed a consensus partition at that parameter value.For the CORTICAL network we observed an optimal resolution parameter of γ = 2.8 (the SUBCORTICAL + CORTICAL network described in the main text achieved its optimal resolution parameter at γ = 2.6).To assess the similarity of the CORTI-CAL and SUBCORTICAL + CORTICAL consensus partitions, we compared them using the z-score of the Rand coefficient (8).Specifically, we removed the sub-cortical regions from the SUBCORTICAL + CORTICAL partition so that both partitions contained the same 1000 nodes, calculated their similarity, and compared this value to a null distribution obtained by randomly and uniformly permuting the module assignments of the CORTICAL partition 10000 times.The observed similarity of the two partitions was z rand = 260.26while the largest similarity in the null distribution was z rand = 3.40 (Fig. S1A).This suggests that the two partitions are similar to each other well beyond what would be expected under this specific null model.
We support this statistical analysis qualitatively by showing both CORTICAL and SUBCOR-TICAL + CORTICAL association matrices side by side with their rows and columns ordered such that nodes assigned to the same module in the SUBCORTICAL + CORTICAL partition were positioned next to each other (Fig. S1B,C).

Robustness to choice of resolution parameter
Multi-scale modularity maximization involves modifying the modularity quality function by including the resolution parameter, γ (9).This modification makes it possible to tune γ to different values and uncover communities of different sizes, thereby partially mitigating the effect of the so-called "resolution limit" (10).Unfortunately, there is no agreed upon method for selecting the optimal value of γ; most studies avoid the issue by setting its value to its default setting of γ = 1 or by using some heuristic rule for choosing its value.In the main text we selected γ to be the value at which the detected partitions were most similar to one another (as measured by the z-score of the Rand coefficient).This resulted in us focusing on communities detected at γ = 2.6.

It is possible that different heuristics could highlight different resolution parameters. While
it is beyond the scope of this paper to systematically compare such heuristics, we wish to demonstrate that our consensus modules are robust to reasonable variations in γ.Accordingly, we compared consensus partitions obtained at γ = 2.2, . . ., 3.0 in increments of 0.1 (nine total values) with those described in the main text.As in the previous section, our comparison involved calculating the z-score of the Rand coefficent and comparing it against a null distribution obtained through permutation testing.In general, we found that all of the consensus partitions detected in this range were much more similar to those detected at γ = 2.6 than would be expected under the null model (the minimum observed z-score over any of the nine γ values was 222.66 while the maximum value obtained in any of the null distributions was 5.43) (Fig. S2).
In addition to this statistical analysis, we show the association matrices calculated from the consensus modules detected over this same range.To facilitate comparision, we order their rows and columns according to the consensus partition obtained at γ = 2.6 (described in the main text).We see that, qualitatively, the partitions are excellent matches and that most non-zero elements are located within modules (the blocks along the diagonal) (Fig. S3A).We also show, for the same values of γ, the association matrices calculated from the "raw partitions" (i.e.those obtained from modularity maximization but before consensus clustering).(Fig. S3B).

Robustness to variation in max curvature angle
Reconstructing connectivity matrices from diffusion imaging data with tractography algorithms is challenging and prone to both false positives and false negatives (11,12).In addition, tractography algorithms include multiple parameters that can influence the resulting matrices.One important parameter is the max curvature angle that determines the largest orientation change that a streamline can exhibit between integration steps before it is terminated.While it is beyond the scope of the present study to systematically test all possible values of this parameter (and others), we nonetheless wished to demonstrate that our results are consistent if we vary its value within some reasonable range.We calculated the similarity (z-score of the Rand coefficient) of each consensus partition with respect to the partition obtained γ = 2.6 and then generated a null distribution by calculating a comparable similarity score but where the community labels in each consensus partition were uniformly permuted at random.The subject-level networks described in the main text were reconstructed using a max curvature angle of 35 • .Accordingly, we repeated our analysis using networks reconstructed with max curvature angles of 30 • and 40 • .From these networks, we followed the methods described in the main text to construct group-representative networks.As, perhaps, was expected, the group-representative networks varied in terms of gross network statistics.For example, the total number of connections increased monotonically with maximum curvature angle: 51994, 80513, and 119204 connections for the 30  Our analysis of these networks was identical to what was described in the main text.Specifically, we maximized the SPTL modularity over a range of resolution parameters and identified the optimal parameter as the one that maximized the mean pairwise similarity of partitions detected at that value.We found that for the 30 • and 40 • cutoffs, the optimal resolution parameters were γ = 3.0 and γ = 2.1, respectively.Once the optimal parameter was selected, we used the consensus clustering procedure to generate consensus modules.
We then calculated the similarity of these new consensus modules with the consensus modules described in the main text (where similarity was measuresd as the z-score of the Rand coefficient).To contextualize these scores, we compared them against a null distribution obtained by randomly and uniformly permuting module assignments.The observed similarity of the consensus modules obtained with maximum curvature angles of 30 • and 40 • to the con- sensus modules described in the main text were 181.76 and 159.77, respectively, which far exceeded the maximum values obtained in either null distribution (4.48 and 4.28) (Fig. S4A,C).
Based on these observations, we conclude that the consensus modules obtained using 30 • and 40 • cutoffs are, at least in this specific statistical sense, similar to those described in the main text.For qualitative assurance that the modules were similar to one another, we also generated association matrices for the consensus modules obtained using 30 • and 40 • cutoffs and ordered their rows and columns according to the module assignments described in the main text (obtained using 35 • cutoff).Our expectation was that most non-zero elements would fall within modules, whereas few non-zero elements would appear between modules.Indeed, this is the case (Fig. S4B,D).

Participation coefficients of structural versus functional partitions
In the main text we described brain regions' participation coefficients with respect to structural partitions obtained from structural connectivity networks.Specifically, we compared structural partitions obtained by maximizing the standard Newman-Girvan modularity with those obtained using a novel space-dependent modularity.This approach -hub classificiation based on structural partitions -has been described in a number of high profile studies and we sought to continue this tradition (13,14).However, a number of other studies have calculated participation coefficients with respect to a functional partition obtained, for example, from an ICA analysis of fMRI BOLD time series (15) or from community structure analysis of functional connectivity brain networks (16).This second approach makes it possible to assess indirectly the roles of individual nodes with respect to the brain's functional systems.
Therefore, in the interest of completeness, we compared participation coefficients calculated using the consensus partition obtained from the SPTL model with those obtained using the functional partition described in the main text (taken from (6)).We followed precisely the same methods as described in the main text: (1) we ranked participation coefficients, (2) subtracted one set of coefficients from the other, (3) and grouped these coefficients by functional system to estimate the mean change difference in participation coefficient for each system (the statistical significance of which we assessed using a permutation test).We also produced surface plots to show the spatial distribution of participation coefficients obtained using the functional partition.
Though the results of this procedure are not in precise agreement with those described in the main text, they are qualitatively quite similar, serving to highlight the same changes in participation.Notably, both comparisons highlight the somatomotor network as exhibiting increased participation.

Figure S2 :
FigureS2: Statistical assessment of community robustness to variation in γ.We tested how similar consensus partitions from γ = 2.2 to γ = 3.0 were to the partition obtained at γ = 2.6.We calculated the similarity (z-score of the Rand coefficient) of each consensus partition with respect to the partition obtained γ = 2.6 and then generated a null distribution by calculating a comparable similarity score but where the community labels in each consensus partition were uniformly permuted at random.

Figure S3 :
Figure S3: Qualitative assessment of community robustness to variation in γ.As a qualitative assessment of similarity, we show association matrices over the same range of γ constructed from consensus partitions and the ensemble of partitions obtained from the initial modularity maximization.(A) Association matrices constructed from consensus partitions.(B) Association matrices constructed from the initial modularity maximization.

Figure S4 :
Figure S4: Comparison of consensus modules obtained using different maximum curvature cutoffs.(A,C) The z-scores of Rand coefficents of consensus modules described in the main text with consensus modules obtained using different tractography parameters -specifically maximum curvature angles of 30 • and 40 • .We compare a null distribution (blue) versus what was observed (red), noting that in both cases the observed value exceeds the maximum values of the null distribution by a wide margin.(B,D) Association matrices for consensus modules obtained using 30 • and 40 • maximum curvature angles.The rows and columns of each matrix have been ordered according to the consensus modules described in the main text.Note that most non-zero elements tend to fall within modules, indicating high similarity.

Figure S5 :
Figure S5: Comparison of participation coefficients.(A) Ranked participation coefficients obtained using a functional connectivity (FC) partition (from (6)).(B) Difference in ranked participation coefficient obtained from the SPTL model and the FC partition.(C) Regional participation coefficients grouped by functional system.

Figure S6 :
Figure S6: Correlation of fiber length and Euclidean distance We show a strong correlation of a connection's Euclidean distance (straight line distance between its endpoints) and its fiber length (curvilinear trajectory through space).

Figure S8 :
Figure S8: All detected consensus SPTL modules for DSI dataset Topographic representation of the 31 consensus modules detected by applying modularity maximization using the SPTL null model to the human DSI datset.

Figure S9 :
Figure S9: NG modules in Human DSI (A) Distribution of the z-score of the Rand indices as a function of γ. (B) Association matrices (fraction of times out of 500 Louvain runs that each pair of nodes were assigned to the same module) clustered according to consensus modules for the two peaks in the curve shown in panel A, corresponding to γ = 1.0 (left) and γ = 2.1 (right).(C) Module assignments on the cortical surface for modules detected at γ = 1.0 (left) and γ = 2.1 (right).