Striatal subdivisions that coherently interact with multiple cerebrocortical networks

Abstract The striatum constitutes the cortical‐basal ganglia loop and receives input from the cerebral cortex. Previous MRI studies have parcellated the human striatum using clustering analyses of structural/functional connectivity with the cerebral cortex. However, it is currently unclear how the striatal regions functionally interact with the cerebral cortex to organize cortical functions in the temporal domain. In the present human functional MRI study, the striatum was parcellated using boundary mapping analyses to reveal the fine architecture of the striatum by focusing on local gradient of functional connectivity. Boundary mapping analyses revealed approximately 100 subdivisions of the striatum. Many of the striatal subdivisions were functionally connected with specific combinations of cerebrocortical functional networks, such as somato‐motor (SM) and ventral attention (VA) networks. Time‐resolved functional connectivity analyses further revealed coherent interactions of multiple connectivities between each striatal subdivision and the cerebrocortical networks (i.e., a striatal subdivision‐SM connectivity and the same striatal subdivision‐VA connectivity). These results suggest that the striatum contains a large number of subdivisions that mediate functional coupling between specific combinations of cerebrocortical networks.


| INTRODUCTION
The cerebral cortex and the striatum constitute a critical part of the cortical-basal ganglia loop, in which information is processed through the cerebral cortex, basal ganglia, and thalamus and returns to the cerebral cortex (Alexander, DeLong, & Strick, 1986;Haber, 2003). The striatum has classically been discriminated into the caudate nucleus, which is connected with the prefrontal cortex, and the putamen, which is connected with the frontal motor cortex, and it is thought to select the most appropriate behavior by interacting with the cerebral cortex (Alexander et al., 1986;Haber, 2003;Middleton & Strick, 2000). Using resting-state functional and structural connectivity (Behrens et al., 2003;Fox & Raichle, 2007;van den Heuvel & Pol, 2010), the precise architecture of the human striatum and its converging input from the cerebral cortex have been revealed by using clustering analyses based on the global similarity of functional/structural connectivity (Barnes et al., 2010;Choi, Yeo, & Buckner, 2012;Di Martino et al., 2008;Draganski et al., 2008;Garcia-Garcia et al., 2018;Janssen, Jylänki, Kessels, & van Gerven, 2015;Jarbo & Verstynen, 2015;Jaspers, Balsters, Kassraian Fard, Mantini, & Wenderoth, 2017;Jung et al., 2014;Verstynen, Badre, Jarbo, & Schneider, 2012). Recent animal tracer studies have also revealed fine organization in the striatum, particularly in its rostral part, in relation to the converging input from the cerebral cortical areas (Averbeck, Lehman, Jacobson, & Haber, 2014;Choi, Ding, & Haber, 2017;Choi, Tanimura, Vage, Yates, & Haber, 2017;Haber, Kim, Mailly, & Calzavara, 2006). However, it remains unclear how the striatum functionally interacts with Akitoshi Ogawa and Takahiro Osada contributed equally to this work. multiple cerebrocortical networks in the temporal domain. The central question of this study is to reveal the temporal aspects of the multiple striatal-cerebrocortical interactions; for this purpose, it is important to analyze such interactions on the basis of the functionally distinct subdivisions.
There are two major approaches for areal parcellation (Schaefer et al., 2017): Clustering analyses reveal cortical areas that represent the global cortical functional architecture well. In contrast, boundary mapping analyses detect boundaries based on local gradients of connectivity profiles. It has been suggested that the local gradient approach is more suitable for delineating cortical areas because detecting abrupt changes in connectivity profiles is similar to the histological delineation of cortical areas (Barnes et al., 2010;Biswal et al., 2010;Cohen et al., 2008;Glasser et al., 2016;Gordon et al., 2016;Hirose et al., 2012Hirose et al., , 2013Hirose et al., , 2016Laumann et al., 2015;Margulies et al., 2007;Nelson et al., 2010;Poldrack et al., 2015;Schaefer et al., 2017;Xu et al., 2016). Therefore, it is expected that boundary mapping analyses are more suitable for the stable detection of smaller striatal subdivisions and that boundary mapping analyses can reveal the temporal aspects of the striatal-cortical interactions more accurately.
In the present functional MRI study, boundary mapping analyses were applied to the striatum to obtain the finer functional architecture of the human striatum. We delineated approximately 100 striatal subdivisions reproducibly across independent data sets. Based on previous studies of converging input from the cortex, it was expected that each striatal subdivision was connected with multiple cerebrocortical networks. We next examined how the striatum interacted with the cerebrocortical networks by analyzing connectivity changes over time (Hutchison et al., 2013;Shine et al., 2015Shine et al., , 2016Zalesky, Fornito, Cocchi, Gollo, & Breakspear, 2014) to reveal whether the multiple interactions between a striatal subdivision and multiple cortical networks are coherent or independent.

