Individualized Functional Subnetworks Connect Human Striatum and Frontal Cortex

Subjects

Data were collected from 10 healthy, right-handed, young adult subjects (5 females; age: 24–34). Two of the subjects are authors (NUFD and SMN), and the remaining subjects were recruited from the Washington University community. Informed consent was obtained from all participants. The study was approved by the Washington University School of Medicine Human Studies Committee and Institutional Review Board. Other findings using these participants have been previously reported in (Gordon, Laumann, Gilmore, et al. 2017c; Gratton et al. 2018; Marek et al. 2018; Greene et al. 2020; Sylvester et al. 2020).

MRI Image Acquisition

Imaging for each subject was performed on a Siemens TRIO 3 T MRI scanner over the course of 12 sessions conducted on separate days, each beginning at midnight. Structural MRI was conducted across two separate days. In total, four T1-weighted images (sagittal, 224 slices, 0.8 mm isotropic resolution, TE = 3.74 ms, TR = 2400 ms, TI = 1000 ms, flip angle = 8 degrees), four T2-weighted images (sagittal, 224 slices, 0.8 mm isotropic resolution, TE = 479 ms, TR = 3200 ms), four MRA (transverse, 0.6 × 0.6 × 1.0 mm, 44 slices, TR = 25 ms, TE = 3.34 ms), and eight MRVs, including four in coronal and four in sagittal orientations (sagittal: 0.8 × 0.8 × 2.0 mm thickness, 120 slices, TR = 27 ms, TE = 7.05 ms; coronal: 0.7 × 0.7 × 2.5 mm thickness, 128 slices, TR = 28 ms TE = 7.18 ms), were obtained for each subject. Analyses of the MRA and MRV scans are not reported here.

On 10 subsequent days, each subject underwent 1.5 h of functional MRI scanning beginning at midnight. In each session, we first collected thirty contiguous minutes of resting-state fMRI data, in which subjects visually fixated on a white crosshair presented against a black background. Each subject was then scanned during performance of three separate tasks, which are not examined here. Across all sessions, each subject was scanned for 300 total minutes during the resting state. All functional imaging was performed using a gradient-echo EPI sequence (TR = 2.2 s, TE = 27 ms, flip angle = 90°, voxel size = 4 mm × 4 mm × 4 mm, 36 slices). In each session, one gradient echo field map sequence was acquired with the same prescription as the functional images. An EyeLink 1000 eye-tracking system (http://www.sr-research.com) allowed continuous monitoring of subjects’ eyes in order to check for periods of prolonged eye closure, potentially indicating sleep. Only one subject (MSC08) demonstrated prolonged eye closures.

Cortical Surface Generation

Generation of cortical surfaces from the MRI data followed a procedure similar to that previously described in (Marek et al. 2018). First, anatomical surfaces were generated from the subject’s average T1-weighted image in native volumetric space using FreeSurfer’s default recon-all processing pipeline (version 5.3). This pipeline first conducted brain extraction and segmentation. After this step, segmentations were hand-edited to maximize accuracy. Subsequently, the remainder of the recon-all pipeline was conducted on the hand-edited segmentations, including generation of white matter and pial surfaces, inflation of the surfaces to a sphere, and surface shape–based spherical registration of the subject’s original surface to the fsaverage surface (Dale et al. 1999; Fischl et al. 1999). The fsaverage-registered left and right hemisphere surfaces were brought into register with each other using deformation maps from a landmark-based registration of left and right fsaverage surfaces to a hybrid left–right fsaverage surface (“fs_LR”; Van Essen et al. 2012). These fs_LR spherical template meshes were input to a flexible multimodal surface matching (MSM) algorithm using sulc features to register templates to the atlas mesh (Robinson et al. 2014). These newly registered surfaces were then down-sampled to a 32 492 vertex surface (fs_LR 32k) for each hemisphere. The various structural metric data (thickness, curvature, etc.) from the original surfaces to the fs_LR 32k surface were composed into a single deformation map allowing for one-step resampling. MSM registration provided a more optimal fit of pial and white surfaces and reduced areal distortion (Glasser et al. 2016). These various surfaces in native stereotaxic space were then transformed into atlas space (711-2B) by applying the previously calculated T1-to-atlas transformation.

fMRI Preprocessing

Functional data were preprocessed to reduce artifacts and to maximize cross-session registration. All sessions underwent correction of odd versus even slice intensity differences attributable to interleaved acquisition, intensity normalization to a whole-brain mode value of 1000, and within run correction for head movement. Atlas transformation was computed by registering the mean intensity image from a single BOLD session to Talairach atlas space (Talairach and Tournoux 1988) via the average high-resolution T2-weighted image and average high-resolution T1-weighted image. All subsequent BOLD sessions were linearly registered to this first session. This atlas transformation, mean field distortion correction (see below), and resampling to 2-mm isotropic atlas space were combined into a single interpolation using FSL’s applywarp tool (Smith et al. 2004). All subsequent operations were performed on the atlas-transformed volumetric time series.

Distortion Correction

A mean field map was generated based on the field maps collected in each subject (Laumann et al. 2015). This mean field map was then linearly registered to each session and applied to that session for distortion correction. To generate the mean field map, the following procedure was used: 1) Field map magnitude images were mutually coregistered. 2) Transforms between all sessions were resolved. Transform resolution reconstructs the n-1 transforms between all images using the n(n − 1)/2 computed transform pairs. 3) The resolved transforms were applied to generate a mean magnitude image. 4) The mean magnitude image was registered to an atlas-representative template. 5) Transforms from individual session magnitude images to atlas space were computed by composing the session-to-mean and mean-to-atlas transforms. 6) Phase images were then transformed to atlas space using the composed transforms, and a mean phase image in atlas space was computed.

