The impact of graph construction scheme and community detection algorithm on the repeatability of community and hub identification in structural brain networks.

A critical question in network neuroscience is how nodes cluster together to form communities, to form the mesoscale organisation of the brain. Various algorithms have been proposed for identifying such communities, each identifying different communities within the same network. Here, (using test-retest data from the Human Connectome Project), the repeatability of thirty-three community detection algorithms, each paired with seven different graph construction schemes were assessed. Repeatability of community partition depended heavily on both the community detection algorithm and graph construction scheme. Hard community detection algorithms (in which each node is assigned to only one community) outperformed soft ones (in which each node can belong to more than one community). The highest repeatability was observed for the fast multi-scale community detection algorithm paired with a graph construction scheme that combines nine white matter metrics. This pair also gave the highest similarity between representative group community affiliation and individual community affiliation. Connector hubs had higher repeatability than provincial hubs. Our results provide a workflow for repeatable identification of structural brain networks communities, based on the optimal pairing of community detection algorithm and graph construction scheme.

INTRODUCTION

The human brain can be modelled as a network (Bassett & Sporns, 2017) and summarised as a graph. In structural networks, the nodes of the graph are small volumes of tissue which are interconnected via white matter tracts (edges). Graph theory can provide novel insights into healthy human brain function (Bassett et al., 2011; Braun et al., 2015) and its alteration in various diseases (Aerts, Fias, & Caeyenberghs, 2016; Baker et al., 2015; Braun et al., 2016; Collin et al., 2016; Drakesmith et al., 2015; Imms et al., 2019; Nelson, Bassett, Camchong, Bullmore, & Lim, 2017; Vidaurre et al., 2018).

An open question in network neuroscience is how neural units cluster together to form inter‐connected groups and provide the coordinated brain activity that gives rise to action, perception and behaviour (Park, & Friston, 2013; Bassett & Mattar, 2017; Pessoa, 2018). Modularity is a quintessential concept in network neuroscience, wherein neural units are densely connected to one another, forming clusters or modules (Meunier, Lambiotte, & Bullmore, 2010). This is an efficient architecture allowing a complex network to integrate information locally, while maintaining its adaptability to any external stimulus. Networks in nature often show hierarchical, modular organisation (Blondel, Guillaume, Lambiotte, & Lefebvre, 2008; Fortunato, 2010; Fortunato & Castellano, 2012; Lancichinetti & Fortunato, 2009; Lancichinetti, Fortunato, & Kertesz, 2009; Newman & Girvan, 2004 Newman, 2012; Meunier, Lambiotte, Fornito, Karen, & Bullmore, 2009). In the brain, such hierarchical modularity could support segregated neuronal interactions and their integration at the global level. Networks with such structure (Fortunato, 2010) are more complex than those with random structure (Sporns, Tononi, & Edelman, 2000), and have been well demonstrated in functional brain networks (Sporns, 2013; Sporns & Betzel, 2016).

Commonly studied graph‐theoretical metrics, such as global and local efficiency, clustering coefficient, shortest path length, and small‐worldness (Rubinov & Sporns, 2010), provide useful information related to the global and local properties of the network, but to investigate the mesoscale network organisation, community (or modularity) detection techniques are more appropriate (Giusti, Ghrist, & Bassett, 2016; Meunier et al., 2010). Following community partition, a frequently used methodology to identify the community structure uses two modular network metrics (Guimera & Amaral, 2005; van den Heuvel & Sporns, 2013). The participation coefficient, Pi, of a node quantifies the distribution of its links among the modules of the network, while the within‐module z‐score, zi, measures the connectedness of a node to other nodes in the module.

Complex network analysis of multimodal structural and functional brain connectivity has identified a subset of brain areas that play a key role for efficient neural signalling and communication (van den Heuvel & Sporns, 2013). These brain areas, called hubs, support dynamic functional coupling within and between functional subnetworks (Friston, 2011). In empirical structural brain networks, the term ‘rich‐club’ characterises brain areas/nodes with high degree that are more densely interconnected between each other compared to the rest of the network. The metrics Pi and zi can be used to classify hubs as either provincial or connector hubs (Guimera &Amaral, 2005). Please refer to Appendix A for more detail on the computation of Pi and zi and how this classification is made.

Various methodologies for structural network generation have been proposed, involving tractography with different algorithms and assigning edge weights using different diffusion MRI‐based metrics. The resulting graphs are quite different from each other and have different levels of robustness and repeatability (Dimitriadis et al., 2017; Owen et al., 2013; Smith, Tournier, Calamante, & Connelly, 2015; Yuan et al., 2019; Zhong, He, & Gong, 2015). We recently explored the repeatability of structural brain graphs, their edge weights and graph‐theoretical metrics, for 21 different edge‐weighting schemes (Messaritaki, Dimitriadis, & Jones, 2019a). We demonstrated that integrating several metrics as edge weights is very good at capturing differences between populations, and is interesting from the perspective of developing biomarkers (Clarke, Messaritaki, Dimitriadis, & Metzler‐Baddeley, 2021; Dimitriadis et al., 2017).

We constructed structural brain networks from a set of test–retest diffusion MRI scan data from the Human Connectome Project (HCP) using the b = 2000 s/mm2 data and the seven most reproducible graph‐construction schemes as derived from our previous study on the same data (Messaritaki et al., 2019a). We then applied thirty‐three community detection algorithms. The “hard” algorithms assign every node to only one community, while the “soft” algorithms can assign a node to multiple communities. For every pair of community detection algorithm and graph construction scheme, we estimated the reproducibility of nodal Pi and zi and of provincial and connector hubs, based on both modular network metrics. Our aim was to identify the combination of graph construction scheme and community detection algorithm with the highest agreement of individual communities between the two repeat scan sessions.

