Functional Connectivity-Based Parcellation of the Thalamus: An Unsupervised Clustering Method and Its Validity Investigation

Abstract

Node definition or delineating how the brain is parcellated into individual functionally related regions is the first step to accurately map the human connectome. As a result, parcellation of the human brain has drawn considerable attention in the field of neuroscience. The thalamus is known as a relay in the human brain, with its nuclei sending fibers to the cortical and subcortical regions. Functional magnetic resonance imaging techniques offer a way to parcellate the thalamus in vivo based on its connectivity properties. However, the parcellations from previous studies show that both the number and the distribution of thalamic subdivisions vary with different cortical segmentation methods. In this study, we used an unsupervised clustering method that does not rely on a priori information of the cortical segmentation to parcellate the thalamus. Instead, this approach is based on the intrinsic resting-state functional connectivity profiles of the thalamus with the whole brain. A series of cluster solutions were obtained, and an optimal solution was determined. Furthermore, the validity of our parcellation was investigated through the following: (1) identifying specific resting-state connectivity patterns of thalamic parcels with different brain networks and (2) investigating the task activation and psychophysiological interactions of specific thalamic clusters during 8-Hz flashing checkerboard stimulation with simultaneous finger tapping. Together, the current study provides a reliable parcellation of the thalamus and enhances our understating of thalamic. Furthermore, the current study provides a framework for parcellation that could be potentially extended to other subcortical and cortical regions.

Introduction

The human brain is a complex system with densely interconnected regions and mapping of this “human connectome” has become a major focus of neuroscience research (Fornito et al., 2013; Sporns, 2013; Wig et al., 2014b). The first crucial step to delineate the topological properties of the human brain is node definition or delineating how the brain is parcellated into individual subunits that can be assessed for functional or structural connections with each other (Fornito et al., 2013). As a result, more and more researchers have lent their attention to parcellation of the human brain, in both cortical and subcortical areas (Bzdok et al., 2013; Craddock et al., 2012; Kim et al., 2013; Shen et al., 2013; Wig et al., 2014a).

The thalamus acts as a gateway between subcortical areas and the cerebral cortex (Liang et al., 2013; Sherman and Guillery, 2002), with nearly all areas of the neocortex receiving afferents from the thalamus (Guillery and Sherman, 2002). Moreover, the thalamus plays an important role in regulating states of sleep and wakefulness (Steriade and Llinas, 1988), and thalamocortical connectivity is reported to be fundamental for the establishment of oscillatory brain waves (Jones, 2001). Dysfunction of the thalamus is implicated in the pathophysiology of neurological and neuropsychiatric disorders (Cauda et al., 2009; Welsh et al., 2010; Woodward et al., 2012). Based on its cytoarchitecture, the thalamus can be divided into several nuclei (Morel et al., 1997) that show distinct anatomical connectivity with different cortical and subcortical regions (Jones, 2007; Sherman and Guillery, 2002). In vivo parcellation of the thalamus into different subdivisions would provide locations of specific thalamic nuclei and, thus, be of great importance for understanding thalamic functioning and revealing thalamic involvement in the pathophysiology of brain disorders.

Magnetic resonance imaging (MRI) techniques are powerful methods for parcellating brain structures in vivo based on either anatomical or functional connectivity to provide insight beyond cytoarchitecture information. For example, there are a few reports that have used a diffusion tensor MRI-based fiber tractography technique (Basser et al., 1994) to parcellate the thalamus into subdivisions according to its anatomical connectivity with seven cortical areas (Behrens et al., 2003; Johansen-Berg et al., 2005; Mastropasqua et al., 2015). Those studies (Behrens et al., 2003; Johansen-Berg et al., 2005; Mastropasqua et al., 2015) were limited to the investigation of connections between the thalamus and large cortical regions. However, as noted by Behrens and colleagues (2003), the nonhuman primate literature provides support for a finer grained topographic mapping between subregions of, for example, the mediodorsal nucleus of the thalamus and smaller regions in the prefrontal cortex (Goldman–Rakic and Porrino, 1985; Kievit and Kuypers, 1977).

Furthermore, while the anatomical structure is thought to be the foundation of brain function, studies have found a discrepancy between anatomical and functional connectivity in the human brain (Honey et al., 2009, 2010). Thus, for functional MRI (fMRI) studies, it is equally essential to parcellate the brain based on its functional connectivity. Over the last decade, resting-state fMRI has emerged as a powerful tool for mapping the spontaneous or intrinsic functional connectivity of the human brain with high sensitivity in a completely noninvasive manner that is easy to implement (Biswal et al., 1995; Kim et al., 2013; Zhang et al., 2008, 2010; Zou et al., 2009). Kim and colleagues (2013) used an ICA-based approach to parcellate the thalamus and basal ganglia into several subdivisions, and Zhang et al. used seed-based resting-state fMRI to parcellate the thalamus into five subdivisions with a Winner-Take-All method to label each thalamic voxel (Zhang et al., 2008, 2010). The number of thalamic subdivisions in ICA-based approach is directly determined by the number of independent components. A different definition of independent components will result in varied thalamic parcellation results. By definition, the approach utilized by Zhang and colleagues is based on a priori knowledge of how the cerebral cortex is segmented. As a result, parcellation of the thalamus with different templates of cortical areas, for example, using the Automated Anatomical Labeling (AAL) template (Tzourio-Mazoyer et al., 2002) or the Brodmann Area template (Brodmann, 1909), will result in a different parcellation of thalamic nuclei, both in the number of parcels and in the spatial distribution of each parcel. Previous studies have shown that using functionally inaccurate region of interests (ROIs), for example, anatomically defined ROIs, for network analysis is extremely damaging to network estimation (Smith et al., 2011).