| MRI procedures
All MRI data were acquired using a 3-T MRI scanner (Siemens Skyra, Erlangen, Germany). T1-weighted structural images were obtained for anatomical reference (resolution = 0.8 × 0.8 × 0.8 mm 3 ). Functional images were obtained using multiband gradient-echo echo-planar sequences (Feinberg et al., 2010;Moeller et al., 2010;Xu et al., 2013; TR = 4.0 s, TE = 41.6 ms, flip angle = 73º, FOV = 160 × 160 mm 2 , matrix size = 128 × 128, 120 contiguous slices, voxel size = 1.25 × 1.25 × 1.25 mm 3 , multiband factor = 4). We acquired 100 volumes in each fMRI run at the resting state and repeated the process for 10 runs in each of the 10 daily sessions. Thus, 10,000 total volumes were collected for each participant. Eight of the 10 participants also participated in a task fMRI experiment and were scanned for six runs using the same scanning parameters.
To increase the multiband factor, a small FOV (160 × 160 mm 2 ) was set. This small FOV did not always cover the whole brain along the anterior-posterior axis. However, the aliasing artifacts in the frontal lobe (from the occipital lobe) were minimal, and the aliasing overlap artifact observed in the frontal lobe was 0.1%-4.9% of the whole cerebral cortex, keeping the striatum intact. It has also been demonstrated that parcellation procedures do not require the entire cerebral cortex for the appropriate detection of the connectivity transitions (Hirose et al., 2012(Hirose et al., , 2013. Therefore, it is unlikely that the FOV influenced either the signal acquisition in the striatum or the data analyses used in the present study.
Parcellation analyses based on boundary mapping (Cohen et al., 2008;Glasser et al., 2016;Gordon et al., 2016;Hirose et al., 2012Hirose et al., , 2013Hirose et al., , 2016Laumann et al., 2015;Osada et al., 2017) were applied to the striatum (Figure 1). Each voxel in the striatum of each participant was used as a seed to calculate its correlations with the voxels in the gray matter of the cerebral cortex. One part of the cerebral cortex that suffered from overlap due to aliasing artifacts was removed from the target voxels. Voxel-wise correlation coefficients in the cerebral cortex were converted to Fisher's z (Fisher z transformation).
The similarity of the spatial patterns of the correlation maps was then evaluated using correlation coefficients, and similarity maps were generated. After minimal spatial smoothing (full width at halfmaximum [FWHM] = 1.25 mm), spatial gradients of the similarity maps were computed for each seed voxel. A three-dimensional watershed algorithm (Vincent & Soille, 1991) was applied to the gradient maps, and the binary watershed maps were averaged across seed voxels after spatial smoothing (FWHM = 1.25 mm) to generate a boundary probability map. The watershed algorithm was again applied to the boundary probability map for each participant. A cluster with three or more contiguous voxels was defined as a striatal subdivision, and the voxel with the smallest boundary probability in the striatal subdivision was defined as a striatal subdivision center (SSC).
The striatum was manually segmented using the functional image of each participant. The border of the striatum was clear for the most part, but the ventral surface of the ventral striatum that was relatively obscure was determined using the T1-weighted structural images as reference. It should be noted that the boundary mapping of the ventral striatum is relatively less sensitive for detecting striatal subdivisions. The striatum analyzed in the present study includes both the dorsal and ventral striatum. However, the ventral striatum is much smaller in size than is the dorsal striatum, and the outermost voxels of the ventral striatum, which are not detected as centers of the subdivisions using boundary mapping analyses, occupy a large portion of the ventral striatum.

| Parcellation of the cerebral cortex
Functional images were preprocessed for resting-state functional connectivity in the same way as the striatum (see section 2.3). For volume to surface mapping, the middle of the gray matter was identified within a functional image for each participant using Caret software (Marcus et al., 2011;Van Essen et al., 2001), and a fiducial surface image was generated from the middle of the gray matter.
The cerebrocortical parcels were calculated in a similar manner, as described previously (Cohen et al., 2008;Gordon et al., 2016;Laumann et al., 2015; Supporting Information Figure S1). Each vertex in the fiducial surface in the cerebral cortex of each participant was used as a seed to calculate its correlations with all of the vertices. The similarity of the spatial patterns of the correlation maps was then evaluated using correlation coefficients, and similarity maps were generated. After spatial smoothing (FWHM = 6.0 mm) (Gordon et al., 2016;Laumann et al., 2015), spatial gradients of the similarity maps were computed for each seed vertex. A two-dimensional watershed algorithm was applied to the gradient maps, and the binary watershed maps were averaged across the seed vertices after spatial smoothing (FWHM = 6.0 mm) to generate a boundary probability map. The watershed algorithm was again applied to the boundary probability map to delineate cerebrocortical parcels for each participant. The voxel with the smallest boundary probability in the cerebrocortical parcel was defined as a cerebrocortical parcel center.

| Quality control of head movement
We evaluated the amount of head motion by using frame-wise displacement (FD) (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012).
FD is a measurement of instantaneous head motion that can be calculated as a locational difference between two successive volumes and is an important measure for the quality control of resting-state data (Gordon et al., 2016;Laumann et al., 2015). Frames with FD > 0.2 mm were censored, as well as uncensored segments of data lasting fewer than five contiguous volumes; all such data were excluded from the subsequent analysis. The included images were 75.1 AE 12.6% (mean AE SD) of the total acquired images, and the resultant FD was 0.084 AE 0.012 mm (mean AE SD).

| Evaluation of the relationship between the striatal subdivisions and cerebrocortical parcels
To evaluate whether the striatal subdivisions and cerebrocortical parcels were related to each other, the voxel-wise variation of connectivity with the cerebrocortical parcels was examined within each striatal subdivision. For each center voxel of the cerebrocortical parcels, Fisher's z value was calculated with all of the voxels in the striatum; then the SD of the z values was calculated within each striatal subdivision.  Overview of boundary mapping analyses of the striatum. functional images were processed, and the correlation, similarity, gradient, boundary probability, and striatal parcel center maps were generated as shown for the analysis flow. Representative examples for these maps are shown in the right plots 2.7 | Replication of larger-scale functional organization using public atlases We replicated previous observations on the larger-scale architecture of the striatal connectivity with cortical gyri. Individual cortical surfaces were converted to the cerebrocortical surface atlas (Yeo et al., 2011) by Freesurfer (Fischl, Sereno, Tootell, & Dale, 1999), and multimodal surface matching (MSM; Robinson et al., 2014) was then applied. The correlations of each SSC with all of the voxels in the precentral gyrus (PCG), middle frontal gyrus (MFG), orbitofrontal gyrus (OFG), and inferior parietal lobule (IPL) (Desikan et al., 2006) were calculated, and the Fisher's z values across all of the voxels within each gyrus were averaged. The gyrus that was connected most strongly with each of the SSCs was determined out of the four gyri, and four SSC binary maps were generated that showed the SSCs that were most strongly connected with each gyrus. The binary SSC maps were averaged across participants for each gyrus, with spatial smoothing by a 4-mm FWHM kernel. The 4-mm kernel was chosen to investigate the striatum on a larger scale, but the kernel was kept smaller than the 6-mm kernel applied to the cerebral cortex. The resultant values for the smoothed images were normalized such that the total values across the entire striatum were 1.
Similarly, correlations of each SSC with the seven networks including the somato-motor (SM), ventral attention (VA), limbic (Lim), and fronto-parietal (FP) (Yeo et al., 2011) were calculated. The network that was connected most strongly with each of the SSCs was determined out of the seven networks, and seven SSC binary maps were generated to show the SSCs that were most strongly connected with each network. The binary SSC maps were averaged across participants for each network, with spatial smoothing by a 4-mm FWHM kernel.

| Dynamic functional connectivity analysis
We determined the two of the seven cerebrocortical networks that were connected most strongly with each SSC. Most of the striatal subdivisions were connected with the specific combinations of cerebrocortical networks, primarily SM-VA and VA-FP combinations. To investigate whether the multiple interactions were temporally coherent, we used a dynamic functional connectivity analysis (Hutchison et al., 2013) using multiplication of temporal derivatives (MTD), a sensitive index of functional connectivity within a small time-window (Shine et al., 2015(Shine et al., , 2016. The MTD allows greater temporal resolution of time-resolved connectivity in BOLD time series data than does the conventional sliding-window Pearson's correlation coefficient (Shine et al., 2015).
Functional images were preprocessed for dynamic functional connectivity analysis (Leonardi & Van De Ville, 2015;Shine et al., 2015Shine et al., , 2016. Images were corrected for slice timing and realigned using SPM8 (www.fil.ion.ucl.ac.uk/spm/). Temporal filters (0.05 Hz < f < 0.125 Hz) were applied to the functional images. A general linear model was used to regress out nuisance signals that correlated with head motion, whole-brain global signals, averaged ventricular signals, and averaged white matter signals. Time series data were extracted from each voxel in the SM, VA, and FP networks (Yeo et al., 2011) and from each SSC that was connected most strongly with the SM-VA or VA-FP network combination. The MTD between the striatal subdivisions and the cerebrocortical networks was calculated based on the sliding window of 7 TRs, with a step size of 1 TR. The MTD was then averaged across the voxels in each cerebrocortical network and further across the voxels of the SSCs that were connected most strongly with the SM-VA or VA-FP network.

| Striatal subdivisions
The parcellation analysis (Figure 1) (Figure 2b). In a group analysis with fine normalization using DARTEL (Ashburner, 2007), the boundary pattern in the striatum was almost lost (Supporting Information Figure S2b). The boundary probability map was highly reproducible between independent data sets, when the total data sets were divided into odd or even runs (Figure 2c Table 1 for all participants). When the high spatial frequency components (≤2.5 mm) of the probability map that constituted the boundary pattern were extracted using Fourier transformation, the average correlation coefficient of the high-frequency component was .37 (Fisher's z = 0.39) [t(9) = 7.2, p < .001 after z transformation]. Figure 2d demonstrates spatial profiles of the probabilistic boundary patterns that exhibited a close match across the two independent data sets. As shown in the distribution of the distance between adjacent striatal subdivisions in Supporting Information Figure S3, the size of the striatal subdivisions was reproducible across participants. The distribution of the size was also similar between the left and right striatum. In almost all of the participants, the most common distance between adjacent subdivisions was 2 mm.

| Relationship of the striatal subdivisions with cerebrocortical parcels
We examined how the striatum was related to the cerebrocortical parcels at the subdivision level by calculating the across-voxel SD of

| Intricate organization of the striatal subdivisions connected with the cerebral cortex
We next examined the functional connectivity of striatal subdivisions with the cerebral cortex. The SSC was defined as the voxel with the smallest boundary probability value in a subdivision, and the functional connectivity was calculated with the SSC voxels as seeds. The functional connectivity with the cerebral cortex was significantly reproducible between the two independent data sets, when the overall data sets were divided into odd or even runs (Figure 4a), with an average correlation coefficient of .46 (Fisher's z = 0.50) [t(9) = 21.0, p < .001 after z transformation] (Table 1). Figure 4b shows a correlation matrix of SSCs in the odd/even runs in the same participant, demonstrating high reproducibility between cerebrocortical correlation maps from the same subdivisions along the diagonal line.
The similarity of cerebrocortical correlation maps was then examined from adjacent (2.36 voxels or 2.95 mm apart, on average) striatal   (Table 1).

| Larger-scale functional organization of the striatum
Previous studies of striatal functional architecture have revealed a well-organized relationship for cerebrocortical-striatal connectivity (Barnes et al., 2010;Choi et al., 2012;Di Martino et al., 2008;Draganski et al., 2008;Garcia-Garcia et al., 2018;Janssen et al., 2015;Jaspers et al., 2017;Jung et al., 2014). We confirmed, as a positive control, that the larger-scale observations could be replicated using the present data set. Supporting Information Figure S4a shows group average maps of the SSCs that were most strongly connected with the PCG, MFG, OFG, or inferior parietal lobule, based on the atlas by Desikan et al. (2006). Consistent with previous studies (Barnes et al., 2010;Choi et al., 2012;Di Martino et al., 2008;Draganski et al., 2008;Garcia-Garcia et al., 2018;Janssen et al., 2015;Jaspers et al., 2017;Jung et al., 2014), the middle part of the putamen was most strongly connected with the PCG, the caudate nucleus was most strongly connected with the MFG, and the caudate nucleus extending to ventral striatum was most strongly connected with the OFG.
Similarly, Supporting Information Figure S4b shows group average maps of SSCs that were connected most strongly with one of the seven cerebrocortical networks, including the SM, VA, Lim, and FP networks, as reported previously (Yeo et al., 2011). Consistent with a previous study (Choi et al., 2012), the ventral part of the putamen was most strongly connected with the SM network, the dorsal part of the putamen was connected most strongly with the VA network, the caudate nucleus was most strongly connected with the FP network, and the caudate nucleus extending to the ventral striatum was most strongly connected with the Lim network. These results confirmed the larger-scale functional organization of the striatum.

| Striatal subdivisions connected with multiple cerebrocortical networks
We then examined how the striatal subdivisions were convergently connected with the two cerebrocortical networks. When the second strongest connectivity was less than 30% of the strongest one, the SSC was excluded from the analysis, resulting in exclusion of 21.2% of the total SSCs. Figure 5a shows the combinations of the two cerebrocortical networks with which the SSCs were connected most strongly.
The most common combination (12.2%) was the SM-VA networks, and the distribution of the combinations was highly biased toward specific combinations, such as the SM-VA and VA-FP networks [χ 2 test, χ 2 (20) = 35.2, p < .019].
It has been demonstrated that the spatial extent of the DM network is variable among individuals (Braga & Buckner, 2017). To test whether the individual variability affected the distribution of the combinations of the connected networks, the regions that were likely to  (Figure 5b), which indicates that the biased distribution was not due to individual differences in the spatial extent of the networks.
To further examine the distribution of the combinations in a more conservative way, 37.4% and 62.2% of the total SSCs were excluded from analysis when the second strongest connectivity was less than 50% and 70% of the strongest one. The distribution of the combinations was almost unchanged (Supporting Information Figure S5a). The distribution was again almost unchanged when the size of the networks was reduced to approximately one half of the original size (Supporting Information Figure S5b).

| Striatal cerebrocortical interaction revealed by dynamic functional connectivity
Striatal subdivisions were likely connected with multiple cerebrocortical networks such as the SM-VA and VA-FP networks. Figure 6a (left) shows the time courses of MTD (Shine et al., 2015(Shine et al., , 2016

| Brain activation during eye movements
We also measured the brain activity when participants performed a standard eye movement task (Supporting Information Text S1 and Figure S6a).
Although a group analysis revealed activations of both the frontal eye field and the striatum, a single-level analysis revealed activations only in the  Figure S6b). The lack of activation in the striatum at this single-level analysis is consistent with a study by the Human Connectome Project Van Essen et al., 2013), which reported a similar phenomenon.

| DISCUSSION
In the present study, approximately 100 subdivisions were revealed using boundary mapping analyses applied to functional images of the striatum. The cerebrocortical parcels were more likely connected uni- It has been demonstrated that the cerebrocortical parcels show homogeneity within their regions (Gordon et al., 2016). The homogeneity is estimated by the greater proportion of the principal component of functional connectivity profiles against a null hypothesis using shifted boundaries. The homogeneity is demonstrated when the number of surface vertices within a parcel is sufficiently great (see figure 3 in Gordon et al., 2016). Thus, the homogeneity analysis was not applicable to the small structure of the striatum, and it is not certain whether the parcellated striatal areas can be considered as "striatal parcels." On the other hand, the striatal subdivisions were highly reproducible when the data sets were divided into two halves. Moreover, the boundary pattern of the striatal subdivisions produced smaller across-voxel SD of functional connectivity with the cerebrocortical parcels. Therefore, although a much higher spatial resolution would be required to demonstrate the "striatal parcels," it is possible that the boundaries of the striatal subdivisions presented in this study reflects functional differences in the areas separated by the boundary. It has been shown that striatum contains chemical compartments called "striosomes" that are approximately 0.5-0.8 mm in size in humans (Graybiel & Ragsdale, 1978).
Although it seems unlikely that the striatal subdivisions reported in the present study simply correspond to the striosomes, future investigations might reveal that multiple striosomes or a complex of a striosome and the surrounding compartment, called the "matrix," are functionally related to striatal subdivisions.
The cerebral  . To investigate individual differences, however, the individual brains have to be aligned to a common atlas. Although the alignment of the cerebral cortex was successful, even fine normalization using DARTEL was not successful in the striatal parcellation results (Supporting Information Figure S2). Thus, normalization using multimodal features such as structure and functional connectivity, as implemented in MSMAll (Glasser et al., 2016;Robinson et al., 2014) for the cerebral cortex, seems to be required for the striatum.
The analyses of the present study can be applied to various striatal-cortical interactions. For example, Zhang, Ide, and Li (2012) examined the functional connectivity of three motor areas (SMA, aPreSMA, and pPreSMA) with the subcortical structures including the striatum and revealed a differential pattern of functional connectivity within the striatum when the seed was moved along the y axis. The analyses of the present study may reveal more refined architecture of the striatal-cortical connections by examining the patterns of functional connectivity between the subdivisions and the cortical parcels in these motor areas. Such analyses may also be used to investigate whether the multiple striatal-cortical interactions from these motor areas are temporally coherent or independent.
Although the similarity of correlation maps from adjacent striatal subdivisions ( Figure 4) was barely significant in the group statistics, the average correlation coefficient was very small (Fisher's z = 0.017), suggesting that the functional organization of the striatum is rather intricate when viewed at the scale of the striatal subdivisions. At the same time, the larger scale organization of the striatum from the same data (Supporting Information Figure S4) was consistent with that found in previous studies due to larger spatial smoothing and group averaging. Animal tracer studies have demonstrated the convergence of projections from multiple cortical areas to the striatum (Averbeck et al., 2014;Choi, Ding, et al., 2017;Choi, Tanimura, et al., 2017;Haber et al., 2006). Human studies using structural and functional connectivity have also provided consistent evidence of converging projections (Barnes et al., 2010;Choi et al., 2012;Di Martino et al., 2008;Draganski et al., 2008;Janssen et al., 2015;Jarbo & Verstynen, 2015;Jaspers et al., 2017;Jung et al., 2014;Verstynen et al., 2012).
The present study further suggests that many of the striatal subdivisions are connected with multiple cerebrocortical networks that consist of multiple cortical areas. The combinations of the cerebrocortical networks connected with the striatal subdivisions were not equally distributed but were restricted to specific networks, excluding the possibility that striatal subdivisions equally cover all of the combinations of cerebrocortical networks, as "hypercolumns" do (Hubel, 1982). For example, it is possible that the sensory and motor processing of the SM network might trigger the bottom-up attentional processing of the VA network, through input from one of the networks and output to both of the networks, via the cortical basal ganglia loop.
Although the precise interpretation of functional connectivity results may require further investigation, the present study revealed approximately 100 striatal subdivisions that may mediate the functional integration of cerebrocortical networks.

ACKNOWLEDGMENTS
We thank T. Kamiya and K. Yamada for technical assistance. This work was supported by JSPS KAKENHI grants 16K16076 to A.O.,