Application of mean field map to individual fMRI sessions: 1) For each session, field map uncorrected data were registered to atlas space, as above. 2) The generated transformation matrix was then inverted and applied to the mean field map to bring the mean field map into the session space. 3) The mean field map was used to correct distortion in each native-space run of resting state and task data in the session. 4) The undistorted data were then reregistered to atlas space. 5) This new transformation matrix and the mean field map then were applied together to resample each run of resting state and task data in the session to undistorted atlas space in a single step.

RSFC Preprocessing

Additional preprocessing steps to reduce spurious variance unlikely to reflect neuronal activity were executed as recommended in (Power et al. 2014; Ciric et al. 2017). First, temporal masks were created to flag motion-contaminated frames. We observed that two subjects (MSC 03 and MSC 10) had a high-frequency artifact in the motion estimates calculated in the phase encode (anterior–posterior) direction that did not appear to reflect biological movement. We thus filtered the motion estimate time courses in this direction only to retain effects occurring below 0.1 Hz in all subjects for consistency (Gratton et al. 2020). Motion-contaminated volumes were then identified by frame-by-frame displacement (FD). Frames with FD > 0.2 mm were flagged as motion-contaminated. Across all subjects, these masks censored 28% ± 18% (range: 6–67%) of the data; on average, subjects retained 5929 ± 1508 volumes (range: 2733–7667), corresponding to 217 ± 55 min (range: 100–281).

After computing the temporal masks for high motion frame censoring, the data were processed with the following steps: 1) demeaning and detrending, 2) linear interpolation across censored frames using so that continuous data can be passed through, 3) a band-pass filter (0.005 Hz < f < 0.01 Hz) without reintroducing nuisance signals (Hallquist et al. 2013) or contaminating frames near high motion frames (Carp 2013).

Next, the filtered BOLD time series underwent a component-based nuisance regression approach (Marek et al. 2018). Nuisance regression using time series extracted from white matter and cerebrospinal fluid (CSF) assumes that variance in such regions is unlikely to reflect neural activity. Variance in these regions is known to correspond largely to physiological noise (e.g., CSF pulsations), arterial pCO₂-dependent changes in T2^*-weighted intensity, and motion artifact; this spurious variance is widely shared with regions of interest in gray matter. We also included the mean signal averaged over the whole brain as a nuisance regressor. Global signal regression (GSR) has been controversial. However, the available evidence indicates that GSR is a highly effective denoising strategy (Power et al. 2015; Ciric et al. 2017).

Nuisance regressors were extracted from white matter and ventricle masks, first segmented by FreeSurfer (Fischl 2012), and then spatially resampled in register with the fMRI data. Voxels surrounding the edge of the brain are particularly susceptible to motion artifacts and CSF pulsations (Satterthwaite et al. 2013; Patriat et al. 2015); hence, a third nuisance mask was created for the extra-axial compartment by thresholding the temporal standard deviation image (SDt > 2.5%), excluding a dilated whole-brain mask. Voxel-wise nuisance time series were dimensionality reduced as in CompCor (Behzadi et al. 2007), except that the number of retained regressors, rather than being a fixed quantity, was determined, for each noise compartment, by orthogonalization of the covariance matrix and retaining components ordered by decreasing eigenvalue up to a condition number of 30 (max eigenvalue/min eigenvalue > 30). The retained components across all compartments formed the columns of a design matrix, X, along with the global signal, its first derivative, and the six time series derived by retrospective motion correction. The columns of X are likely to exhibit substantial colinearity. Therefore, to prevent numerical instability owing to rank-deficiency during nuisance regression, a second-level SVD was applied to XX^T to impose an upper limit of 250 on the condition number. This final set of regressors was applied in a single step to the filtered, interpolated BOLD time series, with censored data ignored during beta estimation. Finally, the data were upsampled to 2 mm isotropic voxels. Censored frames were then excised from the data for all subsequent analyses.

