Abstract
Investigation of the functional macro-scale organization of the human cortex is fundamental in modern neuroscience. Although numerous studies have identified networks of interacting functional modules in the gray-matter, limited research was directed to the functional organization of the white-matter. Recent studies have demonstrated that the white-matter exhibits blood oxygen level-dependent signal fluctuations similar to those of the gray-matter. Here we used these signal fluctuations to investigate whether the white-matter is organized as functional networks by applying a clustering analysis on resting-state functional MRI (RSfMRI) data from white-matter voxels, in 176 subjects (of both sexes). This analysis indicated the existence of 12 symmetrical white-matter functional networks, corresponding to combinations of white-matter tracts identified by diffusion tensor imaging. Six of the networks included interhemispheric commissural bridges traversing the corpus callosum. Signals in white-matter networks correlated with signals from functional gray-matter networks, providing missing knowledge on how these distributed networks communicate across large distances. These findings were replicated in an independent subject group and were corroborated by seed-based analysis in small groups and individual subjects. The identified white-matter functional atlases and analysis codes are available at http://mind.huji.ac.il/white-matter.aspx. Our results demonstrate that the white-matter manifests an intrinsic functional organization as interacting networks of functional modules, similarly to the gray-matter, which can be investigated using RSfMRI. The discovery of functional networks within the white-matter may open new avenues of research in cognitive neuroscience and clinical neuropsychiatry.
SIGNIFICANCE STATEMENT In recent years, functional MRI (fMRI) has revolutionized all fields of neuroscience, enabling identifications of functional modules and networks in the human brain. However, most fMRI studies ignored a major part of the brain, the white-matter, discarding signals from it as arising from noise. Here we use resting-state fMRI data from 176 subjects to show that signals from the human white-matter contain meaningful information. We identify 12 functional networks composed of interacting long-distance white-matter tracts. Moreover, we show that these networks are highly correlated to resting-state gray-matter networks, highlighting their functional role. Our findings enable reinterpretation of many existing fMRI datasets, and suggest a new way to explore the white-matter role in cognition and its disturbances in neuropsychiatric disorders.
Introduction
The functional organization of the gray matter as composed of distinct interacting “modules” is a cornerstone of cognitive neuroscience and clinical neuropsychiatry (Brodmann, 1909; Fodor, 1983; Mesulam, 1998). In recent years, resting-state functional MRI (RSfMRI) has emerged as a powerful technique for investigation of the human cortex's functional architecture. By measuring the correlations between fMRI signal fluctuations in different cortical regions, distinct networks were identified, providing insight into the coordinated activity of these regions (Biswal et al., 1995; Greicius et al., 2003; Yeo et al., 2011; Park and Friston, 2013). However, the white-matter components enabling this cross talk were only explored using structural methods such as diffusion-tensor imaging (DTI; Assaf and Pasternak, 2008; Hagmann et al., 2008; Honey et al., 2009). Although these methods can uncover the fine details of white-matter structure, they do not provide direct evidence of their function, and cannot replace direct in vivo recording of activity inside the human white matter.
Recent works have demonstrated that the white matter also carries functional information, which may be detected by functional MRI, as distinct white-matter tracts were shown to be activated in response to perceptual and motor task performance (Mazerolle et al., 2008; Yarkoni et al., 2009; Fabri et al., 2011; Gawryluk et al., 2011, 2014; Fabri and Polonara, 2013). Moreover, activity in the white matter corresponded to task demands, such that parts of the corpus callosum were active for tasks requiring interhemispheric transfer, and the internal capsule was active for tasks requiring motor activity (Mazerolle et al., 2008; Gawryluk et al., 2011). Although this activity was weaker than that of the gray matter (presumably due to blood vessels in the white matter being five to seven times less dense than in the gray matter; Gawryluk et al., 2014), it was still reliably measured in individual subjects. These studies established the existence of functional brain activity in the white-matter, and showed that fMRI is a promising tool for the investigation of white-matter function in vivo.
Several recent studies further explored white-matter functional activity during RSfMRI. Clustering of resting-state white-matter blood oxygen level dependent (BOLD) signals yielded patterns resembling white-matter tracts (Mezer et al., 2009). However, these clusters were demonstrated only in a single slice, and their relation to white-matter tracts and gray-matter networks were not investigated. Two additional studies demonstrated that fMRI activity along specific white-matter tracts is highly homogenous during the resting-state, and showed that these white-matter signals have characteristics suggestive of hemodynamic (BOLD) changes associated with neuronal activity, and are altered when performing language tasks (Ding et al., 2013, 2016). A final recently published study used independent-components analysis and hierarchical clustering to demonstrate the existence of clusters of correlated activity within the white matter (Marussich et al., 2017). This study further demonstrated that in white-matter clusters in the occipital lobe, connectivity changes occurred when subjects passively watched movies compared with resting-state, providing additional support for the functional character of these signals. However, these studies mostly focused on specific white-matter pathways, used relatively small sample sizes, and did not characterize the correspondence between white-matter signals and gray-matter resting-state networks or structural white-matter tracts across the whole brain.
To characterize large-scale white-matter functional networks, we applied here clustering on RSfMRI data of white-matter regions from a large cohort of subjects (N = 176), and investigated the resulting white-matter functional networks and their relation to both structural white-matter tracts, identified by DTI, and to activity in known gray-matter networks.
Materials and Methods
Subjects and MRI acquisition.
Resting-state data from 176 healthy subjects (120 males, mean age 37.2 ± 20 years) was acquired from the NKI/Rockland sample public dataset (http://fcon_1000.projects.nitrc.org/indi/pro/nki.html; Biswal et al., 2010; Nooner et al., 2012). All subjects were scanned using a Magnetom TrioTim syngo 3T scanner (Siemens Medical Solutions). The data for each subject consisted of a 10 min resting-state EPI scan with eyes open (38 slices, TR = 2.5 s, 3 × 3 × 3 mm isovoxels, flip angle 800, 260 volumes, interleaved slices) and an MPRAGE anatomical scan (1 × 1 × 1 mm resolution). Additional data were acquired from the Cambridge-Buckner dataset from the 1000 connectomes project (https://www.nitrc.org/frs/?group_id=296; Biswal et al., 2010); the data in this dataset (94 subjects, mean age 21.4 ± 3 years) consisted of 6 min resting-state scans with eyes open (47 slices, TR = 3 s, 119 volumes, interleaved slices) and MPRAGE anatomical scans.
DTI: 12 subjects (8 males, age 38.1 ± 13.57) were scanned in the Hadassah Hebrew University Medical Center using a Trio 3T scanner (Siemens Medical Solutions). All subjects had no personal history of neurologic or psychiatric disorders, and had normal structural MRI results. DTI acquisition parameters were as follows: single-shot, spin-echo diffusion protocol, TE = 94 ms, TR = 7127–8224 ms, FOV = 260 × 260 mm, matrix = 128 × 128, 52–60 axial slices, 2 mm thick slices, b = 0 and b = 1000 s/mm2. The high b value was obtained by applying gradients along 64 different diffusion directions, including two averages. All participants gave written informed consent, and the study was approved by the ethical committee of the Hadassah Hebrew University Medical Center.
Functional and anatomical MRI preprocessing.
Preprocessing was performed using SPM8 (www.fil.ion.ucl.ac.uk/spm), DPARSFA (Yan and Zang, 2010), and in-house MATLAB (MathWorks) scripts (Peer et al., 2014, 2016). Anatomical images were segmented into white matter, gray matter, and CSF, using SPM8's New Segment algorithm. Functional time-series were slice-time-corrected, motion-corrected to the mean functional image using a trilinear interpolation with 6 degrees of freedom, and coregistered with the anatomical image. Subjects with maximum motion >3 mm were excluded from further analyses, resulting in a cohort of 176 subjects of the original 204 subjects in the NKI/Rockland database. Additional preprocessing steps included the removal of linear trends to correct for signal drift and filtering with a 0.01–0.15 Hz bandpass filter to reduce non-neuronal contributions to BOLD fluctuations. To minimize the effect of subjects' motion on functional connectivity characteristics (Power et al., 2012; Satterthwaite et al., 2012; Van Dijk et al., 2012), we performed multiple regression of 24 motion parameters (Friston et al., 1996): six rigid-body head motion parameter values—x, y, and z translations and rotations—their value at the previous time point, and the 12 corresponding squared values. To further reduce motion effects, motion “spikes” were also included as separate regressors (identified by framewise displacement of 1 mm), effectively eliminating (censoring) the data at the spike without further changes to correlation values (Power et al., 2012; Satterthwaite et al., 2013). Finally, a regressor for the mean CSF signal was also included as a nuisance source. To avoid elimination of signals of interest, we did not include white-matter and global brain-signal regressors.
Following these preprocessing steps, images were spatially smoothed [4 mm full-width half-maximum (FWHM), isotropic], where smoothing was performed separately on the white matter and the gray matter of each subject, to avoid mixing their signals. For separation of smoothing, we used the segmentation results from each subject to identify white- or gray-matter voxels (identified using a threshold of 0.5 using SPM8's tissue segmentation). We then created two functional series of images from the original functional data, one containing only data of white-matter voxels and one containing only data of gray-matter voxels (by zeroing any voxel not identified as white matter or gray matter, respectively). The smoothing algorithm of SPM8 (FWHM = 4 mm) was then applied on both sets of images. Finally, we combined the two images into full functional images, using only the smoothed data from the white-matter and gray-matter voxels, respectively. Following smoothing, all brain images were normalized to standard anatomical space (Montreal Neurological Institute EPI template, resampling to 3 mm cubic voxels) using the DARTEL algorithm (Ashburner, 2007).
To measure the effect of different preprocessing options on the data, we also attempted full removal (censoring) of high-motion time-points (framewise displacement >1 mm) from the data (Power et al., 2012) instead of using motion spike regressors. We further attempted to use mean global signal regression in the nuisance covariates removal stage using DPARSFA's SPM mean mask.
Creation of group white-matter and gray-matter masks.
To obtain masks for selection of voxels for clustering across the group, we used the segmentation results of each subject. For each voxel, we identified it as white matter, gray matter, or CSF according to the maximum probability from the three segmentation images, for each subject separately. This resulted in a unique designation of tissue type for each voxel. These masks were averaged across subjects to obtain for each voxel the percentage of subjects in which it was classified as white or gray matter. For white matter, voxels identified in >60% of subjects as white matter were used for mask creation. For gray matter, we used a more lenient threshold of >20% of subjects to avoid missing cortical regions with high variability, but excluded any voxels included in the white matter mask. We then compared the resulting masks to the functional data, and removed voxels identified as white or gray matter yet having functional data in <80% of the subjects (such as parts of the medulla and spinal cord). Finally, to correctly classify deep brain structures which are incorrectly assigned as white matter using existing segmentation algorithms (Babalola et al., 2009; Wonderlick et al., 2009; Lorio et al., 2016), we used the Harvard-Oxford Atlas (Desikan et al., 2006) to mark as gray matter (and remove from the white-matter mask) voxels included in the thalamus, caudate nucleus, putamen, globus pallidus, and nucleus accumbens. The resulting white-matter and gray-matter group masks contained 20,059 and 42,044 voxels, respectively.
Generation of subject-specific and group-averaged correlation matrices.
To create data for clustering from each subject, we computed Pearson's correlation coefficients between each white-matter voxel and other white-matter voxels. To reduce computational complexity of the subsequent clustering, we subsampled the white-matter mask using an interchanging grid, taking any second voxel along the image rows and columns, and shifting it by 1 between slices to avoid missing entire columns of data. For each white-matter voxel, correlation was calculated to all of the subsampled mask nodes, resulting in a whole-brain correlation pattern for each voxel (20,059 × 4998 matrix; Yeo et al., 2011; Craddock et al., 2012). Correlation matrices were computed for each subject and then averaged across all 176 subjects to obtain a group-level correlation matrix. Importantly, for each voxel-to-voxel correlation, data used for averaging across subjects were taken only from subjects for whom both voxels were identified as white matter during segmentation, therefore excluding contributions from gray-matter signals. The same process was performed to obtain subject-specific gray-matter correlation matrices and a group-averaged gray-matter correlation matrix (42,044 × 4645 matrix).
Clustering of voxels according to their connectivity patterns.
To identify white- and gray-matter functional networks, we used a clustering approach on the resting-state correlation matrices (Yeo et al., 2011; Craddock et al., 2012; Blumensath et al., 2013; Moreno-Dominguez et al., 2014; Peer et al., 2014). K-means clustering (implemented in MATLAB) was performed on the rows of the average group-level correlation matrix (distance metric–correlation, 10 replicates). This resulted in identification of clusters of white-matter voxels with similar connectivity patterns to the rest of the voxels. Clustering was performed for all numbers of clusters (K) ranging between 2 and 22, to obtain coarse and fine-grained parcellation results. Clustering was additionally performed on connectivity matrices obtained from half of the subjects, for stability analysis, and from random subgroups of different sizes (n = 80, 70, 60, 50, 40, 30, 20, 10, 5 and 1 subjects; 10 randomly selected groups of each size).
Identification of the most stable numbers of networks.
To identify the most stable number of networks, we measured the stability of the clustering solution for each number of clusters, from 2 to 22 (Lange et al., 2004; Yeo et al., 2011). The full connectivity matrix was divided randomly into four folds, each containing the correlation of all voxels to a fourth of the subsampled masks (20,059 × 1249 voxels per fold for white matter; 42,044 × 1161 voxels for gray matter). For each number of clusters, clustering was separately performed on each of the four folds; clustering according to different features should provide approximately similar results for numbers of clusters which represent stable clustering solutions (Lange et al., 2004; Yeo et al., 2011). As clustering gives random labels to clusters, it is difficult to directly measure similarity between clustering solutions using labels. Therefore, to measure this similarity between clustering solutions, we computed an adjacency matrix for each solution (where each cell represents whether 2 voxels belong to the same cluster in this solution), and these adjacency matrices were compared to each other using Dice's coefficient. For each number of clusters, all four adjacency matrices were compared to each other and an average Dice's coefficient was obtained. Although clustering similarity naturally decreases with larger numbers of clusters (Yeo et al., 2011), plotting the Dice's coefficient allows identification of highly stable solutions as local peaks on the graph (Figs. 1A,B, 7A).
In addition to evaluation of clustering solutions according to stability of clustering for different sets of features, we also evaluated stability of solutions for different sets of subjects. Subjects were randomly divided into two groups of 88 subjects each, and full group-average correlation matrices were computed for each group (20,059 × 4998 matrix for each group). Clustering was separately performed on the full connectivity matrix from each subject group. The two solutions were again compared for each number of clusters using Dice's coefficient as described above.
Our results of white-matter clustering demonstrated several stable clustering solutions for the white matter, at 2, 5, 8, and 12 clusters (Figs. 1, 2). In this paper, we chose to focus on the most detailed level of description of 12 functional networks. Despite this, other stable numbers of clusters might also be of interest, and may be characterized in future works.
Symmetry calculation.
To measure the symmetry of the identified white-matter functional networks, we used the clustering solutions, and flipped half of the data along the midsagittal plane to obtain clustering solutions for the two hemispheres. We then computed an adjacency matrix for each hemisphere clustering solution separately (as described in “Identification of the most stable numbers of networks”). The clustering solutions were compared for the two hemispheres using Dice's coefficient, using only voxels designated as white matter in the two hemispheres. To measure the significance of the symmetry values, we used a permutation test, where the adjacency matrix for each hemisphere was randomly permuted 100 times. Dice's coefficient was computed for the two permuted adjacency matrices at each iteration, and the nonpermuted result was then compared with the resulting permutations distribution.
Verification of network boundaries using correlation transitions.
To support the validity of the clustering results, we searched for sharp transitions in correlation patterns for identification of boundaries between functional networks (Cohen et al., 2008; Hirose et al., 2012). For each white-matter voxel, we calculated the correlation of its signal across time to that of all white-matter voxels in the 5 × 5 × 5 voxels box surrounding it (using each subject's own white-matter mask for voxel selection, to avoid mixing with gray-matter signals). Voxels lying at the border between different functional networks are expected to show lower average correlations to surrounding voxels (Hirose et al., 2012). To compare these results to the 12-networks clustering solution, we calculated how many different networks (clusters) are included in each voxel's surrounding 5 × 5 × 5 voxels box, and correlated these values with the correlation to neighboring voxels values calculated earlier.
DTI preprocessing and identification of white-matter tracts.
To obtain average white-matter tracts, we used data from a separate group of subjects (see above “Subjects”). DTI preprocessing was performed using the MrVista software (https://github.com/vistalab/vistasoft), and included removing Eddy current distortions and subject motion. Data were registered and aligned to the T1 image using mutual information algorithm.
To identify average white-matter tracts, we used the recently developed AFQ algorithm, which automatically identifies 20 major fiber tracts in each subject based on the JHU white-matter tractography atlas (Mori et al., 2005; Hua et al., 2008; Yeatman et al., 2012). White-matter tracts included the bilateral corticospinal tracts, anterior thalamic radiations, cingulum (cingulate gyrus and hippocampal parts), forceps major, forceps minor, inferior fronto-occipital fasciculus, inferior longitudinal fasciculus, superior longitudinal fasciculus (main and temporal parts), and uncinate fasciculus.
Following tracts identification, we used the anatomical T1-weighted images of each subject to normalize his anatomical and tract data into MNI space using SPM8. For each subject, voxels where a tract was identified in >2 subjects, and included in the group white-matter mask, were designated as belonging to the tract.
Correspondence between white-matter functional networks and DTI tracts.
Similar anatomical DTI tracts from the two hemispheres were combined, to compare them to the bilateral functional networks. To measure the similarity between each white-matter functional network and each anatomical DTI tract, we calculated the percentage of voxels in the functional network identified as belonging to that tract. It is important to note that this comparison is not complete, as all identified white-matter tracts together cover only 45% of the subjects' average white-matter mask, whereas functional connectivity clustering was applied to the full mask.
Comparison of signals between gray-matter and white-matter functional networks.
To measure the relations between gray- and white-matter functional networks, we extracted the average signal time courses from white- and gray-matter networks by averaging across all voxels belonging to each network, separately for each subject. To avoid mixing of gray-matter and white-matter signals in each subject, voxels used for signal averaging in each subject included only those belonging to the required network and identified as white or gray matter, respectively, in this specific subject, by tissue of maximum probability in the tissue segmentation mask. Signal amplitudes in each frequency were extracted using the Fourier transform (MATLAB's FFT function) from each white- and gray-matter network. The resulting frequency graphs were averaged across participants to produce an average frequency-power graph for each network, and linear fit was computed for the mean graphs. We furthermore computed Pearson's correlation between all white-matter and gray-matter networks' time courses for each subject. These correlation values were averaged across subjects to obtain a group-level matrix representing the relation between gray-matter and white-matter networks. To measure the consistency of these results, these steps were also performed in two independent samples (88 subjects each) from the whole group, and the absolute difference between the resulting correlation matrices was measured for each cell (correlation value).
Seed-based analysis.
For each of the 12 identified white-matter functional networks, we defined 5 random seeds of 10 spatially contiguous voxels each. Seeds were selected by setting a box of size 6 × 6 × 6 voxels and moving it randomly across the image until finding a location containing at least 10 voxels belonging to the required network (as defined on the whole group: N = 176 subjects). For each of 10 randomly selected subjects, we extracted the seed's time course and correlated it with the time courses of all other white-matter voxels of the corresponding subject. To measure whether the seed-based analysis indeed identified the required network, we calculated the percentage of voxels identified as correlated to the seed (using a threshold of r > 0.4) and belonging to the same network that was used to define the seed. We also performed this analysis on groups of 5 and 10 subjects, by averaging the correlation coefficients of each seed to each white-matter voxel across subjects (using for each voxel only subjects in which it was defined as white matter).
Visualization of functional white-matter networks.
Visualization of functional networks was performed using the ITK-SNAP (Yushkevich et al., 2006) and ParaView (http://www.paraview.org/) software packages (Madan, 2015). 3D images were smoothed for visualization purposes using a Laplacian filter (iterations = 100).
Experimental design and statistical analysis.
Similarity between different clustering solutions (data from different subject groups or with different processing options) was calculated using Dice's similarity coefficient (see “Identification of the most stable numbers of networks”). To calculate the significance of these similarity coefficients, we used a permutation test and compared the coefficients to a distribution of similarity coefficients between randomly assigned clustering results and the whole-group clustering result. To assess the symmetry of clustering solutions, we used a permutation test, where adjacency matrices were randomly permuted 100 times to create a null hypothesis Dice's coefficient distribution (see “Symmetry calculation”). Significance of correlation coefficients was obtained from MATLAB's “corr” function, which compares the observed correlation to a Student's t distribution under the null hypothesis that the correlation is zero. The significance of linear fits to the Fourier transform results from each network was assessed using MATLAB's “regress” function, which performs an F test on the regression model.
Results
Clustering analysis of RSfMRI reveals symmetrical functional networks in the white matter
To identify white-matter functional networks, we used a clustering approach on the white-matter resting-state voxelwise correlation matrices obtained from 176 healthy subjects (Yeo et al., 2011; Craddock et al., 2012; Blumensath et al., 2013; Moreno-Dominguez et al., 2014; Peer et al., 2014). Only voxels identified as white matter were used from each subject (see Materials and Methods). K-means clustering of white-matter functional connectivity revealed an intricate pattern of interlaced, symmetrical functional networks inside the white matter (Figs. 1, 2). The numbers of clusters that represented the most stable functional segregations were 2, 5, 8, and 12 (Fig. 1A); we focus here on the detailed clustering into 12 networks (Fig. 2; see Materials and Methods). The resulting networks can be divided into three layers: superficial (outer), middle, and deep (inner; Fig. 3), where superficial networks mostly correspond to functional subdivisions of the overlying gray matter. The functional networks were found to be highly symmetrical across the two cerebral hemispheres (permutation test, p < 0.01; Fig. 2), and several of them further demonstrated fine connecting commissural bridges between hemispheres (Fig. 2A; networks 1,4,5,11,12), suggesting a nonrandom, functionally related origin. Areal boundary mapping by local correlation strength (Cohen et al., 2008; Hirose et al., 2012) supported the clustering-based division to functional networks, as there was a strong negative correlation between the number of networks included in each local white-matter region to the signal homogeneity (correlation) between voxels in this region (r = −0.4, t(20,057) = −62.6, p < 1020).
Functional networks in the brain's white-matter: clustering solutions for different numbers of clusters. K-means clustering was applied on resting-state BOLD fluctuations from the white matter of 176 subjects. A, Stability of clustering for different numbers of clusters, measured by similarity of clustering solutions obtained using different subsets of features. B, Stability of clustering for different numbers of clusters, measured by similarity of clustering solutions between two independent subject groups. Asterisks represent local graph maxima. Both analyses indicate stable solutions at 2, 5, 8, and 12 clusters. C, Results of clustering of the white matter into 2, 5, and 8 functional networks, revealing highly symmetrical functional networks in the white matter. For the 12-clusters solution see Figure 2.
Functional networks in the brain's white-matter: 12 networks. K-means clustering into 12 clusters revealed highly symmetrical functional networks in the white matter, closely resembling known gray-matter networks and white-matter tracts. A, The 12 identified networks. B, Whole-brain representation of all networks in the white-matter.
White-matter functional networks are organized in three layers. Networks can be classified as superficial (1), middle (2), and deep (3), each layer directly overlying the prior one. A, Schematic representation of the relations between the different layers. B, The networks belonging to each layer.
Clustering stability with different processing and acquisition parameters
To measure stability, we performed the same clustering on data of two independent groups of 88 subjects from the original dataset. Clustering indicated a remarkable consistency between solutions for the two groups (Fig. 4A,B; Dice's similarity coefficient = 0.82, permutation test, p < 0.001), indicating that the clustering results are robust and do not depend on random noise or on the characteristics of specific subjects.
Stability of clustering across different subjects, processing and acquisition parameters. A, B, Results of clustering of two group of subjects with similar acquisition and processing parameters show high correspondence (Dice's coefficient = 0.82), indicating that the clustering is not random or influenced by noise sources. C, Clustering of subjects after application of global-signal regression. D, Clustering of a separate subject group with different acquisition parameters and demographics (TR, number of slices, number of time-points, age distribution).
We further measured the effects of several preprocessing stages on the results. Censoring (removal) of time points (Power et al., 2012) instead of adding regressors to eliminate motion spikes (Satterthwaite et al., 2013) had no effect on the clustering results (clustering to 12 networks, Dice's similarity coefficient between results with either censoring or spike regressors = 0.9918). Applying global signal regression resulted in symmetrical networks that were similar to the networks derived from the original clustering; however, several deep and middle networks (networks 5, 7, and 11) were merged into a single network, whereas networks 1 and 12 were subdivided into anterior and posterior subnetworks (Fig. 4C). Performing the clustering on a separate dataset (Cambridge-Buckner dataset, 1000 functional connectomes project; Biswal et al., 2010), with different acquisition parameters (TR, number of slices, number of time-points, subjects' age distribution; see Materials and Methods), resulted in clustering into generally similar and highly symmetrical networks, yet again exhibited merging and subdivisions of specific networks (merging of networks 1 and 4, and networks 7 and 11; subdivision of networks 10 and 12; Fig. 4D).
We next measured the stability of clustering in smaller groups of subjects. Results for all group sizes (down to the single-subject level) were significantly more similar to the whole-sample clustering than chance level (permutation test, p < 0.01), and clusters' symmetry was preserved (Fig. 5). However, for group sizes smaller than n = 20, clustering progressively became less stable. These results were comparable for gray and white matter.
Clustering stability in different sample sizes. Clustering was performed on randomly selected subgroups of subjects of differing sizes (10 subgroups of each size). A, Similarity between subgroups clustering and whole-group (N = 176) clustering (mean ± SD) for gray matter (red) and white matter (blue). Dashed line indicates similarity of random label assignment to the whole-group clustering. B, Examples of clustering results on subgroups of different sizes. Note that results in the n = 60 and n = 30 subgroups are similar to the whole-sample clustering (Fig. 2), unlike the smaller subgroups.
Functional white-matter networks correspond to structural white-matter tracts identified by DTI
To investigate the similarity between the identified networks and known white-matter structures we compared the functional white-matter networks to structural white-matter tracts, identified using DTI in 11 independent subjects. This analysis revealed complex interactions between functional networks and anatomical tracts. Although six of the 12 networks showed good (>40%) spatial correspondence with specific anatomical tracts (e.g., network 7 and the inferior longitudinal fasciculus; Fig. 6A), many of these networks extended further than the anatomical tract (e.g., network 2 extends beyond the uncinate fasciculus, network 4 extends beyond the forceps minor; Fig. 6B,C), and the remaining six networks simultaneously corresponded to several tracts (Table 1). The functional-anatomical overlap was stable when computed on clustering from two independent subject groups (n = 88 subjects each; correlation between the 2 functional-anatomical similarity matrices = 0.99). These findings support the reliability of our functional analyses and indicate that during resting state, white-matter tracts demonstrate homogenous functional activity, and several white-matter tracts may act in concert, presumably to support coordinated activity in distributed cortical networks.
Comparison of functional resting-state white-matter networks with known white-matter tracts identified using DTI. A–C, Examples of networks demonstrating extensive overlap with DTI tracts. Note that in many cases functional networks overlap with several DTI tracts. D, Systematic quantification of spatial correspondence between all functional networks and DTI tracts. Matrix cells' values represent the percentage of each anatomical tract occupied by the functional network. Functional networks show a high overlap with structural tracts, although several functional networks extend across several tracts. See Table 1 for full details of overlap. Sup., Superior; inf., inferior; long., longitudinal.
Overlap of white-matter functional networks with anatomical tracts, correlation with functional gray-matter networks, and putative designation
White-matter functional networks act in concert with gray-matter networks
Application of the same clustering procedure on signals from the gray-matter resulted in identification of nine gray-matter networks (Fig. 7A), closely corresponding to known cortical network subdivisions of the cerebrum and cerebellum (Buckner et al., 2011; Power et al., 2011; Yeo et al., 2011). To examine the relations between gray- and white-matter networks, we quantified the correlation between the average resting-state BOLD signals from different networks, for each subject separately. Across subjects, specific white-matter and gray-matter networks were highly correlated in their activity, with correlations between several network pairs rising above r = 0.8 (Fig. 7B,C; Table 1). In particular, superficial frontotemporal and frontoparietal white-matter networks were highly correlated (r > 0.7) to distributed networks, such as the default-mode, dorsal attention, ventral attention, and frontoparietal control networks, suggesting a role for white-matter networks in communication between distant regions of these gray-matter networks (Fig. 7C; Table 1). Other superficial networks were mostly correlated to their overlying gray-matter networks, possibly mediating close-range communications within these networks (Fig. 7; Table 1). In contrast to these results, the signals of deep white-matter networks were in general only weakly correlated to those of gray-matter networks (average r < 0.4; Fig. 7B; Table 1). The white- and gray-matter signals correlation results were consistent when measured in two independent samples from the data of 88 subjects each (average absolute difference in correlation values between the samples = 0.02, maximum difference between groups = 0.06).
Relation between white-matter and gray-matter functional networks. Gray-matter networks were identified using K-means clustering. A, Stability of clustering of the gray-matter signals for different numbers of clusters, measured by similarity of clustering solutions obtained using different subsets of features (asterisks represent local graph maxima). The stable solution at K = 9 corresponds to known parcellations of resting-state networks from the literature. B, Average correlation between signals from white-matter and gray-matter functional networks. Many white- and gray-matter networks show a high degree of correlation in their functional BOLD signal, suggesting a role for white-matter networks in mediating communication within gray-matter resting-state networks. C, White-matter networks showing high connectivity to gray-matter resting-state networks. White-matter and gray-matter networks are presented within and around the brain, respectively, and lines indicate pairs of white- and gray-matter networks with average signal correlation between them surpassing 0.7. Functional white-matter networks with lower levels of connectivity to gray-matter networks are not presented. Note that frontoparietal white-matter networks (blue, purple) are highly correlated to gray-matter networks with separate frontal and parietal parts (r > 0.7).
To further compare the signals of gray- and white-matter networks, Fourier transform was applied to the signals. The amplitudes of white-matter networks were lower than amplitudes of gray-matter networks for all frequencies (Fig. 8). Furthermore, all white-matter and gray-matter networks showed a gradual decrease in amplitude with increased frequency (linear fits for all networks, F(87) = 49.5–167.3, p values = 4e−10–5e−22), indicating greater activity at low frequencies. All gray-matter networks and superficial white-matter networks showed maximum activity at the lowest frequencies of ∼0.01 Hz, whereas the middle and deep white-matter networks (networks 5, 7, 10, and 11) showed maximal amplitude at ∼0.07 Hz (Fig. 8).
Power spectral analysis of signals in white- and gray-matter networks. A, Power-frequency graphs for the 12 white-matter networks. B, Power-frequency graphs for the nine gray-matter networks. Power is greater at lower frequencies in all networks. Note the different peak amplitude for deep white-matter networks (5, 7, 10, and 11; for network numbers see Fig. 2).
White-matter seed-based analysis in individual subjects
The clustering approach was found to be reliable across large groups of subjects. To test whether white-matter functional networks can be investigated in studies with smaller numbers of subjects, we measured seed-based connectivity to randomly defined seeds from each network in individual subjects and groups of 5 and 10 subjects. We found that on average 748 ± 334 voxels were correlated to the seed in each subject, of which 42% on average belonged to the white-matter network in which the seed was located (compared to the expected 8% given the null hypothesis that seed-based correlation was random). When averaging across 5 or 10 subjects the number of correlated voxels belonging to the network rose to 59% and 64%, respectively (examples of successful seed-based analyses are shown in Fig. 9). These data suggest that the seed-based approach can be used to investigate white-matter networks in studies with smaller subject cohorts.
Seed-based analysis in a small group identifies the functional networks. Examples of successful seed-based analysis in a group of five subjects. Seeds are marked in blue and correlated voxels in red. The analysis identified extensive networks around the seed with symmetric components in the contralateral hemisphere.
Discussion
Our investigation of functional connectivity in the white-matter of a large group of subjects revealed that it is composed of distinct symmetrical functional networks with commissural components between the hemispheres. These functional networks are related to specific, anatomically established white-matter tracts. White-matter networks were found to be organized in three layers. Superficial white-matter networks are generally correlated with established gray-matter networks, whereas deep networks are relatively independent of the gray-matter networks. Notably, these findings were replicated in a large independent group of subjects and remained relatively similar with different acquisition and processing parameters, and therefore cannot be attributed to noise or artifacts. Seed-based analysis in individual subjects and areal boundary mapping supported the clustering into distinct functional networks and may provide additional ways for their investigation. We provide the identified white-matter functional atlases and software for identification of white-matter functional networks in http://mind.huji.ac.il/white-matter.aspx.
Previous studies have already noted fMRI activity and connectivity inside the white matter. Task-related activations were reported in the corpus callosum and internal capsule (Mazerolle et al., 2008; Fabri et al., 2011; Gawryluk et al., 2011, 2014; Fabri and Polonara, 2013). Further investigations demonstrated correlated activity along white-matter tracts during rest and prolonged task performance in small groups of subjects (Ding et al., 2013, 2016; Marussich et al., 2017). Our findings replicate and extend these results as symmetrical white-matter networks were found to exist throughout the brain, and to correspond to combinations of known white-matter tracts. RSfMRI has been extensively used to identify large-scale networks or regions in the gray matter (Buckner et al., 2011; Power et al., 2011; Yeo et al., 2011; Craddock et al., 2012), and here we demonstrate that similar networks exist within the white matter. In view of these findings, future works should further characterize the roles of the identified white-matter networks during task performance. Importantly, connectivity disturbances in gray-matter resting-state is frequently used as a biomarker for neuropsychiatric disorders (Buckner et al., 2005; Biswal et al., 2010; Zhang and Raichle, 2010; Salomon et al., 2011; Bick et al., 2012; Finke et al., 2013; Peer et al., 2014; Reuveni et al., 2016). Investigation of differences in white-matter resting-state correlations between health and disease can add an additional layer of information in future studies.
An important issue is whether the observed white-matter signals reflect neuronal-related activity in the white matter. Concerns may relate to three issues: random noise artifacts (such as head motion), nonrandom factors affecting the signals such as differential distribution of blood vessels across the white-matter, and contribution of gray-matter signals to the observed networks due to smoothing and partial-volume effects. The possibility that the observed correlations arise from random noise can be contrasted by the fact that networks are replicated with high fidelity across independent subject groups, as well as their highly symmetric character and the existence of fine interhemispheric connecting bridges. The correspondence of signals from white-and gray-matter networks, which are known to be related to task-related activity (Smith et al., 2009; Yeo et al., 2011; Cole et al., 2014; Tavor et al., 2016), further suggests that the white-matter signals correspond to functional activity. Previous findings linking white-matter activity and connectivity to task-related changes also strengthen this conclusion (Gawryluk et al., 2011, 2014; Fabri and Polonara, 2013; Ding et al., 2016; Marussich et al., 2017). Finally, we took several measures to ensure that gray-matter signals do not interfere with white-matter signals, by smoothing the white-matter and gray-matter separately, and using from each subject only voxels which were identified as white-matter in this specific subject. These measures, together with the fact that all identified networks extend deeper than the white-gray matter border, ensure that these signals arise from activity within the white-matter itself.
Most of the current knowledge of white-matter tracts is obtained from structural measures such as diffusion imaging. We report here a novel method for in vivo observation of the functional white matter in humans. Our findings extend the existing anatomical knowledge of white-matter tracts, as they show functional networks composed of several tracts acting in concert, similarly to the discovery of interacting groups of gray-matter regions during the resting state. Furthermore, our results identified superficial white-matter systems, locally connecting adjacent gray-matter regions, which are harder to describe using DTI as they combine many small fibers with different orientations (Oishi et al., 2008). Finally, whole-brain functional investigation enables to directly link white-matter networks to their corresponding gray-matter sources and destinations.
The identification of white-matter functional networks may further complement our understanding of gray-matter functional networks, and provide clues as to how these fragmented networks are connected. Past studies attempted to identify white-matter tracts connecting disconnected nodes of cortical networks (Damoiseaux and Greicius, 2009; Greicius et al., 2009; van den Heuvel et al., 2009). Here we directly measured the functional correlations between white- and gray-matter voxels to identify the linking white-matter tracts, including the commissural fibers. In addition, we identified homogenous functional networks in the deep white matter, which are less correlated in their activity to specific gray-matter networks. The deep white-matter networks also exhibit a different spectral profile compared with surface networks, which is less similar to the gray matter. It is plausible that these deep networks may provide the means of communication between different gray-matter networks.
Despite the large amount of research in the field of RSfMRI, findings of resting-state correlations inside the white matter have only recently been reported. This may be related to several factors: first, the white matter contains a much smaller number of blood vessels (Gawryluk et al., 2014), resulting in low signal-to-noise ratio and weaker correlations within the white matter than those within gray-matter networks. Additionally, white-matter tracts intersect with each other (Alexander et al., 2001, 2002; Tuch et al., 2002; Tournier et al., 2004), possibly lowering the signal correlations even more. Moreover, the white-matter is often overlooked during functional MRI analyses, either implicitly or explicitly by application of gray-matter masking (Yarkoni et al., 2009; Gawryluk et al., 2014). Finally, white-matter signals are often explicitly removed during RSfMRI preprocessing by their use as regressors of no interest (Behzadi et al., 2007; Power et al., 2012; Satterthwaite et al., 2013). In light of our findings, this practice may be problematic because the white matter contains a multitude of signals correlated to those of gray-matter networks.
BOLD signals have not been thoroughly characterized in the white matter. As these signals depend on cerebral blood flow and volume, they are generally lower in magnitude than gray-matter signals but have been reported to be of similar shape in response to stimuli (Fraser et al., 2012). Moreover, white-matter signals have a similar frequency range as gray-matter signals (as we have also observed; Fig. 8), occur ∼TE = T2* (Ding et al., 2016), and are influenced by anesthesia (Wu et al., 2016), further supporting their functional origin. The precise source of the white-matter fMRI signals is unclear, as BOLD signals are thought to arise mainly from postsynaptic potentials (Logothetis et al., 2001), which are largely absent in the white-matter. Gawryluk et al. (2014) discuss two potential sources for white-matter BOLD signals: spiking-related metabolic demands (Smith et al., 2002) and activity of astrocytes and NO-producing neurons (Petzold and Murthy, 2011; Barbaresi et al., 2014). Moreover, as these authors pointed out, previous investigations of the origin of BOLD signals were performed exclusively in the gray matter, and therefore further investigations in the white matter may shed light on the source of this activity.
Despite our systematic findings, this study is not free of limitations. First, in contrast to gray-matter, white-matter tracts can cross each other, thereby mixing the signal from different functional systems in the same locations (Alexander et al., 2001, 2002; Tuch et al., 2002; Tournier et al., 2004). As the K-means clustering algorithm we used is a winner-takes-all approach, it can only assign one label to each brain voxel, therefore possibly neglecting fiber crossings and missing several tract continuations. Second, the AFQ algorithm we used to identify white-matter tracts from DTI data identifies only 20 major tracts, which span approximately half of the white-matter; future works may provide more accurate estimates of the overlap between functional networks and anatomical tracts across the whole brain. Third, we found that preprocessing and acquisition parameters affect the clustering, creating finer separation to networks in specific brain parts and coarser separation in others; further research is necessary to characterize the optimal acquisition and processing parameters for white-matter networks investigation. Finally, our study involved a large group of subjects; the clustering approach was found to be less reliable for functional network identification in smaller groups (N < 20) as well as individual subjects. Nevertheless, works in individual subjects may be done by averaging signals within networks using an existing atlas of functional white-matter networks (as we provide here). Seed-based analysis provides an additional approach for investigation of these signals, setting the stage for application in studies with smaller subject cohorts and individual patients.
To conclude, we have shown here that distinct functional networks in the human white matter can be identified in vivo using RSfMRI correlations, and are related to both gray-matter functional networks and structural white-matter tracts. Our findings were highly replicable across large groups of subjects, demonstrating that they cannot be attributed to noise sources. These white-matter networks add a new dimension to research of RSfMRI, and may inform investigations of white-matter function and structure and provide a missing link for understanding communication within and across whole-brain functional networks. We provide the resulting atlases and software for clustering at http://mind.huji.ac.il/white-matter.aspx; these can be used to identify white-matter functional networks in existing and newly recorded datasets of resting-state and task fMRI. Future investigations may uncover the functional roles of each of these networks, and their function in health and disease.
Notes
Supplemental material for this article is available at http://mind.huji.ac.il/white-matter.aspx, including white-matter functional network atlas, and codes for white-matter functional analyses. This material has not been peer reviewed.
Footnotes
This work was supported by the Israel Science Foundation, the Ministry of Science and Technology of Israel (Grant 3-10789 to M.P.), and the Azrieli Foundation for the Azrieli Fellowship to M.N. We thank Roey Schurr and Rachel Fried for helpful discussions.
The authors declare no competing financial interests.
- Correspondence should be addressed to Michael Peer, Computational Neuropsychiatry Lab, Department of Neurology, Hadassah Hebrew University Medical Center, Jerusalem 91120, Israel. michael.peer{at}mail.huji.ac.il
