Abstract
Decades of research have greatly improved our understanding of intrinsic human brain organization in terms of functional networks and the transmodal hubs within the cortex at which they converge. However, substrates of multinetwork integration in the human subcortex are relatively uncharted. Here, we leveraged recent advances in subcortical atlasing and ultra-high field (7 T) imaging optimized for the subcortex to investigate the functional architecture of 14 individual structures in healthy adult males and females with a fully data-driven approach. We revealed that spontaneous neural activity in subcortical regions can be decomposed into multiple independent subsignals that correlate with, or “echo,” the activity in functional networks across the cortex. Distinct subregions of the thalamus, striatum, claustrum, and hippocampus showed a varied pattern of echoes from attention, control, visual, somatomotor, and default mode networks, demonstrating evidence for a heterogeneous organization supportive of functional integration. Multiple network activity furthermore converged within the globus pallidus externa, substantia nigra, and ventral tegmental area but was specific to one subregion, while the amygdala and pedunculopontine nucleus preferentially affiliated with a single network, showing a more homogeneous topography. Subregional connectivity of the globus pallidus interna, subthalamic nucleus, red nucleus, periaqueductal gray, and locus coeruleus did not resemble patterns of cortical network activity. Together, these finding describe potential mechanisms through which the subcortex participates in integrated and segregated information processing and shapes the spontaneous cognitive dynamics during rest.
SIGNIFICANCE STATEMENT Despite the impact of subcortical dysfunction on brain health and cognition, large-scale functional mapping of subcortical structures severely lags behind that of the cortex. Recent developments in subcortical atlasing and imaging at ultra-high field provide new avenues for studying the intricate functional architecture of the human subcortex. With a fully data-driven analysis, we reveal subregional connectivity profiles of a large set of noncortical structures, including those rarely studied in fMRI research. The results have implications for understanding how the functional organization of the subcortex facilitates integrative processing through cross-network information convergence, paving the way for future work aimed at improving our knowledge of subcortical contributions to intrinsic brain dynamics and spontaneous cognition.
Introduction
A large body of research in the past decades has focused on descriptions of the macroscopic organization of the human brain in terms of intrinsic functional connectivity (FC) and its role in orchestrating cognition and behavior (Damoiseaux et al., 2006; W. H. Lee et al., 2019; Liégeois et al., 2019). The integration of distributed, functionally specialized brain networks is thought to be essential, especially for higher-level cognition and consciousness (Senden et al., 2014; Bell and Shine, 2016). With a variety of methods, specific sites for network convergence have been identified in the posterior cingulate cortex (PCC), anterior cingulate cortex, and the posterior parietal cortices (Tomasi and Volkow, 2011; Bell and Shine, 2015; Lyu et al., 2021), revealing an ensemble of transmodal regions in the cortex that enable efficient global communication (Van der Heuvel and Sporns, 2011; Grayson et al., 2014). With a novel multivariate approach, it was revealed that subtle signals from functionally specialized subdivisions within these regions have connectivity profiles that mirror, or “echo,” the activity of different networks, potentially indicating a mechanism through which they facilitate cross-network information integration (Leech et al., 2012; Braga et al., 2013; Braga and Leech, 2015).
Although this work has provided important insights, the dominating corticocentric view overlooks potential contributions from the highly diverse and interconnected structures in the subcortex (Bell and Shine, 2016; Forstmann et al., 2017; Tian et al., 2020). This knowledge gap is likely related to the challenges associated with visualizing the subcortex using conventional MRI because of the varied magnetic tissue properties and generally weaker signal-to-noise ratio (SNR) compared with the cortex (De Hollander et al., 2017; Keuken et al., 2018). Nonetheless, many subcortical structures are part of extensive cortico-subcortical circuitry and demonstrate widespread FC to networks, including the default mode network (Haber, 2003; Bär et al., 2016; T. W. Lee and Xue, 2018; Ji et al., 2019; Li et al., 2021). Compared with the smaller subcortical nuclei in the deep brain, larger structures, such as the thalamus (Tha) and striatum (Str), have received a relatively high amount of attention, establishing their hub-like properties and roles in integrative processing (Choi et al., 2012; Jarbo and Verstynen, 2015; Hwang et al., 2017; Greene et al., 2020; Seitzman et al., 2020; Cheng and Liu, 2021). However, most of the subcortex remains underrepresented in human fMRI studies, and the majority of available evidence is based on lower field strength (3 Tesla), often combined with extensive spatial smoothing, both of which limit the spatial resolution needed to resolve smaller nuclei and increase the risk for signal blurring (De Hollander et al., 2015; Forstmann et al., 2017).
Because of these shortcomings, the functional architecture of the subcortex and its role in integrative processing remain poorly understood. Given that subcortical dysfunction is heavily implicated in a wide range of neuropsychiatric diseases, advancing this knowledge may be vital for our understanding of healthy cognitive functioning as well as improving disease models. Charting the topography of network echoes within the subcortex provides a compelling approach to accomplish new insights into the subcortical contributions to whole-brain communication and higher-level cognition. Following previous work (Leech et al., 2012; Braga et al., 2013), we define an echo as a unique subregional connectivity profile that traces the activity pattern of a functional network. By leveraging recent advances in automated parcellation algorithms and sensitive fMRI protocols for the subcortex at ultra-high field (Bazin et al., 2020; Miletic et al., 2020), we aim to extend the previously established multivariate echo analysis to a large set of subcortical structures, including those rarely studied with human fMRI: the Tha, Str, globus pallidus externa (GPe), globus pallidus interna (GPi), subthalamic nucleus (STN), claustrum (Cl), hippocampus (HPC), amygdala (Amg), substantia nigra (SN), red nucleus (RN), ventral tegmental area (VTA), locus coeruleus (LC), periaqueductal gray (PAG), and pedunculopontine nucleus (PPN). Similar to findings for the cortex, we expect that subcortical structures organized to facilitate multinetwork integration demonstrate a heterogeneous subregional topography of intrinsic echoes from separate functional networks, which are likely hidden with previous univariate connectivity analyses.
Materials and Methods
Participants
The study was approved by the Ethics Review Board of the University of Amsterdam and the Regional Committees for Medical and Health Research Ethics in Norway. Forty healthy adults between 19 and 39 years old (21 female, mean age = 26.5 years, SD = 5.5 years) were recruited from the general population in Norway and screened for MRI compatibility. Exclusion criteria were self-reported (history of) neurologic or psychiatric disease, impaired vision, or any contraindications for MRI, such as metal implants. Written informed consent was obtained from all participants before data collection. All materials, code, and unthresholded group-level statistical maps from multivariate as well as (supplementary) univariate connectivity analyses are publicly available in an Open Science Framework repository at https://osf.io/wt3uc.
fMRI acquisition and preprocessing
Neuroimaging data were collected with a Siemens MAGNETOM Terra 7 Tesla (7 T) system with a 32-channel phased-array head coil. Structural images were obtained with a MP2RAGE sequence (Marques et al., 2010) in 224 sagittal slices at 0.75 mm isotropic voxel resolution (TR = 4300 ms; TI1,2 = 840, 2370 ms; flip angles1,2 = 5, 6°; TE = 1.99 ms; FOV = 240 × 240 × 168 mm). Functional images were acquired using a gradient EPI sequence with a voxel resolution of 1.5 mm isotropic (82 transverse slices per volume; TR = 1380 ms; TE = 14 ms; flip angle = 60°; in-plane acceleration factor (GRAPPA) = 3; multiband acceleration factor = 2; partial Fourier = 6/8). An additional EPI sequence with opposite phase-encoding direction was performed for susceptibility distortion correction purposes. Heart rate and respiratory data were acquired with a fingerclip and waistband, respectively, to correct for physiological noise, which is especially prominent in the subcortex.
MR images were preprocessed with fMRIPrep (version 20.2.6) (Esteban et al., 2019) in the Nipype framework (Gorgolewski et al., 2011). The structural (T1-weighted) scan was corrected for intensity nonuniformity with N4BiasFieldCorrection (ANTs version 2.3.3) (Tustison et al., 2010) and skull-stripped with antsBrainExtraction using the OASIS30ANTs target template. Brain tissue segmentation of CSF, white matter, and gray matter was performed with FAST (FSL version 5.0.9) (Zhang et al., 2001). For each of the two resting-state runs, a reference volume and its skull-stripped version were generated. A fieldmap based on the EPI references with opposing phase-encoding directions was calculated with 3dQwarp (AFNI) (Cox, 1996), and susceptibility distortion correction was applied to the EPI reference before coregistration to the T1-weighted reference using the boundary-based registration cost-function in bbregister with 6 degrees of freedom (FreeSurfer) (Greve and Fischl, 2009). Head-motion parameters (rotation and translation) were estimated with MCFLIRT (FSL version 5.0.9) (Jenkinson et al., 2002), and slice-time correction to half of the acquisition range (0.674 s) was performed with AFNI’s 3dTshift. Following fMRIPrep, data were spatially smoothed with a FWHM Gaussian kernel of 1.5 mm using SUSAN (S. M. Smith and Brady, 1997) and denoised with a first-level GLM in FEAT (Woolrich et al., 2001) that included fMRIPrep-derived confound regressors, including the following: mean signal in CSF and white matter, framewise displacement, six rotation and translation parameters, and discrete-cosine transform basis functions to model low-frequency scanner drifts. In addition, cardiac and respiratory sources of nuisance were based on acquired physiological data and modeled with RETROICOR (Glover et al., 2000) using the MATLAB PhysIO toolbox (Kasper et al., 2017) in TAPAS (Frässle et al., 2021). For one subject with missing physiological data, the same number of fMRIPRrep’s anatomic component-based noise correction (aCompCor; Behzadi et al., 2007) regressors were entered in the model instead. The modeled data were obtained via linear regression and normalized. Finally, the two residual runs were concatenated and registered to the ICBM 152 Nonlinear Assymetrical template version 2009c (MNI152Nlin2009cAsym) (Fonov et al., 2009) using the nonlinear registration tool in antsRegistration (Avants et al., 2008) with the transformation parameters provided by fMRIPrep.
Experimental design and ROIs
Two runs of 15 min eyes-open wakeful rest (fixation on centered cross) were collected together with structural (or anatomical) scans during the first of four sessions that were part of a larger multisession 7 T study. The anatomical and experimental data acquired during the other sessions are not part of this study. Figure 1 provides an overview of the analysis, extending the data-driven echo approach (Leech et al., 2012; Braga et al., 2013) to the subcortex. With this multivariate technique, unique FC patterns are estimated while controlling for other subsignals within a region, revealing a more subtle subregional functional organization beyond a region’s global connectivity profile that remains concealed with univariate analyses (Leech et al., 2012).
Overview of the data analysis.
Fourteen subcortical ROIs were defined based on open-source parcellations (Table 1; Fig. 2a). Binary ROI masks were computed from the Multi-contrast Anatomical Subcortical Structure Parcellation algorithm (MASSP) (Bazin et al., 2020) that is based on quantitative MRI data (N = 105, ages 18-80) from the 7 T Amsterdam ultra-high field adult lifespan database (Alkemade et al., 2020) in high-resolution MNI space (MNI152Nlin2009bAsym) (Fonov et al., 2009). The MASSP parcellations include the Tha, Str, Cl, GPe, GPi, SN, STN, VTA, RN, Amg, PAG, and PPN. The LC was defined with the 7 T Probabilistic LC Atlas based on 53 healthy adults aged 52-84 years (Ye et al., 2021). In addition, the 17-network cortical parcellation (Yeo et al., 2011) was used for extracting a mask of the HPC, which was taken from the Default C network. To validate the results for noncortical structures, we also assessed whether we could reproduce the pattern of echoes within various cortical regions, including the PCC, medial prefrontal cortex (mPFC), and visual cortex (Braga et al., 2013). We used the same cortical network parcellation to derive masks for the striate and extrastriate cortex (Visual Central network) and the PCC and mPFC (Default A network). For bilateral ROIs, left and right hemispheres were combined into a single binary mask, and all masks were resampled to the resolution of the functional data with FLIRT using nearest-neighbor interpolation (version 6.0) (Jenkinson and Smith, 2001). The probabilistic LC mask was thresholded liberally so that voxels that overlapped 1% or more were included in the resampled mask.
Parcellation details for ROIs
Parcellations of subcortical ROIs and reference networks. a, Subcortical ROIs defined with open-source atlases and b data-driven reference networks from a whole-brain canonical ICA on the resting-state timeseries, labeled according to their maximum spatial correlation with a 17-network cortical parcellation. Corresponding whole-brain tSNR maps are shown in Extended Data Figure 2-1. SomA/B, Somatomotor A/B; ConA/B/C, Control A/B/C; TemPar, Temporal Parietal; DefA/B, Default A/B; VisC, Visual Central; VisP, Visual Peripheral; LimA/B, Limbic A/B; SalA/B, Salience/Ventral Attention A/B.
Figure 2-1
Whole-brain temporal signal to noise ratio (tSNR). For each of the two fMRI runs, voxel-wise tSNR values were calculated as the ratio of the mean and standard deviation of the resting-state timeseries after temporal high-pass filtering (1/128s) to remove low-frequency signal drifts. Individual tSNR maps (n = 40) were registered to standard MNI space (MNI152Nlin2009cAsym) with ANTs before voxel-wise tSNR values were averaged across subjects and runs to create the group-level map. The black contours outline the regions of interest that were included in the study. Download Figure 2-1, TIF file.
Statistical analysis
The individual preprocessed resting-state timeseries were masked with each of the binary ROIs and decomposed into 10 spatiotemporal independent subregions with a spatially restricted group canonical independent component analysis (canICA) as implemented in Nilearn. Although the temporal concatenation ICA approach is a popular technique in combination with dual regression, biases in the estimation of group-level networks may arise with varying degrees of interindividual variability (Hu and Yang, 2021). Instead, canICA applies a hierarchical approach in which individual data are decomposed before canonical correlation analysis to identify group commonalities (Varoquaux et al., 2010). The ROI-wise canICAs were restricted to find 10 independent components. Model order selection constitutes a main challenge in ICA, and the exact number of underlying signals in the diverse subcortical structures remains unknown. While prior analyses on the PCC demonstrated qualitatively similar outcomes for various model orders (Leech et al., 2012), conducting such comprehensive comparisons for all included structures was beyond the scope of this study. Instead, we opted to follow previous approaches and fix the number of components, addressing interregional differences in network echoes rather than precise dimensionality of individual structures.
Following spatiotemporal decomposition, the unique whole-brain FC of each independent component (subregion) was then investigated with dual regression (Beckmann et al., 2009; Zuo et al., 2010). First, the 10 spatial maps from the canICA were regressed onto every individual’s whole-brain resting-state data to estimate the subject-specific time course for each subregion. By simultaneously entering all 10 spatial maps as design matrix, the time course for each subregion was estimated while statistically controlling for the variance in the other subregions’ time courses. Second, the 10 subject-specific independent time courses were regressed onto the subject’s resting-state data to obtain spatial maps corresponding to the whole-brain, voxel-wise unique FC of each subregion. These subject-level FC maps were then combined in a nonparametric group-level analysis using random permutation testing (5000 permutations) with threshold-free cluster enhancement. This resulted in one group-level t-statistical map for each of the 10 subregions within each individual ROI that was thresholded with family-wise error correction at p < 0.05.
To quantify the presence of echoes from canonical resting-state networks within subcortical regions, the thresholded group-level FC maps were spatially correlated with data-driven reference networks obtained from a canICA on the whole-brain timeseries restricted to find 20 independent components. Based on visual inspection and low spatial Pearson product-moment correlation coefficients with an established 17-network cortical parcellation (Yeo et al., 2011), four independent components (r = 0.05, r = 0.04, r = 0.13, r = 0.04) were identified as artifactual and removed from further analysis. The resulting 16 reference networks were masked with the cortical network parcellation to remove any voxels located outside cortical gray matter (e.g., cerebral white matter, subcortex, CSF).
The extent of the spatial correlation between the FC map for each subregion and the reference networks was used to identify whether patterns of cortical network activity were mirrored, or echoed, in the unique subregional time courses.
Results
Data-driven networks correspond to existing cortical network parcellations
The 16 data-driven reference networks were labeled automatically according to their maximum spatial correlation with the well-established 17-network cortical parcellation (Yeo et al., 2011; Fig. 2b), which is based on rs-fMRI data from 1000 individuals. Despite large differences in field strength, data resolution, and parcellation method, we found correlation coefficients ranging from 0.21 to 0.67 (mean r = 0.44, SD = 0.14), generally indicating moderate to good spatial overlap with their reference network counterparts (Fig. 2b, bottom right): Somatomotor A (r = 0.66), Somatomotor B (r = 0.30), Control A (r = 0.46), Control B (r = 0.51), Control C (r = 0.57), Salience/Ventral Attention A (r = 0.34), Salience/Ventral Attention B (r = 0.54), Temporal Parietal (r = 0.21), Dorsal Attention A (DorA, r = 0.46), Dorsal Attention B (DorB, r = 0.25), Default A (r = 0.42), Default B (r = 0.47), Limbic A (r = 0.30), Limbic B (r = 0.41), Visual Central (r = 0.49), and Visual Peripheral (r = 0.67). The data-driven Temporal Parietal network also partially overlapped with the Control A network parcellation (r = 0.15).
Together, the reference networks covered 66% of cortical gray matter defined in the parcellation by Yeo et al. (2011). The strongest deviation was observed in the anterior temporal cortex, which was not remedied by increasing model order (40 or 100 independent components) or a cortically restricted canICA. To assess corresponding variations in temporal SNR (tSNR), we calculated voxel-wise tSNR values as the ratio of the mean and SD of the resting-state timeseries after temporal high-pass filtering (1/128 s). Individual tSNR maps were registered to standard MNI space and averaged (voxel-wise) across subjects and runs. Compared with other cortical areas, reduced tSNR in the temporal lobe was observed; and as a consequence, temporal networks were underrepresented in the analysis (Extended Data Fig. 2-1).
Subcortical structures echo signals from different resting-state networks
The 10 thresholded FC maps for each ROI, representing the unique whole-brain FC of each subregion at group level, were spatially correlated with the 16 unthresholded spatial maps of the data-driven reference networks. Figure 3a summarizes the degree of network echoes for the nine ROIs that demonstrated at least one spatial correlation with any reference network above a threshold that was arbitrarily set at the 97th percentile of all spatial correlations (r = 0.16). Echoes were summarized by counting above-threshold spatial correlations in terms of (1) the number of reference networks represented in each ROI and (2) the number of subregions that echoed a reference network. For example, six distinct striatal subregions displayed FC profiles that spatially correlated above-threshold with in total 10 different resting-state networks. Figure 3b presents the actual maximum spatial correlations between each ROI and each reference network, independent of subregion. The reference network that was represented most often was the Salience B network, correlating above-threshold with seven ROIs, followed by Default A, Control C, and Visual Peripheral, each with at least one above-threshold spatial correlation with six different ROIs.
Echoes of intrinsic connectivity networks in the subcortex. a, The number of distinct subregions within an ROI with an FC profile that resembled a reference network (Subregions) and the number of different reference networks that were echoed within a region (Networks) both defined by counting above-threshold spatial correlations. b, The maximum spatial correlation between each ROI and each reference network, independent of subregion, for nine ROIs that demonstrated at least one above-threshold spatial correlation to any reference network. Subregional FC profiles for a subset of structures and their spatial correlation with reference networks are illustrated in Extended Data Figure 3-1. The same analysis was repeated with reference networks taken from the 17-network cortical parcellation (Yeo et al., 2011) shown in Extended Data Figure 3-2 as well as for three cortical ROIs (Extended Data Fig. 3-3). SomA/B, Somatomotor A/B; ConA/B/C, Control A/B/C; TemPar, Temporal Parietal; DefA/B, Default A/B; VisC, Visual Central; VisP, Visual Peripheral; LimA/B, Limbic A/B; SalA/B, Salience/Ventral Attention A/B.
Figure 3-1
Functional connectivity patterns of subcortical subregions and their spatial overlap with intrinsic connectivity networks. Diversity in whole-brain functional connectivity (FC) of distinct subregions of subcortical structures plotted on cortical surface meshes and the maximum spatial correlation with data-driven reference networks (four out of sixteen networks shown for illustration). Although the spatial correlations are calculated from the unthresholded spatial maps, the reference networks were thresholded by assigning each voxel to its most strongly associated network based on the group canICA (i.e., every voxel is assigned to only one network and networks are non-overlapping) for illustration purposes. The subregion-specific FC maps are the group-level results of a dual regression analysis on the time course for each subregion while controlling for the variance in the other subregions, statistically tested with random permutation testing and thresholded at p < .05. Labels: thalamus (Tha), striatum (Str), claustrum (Cl), hippocampus (HPC), subsantia nigra (SN), globus pallidus externa (GPe), Default A (DefA), Default B (DefB), Somatomotor A (SomA), Salience/Ventral Attention A (SalA). Download Figure 3-1, TIF file.
Figure 3-2
Echoes of well-established cortical intrinsic connectivity networks in the subcortex. (a) The number of distinct subregions within a region of interest (ROI) with a functional connectivity profile that resembled a reference network (‘Subregions’) and the number of different reference networks that were echoed within a region (‘Networks’) both counted as the number of above-threshold spatial correlations. Reference networks were taken from the 17-network cortical parcellation (Yeo et al, 2011). (b) The maximum spatial correlation between each ROI and each reference network, independent of subregion, for nine ROIs that demonstrated at least one above-threshold spatial correlation. Labels: thalamus (Tha), striatum (Str), globus pallidus externa (GPe), claustrum (Cl), hippocampus (HPC), amygdala (Amg), substantia nigra (SN), ventral tegmental area (VTA), pedunculopontine nucleus (PPN), Somatomotor A (SomA), Somatomotor B (SomB), Control A (ConA), Control B (ConB), Control C (ConC), Temporal Parietal (TemPar), Dorsal Attention A (DorA), Dorsal Attention B (DorB), Default A (DefA), Default B (DefB), Default C (DefC), Visual Central (VisC), Visual Peripheral (VisP), Limbic A (LimA), Limbic B (LimB), Salience/Ventral Attention A (SalA), Salience/Ventral Attention B (SalB). Download Figure 3-2, TIF file.
Figure 3-3
Echoes of intrinsic connectivity networks in cortical regions of interest. (a) The maximum spatial correlation between the whole-brain functional connectivity (FC) of each cortical ROI with data-driven reference networks. The results demonstrate greater functional heterogeneity within posterior cingulate cortex (PCC) and medial prefrontal cortex (mPFC), as evident in more distributed patterns of FC with default mode, control, and salience networks compared to the visual cortex (VC), which showed a more uniform organization dominated by a preferential connection with the visual peripheral network. This is consistent with previous work (Braga et al, 2013) and provides a validation for our novel application of the multivariate analysis within subcortical regions of interest. (b) The results of an identical analysis but with the 17-network cortical parcellation (Yeo et al, 2011) as reference networks, revealing a less pronounced but qualitatively similar pattern of results compared to the data-driven networks. Labels: Somatomotor A (SomA), Somatomotor B (SomB), Control A (ConA), Control B (ConB), Control C (ConC), Temporal Parietal (TemPar), Dorsal Attention A (DorA), Dorsal Attention B (DorB), Default A (DefA), Default B (DefB), Default C (DefC), Visual Central (VisC), Visual Peripheral (VisP), Limbic A (LimA), Limbic B (LimB), Salience/Ventral Attention A (SalA), Salience/Ventral Attention B (SalB). Download Figure 3-3, TIF file.
Seven subcortical ROIs echoed signals from more than one network, including the following: Tha, Str, HPC, Cl, GPe, SN, and VTA. The former four ROIs furthermore showed that the echoes from different reference networks were distributed among multiple subregions, indicating evidence for a heterogeneous functional organization. In contrast, both the Amg and PPN showed medium and small spatial correlations, respectively, with only one reference network (Amg: r = 0.37 [DefA]; PPN: r = 0.19 [SalB]). The GPi, STN, RN, PAG, and LC failed to show evidence of echoes as none of their subregions demonstrated a connectivity pattern that resembled the pattern of an intrinsic connectivity network. In some cases, a subregion’s FC profile was widespread and shared spatial similarity with more than one reference network. Extended Data Figure 3-1 presents a few FC maps to illustrate the diversity and similarity in connectivity profiles to different reference networks across a subset of subcortical structures.
The FC maps of each subregion were also spatially correlated with the 17-network cortical parcellation (Yeo et al., 2011), which yielded generally lower spatial correlations but a qualitatively similar pattern of results (Extended Data Fig. 3-2). To validate these novel results for the subcortex, we repeated the analyses for three cortical regions that were previously investigated. Results for the PCC, mPFC, and visual cortex are presented in Extended Data Figure 3-3 and are largely consistent with previous findings (Leech et al., 2012; Braga et al., 2013).
Topographic organization of functionally heterogeneous subcortical structures
Figure 4 shows the topographic pattern of network echoes in the subregions of the seven ROIs with more than one above-threshold spatial correlation. Subregions are color-coded according to the reference network they echoed most strongly, whereas subregions with a maximum spatial correlation below threshold (r < 0.16) are translucent. For every ROI, there were several subregions that did not mirror the activity in any intrinsic connectivity network, because they were predominantly functionally connected to other subcortical structures or because their signal largely reflected noise on visual inspection.
Topography of network echoes within heteromodal subcortical structures. Spatiotemporal decomposition of subcortical structures into independent subregions, color-coded according to their strongest network echo or made translucent if their maximum spatial correlation with any reference network did not reach threshold. SomA/B, Somatomotor A/B; ConA/B/C, Control A/B/C; TemPar, Temporal Parietal; DefA/B, Default A/B; VisC, Visual Central; VisP, Visual Peripheral; LimA/B, Limbic A/B; SalA/B, Salience/Ventral Attention A/B.
Five thalamic subregions echoed signals from various reference networks, demonstrating a heterogeneous organization that was mostly symmetrically distributed in bilateral subdivisions. Left and right ventromedial subregions were both most strongly correlated to the Somatomotor A network (left: r = 0.26, right: r = 0.20), although the right subregion’s connectivity profile also spatially overlapped with Salience B (r = 0.20). A more dorsomedial bilateral subregion displayed a connectivity pattern that correlated with the pattern of multiple reference networks, including Default A (r = 0.38), Default B (r = 0.32), and Control A (r = 0.25). Another bilateral subregion, more dorsolaterally located, correlated most strongly with the DorA network (r = 0.31), although there was also spatial overlap with Somatomotor A (r = 0.29), DorB (r = 0.25), and Visual Peripheral (r = 0.24) networks. Finally, the Default B network was represented in the posterior part of the left-sided Tha (r = 0.22).
Within the Str, there were six different subregions that echoed one or more reference networks, located mostly within the caudate nucleus. A subregion primarily in the left tail of the caudate nucleus spatially correlated with the Default B network (r = 0.21), whereas a subregion covering more of the right tail of caudate nucleus most strongly echoed Control B (r = 0.26), although its widespread connectivity pattern also overlapped with Temporal Parietal (r = 0.23) and Salience A (r = 0.22) networks. A bilateral subregion covering the NAc correlated most strongly with Default A (r = 0.40), whereas another bilateral subregion in the mediodorsal part of the caudate head was functionally connected with Control A (r = 0.26) and Default B (r = 0.21) networks. Subregions that most strongly echoed the Salience A network included a division in the posterior parts of the left caudate tail and left putamen (r = 0.20) as well as a bilateral region in the lateral NAc (r = 0.19).
For the HPC, we observed that different intrinsic connectivity networks were echoed within four different subregions. In the left hemisphere, a posterior dorsal subregion correlated most strongly with Default A (r = 0.34), whereas a more ventrally located subregion correlated exclusively with the Limbic A network (r = 0.24). A bilateral anteromedial subregion was functionally connected to the Visual Central network (r = 0.21), whereas a posterior dorsal subregion in the right hemisphere echoed the Visual Peripheral (r = 0.30) as well as the Dorsal Attention networks (DorA: r = 0.28, DorB: r = 0.30).
Five subregions of the Cl showed an FC profile that correlated with different reference networks. A small, bilateral subregion in the ventral Cl had a widespread cortical connectivity that had the strongest spatial similarity with DorA (r = 0.23), but also Somatomotor A (r = 0.20), Dorsal attention B (r = 0.19), and Salience B (r = 0.19) networks. Left and right subdivisions in the posterior part both echoed the Salience A network (r = 0.26 and r = 0.21, respectively). In addition, an exclusive functional connection with the Default B network was observed in an anterior subregion of the left Cl (r = 0.32) and with the Somatomotor A network in a more posterior subregion of the right Cl (r = 0.35).
The GPe and SN each had one subregion with a widespread connectivity profile comprising seven and three reference networks, respectively (Fig. 3b). In the GPe, a bilateral dorsolateral subdivision most strongly echoed the Somatomotor A network (r = 0.26), but its signal also correlated with activity in DorA (r = 0.23) and Control networks A and B (r = 0.20 and r = 0.22, respectively). The most pronounced network echo within the SN was from Default A (r = 0.24) and came from a bilateral subregion in the medial anterior SN. The same subregion also showed traces from Salience B (r = 0.22) and Control C (r = 0.16) networks. For the VTA, a large inferomedial subdivision in the right hemisphere was most strongly connected to Salience B (r = 0.19) and just below threshold to Visual Peripheral (r = 0.15) networks. Echoes from the Default A network were furthermore present in two other subregions of the VTA, but spatial correlations were weaker (r = 0.15 and r = 0.13).
Discussion
Despite accumulating insights into the mechanisms of functional integration within the cortex, subcortical substrates of cross-network convergence are largely unexplored. Nonetheless, the subcortex is embedded within an extensive cortico-subcortical architecture that is thought to serve integrative rather than purely segregated functions (Haber, 2003). Here, we aimed to more closely examine the underlying functional organization of subcortical nuclei and their subregional connectivity to functional networks across the cortex.
Consistent with our expectations, we show that individual subcortical structures contain a composite of neural signals that can be decomposed into activity traces of intrinsic network activity. In their study, Braga et al. (2013) showed that activity in multiple networks converges at specific transmodal zones in the cortex, as reflected in a mixture of signals that partially correlate with different networks. We demonstrate that this property is not limited to cortical regions by revealing potential mechanisms for multinetwork integration in the subcortex. The results provide the strongest evidence for functional heterogeneity within the Tha, Str, Cl, and HPC, for which we observed a complex pattern of subregional whole-brain FC that resembled spontaneous activity in distinct functional networks. Subregions in left and right hemispheres had similar spatiotemporal signatures that echoed the same functional networks, showing a symmetrical bilateral topography that is consistent with prior work (Cheng and Liu, 2021).
The Tha and Str are the most commonly represented noncortical structures in studies of global brain connectivity, providing support for their putative role as hub regions (Van der Heuvel and Sporns, 2011; Bell and Shine, 2015, 2016). Whereas several studies report an amalgamation of primarily sensory information within thalamic subregions consistent with its gating function (Tomasi and Volkow, 2011; Ji et al., 2019), we observed traces of somatomotor as well as default mode and dorsal attention networks. The somatomotor subdivisions also spatially overlapped with cingulo-opercular regions of the salience network, which aligns with findings of a “motor integration zone” within ventral thalamic nuclei (Greene et al., 2020). Additionally, dorsal attention, somatomotor, and visual networks converged in a dorsolateral subregion, similar albeit slightly less posterior to the “visual integration zone” in the pulvinar nucleus reported earlier (Greene et al., 2020). For the Str, we observed signal echoes from default mode, control, and salience networks predominantly within the caudate head and left tail, right tail, and left putamen, respectively. Despite large methodological differences across studies, these findings are consistent with prior evidence for “cognitive” integration within the Str (Choi et al., 2012; Greene et al., 2020; Seitzman et al., 2020) and supports thalamic and striatal roles in information integration and higher-level cognitive functioning (Haber, 2003; Hwang et al., 2017).
Although organizational principles may broadly concur, precise functional boundaries and network connections diverge across studies. For example, the subregional profiles identified here partially deviate from another data-driven co-partitioning (Cheng and Liu, 2021) and a voxel-wise winner-take-all approach (Seitzman et al., 2020) for the Tha, as well as the from the striatal architecture reported by Choi et al. (2012). Additionally, we found interhemispheric differences in the HPC (i.e., visual and dorsal attention network echoes in the right and default mode and limbic in the left side) that are inconsistent with reports of lateralized subdivisions along an anterior-posterior axis, as well as the location along this axis of the preferential connection to the default mode network (Blessing et al., 2016; Cheng et al., 2020; Ezama et al., 2021). Given differences in connectivity with entorhinal and parahippocampal cortex (Qin et al., 2016; Seoane et al., 2022), it is possible that the extent of hippocampal and surrounding voxels included in the analysis explains some of the discrepancies across studies, which might be further exacerbated by the effects of spatial smoothing. Furthermore, high degrees of individual variability in subcortical anatomy and FC may result in distortions of group-level estimations (De Hollander et al., 2015; Greene et al., 2020; Sylvester et al., 2020; Tian et al., 2020; Marek and Greene, 2021).
Similar to previous observations for the cortex (Braga et al., 2013), we demonstrate that functional heterogeneity is not ubiquitously present throughout the subcortex. Within the GPe, SN, and VTA, only one subregion’s connectivity profile resembled patterns of functional network activity. A region in the dorsolateral GPe echoed somatomotor as well as dorsal attention and control networks, indicating an integrative site that may support its known role in voluntary, planned movement. Both the SN and VTA showed a pattern of converging signals from default mode and salience networks, although less evident in the VTA. Whereas this association with the default mode network is more established (Bär et al., 2016; Zhang et al., 2016; Edlow, 2021; Li et al., 2021), connectivity to the salience network is less known and may indicate involvement in attention and spontaneous cognition (O’Callaghan et al., 2021).
No clear evidence for functional integration was observed for the Amg and PPN. Whereas the PPN likely takes part in more specialized subcortical circuitry involved in arousal and locomotion (Martinez-Gonzales et al., 2011; Bennarroch, 2013), the Amg was previously proposed as hub structure (Tomasi and Volkow, 2011) and showed dissociable FC profiles from its separate nuclei (Kerestes et al., 2017). Although we did not find evidence for such heterogeneity when controlling for other subregional time courses, we observed an intact connection with the default mode network, which is supported by other work (Kerestes et al., 2017; Sylvester et al., 2020; Harrison et al., 2021). For the remaining structures (i.e., GPi, STN, RN, PAG, and LC), we failed to find network echoes. Although previous univariate FC studies have indicated correlations with widespread cortical activity for some of these structures (e.g., Zhang et al., 2016; Anteraper et al., 2018), the multivariate analysis here did not result in a clear group-level pattern of cortical connectivity. Similar to the PPN, these structures may be less involved in integrating spontaneous signals from distributed functional processes across the cortex, but are likely more strongly embedded in local networks to support segregated functional processing (Singh et al., 2022). Recent findings suggest that neuromodulatory nuclei for dopaminergic and noradrenergic systems are driving systems-level integration and cognition (Zhang et al., 2016; De Gee et al., 2017; Liu et al., 2017). However, not all findings converge. For example, Bär et al. (2016) showed that LC connectivity to the default mode network disappeared when controlling for adjacent neural signals and that hub-like features of midbrain nuclei were not supported by a graph theory analysis. The results presented here align with this observation and emphasize that integrative properties of these structures, among which the LC, remain somewhat elusive. Given proposed roles of the LC in mediating the dynamics of cortical connectivity and neural gain (Aston-Jones and Cohen, 2005; Munn et al., 2021), it is perhaps not surprising that no dissociable traces of functional network activity are observed. That is, the LC may drive global states of network integration and segregation rather than serving as a convergence zone in itself.
In conclusion, our results suggest that subcortical structures exhibit varying degrees of functional heterogeneity. This characteristic might be expressed along a gradient, where structures adjacent to the cortex seem more likely to support multinetwork integration compared with deep brain nuclei. However, several factors may confound interpretations of interregional differences in the subcortex. For example, deep brain nuclei are generally smaller in size and have weaker SNR, while subcortex near the cortex is susceptible to signal bleeding from adjacent cortical voxels, to which they are also reciprocally connected (Choi et al., 2012). This issue might be especially prominent in the Cl, which is a thin sheet-like structure situated directly between the Str and insula. In a recent study, Krimmel et al. (2019) used a novel regression technique on similar high-resolution fMRI data (1.5 mm isotropic voxels) to isolate the signal in the Cl from nearby cortical and striatal voxels, which preserved the widespread FC with cortical networks involved in attention and cognitive control. Although we did not correct for potential signal bleeding beyond limiting the amount of spatial smoothing, our finding of functionally heterogeneous network echoes within the Cl’s subdivisions coincides with this work and its postulated role in attention and cognition (Bell and Shine, 2015; Krimmel et al., 2019; J. B. Smith and Jackson, 2020).
It should be noted that recent work highlights the difference in FC between eyes-open and eyes-closed resting-state conditions, particularly with regard to internetwork connectivity of visual and sensorimotor networks to default mode and salience networks (Agcaoglu et al., 2019; Costumero et al., 2020; Han et al., 2023). While a large portion of studies on subcortical connectivity cited here are correspondingly based on eyes-open resting-state fMRI (e.g., Choi et al., 2012; Blessing et al., 2016; Hwang et al., 2017; Greene et al., 2020; Seitzman et al., 2020; Sylvester et al., 2020), future efforts could contrast our results to potential reconfigurations during other resting-state and experimental conditions. Investigating changes in the pattern of echoes according to external factors, such as cognitive demand, and internal state are likely necessary to illuminate their functional relevance (e.g., Leech et al., 2012).
Although the precise significance of network echoes for cognition and behavior is not resolved, we strengthen the evidence that the subcortex participates in cross-network integration through echoing intrinsic network activity. These results may ignite new intriguing hypotheses on the mechanisms of spontaneous cognitive processes, such as mind wandering (Mittner et al., 2016; Zuberer et al., 2021). Previous work has shown that mind wandering correlates with activity and connectivity in the default mode and frontoparietal control networks as well as the subcortex (Mittner et al., 2014; Kucyi et al., 2017; Groot et al., 2022). Given that both subtle and pronounced reorganizations in FC occur with changes in task demand (Leech et al., 2012; Braga et al., 2013; Tian et al., 2020), investigations of how the complex pattern of echoes in the subcortex is perturbed by attentional changes may reveal novel insights into the mechanisms that drive mind wandering.
Footnotes
This work was supported by The Netherlands Organization of Scientific Research Grant 016.Vici.185.052 to B.U.F.
The authors declare no competing financial interests.
- Correspondence should be addressed to Matthias Mittner at matthias.mittner{at}uit.no