Surface Processing and CIFTI Generation of BOLD Data

Surface processing of BOLD data proceeded through the following steps. First, the BOLD fMRI volumetric time series were sampled to each subject’s original midthickness left- and right-hemisphere surfaces (generated as the average of the white and pial surfaces) using the ribbon-constrained sampling procedure available in Connectome Workbench 1.0. This procedure samples data from voxels within the gray matter ribbon (i.e., between the white and pial surfaces) that lie in a cylinder orthogonal to the local midthickness surface weighted by the extent to which the voxel falls within the ribbon. Voxels with a time series coefficient of variation 0.5 standard deviations higher than the mean coefficient of variation of nearby voxels (within a 5-mm sigma Gaussian neighborhood) were excluded from the volume to surface sampling, as described in (Glasser et al. 2013). Once sampled to the surface, timecourses were deformed and resampled from the individual’s original surface to the 32-k fs_LR surface in a single step using the deformation map generated above (in “Cortical Surface Generation”). This resampling allows point-to-point comparison between each individual registered to this surface space.

These surfaces were then combined with volumetric subcortical and cerebellar data into the CIFTI format using Connectome Workbench (Marcus et al. 2011), creating full brain timecourses excluding non–gray matter tissue. Subcortical (including accumbens, amygdala, caudate, hippocampus, pallidum, putamen, and thalamus) and cerebellar voxels were selected based on the FreeSurfer segmentation of the individual subject’s native-space average T1, transformed into atlas space, and manually inspected. Finally, the BOLD timecourses were smoothed with a geodesic 2D (for surface data) or Euclidean 3D (for volumetric data) Gaussian kernel of σ = 2.55 mm.

Regression of Adjacent Cortical Tissue from RSFC BOLD

Some striatal regions (e.g., lateral putamen) are in close anatomical proximity to cortex, potentially resulting in spurious functional coupling between the cortical vertices and adjacent subcortical voxels due to local autocorrelation. To reduce this effect, RSFC BOLD time series from all vertices falling within 20-mm Euclidean distance of a source voxel were averaged and then regressed from the voxel time series in order to remove the autocorrelation, following (Buckner et al. 2011; Marek et al. 2018; Greene et al. 2020).

Mapping Multiscale Network Structure

The network organization of each subject’s brain was delineated following (Gordon et al. 2020) using the graph-theory-based Infomap algorithm for community detection (Rosvall and Bergstrom 2008). In this approach, we calculated the cross-correlation matrix of the time courses from all brain vertices (on the cortical surfaces) and voxels (in subcortical structures), concatenated across sessions. Correlations between vertices/voxels within 30 mm of each other were set to zero in this matrix to avoid basing network membership on correlations attributable to spatial smoothing. Geodesic distance was used for within-hemisphere surface connections and Euclidean distance for subcortical-to-cortical connections. Connections between subcortical structures were disallowed, as we observed extremely high connectivities within nearly the entire basal ganglia that would prevent network structures from emerging. Interhemispheric connections between the cortical surfaces were retained, as smoothing was not performed across the midsagittal plane.

We observed that connectivity patterns within regions known to have low BOLD signal due to susceptibility artifact dropout (e.g., ventral anterior temporal lobe and portions of orbitofrontal cortex) were effectively random. To avoid having the delineated network structures distorted by such random connections, we first calculated a set of common low-signal regions as the vertices in which the average mode-1000 normalized BOLD signal across subjects and timepoints was less than 750 (as in Wig et al. 2014; Gordon et al. 2016). All connections to these low-signal regions were set to zero. These vertices converged with regions known to suffer severe signal loss due to susceptibility artifact and to the presence of veinous structures (see Supplementary Fig. 1 for excluded vertices).

This matrix was then thresholded to retain at least the strongest 0.1% of connections to each vertex and voxel (Gordon et al. 2020). Note that this thresholding approach differs from previous approaches for forming brain graphs from functional connectivity data (e.g., Power et al. 2011; Gordon, Laumann, Gilmore, et al. 2017c). The typical procedure applies a uniform edge density threshold to all functional connectivity values in the brain. The weakness of the uniform threshold approach is that subcortical structures generally have decreased BOLD signal-to-noise relative to cortex due to their greater distance from the MR head coil. The result is that functional connectivity patterns seeded from striatal voxels have weaker peak connectivity strengths, even though they may appear well organized and coherent with known cortical networks. Thus, with a uniform threshold, these regions are frequently not identified as being networked with cortical regions.