To obviate those mentioned issues, in the present study, a recently proposed approach for parcellating brain regions based on the similarity of their functional profiles with other regions of the brain (Craddock et al., 2012; Kahnt et al., 2012; Kelly et al., 2012; van den Heuvel et al., 2008) was adopted to parcellate the thalamus into different spatially continuous subdivisions without any spatial restrictions. An optimal solution of the number of thalamic subdivisions was obtained by evaluating the stability of diverse parcellation results. The validity of this approach was investigated by assessing the following: (1) specific functional connectivity fingerprints of the thalamic parcels with known brain networks and (2) brain activation and psychophysiological interactions (PPIs) of specific thalamic nuclei during a visuomotor task (8-Hz flashing checkerboard stimulation with simultaneous finger tapping). The visuomotor task allows us to examine the validity of our parcellation because it activates specific thalamic subdivisions (i.e., lateral dorsal and ventral posterior nuclei (Purushothaman et al., 2012; Stepniewska et al., 2007; Strick, 1976)) and it can be used to assess PPIs, which characterize direct connections between thalamic nuclei and other regions during task performance. For example, Purushothaman and colleagues found that the posterior part of thalamus, which sends fibers to primary visual areas, was involved in visual signal processing.

Materials and Methods

Data acquisition

Both resting-state fMRI and visuomotor task fMRI datasets were included in this study.

Resting-state fMRI dataset

Resting-state fMRI data that were acquired at the Beijing Normal University and have been made freely available on the website of the 1000 Functional Connectomes Project (www.nitrc.org/projects/fcon_1000) were used for this study. The data from 60 subjects (30 females/30 males, 20.8±2.6 years) with complete coverage of cortical regions were selected from the available scans. All subjects provided written informed consent before the experiment, and the study was approved by the Institutional Review Board of Beijing Normal University.

MR images from this dataset were acquired using a Siemens 3T Trio MRI scanner. Resting-state functional images were obtained using a T₂*-weighted echo-planar imaging sequence with the following parameters: 33 axial slices, thickness/gap=3.0/0.6 mm, FOV=200×200 mm² with in-plane matrix=64×64, TR/TE/θ=2000 ms/30 ms/90^o, and scanning duration=7 min and 30 sec with 225 functional volumes. A high-resolution structural 3D MRI was acquired using a T₁-weighted magnetization-prepared rapid gradient echo (MPRAGE) imaging sequence (spatial resolution=1×1×1.33 mm, TR/TE/θ=2530 ms/3.39 ms/7^o, inversion time (TI)=1100 ms). During the resting-state scanning session, subjects were asked to keep their eyes closed, but not fall asleep.

Visuomotor task fMRI dataset

Visuomotor task data were acquired in eighteen healthy subjects (9 females/9 males; 22.8±1.8 years) at the Chinese Academy of Science. This study was approved by the Institutional Review Board of Chinese Academy of Science.

MR images of this dataset were also collected on a Siemens 3T Trio scanner. Task fMRI parameters of this dataset are as follows: 22 axial slices without gap, slice thickness=5 mm, FOV=220×220 mm² with in-plane matrix=64×64, and TR/TE/θ=2500 ms/30 ms/90^o. A high-resolution structural 3D MRI was acquired using a T₁-weighted MPRAGE sequence (spatial resolution=1 mm³ isotropic, TR/TE/θ=2530 ms/3.39 ms/7^o, TI=1100 ms). In the visuomotor task scans, a block design was used. It began with a 30-s off block and consisted of twelve cycles of alternated 20-s on and 20-s off blocks. During the on block, subjects were asked to view a screen with 8-Hz flashing checkerboard stimulation and tap their finger simultaneously.

Resting-state data processing

The resting-state fMRI were processed using SPM8 (www.fil.ion.ucl.ac.uk/spm/) and REST (Song et al., 2011; Yan and Zang, 2010) after discarding the first ten functional volumes of each subject to avoid transient signal changes before the longitudinal magnetization reached a steady state. Preprocessing steps are shown in a flow chart (Fig. 1). Steps in SPM8 included slice timing correction, head motion correction, and spatial normalization to the Montreal Neurological Institute (MNI) space (Collins et al., 1994). The voxel size of the normalized functional images was 3 mm isotropic. No subject had head motion larger than 2 mm in translation and/or 2° in rotation. As a result, no subject was excluded for large head motion. Linear trend removal and temporal band-pass filtering (0.01–0.08 Hz) (Biswal et al., 1995) were then performed using REST. Spatial smoothing with a Gaussian kernel of 4 mm full width at half maximum (FWHM) was applied to our imaging data. To reduce to potential contamination of non-neuronal sources of noise, the anatomical CompCor (aCompCor) method was adopted for noise estimation and removal of noise signals (Behzadi et al., 2007).

FIG. 1.

Flowchart. Procedures of resting-state functional MRI (fMRI) data analysis and thalamic parcellation.

To further reduce head motion effects, a temporal scrubbing method was applied to the preprocessed fMRI data with the threshold for framewise displacement (FD) set to 0.5 mm (Power et al., 2012). After the scrubbing procedure, most subjects (44/60) had no volumes scrubbed and the most volumes removed were 35 (in one subject). Scrubbing has been shown to improve functional connectivity maps, and in this study, this procedure will have minimal impact from slight differences in degrees of freedom from subject to subject (since most subjects did not have any volumes scrubbed and this is a single-group study). If the scrubbed data are less than 5 min, then exclude that data. In this study, no subjects are excluded in this procedure.

Resting-state functional connectivity-based parcellation