The quality criterion for the estimated community partitions was also important in our study. To this end, we compared the quality index of the community partitions estimated over the original graphs with the quality indices of the community partitions computed over surrogate null versions of the original graph (Guimerà, Sales‐Pardo, & Amaral, 2004). We previously reported a statistical procedure for performing condition and group comparisons in terms of brain communities (Dimitriadis et al., 2012). Here, we applied a similar approach to assess between‐scan pairwise community similarity for every pair of graph construction schemes and community detection algorithms. We adopted a proper community partition distance metric, the Normalised Mutual Information (NMI) (Lancichinetti & Fortunato, 2009; Alexander‐Bloch et al. 2012). Finally, we derived a consensus cluster across participants and repeat scans (Dong, Frossard, Vandergheynst, & Nefedov, 2014; Ozdemir, Bolanos, Bernat, & Aviyente, 2015). The agreement of consensus cluster with individual communities adopting NMI was also used as an objective criterion of the optimal combination of graph construction scheme and community detection algorithm.

We note that the analysis presented here does not aim to assess how well these structural networks represent the functional organisation of the human brain. The accuracy of these networks and of the metrics used as edge‐weights in representing the functional organisation of the brain has been validated in recent work by Messaritaki et al. (2021). Additionally, the metrics used as edge‐weights are routinely used in network analyses in the literature (e.g., Caeyenberghs, Metzler‐Baddeley, Foley, & Jones, 2016; Nigro et al., 2016; Taylor, Cheol, Han, Weber, & Kaiser, 2015). Our analysis does, however, address one aspect of the accuracy of the partition of the structural connectome. If a partition of the structural connectome is not repeatable in the absence of changes resulting from maturation or intervention, then that partition is not an accurate representation of the modular organisation of the structural connectome. Only partitions that are repeatable can convey reliable information about the structural organisation of the human brain. In other words, even though the repeatability of a partition is not a sufficient condition for it to be representative of the brain's structural organisation, it is a necessary one.

The rest of this manuscript is organised as follows: Section 2 describes the graph‐construction schemes, community detection algorithms, community partition similarity, the methodology for detecting connector and provincial hubs, and their repeatability. Section 3 reports our results in terms of repeatable community partitions across the 2D space of graph‐construction schemes/community detection algorithms, the repeatability of nodal Pi/ zi and the detection of connector/provincial hubs. The Discussion summarises the main outcome of our study explaining its advantages, limitations, and suggestions for future directions.

METHODS

All analyses were performed using MATLAB (2019a; The Mathworks, Inc., MA).

Data

We analysed the test–retest MRI and diffusion‐MRI data set from the multimodal neuroimaging database of the Human Connectome Project (HCP) (Glasser et al., 2013; Sotiropoulos et al., 2013; van Essen et al., 2013). We used the data from the 37 participants for whom there were 90 gradient directions for each b‐value. The participants on this test–retest data set were scanned twice with the between‐scan time interval ranging between 1.5 and 11 months. The age‐range of the participants was 22–41 years. The test–retest time interval is shorter than the expected time over which maturation‐induced structural changes can be measured with diffusion MRI (dMRI).

The diffusion‐weighted images (DWIs) had a resolution of (1.25 × 1.25 × 1.25) mm3 and were acquired for three different diffusion weightings (b‐values: 1000, 2000 and 3,000 s/mm2) across 90 gradient directions for each b‐value. The HCP acquisition details and pre‐processing are described in Sotiropoulos et al. (2013) and Sotiropoulos et al. (2013), Feinberg et al. (2010), Moeller et al. (2010), Setsompop et al. (2012), Xu et al. (2012), Glasser et al. (2013). Specifically, the diffusion images were corrected for EPI distortions, eddy‐current distortions, participant movement and gradient nonlinearities. The diffusion data were also registered to the structural data. We performed the following analyses using the b = 2000 s/mm2 data.

Tractography

We performed whole‐brain tractography using ExploreDTI‐4.8.6 (Leemans, Jeurissen, Sijbers, & Jones, 2009), estimating the fibre orientation distribution function (fODF) using constrained spherical deconvolution (CSD) (Tournier et al., 2004). Tracking was initiated on a 2 × 2 × 2 mm grid, with a 1 mm step size, angular threshold of 30∘ and fibre length range of 50–500 mm.

Graph generation

Node definition

The Automated Anatomical Labelling (AAL) atlas (TzourioMazoyer et al., 2002) was used to define 90 cortical and subcortical areas (45 areas per hemisphere) as nodes of the structural brain graphs. Structural brain networks (SBN) were generated for each participant using Explore DTI‐4.8.6 (Leemans et al., 2009).

Edge weights

Edges were weighted using the seven most reproducible graph‐construction schemes identified previously with the same data set (Messaritaki et al., 2019a), and which were based on different combinations of the nine metrics listed in Table 1 (see Section 2.3.4). Each graph was normalised to have a maximum edge weight of 1, while the elements in the main diagonal were set to zero (see Figure 1).

TABLE 1

Metrics used in connectivity matrices

Metric	Abbreviation
Fractional anisotropy	FA
Mean diffusivity	MD
Radial diffusivity	RD
Number of streamlines	NS
Percentage of streamlines	PS
Streamline density	SLD
Tract volume	TV
Tract length	TL
Euclidean distance between nodes	ED

FIGURE 1

Flowchart of the construction of a structural brain network based on tractography and diffusion metrics (see Table 1)

Integrated edge‐weights

Combining multiple metrics into an integrated edge weight is supported by the fact that each metric conveys information about different tissue properties, while at the same time topological properties of SBNs are affected by more than one metric. Here, using the data‐driven algorithm described in our previous work (Dimitriadis et al., 2017,c) the nine metrics in Table 1 were used to form an integrated SBN for each participant and scan session.

The algorithm down‐weights the more similar metrics and up‐weights the most dissimilar metrics to enhance the integration of complementary topological information across the nine metrics. An orthogonal‐minimal‐spanning‐tree (OMST) algorithm was then applied to the resulting networks, selecting the edges that preserve connectivity between nodes, while guaranteeing that the overall network efficiency is maximised. More details on the OMST algorithm and its implementation can be found in our previous work (Dimitriadis et al., 2017,b,c; Messaritaki, Dimitriadis, & Jones, 2019b) and the related code is freely available at https://github.com/stdimitr/multi-group-analysis-OMST-GDD.

Graph construction schemes

Seven graph construction schemes were used in this study, summarised in Table 2 and falling broadly into two categories. We briefly explain their construction methodologies here.

TABLE 2