By contrast, the current approach thresholds the connectivity maps seeded from each cortical or subcortical point in the brain separately, always retaining at least the 0.1% strongest connections. The result is that the connectivity of all striatal voxels are represented in the brain graph used for community detection. Previous validation of this procedure (Gordon et al. 2020) has shown that the 0.1% density threshold identifies subnetwork divisions in resting data that best explain task activations in these subjects.

The thresholded matrices were used as inputs for the Infomap algorithm, which calculated community assignments separately for each threshold. The resulting communities represent subnetworks in the brain. Small networks with 10 or fewer vertices/voxels were considered unassigned and removed from further consideration. The above analysis was conducted in each individual subject.

Identifying Matched Striatal Subnetworks in Individuals

For each subject, we considered only subnetworks that had at least some representation in the striatum (the individually defined caudate, putamen, nucleus accumbens, and globus pallidus). Following (Gordon et al. 2020), we visually examined the cortical and subcortical topographies of each subnetwork with striatal representation, as well as their topological arrangement relative to each other. After careful consideration of the subnetworks observable in this population, matched subnetworks were identified based on the following heuristic rules:

A subnetwork in inferior ventral striatum usually also had representation in subgenual cingulate and/or orbitofrontal cortex.

A subnetwork in anterior ventral striatum had representation in bilateral pregenual cingulate cortex. This subnetwork was previously described in (Gordon et al. 2020).

A subnetwork in superior ventral striatum had representation in anterior cingulate cortex and middle insula.

A subnetwork that was spatially variable within striatum, with representation most often in medial putamen, but was consistently present in dorsal anterior cingulate cortex and dorsal insula.

A subnetwork in lateral middle putamen and posterior dorsomedial prefrontal cortex.

A subnetwork in dorsal caudate that had consistent representation as the dorsomedial cluster of the frontoparietal network and inconsistent representation in lateral frontal and parietal cortex.

A subnetwork in medial anterior putamen (especially left putamen) that had consistent cortical representation converging with the language network (Braga et al. 2020).

A subnetwork in lateral anterior caudate (especially left caudate) that had representation to a set of cortical regions directly anterior to the language network above.

A subnetwork (sometimes multiple subnetworks) in ventral posterior putamen that had representation somewhere in the central sulcus, with the location of the representation varying substantially by subject.

Additionally, we commonly observed subnetworks that were generally in posterior lateral putamen (but with variable representation across subjects) that had cortical representation only in middle insular cortex, physically adjacent to the striatal representation, or (in one subject) in unusual locations such as medial occipital lobe. These subnetworks also often cut implausibly across disjoint anatomical structures such as putamen and thalamus. These subnetworks were tentatively classified as noise-related and were ignored for future analyses.

Determining Amount of Data Required to Reliably Identify Individual-Specific Subnetworks

To determine the quantity of data required to reliably identify each subject’s frontostriatal subnetworks, we iteratively subsampled varying amounts (5–100 min) of contiguous data from each subject’s concatenated cross-session fMRI timecourses. Given the computational demands of this analysis, only 10 iterations were conducted for each data amount within each subject. For each iterative subsampling, we replicated the community detection approach employed for the whole dataset. The communities derived from the subsampled data were matched to that subject’s full-data frontostriatal subnetworks based on the Dice coefficient of overlap. We allowed a combination of multiple communities to match the subnetwork if that combination produced the optimal Dice overlap. The overall overlap for each iteration was computed as the average Dice coefficient across all frontostriatal subnetworks.

Determining Individual-Specific Large-Scale Network Identities

We identified the set of canonical large-scale networks using the individual-specific network matching approach described in (Gordon, Laumann, Gilmore, et al. 2017c). Briefly, the Infomap algorithm was applied to each subject’s correlation matrix thresholded at a range of edge densities spanning from 0.01% to 5%. At each threshold, the algorithm returned community identities for each vertex and voxel. Communities were labeled by matching them at each threshold to a set of independent group average networks described in (Gordon, Laumann, Gilmore, et al. 2017c). The matching approach proceeded as follows: 1) At each density threshold, all identified communities were compared with the independent group networks using the Jaccard Index of spatial overlap. 2) The community with the best match (highest overlap) to one of the independent networks was assigned that network identity, and then not considered for further comparison with other independent networks within that threshold. Matches lower than Jaccard = 0.1 were not considered (to avoid matching based on only a few vertices). Matches were first made with the large, well-known networks (default, lateral visual, motor hand, motor mouth, cingulo-opercular, frontoparietal, and dorsal attention), and then to the smaller, less well-known networks (ventral attention, salience, parietal memory, contextual association, medial visual, and motor foot). 3) In each individual and in the average, a “consensus” network assignment was derived by collapsing assignments across thresholds, giving each node the assignment it had at the sparsest possible threshold at which it was successfully assigned to one of the known group networks.