The first goal of this study was to segment the entire thalamus (both left and right parts together) into distinct subdivisions based on their resting-state functional connectivity patterns with the whole brain. The template of the thalamus was extracted from the AAL template (Tzourio-Mazoyer et al., 2002). The parcellation approach used here was modified from previous studies (Craddock et al., 2012; Kelly et al., 2012; van den Heuvel et al., 2008). A flow chart of parcellation procedure is shown in Figure 1. For each subject, functional connectivity between each thalamic voxel and every voxel of the whole brain was calculated (Fig. 1). Then, an individual similarity matrix was computed using Pearson's correlation coefficient to evaluate similarity between those spatial functional connectivity maps of each pair of thalamic voxels. A K-means clustering algorithm was applied to this similarity matrix to obtain individual thalamic parcellations, generating an individual adjacency matrix (Craddock et al., 2012). The group adjacency matrix was obtained by averaging all individual adjacency matrices. Then, the group adjacency matrix was input to the K-means clustering algorithm to obtain the group-level clustering result. For each clustering procedure (both single-subject and group levels), the number of clusters was chosen from 2 to 15.

Determination of the optimal solution

Since the cluster number, K, is a freely chosen parameter, we further assessed whether an optimal solution, that is, stable solution, existed. The variation of information (VI) index (Kahnt et al., 2012; Meilă, 2007) and Dice's coefficient (DC) (Craddock et al., 2012; Zhang et al., 2014b) were calculated here to assess the stability of parcellation procedure and determine an optimal solution. The computation of VI and DC is shown in Supplementary Materials (Supplementary Data are available online at www.liebertpub.com/brain). To obtain those indices, a split-half comparison approach was used, which divided subjects randomly into two subgroups (30 subjects in each group) for 100 instances. For each instance, the above-mentioned parcellation approach was applied to each subgroup, and VI and DC values were computed according to the parcellation results of the two subgroups. The K value, whose values of VI were statistically indistinguishable from that of the K −1 solution (Kahnt et al., 2012; Kelly et al., 2010) or with local maximum DC value, is defined as a stable solution. Combining the values of VI and DC, we can obtain the optimal solution.

Functional connectivity fingerprints of thalamic subdivisions

To define the functional connectivity fingerprints of the estimated thalamic subdivisions from the optimal solution, we calculated the Pearson correlation coefficient between mean time courses of each thalamic subdivision and those of several common resting-state brain networks (Smith et al., 2009). Nine of ten functional networks were selected, including three visual networks, the default mode network, motor network, auditory network, executive control network, and left and right frontal–parietal networks (Supplementary Fig. S1). The tenth network reported by Smith and colleagues (2009), which covered the cerebellum, was not assessed in the current study because the fMRI data from resting-state fMRI dataset did not fully cover the cerebellum. The Pearson correlation coefficients were then transformed to Fisher Z-scores and then a series of one-sample t-tests were used to assess the statistical significance of each correlation. The statistical threshold was FWE corrected at p<0.05.

Specific functional response of the thalamic subdivisions to visuomotor task

For the task fMRI dataset, preprocessing steps were as follows: (1) slice timing; (2) head motion correction; (3) spatial normalization; (4) linear trend removal; and (5) spatial smoothing with 4 mm FWHM. For the first-level analyses, a general linear model (GLM) was constructed to detect task-based activation areas in each subject. Specifically, the block design convolved with a hemodynamic response function and six head motion parameters were used as predictors in the GLM. The spatial maps of regression coefficients for the block design predictor were used as contrast maps. Group activation maps were generated from a higher level GLM analysis implementing a one-sample t-test of the contrast maps from the first-level GLM. Whole-brain results were corrected for multiple tests to control for familywise error (FWE) using Monte Carlo simulations (FWE corrected at p<0.01, with uncorrected p<0.001 and minimum cluster size=270 mm³).

To investigate whether specific thalamic clusters responded during the visuomotor task, we averaged the raw regression coefficients in each cluster from the optimal solution. Then, a one-sample t-test of average activation amplitude (computed from raw regression coefficients) for each thalamic cluster was conducted across subjects with a significance threshold of p<0.05. To investigate the functional connectivity profiles of specific clusters during task performance, each thalamic cluster was used as a seed region to perform whole-brain voxelwise PPI analysis (Friston et al., 1997). As running in SPM, the centered block design convolved with the hemodynamics response function; the first eigenvalue of BOLD signals of given thalamic subdivisions and their interactions were used as regressors in PPI analysis. Only the element for beta value of interaction term was set to one in the contrast matrix. A series of two-tailed one-sample t-tests were then conducted to test for significance of previous results with a whole-brain multiple comparison correction based on Monte Carlo simulations (corrected p<0.01 with uncorrected p<0.001, and minimum cluster size was 270 mm³). In the case of a simple first-level GLM, the averaged activation amplitudes can be used to test hypotheses related to the simple question of “what regions in the brain have activity that is being modulated by the task?” In the case of PPI, the hypotheses being tested are asking the question “what regions of the brain show increased connectivity with the thalamic nucleus during the visuomotor stimulation itself (e.g., relative to the control condition)?” These are very different hypotheses that inform different aspects of thalamic nuclei behavior during visuomotor stimulation. Thus, a PPI analysis provides strong evidence for testing our parcellation scheme, as direct connections of thalamic nuclei with visual and motor regions should be separately disentangled (whereas the GLM activation results will just show nuclei involved in task response, but not their direct connections during the task).

Results

Thalamus parcellation results