Summary of the graph‐construction schemes

Abbreviation	Initial edge weights	Topology	Final edge weights	Symbol
NS–OMST	NS	OMST	Unchanged	A
NS + FA OMST	Lin. Comb. of NS and FA	OMST	Unchanged	B
9‐m OMST	Lin. Comb. of all 9 metrics in Table 1	OMST	Unchanged	C
NS‐thr	NS	Keep highest‐NS edges	Unchanged	D
NS‐t/FA‐w	NS	Keep highest‐NS edges	Re‐weight with FA	E
NS‐t/MD‐w	NS	Keep highest‐NS edges	Re‐weight with MD	F
FA‐t/NS‐w	FA	Keep highest‐FA edges	Re‐weight with NS	G

The first category includes graphs constructed via the data‐driven algorithm (Dimitriadis et al., 2017, c). (a) NS‐OMST: apply the OMST filtering algorithm (Dimitriadis et al., 2017,c) to the NS‐weighted matrix. (b) NS + FA−OMST: Integrate the NS‐weighted and FA‐weighted matrices with the data‐driven algorithm. (c) 9 m‐OMST: Integrate all nine diffusion metrics (as originally reported in Dimitriadis, Drakesmith, et al., 2017, see Table 2).

The second category includes SBNs with edges weighted by the NS or the FA and applying a threshold to remove edges with the lowest weights. The threshold was determined by imposing the constraint that the graphs exhibit the same sparsity as the OMST graphs that exhibited the highest reproducibility (Messaritaki et al., 2019a). Once the topology of each of those graphs was specified, the weights of the edges were either kept as they were or re‐weighted with one of the remaining two metrics. These graphs are as follows (see Table 2). (d) NS‐thr: Keep the highest‐NS edges. (e) NS‐t/FA‐w: Threshold to keep the highest‐NS edges, then reweight those edges with their FA. (f) NS‐t/MD‐w: Keep the highest‐NS edges, then reweight those edges with their MD. (g) FA‐t/NS‐w: Keep the highest‐FA edges, then reweight those edges with their NS.

As we have shown previously, these seven schemes exhibit different values of similarity between them, from 0.99 to 0.42 (Messaritaki et al., 2019a, Table 2), motivating their inclusion in a study on the repeatability of community detection.

Community detection algorithms

Communities or modules are defined as subgroups of nodes that are more interconnected with each other compared to the rest of the network (Newman & Girvan, 2004; Radicchi, Castellano, Cecconi, Loreto, & Parisi, 2004). In the present study, we compared thirty‐three different community detection algorithms, comprising twenty‐six with hard clustering and seven with soft clustering. (see Figure 2). In hard clustering, community membership can be represented as a vector that encapsulates the assignments of every brain area to every detected graph cluster (community). In our case, clustering has a dimension of 1 x 90, equalling the number of brain regions in the AAL parcellation. In soft clustering, the outcome is a matrix that encapsulates how many soft clusters a given node (brain area) belongs to. A more detailed description of the adopted community detection algorithms is provided in Appendix B.

FIGURE 2

An example of hard and soft clustering in a toy example containing 7 nodes. (a) Hard clustering: A node can only belong to one cluster. The table shows the community assignment to every node. (b) Soft clustering: Five out of seven nodes are clustered in a single cluster/community {nodes 1,2,3,4,6} while nodes {5 and 7} belong to two communities: node 5 belongs to communities 1 and 3 while node 7 belongs to communities 2 and 3. The table shows the community assignment to every node

In the present study, we considered for the very first time, to the best of our knowledge, a large number of community detection algorithms. We adopted thirty‐three graph partition algorithms further divided into twenty‐six hard clustering algorithms and seven soft clustering algorithms.

Permutation test on quality modular indices

For every participant, scan and graph construction method, we produced 1,000 surrogate null graph models by randomising the weighted connections while preserving both the degree and strength of every node and the overall connectedness of the network (Rubinov & Sporns, 2010).

All of the hard clustering algorithms (no.s 1–26) involved a Q quality index for the communities detected. For further details of Q quality indices see Le Martelot and Hankin (2011, 2012a, 2012b).

For the soft clustering algorithms (no. 27–33), we estimated the normalised mutual information (NMI; see Appendix C) between the original community affiliation and the surrogate null communities produced via the application of every algorithm to the surrogate graph model.

Between‐scan community detection agreement

We quantified the graph‐partition distance with the normalised mutual information (NMI; see Appendix C).

Consensus clustering

A consensus matrix was constructed for every pair of {graph construction scheme—community detection algorithm} that showed high test–retest reproducibility across the cohort (group‐averaged NMI > 0.9). We quantified how many times two nodes across the 74 SBNs (37 participants x 2 scans) were classified in the same community and this entry, ti,j, was assigned to the relevant pair of nodes. The consensus matrix has the same dimensions as the original SBN (90 × 90 in our case) with entries assuming integer values between 0 and 74 {37 participants x 2 scans}, which were then transformed to denote the probability of a pair of nodes (brain areas) being classified as belonging to the same community across the cohort and scan sessions. We converted the consensus matrix into a probability one by dividing each entry by 74.

In order to get a consensus or group representative community per graph‐construction scheme and community detection algorithm, consensus matrices should be iteratively thresholded and clustered with a community detection algorithm (Lancichinetti & Fortunato, 2012). This algorithm uses an absolute arbitrary threshold to eliminate weak connections and iteratively apply a graph partition technique. Instead of an arbitrary filtering scheme, we adopted our OMST algorithm (Dimitriadis et al., 2017,b,c) to topologically‐filter the consensus matrix in a data‐driven way. We then extracted the consensus—group representative community by applying the community detection algorithms across the graph construction schemes (Newman, 2006). See Figure 3e for an example of a consensus matrix.

FIGURE 3