Computing Lobe-Wise Striatal Subnetwork Representation

For each individual subject, we defined the frontal, insular, parietal, temporal, and occipital lobes based on the lobe identities of the regions of the subject-specific FreeSurfer-generated Destrieux atlas parcellation, which was deformed into fs_LR_32k space to match the functional data. Lobe-wise representation was computed within each subject as the percent of cortical vertices among all nonnoise subnetworks that were within each lobe.

Visualizing Subnetwork Overlap across Subjects

For each matched nonnoise subnetwork, the overlap across individuals was visualized using Connectome Workbench. For cortex, the number of individuals with the subnetwork present was calculated in each vertex. For striatum, this was less straightforward, as the subcortical structures were individualized, and so were different shapes and sizes across people. To represent overlap in striatum, we first mapped the subcortical voxels of each individual among MSC02–10 to the subcortical structures in MSC01. This was done by assigning each voxel in the MSC01 structures the network ID of the nearest (by Euclidean distance) subcortical voxel in that individual. Overlap was then calculated across MSC01-space subject subnetworks as the number of individuals with the subnetwork present in each MSC01-space voxel. It should be noted that this procedure was used only for visualization of cross-subject overlap.

Striatal Subnetworks Are Individually Specific but Spatially Similar across Subjects

We delineated individual-specific subnetworks connecting striatum to cortex and matched them across subjects following procedures described in (Gordon et al. 2020). Ten separate subnetworks with both striatal and cortical representation were identified that were spatially and topologically consistent across subjects, with some interindividual variability in the localization of each subnetwork. The spatial representations of these subnetworks spanned from posterior putamen to nucleus accumbens in striatum, and from central sulcus and parietal/temporal lobes in cortex to subgenual cingulate/orbitofrontal cortex. See Figure 1 for an example subject; see Supplementary Figure 2 for all subjects; see Figure 2 for the overlap of each subnetwork across subjects; and see Supplementary Table 1 for information about the consistency of connections. These 10 subnetworks existed as substructures within several known large-scale brain networks, including the default mode, salience, cingulo-opercular, frontoparietal, language, and somatomotor networks (see Fig. 3 for an example subject; see Supplementary Table 1 for mode network memberships across subjects).

Corticostriatal subnetworks could not be reliably detected if low quantities of data were used for subnetwork identification (Dice overlap < 0.5). The reliability of identified subnetworks increased in all subjects when increasing data quantities were analyzed and reached an average Dice overlap of 0.75 when at least 90 min of data were utilized (Supplementary Fig. 3). Reliability was lower than for cortex-only functional network for similar quantities of data (Gordon, Laumann, Gilmore, et al. 2017c) likely due to the lower BOLD signal SNR in subcortical structures.

Across subjects, an average of 95% (range: 91–99%) of striatal voxels were matched to one of the ten subnetworks. Thus, the subnetworks described here represent the large majority of strong functional corticostriatal connections, with very little of the striatum containing idiosyncratic subnetwork connections that cannot be matched across subjects.

These subnetworks showed significant overlap with the nonhuman primate-based model of corticostriatal connections advanced by Haber (2016). First, these subnetworks, which represent very strong corticostriatal functional connections, most commonly connected striatum to frontal cortex. In each subject, we calculated the percent of cortical vertices with any subnetwork representation that were within each cortical lobe (Fig. 4A). Across subjects, we found that 70% ± 6% (mean ± SD) of these subnetworks’ cortical vertices were within the frontal lobe, which was higher in every subject than in the insula (10% ± 2%), parietal lobe (9% ± 3%), temporal lobe (10% ± 6%), or occipital lobe (1% ± 1%); all paired ts(9) > 17.0, all Ps < 10⁻⁷. Similarly, we evaluated the strength of within-subnetwork corticostriatal RSFC for each lobe in each subject (Fig. 4B). Across subjects, within-subnetwork RSFC between striatal and frontal regions was 0.25 ± 0.02 (mean ± SD), which was higher in every subject than RSFC between striatum and insula (0.15 ± 0.03), parietal lobe (0.11 ± 0.04), temporal lobe (0.10 ± 0.02), or occipital lobe (0.12 ± 0.06); all paired ts(9) > 6.8, all Ps < 0.0002.

The spatial distributions of these subnetworks in large parts converged with the Haber (2016) model (Fig. 5), with several important exceptions (see below). Projections to posterior and middle putamen originated in motor and premotor cortex, respectively; projections to dorsal caudate originated from dorsomedial and dorsolateral prefrontal cortex; and projections to ventral anterior caudate and medial putamen originated from the dorsal anterior cingulate cortex. Projections to nucleus accumbens originated from medial and orbitofrontal cortex; though notably, the targets of the medial and orbitofrontal cortex projections within nucleus accumbens (to the shell and core, respectively) appeared reversed relative to the nonhuman primate model.