Parcellation results with different cluster solutions obtained from the resting-state data are shown in Figure 2. The K=2 cluster solution revealed two spatially symmetric subdivisions of the thalamus. The larger cluster (Cluster 1) stretched from the ventral anterior part to the medial and almost covered the complete dorsal part of the thalamus. The other cluster (Cluster 2) covered only the ventral posterior and ventral lateral part of the thalamus. Interestingly, the two subdivisions appeared in both the left and right hemispheres and showed excellent symmetry. As we parcellated the entire thalamus together, the symmetric results indicated that the parcellation was dependent on the functional connectivity profiles and not constrained by the distance relationship (Mars et al., 2012; Tomassini et al., 2007). Compared with the K=2 cluster solution, Cluster 3 from the K=3 solution mostly stemmed from Cluster 2 (the ventral posterior cluster) of the K=2 solution. Cluster 4 from the K=4 solution was spilt from the medial cluster from K=3 and it covered the posterior part of the medial thalamus. Cluster solutions with K=5 to K=8 resulted in further splitting of these clusters into subdivisions to show a clear hierarchical structure of the parcellations. As it is shown in Figure 2, Cluster 9 from the K=9 solution was lateralized to the left part of the thalamus. Similarly, when K>9, one or more clusters were lateralized to a single hemisphere, thus those results are not shown here.

FIG. 2.

The resting-state functional connectivity-based thalamic parcellation. Resting-state functional connectivity-based parcellation of the thalamus with different cluster solutions (K=2, 3…9). Seven slices (Z values from −1 mm to 17 mm, slice thickness: 3 mm) covering the entire thalamus are shown. Different colors represent different thalamic subdivisions from the parcellation. L and R denote the left and right hemispheres, respectively. Color images available online at www.liebertpub.com/brain

Determination of the optimal cluster solution

To determine the optimal solution, a split-half procedure with 100 randomly generated split-half groups was used and VI and DC values were computed (Fig. 3) to assess the stability of the parcellation and the similarity between the adjacency matrices calculated for each group in the split-half procedure. VI increased as a function of K, indicating that the stability decreased as K increased (Fig. 3A). After performing a series of paired t-test statistics, VI was found to be significantly different with K (p<0.001), except for K=7. VI indices for the K=6 and K=7 parcellation were statistically undistinguishable (p=0.988), moreover, the parcellation with K=7 corresponded to a local maximum of the DC values (Fig. 3B); thus, K=7 was selected as the optimal solution.

FIG. 3.

VI and DC values for each cluster solution (K). (A) Values of VI as a function of K. Values of VI for each K cluster solution are significantly different from the K+1 solution (p<0.001), except for the difference between K=6 and K=7 (p=0.988), (B) Values of DC as a function of K. K=7 is a local maximum and, thus, considered to be a relatively stable solution. VI, variation of information; DC, Dice's coefficient.

Functional connectivity fingerprints between the thalamic parcellation and resting-state brain networks

The functional connectivity fingerprints of the thalamic subdivisions from the optimal thalamic parcellation (K=7) are shown in Figure 4, and the connectivity strength computed as values of the group averaged correlation coefficients between mean time courses of each thalamic subdivision and each resting-state brain network is listed in Table 1. It can be seen that different subdivisions of the thalamus showed distinct functional connectivity fingerprints. Table 1 shows that all statistically significant connections are positive. Cluster 1 was significantly correlated with default mode network, motor network, auditory network, executive control network, and bilateral frontal–parietal networks (Fig. 4). Cluster 2 demonstrated significant functional connectivity with all networks, except for bilateral frontal–parietal networks. Among these connections, the one with the motor network has the greatest connectivity strength (r=0.341). Cluster 3 showed significant but weak correlations with motor, auditory, and executive control networks. Cluster 4 and Cluster 6, which are almost the same cluster in the K=5 solution, showed a similar connectivity fingerprint, with very significant connections with the default mode network and executive control network. Cluster 5 and Cluster 7, although with distinct connectivity patterns, showed significant connectivity with almost all networks except for occipital pole visual and lateral visual networks.

FIG. 4.

Functional fingerprints of thalamic subdivisions from the optimal solution (K=7). Fingerprints of Clusters 1–7 from K=7 solution. Each value of connectivity strength in fingerprint is the group averaged correlation coefficient between mean time course of each thalamic subdivision and each functional network. The MNI axial slice for the overlay corresponds to Z=11 mm. VN1, medial visual network; VN2, occipital pole visual network; VN3, lateral visual network; DMN, default mode network; MN, motor network; AN, auditory network; ECN, executive control network; rFPN, right frontal–parietal network; lFPN, left frontal–parietal network. Color images available online at www.liebertpub.com/brain

Table 1.

Mean Functional Connectivity Strength Between Thalamic Subdivisions and Brain Networks

	Cluster number
Brain networks	1	2	3	4	5	6	7
VN1	0.034	0.196	0.011	0.117	0.134	0.209	0.150
VN2	−0.006	0.139	0.000	0.125	0.087	0.206	0.101
VN3	−0.036	0.175	0.030	−0.026	0.093	0.018	0.030
DMN	0.172	0.176	0.013	0.335	0.163	0.392	0.310
MN	0.127	0.341	0.157	0.117	0.294	0.063	0.151
AN	0.164	0.225	0.134	0.009	0.313	0.153	0.249
ECN	0.355	0.136	0.141	0.299	0.289	0.543	0.459
rFPN	0.278	0.073	0.067	0.171	0.160	0.228	0.283
lFPN	0.261	0.069	0.048	0.147	0.120	0.192	0.247

Thalamic subdivisions are those from the optimal solution (K=7) using the resting-state fMRI data. Functional connectivity values in bold indicate statistically significant connectivity across subjects (corrected p<0.05).

VN1, medial visual network; VN2, occipital pole visual network; VN3, lateral visual network; DMN, default mode network; MN, motor network; AN, auditory network; ECN, executive control network; rFPN, right frontal–parietal network; lFPN, left frontal–parietal network.

Activations and PPI of the thalamic subdivisions during the visuomotor task