Outline of the presented methodology. The demonstration based on 9‐m OMST graph‐construction scheme and gso‐discrete mode community detection algorithm. (a) Repeat—Scan Sessions. (b) Structural brain networks from participant 1 from both sessions using 9‐m OMST graph‐construction scheme. (c) Individual community affiliation of participant 1 from scan session 1. Each colour represents one community. (d) Vectorised community affiliations of the whole cohort from scan sessions 1 and 2 separated with a red line. Every module is coded with a different colour. (e) Consensus matrix is built over group community affiliations across both sessions as presented in (d). Weights in the consensus matrix refer to the total number of times two brain areas are grouped together across the cohort and scan sessions with the maximum value being (number of participants) x (scan sessions) = 74. (f) Representative community affiliation after graph partitioning the consensus matrix presented in (e). Each community is encoded to a different colour. Similarity NMI distance has been estimated between representative community affiliation presented in (f) and individual community affiliations presented in (c)

Agreement of consensus representative community with individual community structures

An important criterion of our analysis is the high similarity between the consensus clustering and individual clustering for every graph construction scheme that showed high group‐averaged community similarity (NMI > 0.9). To this end, we estimated this community similarity for every case.

Figure 3 illustrates the various steps of the analysis.

Evaluating the combined graph construction schemes—Community detection algorithms

As mentioned previously, we first identified the combinations of graph construction schemes and community detection algorithms with higher group‐averaged between‐scan community affiliation agreement (NMI > 0.9) with a p < .05 based on the bootstrapping procedure. We then adopted a criterion of highest community similarity between the consensus clustering with individual community affiliation (clustering). It is important that consensus clustering expresses the inter‐subject variability and acts as a vector median for the whole group (Dimitriadis et al., 2012). The final ranking of pairs of graph construction schemes and community detection algorithms was based on: (a) high between‐scan group‐averaged community similarity quantified with NMI, and supported by a p < .05 ((bootstrapping, see Appendix C), (b) Q quality index with p < .05 based on surrogate null brain models and (c) high community similarity between consensus clustering and individual community affiliations (clusterings) assessed via a two‐way analysis of variance (ANOVA) (see Section 2.7).

Modular driven structural brain hub detection

Appendix A describes in detail the computation of Pi and zi and how hubs are classified as either provincial or connector hubs (Guimera & Amaral, 2005). Here, we applied the aforementioned hub detection methodology solely on the graph construction schemes and community detection algorithm that fulfil the evaluation criteria of Section 2.3.10.

Reliability of nodal participation coefficient Pi and within‐module Z‐score zi

We also explored the intra‐class correlation coefficient (ICC) of nodal participation coefficient P and within‐module z‐score z . As a main outcome of this hub detection approach, we quantified the consistency of connector/provincial hub detection first within participant between scans, and secondly across the cohort.

Assessing a reproducible structural core of the human brain

We detected structural hubs for every participant, scan session and graph construction scheme by applying an absolute threshold to the participation coefficient and within‐module z‐score (Guimera & Amaral, 2005; Hagmann et al., 2007, 2008); van den Heuvel & Sporns, 2013). We estimated an agreement index that quantifies the percentage of connector/provincial hubs that were detected in both scans. This agreement index is defined as:where CH1,2 are two vectors of size 1 ×  37 (number of subjects) with ones in positions where a brain area is detected as connector or provincial hub in specific subjects. This agreement index is normalised by the total number of participants and takes the absolute value of 1 when a node/ROI is detected as either connector or provincial hub across all participants and in both scans. We characterized an ROI as either provincial or connector hub only if the related Agreement index equals 1.

Statistical analysis

Firstly, we determined pairs of graph construction schemes and community detection algorithms that fulfilled the evaluation criteria presented in Section 2.3.10. The first two criteria were evaluated via bootstrapping and surrogate null models, while for the third, based on the similarity of individual community partitions with consensus community partition, we run a two‐way ANOVA over the pairs of graph construction schemes and community detection algorithms (p < .05). We run also a two‐way ANOVA over the pairs of graph construction schemes and community detection algorithms to assess the effect of both factors and their synergy to the repeatability of network topology quantified with NMI.

The detection of reproducible brain structural hubs using community‐based hub detection network metrics require reproducible communities. For that reason, we followed hub detection analysis over the best pairs of graph construction schemes and community detection algorithms. Then, we adopted a two‐way ANOVA (p < .05) to detect the best pair of graph construction scheme and community detection algorithms over hub detection analysis, using as input the ICC for the participation index and the within‐module z‐score across nodes and the agreement‐index.

RESULTS

Quality of the detected communities

Q original values were transformed into p‐values by comparing them with the 1,000 surrogate (permuted) Q values. p‐values ranged between .013 and .021 across participants, scans and graph construction schemes. For the soft clustering algorithms (no. 27–33) that do not include a Q quality index, we measured the NMI between the original community affiliation and the 1,000 surrogate‐based communities. The group‐averaged NMI ranged between 0.08 and 0.12 which supported the quality of the derived soft communities with the whole set of soft clustering algorithms. These findings support the quality of the extracted graph partitions and allow us to include all the participants, scans, graph construction schemes and community detection algorithms in our analysis.

Group‐averaged between‐scan agreement of communities affiliations

Figure 4 shows group‐averaged between‐scan agreement of community affiliations across graph constructions schemes and community detection algorithms. Based on the highest group‐averaged NMI values and the detected p‐values (p = .0001) derived from the permutation test, we detected the following four community detection algorithms across the seven graph construction schemes as having the highest between‐scan agreement (NMI > 0.9):

FIGURE 4

Between‐scan agreement of communities affiliations across graph construction schemes and community detection algorithms. Every subplot refers to one of the seven graph‐construction schemes. The bars define the group‐averaged between‐scan agreement of community affiliations. Numbers below the plot in A refer to the number list of community detection algorithms represented in Section 2.3.5 and in Appendix B. Community detection algorithms with the highest agreement between the two scans (NM1 > 0.9) were: mscd_afg,mscd_rb, mscd_rn and mscd_so. For the abbreviations and numbering of the community detection algorithms please see Appendix B

(mscd_afg): Fast multi‐scale community detection algorithm using the criterion from Arenas, Fernandez, and Gomez (2008)

(mscd_rb): Fast multi‐scale community detection algorithm using the criterion from Reichardt & Bornholdt, (2006).

(mscd_rn): Fast multi‐scale community detection algorithm using the criterion from Ronhovde and Nussinov (2009)

