Participants

The data used for this study were acquired by the Nathan Kline S. Institute for Psychiatric Research and made available by the 1000 Functional Connectomes Project INDI (Enhanced NKI Rockland Sample24, http://fcon_1000.projects.nitrc.org/indi/enhanced/). Institutional review board approval for this project was obtained at the Nathan Kline Institute (#239708). All methods were carried out in accordance with these guidelines and all participants gave written informed consent. We used a subsample of 309 participants, for which complete neuroimaging data were available (age: 18–60 years, M = 38.93, SD = 13.94; gender: 199 female, 110 male; handedness: 262 right, 22 left, 25 ambidextrous). For this sample, the Full Scale Intelligence Quotient (FSIQ), as assessed with the Wechsler Abbreviated Scale of Intelligence (WASI)25, ranged from 67 to 135 (M = 99.12, SD = 13.23).

fMRI Data Acquisition

MRI data were acquired on a 3 Tesla whole-body MRI scanner (MAGNETOM Trio Tim, Siemens, Erlangen, Germany). Functional resting-state data were obtained using a T2*-weighted BOLD-sensitive gradient-echo EPI sequence with 38 transversal axial slices of 3 mm thickness (120 volumes; field of view [FOV] 216 × 216 mm; repetition time [TR] 2500 ms; echo time [TE] 30 ms; flip angle 80°; voxel size 3 × 3 × 3 mm; acquisition time 5.05 min). For coregistration, three-dimensional high-resolution anatomical scans were obtained with a sagittal T1-weighted, Magnetization Prepared-Rapid Gradient Echo sequence scan (176 sagittal slices; FOV 250 × 250 mm; TR 1900 ms; TE 2.5 ms; flip angle 9°; voxel size 1 × 1 × 1 mm; acquisition time 4.18 min).

Preprocessing

Data was preprocessed using FSL (http://www.fmrib.ox.ac.uk/fsl/) and AFNI (http://afni.nimh.nih.gov/afni) with the scripts released by the 1000 Functional Connectomes Project (http://www.nitrc.org/projects/fcon_1000), comprising: 1. Discarding the first four EPI volumes to allow for signal equilibration, 2. Slice-time correction, 3. Three-dimensional motion correction, 4. Time-series despiking, 5. Spatial smoothing (6 mm full-width half-maximum Gaussian kernel), 6. Four-dimensional mean-based intensity normalisation, 7. Bandpass temporal filtering (0.005–0.1 Hz), 8. Removing linear and quadratic trends, 9. Normalisation of the individual EPI volumes to MNI152 space (3 × 3 × 3 mm) via nonlinear transformation and by the use of each subject’s anatomical scan, 10. Elimination of nine nuisance signals (white matter, cerebrospinal fluid, global mean, six motion parameters) by regression.

Graph Analyses

Graph analyses were performed with the open source python package network-tools 26, specifically developed for the analysis of functional and structural brain network graphs.

Graph construction

As nodes, we used those 5,411 voxels that covered all grey matter in the EPI images down-sampled to 6 × 6 × 6 mm. For each subject separately, edges were assumed between nodes showing high positive correlations of BOLD signal time series. Edges of physically short distance (< 20 mm) were excluded, due to their increased susceptibility to motion artefacts and potential correlations arising from shared nonbiological signal27. Most graph metrics are strongly influenced by the density of the graph28. This has specifically been shown for modularity29. To avoid biases due to individual differences in graph density, the main analyses were performed on thresholded and binarised graphs (as recommended for the study of individual differences in graph topology28,30). In contrast, weighted graphs usually vary in density (i.e., the mean weight of edges), even if the number of edges is held constant across individuals. For the purpose of comparison, we also conducted all analyses on weighted graphs (see Supplementary Tables S4, S5). Discussion of results, however, will rely on the results for the binarised graphs. We applied five different thresholds to the correlation matrix, retaining the strongest 10, 15, 20, 25, or 30% edges, thereby also excluding all negative edges31. This resulted in five graphs of different density per person. Community detection and the calculation of graph metrics were performed separately for the five graphs, and resulting graph metrics were averaged for each participant. This averaging procedure was applied to enhance the reliability of findings, as the resulting measures of graph properties are robust across a wider range of thresholds.

Measures of modular network organisation

To study the modular organisation of the functional brain network graphs, we applied the Louvain algorithm32. It finds the optimal modular partition in an iterative procedure maximizing the global modularity Q 33:

$$Q={\sum }_{s=1}^{m}[\frac{{l}_{ins}}{L}-\,{(\frac{{k}_{s}}{2L})}^{2}]\,,$$ (1)

where m is the number of modules, l ins is the number of edges inside module s, L is the total number of edges in the network, and k s is the total degree of the nodes in module s. Thus, the actual fraction of within-module edges is represented by the first term, whereas the expected fraction of within-module edges is represented by the second term. If the first term (actual within-module edges) is much larger than the second term (expected within-module edges), there are many more edges inside module s than expected by chance. In that case, s can be defined as a module and the global modularity Q, which results from summing up these differences (actual – expected within-module edges) over all modules m in the network, is increased. Usually, modularity values of Q > 0.3 indicate a modular network structure34.

The Louvain algorithm starts by assigning a different module to each node. Then, the first step is a greedy optimisation, where nodes adopt the modules of one of their neighbour nodes, if this reassignment increases the global modularity Q (see above). In the second step, a meta-network is built, whose nodes are the modules found in the first step. Both steps are repeated until no improvement in global modularity Q is possible and the optimal partition is found35,36,37. In addition to global modularity Q, for each participant, we calculated three further whole-brain measures of modular network organisation for the final module partition: number of modules, average module size, and the variability in module size.

The embedding of each node within the modular partition can be described by two graph-theoretical metrics:

(i) The participation coefficient p i represents between-module connectivity and is defined as:

$${p}_{i}=1-\,{\sum }_{m{\epsilon }M}{(\frac{{k}_{i}(m)}{{k}_{i}})}^{2}\,,$$ (2)

where k i is the degree of node i (i.e., the number of edges directly attached to node i) and k i (m) is the subset of edges that connect node i to other nodes within the same module13,38. The participation coefficient p i of a node is 0 when all of its edges are within its own module, and close to 1 when its edges are uniformly distributed among its own and other modules38.

(ii) The within-module degree z i represents within-module connectivity and is defined as:

$${z}_{i}=\,\frac{{k}_{i}\,({m}_{i})-\,\overline{k\,}\,({m}_{i})}{{\sigma }^{k({m}_{i})}},$$ (3)

where m i is the module of node i, k i (m i ) is the within-module degree of node i, \(\bar{k}\) (m i ) and σ k(mi) are the mean and standard deviation of the within-module degree distribution of module m i 38. Positive values indicate that a node is highly connected to nodes within its own module, whereas negative values indicate low levels of connectivity within the same module15.

These two graph metrics, participation coefficient p i and within-module degree z i , allow the characterisation of a node’s embedding within the modular brain network free of any biases due to different module sizes18 (as would be the case when simply comparing numbers of between- and within-module connections). The distributions of participation coefficient p i and within-module degree z i were visualised by averaging the individual mean p i - and z i -values of each node across participants and projecting them to the surface of the brain (Fig. 1).

Figure 1 Values of participation coefficient p i (top row) and within-module degree z i (bottom row) averaged across all participants. Higher values are shown in warm colours (green to red), lower values are shown in cool colours (blue to pink). Graph metrics (p i and z i ) were calculated for binarised and proportionally thresholded graphs using five different cut-offs (i.e., graphs were defined by the 10%, 15%, 20%, 25%, or 30% strongest edges). For each participant, individual mean maps for the graph metrics were calculated by averaging across the five thresholds. Displayed here are the group average maps for p i and z i that resulted from averaging all individual mean maps across participants (see Methods section for details on the procedure). For the lateral view, values were projected to the surface of the brain. The medial view displays graph values in the x-plane. L, left; R, right. Full size image

Node-type analysis

Functional cartography38 uses the above-described metrics (participation coefficient p i , within-module degree z i ) to classify network nodes into seven different types regarding their roles in within- and between-module communication (see Fig. 2A). As proposed in previous work38,39, we classified nodes with within-module degree z i ≥ 1 as hubs (18% of all nodes) and nodes with z i < 1 as non-hubs. Depending on the participation coefficient p i , non-hubs were further divided into ultra-peripheral (p i ≤ 0.05), peripheral (0.05p i ≤ 0.62), non-hub connector (0.62 < p i ≤ 0.80), and non-hub kinless nodes (p i > 0.80), whereas hubs were classified into provincial (p i ≤ 0.30), connector (0.30 < p i ≤ 0.75), or kinless hubs (p i > 0.75).

Figure 2 Node-type analysis. (A) Definition of node types as a function of participation coefficient p i and within-module degree z i . Adapted from Guimerà and Amaral (2005; method also known as functional cartography; cf. Methods). (B) Group-average proportions of node types across the entire cortex. Proportions of node types were calculated for each individual subject separately and then averaged across all subjects. (C) Illustration and anatomical distribution of node types within the human brain for one exemplary participant. Non-hub nodes (z i ≤ 1) are shown in cool colours (green to blue), hubs (z i > 1) are shown in warm colours (yellow to red). Graph metrics (p i and z i ) and the respective node-type proportions were calculated for binarised and proportionally thresholded graphs using five different cut-offs (i.e., graphs were defined by the 10%, 15%, 20%, 25%, or 30% strongest edges). For each participant, individual node-type proportions were calculated by averaging across the five thresholds. Displayed in B are the group average proportions of node types resulted from averaging all individual node-type proportions across participants (see Methods section for details on the procedure). The x- and z-coordinates represent coordinates of the Montreal Neurological Institute template brain (MNI152). Full size image

Intelligence-related differences in modular network organisation

All individual-differences analyses were done after exclusion of outliers, i.e., subjects with values > 3 SD above/below the mean of the respective variable of interest. For all whole-brain measures and proportions of node types we used SPSS22 (IBM Corp., Armonk, NY) to calculate partial correlations with WASI FSIQ, including potential confounding effects of age, sex, and handedness as covariates. Effects with p < 0.05 were considered statistically significant (Bonferroni-adjusted p-values are .013 for global modularity measures and .007 for node-type proportions). To quantify evidence for the null hypothesis (i.e., absence of an association) we calculated Bayes Factors40,41 (BF 01 ), using Bayesian linear regression and the default prior42 as implemented in JASP (https://jasp-stats.org). In accordance with Jeffreys40, we interpret BF 01 > 3 as substantial evidence for the null hypothesis.

To investigate the association between intelligence and whole-brain aspects of modular network organisation, we calculated partial correlations between WASI FSIQ and global modularity Q, number of modules, average module size, and the variability in module size. Furthermore, we tested for associations between intelligence and the whole-brain proportions of each node type as determined in the node-type analysis. To study the association between intelligence and node-specific aspects of modular organisation (i.e., between- and within-module connectivity), we set up two separate regression models in SPM8 (Statistic Parametric Mapping, Welcome Department of Imaging Neuroscience, London, UK), one for predicting the individual maps of participation coefficient p i , and one for predicting the individual maps of within-module degree z i (both maps upsampled to 3 × 3 × 3 mm) by WASI FSIQ. To control for the potential confounding effects of age, sex, and handedness, these variables were included as covariates of no interest in all regression models. P-values were corrected for multiple comparisons using a Monte Carlo-based cluster-level thresholding procedure43. An overall threshold of p < 0.05 (FWE-corrected) was applied by combining a voxel-level threshold of p < 0.005 with a cluster-level threshold of k > 26 voxels (3dClustSim; AFNI version August 2016; 10,000 permutations; voxel size: 3 × 3 × 3 mm)44. As ultimately the modular organisation of the whole brain network is always defined by both, between-module and within-module connectivity, we also tested for an overlap of intelligence-related effects in both measures, i.e., participation coefficient p i and within-module degree z i .

Data availability statement

The data used in the present work can be accessed under the following link: http://fcon_1000.projects.nitrc.org/indi/enhanced/.