Abstract
The thalamus is globally connected with distributed cortical regions, yet the functional significance of this extensive thalamocortical connectivity remains largely unknown. By performing graph-theoretic analyses on thalamocortical functional connectivity data collected from human participants, we found that most thalamic subdivisions display network properties that are capable of integrating multimodal information across diverse cortical functional networks. From a meta-analysis of a large dataset of functional brain-imaging experiments, we further found that the thalamus is involved in multiple cognitive functions. Finally, we found that focal thalamic lesions in humans have widespread distal effects, disrupting the modular organization of cortical functional networks. This converging evidence suggests that the human thalamus is a critical hub region that could integrate diverse information being processed throughout the cerebral cortex as well as maintain the modular structure of cortical functional networks.
SIGNIFICANCE STATEMENT The thalamus is traditionally viewed as a passive relay station of information from sensory organs or subcortical structures to the cortex. However, the thalamus has extensive connections with the entire cerebral cortex, which can also serve to integrate information processing between cortical regions. In this study, we demonstrate that multiple thalamic subdivisions display network properties that are capable of integrating information across multiple functional brain networks. Moreover, the thalamus is engaged by tasks requiring multiple cognitive functions. These findings support the idea that the thalamus is involved in integrating information across cortical networks.
Introduction
The mammalian brain can be conceptualized as a thalamocortical system. Every cortical region receives projections from the thalamus and in turn sends outputs to one or multiple thalamic nuclei (Jones, 2001). Thalamocortical projections relay nearly all incoming information to the cortex as well as mediate corticocortical communication (Sherman and Guillery, 2013). Thus, full insight into brain function requires knowledge of the organization and properties of thalamocortical interactions.
The thalamus can be divided into the following two types of nuclei: first-order and higher-order thalamic nuclei (Sherman, 2007). First-order thalamic nuclei, such as the lateral geniculate nucleus (LGN) and the ventral posterior (VP) nuclei, receive inputs from ascending sensory pathways or other subcortical brain regions. In contrast, higher-order thalamic nuclei, such as the mediodorsal and pulvinar nuclei, receive inputs predominately from the cortex. Higher-order thalamic nuclei have both reciprocal and nonreciprocal connections with multiple cortical regions (Giguere and Goldman-Rakic, 1988; Selemon and Goldman-Rakic, 1988; Haber and Knutson, 2010), suggesting that, in addition to relaying external sensory information to the cortex through first-order nuclei, another principle function of the thalamocortical system is to mediate corticocortical information transfer within the cortex (Sherman, 2016).
Graph-theoretic network analysis of resting-state functional MRI (rs-fMRI) data is well suited for exploring the organization and functional properties of the thalamocortical system. Previous functional connectivity analyses of rs-MRI data have demonstrated that the human cerebral cortex can be decomposed into several modular functional networks (Power et al., 2011; Yeo et al., 2011), each of which is potentially involved in executing a discrete set of cognitive functions (Smith et al., 2009; Bertolero et al., 2015). Graph-theoretic measures can be further used to characterize the topological properties of each region (Sporns et al., 2007). For instance, a brain region with many within-network connections has a strong “provincial hub” property (Guimerà and Nunes Amaral, 2005), presumably to promote within-network interactions for executing specialized functions of the network. In contrast, a brain region with many between-network connections has a strong “connector hub” property, presumably to mediate interactions between functional networks. Cortical connector and provincial hubs have distinct contributions to the modular organization of cortical networks (Gratton et al., 2012), and cortical connector hubs have been found to be involved in multiple cognitive tasks (Cole et al., 2013; Warren et al., 2014; Bertolero et al., 2015; Yeo et al., 2015).
The topographic properties of the thalamocortical system are largely unknown. The thalamus and its constituent nuclei have global and extensive connectivity with multiple brain structures and are likely major hub regions for functional brain networks (Cole et al., 2010; Crossley et al., 2013). Thalamic nuclei are traditionally hypothesized to function as modality-specific relays. This view of thalamic function would predict that thalamic nuclei should exhibit strong provincial hub properties to support information communication within relatively encapsulated processing channels (Alexander et al., 1986). However, it is not known whether all thalamic nuclei serve as provincial hubs or whether some nuclei, such as the higher-order nuclei, further support between-network interactions, serving as connector hubs. It has been shown that a single thalamic nucleus could have connections with multiple anatomically segregated cortical regions (Giguere and Goldman-Rakic, 1988; Selemon and Goldman-Rakic, 1988). If these cortical regions innervated by a single thalamic nucleus belong to the same functional network, then this thalamic nucleus should exhibit strong provincial hub properties for supporting modality-selective processes (Fig. 1A). Alternatively, if these cortical regions belong to different functional networks, then this thalamic nucleus is capable of integrative processes mediating interactions across multiple cortical networks and should exhibit strong connector hub properties (Fig. 1B). These hypotheses are not mutually exclusive—the thalamus could contain nuclei that act as “global kinless” hubs (Guimerà and Nunes Amaral, 2005; Guimerà et al., 2007), which are strongly associated with not only one but multiple functional networks involved in both modality-selective and multimodal integrative processes.
The goal of this study was to elucidate the network topological roles of thalamic nuclei in functional brain networks. To measure network properties of thalamocortical functional connectivity, we performed graph-theoretic network analyses on rs-fMRI data collected from healthy human participants. To relate network topology to cognitive functions, we analyzed the task-related activity using a meta-analysis of 10,449 functional neuroimaging experiments from the BrainMap database (Laird et al., 2005; Yeo et al., 2015). Finally, we examined the contribution of thalamic nuclei to cortical network organization by analyzing rs-fMRI data from human patients with focal thalamic lesions.
Materials and Methods
Datasets.
For the main analyses, we analyzed publically available rs-fMRI data from 303 subjects (mean age = 21.28 years, SD = 2.64; age range = 19–27 years; 128 males) that were acquired as part of the Brain Genomics Superstruct dataset (Holmes et al., 2015). For each subject, two runs (6.2 min each) of rs-fMRI data were collected using a gradient-echo echo-planar imaging sequence with the following parameters: relaxation time (TR) = 3000 ms; echo time (TE) = 30 ms; flip angle = 85°; and 3 mm3 isotropic voxels with 47 axial slices. Structural data were acquired using a multi-echo T1-weighted magnetization-prepared rapid acquisition gradient-echo (MPRAGE) sequence (TR = 2200 ms; TE = 1.54 ms for image 1 to 7.01 ms for image 4; flip angle = 7°; 1.2 mm3 isotropic voxel). We replicated our main analyses with another publically available rs-fMRI dataset from 62 healthy adults (mean age = 23.96 years, SD = 5.24; age range = 18–37 years; 27 males) that was acquired as part of the NKI-Rockland sample (Nooner et al., 2012). For each subject, 9 min and 35 s of rs-fMRI data were acquired using a multiband gradient-echo echo-planar imaging sequence with the following parameters: TR = 1400 ms; TE = 30 ms; multiband factor = 4; flip angle = 65°; 2 mm3 isotropic voxels with 64 axial slices. Structural data were acquired using an MPRAGE (TR = 1900 ms; TE = 2.51 ms; flip angle = 9°; 1 mm3 isotropic voxel). For both datasets, subjects were instructed to stay awake and keep their eyes open.
For the lesion analyses, we analyzed rs-fMRI data from four patients with focal thalamic lesions (ages of the four patients: S1 = 81 years; S2 = 50 years; S3 = 55 years; S4 = 84 years; all males; all were scanned at least 6 months after their stroke). Patient S1 has bilateral lesions, and all other patients have unilateral lesions. Two runs of rs-fMRI data were collected (10 min each; TR = 2000 ms; TE = 30 ms; flip angle = 72°; 3.5 mm2 in plane resolution with 34 axial 4.2 mm slices). Structural images were acquired using an MPRAGE sequence (TR = 2300 ms; TE = 2.98 ms; flip angle = 9°; 1 mm3 voxels). Patients were instructed to stay awake and keep their eyes open. Informed consent was obtained from all patients in accordance with procedures approved by the Committees for Protection of Human Subjects at the University of California, Berkeley.
Functional MRI data preprocessing.
Image preprocessing was performed with the Configurable Pipeline for the Analysis of Connectomics software (Sikka et al., 2014). First, brain images were segmented into white matter (WM), gray matter, and CSF. Rigid body motion correction was then performed to align each volume to a temporally averaged volume, and a boundary-based registration algorithm was used to register the functional data to the anatomical image. Advanced Normalization Tools was used to register the images to the MNI152 template using a nonlinear normalization procedure (Avants et al., 2008). All images were spatially resampled to a 2 mm voxel resolution. We then performed nuisance regression to further reduce non-neural noise and artifacts. To reduce motion-related artifacts, we used the Friston-24 regressors model during nuisance regression (Friston et al., 1996). WM and CSF signals were regressed using the anatomical CompCor approach with five components for each tissue class (Behzadi et al., 2007). Linear and quadratic drifts were also removed. Because the physical proximity between the thalamus and the ventricles could result in the blurring of the fMRI signal, we further regressed out the mean signals from CSF, WM, and gray matter that were within 5 voxels (10 mm) from the thalamus. After regression, data were bandpass filtered from 0.009 to 0.08 Hz and scaled to a whole-brain mean value of 10,000. Given that the thalamus is a relatively small structure, to avoid signal blurring we did not perform any spatial smoothing.
Identifying cortical functional networks.
Following preprocessing, mean rs-fMRI time series were extracted from 333 cortical regions of interest (ROIs; Gordon et al., 2016) and concatenated across runs for subjects with multiple rs-fMRI scans. Cortico-cortical functional connectivity was assessed in each subject by computing Pearson correlations between all pairs of cortical ROIs, resulting in a 333 × 333 correlation matrix. This correlation matrix was then thresholded to retain the strongest functional connections. To identify the modular structure of cortico-cortical functional connectivity, we performed an InfoMap algorithm (Rosvall and Bergstrom, 2007) to partition the correlation matrix into putative modular functional networks. InfoMap is one of the best performing subnetwork detection algorithms available and has been successfully used to identify cortical network organization (Power et al., 2011; Bertolero et al., 2015). Specifically, for each subject, we first identified functional networks using the InfoMap algorithm at a threshold of density of 0.15 (keeping 15% of connections of all possible cortico-cortical connections). Based on this partition result, we then constructed a consensus matrix that consisted of values of 1 where the two ROIs are in the same module and values of 0 elsewhere (Lancichinetti and Fortunato, 2012). We then rethresholded the correlation matrix by decreasing the density in steps of 0.001 and ran the InfoMap algorithm again to obtain a new partition. The consensus matrix was then updated for the new partition, except for rows and columns for which the ROI had no connections that survived the new threshold or the ROI was in a fragmented module of fewer than five ROIs. The consensus matrix continued to be updated to the minimal threshold of density (0.01). Thus, this subject-specific consensus matrix represents the modular network assignment for each pair of nodes at their sparsest level possible (i.e., before they become disconnected from the graph). We then aggregated individual subjects' module organization by averaging consensus matrices across subjects. This averaged consensus matrix was then submitted to the same recursive method described above to identify group-level cortical functional networks. Following methods reported in previously published studies (Power et al., 2011; Gordon et al., 2016), networks with five or fewer ROIs were eliminated from further analyses, and ROIs that were not clustered into networks were excluded from further analyses.
Thalamus parcellation.
Previous studies of functional brain networks using graph theory often excluded subcortical structures or examined the thalamus with gross or no subdivisions. Given the complex structure of the thalamus with multiple distinct nuclei, the thalamus is likely not uniformly interacting with the cortex. Instead, different thalamic subdivisions have distinct structural connectivity with the cortex and thus functional connectivity with the cortex (Behrens et al., 2003; Arcaro et al., 2015; Yuan et al., 2016; Kumar et al., 2017). We therefore examined thalamocortical network properties of using three different thalamic parcellations. To localize the thalamus, the Morel Atlas (Krauth et al., 2010) was used to define its spatial extent (2227 2 mm3 voxels included in the atlas, registered to the MNI152 template). To identify thalamic subdivisions, three different thalamic atlases were used. We first performed a custom winner-take-all functional parcellation using rs-fMRI data. We calculated partial correlations between the mean BOLD signal of each cortical functional network and the signal in each thalamic voxel, while removing signal variance from other functional networks. Partial correlations were then averaged across subjects, and each thalamic voxel was labeled according to the cortical network with the highest partial correlation. The Morel atlas identified thalamic nuclei based on cyto- and myelo-architecture in stained slices of postmortem tissue collected from five postmortem brains (Morel et al., 1997) and subsequently transformed to MNI space (Krauth et al., 2010). The Oxford-FSL thalamic structural connectivity atlas defined thalamic subdivisions based on its structural connectivity with different cortical regions estimated from diffusion imaging data (Behrens et al., 2003). Each of these atlases is sensitive to a specific type of anatomical or functional information; therefore, using all three atlases allows us to derive a more complete characterization of thalamocortical network properties.
Thalamic and cortical nodal properties.
To formally quantify the network properties of thalamocortical functional connectivity for each subject, we first extracted the signal from the thalamus by using either the preprocessed signal of each voxel or the averaged voxelwise BOLD signal within each thalamic subdivision. We then calculated the partial correlation between the mean BOLD signal of each cortical ROI and thalamic voxel or subdivision. Partial correlation was calculated by removing signal variance from all other cortical ROIs. Given the large number of cortical ROIs, a dimension reduction procedure using principal component analysis was performed based on signals from cortical ROIs not included in the partial correlation calculation, and eigenvectors that explained 95% of the variance were entered as additional nuisance regressors in the model. Note that no correlations were calculated between thalamic voxels. This resulted in two thalamocortical connectivity matrices: voxelwise and subdivision-wise. For voxelwise matrices, 2227 thalamic voxels were included. For the functional parcellation, Oxford-FSL and Morel atlases, 9, 7, and 15 thalamic subdivisions were included, respectively. Note that each thalamic voxel and thalamic subdivision all had the same number of functional connections with cortical ROIs (total = 333). Thalamocortical connectivity matrices were then averaged across subjects. For cortical ROIs, we calculated Pearson correlations between all pairwise cortical ROIs to obtain cortico-cortical connectivity matrices (matrix size, 333 × 333). Matrices were then averaged across subjects. Following recommendations from a previous study (Power et al., 2011), we explored network properties across a range of network density thresholds (density = 0.01–0.15) and submitted these to graph analyses. Note that, because we did not perform global signal regression, no negative correlations entered our graph analyses.
For estimating thalamocortical functional connectivity, we chose partial correlations over full correlations because past studies have shown that detailed thalamocortical connectivity patterns could be obscured without accounting for the shared variance between cortical regions (Zhang et al., 2008). We compared graph matrices of thalamic voxels calculated using partial or full correlations and found that results were independent of whether full or partial correlations were used. We also calculated graph metrics at the level of each thalamic subdivision by averaging voxelwise signals within each thalamic subdivision before graph analyses and found consistent results when compared with voxelwise metrics.
For each thalamic voxel, subdivision, and cortical ROI, we calculated the participation coefficient (PC), a measure of the connector hub property, and the within-module degree z-score (WMD), a measure of the provincial hub property (Guimerà and Nunes Amaral, 2005). To ensure that our results were not biased by a single specific threshold, all graph metrics were calculated across a range of thresholds (density = 0.1–0.15). We report results that were averaged across thresholds.
To calculate the WMD, correlation matrices were first binarized by setting weights above the density threshold to 1. WMD values were calculated across a range of density thresholds (edge density, 0.1–0.15) and averaged across thresholds. Weights were binarized to equate the connectivity weights between thalamocortical and corticocortical networks. WMD is calculated as WMD = , where C̅W̅s̅ is the average number of connections among all cortical ROIs within cortical network s, and σCWs is the SD of the number of connections of all ROIs in network s. Ki is the number of connections between i and all cortical ROIs in network s. Because our goal was to understand the contribution of the thalamus to cortical network organization, WMD scores of thalamic subdivisions were calculated using the mean and SD of the within-network degree (number of intra-network connections) calculated from each cortical functional network. For comparison, we also calculated WMD using fully weighted matrices without thresholding, and found similar results. Note that because WMD is normalized by the number of regions within the associated functional network, unlike degree-based centrality measures, it is not influenced by the size of functional networks. Given that we assigned each thalamic voxel to a cortical functional network in the functional parcellation atlas, higher WMD values reflect more within-network connections of the thalamic voxel with the network it was assigned to.
The PC value of region i is defined as follows: where Ki is the sum of the connectivity weight of i, Kis is the sum of the connectivity weight between i and the cortical network s, and NM is the total number of networks. If a region has connections uniformly distributed to all cortical networks, then its PC value will be close to 1; otherwise, if its connectivity is concentrated within a specific cortical network, its PC value will be close to 0. Given that we identified nine cortical functional networks, the maximum possible PC value would be 0.89. We therefore further divided PC values by this theoretical upper limit so that the highest possible PC value given the network architecture would be 1. For comparison, we also calculated the PC using binary networks by setting weights above the threshold to 1 and found similar results. Because PC is independent of the node degree, it is less biased by the number of ROIs within each functional network, resulting in a superior measure of hub properties (Power et al., 2013). It is also important to note that a PC value is independent of the functional network that an ROI belongs to; therefore, our thalamic PC calculations are not influenced by which cortical network thalamic voxels or subdivisions are assigned to it.
Patients with focal thalamic lesions.
Lesion masks were manually traced in the native space according to visible damage on a T1-weighted anatomical scan and further guided by hyperintensities on a T2-weighted fluid-attenuated inversion recovery image. Lesion masks were then warped into the MNI space using the same nonlinear registration parameters calculated during preprocessing. For comparing the effects of local thalamic lesions, rs-fMRI volumes with framewise displacement that exceeded 0.5 mm were removed from further analysis (scrubbed) before bandpass filtering. More volumes were scrubbed for patients when compared with healthy control subjects (mean percentage of frames scrubbed for patients = 11.04%, SD = 10.63%; mean percentage of frames scrubbed for healthy control subjects = 1.02%, SD = 1.62%; no subjects were excluded).
Modularity.
Modularity can be measured by Newman's modularity Q (Newman, 2006), which was defined as follows: where eii is the fraction of the connectivity weight connecting ROIs within a cortical functional network i, ai is the fraction of the connectivity weight connecting ROIs in cortical functional network i to other cortical networks, and m is the total number of cortical functional networks. The network partition from Figure 2A was used. Modularity was calculated for the whole brain (including all cortical ROIs). Note that no thalamic subdivision was included into this analysis. For lesions analyses, the mean and SD of Q calculated from all healthy control subjects were used to normalize Q values from each patient (Gratton et al., 2012), as follows: where Q is the patient's modularity score, Q̄ is the average modularity score of healthy control subjects, and δ is the SD in modularity of healthy control subjects. Modularity calculations were performed separately for each density threshold (density = 0.01–0.15), and results are presented across thresholds.
Within-network and between-network connectivity strength.
We further analyzed how focal thalamic lesions could affect connectivity between and within cortical functional networks. For each patient, we first mapped the cortical functional networks their thalamic lesions were connected to according to the functional parcellation atlas from healthy subjects (see Thalamic parcellation section). For each patient and healthy control subject, we further normalized each functional connection by subtracting the mean and divided by the SD of all connections within each subject. For the affected cortical networks, we then summed its total connectivity strength with other cortical networks (between-network) and with its own respective networks (within-network). The mean and SD of connectivity strength for healthy control subjects were then used to normalize the patient's values to z-scores.
Meta-analysis of functional neuroimaging experiments in the BrainMap database.
We reanalyzed data presented in a previously published meta-analysis of the BrainMap database (Yeo et al., 2015). In the meta-analysis, a hierarchical Bayesian model was used to derive a set of 12 cognitive components that best describe the relationship between behavioral tasks and patterns of brain activity in the BrainMap database (Laird et al., 2005). Specifically, each behavioral task (e.g., Stroop, stop-signal task, finger tapping) engages multiple cognitive components, and in turn each cognitive component is supported by a distributed set of brain regions. To determine whether or not a thalamic voxel is recruited by a cognitive component, a threshold of p = 10−5 was used. This is an arbitrary yet stringent threshold that was used in two prior studies (Bertolero et al., 2015; Yeo et al., 2015). Critically, there is potential spatial overlap between components. Therefore, brain regions that can flexibly participate in multiple cognitive components could be identified by calculating the number of cognitive components that each brain region engages. The number of cognitive components was summed for each voxel and cortical ROI and was defined as a “cognitive flexibility” score (Yeo et al., 2015).
Results
Identification of cortical networks
To identify cortical functional networks, we first measured functional connectivity matrices among 333 cortical ROIs (Gordon et al., 2016) then performed a network partition analysis to estimate cortical network organization (see Materials and Methods). We found that the cerebral cortex can be decomposed into nine functional networks (Fig. 2A), an organization scheme largely similar to those in previous studies (Power et al., 2011; Yeo et al., 2011; Bertolero et al., 2015; Gordon et al., 2016).
Parcellation of the thalamus
Given that the thalamus can be subdivided using different approaches, we performed our analyses using three different atlases based on data from rs-fMRI, diffusion tensor imaging (DTI), and postmortem histology (Fig. 2A–C; for details, see Materials and Methods). Using resting-state functional connectivity MRI data, we identified thalamic subdivisions that demonstrated the strongest functional connectivity with the different cortical functional networks reported above (Fig. 2A; henceforth referred to as the functional parcellation atlas). We further replicated these results with an independent dataset and found high correspondence among datasets (Fig. 2A; normalized mutual information = 0.69; z-scored Rand coefficient = 144.13; p < 10−5). The Oxford-FSL thalamocortical structural connectivity atlas (Fig. 2B) subdivides the thalamus based on structural connectivity (estimated using probabilistic diffusion tractography on DTI data) into the following seven large cortical areas: primary motor, primary somatosensory, occipital, premotor, prefrontal (including medial and orbitofrontal cortices), parietal, and temporal cortices (Behrens et al., 2003). The Morel atlas (Fig. 2C) subdivides the thalamus into smaller nuclei based on cyto- and myelo-architecture information from five postmortem brains (Morel et al., 1997; Krauth et al., 2010). We further classified each thalamic nucleus from the Morel atlas into first-order or higher-order thalamic nuclei (Sherman and Guillery, 2013; Sherman, 2016). Our parcellation results and numerous replicable human MRI-based parcellations (Behrens et al., 2003; Johansen-Berg et al., 2005; Zhang et al., 2008; Kim et al., 2013; Arcaro et al., 2015; Kumar et al., 2015, 2017; Ji et al., 2016; Yuan et al., 2016) suggest that distinct signals from different thalamic subdivisions can be reliably detected with a conventional 2–4 mm voxel resolution.
Network properties of thalamocortical functional connectivity
To determine the network property of each thalamic subdivision, we estimated functional connectivity between each thalamic voxel and every cortical ROI (see Materials and Methods) to generate a thalamic voxel-to-cortical ROI thalamocortical graph. Graph metrics were calculated for every thalamic voxel and were pooled across voxels for different categories of thalamic subdivisions (i.e., first-order and higher-order nuclei or subdivisions within the three different thalamic atlases). Our goal was to examine the topological properties of the thalamus network in the context of functional brain networks. For example, if a thalamic subdivision is found to have strong connector hub properties, how does it compare to cortical ROIs that are connector hubs for cortical functional network? Therefore, we further calculated graph metrics for each cortical ROI by analyzing patterns of corticocortical functional connectivity.
Within-module degree analyses
Provincial hub properties can be measured by the graph theory metric WMD, which measures the number of within-network connections of a region, z-scored by the mean and SD of within-network connections of all regions in that network (Guimerà and Nunes Amaral, 2005). Higher values reflect more within-network connections. We found that thalamic voxels in both first-order and higher-order thalamic nuclei exhibited higher WMD values than most cortical ROIs (Fig. 3A). A two-sample Kolmogorov–Smirnov test showed that on average thalamic voxels had significantly higher WMD values when compared with cortical ROIs (mean WMD for thalamic voxels = 1.62, SD = 1.9; mean WMD for cortical ROIs = 0, SD = −0.89; D(333,2227) = 0.45, p < 10−10). We further compared the WMD values of thalamic voxels with those of cortical provincial hubs. Because both cortical and thalamic WMD values exhibited a unimodal distribution, we arbitrarily defined cortical provincial hubs as cortical ROIs with WMD values >90% of all cortical ROIs (threshold: WMD = 1.04). We found that on average thalamic voxels exhibited WMD values that were comparable to those of cortical provincial hubs (mean WMD for cortical provincial hubs = 1.35, SD = 0.26; mean WMD for voxels within first-order thalamic nuclei = 1.37, SD = 1.79; mean WMD for voxels within higher-order thalamic nuclei mean WMD = 1.75, SD = 1.94).
Participation coefficient analyses
Connector hub properties can be measured by the graph theory metric called PC, which is a measure of the strength of inter-network connectivity for each region normalized by their expected value (Guimerà and Nunes Amaral, 2005). Higher values reflect more inter-network connections. For each thalamic voxel, we calculated its PC value based on its inter-network connectivity pattern to all cortical ROIs. We found that thalamic voxels in both first-order and higher-order thalamic nuclei exhibited higher PC values than most cortical ROIs (Fig. 3C). A two-sample Kolmogorov–Smirnov test showed that thalamic voxels had significantly higher PC values compared with cortical ROIs (mean PC for thalamic voxels = 0.76, SD = 0.12; mean PC for cortical ROIs = 0.36, SD = 0.22; D(333,2227) = 0.65, p < 10−10). We further compared the PC values of thalamic voxels to those of cortical connector hubs. Because PC values also exhibited a unimodal distribution, we arbitrarily defined cortical connector hubs as cortical ROIs with PC values >90% of all cortical ROIs (threshold: PC = 0.64). We found that, on average, voxels in both first-order and higher-order thalamic nuclei exhibited high PC values that were comparable to those of cortical connector hubs (mean PC for cortical connector hubs = 0.72, SD = 0.05; mean PC for voxels in first-order thalamic nuclei = 0.74, SD = 0.13, mean PC for voxels in higher-order thalamic nuclei = 0.77, SD = 0.11).
Spatial distribution of PC and WMD values
We found that cortical ROIs in the precuneus, and medial frontal, inferior parietal, insular, and middle frontal cortices exhibited high WMD values, whereas ROIs in anterior and inferior frontal, superior precentral sulcus, intraparietal sulcus, and lateral occipital cortices exhibited high PC values (Fig. 4A). High PC and WMD values were found throughout the thalamus (Fig. 4B,C). To determine the differences in the spatial distribution of connector and provincial hub properties in the thalamus, we identified thalamic voxels that exhibited WMD or PC values greater than the cortical connector and provincial hubs (cortical provincial hub threshold WMD = 1.04; cortical connector hub threshold PC = 0.64; hubs were arbitrarily defined as cortical ROIs with PC or WMD values >90% of all cortical ROIs). We found that anterior, medial, posterior, and dorsal parts of the thalamus exhibited both high PC and WMD values. In addition, portions of the lateral thalamus also exhibited only high PC values, and a small portion of the ventral thalamus exhibited only strong high WMD values (Fig. 4D).
PC and WMD values of each thalamic subdivision
We calculated the median WMD and PC values across voxels for each thalamic subdivision and compared those values to cortical connector and provincial hubs (see definition above). Based on the functional parcellation atlas, thalamic subdivisions that showed dominant functional coupling with cingulo-opercular (CO), default mode (DM), frontoparietal (FP), medial temporal (mT), and superior FP (sFP) networks exhibited high WMD values numerically comparable to those of cortical provincial hubs (Fig. 5A). Based on the Oxford-FSL thalamocortical structural connectivity atlas, thalamic subdivisions with dominant structural connectivity with the prefrontal cortex and temporal cortices showed high WMD values comparable to cortical provincial hubs (Fig. 5B). Based on the Morel histology atlas, thalamic subdivisions with high WMD values comparable to cortical provincial hubs included the anterior nucleus (AN), LGN, the ventral lateral nucleus (VL), intralaminar nuclei (IL), lateral posterior nucleus (LP), medial dorsal nucleus (MD), medial pulvinar (PuM), and ventral anterior nucleus (VA; Fig. 5C). We found that all thalamic subdivisions exhibited high PC values that were comparable to or higher than those of cortical connector hubs (Fig. 5D–F).
Connectivity patterns of specific thalamic nuclei
Based on the Morel histology atlas, we found that AN, LGN, VL, IL, LP, MD, PuM, and VA exhibited both high PC and WMD values comparable to those of cortical connector and provincial hubs. To further probe their connectivity patterns, for each nucleus we calculated the mean functional connectivity strength with each of the nine cortical functional networks (using partial correlations; see Materials and Methods), and divided by the summated total connectivity strength of the nucleus with all networks. If a nucleus is diffusely interacting with all functional networks, then it should devote ∼11% (1/9 = 0.11) of its total connectivity for each network. In contrast, if a nucleus interacts only with a selective network, the majority of its connectivity strength should be devoted to that network, while connectivity with other networks should be considerably lower. Consistent with the high PC values we observed in these nuclei, we found that each of these thalamic nuclei exhibited a diffuse functional connectivity pattern, with strong connectivity (>11% of its total connectivity strength) with multiple cortical functional networks (Fig. 6). Specifically, we found that multiple nuclei showed strong connectivity with both the CO and DM networks (AN, MD, VA, IL, LP, and VL) and that most of these nuclei had strong connectivity with at least one other cortical functional network (e.g., AN, VA, LP, VL). Altogether, every nucleus had strong functional connectivity with at least three or more cortical functional networks, and many cortical functional networks had strong functional connectivity with multiple thalamic nuclei.
Meta-analysis of the BrainMap database
We found that multiple thalamic subdivisions exhibited both strong WMD and PC values, suggesting that the thalamus is capable of mediating information communication both within and between multiple functional brain networks. Given that individual cortical functional networks are putatively associated with distinct cognitive functions (Smith et al., 2009; Bertolero et al., 2015), it is likely that any individual thalamic nucleus with distributive connectivity with multiple cortical functional networks is involved in multiple cognitive functions. We tested this hypothesis by analyzing the results from a published meta-analysis of 10,449 functional neuroimaging studies (Yeo et al., 2015). This published meta-analysis derived latent variables—an ontology of cognitive functions or “cognitive components”—that best described the relationship between 83 behavioral tasks and the corresponding brain activity. From these data, a “cognitive flexibility score” can be estimated by summing the number of cognitive components that are engaged by every brain region (Yeo et al., 2015). Here, a brain region with a high cognitive flexibility score is assumed to be involved in more cognitive functions. We have previously demonstrated that cortical connector hubs exhibit high cognitive flexibility scores (Bertolero et al., 2015). We hypothesize that thalamic subdivisions that are recruited by multiple cognitive components likely serve as integrative connector hubs, whereas thalamic subdivisions that are recruited by a limited number of specific functions are likely involved in domain general processes.
Consistent with our hypothesis, both first-order and higher-order thalamic nuclei were found to be involved in multiple cognitive components (Fig. 7A). A two-sample Kolmogorov–Smirnov test showed that thalamic voxels have higher cognitive flexibility scores when compared with cortical ROIs (mean for thalamic voxels = 3.67, SD = 1.71; mean for cortical ROIs = 0.36, SD = 0.22; D(333,2227) = 0.54, p < 10−10). We further examined the specific cognitive components (C1–C12; for details, see Yeo et al., 2015) that recruited each of the thalamic subdivisions. As an example, VL, with projections to motor and premotor cortices (Alexander et al., 1986), is recruited by components C1 and C2 that predominately recruit motor cortices (Fig. 7B). However, VL also participates in other cognitive components that recruit lateral prefrontal, medial prefrontal, and parietal cortices (C8, C9, and C12; Fig. 7B). Note that most of these cognitive components not only recruit VL, but also engage other parts of the thalamus (Fig. 7B). This observation is consistent with the results we presented in Figure 6, showing that multiple thalamic subdivisions could be associated with the same cognitive component.
Thalamic lesions have global and distal effects on cortical network organization
Modularity is a metric that quantifies the extent to which the brain is differentiated into separable subnetworks and is an essential property found in many complex systems (Sporns and Betzel, 2016). Previous studies have suggested that connector and provincial hubs have distinct contributions to modular organization. For example, a lesion study showed that damage to connector hubs, but not to provincial hubs, causes more severe disruption of the modular organization of the network (Gratton et al., 2012), suggesting that focal lesions to connector hubs can have a widespread impact on network organization when between-network connections are disrupted. Based on these findings, we predict that if thalamic subdivisions serve as connector hubs for functional brain networks, lesions to those subdivisions should reduce the modular organization of these networks. Modularity can be measured by Newman's modularity Q (Newman, 2006), a comparison between the number of connections within a module to the number of connections among modules. In four patients with focal thalamic lesions [Fig. 8A; one patient (S1) has bilateral lesions, and three patients (S2–S4) have unilateral lesions], we examined the effect of a thalamic lesion on cortical modularity across the whole cerebral cortex (Fig. 8B). Each patient's Q score was converted to a z-score using the mean and SD of healthy control subjects (see Materials and Methods). In all four patients, whole-brain modularity was lower (as indicated by negative z-scores).
Reduction in modularity could be caused by increased between-network connectivity or decreased within-network connectivity. Thus, we further quantified how focal thalamic lesions affect between- and within-network connectivity strength. For each patient, we first mapped the lesioned thalamic voxels, and then we identified cortical networks that are connected to these lesioned voxels using data from healthy control subjects (Figs. 2A, 8A, right). For the connected networks of each lesion, we calculated the change of between- and within-connectivity strength relative to those of healthy control subjects (Fig. 8C; see Materials and Methods for the z-scoring procedure). In all four patients with thalamic lesions, the between-network connectivity strength was higher than that of healthy control subjects (as indicated by positive z-scores); in three patients, within-network connectivity strength was lower than that of healthy control subjects (as indicated by negative z-scores). Overall, increased between-network connectivity strength (mean z-score = 2.79, SD = 1.73) was more prominent than decreased within-network connectivity (mean z-score = −0.71, SD = 1.10), suggesting that focal thalamic lesions increase between-network connectivity.
Secondary analyses
Replication
We replicated the WMD and PC analyses using an independent rs-fMRI dataset with higher native voxel resolution (2 mm3); the spatial correlation values across both cortical ROIs (333 ROIs) and thalamic voxels (2227 voxels) between the test and replication datasets for PC and WMD scores were 0.74 and 0.78, respectively. We also replicated our results using a different cortical ROI definition template that consists of 320 cortical ROIs (Craddock et al., 2012), and the thalamic voxelwise spatial correlation values for PC and WMD scores were 0.63 and 0.78, respectively.
InfoMap versus spectral modularity detection
To ensure that our lesion analysis was not specific to calculating Q with a modular partition derived using the InfoMap algorithm, we recalculated Q values for each healthy control subject using the spectral modularity detection algorithm that maximizes Q values to derive a modular partition (Newman, 2006), and found it to be highly correlated with and slightly lower than Q values calculated using the InfoMap algorithm (at density = 0.05, r = 0.86; mean Q for InfoMap algorithm = 0.57, SD = 0.06; mean Q for spectral algorithm = 0.53, SD = 0.06).
Discussion
In this study, we provide evidence suggesting that the human thalamus is a critical integrative hub for functional brain networks. First, we found that all thalamic subdivisions exhibited strong between-network connectivity, indicating that a single thalamic subdivision connects not only with multiple cortical regions, but also with multiple cortical functional networks. From a meta-analysis of 10,449 neuroimaging experiments, we further found that the thalamus is engaged by multiple cognitive functions, supporting a role in multimodal information processing. Finally, we found that focal thalamic lesions cause a disruption of the modular structure of cortical functional networks, further underscoring the critical contribution of thalamic function to brain network organization.
The human brain is composed of modular functional networks (Crossley et al., 2013; Bertolero et al., 2015), which comprise provincial hubs—brain regions important for within-network communication—and connector hubs—brain regions important for communication between networks. Here, we used graph-theoretic measures to estimate provincial and connector hub properties of the thalamus. Consistent with traditional interpretations of the relay of information to the cortex from the thalamus, we found that multiple thalamic subdivisions exhibited strong provincial hub properties (high WMD values), supporting modality-selective processes. However, both first-order and higher-order thalamic nuclei also exhibited prominent connector hub properties (high PC values), an indication of their involvement in integrative, multimodal processes (Figs. 3, 4, 5, 6, 7). For example, we found that both first-order (Fig. 6, LGN) and higher-order (Fig. 6, PuM) nuclei showed strong functional connectivity with more than three cortical functional networks. We also found that many thalamic nuclei showed strong connectivity with the CO and DM networks, a finding that is consistent with previous studies that have classified the thalamus as part of the CO, saliency, or DM networks (Dosenbach et al., 2007; Seeley et al., 2007; Power et al., 2011). However, these nuclei also exhibited strong functional connectivity with at least one additional cortical network. Based on these results, both first-order and higher-order thalamic nuclei should be more accurately described as “global kinless” hubs (Guimerà and Nunes Amaral, 2005; Guimerà et al., 2007), with strong connectivity with not only one functional network but homogenously interacting with multiple networks. The global hub property of thalamic nuclei could potentially allow the thalamus to send and access information across diverse cortical functional networks. Via the convergence of information, the thalamus may serve as a global hub that subserves multiple cognitive functions. Although it has previously been proposed that the thalamus is not simply a relay station and that higher-order thalamic nuclei could further serve to mediate cortical-to-cortical communication within its anatomical projection zone (Sherman, 2016), the notion that both the first-order and higher-order nuclei also play an integrative role interacting with multiple cortical functional networks has not been demonstrated in humans. In a rodent fMRI study (Liska et al., 2015), using a measure that was analogous to PC, it was demonstrated that the rodent thalamus has the most diverse connectivity pattern with multiple functional subnetworks, more diverse than the hippocampus, basal forebrain, DM, and lateral cortical networks.
Higher-order thalamic nuclei, which receive inputs predominately from the cortex, are hypothesized to provide transthalamic routes to support cortico-cortical interactions within a functional network that receives its projections (Saalmann et al., 2012; Sherman, 2016). For example, the posterior nucleus transfers information from primary to secondary somatosensory areas (Theyel et al., 2010). Likewise, the pulvinar has extensive reciprocal connections with striate and extrastriate visual cortices (Adams et al., 2000) and is thought to modulate information communication between visual areas (Saalmann and Kastner, 2011; Saalmann et al., 2012). In contrast, first-order thalamic nuclei, which receive projections from peripheral sensory organs or other subcortical structures, have projections to primary cortices and are thought to act as modality-selective relays for a limited type of afferent signal to the cortex. Our graph-theoretic analyses of thalamocortical functional connectivity provide evidence suggesting that both first-order and higher-order thalamic nuclei participate in information exchange with multiple cortical functional networks.
How can an individual thalamic nucleus interact with multiple functional cortical networks? One possibility is that thalamic nuclei could have dense reciprocal connections with cortical connector hubs that in turn are connected with multiple cortical functional networks. For example, from our analyses of the Oxford-FSL and functional atlases, we found that thalamic subdivisions that were most strongly connected with the frontal cortex, the temporal cortex and frontoparietal functional networks (e.g., FP, sFP, CO, DF), showed kinless global hub properties (prominent connector plus provincial hub properties). Previous studies have found that cortical connector hubs are primarily located in frontoparietal and temporal cortices (Power et al., 2013; Bertolero et al., 2015). Our findings suggest that, in the context of functional brain networks, many thalamic subdivisions could share similar information-processing roles with associated heteromodal regions. In addition, both first-order and higher-order nuclei are known receive nonreciprocal corticothalamic innervations from multiple cortical regions (McFarland and Haber, 2002). Higher-order thalamic nuclei are also known to project to more than one area of the cerebral cortex (Jones, 1998). For example, MD projects to the lateral frontal cortex, the insula, the anterior cingulate cortex, the temporal cortex, and the supplementary motor area (Markowitsch et al., 1985; Giguere and Goldman-Rakic, 1988), and these regions likely span multiple different cortical networks (e.g., FP, DF, CO, and SM networks). Also, thalamic nuclei not only comprise “core” thalamocortical cells that project to middle layers of the cerebral cortex that are constrained by the borders of a functional area, but also comprise “matrix” thalamocortical cells that diffusely project to superficial layers of the cerebral cortex, and are unconstrained by the boundaries of cortical topographic representations (Jones, 1998). The inhibitory thalamic reticular nucleus also receives projections from multiple cortical and subcortical regions and further projects to multiple thalamic nuclei (Sherman and Guillery, 2013). Thus, the reticular nucleus could receive information from different cortical networks and modulate activity in both first-order and higher-order thalamic nuclei. To summarize, individual thalamic nucleus receive projections and project to multiple cortical regions that belong to different cortical functional networks. Thus, the anatomical and functional architectures of the thalamus are capable of simultaneously receiving and transmitting signals among multiple brain regions that belong to different cortical functional networks. Our results further suggest a many-to-many relationship of the thalamocortical system. Specifically, a single thalamic subdivision could be functionally connected with multiple cortical functional networks, and a single cortical functional network could in turn be functionally connected with multiple thalamic subdivisions.
Consistent with its extensive connectivity with multiple cortical functional networks, we found that the thalamus is one of the most “cognitively flexible” brain regions, indicating that the thalamus is involved in a diverse range of behavioral tasks. This observation derived from our meta-analysis of the BrainMap database is further supported by several representative empirical studies demonstrating that the thalamus mediates interactions between high-level cognitive processes (e.g., attention and working memory) and more elementary sensorimotor functions (Saalmann et al., 2012; de Bourbon-Teles et al., 2014; Zhou et al., 2016). For example, a monkey electrophysiology study found that deactivating the pulvinar reduced the attentional effects on sensory-driven evoked responses recorded in V4 (Zhou et al., 2016). Also, optogenetically perturbing thalamic activity in rodents impaired the ability of the animals to select between conflicting visual and auditory stimuli (Wimmer et al., 2015). Finally, VL lesions in humans impair their ability to use a memorized cue during working memory necessary for guiding a visual search of multiple visual stimuli (de Bourbon-Teles et al., 2014). Together, results from our graph analyses of thalamocortical functional connectivity and meta-analysis of task-related thalamic activity patterns suggest that the thalamus participates in interactions among multiple functional cortical networks, networks that are putatively involved in distinct cognitive functions. Based on our empirical results, future task-based studies can test this hypothesis regarding the role of the thalamus in integrative functions.
Previous studies (Gratton et al., 2012) have suggested that connector hubs are critical for maintaining the modular architecture of functional brain networks. Mathematically, whole-brain modularity is inversely related to between-network connectivity; therefore, a loss of connector hub should increase modularity. However, our previous work with patients with focal cortical lesions found that damage to cortical connector hubs decreased modularity (Gratton et al., 2012). Furthermore, a transcranial magnetic stimulation study (Gratton et al., 2013) further found that disruption of the connector hub function increased between-network connectivity. This suggests that, in addition to mediating internetwork interaction, connector hubs could be further involved in maintaining the modular structure of functional networks through decreasing between-network connectivity. Consistent with these prior studies, we found that thalamic lesions reduce the modularity of cortical functional networks and increase between-network connectivity strength. This suggests that thalamic connector hub damage weakens the modular organization of cortical networks by increasing the strength of between-network connectivity. Moreover, since the effect of a thalamic lesion is not constrained only to cortical regions it directly projects to, it can be considered an example of “connectomal diaschisis” (Carrera and Tononi, 2014).
One possible interpretation of our findings from patients with thalamic lesions is that through its extensive between-network connectivity with multiple cortical regions, the thalamus acts as a “gate” to inhibit cortico-cortical connectivity while maintaining the modular structure of cortical networks. A weakening of the “gating” function following a thalamic lesion would disinhibit cortico-cortical connectivity, causing an increase in between-network connectivity and a decrease in modularity. Alternatively, if functional brain networks must maintain a stable between-network connectivity strength for optimal between-network interactions, a compensatory mechanism may lead to increased direct cortco-cortical functional connectivity strength in response to decreased cortico-thalamo-cortical connectivity. To pinpoint the underlying mechanism would require techniques that can simultaneously manipulate and record the activities of the thalamocortical system.
Note Added in Proof: During the proof stage, the “replication dataset” in Figure 2A was replaced with an updated parcellation file to correct for a figure error in the Early Release version published 27 April 2017. Figure 2 has now been corrected.
Footnotes
This work was supported by National Institutes of Health (NIH) Grants RO1-NS-79698 and F32-NS-090757. The Morel atlas was obtained by written consent from Professor G. Szekely from the Computer Vision Laboratory of ETH Zurich, Zurich, Switzerland. This project was made possible by a collaborative agreement allowing comprehensive access to the BrainMap database, a copyrighted electronic compilation owned by the University of Texas Board of Regents. BrainMap is supported by the NIH, National Institute of Mental Health Award R01-MH-074457.
The authors declare no competing financial interests.
- Correspondence should be addressed to Dr. Kai Hwang, 132 Barker Hall, MC 3190, University of California Berkeley, California, CA 94720. kai.hwang{at}berkeley.edu