Visuomotor task activation of thalamic regions is shown in Figure 5A. The activated regions are located mainly in the bilateral ventral posterior part and left ventral lateral part of the thalamus, which overlap primarily with Cluster 2, Cluster 3, and Cluster 5 for the K=7 solution. The percentage overlap between visuomotor activated regions and the thalamic subdivisions is listed in Table 2. The activation amplitude of each thalamic cluster is shown in Figure 5B. Clusters 2, 3, and 5 showed statistically significant responses to visuomotor stimulation.

FIG. 5.

Response of specific thalamic clusters from the optimal solution (K=7) to visuomotor stimulation. (A) Task activation results are shown in the upper panel (corrected p<0.01) with the corresponding thalamic parcellation at the same Z locations. (B) Bar and error bar plots of activation amplitude for each thalamic subdivision, averaged across subjects. Clusters 2, 3, and 5 show significant activation during the task. (C) Regions showing significant PPI effects with thalamic clusters during the visuomotor task period relative to the control period. Color images available online at www.liebertpub.com/brain

Table 2.

Percentage Overlap Between Thalamic Subdivisions and Visuomotor Task Activations

	Cluster number
	1	2	3	4	5	6	7
Task activations	0%	50.8%	15.9%	3.4%	25.9%	0%	4.0%

Thalamic subdivisions are those from K=7 parcellation using the resting-state fMRI data (with global signal regression).

PPIs during the visuomotor task condition relative to baseline were investigated between every subdivision of the optimal solution (K=7) and the whole brain. Only clusters (Cluster 2, Cluster 3, and Cluster 5) activated by the visuomotor task showed significant negative PPIs with the visual cortex or/and motor cortex (Fig. 5C). Clusters 2 and 5 showed negative PPIs with visual and motor cortices, while Cluster 3 only showed a negative PPI with the visual cortex. There were no significant positive PPIs of thalamic subdivisions with other regions of the brain.

Discussion

In the present study, we parcellated the thalamus into different subdivisions using an unsupervised approach based on their resting-state functional connectivity profiles with the whole brain. The validity of our thalamic parcellation approach was tested by assessing the following: (1) the functional connectivity fingerprints of our thalamic parcels with known brain networks and (2) the thalamic responses and PPIs with visual and motor regions during visuomotor stimulation.

Connectivity-based parcellation of the thalamus

Behrens and colleagues proposed a novel method to parcellate the thalamus based on its anatomical connectivity with a priori segmented cortical areas (Behrens et al., 2003; Johansen-Berg et al., 2005; Mastropasqua et al., 2015). Depending on the number of cortical regions used for the parcellation, they obtained different parcellation results, that is, parcellation of the thalamus into five or seven subdivisions. In the human brain, the anatomical structure is assumed to be the foundation of its function. With equal importance, the underlying anatomical structure is continually modulated by function in relation to experience (Damoiseaux and Greicius, 2009; Zhang et al., 2010). It is commonly assumed that functional connectivity between brain regions reflects the anatomical connectivity between them (Damoiseaux and Greicius, 2009; Honey et al., 2009, 2010). However, the relationship between functional and structural connectivity appears not be so straightforward. Studies have shown that the strength of functional and structural connectivity is positively correlated when structural connectivity is present; functional connectivity is also observed between brain regions without structural connectivity (Honey et al., 2009). Thus, a thalamic parcellation based on functional connectivity properties is needed for functional MRI studies.

To that end, a few studies have defined thalamic subdivisions based on functional connectivity with a template of brain cortical regions (Zhang et al., 2008, 2010). Kim and colleagues (2013) proposed a multilevel independent component analysis (ICA) technique to parcellate the thalamus and basal ganglia. However, ICA has some intrinsic methodological limitations (Hyvarinen, 2011), specifically, ICA-based schemes will depend on the choice of dimensionality of the ICA. In the current study, we adopted a technique based on the intrinsic functional connectivity profiles of thalamus. Compared to previous seed-based connectivity parcellation methods, it does not rely on a priori definition of cortical areas. Compared to ICA-based techniques, our approach is more hypothesis driven in the sense of the determination of number of nuclei. Moreover, the ICA will separate out the connectivity of each parcel into different components according to the networks they belong to. This makes interpretation a little trickier as one has to worry about mixing in many different components. However, the seed-based connectivity approach, as what we do, can assess the whole connectivity profile of each parcel. Thus, the results of the proposed methods are perhaps easier to interpret.

Resting-state functional connectivity with specific brain networks

Distinct and anatomically relevant functional connectivity fingerprints of each subdivision with brain networks provided evidence for the validity of our parcellation approach. For instance, the Cluster 6 mainly overlapped with the mediodorsal part of the thalamus. Significant functional connectivity was observed between this subdivision and the executive control and default mode networks (Fig. 4). Previous studies have shown that the mediodorsal nucleus anatomically projects to the prefrontal cortex and posterior cingulate cortex (Jones, 2007), which are key regions in the executive control network and the default mode network, respectively. Thus, our findings suggest a close correspondence between the functional and anatomical connectivity of Cluster 6. Notably, Clusters 2 and 3 are primarily located in the ventral lateral and ventral posterior parts of the thalamus, which sends fibers to motor and sensory cortices. Figure 4 shows that these subdivisions correlated with the motor network. Moreover, Cluster 4 mainly overlapped with the medial part of the pulvinar nucleus. This part of thalamus is anatomically connected with the prefrontal cortex, cingulate cortex, and visual and auditory cortices (Jones, 2007). Thus, the significant functional connectivity observed between Cluster 4 and these cortical areas (Fig. 4 and Table 1) is consistent with the known anatomical connections of the pulvinar.

Response of specific thalamic clusters to visuomotor stimulation

Visuomotor task activation