(mscd_so): Fast multi‐scale community detection algorithm using stability optimisation (Le Martelot & Hankin, 2012a). Two‐way ANOVA reported an effect of the graph construction scheme on the repeatability of network topologies between scans assessed with NMI (df = 6, F = 199.64, p = .124 × 10−24, corrected for multiple comparisons), an effect of the community detection algorithm (df = 32, F = 980.07, p = .345 × 10−49, corrected for multiple comparisons), and an interaction effect (df = 192, F = 28.18, p = .226 × 10−9, corrected for multiple comparisons).

Similarity of individual community partitions with consensus community partition

The highest similarity between individual community partitions and consensus community partition was detected for the combination of 9‐m‐OMST graph construction scheme and mscd_so community detection algorithm. The second highest similarity was detected for mscd_afg and 9‐m OMST. Mscd with rb and rn criterions failed to produce an acceptable community similarity between a group representative community estimated via consensus clustering and the individual community affiliations (see Table 3). Two‐way ANOVA reported an effect of group‐averaged similarity of individual community partitions with consensus community partition across the seven graph construction schemes (df = 6, F =104.96, p = .783 × 10−17, corrected for multiple comparisons), an effect of the four community detection algorithms (df = 3, F = 2481.91, p = .641 × 10−36, corrected for multiple comparisons), and an interaction effect (df = 18, F = 28.94, p = .132 × 10‐9 corrected for multiple comparisons). Based on the mean group‐averaged similarity across the seven graph construction schemes, we ranked the four community‐detection algorithms. The mscd_so community detection algorithm produced the highest group‐averaged similarity of individual community partitions with consensus community partition compared to the four community detection algorithms.

TABLE 3

Group‐averaged similarity of individual community partitions with consensus community partition. Similarities are expressed in NMI scale

	NS‐OMST	NS + FA OMST	9‐m OMST	NS‐thr	NS‐t/FA‐w	NS‐t/MD‐w	FA‐t/NS‐w
mscd_afg	0.64	0.57	0.70	0.61	0.61	0.61	0.39
mscd_rb	0.53	0.39	0.39	0.32	0.31	0.30	0.18
mscd_rn	0	0	0.36	0	0	0	0
mscd_so	0.70	0.65	0.75	0.68	0.68	0.68	0.64

Note: We assigned with bold the top ranked values.

Evaluation of the best combination of graph construction scheme and community detection algorithm

Based on the evaluation criteria presented in Sections 2.3.6, 2.3.7, 2.3.8, 2.3.9 for the final ranking of the pairs of graph construction schemes and community detection algorithms, and taking into account the individual scores, we ranked as the best pair the combination of mscd_so with 9‐m OMST. This pair showed a high between‐scan group‐averaged community similarity estimated with NMI and evaluated via bootstrapping (p = .0001). The Q quality index has a value of 0.67 with p = .001. The mscd_so community detection algorithm also produced the highest group‐averaged similarity of individual community partitions with consensus community partition compared to the four community detection algorithms supported by two‐way ANOVA.

Figure 5 illustrates the topology of the nine communities across the 90 AAL brain areas based on the combination of mscd_so with 9‐m OMST. Importantly, modules number 1, 7 and 9 group together brain areas located within the left hemisphere, modules number 2, 6 and 8 group together brain areas located within the right hemisphere, while modules 3, 4 and 5 involve areas from both hemispheres. Specifically, they involve bilateral ROIs from the fronto‐parietal, cingulo‐opercular and default mode networks like rectus, anterior and middle gyrus, frontal superior gyrus, frontal superior medial gyrus, supplementary motor area, precuneus, cuneus, calcarine and occipital superior gyrus. The bilateral superior temporal gyrus, superior temporal pole and middle temporal pole play an inter‐hemispheric modular connector role (see * in Figure 5). Five out of thirteen consistent connector hubs are located within inter‐hemispheric modules. Interestingly, eight homologous brain areas are grouped together in either left (Module 8) or right hemisphere (Module 9). These areas are: hippocampus, parahippocampal gyrus, amygdala, inferior temporal gyrus, thalamus, pallidum, fusiform and lingual gyrus.

FIGURE 5

Topological Layout of Modular Assignment into the 90 AAL brain areas based on the community affiliation extracting from the consensus matrix related to 9‐m OMST graph construction scheme and mscd_so community detection algorithm. With ‘*’, we denoted the connector hubs detected consistently across participants and repeat scans from the same combination of {mscd‐so, 9‐m OMST} (see Section 3.5). This circular plot illustrates the 90 AAL brain areas into 45 of the left hemisphere on the left semi‐circle and 45 of the right hemisphere on the right semi‐circle. Our analysis gave nine communities/modules where each one is encoded with a different colour

ICC of nodal participation coefficient index and within‐module z‐score

Table 4 shows the group‐averaged ICC of nodal Participation Coefficient, P i. ICC values were first estimated per node and were then averaged across the 90 nodes for every pair of graph construction scheme and community detection algorithms. The highest ICC score was detected for the mscd_so algorithm combined with both the 9‐m OMST and NS‐thr graph construction schemes. On average across the seven graph construction schemes, the mscd_so algorithm also demonstrated the highest average ICC score. Two‐way ANOVA revealed a main effect of network‐averaged ICC of Participation Coefficient index Pi across the four community detection algorithms (df = 3, F = 34.63, p = 1.09 x 10‐7, corrected for multiple comparisons), a main effect of the graph construction scheme (df = 6, F = 45.49, p = .651 × 10−10, corrected for multiple comparisons), and no significant interaction effect. The mscd_so community detection algorithm produced the highest group‐averaged ICC of nodal Participation Coefficient Pi.

TABLE 4

ICC of nodal participation coefficient index Pi across every combination of graph construction scheme with community detection algorithms

	NS‐OMST	NS + FA OMST	9‐m OMST	NS‐thr	NS‐t/FA‐w	NS‐t/MD‐w	FA‐t/NS‐w
mscd_afg	0.54	0.50	0.62	0.67	0.58	0.60	0.38
mscd_rb	0.61	0.55	0.68	0.73	0.67	0.65	0.46
mscd_rn	0.58	0.55	0.68	0.67	0.61	0.61	0.43
mscd_so	0.69	0.70	0.80	0.80	0.77	0.76	0.43

