Human Fronto-Striatal Connectivity is Organized into Discrete Functional Subnetworks

The striatum is interconnected with the cerebral cortex via multiple recurrent loops that play a major role in many neuropsychiatric conditions. Primate cortico-striatal connections can be precisely mapped using invasive tract-tracing. However, noninvasive human research has not mapped these connections with anatomical precision, limited by the practice of averaging neuroimaging data across individuals. Here we utilized highly-sampled resting-state functional connectivity MRI for individually-specific precision functional mapping of cortico-striatal connections. We identified ten discrete, individual-specific subnetworks linking cortex—predominately frontal cortex—to striatum. These subnetworks included previously unknown striatal connections to the human language network. The discrete subnetworks formed a stepped rostral-caudal gradient progressing from nucleus accumbens to posterior putamen; this organization was strongest for projections from medial frontal cortex. The stepped gradient organization fit patterns of fronto-striatal connections better than a smooth, continuous gradient. Thus, precision subnetworks identify detailed, individual-specific stepped gradients of cortico-striatal connectivity that include human-specific language networks.

1993; Pandya, 1998, 1993;Yeterian and Van Hoesen, 1978), though projections from 1 frontal cortex are thought to be the main input driving striatal function (Haber, 2003). Non-human 2 primate fronto-striatal circuits have an ordered organization that is mirrored in striatum and cerebral 3 cortex (Haber, 2016). Specifically, rostral portions of frontal cortex project to punctate regions within 4 rostral striatum (Haber et al., 1995;Kunishio and Haber, 1994), while progressively more caudal 5 frontal cortex projects to progressively caudal striatal targets (Averbeck et al., 2014;Calzavara et 6 al., 2007;Flaherty andGraybiel, 1994, p. 1975;Künzle, 1975;Selemon and Goldman-Rakic, 1985). 7 The fronto-striatal organization observed in non-human primates is largely assumed to have 8 been conserved in the human brain (Haber, 2016). However, the complex cognitive and psychiatric 9 symptoms caused by striatal dysfunction point towards human-specific aspects of fronto-striatal 10 connectivity. Further, relative to non-human primates, the prefrontal cortex is expanded in humans 11 (Smaers et al., 2011), and pathways for language processing have developed (Dick and Tremblay,12 2012; Friederici, 2017). It is unknown how human-specific circuits integrate into fronto-striatal 13 organization. 14 Precisely mapping human fronto-striatal connections with non-invasive imaging has been 15 challenging. Techniques such as diffuse optical tomography, electroencephalography, and 16 magnetoencephalography cannot obtain signal from the brain center. MRI-based techniques such 17 as resting-state functional connectivity (RSFC), which involves identifying correlated patterns of 18 activity in spatially distant regions of the brain (Biswal et al., 1995), have proven successful in 19 describing the brain's network organization (Power et al., 2011;Yeo et al., 2011). However, these 20 techniques suffer rapid drop-off of signal-to-noise (SNR) ratios as distance increases from the head 21 coil. Standard approaches for overcoming low SNR in the brain center require registering study 22 participants to a common atlas and averaging their data. Such group-averaged neuroimaging 23 studies have been limited to describing broad, nonspecific connections between large swaths of 24 striatum and widespread regions of cortex, distributed large-scale networks, or entire cortical lobes 25 ( scale networks that comprise strong, functionally differentiated connections between spatially 17 specific cortical and subcortical regions (Gordon et al., 2020). The high spatial specificity of 18 connections that can be identified using this subnetwork-PFM technique makes it ideal for 19 characterizing precise, spatially-specific fronto-striatal connections. 20 Precise characterization of individual-specific human cortico-striatal organization would allow us 21 to address fundamental questions in neuroscience, such as "is brain organization continuous, or 22 discrete?" The rostral-caudal organization of fronto-striatal projections has been termed a 23 "functional gradient" (Haber, 2003), in that it describes a progression along a single rostral-caudal 24 axis. Group-averaged neuroimaging data have been used to present functional cortico-striatal 25 gradients as a continuous transition along a primary organizational axis (Jarbo and Verstynen, 26

Striatal Subnetworks are Individually Specific but Spatially Similar Across Subjects 3
We delineated individual-specific subnetworks connecting striatum to cortex and matched them 4 across subjects following procedures described in (Gordon et al., 2020). Briefly, in each subject we 5 calculated RSFC strength between each pair of points in the brain, and then used a data-driven 6 subnetwork detection procedure designed to identify and operate upon the strongest functional 7 connections of every point in the brain. We then matched subnetworks with striatal and cortical 8 representation across subjects. See Methods for details. 9 Ten separate subnetworks with striatal and cortical representation were identified that were 10 spatially and topologically consistent across subjects, with some inter-individual variability in the 11 localization of each subnetwork. See  Table S1 for information about 13 the consistency of connections. These ten subnetworks existed as substructures within several 14 known large-scale brain networks, including the Default-Mode, Salience, Cingulo-opercular, Fronto-15 parietal, Language, and Somatomotor networks (see Fig 3 for an example subject; see Table S1 for 16 mode network memberships across subjects). The spatial representations of these subnetworks 17 spanned from posterior putamen to nucleus accumbens in striatum, and from central sulcus and 18 parietal/temporal lobes in cortex to subgenual cingulate/orbitofrontal cortex. vertices/voxels are colored based on the most common subnetwork present across the ten 7 subjects. A putative Language network (bolded in the inset legend) was not identified in the Haber 8 model. Note that a Motor subnetwork is present posterior of the displayed slice in putamen, and a 9 second Language network is present anterior of the displayed slice in left caudate. 10 11 12