Activity of the thalamic clusters from the K=7 parcellation showed responses in certain clusters to visuomotor stimulation that are consistent with previous reports of visuomotor task activation (Calhoun et al., 2001; Miao et al., 2014; Witt et al., 2008), providing validity of our thalamic parcellation results. Specifically, Clusters 2, 3, and 5 were activated during the visuomotor task. Cluster 2 primarily comprised the posterior lateral part of the thalamus, a region that has been shown to project to both primary and higher order visual areas and to contribute to visual perception and attention (Fischer and Whitney, 2012; Purushothaman et al., 2012). The ventral lateral and ventral posterior parts of the thalamus, which were in Cluster 3 and Cluster 5, relay somatosensory and motor information to and among the brain (Morel et al., 1997; Stepniewska et al., 2007; Strick, 1976) and these regions were activated during the finger-tapping task when the flashing checkerboard stimulation was on.

PPI during the visuomotor task

Interestingly, only those thalamic clusters, which were activated by the visuomotor task, showed a significant PPI with visual and motor cortices. However, the PPI analysis revealed that in fact Clusters 2, 3, and 5 showed strong negative connectivity with visual and/or motor cortices. Although this may seem counterintuitive, theoretically, it can be explained. The statistical model used to detect brain activation (which is a GLM with a predictor that corresponds to the task plus any nuisance predictors) detects regions with increased or decreased BOLD signal responses to a given task (e.g., a main effect of the task), while a PPI analysis, which is also a GLM, but includes a ROI time course and an interaction regressor, investigates the interaction between the seed ROI and other regions of the brain during task performance. In other words, the PPI helps distinguish whether two regions activated by a task become more strongly functionally connected during task performance (e.g., have a stronger positive or negative correlation) or whether the regions act independently with both becoming more active during the task. Several previous studies have reported negative PPIs of thalamus with other brain regions during task performance. For example, Chang and colleagues (2011) reported a negative PPI between the ventral lateral thalamus and primary auditory cortex during a speech production task, and Zhang and colleagues (2014a) reported a negative PPI between the thalamus and medial prefrontal cortex during a stop signal task. Second, a recent animal electrophysiology study found that the increases in the thalamic action potential firing rate during sensory stimulation correlated with decreases in the amplitude of low-frequency local field potential (LFP) fluctuations of neurons from the barrel cortex (Poulet et al., 2012). BOLD fMRI signals are closely related to LFP signals (Logothetis et al., 2001), therefore an increase in thalamic action potential firing leading to a decrease in the amplitude of the LFP could be one possible mechanism for the strong negative PPIs we observed between these two regions during visuomotor stimulation, although clearly this is speculative. Thus, the presence of a strong PPI effect between those activated clusters and visual and motor cortices provides further evidence to the validity of our parcellation scheme. For example, it is consistent that thalamic nuclei would interact with visual and motor cortices during a visuomotor task. Our findings reveal that this occurs through negative PPIs with the thalamic nuclei, which sheds insight into mechanisms underlying thalamocortical interactions. Future work to understand the exact processes giving rise to these negative PPIs is warranted.

Limitations

There are potential limitations of our parcellation procedure that should be considered. First, the parcellation results were performed on group-averaged data. More work is needed to obtain stable individual parcellation results. Second, considering the convergence and divergence between functional connectivity and structural connectivity (Buckner et al., 2008; Honey et al., 2010), a multimodal diffusion-based anatomical scan (e.g., diffusion spectrum imaging) to assess structural connectivity combined with functional connectivity measurements would shed great light on the correspondence between structural and functional connectivity and perhaps result in a more robust thalamic parcellation. Third, improvements in the signal to noise ratio and spatial resolution would result in a parcellation of the thalamus with higher precision. Fourth, considering that thalamus is surrounded by the ventricular system, spatial smoothing might be a potential source of contamination of thalamic signals. In this study, to minimize the effect of smoothing, a 4-mm FWHM smooth kernel was used, with a 3 mm voxel size of normalized functional images. In addition, the seed-based functional connectivity analysis method has intrinsic shortcomings. For example, it is less sensitive in detecting interindividual variability than ICA-based analyses (Smith et al., 2014) and one must contend with the global signal regression problem. In this study, we used the aCompCor approach for removing physiological noise signals that contribute to global signals and, thus, obviated the issue of GSR.

Conclusion

In the current study, we used an unsupervised clustering approach to parcellate the thalamus based on resting-state functional connectivity profiles. Our optimal thalamic parcellation scheme resulting in five symmetric thalamic clusters was validated by assessing the distinct fingerprints of functional connectivity profiles of the resulting thalamic clusters with brain networks and by assessing the thalamic responses and PPIs during a visuomotor task.

Footnotes

Acknowledgments

This work was supported by China's National Strategic Basic Research Program (973) (2012CB720700 and 2015CB856400) and the Natural Science Foundation of China (81227003, 81430037, 31421003, 31200761, and 81201142).

Author Disclosure Statement

No competing financial interests exist.

References

Basser

, Mattiello

, LeBihan

. 1994. MR diffusion tensor spectroscopy and imaging. Biophys J, 66:259–267.

Behrens

, Johansen-Berg

, Woolrich

, Smith

, Wheeler-Kingshott

, Boulby

, et al. 2003. Non-invasive mapping of connections between human thalamus and cortex using diffusion imaging. Nat Neurosci, 6:750–757.

Behzadi

, Restom

, Liau

, Liu

. 2007. A component based noise correction method (CompCor) for BOLD and perfusion based fMRI. NeuroImage, 37:90–101.

Biswal

, Yetkin

, Haughton

, Hyde

. 1995. Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn Reson Med, 34:537–541.

Brodmann