Note: We denote the top ranked values in bold letters.

Figure 6 shows the nodal ICC for the {mscd_so,9‐m OMST} and {mscd_afg, NS‐thr} pairs. Applying a Wilcoxon Rank‐Sum‐test between the two sets of 90 ICCs, we detected a significantly higher ICC for the mscd_so – 9‐m OMST pair (p‐value = .041).

FIGURE 6

ICC of nodal participation coefficient index Pi for the best combinations of graph construction scheme and community detection algorithm. (a)mscd_so – 9‐m OMST. (b) mscd_afg – NS‐thr

Table 5 shows the group‐averaged ICC of nodal within‐module Z‐scores, z i. ICC values were first estimated per node and were then averaged across the 90 brain areas for every pair of graph construction scheme and community detection algorithms. The highest ICC scores were found for the mscd_rb and mscd_so mode algorithms combined with the NS‐thr graph construction scheme. On average across the seven graph construction schemes, the mscd_so mode algorithm demonstrated the highest average ICC score. These findings reached statistical significance. Two‐way ANOVA revealed a main effect of network‐averaged ICC of zi across the four community detection algorithms (df = 3, F = 10.16, p = .00004, corrected for multiple comparisons), a main effect of the graph construction scheme (df = 6, F = 3.9, p = .00113, corrected for multiple comparisons), and no significant interaction effect. The mscd_rb community detection algorithm produced the highest group‐averaged ICC of nodal within‐module Z‐scores, zi.

TABLE 5

ICC of nodal within‐module z‐score across every combination of graph construction scheme with community detection algorithms

	NS‐OMST	NS + FA OMST	9‐m OMST	NS‐thr	NS‐t/FA‐w	NS‐t/MD‐w	FA‐t/NS‐w
mscd_afg	0.71	0.66	0.65	0.79	0.69	0.69	0.69
mscd_rb	0.74	0.66	0.64	0.81	0.72	0.72	0.73
mscd_rn	0.56	0.57	0.60	0.68	0.55	0.57	0.61
mscd_so	0.73	0.59	0.73	0.75	0.60	0.62	0.52

Note: We denote the top ranked values in bold letters.

Figure 7 shows the nodal ICC for the {mscd_rb, NS‐thr} and {mscd_so, NS‐thr}. Applying a Wilcoxon Rank‐Sum‐test between the two sets of 90 ICCs, we found the mscd_rb algorithm had significantly higher ICC (p‐value = .0335 × 10−9).

FIGURE 7

ICC of nodal within‐module z‐score for the best combinations of graph construction scheme and community detection algorithm. (a) mscd_rb – NS‐thr. (b) mscd_so –NS‐thr

Reproducibility of structural hubs detection based on participation coefficient index and within‐module z‐score

We estimated the Agreement index of both connector and provincial hub detection across the cohort. The Agreement index of provincial hub detection across the cohort was higher than the Agreement index for connector hubs (Table 6 versus Table 7).The highest Agreement index for provincial hub detection was found for {mscd_so, 9‐m OMST} (Table 6). On average across the seven graph construction schemes, the mscd_so algorithm demonstrated the highest average Agreement index for provincial hub detection. Two‐way ANOVA revealed a main effect of network‐averaged Agreement index of provincial hubs across the four community detection algorithms (df = 3, F = 54.19, p = .542 × 10−12, corrected for multiple comparisons), a main effect of the graph construction scheme (df = 6, F = 3.1, p = .00291, corrected for multiple comparisons), and no significant interaction effect. The mscd_so community detection algorithm produced the highest Agreement index for provincial hub detection. However, no ROI is detected as provincial hub consistently across the cohort (Agreement < 1).

TABLE 6

Agreement indices of provincial hubs across every combination of graph construction scheme and community detection algorithm

	NS‐OMST	NS + FA OMST	9‐m OMST	NS‐thr	NS‐t/FA‐w	NS‐t/MD‐w	FA‐t/NS‐w
mscd_afg	0.69	0.72	0.75	0.60	0.58	0.51	0.52
mscd_rb	0.46	0.23	0.28	0.30	0.31	0.30	0.21
mscd_rn	0.00	0.00	0.27	0.00	0.00	0.00	0.00
mscd_so	0.69	0.73	0.83	0.80	0.59	0.61	0.24

Note: We assigned with bold the top ranked values.

TABLE 7

Agreement indices of connector hubs across every combination of graph construction scheme and community detection algorithm

	NS‐OMST	NS + FA OMST	9‐m OMST	NS‐thr	NS‐t/FA‐w	NS‐t/MD‐w	FA‐t/NS‐w
mscd_afg	0.40	0.13	0.23	0.38	0.32	0.31	0.43
mscd_rb	0.65	0.42	0.45	0.67	0.57	0.58	0.80
mscd_rn	0.13	0.05	0.09	0.19	0.15	0.14	0.49
mscd_so	0.17	0.03	0.05	0.12	0.12	0.12	0.17

Note: We assigned with bold the top ranked values.

The highest Agreement index for connector hubs detection was found for {mscd_rb, FA‐t/NS‐w} (Table 7). On average across the seven graph construction schemes, the mscd_rb algorithm demonstrated the highest average Agreement index of connector hub detection. These findings are supported statistically. Two‐way ANOVA revealed a main effect of network‐averaged Agreement index of connector hubs across the four community detection algorithms (df = 3, F = 89.75, p = .510 × 10−14 corrected for multiple comparisons), an effect of the graph construction scheme (df = 6, F = 11.77, p = .215 x 10‐5, corrected for multiple comparisons), and no significant interaction effect. The mscd_rb community detection algorithm produced the highest agreement index for connector hub detection. Figure 8 shows the Agreement index for connector hubs for the pair {mscd_rb, FA‐t/NS‐w}. It is clear that 13 ROIs are consistently detected as connector hubs across the cohort (Agreement = 1).

FIGURE 8

Agreement index of connector hubs for the best pair of {mscd_rb, FA‐t/NS‐w}