Fronto-Striatal Organization is Better Explained as a Stepped than a Continuous Gradient 13
Fronto-striatal connectivity is often presented as exhibiting a smooth and continuous rostral-1 caudal gradient (e.g., model in Fig 5A, left). Visual examination of the subnetworks identified here 2 (see Fig 1) does support the presence of a rostral-caudal gradient; however, it also suggests that 3 this gradient may be not smooth, but comprised of a progression of discontinuous cortical areas 4 with discrete connectivity profiles (model in Fig 5A, right). Such a gradient may be conceptualized 5 as "stepped" rather than smooth and continuous. A stepped gradient would be consistent with both 6 the postulated rostral-caudal gradient of fronto-striatal projections (Vogelsang and D'Esposito, 7 2018) as well as with the areal model of cortex (Brodmann, 1909;Churchland and Sejnowski, 8 1988). As a third possibility, both the continuous and stepped models could simultaneously be 9 present, such that regions exhibit categorically distinct fronto-striatal connections from each other, 10 but a rostral-caudal gradient also exists overlaid on top of the discrete regions. Here, we used linear 11 regression models to test whether the locations of strong fronto-striatal connections were best 12 explained by continuous rostral-caudal position, or by discrete subnetwork identities, or by both. 13 We first had to define the rostral-caudal position for each brain location of interest. We defined a 14 separate rostral-caudal axis for the caudate, the left and right putamen, and the prefrontal cortex 15 ( Fig 5B), each with its own unique geometry. We then labeled each striatal voxel and each frontal 16 cortex vertex according to its position along that fit line (i.e., the point along the structure-specific fit 17 line that is closest to that voxel / vertex). The calculated rostral-caudal position of each voxel / 18 vertex can be seen in Fig 5C. Note that we elected to use this anatomically-based rostral-caudal 19 orientation rather than explicitly fitting a gradient to the RSFC data because we did not want the 20 continuous gradient overfit to the data when comparing against a discrete organizational model. 21 for all points within frontal cortex plus anterior insula (left) and striatum (right). These rostral-caudal 7 positions were used to compare models of organization described in (A). Note that while striatal 8 connections with non-frontal sources were observed (Figs 1-3, Fig S2), rostral-caudal modeling 9 here is restricted to the frontal cortex. 10 11 We then conducted two separate, symmetric series of regressions in every subject. First, for 1 every frontal/insular cortical vertex with a subnetwork identity, we calculated the voxel within the 2 striatum that had the strongest functional connectivity with that vertex. As stronger functional 3 connectivity is arguably likely to represent more direct connections, we take this voxel of maximal 4 connectivity to be the striatal target of that cortical vertex. 5 To test the continuous gradient hypothesis, which holds that the rostral-caudal positions of 6 connected cortical and subcortical regions should vary continuously with each other, we regressed 7 the rostral-caudal position of the cortical source vertex against the continuous rostral-caudal 8 position of the striatal target voxel. To test the discrete, stepped gradient hypotheses, we conducted 9 a one-way ANOVA examining whether the rostral-caudal position of the striatal target voxel was 10 explained by the subnetwork identity of the cortical source vertex. In each case, the adjusted R 2 11 value was calculated as the measure of variance explained, in order to allow direct comparisons 12 between these two regressions which had differing numbers of independent variables. It is 13 important to note that the dependent variable in both of these regressions was the rostral-caudal 14 position of the subcortical target voxel. As a result, the adjusted R 2 values are directly comparable 15 between the two; and further, this rostral-caudal position is not a priori known from either the 16 cortical vertex's position or its subnetwork identity. Finally, to test whether both the continuous and 17 discrete models might be simultaneously true, we entered both the cortical position and subnetwork 18 identity into an ANCOVA model testing against subcortical voxel position. 19 We found that both the continuous and the discrete, stepped gradient hypotheses successfully 20 fit the data, but that the discrete stepped hypothesis explained more variance in the data (Fig 6,  When both factors were entered into the same ANCOVA model, only the discrete stepped 3 gradient hypothesis explained unique variance in the data (Fig 6E, bottom). Across subjects, the 4 rostral-caudal position of a cortical vertex explained only 0.6% ± 0.8% of unique variance in the 5 position of its striatal target, after controlling for subnetwork identity. By contrast, the subnetwork 6 identity of the cortical vertex still explained 40% ± 7.5% of unique variance in the position of its 7 striatal target, after controlling for rostral-caudal position. 8 The above regressions were asymmetric: they explained the rostral-caudal position of a 9 subcortical target based on the position and/or subnetwork identity of its cortical source. To ensure 10 symmetry of these effects, we repeated all of the regressions in the reverse direction, in which the 11 most strongly connected cortical source vertex was calculated for each striatal voxel, and the 12 rostral-caudal position of a cortical source vertex was explained based on the position and/or 13 subnetwork identity of the striatal voxel. 14 In these regressions, once again both the continuous and the discrete, stepped gradient 15 hypotheses successfully fit the data, but the discrete stepped hypothesis explained more variance 16 in the data (Fig 6 identity explaining more variance than rostral-caudal position in 9/10 subjects. 21 When both factors were entered into one model, again only the discrete stepped gradient 22 hypothesis explained unique variance in the data (Fig 6F, bottom), with the rostral-caudal position 23 of a subcortical voxel explaining only 3% ± 4% of unique variance after controlling for subnetwork 24 identity, and the subnetwork identity of the voxel still explaining 22% ± 12% of unique variance. 25 These findings indicate that the rostral-caudal progression of strong fronto-striatal connections is 1 organized as a stepped gradient composed of a series of discrete subnetworks. of its most strongly connected striatal target (y-axis), and vice versa (B). However, the discrete 5 subnetwork identity of a cortical vertex (delineated by color and by x-axis position) explained more 6 variance in the rostral-caudal position of its striatal target (y-axis), and vice versa (D). In A-D, each 7 dot represents one cortical vertex (A,C) or striatal voxel (B,D). Across subjects, for both the 8 cortex→striatum regressions (E) and the striatum→cortex regressions (F), both continuous rostral-9 caudal position and discrete subnetwork identity explained substantial variance, but subnetwork 10 identity explained more variance in every subject (top). When both factors were employed in the 11 same model (bottom), the discrete subnetwork identity still explained substantial unique variance, 12 but the continuous gradient factor explained almost no unique variance. 13 14 15 The progression of subnetworks along the rostral-caudal axis is more striking on the medial 16 surface of the frontal cortex than on the lateral surface (e.g. Fig 1). It is possible that the relative 17 explanatory power of continuous vs stepped gradients might vary based on whether projections are 18 from the medial or lateral surface. To explore this possibility, we repeated the regressions above 19 twice: once while restricting cortical source vertices to only the medial, and then to the lateral, 20 frontal cortical surfaces. We found that on the medial surface, adjusted R 2 values were much higher 21 for both the continuous (cortex→striatum: 33% ± 17%; striatum→cortex: 41% ± 12%) and the 22 stepped (cortex→striatum: 61% ± 13%; striatum→cortex: 66% ± 12%) gradient models, but the 23 stepped model explained more variance than the continuous model in every subject, in both 24 directions. On the lateral surface, adjusted R 2 values for the continuous gradient (cortex→striatum: 25 18% ± 24%; striatum→cortex: 24% ± 14%) were similar to those observed across the whole frontal 26 cortex, but adjusted R 2 values for the stepped gradient were slightly lower (cortex→striatum: 50% ± 27 13%; striatum→cortex: 38% ± 14%); however, the stepped model explained more variance than the 28 continuous model in every subject in the striatum→cortex direction, and in 9/10 subjects in the 29 cortex→striatum direction. These results suggest that the rostral-caudal organization of fronto-30 striatal connections, which is organized as a stepped gradient, is more prominent in the medial than 31 the lateral frontal cortex. 32

DISCUSSION 1
In this work, we identified ten subnetworks representing very strong, anatomically precise, 2 individually-specific connections between cortex and striatum that could be matched across 3 subjects. The identification of these connections within individual humans represents a substantial 4 breakthrough in our ability to understand these circuits that are both critical for complex behaviors 5 and implicated in many neuropsychiatric disorders. The neurobiological validity of the observed 6 cortico-striatal subnetwork structures is supported by their general convergence with findings from 7 invasive, anatomically precise tract-tracing studies in non-human primates. Critically, the identified 8 subnetworks also diverge from non-human primate-based models of fronto-striatal connectivity in 9 several ways, most prominently by identifying left-lateralized striatal elements of the human 10 language network. 11 12

Most Cortico-striatal Subnetworks Converge with Primate Tract-tracing 13
We identified striatal subnetworks composed of only very strong functional connections in single 14 individuals, and we found that 1) the majority of these very strong cortico-striatal connections were 15 from frontal cortex, and 2) frontal regions of these subnetworks had stronger connectivity with 16 striatum than non-frontal regions, suggesting that the non-frontal regions we did observe were less 17 strongly associated with striatal activity. This converges with non-human primate work arguing that, The identified fronto-striatal circuits were highly ordered as a topographically mirrored rostral-3 caudal gradient of cortico-striatal projections. This organization converged with the organization of 4 these circuits proposed based on findings from non-human primates (Haber, 2016(Haber, , 2003; see anterior cingulate to anterior caudate and putamen, and particularly to sites more dorsal and lateral 26 than the orbitofrontal cortex (OFC) or ventromedial prefrontal cortex (vmPFC) projections described 1 below (Haber et al., 2006). This set of cortical regions (dorsal anterior cingulate, anterior insula, 2 anterior striatum) is considered the core of the adjacent "Cingulo-opercular" and "Salience" 3 networks, which are critical for a wide variety of cognitive control operations 4 5) We found functional connections between nucleus accumbens and pregenual medial 5 prefrontal cortex (mPFC), and a separate set with OFC / subgenual cingulate cortex ( respectively. The nucleus accumbens is anatomically segregated into a more dorsal shell and a 9 more ventral core (Salgado and Kaplitt, 2015). Here, the nucleus accumbens nodes of these two 10 subnetworks exhibited a similar dorsal/ventral distribution; as such, we tentatively label these nodes 11 the accumbens shell and core, respectively. 12 Interestingly, these connections only partially converged with connectivity identified in non-13 human primates, in which the nucleus accumbens core preferentially receives projections from 14 pregenual mPFC, while the shell receives projections from OFC-the reverse of the arrangement 15 we observed here. It is possible that this divergence from non-human primate work is driven by 16 fMRI data quality issues, such as the known presence of susceptibility artifact within ventral 17 prefrontal and orbitofrontal cortex in human fMRI. However, human diffusion tractography imaging 18 has also demonstrated that the accumbens core is more strongly connected to OFC than the shell 19 consider an alternate possibility that this divergence represents a true organizational difference 21 between humans and non-human primates related to reward processing. 22 Several of the connections identified here, including to the posterior and middle putamen, dorsal 23 caudate, and nucleus accumbens core, were restricted to a single target subcortical structure. This 24 represents an increase in target specificity relative to the (Haber, 2016) model, in which similar 25 fronto-striatal projections span the internal capsule (Fig 4A,B, green and orange) or extend out of 26 nucleus accumbens into caudate (Fig 4A,B, red). We consider it most likely that frontal projections 1 to multiple striatal structures do exist in humans, but that projections to secondary structures do not 2 drive striatal function as strongly as projections to the primary target structure. 3 4 Two Cortico-Striatal Subnetworks Diverged from Primate Models but Converged with Human 5

Language Networks 6
We observed two subnetworks that did not converge with proposed models of non-human 7 primate fronto-striatal connectivity (Haber, 2016). These consisted of two adjacent sets of 8 Prior work has demonstrated that fronto-striatal connections are organized as mirrored cortical 9 and striatal gradients spanning from primary motor cortex / posterior putamen to nucleus 10 accumbens / ventromedial prefrontal cortex. This finding is consistent across humans and non-11 human primates, and across methodologies ranging from invasive tract-tracing to diffusion MRI to 12 RSFC to task-based fMRI (Haber, 2003 Fronto-striatal gradients are often presented, for visual or mathematical convenience, as if they 20 were smooth and continuous. However, the concept of a smooth, continuous fronto-striatal gradient 21 conflicts with evidence that the brain is organized as discrete, discontinuous areas networked Importantly, a progressive rostral-to-caudal series of such discontinuous cortical areas with 5 topographically ordered striatal projections-a stepped gradient-can visually ( Fig 5A) and 6 statistically (Fig 6A,B) appear continuous, particularly if the categorical divisions are unknown. 7 Here, direct comparison of continuous gradient models against discontinuous, stepped gradient 8 models demonstrated that the stepped gradient models explained the data much better (Fig 6). This 9 indicates that there is little evidence for a smooth, continuous cross-areal gradient, either as a 10 primary model of human fronto-striatal organization or as a secondary effect superimposed on top 11 of a stepped gradient. Instead, the data is best explained as a stepped rostral-caudal gradient 12 composed of a progressive series of discrete fronto-striatal subnetworks. Interestingly The present approach for comparing continuous and stepped gradient models does have 5 limitations. First, the infomap technique, while excellent at recovering true networked connections, 6 is not always stable in determining when subnetworks should be combined or split. As a result, the 7 categorical factor entered into the regression model may sometimes have too many or too few 8 categories. Second, using the rostral-caudal location of the strongest functional connection as a 9 dependent variable ignores the possibility that even very strong cortico-striatal connections are 10 many-to-one (as is suggested by the identified subnetworks). Thus, the model may not accurately 11 represent the brain's true connectional organization. Importantly, both of these issues would be 12 expected to reduce the variance explained of the stepped gradient model, but would not affect the 13 continuous gradient model. The fact that the stepped model outperforms the gradient model despite 14 these limitations is strong evidence that it more accurately represents the organization of fronto-15 striatal connections. 16 17

Integration 19
It is well established in rodents and non-human primates that cortico-striatal circuits are further 20 projected to the thalamus, and then back to cortex to form feedback loops (Haber, 2003;Nakano et 21 al., 2000). Thus, in principle the cortico-striatal subnetwork structures we observe here should also 22 have representation in the thalamus. While we did observe this representation to some degree, the 23 thalamus did not consistently contain representations of all cortico-striatal subnetworks; and further, 24 the thalamus tended to have extensive representation of subnetworks characterized as noise-25 related. This suggests that the quality of the MSC fMRI data, while adequate for winner-take-all 26 style whole-network mapping in thalamus (Greene et al., 2020), may not be adequate for the data-1 driven mapping of highly specific cortico-striato-thalamo-cortical circuits in thalamus. Future work 2 may be able to more broadly conduct such mapping using such higher-quality data collected with 3 advanced multiband sequences. 4 While the presented delineations of cortico-striatal circuits are neurobiologically compelling and 5 individual-specific, they are incomplete. In examining only the very strongest functional connections 6 to each striatal voxel, we focus on the likely driving inputs to the striatum, but we also exclude the 7 wealth of multi-lobe, modulatory inputs to striatum that have been described in humans (Choi et al., demonstrated that fronto-striatal circuits are only partially segregated, with some degree of cross-12 circuit connections (Haber, 2003). Similarly, at the large-scale network level, human RSFC work 13 has also demonstrated cross-network integration zones (Greene et al., 2020). In reproducing 14 analyses from (Greene et al., 2020) at the subnetwork level rather than the network level, we 15 similarly observed extensive cross-subnetwork connections (Fig S4). Untangling the substantial 16 complexity of these multi-subnetwork connections for comprehensive evaluation and categorization 17 must be conducted in future work. 18 19

Mapping Discrete Cortico-Striatal Circuits in Individual Humans 20
In describing cortico-striatal organization as a stepped gradient composed of discrete 21 subnetworks, the present results have important implications for the future study of the human 22 striatum. The discrete subnetworks present within the striatum likely enable a variety of different 23 cognitive and motor processes that can be localized to each separate, individual-specific striatal 24 network structure. With such localization, we can generate strong hypotheses about how individual-25 specific symptoms observed in neurosychiatric disorders might arise from disruption of those 26 specific fronto-striatal circuits, or about how neurodegenerative diseases or brain injuries that 1 damage specific fronto-striatal circuits would impair associated cognitive and motor functions. We 2 can then explore whether symptoms can be alleviated or impaired function restored via brain 3 stimulation approaches, including noninvasive TMS targeting the cortical sources of these fronto-4 striatal projections, or deep electrode implantation targeting their striatal targets. 5 6

2
Subjects 3 Data were collected from ten healthy, right-handed, young adult subjects (5 females; age: 24-4 34). Two of the subjects are authors (NUFD and SMN), and the remaining subjects were recruited 5 from the Washington University community. Informed consent was obtained from all participants. 6 The study was approved by the Washington University School of Medicine Human Studies 7 Committee and Institutional Review Board. Other findings using these participants have been 8 previously reported in (Gordon et  images. An EyeLink 1000 eye-tracking system (http://www.sr-research.com) allowed continuous 30 monitoring of subjects' eyes in order to check for periods of prolonged eye closure, potentially 31 indicating sleep. Only one subject (MSC08) demonstrated prolonged eye closures. 32 33 Cortical surface generation 34 Generation of cortical surfaces from the MRI data followed a procedure similar to that previously Matching (MSM) algorithm using sulc features to register templates to the atlas mesh (Robinson et  47 al., 2014). These newly registered surfaces were then down-sampled to a 32,492 vertex surface 48 (fs_LR 32k) for each hemisphere. The various structural metric data (thickness, curvature, etc.) 49 from the original surfaces to the fs_LR 32k surface were composed into a single deformation map 50 allowing for one step resampling. MSM registration provided a more optimal fit of pial and white 51 surfaces and reduced areal distortion (Glasser et al., 2016b). These various surfaces in native 1 stereotaxic space were then transformed into atlas space (711-2B) by applying the previously 2 calculated T1-to-atlas transformation. 3 4 fMRI Preprocessing 5 Functional data were preprocessed to reduce artifacts and to maximize cross-session 6 registration. All sessions underwent correction of odd vs. even slice intensity differences attributable 7 to interleaved acquisition, intensity normalization to a whole brain mode value of 1000, and within 8 run correction for head movement. Atlas transformation was computed by registering the mean 9 intensity image from a single BOLD session to Talairach atlas space (Talairach and Tournoux,  10 1988) via the average high-resolution T2-weighted image and average high-resolution T1-weighted 11 image. All subsequent BOLD sessions were linearly registered to this first session. This atlas 12 transformation, mean field distortion correction (see below), and resampling to 2-mm isotropic atlas 13 space were combined into a single interpolation using FSL's applywarp tool (Smith et al., 2004). All 14 subsequent operations were performed on the atlas-transformed volumetric time series. 15 16 Distortion correction 17 A mean field map was generated based on the field maps collected in each subject (Laumann 18 et al., 2015). This mean field map was then linearly registered to each session and applied to that 19 session for distortion correction. To generate the mean field map the following procedure was used: 20 (1) Field map magnitude images were mutually co-registered. (2) Transforms between all sessions 21 were resolved. Transform resolution reconstructs the n-1 transforms between all images using the 22 n(n-1)/2 computed transform pairs. (3) The resolved transforms were applied to generate a mean 23 magnitude image. (4) The mean magnitude image was registered to an atlas representative 24 template. (5) Individual session magnitude image to atlas space transforms were computed by 25 composing the session-to-mean and mean-to-atlas transforms. (6) Phase images were then 26 transformed to atlas space using the composed transforms, and a mean phase image in atlas 27 space was computed. 28 Application of mean field map to individual fMRI sessions: (1) For each session, field map 29 uncorrected data were registered to atlas space, as above. (2) The generated transformation matrix 30 was then inverted and applied to the mean field map to bring the mean field map into the session 31 space. (3) The mean field map was used to correct distortion in each native-space run of resting 32 state and task data in the session. (4) The undistorted data were then re-registered to atlas space. 33 (5) This new transformation matrix and the mean field map then were applied together to resample 34 each run of resting state and task data in the session to undistorted atlas space in a single step. 35  After computing the temporal masks for high motion frame censoring, the data were processed 49 with the following steps: (i) demeaning and detrending, (ii) linear interpolation across censored 50 frames using so that continuous data can be passed through (iii) a band-pass filter (0.005 Hz < f < 51 0.01 Hz) without re-introducing nuisance signals (Hallquist et al., 2013) or contaminating frames 1 near high motion frames (Carp, 2013 and cerebrospinal fluid (CSF) assumes that variance in such regions is unlikely to reflect neural 5 activity. Variance in these regions is known to correspond largely to physiological noise (e.g., CSF 6 pulsations), arterial pCO2-dependent changes in T2*-weighted intensity and motion artifact; this 7 spurious variance is widely shared with regions of interest in gray matter. We also included the 8 mean signal averaged over the whole brain as a nuisance regressor. Global signal regression 9 (GSR) has been controversial. However, the available evidence indicates that GSR is a highly 10 effective de-noising strategy (Ciric et al., 2017;Power et al., 2015). 11 Nuisance regressors were extracted from white matter and ventricle masks, first segmented by 12 FreeSurfer (Fischl, 2012), then spatially resampled in register with the fMRI data. Voxels  13 surrounding the edge of the brain are particularly susceptible to motion artifacts and CSF pulsations 14 (Patriat et al., 2015;Satterthwaite et al., 2013); hence, a third nuisance mask was created for the 15 extra-axial compartment by thresholding the temporal standard deviation image (SDt > 2.5%), 16 excluding a dilated whole brain mask. Voxel-wise nuisance time series were dimensionality reduced 17 as in CompCor (Behzadi et al., 2007), except that the number of retained regressors, rather than 18 being a fixed quantity, was determined, for each noise compartment, by orthogonalization of the 19 covariance matrix and retaining components ordered by decreasing eigenvalue up to a condition 20 number of 30 (max eigenvalue / min eigenvalue > 30). The retained components across all 21 compartments formed the columns of a design matrix, X, along with the global signal, its first 22 derivative, and the six time series derived by retrospective motion correction. The columns of X are 23 likely to exhibit substantial co-linearity. Therefore, to prevent numerical instability owing to rank-24 deficiency during nuisance regression, a second-level SVD was applied to XX T to impose an upper 25 limit of 250 on the condition number. This final set of regressors was applied in a single step to the 26 filtered, interpolated BOLD time series, with censored data ignored during beta estimation. Finally, 27 the data were upsampled to 2mm isotropic voxels. Censored frames were then excised from the 28 data for all subsequent analyses. 29 30 Surface processing and CIFTI generation of BOLD data 31 Surface processing of BOLD data proceeded through the following steps. First, the BOLD fMRI 32 volumetric timeseries were sampled to each subject's original mid-thickness left and right-33 hemisphere surfaces (generated as the average of the white and pial surfaces) using the ribbon-34 constrained sampling procedure available in Connectome Workbench 1.0. This procedure samples 35 data from voxels within the gray matter ribbon (i.e., between the white and pial surfaces) that lie in a 36 cylinder orthogonal to the local mid-thickness surface weighted by the extent to which the voxel falls 37 within the ribbon. Voxels with a timeseries coefficient of variation 0.5 standard deviations higher 38 than the mean coefficient of variation of nearby voxels (within a 5 mm sigma Gaussian 39 neighborhood) were excluded from the volume to surface sampling, as described in (Glasser et al., 40 2013). Once sampled to the surface, timecourses were deformed and resampled from the 41 individual's original surface to the 32k fs_LR surface in a single step using the deformation map 42 generated above (in "Cortical surface generation"). This resampling allows point-to-point 43 comparison between each individual registered to this surface space. 44 These surfaces were then combined with volumetric subcortical and cerebellar data into the 45 CIFTI format using Connectome Workbench (Marcus et al., 2011), creating full brain timecourses 46 excluding non-gray matter tissue. Subcortical (including accumbens, amygdala, caudate, 47 hippocampus, pallidum, putamen, and thalamus) and cerebellar voxels were selected based on the 48 FreeSurfer segmentation of the individual subject's native-space average T1, transformed into atlas 49 space, and manually inspected. Finally, the BOLD timecourses were smoothed with a geodesic 2D 50 (for surface data) or Euclidean 3D (for volumetric data) Gaussian kernel of σ = 2.55 mm. 51 1 Regression of adjacent cortical tissue from RSFC BOLD 2 Some striatal regions, particularly including lateral putamen, are in close anatomical proximity to 3 cortex, resulting in spurious functional coupling between the cortical vertices and adjacent 4 subcortical voxels. To reduce this artifact, RSFC BOLD time series from all vertices falling within 5 20mm Euclidean distance of a source voxel were averaged and then regressed from the voxel time timeseries were used for all subsequent analyses. 8 9 Mapping multi-scale network structure 10 The network organization of each subject's brain was delineated following (Gordon et al., 2020) 11 using the graph-theory-based Infomap algorithm for community detection (Rosvall and Bergstrom, 12 2008). In this approach, we calculated the cross-correlation matrix of the time courses from all brain 13 vertices (on the cortical surfaces) and voxels (in subcortical structures), concatenated across 14 sessions. Correlations between vertices/voxels within 30 mm of each other were set to zero in this 15 matrix to avoid basing network membership on correlations attributable to spatial smoothing. 16 Geodesic distance was used for within-hemisphere surface connections and Euclidean distance for 17 subcortical-to-cortical connections. Connections between subcortical structures were disallowed, as 18 we observed extremely high connectivities within nearly the entire basal ganglia that would prevent 19 network structures from emerging. Inter-hemispheric connections between the cortical surfaces 20 were retained, as smoothing was not performed across the mid-sagittal plane. 21 We observed that connectivity patterns within regions known to have low BOLD signal due to 22 susceptibility artifact dropout (e.g. ventral anterior temporal lobe, portions of orbitofrontal cortex) 23 were effectively random. To avoid having the delineated network structures distorted by such 24 random connections, we first calculated a set of common low-signal regions as the vertices in which 25 the average mode-1000 normalized BOLD signal across subjects and timepoints was less than 750 26 (as in Gordon et al., 2016;Wig et al., 2014). All connections to these low-signal regions were set to 27 zero. 28 This matrix was then thresholded to retain at least the strongest .1% of connections to each 29 vertex and voxel, following observations in (Gordon et al., 2020) that this level of network division 30 represents divisions that, across these subjects, optimally divide the resting state data and explain 31 task activation patterns. 32 These thresholded matrices were used as inputs for the Infomap algorithm, which calculated 33 community assignments (representing brain network structures) separately for each threshold. 34 Small networks with 10 or fewer vertices / voxels were considered unassigned and removed from 35 further consideration. The above analysis was conducted in each individual subject. 36 37 Identifying matched striatal subnetworks in individuals 38 For each subject, we considered only subnetworks that had at least some representation in the 39 striatum (the individually defined caudate, putamen, nucleus accumbens, and globus pallidus). 40 Following (Gordon et al., 2020), we visually examined the cortical and subcortical topographies of 41 each subnetwork with striatal representation, as well as their topological arrangement relative to 42 each other. After careful consideration of the subnetworks observable in this population, matched 43 subnetworks were identified based on the following heuristic rules: 44 A subnetwork in inferior ventral striatum usually also had representation in subgenual cingulate 45 and/or orbitofrontal cortex. 46 A subnetwork in anterior ventral striatum had representation in bilateral pregenual cingulate 47 cortex. This subnetwork was previously described in (Gordon et al., 2020). 48 A subnetwork in superior ventral striatum had representation in anterior cingulate cortex and 49 middle insula. 50 A subnetwork that was spatially variable within striatum, with representation most often in 1 medial putamen, but was consistently present in dorsal anterior cingulate cortex and dorsal insula. 2 A subnetwork in lateral middle putamen and posterior dorsomedial prefrontal cortex. 3 A subnetwork in dorsal caudate that had consistent representation as the dorsomedial cluster of 4 the Fronto-Parietal network and inconsistent representation in lateral frontal and parietal cortex. 5 A subnetwork in medial anterior putamen (especially left putamen) that had consistent cortical 6 representation converging with the Language network (Braga et al., 2020). 7 A subnetwork in lateral anterior caudate (especially left caudate) that had representation to a set 8 of cortical regions directly anterior to the Language network above. 9 A subnetwork (sometimes multiple subnetworks) in ventral posterior putamen that had 10 representation somewhere in the central sulcus, with the location of the representation varying 11 substantially by subject. 12 Additionally, we commonly observed subnetworks that were generally in posterior lateral 13 putamen (but with variable representation across subjects) that had cortical representation only in 14 middle insular cortex, physically adjacent to the striatal representation, or (in one subject) in 15 unusual locations such as medial occipital lobe. These subnetworks also often cut implausibly 16 across disjoint anatomical structures such as putamen and thalamus. These subnetworks were 17 tentatively classified as noise-related and were ignored for future analyses. 18 19 Computing lobe-wise striatal subnetwork representation 20 For each individual subject, we defined the Frontal, Insular, Parietal, Temporal, and Occipital 21 lobes based on the lobe identities of the regions of the subject-specific Freesurfer-generated 22 Destrieux atlas parcellation, which was deformed into fs_LR_32k space to match the functional 23 data. Lobe-wise representation was computed within each subject as the percent of cortical vertices 24 among all non-noise subnetworks that were within each lobe. 25 26 Visualizing subnetwork overlap across subjects 27 For each matched non-noise subnetwork, the overlap across individuals was visualized using 28 Connectome Workbench. For cortex, the number of individuals with the subnetwork present was 29 calculated in each vertex. For striatum, this was less straightforward, as the subcortical structures 30 were individualized, and so were different shapes and sizes across people. To represent overlap in 31 striatum, we first mapped the subcortical voxels of each individual among MSC02-10 to the 32 subcortical structures in MSC01. This was done by assigning each voxel in the MSC01 structures 33 the network ID of the nearest (by Euclidean distance) subcortical voxel in that individual. Overlap 34 was then calculated across MSC01-space subject subnetworks as the number of individuals with 35 the subnetwork present in each MSC01-space voxel. 36 37 Defining position along a rostral-caudal axis in striatum and frontal cortex 38 For each subject, we calculated four separate curves to fit the rostral-caudal orientation of 1) 39 bilateral caudate, 2) left putamen, 3) right putamen, and 4) frontal cortex. 40 For caudate, the structure is oriented vertically, and the rostral-caudal organization is best 41 represented as a curve running from ventral striatum up through caudate head and back to caudate 42 tail-that is, a curve in the y and z dimensions, with the x dimension irrelevant. We thus created a 43 rostral-caudal axis as a lowess curve fit to the y-and z-coordinates of the caudate voxels. This 44 curve represents the approximate geometric center of the structure in the y and z dimensions. In 45 the dorsal posterior portion of the caudate (body and tail), the rostral-caudal organization is simply 46 anterior-to-posterior. Therefore, once the curve reached the maximum z coordinate of the structure, 47 the lowess fit was terminated and the axis was continued by projecting straight backwards along the 48 y axis. For this fitting, nucleus accumbens was combined with caudate. Additionally, left and right 49 caudate were combined, as the difference between the left and right structures (in the x-axis) was 50 irrelevant for the fit. 51 For putamen, the structure is oriented horizontally, and thus the rostral-caudal organization is 1 best represented as a curve in the x and y dimensions, with the z dimension irrelevant. Thus, as 2 with caudate, we created a rostral-caudal axis as a lowess curve fitted to the x-and y-coordinates 3 of the putamen voxels, representing the approximate geometric center of the structure in the x and 4 y dimensions. Once the curve reached the maximum x coordinate of the structure, the lowess fit 5 was terminated and the axis was continued by projecting straight backwards along the y axis. 6 Because the left and right putamen are mirrored in the x dimension (rather than oriented 7 equivalently, as with the caudate z dimension), the fit was necessarily conducted separately for the 8 left and right structures. For this analysis, globus pallidus was not included in the fit. 9 For frontal cortex, the structure is oriented vertically, as with caudate. Therefore, as with 10 caudate, we created a rostral-caudal axis as a lowess curve fit to the y-and z-coordinates of the 11 cortical vertices in the frontal cortex. For each subject, frontal cortex was defined using the subject-12 specific Destrieux atlas parcellation, as above, except that we also included anterior insula regions, 13 which commonly included substantial representation of striatal networks (excluding insula regions 14 did not change the results; see Fig S5B). As with caudate, we terminated the curve and projected 15 backwards along the y axis once the curve reached the maximum z coordinate. 16 Finally, we calculated the coordinate of each striatal voxel and each cortical vertex along the 17 rostral-caudal axis. First, the distances along the four rostral-caudal fit curve axes were converted 18 to rank value (rank within each curve) to ensure comparability across structures. Then, for each 19 voxel / vertex, the rostral-caudal coordinate was taken to be the rank of the point on the relevant 20 structure's rostral-caudal axis closest (i.e. minimum Euclidean distance) to that voxel/vertex. Note 21 that this calculation-the closest point on the relevant axis-is identical to how x/y/z axis 22 coordinates are defined in a standard 3D space. 23 In each case, rostral-caudal positions were converted to rank order along the rostral-caudal 24 gradient to account for nonuniform effects of position. 25 26 Examining continuous and discrete organizations of strong fronto-striatal connections 27 We conducted two separate sets of analyses: One in which we identified the most strongly 28 connected striatal target voxel for each frontal cortex source vertex, and the second in which we 29 identified the most strongly connected frontal source vertex for each striatal target voxel. 30 Regression Set 1: Cortex→Striatum 31 For each frontal cortex vertex, we identified the striatal voxel with which it had the strongest 32 functional connection. This was done by calculating correlations between the timecourses of all 33 subcortical voxels and all frontal cortex vertices (frontal cortex defined as above). For each frontal 34 vertex, the striatal voxel with the most strongly correlated timecourse was taken to be that vertex's 35 "target"-the voxel most likely to be directly anatomically connected to that voxel. Voxels with 36 correlations to their sources that were weaker than R=.2 were considered too weak to be reliable 37 and were excluded from subsequent analyses. However, inclusion of these voxels does not change 38 the results; see Fig S5A. 39 We tested which factors of a frontal source vertex explain the rostral-caudal position of its 40 striatal target. We first tested whether fronto-striatal connections could appear to be a continuous 41 rostral-caudal gradient, in which continuous movement from vertex to vertex along the frontal cortex 42 rostral-caudal axis is mirrored by a continuous rostral-caudal movement of those vertices' striatal 43 targets. Specifically, separately in each subject, we tested whether the rostral-caudal position of 44 cortical vertices was correlated with the rostral-caudal position of their striatal targets. We evaluated 45 this effect as the adjusted square of the correlation value (R 2 ), representing the variance in striatal 46 voxel position explained by cortical target position, corrected for the use of the single regressor. 47 We then tested whether fronto-striatal connections could appear as a stepped gradient 48 composed of discrete networks, in which the location of a striatal target voxel is explained by the 49 network of its frontal source vertex. Specifically, separately in each subject, we conducted a one-50 way ANOVA with the network identity of frontal vertices as an independent variable and the rostral-51 caudal position of their striatal targets as the dependent variable. We evaluated this effect as the 1 adjusted R 2 of the main factor of the ANOVA, representing the variance in striatal target position 2 explained by frontal vertex network ID, corrected for the number of levels in the ANOVA. Note that 3 this adjusted R 2 value is directly comparable to the adjusted R 2 value from the correlation above. 4 Finally, we tested whether fronto-striatal connections could appear as a continuous or a stepped 5 rostral-caudal gradient when accounting for the other organization, as well as which organization 6 explained more variance. Separately in each subject, we conducted a one-way ANCOVA with the 7 network identity of frontal vertices as one discrete independent variable, the rostral-caudal position 8 of frontal vertices as a second, continuous variable, and the rostral-caudal position of the striatal 9 targets as the dependent variable. The explanatory power of each factor (continuous gradient vs 10 discrete network) was compared against each other by calculating and computing the adjusted R 2 11 of each factor within the ANCOVA. 12 13 Regression Set 2: Striatum→Cortex 14 These regressions were conducted using the same logic as above, but using the striatal voxels 15 as independent variables and the frontal vertices as dependent variables. Specifically, for each 16 striatal voxel, we identified the frontal cortex vertex with which it had the strongest functional 17 connection. This was done as above by calculating correlations between the timecourses of all 18 subcortical voxels and all frontal cortex vertices. For each voxel, the cortical vertex with the most 19 strongly correlated timecourse was taken to be its "source". Voxels with correlations to their sources 20 that were weaker than R=.2 were excluded from subsequent analyses (again, inclusion of these 21 voxels does not change results; Fig S5A). 22 We tested which factors of a striatal target voxel explain the rostral-caudal position of its frontal 23 source. We first tested whether the rostral-caudal position of striatal voxels was correlated with the 24 rostral-caudal position of their frontal sources. We evaluated this effect as the adjusted square of 25 the correlation value (R 2 ), representing the variance in fontal vertex position explained by striatal 26 voxel position, corrected for the use of the single regressor. 27 We then tested whether the location of a frontal source vertex is explained by the network of its 28 striatal target voxel by conducting a one-way ANOVA with the network identity of striatal voxels as 29 an independent variable and the rostral-caudal position of their frontal sources as the dependent 30 variable. We evaluated this effect as the adjusted R 2 of the main factor of the ANOVA, representing 31 the variance in frontal source position explained by striatal voxel network ID, corrected for the 32 number of levels in the ANOVA. This adjusted R 2 value is directly comparable to the adjusted R 2 33 value from the correlation above. 34 Finally, we tested whether fronto-striatal connections could appear as a continuous or a stepped 35 rostral-caudal gradient when accounting for the other organization. Separately in each subject, we 36 conducted a one-way ANCOVA with the network identity of striatal voxels as one discrete 37 independent variable, the rostral-caudal position of striatal voxels as a second, continuous variable, 38 and the rostral-caudal position of the frontal source vertices as the dependent variable. The 39 explanatory power of each factor was compared against each other by calculating and computing 40 the adjusted R 2 of each factor within the ANCOVA. 41 42 43 44 DATA AND SOFTWARE AVAILABILITY 45 Data: Raw MRI data from the MSC Dataset, as well as segmented cortical surfaces, 46 preprocessed volumetric and cifti-space RSFC timecourses, and preprocessed task timecourses 47 and contrasts, have been deposited in the Openneuro data repository 48 (https://openneuro.org/datasets/ds000224/versions/1.0.2) under the label "Midnight Scan Club". 49 Code: Code to perform all preprocessing and analysis is available at 50 https://github.com/MidnightScanClub. 51    Figure S4: Cross-subnetwork integration in subject MSC01. Voxels with preferential RSFC to one network (network-specific) are represented by solid colors, and voxels functionally connected to multiple networks (integrative) are represented by cross-hatching. Analyses performed following procedures in Greene et al. (2020), except that we evaluated voxel RSFC to matched fronto-striatal subjects rather than to whole networks. Figure S5: Fronto-striatal organization explained by continuous and stepped gradients after alternate analysis procedures. A) With no minimum RSFC thresholding (minimum RSFC set to Z(r)=.2 in main text analyses), discrete network identity explained more variance than continuous rostral-caudal position for both the Cortex→Striatum (left) and the Striatum→Cortex (right) analyses. B) After excluding anterior insula voxels from consideration, discrete network identity explained more variance than continuous rostral-caudal position for both the Cortex→Striatum (left) and the Striatum→Cortex (right) analyses.