Importantly, we discovered two subnetworks missing from the nonhuman primate-derived model of corticostriatal connections advanced by Haber (2016) (Fig. 5): 1) a subnetwork (Figs 1–3, #3, Med Ant Putamen, pink) with representation primarily in medial putamen and sometimes lateral caudate head, as well as in a distributed set of cortical regions including dorsomedial prefrontal cortex, inferior frontal gyrus, posterior middle frontal gyrus, and superior temporal sulcus; and 2) a subnetwork (Figs 1–3, #6, Ant Caudate, dark green) with representation in anterior lateral caudate and distributed cortical regions immediately adjacent to the medial putamen subnetwork. See Figure 6A for these subnetworks in all subjects. Notably, the cortical distributions of the medial putamen subnetwork converges with the known distribution of the human language network (Braga et al. 2020), while the distribution of the lateral anterior caudate subnetwork converges with a previously described default mode subnetwork known to link core default regions to the language network (Gordon et al. 2020).

Like the human language network, these two subnetworks were left-lateralized, with more extensive left than right hemisphere representation in both cortical and subcortical structures. To quantify this observation, for each subnetwork, we calculated a laterality index as ([# left hemisphere vertices/voxels] − [# right hemisphere vertices/voxels])/[# total vertices/voxels], separately in cortex and in striatum. With this index, positive values indicate left lateralization, while negative values indicate right lateralization. In every subject, these two subnetworks were left-lateralized in both cortex and striatum (medial putamen subnetwork: cortex mean laterality = 0.62 ± 0.24, one-sample t(9) = 7.9, P < 10⁻⁴; striatum mean laterality = 0.73 ± 0.32, t(9) = 7.2, P < 10⁻⁴; lateral caudate subnetwork: cortex mean laterality = 0.67 ± 0.34, one-sample t(9) = 5.7, P < 0.001; striatum mean laterality = 0.72 ± 0.36, t(9) = 5.4, P < 0.002); see Figure 6B. By contrast, no other subnetwork exhibited consistent lateralization across subjects (all ts(9) < 2.2, all Ps > 0.05 in both cortex and striatum).

Most Corticostriatal Subnetworks Converge with Primate Tract-Tracing

We identified striatal subnetworks composed of only very strong functional connections in single individuals, and we found that 1) the majority of these very strong corticostriatal connections were from frontal cortex, and 2) frontal regions of these subnetworks had stronger connectivity with striatum than nonfrontal regions, suggesting that the nonfrontal regions we did observe were less strongly associated with striatal activity. This converges with nonhuman primate work arguing that, while striatum does receive inputs from diverse cortical sources across frontal, parietal, insular, and temporal cortex (Yeterian and Van Hoesen 1978; Van Hoesen et al. 1981; Selemon and Goldman-Rakic 1985, 1988; Cavada and Goldman-Rakic 1991; Steele and Weller 1993; Webster et al. 1993; Yeterian and Pandya 1993, 1998; Cheng et al. 1997; Chikama et al. 1997), the strongest, driving projections to striatum arise from frontal lobe (Haber 2003). Thus, the present subnetwork connections represent a substantial increase in specificity over prior human neuroimaging work that identified corticostriatal connections originating from large swaths of cortex and from distributed networks across multiple cortical lobes (Choi et al. 2012; Greene et al. 2014; Tziortzi et al. 2014; Morris et al. 2016; Parkes et al. 2017).

The identified frontostriatal circuits were highly ordered as a topographically mirrored rostral–caudal progression of corticostriatal projections. This organization converged with the organization of these circuits proposed based on findings from nonhuman primates (Haber 2003, 2016); see Figure 5. Specifically, moving from caudal to rostral striatum:

1. We found functional connections between somatomotor regions in the central sulcus and the posterior putamen (Fig. 5, magenta), consistent with previous observations of projections from precentral motor cortex to posterior putamen in humans (Zeharia et al. 2015) and nonhuman primates (Künzle 1975; Flaherty and Graybiel 1994). These connections are likely related to striatal control over motor function.

2. We found functional connections between regions just anterior to the somatomotor strip, particularly in dorsomedial prefrontal cortex, and the middle lateral putamen (Fig. 5, green), consistent with previous observations in nonhuman primates of projections from pre-/supplementary motor cortex to middle lateral putamen (Künzle 1978; Haber 2003). These circuits are engaged during performance of complex (nonautomatized) motor skills (Boecker et al. 1998).

3. We found functional connections between dorsomedial and dorsolateral prefrontal cortex, within the frontoparietal network, and a large extent of dorsal caudate (Fig. 5, orange). This observation is consistent with previous accounts of projections from dorsolateral prefrontal cortex to terminals in dorsal caudate that span a significant rostral–caudal extent (Selemon and Goldman-Rakic 1985; Arikuni and Kubota 1986; Haber et al. 2006), which have been argued to play a critical role in working memory (Partiot et al. 1996).

4. We found several sets of functional connections between dorsal anterior cingulate cortex/anterior insula and the anterior caudate/putamen (Fig. 5, yellow/maroon/brown). These observations align with known projections from primate dorsal anterior cingulate to anterior caudate and putamen, and particularly to sites more dorsal and lateral than the orbitofrontal cortex (OFC) or ventromedial prefrontal cortex (vmPFC) projections described below (Haber et al. 2006). This set of cortical regions (dorsal anterior cingulate, anterior insula, and anterior striatum) is considered the core of the adjacent “cingulo-opercular” and “salience” networks, which are critical for a wide variety of cognitive control operations.

5. We found functional connections between nucleus accumbens and pregenual medial prefrontal cortex (mPFC), and a separate set with OFC/subgenual cingulate cortex (Fig. 5, blue/red). These connections match circuitry known to be critical for reward processing/addiction (Haber and Knutson 2010; Koob and Volkow 2016) and mood disorders (Drevets et al. 2008), respectively. The nucleus accumbens is anatomically segregated into a more dorsal shell and a more ventral core (Salgado and Kaplitt 2015). Here, the nucleus accumbens nodes of these two subnetworks exhibited a similar dorsal/ventral distribution; as such, we tentatively label these nodes the accumbens shell and core, respectively.

Interestingly, these connections only partially converged with connectivity identified in nonhuman primates, in which the nucleus accumbens core preferentially receives projections from pregenual mPFC, while the shell receives projections from OFC—the reverse of the arrangement we observed here. It is possible that this divergence from nonhuman primate work is driven by fMRI data quality issues, such as the known presence of susceptibility artifact within ventral prefrontal and orbitofrontal cortex in human fMRI. However, human diffusion tractography imaging has also demonstrated that the accumbens core is more strongly connected to OFC than the shell (Baliki et al. 2013; Xia et al. 2017; Cartmell et al. 2019), convergent with our present findings. Hence, this divergence may represent a true organizational difference between humans and nonhuman primates related to reward processing.

Several of the connections identified here, including to the posterior and middle putamen, dorsal caudate, and nucleus accumbens core, were restricted to a single target subcortical structure. This represents an increase in target specificity relative to the (Haber 2016) model, in which similar frontostriatal projections span the internal capsule (Fig. 5A,B, green and orange) or extend out of nucleus accumbens into caudate (Fig. 5A,B, red). We consider it most likely that frontal projections to multiple striatal structures do exist in humans, but that projections to secondary structures do not drive striatal function as strongly as projections to the primary target structure.

Interestingly, while every subnetwork had representation in the medial frontal cortex in every subject, only a few (e.g., dorsal caudate and posterior putamen) were consistently represented on the lateral surface of the frontal cortex. The reasons for these apparent preferential connections between medial frontal cortex and striatum are not immediately clear. However, we hypothesize that they may be driven by the tendency of evolution to optimize connections between brain regions while minimizing “wiring cost”—particularly the volume of interregional axons—within the constraints of the volume of the skull (Klyachko and Stevens 2003; Bullmore and Sporns 2012). Notably, the medial cortical elements of the subnetworks appear to be arranged so as to globally minimize distances to their striatal targets (see Fig. 1). Further work is needed to test this possibility.

Language-Related Corticostriatal Subnetworks Specific to Humans

We observed two subnetworks that did not converge with proposed models of nonhuman primate frontostriatal connectivity (Haber 2016). These consisted of two adjacent sets of distributed frontotemporal regions connected to anterior medial putamen and anterior lateral caudate (Fig. 5B, pink/dark green). Based on several lines of evidence, we argue that these connections represent circuitry for language processing. First, the cortical elements of these connections are highly convergent with previous descriptions of language-related networks in group-average (Friederici and Gierhan 2013; Ji et al. 2019) and individual (Braga et al. 2020; Gordon et al. 2020) human neuroimaging data. Specifically, cortical aspects of the anterior medial putamen (pink) subnetwork converge with the “core” language network (Braga et al. 2020), while those of the anterior lateral caudate (dark green) subnetwork correspond to a subnetwork we have previously described as linking the language network to the default mode network (Gordon et al. 2020). Second, as with human language function (Friederici 2011; Corballis 2012; Ocklenburg et al. 2014), connections of these two subnetworks were consistently lateralized, with greater extent in the left hemisphere than in the right hemisphere in both cortex and striatum. No other corticostriatal network structure exhibited consistent lateralization across subjects. Third, the striatal elements of these connections converge well with previous work reporting activation in anterior medial putamen (Friederici et al. 2003) and anterior lateral caudate (Mestres-Missé et al. 2012) in response to specifically effortful, nonautomatic language processing. These findings suggest that the striatal elements of these language subnetworks may be critical for learning and executing new language-related tasks, similar to the learning and optimization role ascribed to the striatum in other behavioral domains (Cohen and Frank 2009).

Interestingly, this network is not wholly novel in humans; some connections of this human corticostriatal language network do exist in nonhuman primates. Lateral prefrontal cortex and superior temporal sulcus have been shown to project to each other (Yeterian et al. 2012), and superior temporal cortex projects to lateral caudate and medial putamen (Selemon and Goldman-Rakic 1985; Yeterian and Pandya 1998). However, the multiple elements of the language network (inferior frontal gyrus, posterior middle frontal gyrus, medial superior frontal gyrus, superior temporal sulcus, and striatum) are not known to be all interconnected in nonhuman primates, as our work suggests here in humans. This supports the idea that the human language network is an evolutionarily novel network scaffolded in part upon existing circuitry present in prehuman primates (Dick and Tremblay 2012; Friederici 2017).

Toward a Complete Description of Integrated Cortico-Striato-Thalamic Circuits

It is well established in rodents and nonhuman primates that corticostriatal circuits are further projected to the thalamus, and then back to cortex, forming feedback loops (Nakano et al. 2000; Haber 2003). Thus, in principle, the corticostriatal subnetwork structures we observe here should also be represented in the thalamus. While we did observe this representation to some degree, the thalamus did not consistently contain representations of all corticostriatal subnetworks, and further, the thalamus tended to have extensive representation of subnetworks characterized as noise-related. This suggests that the quality of the MSC fMRI data, while adequate for winner-take-all style whole-network mapping in thalamus (Greene et al. 2020), may not be adequate for the data-driven mapping of highly specific cortico-striato-thalamo-cortical circuits in thalamus. Future work may be able to more broadly conduct such mapping using higher-quality data collected with advanced multiband and multiecho sequences (Lynch et al. 2020).

While the presented delineations of corticostriatal circuits are neurobiologically compelling and individual-specific, they are incomplete. The present analyses focus on the very strongest functional connections to each striatal voxel, representing the likely driving inputs to the striatum. However, frontostriatal circuits are only partially segregated and do exhibit some degree of cross-circuit connectivity (Haber 2003). At the large-scale network level, human RSFC work has also demonstrated cross-network integration zones (Greene et al. 2020). We reproduced analyses from (Greene et al. 2020) at the subnetwork level rather than the network level, and we similarly observed cross-subnetwork connections (Supplementary Fig. 4), suggesting that weaker connections may enable interactions between frontostriatal loops above and beyond the specific, driving inputs that are the focus of the current paper. Untangling the complexity of these multisubnetwork connections for comprehensive evaluation and categorization will be conducted in future work.

Functional connectivity techniques cannot distinguish between mono- and polysynaptic connections, so it is unclear whether the anterior caudate and putamen actually receive direct projections from the language network. However, we note that when only the very strongest functional connections are considered, resulting subnetworks recapitulate known driving monosynaptic frontal cortex inputs to the striatum fairly closely (Haber 2003, 2016). This suggests that while functional connectivity in general may reflect a combination of mono- and polysynaptic projections, the strongest corticostriatal functional connections are more likely to primarily represent monosynaptic projections.

The Utility of Mapping Corticostriatal Circuits in Individual Humans

The present results have important implications for the future study of the human striatum. The discrete subnetworks present within the striatum likely enable a variety of different cognitive and motor processes that can be localized to each separate, individual-specific striatal subnetwork. With such localization, we can generate stronger hypotheses about how individual-specific symptoms observed in neuropsychiatric disorders, neurodegenerative disease, or brain injuries might arise from disruption of those specific frontostriatal circuits. Such symptoms might be optimally alleviated restored via brain stimulation approaches targeting those individual-specific circuits, including noninvasive TMS applied to the cortical sources of frontostriatal projections, or deep electrode implantation localized to their striatal targets.

Notably, the present results illustrate why it is critical to localize functional networks in repeatedly sampled individual humans, rather than in group-average data. The frontal–striatal subnetworks exhibited significant interindividual spatial variability, as illustrated by the limited cross-subject overlap observed in Figure 2. This result is consistent with prior work illustrating that poor cross-subject overlap is more apparent when comparing small network structures than larger networks (Gordon et al. 2020; Harrison et al. 2020). This variability highlights the need for a high degree of precision when attempting to target specific circuits for intervention (e.g., with TMS or DBS).