The group of connector hubs is indicated alongside modular representation of consensus modules illustrated in Figure 5 and also in Table 8. Interestingly, 5 of 13 consistent connector hubs are located within the inter‐hemispheric modules (see Figure 5). Our conclusion is that the combination of modular network metrics Pi and zi succeeded in uncovering a consistent core of connector hubs but failed to detect provincial hubs consistently.

TABLE 8

Consistent connector hubs aligned with the detected module number illustrated in Figure 5

Connector hubs	Module number
FrontSup_L	3
FrontSup_R	3
FrontInfOrb_L	2
RolOper_L	2
Lingual_L	6
Lingual_R	7
OccMid_L	8
OccMid_R	9
Precuneus_R	5
TempSup_R	9
TempPoleSup_L	2
TempMid_L	8
TempInf_R	7

DISCUSSION

We have presented the first extensive study in the literature on the robustness of community detection in structural brain networks by exploring different graph‐construction schemes (previously shown to exhibit high repeatability themselves) and various community detection algorithms. Our main findings have direct implications for longitudinal studies and studies comparing healthy controls versus diseased populations.

The key findings of our analysis can be summarised as follows:

The repeatability of community affiliations depends heavily on the combination of graph‐construction scheme and community detection algorithm. All previously reported studies of network communities adopted a specific pair of graph‐construction scheme and community detection algorithm, with the majority of them focused on Newman's modularity objective criterion (Betzel et al., 2017; Newman, 2006; Sporns & Betzel, 2016). Based on our first criterion of high repeatability of community affiliation between the two scans and across the cohort (NMI > 0.9) supported statistically via bootstrapping (p = .0001), we identified four community detection algorithms as the best choices:

A. mscd_agb.

B. mscd_rb.

C. mscd_rn.

D. mscd_so.

These four algorithms gave excellent repeatability across the entire set of graph‐construction schemes (see Table 2 and Figure 4). Two‐way ANOVA showed an effect of the repeatability of network topologies between scans assessed with NMI across the graph construction schemes, an effect of the community detection algorithm, and an interaction effect of both factors. These findings support our hypothesis that repeatable identification of structural brain networks communities can be derived from the optimal pairing of community detection algorithm and graph construction scheme.

Complementary to the repeatability of the community partitions, we assessed their quality (second criterion). The majority of the graph partition algorithms used here (28 of 33) estimated a quality index in combination with the community partition. By comparing the quality index of the original community partition to that derived from 1,000 surrogate null graph versions, we computed their significance values to be in the range 0.013 ≤  < .021. For the soft clustering algorithms, where this quality index was absent, we estimated the mean clustering distance (mean NMI) between the original partition and partitions derived from 1,000 surrogate null graph versions of the original graph. The mean NMI was below 0.12. Our findings indicated community detection at the mesoscale is repeatable and, as the Q quality for the four selected algorithms reached statistically significant levels, is of high quality.

Our third criterion was the high similarity of consensus community affiliation with individual community affiliations. The majority of the neuroimaging studies that employed a single graph‐construction scheme and one community detection algorithm did not evaluate their findings under this framework (Akiki & Abdallah, 2019; Rasero et al., 2017; Ryali, Chen, Padmanabhan, Cai, & Menon, 2015). The consensus matrix integrates community affiliations across the entire cohort counting the number of times two nodes are assigned to the same community. To reveal a representative community affiliation, we have to apply a community detection algorithm. However, almost all previous studies report this representative community affiliation was determined without estimating its similarity with individual community affiliations. Following a two‐way ANOVA, we observed an effect of group‐averaged similarity of individual community partitions with consensus community partition across the four best community detection algorithms, an effect of the graph construction scheme, and a significant interaction effect of both factors. The mean group‐averaged similarity across the seven graph construction schemes allowed us to rank the community detection algorithms. Of the four community‐detection algorithms, the mscd_so community detection algorithm produced the highest group‐averaged similarity of individual community partitions with consensus community partition. We found the highest similarity between individual community partitions and representative community partition for the {mscd_so, 9‐m OMST} pair. The second highest similarity was detected for {mscd_afg, 9‐m OMST} pair.

An important result of our analysis is that soft clustering community detection algorithms gave the least repeatable results. Therefore, we recommend the use of hard‐clustering algorithms for the detection of brain communities, at least when using the AAL template.

The best combination of graph‐construction scheme (9‐m OMST) and community detection algorithm (mscd_so) revealed nine distinct modules (as illustrated topologically in Figure 5) with interesting findings:

A. Modules 1, 7 and 9 group together brain areas located exclusively within the left hemisphere while modules 2, 6 and 8 group together brain areas located exclusively within the right hemisphere. Modules 3, 4 and 5 group brain areas from both hemispheres together.

B. Modules 3–5 that integrate brain areas from both hemispheres involve bilateral ROIs from fronto‐parietal, cingulo‐opercular and default mode network like rectus, anterior and middle gyrus, frontal superior gyrus, frontal superior medial gyrus, supplementary motor area, precuneus, cuneus, calcarine and occipital superior gyrus. A multi‐tasking fMRI study has previously suggested the flexible role of the fronto‐parietal network in cognitive control and adaptive demands of cognitive tasks (Cole et al., 2013).

C. Five out of thirteen consistent connector hubs are located within inter‐hemispheric modules 3–5 supporting their inter‐connecting role (Figure 5).

D. Interestingly, eight homologous brain areas were grouped together in either left (module 8) or right hemisphere (module 9). These areas are: hippocampus, parahippocampal gyrus, amygdala, inferior temporal gyrus, thalamus, pallidum, fusiform and lingual. Lesions of hippocampus, parahippocampal gyrus, amygdala and fusiform gyrus in participants with temporal lobe epilepsy caused an impaired associative memory in learning tasks that require learning and recall of objects and faces (Weniger, Boucsein, & Irle, 2004). These four brain areas plus the thalamus are those most consistently implicated in neurodegenerative dementias, especially in Alzheimer's Disease, even at an early stage (Manuello et al., 2018).