. 1909. Vergleichende Lokalisationslehre der Grosshirnrinde in ihren Prinzipien dargestellt auf Grund des Zellenbaues. Barth: Leipzig, Germany.

Buckner

, Andrews-Hanna

, Schacter

. 2008. The brain's default network: anatomy, function, and relevance to disease. Ann N Y Acad Sci, 1124:1–38.

Bzdok

, Laird

, Zilles

, Fox

, Eickhoff

. 2013. An investigation of the structural, connectional, and functional subspecialization in the human amygdala. Hum Brain Mapp, 34:3247–3266.

Calhoun

, Adali

, McGinty

, Pekar

, Watson

, Pearlson

. 2001. fMRI activation in a visual-perception task: network of areas detected using the general linear model and independent components analysis. NeuroImage, 14:1080–1088.

Cauda

, Sacco

, D'Agata

, Duca

, Cocito

, Geminiani

, et al. 2009. Low-frequency BOLD fluctuations demonstrate altered thalamocortical connectivity in diabetic neuropathic pain. BMC Neurosci, 10:138.

10.

Chang

S-E

, Horwitz

, Ostuni

, Reynolds

, Ludlow

. 2011. Evidence of left inferior Frontal–Premotor structural and functional connectivity deficits in adults who stutter. Cereb Cortex, 21:2507–2518.

11.

Collins

, Neelin

, Peters

, Evans

. 1994. Automatic 3D intersubject registration of MR volumetric data in standardized Talairach space. J Comp Assist Tomogr, 18:192–205.

12.

Craddock

, James

, Holtzheimer

, Hu

, Mayberg

. 2012. A whole brain fMRI atlas generated via spatially constrained spectral clustering. Hum Brain Mapp, 33:1914–1928.

13.

Damoiseaux

, Greicius

. 2009. Greater than the sum of its parts: a review of studies combining structural connectivity and resting-state functional connectivity. Brain Struct Funct, 213:525–533.

14.

Fischer

, Whitney

. 2012. Attention gates visual coding in the human pulvinar. Nat Commun, 3:1051.

15.

Fornito

, Zalesky

, Breakspear

. 2013. Graph analysis of the human connectome: promise, progress, and pitfalls. NeuroImage, 80:426–444.

16.

Friston

, Buechel

, Fink

, Morris

, Rolls

, Dolan

. 1997. Psychophysiological and modulatory interactions in neuroimaging. NeuroImage, 6:218–229.

17.

Goldman–Rakic

, Porrino

. 1985. The primate mediodorsal (MD) nucleus and its projection to the frontal lobe. J Comp Neurol, 242:535–560.

18.

Guillery

, Sherman

. 2002. Thalamic relay functions and their role in corticocortical communication: generalizations from the visual system. Neuron, 33:163–175.

19.

Honey

, Sporns

, Cammoun

, Gigandet

, Thiran

, Meuli

, et al. 2009. Predicting human resting-state functional connectivity from structural connectivity. Proc Natl Acad Sci U S A, 106:2035–2040.

20.

Honey

, Thivierge

, Sporns

. 2010. Can structure predict function in the human brain?. NeuroImage, 52:766–776.

21.

Hyvarinen

. 2011. Testing the ICA mixing matrix based on inter-subject or inter-session consistency. NeuroImage, 58:122–136.

22.

Johansen-Berg

, Behrens

, Sillery

, Ciccarelli

, Thompson

, Smith

, et al. 2005. Functional-anatomical validation and individual variation of diffusion tractography-based segmentation of the human thalamus. Cereb Cortex, 15:31–39.

23.

Jones

. 2007. The Thalamus. Cambridge: Cambridge Univ. Press.

24.

Jones

. 2001. The thalamic matrix and thalamocortical synchrony. Trends Neurosci, 24:595–601.

25.

Kahnt

, Chang

, Park

, Heinzle

, Haynes

. 2012. Connectivity-based parcellation of the human orbitofrontal cortex. J Neurosci, 32:6240–6250.

26.

Kelly

, Toro

, Di Martino

, Cox

, Bellec

, Castellanos

, et al. 2012. A convergent functional architecture of the insula emerges across imaging modalities. NeuroImage, 61:1129–1142.

27.

Kelly

, Uddin

, Shehzad

, Margulies

, Castellanos

, Milham

, et al. 2010. Broca's region: linking human brain functional connectivity data and non-human primate tracing anatomy studies. Eur J Neurosci, 32:383–398.

28.

Kievit

, Kuypers

. 1977. Organization of the thalamo-cortical connexions to the frontal lobe in the rhesus monkey. Exp Brain Res, 29:299–322.

29.

Kim

, Park

. 2013. Functional connectivity–based identification of subdivisions of the basal ganglia and thalamus using multilevel independent component analysis of resting state fMRI. Hum Brain Mapp, 34:1371–1385.

30.

Liang

, Li

, King

, Zhang

. 2013. Mapping thalamocortical networks in rat brain using resting-state functional connectivity. NeuroImage, 83:237–244.

31.

Logothetis

, Pauls

, Augath

, Trinath

, Oeltermann

. 2001. Neurophysiological investigation of the basis of the fMRI signal. Nature, 412:150–157.

32.

Mars

, Sallet

, Schuffelgen

, Jbabdi

, Toni

, Rushworth

. 2012. Connectivity-based subdivisions of the human right “temporoparietal junction area”: evidence for different areas participating in different cortical networks. Cereb Cortex, 22:1894–1903.

33.

Mastropasqua

, Bozzali

, Spanò

, Koch

, Cercignani

. 2015. Functional anatomy of the thalamus as a model of integrated structural and functional connectivity of the human brain in vivo. Brain Topogr, 28:548–558.

34.

Meilă

. 2007. Comparing clusterings—an information based distance. J Multivariate Anal, 98:873–895.