E. The bilateral superior temporal gyrus, superior temporal pole and middle temporal pole play an inter‐hemispheric integration role. Inter‐hemispheric functional connections between temporal lobes predict language impairment in adolescents born preterm (Northam et al., 2012). Phonological awareness, a key factor in reading acquisition was positively correlated with radial diffusivity of the interhemispheric pathways connecting temporal lobes (Dougherty et al., 2007). This bilateral temporal module could play a key role in many functions and dysfunctions.In our previous study on the same cohort, we focused on the repeatability of network topologies focusing on edge weights and graph theoretical metrics. We demonstrated that network topology and edge weights are repeatable, but the repeatability depends on the graph‐construction scheme (Messaritaki et al., 2019b). The important finding in this work is that the repeatability of network topologies and edge weights does not guarantee the repeatability of community detection at the mesoscale. In the present study, we focused on this important tool for mesoscale network topological investigations, and the detection of robust communities in structural brain networks over the same participants. To the best of our knowledge, this is the first study in the literature that explores the robustness of community detection over a large set of graph‐construction schemes (seven) and community detection algorithms (thirty‐three). Our analysis detected an optimal pair of {mscd_so, 9‐m OMST} that fulfils the three basic criteria: high repeatability of community affiliations between the two scan‐sessions, quality over surrogate null graph partitions and high similarity of group community affiliation with the individual community affiliations. To the best of our knowledge, this is the first time that the second and third criterion were used for the validation of representative consensus community affiliation (this includes studies using a single graph‐construction scheme and community detection algorithm).

The core of our study was an extensive analysis to identify the optimal pair of graph‐construction scheme and community detection algorithm. The choice of this pairing will also affect repeatability of connector and provincial hub detection based on the participation coefficient score Pi and the within‐module z‐score zi. Our results revealed a high repeatability of nodal Pi with the mscd_so algorithm across the seven graph construction schemes. The highest ICC score was reached for the {mscd_so, 9‐m OMST} pair. A significantly higher repeatability of nodal zi was found for mscd_so algorithm compared to the rest of the community detection algorithms. The highest ICC was achieved for the pairs {mscd_rb, NS‐thr} and {mscd_so, NS‐thr}. Statistical analysis revealed a main effect of network‐averaged ICC of zi and Pi across the four community detection algorithms, a main effect of the graph construction schemes, and no significant interaction effect.

The Agreement index revealed consistent connector hubs across the cohorts, but no group consistent provinicial hubs using both modular network metrics. Statistical analysis revealed a main effect of network‐averaged Agreement index of both provincial and connector hubs across the four community detection algorithms, an effect of the graph construction scheme, and no interaction effect. The pair {mscd_so, 9‐m OMST} demonstrated the highest Agreement index for provincial hubs, and the pair {mscd_rb, FA‐t/NS‐w} showed the highest Agreement index for connector hubs. We detected a group of 13 repeatable connector hubs across the cohort (Agreement = 1), but no group consistent provincial hubs (Agreement < 1). Based on our results, we therefore recommend to not use these modular network metrics for the detection of provincial hubs, at least when using the AAL atlas. The designation of a brain node as a hub depends also on the scale at which brain networks are constructed. Many brain areas in a basic atlas template group together functionally heterogenous subareas and it is possible that a finer‐grained parcellation may affect the nodes' classification as a hub or not. For example, the thalamus, despite comprising 50–60 specialised sub‐nuclei (Herrero, Barcia, & Navarro, 2002) is in many studies, including ours, treated as a single node.

Running the comparison study for the whole set of thirty‐three graph partition algorithms (including graph partition of the original graph and 1,000 surrogate null models) takes a few hours on a personal computer. We suggest to the neuroscience community to always run such an analysis over an in‐house test–retest data set acquired with the same settings as in the targeted data set. Optimising the set of algorithms over the test–retest study will increase the chance of repeatability of findings over the single‐scan data set. This process will increase the reproducibility of research findings, especially important for cross‐sectional studies (Welton, Kent, Auer, & Dineen, 2015).

Our study has a few limitations. This data set involves a specific data acquisition protocol and a specific tractography algorithm. We recommend following our analysis for every study because such an investigation could improve the repeatability and reproducibility of the findings at the mesoscale while also increasing the power of the study at the nodal and network level (Messaritaki et al., 2019a). Additionally, in our study we used only one of the three available b‐values to perform the tractography. This was mainly done in order to supplement our previous work (Messaritaki et al., 2019b), and we chose the b‐value of 2000 s/mm2, because this value provides a balance between the b‐value being sufficiently high to resolve crossing fibres with CSD, while at the same time ensuring sufficient SNR in the signal for robust measurements, and that higher‐order effects of the diffusion do not need to be taken into account when calculating the diffusion metrics. Using one b‐value also reflects acquisition protocols routinely used in other studies, and therefore makes our work more applicable to the general literature. At the same time, tractography results could be improved by combining data from all avaialble b‐values, and the implications of using different community detection algorithms in those cases should be explored as well. Moreover, three or more scan sessions would be also desirable to get a more robust assessment of repeatability. Scanning the same participant on different scanners and /or with different protocols would also allow assessment of reproducibility as well as repeatability. Lastly, the reproducibility of estimates of structural brain networks is affected by the resolution of the MR data (Vaessen et al., 2010), the parcellation scheme used (Bassett, Brown, Deshpande, Carlson, & Grafton, 2011), the interval time between the scan sessions and others.

CONCLUSIONS

In this study, we compared several graph‐construction schemes and thirty‐three community detection algorithms for the detection of reproducible communities in structural brain networks. Our extensive analysis showed that every choice in both groups of algorithms exhibits different reproducibility in community detection algorithms, as well as in connector/provincial hubs detection. Our analysis indicates that our analytic pathway should be adopted and performed in every study in order to extract reliable results at the mesoscale of structural brain networks.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

AUTHOR CONTRIBUTIONS

Stavros I. Dimitriadis: conceptualization, methodology, software, validation, formal analysis. investigation, data curation, roles/writing ‐ original draft and funding acquisition. Eirini Messaritaki: Formal analysis, data curation, roles/writing ‐ original draft, funding acquisition. Resources: The original diffusion MRI are free available from the Human Connectome Project. The analysis of tractography and the construction of structural brain networks has been realised by Eirini Messaritaki. Derek K. Jones: Writing–review and editing