35.

Miao

, Gu

, Yan

, Lu

, Wang

, Zhou

, et al. 2014. Detecting resting-state brain activity by spontaneous cerebral blood volume fluctuations using whole brain vascular space occupancy imaging. NeuroImage, 84:575–584.

36.

Morel

, Magnin

, Jeanmonod

. 1997. Multiarchitectonic and stereotactic atlas of the human thalamus. J Comp Neurol, 387:588–630.

37.

Poulet

, Fernandez

, Crochet

, Petersen

. 2012. Thalamic control of cortical states. Nat Neurosci, 15:370–372.

38.

Power

, Barnes

, Snyder

, Schlaggar

, Petersen

. 2012. Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion. NeuroImage, 59:2142–2154.

39.

Purushothaman

, Marion

, Li

, Casagrande

. 2012. Gating and control of primary visual cortex by pulvinar. Nat Neurosci, 15:905–912.

40.

Shen

, Tokoglu

, Papademetris

, Constable

. 2013. Groupwise whole-brain parcellation from resting-state fMRI data for network node identification. NeuroImage, 82:403–415.

41.

Sherman

, Guillery

. 2002. The role of the thalamus in the flow of information to the cortex. Philos Trans R Soc Lond B Biol Sci, 357:1695–1708.

42.

Smith

, Utevsky

, Bland

, Clement

, Clithero

, Harsch

, et al. 2014. Characterizing individual differences in functional connectivity using dual-regression and seed-based approaches. NeuroImage, 95:1–12.

43.

Smith

, Fox

, Miller

, Glahn

, Fox

, Mackay

, et al. 2009. Correspondence of the brain's functional architecture during activation and rest. Proc Natl Acad Sci U S A, 106:13040–13045.

44.

Smith

, Miller

, Salimi-Khorshidi

, Webster

, Beckmann

, Nichols

, et al. 2011. Network modelling methods for fMRI. NeuroImage, 54:875–891.

45.

Song

, Dong

, Long

, Li

, Zuo

, Zhu

, et al. 2011. REST: a toolkit for resting-state functional magnetic resonance imaging data processing. PloS One, 6:e25031.

46.

Sporns

. 2013. The human connectome: origins and challenges. NeuroImage, 80:53–61.

47.

Stepniewska

, Preuss

, Kaas

. 2007. Thalamic connections of the dorsal and ventral premotor areas in New World owl monkeys. Neuroscience, 147:727–745.

48.

Steriade

, Llinas

. 1988. The functional states of the thalamus and the associated neuronal interplay. Physiol Rev, 68:649–742.

49.

Strick

. 1976. Anatomical analysis of ventrolateral thalamic input to primate motor cortex. J Neurophysiol, 39:1020–1031.

50.

Tomassini

, Jbabdi

, Klein

, Behrens

, Pozzilli

, Matthews

, et al. 2007. Diffusion-weighted imaging tractography-based parcellation of the human lateral premotor cortex identifies dorsal and ventral subregions with anatomical and functional specializations. J Neurosci, 27:10259–10269.

51.

Tzourio-Mazoyer

, Landeau

, Papathanassiou

, Crivello

, Etard

, Delcroix

, et al. 2002. Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. NeuroImage, 15:273–289.

52.

van den Heuvel

, Mandl

, Pol

. 2008. Normalized cut group clustering of resting-state fMRI data. PloS One, 3:e2001.

53.

Welsh

, Chen

, Taylor

. 2010. Low-frequency BOLD fluctuations demonstrate altered thalamocortical connectivity in schizophrenia. Schizophr Bull, 36:713–722.

54.

Wig

, Laumann

, Cohen

, Power

, Nelson

, Glasser

, et al. 2014a. Parcellating an individual subject's cortical and subcortical brain structures using snowball sampling of resting-state correlations. Cereb Cortex, 24:2036–2054.

55.

Wig

, Laumann

, Petersen

. 2014b. An approach for parcellating human cortical areas using resting-state correlations. NeuroImage, 93:276–291.

56.

Witt

, Laird

, Meyerand

. 2008. Functional neuroimaging correlates of finger-tapping task variations: an ALE meta-analysis. NeuroImage, 42:343–356.

57.

Woodward

, Karbasforoushan

, Heckers

. 2012. Thalamocortical dysconnectivity in schizophrenia. Am J Psychiatry, 169:1092–1099.

58.

Yan

C-G

, Zang

Y-F

. 2010. DPARSF: A MATLAB Toolbox for “Pipeline” Data Analysis of Resting-State fMRI. Front Syst Neurosci, 4:13.

59.

Zhang

, Snyder

, Fox

, Sansbury

, Shimony

, Raichle

. 2008. Intrinsic functional relations between human cerebral cortex and thalamus. J Neurophysiol, 100:1740–1748.

60.

Zhang

, Snyder

, Shimony

, Fox

, Raichle

. 2010. Noninvasive functional and structural connectivity mapping of the human thalamocortical system. Cereb Cortex, 20:1187–1194.

61.

Zhang

, Hu

, Bednarski

, Erdman

, Li

C-sR

. 2014a. Error-related functional connectivity of the thalamus in cocaine dependence. NeuroImage Clin, 4:585–592.

62.

Zhang

, Fan

, Zhang

, Wang

, Zhu

, Zhang

, et al. 2014b. Connectivity-based parcellation of the human posteromedial cortex. Cereb Cortex, 24:719–727.

63.

Zou

, Long

, Zuo

, Yan

, Zhu

, Yang

, et al. 2009. Functional connectivity between the thalamus and visual cortex under eyes closed and eyes open conditions: a resting-state fMRI study. Hum Brain Mapp, 30:3066–3078.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.19 MB