Organization of the human cerebral cortex estimated within individuals: networks, global topography, and function

1The primate cerebral cortex is populated by specialized networks that support sensory, motor, and higher-order cognitive and affective functions. Characterizing how the networks and their interconnected regions are organized on the cortical surface began more than a century ago with landmark studies of myelogenetic and architectonic patterns (e.g., Refs. –) and continued with modern systems neuroscience integration of anatomical projection data and insights from study of brain lesions (e.g., Refs. –). Over the past decades our laboratory, and the field more broadly, has undertaken data collection efforts and analyses of neuroimaging data with the goal of improving understanding of human network organization and providing noninvasive approaches to measure brain organization for clinical use.

2It is beyond the present scope to cover the extensive literature that has evolved, but it is important to interpret the current effort with awareness that new details and revisions emerge incrementally as the methods and data quality progress. Our efforts presented here reflect another step in that progression. The specific challenge in examining the details of network organization in humans is that the methods are indirect and limited and often noisy within individuals. Despite limitations, advances in structural, diffusion, and functional MRI (fMRI) provide valuable information about human cortical organization, albeit with ambiguities consistent with the complexity of cortical architecture and the low resolution of the techniques.

3Resting-state functional connectivity MRI (fcMRI), based on measuring spontaneous correlated fluctuations between brain regions, has been especially useful for estimating networks (Ref.; see also Refs. –). Explorations in group-averaged fcMRI data, with sample sizes that range from 10 to thousands of participants, reveal network estimates that are consistent across analytical approaches and datasets (e.g., Refs. –). Moreover, estimated networks show similarities to directly observed anatomical projection patterns from tracer injections in the monkey, providing support that they reflect, to a first approximation, anatomically connected networks. Correspondence is far from perfect, and there are unresolved aspects to how indirect human network estimates link to anatomy, a theme that we return to in discussion.

4A recent advance in the field is to use within-individual estimates of networks without recourse to averaging across participants. Architectonic fields tile the cortical mantle with variability in their exact locations, sizes, and borders between individuals. Spatial blurring, inherent in group averaging, impedes the ability to estimate details of network organization. Precision neuroimaging, involving intensive sampling and analysis of data within the individual, preserves idiosyncratic anatomical features.

5Within-individual approaches have been the mainstay in human neuroimaging studies of sensory and motor systems (e.g., Refs. –) and emerged later as a viable method to estimate task-based responses in higher-order association cortex (e.g., Refs. –). Within-individual precision mapping using fcMRI only became emphasized recently, even though the first report was within individuals. Following a landmark demonstration that intensive repeat scanning is possible, multiple groups have pursued within-individual characterization of network organization (e.g., Refs. –; for further discussion see Refs., ).

6Here we continue the investigation of the detailed organization of the cerebral cortex, using within-individual approaches. There are multiple goals and methodological innovations that steer this work. First, we employ deep, intensive imaging to boost the signal-to-noise ratio (SNR) within individual participants. Each new participant was scanned on at least eight separate occasions and often more. Second, we applied a novel multisession hierarchical Bayesian model (MS-HBM; see Ref. ) to automatically estimate networks in the intensively sampled participants. Specifically, the number of networks estimated was set at 15 to capture established networks sometimes missed in simpler network parcellations, as detailed within methods. Third, to enable clinical translational research, we developed an empirical method and projected all network estimates from the surface back into the native-space volume of individual participants, as is needed for presurgical planning and neuromodulation. Fourth, inspired by the possibility to chart global spatial relations between networks (e.g., Ref. ), we also plotted the resulting network estimates on the fully flattened cortical surface. As the results reveal, there are repeating patterns of spatial juxtapositions among networks that provide insight into their evolutionary and developmental origins. Finally, we collected and examined task data within the same intensively sampled participants to test whether within-individual network estimates predict functional response patterns and also to explore between-network functional dissociations.

7The raw data and our network parcellations generated through this research effort are provided to the community as an open resource.

8We sought to estimate networks within individuals with high precision. The analyses proceeded in three stages: 1) a refinement stage established the methods for estimating networks; 2) an implementation stage applied the methods prospectively to 15 new participants; and 3) a functional testing stage explored functional response properties and dissociations between networks.

9In the refinement stage, previously reported datasets (N = 2; see Refs., ) were analyzed to establish novel MS-HBM network estimates that incorporated priors for 15 distinct networks (in contrast to 10 networks used in earlier work). Each of the participants performed 31 independent MRI sessions, allowing considerable data to test for within-individual reliability.

10In the implementation stage, the 15-network MS-HBM was prospectively applied to 15 new participants who were each scanned over 8–11 sessions. The model was estimated for each participant in a fully automated fashion, and the networks were confirmed by model-free seed region-based functional connectivity. After network estimation, the overlap and variability of each network across individuals were examined and atlases constructed for open use.

11In the final functional testing stage, an extensive battery of tasks was administered to the same individuals and analyzed to explore whether the estimated networks predicted functional responses.

12Seventeen native English-speaking volunteers participated for payment. History of a neurological or psychiatric illness was an exclusion criterion. Participants provided written informed consent with protocols approved by the Institutional Review Board of Harvard University. For the refinement stage data, two right-handed adult women ages 22 and 23 yr participated [data previously reported in Braga et al. and Xue et al.]. The refinement stage data participants are labeled S1 and S2 to match Xue et al.. For the implementation stage data, 15 right-handed adults ages 18–34 yr participated (mean = 22.1 yr, SD = 3.9 yr; 9 women). Participants came from diverse racial and ethnic backgrounds (9 of the 17 individuals self-reported as non-White and/or Hispanic). A subset of the participants contributing implementation stage data also enrolled in a study of motor movement mapping. The implementation stage participants are labeled P1–P15.

13Data were acquired at the Harvard Center for Brain Science with a 3-T Prismafit MRI scanner (Siemens Healthineers, Erlangen, Germany). A 64-channel phased-array head-neck coil (Siemens Healthineers) was used in the refinement stage and for a subset of motor task sessions in the implementation stage. A 32-channel phased-array head coil (Siemens Healthineers) was used to acquire all other data in the implementation stage. For functional neuroimaging, the differences between these two coils are minimal, and the data were treated as comparable. Foam and inflated padding mitigated head motion. Participants were instructed to remain still and alert and to look at a rear-projected display through a mirror attached to the head coil. The display had a resolution of 1,280 × 1,024 pixels and a screen width of 43 cm, with an effective viewing distance of 104 cm (54 pixels per visual degree). Eyes were video recorded with an EyeLink 1000 Plus with Long-Range Mount (SR Research, Ottawa, ON, Canada), and alertness was scored during each functional run. MRI data quality was monitored during the scan with Framewise Integrated Real-time MRI Monitoring (FIRMM; see Ref. ).

14Each participant (S1 and S2) was scanned across 31 MRI sessions over 28–40 wk, with no sessions on consecutive days. Each session involved multiple resting-state fixation runs that were reanalyzed here for functional connectivity analysis, for a total of 63 functional MRI (fMRI) runs obtained for each individual. fMRI data were acquired with blood oxygenation level-dependent (BOLD) contrast. A custom multiband gradient-echo echo-planar pulse sequence provided by the Center for Magnetic Resonance Research (CMRR) at the University of Minnesota was used (Refs.,; see also Ref. ): voxel size = 2.4 mm, repetition time (TR) = 1,000 ms, echo time (TE) = 32.6 ms, flip angle = 64°, matrix 88 × 88 × 65, anterior-to-posterior (AP) phase encoding, multislice 5× acceleration, fully covering the cerebrum and cerebellum. Signal dropout was minimized by automatically selecting a slice 25° from the anterior-posterior commissural plane toward the coronal plane. Each run lasted 7 min 2 s (422 frames, with the first 12 frames removed for T1 equilibration). A dual-gradient-echo B0 field map was acquired to correct for spatial distortions: TE = 4.45 and 6.91 ms, with slice prescription/spatial resolution matched to the BOLD sequence. During BOLD scanning, participants fixated a centrally presented plus sign (black on a gray background). The scanner room was illuminated.

15A rapid T1-weighted (T1w) structural scan was obtained with a multiecho magnetization prepared rapid acquisition gradient echo (ME-MPRAGE) three-dimensional sequence: voxel size = 1.2 mm, TR = 2,200 ms, TE = 1.57, 3.39, 5.21, 7.03 ms, TI = 1,100 ms, flip angle = 7°, matrix 192 × 192 × 176, in-plane generalized autocalibrating partial parallel acquisition (GRAPPA) acceleration = 4.

16Each participant (P1−P15) was scanned across 8–11 sessions, most often over 6–10 wk. A few participants had longer gaps between the first and last MRI sessions up to 1 yr. Each session involved multiple fMRI runs to be used for functional connectivity analysis, for a total of 17–24 resting-state fixation runs obtained for each individual. BOLD acquisition parameters were similar to the refinement stage data: voxel size = 2.4 mm, TR = 1,000 ms, TE = 33.0 ms, flip angle = 64°, matrix 92 × 92 × 65 [field of view (FOV) = 221 × 221], 65 slices covering the full cerebrum and cerebellum. Each resting-state fixation run again lasted 7 min 2 s (422 frames, with the first 12 frames removed for T1 equilibration). Dual-gradient-echo B0 field maps were also acquired, with parameters matched to the refinement stage. The first two sessions of P12 were acquired in a different FOV (211 × 211); therefore, the matrix for both BOLD runs and field maps was 88 × 88 × 65 and BOLD TE = 32.6 ms, matching S1 and S2. The change in FOV did not affect the quality of registration or impact the analyses in any way we could detect.

17High-resolution T1w and T2-weighted (T2w) scans were acquired for the implementation stage data based on the Human Connectome Project (HCP; Ref. ): T1w MPRAGE parameters: voxel size = 0.8 mm, TR = 2,500 ms, TE = 1.81, 3.60, 5.39, and 7.18 ms, TI = 1,000 ms, flip angle = 8°, matrix 320 × 320 × 208, 144, in-plane GRAPPA acceleration = 2; T2w sampling perfection with application-optimized contrasts using different flip angle evolution sequence (SPACE) parameters: voxel size = 0.8 mm, TR = 3,200 ms, TE = 564 ms, 208 slices, matrix = 320 × 300 × 208, in-plane GRAPPA acceleration = 2. Rapid T1w structural scans were also obtained as backup with the refinement stage sequence but with matrix 192 × 192 × 144.

18To explore functional response properties, extensive task-based BOLD fMRI data were collected on participants P1–P15. Task runs used the same sequence as the resting-state fixation runs, ensuring that the estimated networks would be spatially aligned to the task-based data. Details of the task designs, stimuli, and run structure are described in Task Paradigms.

19Each BOLD fMRI run was examined for quality. Exclusion criteria generally consisted of the parameters reported in Xue et al., including 1) maximum absolute motion >1.8 mm and 2) slice-based SNR <130. Runs with SNR >100 but also SNR <130 were retained if motion and visual inspection indicated adequate quality. For the functional testing stage data, the maximum absolute motion for the Episodic Projection task was >2.5 mm given the long duration. One borderline motor run (P2) was included with motion of 1.9 mm as the motion was largely due to a linear drift. For the refinement stage, usable resting-state runs were 62 (S1) and 61 (S2) runs. For the implementation stage, usable resting-state runs ranged from 15 (P11) to 24 (P12) runs. For the functional testing stage, usable task runs ranged from 18 (P5) to 70 (P12) runs (see Table 1). All data exclusions were finalized before functional connectivity and task response analyses.

20Data were processed with an in-house preprocessing pipeline (“iProc”) that preserves spatial details by minimizing blurring and multiple interpolations [described in detail in Braga et al.]. For the refinement stage data (S1 and S2), the processed data were taken directly from Xue et al., with relevant method description repeated here. For the implementation stage data (P1–P15), the changes in processing included the use of the high-resolution T1w and T2w structural images. For one participant (P12), the registration failed with the 0.8-mm T1w structural image and their 1.2-mm image was used as a backup. For another participant (P1), only the 0.8-mm T1w structural image was used, without a paired T2w image.

21Data were interpolated to a 1-mm isotropic T1w native-space atlas (with all processing steps composed into a single interpolation) that was then projected with FreeSurfer v6.0.0 to the fsaverage6 cortical surface (40,962 vertices per hemisphere; see Ref. ). Four transformation matrices were calculated: 1) a motion correction matrix for each volume to the run’s middle volume [linear registration, 6 degrees of freedom (DOF); MCFLIRT, FSL], 2) a matrix for field map-unwarping the run’s middle volume, correcting for field inhomogeneities caused by susceptibility gradients (FUGUE, FSL), 3) a matrix for registering the field map-unwarped middle BOLD volume to the within-individual mean BOLD template (12 DOF; FLIRT, FSL), and 4) a matrix for registering the mean BOLD template to the participant’s T1w native-space image, which was resampled to 1.0-mm isotropic resolution (6 DOF; using boundary-based registration, FreeSurfer). The individual-specific mean BOLD template was created by averaging all field map-unwarped middle volumes after being registered to an upsampled 1.2 mm and unwarped midvolume template (an interim target, selected from a low-motion run, typically acquired close to a field map).

22For resting-state fixation runs, confounding variables including six head motion parameters, whole brain signal, ventricular signal, deep cerebral white matter signal, and their temporal derivatives were calculated from the BOLD data in the T1w native-space volume. The signals were regressed out from the BOLD data with 3dTproject (AFNI; Refs., ). The residual BOLD data were then band-pass filtered at 0.01–0.1-Hz using 3dBandpass (AFNI; Refs., ). For task data runs, only whole brain signal was regressed out (see Ref. ). Data registered to the T1w native-space atlas were resampled to the fsaverage6 standardized cortical surface mesh with trilinear interpolation (featuring 40,962 vertices per hemisphere; Ref. ) and then surface-smoothed with a 2-mm full-width-at-half-maximum (FWHM) Gaussian kernel. The iProc pipeline thus allowed for high-resolution and robustly aligned BOLD data, with minimal interpolation and signal loss, output to two relevant spaces: the native-space volume and the fsaverage6 cortical surface. Analyses were performed on the fsaverage6 cortical surface, but the network estimates (parcellations) were projected back into the individual participant’s native-space volume, allowing both surface-based and volume visualization.

23The MS-HBM was implemented to estimate cortical networks. The MS-HBM was independently implemented for the refinement stage data (S1 and S2) and then subsequently for the implementation stage data (in 3 separate groups: P1–P5, P6–P10, and P11–P15). Estimating the model separately for multiple small groups allowed for prospective replication. As the results reveal, the procedure was robust.

24First, the connectivity profile of each vertex on the fsaverage6 cortical surface was estimated as its functional connectivity to 1,175 regions of interest (ROIs) uniformly distributed across the fsaverage surface. For each run of data, the Pearson’s correlation coefficients between the fMRI time series at each vertex (40,962 vertices/hemisphere) and the 1,175 ROIs were computed. The resulting 40,962 × 1,175 correlation matrix was then binarized by keeping the top 10% of the correlations to obtain the functional connectivity profiles.

25Next, the expectation-maximization (EM) algorithm for estimating parameters in the MS-HBM was initialized with a group-level parcellation from a subset of the HCP S900 data release (that itself used the clustering algorithm from our previous study; Ref. ). It is important to note that the goal of applying the model in this study was to obtain the best estimate of networks within each individual participant’s dataset, not to train parameters and apply them to unseen data from new participants (see Ref. ). In this analysis, as with our other studies using this approach, we did not include the validation step described in Kong et al., so no spatial smoothness prior was applied. Only the training step described in Kong et al. was conducted here. A network label assignment for each vertex was obtained for each participant within the training step.

26Data from the two participants (S1 and S2) were analyzed together with the same MS-HBM. The data were used to estimate and compare the 15-network MS-HBM and the 10-network MS-HBM, as well as to explore the reliability of model estimates when subsets of data were analyzed. The results from these initial two participants guided the subsequent processing of the implementation stage data.

27The specific impetus for exploring a 15-network model was that networks at or near to the insula did not distinguish multiple networks that had been reported in the literature, variably labeled the cingulo-opercular network and the salience network (Ref.; see also Ref. ), as well as established distinctions at or around primary visual and somatomotor cortex.1 The 15 candidate networks explored here are labeled2 Somatomotor-A (SMOT-A), Somatomotor-B (SMOT-B), Premotor-Posterior Parietal Rostral (PM-PPr), Cingulo-Opercular (CG-OP), Salience/Parietal Memory Network (SAL/PMN), Dorsal Attention-A (dATN-A), Dorsal Attention-B (dATN-B), Frontoparietal Network-A (FPN-A),3 Frontoparietal Network-B (FPN-B), Default Network-A (DN-A), Default Network-B (DN-B), Language (LANG), Visual Central (VIS-C), Visual Peripheral (VIS-P), and Auditory (AUD).

28A key aspect of our methods is generalization and replication. The 15 participants in the implementation stage data were divided into discovery, replication, and triplication datasets of 5 participants each.4 The MS-HBM, initialized with a 15-network group-level parcellation obtained from a subset of the HCP S900 data (see Ref. ), was applied independently to the three separate datasets.

29When employing the MS-HBM, there are assumptions about the organization of the brain from the group prior, how many networks should be estimated, and assignment of vertices to only a single network. The idiosyncratic patterns of estimated networks thus could be distorted or fail to capture features of the underlying correlation matrix. To confirm that the individual network estimates were not obligated by the assumptions, a model-free seed region-based analysis was conducted using the same data as the MS-HBM, mirroring the procedures outlined by Braga and Buckner. The results were expected to converge if the model did not overly bias network assignments and to diverge if the assignments mismatched the underlying data patterns. Model-free seed region-based confirmation thus served as a check to ensure that network estimates properly captured individual correlation patterns.

30For this control check, the pairwise Pearson correlation coefficients between the fMRI time courses at each surface vertex were calculated for each resting-state fixation run, yielding an 81,924 × 81,924 matrix (40,962 vertices/hemisphere). The matrix was then Fisher r-to-z transformed and averaged across all runs to yield a single best estimate of the within-individual correlation matrix. This averaged matrix was used to explore network organization. The mean correlation maps were assigned to a cortical template combining left and right hemispheres of the fsaverage6 surface into the CIFTI format to interactively explore correlation maps with the Connectome Workbench’s wb_view software. Seed regions with robust functional connectivity correlation maps were manually selected within MS-HBM network boundaries. Anterior and posterior seed regions were recorded and visualized for each network in all the participants. Thresholds were set at z(r) > 0.2 for all seed regions. The color scales of correlation maps were set to range between 0.2 and 0.6, with the Jet lookup table (color bar) for visualization.

31Networks were first estimated and analyzed for each individual on the normalized fsaverage6 surface of FreeSurfer. Surface-based analyses allowed comparisons across individuals and utilization of the group-based priors for initialization of the MS-HBM. However, many applications require network assignments to be utilized within the native-space anatomy of the individual’s own volume (e.g., for presurgical planning and neuromodulation targeting). Given these needs, we devised a robust empirical procedure to project the network estimates back into each individual’s native-space T1w anatomical volume, as well as a procedure to verify that the transformation was valid.

32We constructed three separate images within the native-space volume that each varied from 0 to 255 in one of the three Cartesian x-, y-, and z-coordinate axes (e.g., the x-coordinate image possessed a volume that linearly varied in the x-dimension, going from 0 to 255 with no other variation across the image volume). Each separate axis volume was then projected to the fsaverage6 surface with mri_vol2surf and mri_surf2surf (FreeSurfer v6.0.0) with the same spatial transformation used for the projection of the participant’s BOLD fMRI data onto the fsaverage6 surface. Nearest-neighbor interpolation was used. The matrices for this projection were taken from each participant’s processing pipeline (iProc).

33In this manner, x, y, and z volume coordinates were obtained on the surface using the exact same spatial transformation matrix as originally applied to the BOLD data. We assigned each surface network label to its corresponding x, y, and z coordinates in the native-space volume. This resulted in a sparse 256 × 256 × 256 matrix in the volume, which was filled with nearest-neighbor interpolation (MATLAB knnsearch). We then masked this with the individual’s FSL-reoriented and binarized cortical ribbon generated by FreeSurfer during preprocessing. As a validity check, the final native-space network estimates were projected back to the surface and compared to the original MS-HBM surface estimates for each participant to ensure no spatial distortions.

34The parcellation results from MS-HBM were overlaid onto each individual’s T1w structural image. Sagittal, axial, and coronal slices were chosen to show common landmarks in each individual (midline, left and right insula, anterior commissure, primary sensory and motor cortices).

35Data acquired with BOLD-contrast (T2* images) and echo-planar imaging possess variable distortion and signal dropout due to magnetic susceptibility artifacts, especially near the sinus and ear canals (e.g., Ref. ). Vertex-based SNR maps were computed by taking the preprocessed time series from each resting-state fixation run (before regressing out confounding variables) and dividing the mean signal at each vertex by its standard deviation over time. The SNR maps were then averaged across the runs, resulting in an aggregate within-individual SNR map on the fsaverage6 surface. To visualize these effects in the native anatomy, surface maps were projected to the native-space volume with the procedure described above. The only difference is that linear interpolation (MATLAB scatteredInterpolant) was used to fill in the sparsely filled 256 × 256 × 256 matrix.

36To measure spatial variability across individuals, overlap maps of network assignments were computed. For each individual, the spatial extent of their estimated network was plotted simultaneously with all other participants and the percentage of overlap computed. In addition, the individual networks were plotted next to one another to appreciate the commonalities across individuals as well as the idiosyncratic features of each individual’s estimate (available in the Supplemental Materials).

37Overlap maps were also computed for the model-free seed region correlation maps. These maps make no assumption of a winner-take-all network assignment, so they provide a different view of network consistency or inconsistency across participants. For this final analysis, each individual’s seed region correlation map for each network was thresholded at z(r) > 0.2 and the overlap across participants plotted. The analysis was performed separately for both the anterior and posterior seed regions for each network.

38Finally, the estimates were used to construct consensus and agreement atlases intended to provide the community resources to integrate the present observations into their own processing pipelines.

39The human cerebral cortex is a complex structure with numerous sulci and gyri that can make it difficult to appreciate topographic patterns, including patterns that evolve over medial to lateral views and through complex structures like the insula. To appreciate global topographic relations, a flattened surface was created by editing the inflated surface file with the “TKSurfer” tool of FreeSurfer v6.0.0. The flattened surface was constructed to especially focus on the anatomy at and around the central sulcus extending into the insula. Five linear cuts were made on the midline of the inflated cortical surface, including one along the calcarine sulcus and four roughly equally spaced cuts radiating out from the medial wall. Next, a circular cut was made on the midline to allow the surface to unravel. Finally, the “mris_flatten” tool of FreeSurfer v6.0.0 was employed to create the flattened surface. This procedure was performed separately for the left and right hemispheres.

40After estimation of within-individual networks, functional response properties were explored in independent task-based data collected on the same individuals. The task paradigms were chosen based on literature review and our prior studies because of their ability to differentially activate distinct networks, and to do so robustly. A second feature of the selected task paradigms is that they were amenable to repeat testing either because extensive novel stimuli could be constructed (e.g., sentences, question probes) or because, by their nature, they were resilient to habituation even after many repetitions (e.g., flickering visual stimuli). Task details are described below.

41The motor task was extended from Buckner et al. to examine the organization of the foot, glute, hand, and tongue representations. Novel targeting of the glute representation allowed an intermediate body position to be mapped between the hand and foot [as reported earlier in Saadon-Grosman et al.]. The goal of this task paradigm was to activate somatotopic portions of SMOT-A and SMOT-B while minimizing nonspecific shared responses.

42After extensive prescan training, participants performed six types of active movements in the scanner: 1 and 2) left and right finger taps (thumb to index and thumb to middle); 3 and 4) left and right toes plantar flexion and dorsiflexion; 5) tongue movements from right to left (touching the premolar upper teeth); and 6) contraction and relaxation of their gluteal muscles. Each movement type was performed repeatedly across 10-s movement blocks. Before each movement block, a 2-s visual cue of a drawn body part with a text label informed the participant to initiate one of the six movement types. The fixation crosshair then changed to a slow-flickering black circle to pace the movements. The onset of the black circle cued movement of thumb to index finger, toe plantarflexion, tongue to the right, and glute contraction. The offset of the black circle cued movement of thumb to middle finger, toe dorsiflexion, tongue to the left, and glute relaxation. After five cycles, the word “END” instructed movement cessation. Twenty-four movement blocks (4 per movement type) occurred within each run, with 16-s blocks of passive fixation following each set of six movement blocks. Runs began and ended with fixation, yielding five fixation blocks per run.

43Each run lasted 7 min 8 s (428 frames, with the first 12 frames removed for T1 equilibration). Six motor runs were collected with full counterbalanced orders of movement conditions on each day. Runs were excluded from analysis if participants missed or failed to respond to cues, as confirmed by operators observing their alertness and movements from the control room.

44A visual retinotopic stimulation task was used to map visual cortex (similar to Refs., ). Our design had three levels of eccentricity stimuli (to map eccentricity gradients that span the V1, V2, V3 cluster) and separate vertical versus horizontal meridian stimuli (to map polar angle reversals that separate the borders of V1, V2, and V3; Ref., see also Ref. ). The goal of this task was to activate retinotopic portions of VIS-C and VIS-P.

45The basic stimulus consisted of a circular checkerboard that expanded outward from the central fixation point to approximate cortical expansion in visual cortex. Moving from center, the radius ring of the checkerboard became larger by a log step of 0.29. The resulting checkerboard was rendered out to 36 even rings cropped to a resolution of 1,024 × 1,024 pixels. To localize the meridians, two wedges masked the checkerboard, each covering 0.5° to 16.2° of eccentricity and 11.2° of polar angle. Horizontal wedges were centered at polar angles 360° and 180° and vertical wedges at 0° and 90°. To localize polar angle, the checkerboard was masked with a circular ring, which increased in size with increasing eccentricity. The center ring covered 0.5° to 1.6°, the middle ring 1.6° to 5.1°, and the peripheral ring 5.1° to 16.2°.

46Each run consisted of ten 10-s blocks of visual stimulation (2 blocks of each of the 5 conditions). The beginning, middle, and end of each run included a 20-s block of extended fixation. During stimulation, the checkerboard changed 6 times/s in the order white/black, color, black/white, color, white/black, black/white. The black center fixation dot unpredictably changed to gray (every 1–5 s). To ensure continuous fixation, participants pressed a button every time the fixation dot changed to gray. The primary contrasts of interest were horizontal versus vertical meridian blocks and, separately, the three eccentricity blocks versus each other.

47Each run lasted 4 min 30 s (270 frames, with the first 6 frames removed for T1 equilibration). Five runs were collected for each participant. Runs were excluded from analysis if participants missed trials and the eye video recordings indicated drowsiness. Lights within the scanner room were off during visual topography mapping, and a black occluding board was inserted into the scanner to prevent any light reflections.

48The oddball task explored detection of transient responses to salient, visual oddball targets that were uncommon relative to irrelevant nontargets and distracting nontargets [similar to Wynn et al.]. The goal of the task was to activate the SAL/PMN and CG-OP networks. Both networks have regions at or near the anterior insula and have been variably associated with response to task-relevant transients (see Refs.,, for discussion).

49Participants viewed a sequence of uppercase letters O and K in either black or red. Participants pressed a button with their right index finger when a red K appeared and withheld their responses to all the other letter-color combinations. The random trial ordering was set with Optseq. In each run, 10% of the trials were target red Ks, 10% were lure red Os, 40% were distractor black Ks, and 40% were distractor black Os. The contrast of interest was the target red Ks versus all other trials coded as the implicit baseline.

50Each run lasted 5 min 50 s (350 frames, with the first 6 frames removed for T1 equilibration). After 6 s of fixation overlapping the initial stabilization frames, an initial 20-s block of fixation was followed by a continuous extended block of 300 1-s trials (0.15-s presentation of the letter followed by 0.85 s of fixation) and then a final 20-s block of extended fixation. Before the first trial a 2-s start cue (1-s “Begin,” 1-s fixation) was presented, as well as a similar “End” cue after the final trial. Thus, the design was a rapid, event-related paradigm sandwiched between blocks of extended fixation. Five runs were collected for each participant. Runs were excluded from the analysis if participants missed more than six targets within a task run, which accounted for 20% of the total targets.

51The working memory (N-Back) task was extended from Cohen et al. and Braver et al. to explore demands on cognitive control under varied levels of memory load. Specifically, the N-Back task utilized a 2-back versus 0-back comparison to target FPN-A and FPN-B. In addition, following the design of the HCP N-Back task paradigm, multiple stimulus types/matching rules were included to explore whether the load effect was domain flexible or domain specialized.

52Stimuli were presented sequentially in the center of the computer screen. Participants decided whether the current stimulus matched a consistent template target (the 0-Back or low-load condition) or whether the current stimulus matched the stimulus presented two stimuli back in the sequence (the 2-Back or high-load condition). Participants maintained fixation on a central crosshair throughout the run.

53The stimuli varied across four conditions (Face, Word, Scene, and Letter) that were each presented in separate blocks. Faces and scenes were color images, with scenes showing both indoor and outdoor spaces and chosen not to feature people (faces from HCP, Ref.; scenes generously provided by the Konkle laboratory, Refs., ). Letters included subsets of consonants, and words featured one-syllable words from 10-word sets matched for length and frequency with the Corpus of Contemporary English (Ref.; December 2015 version). In all but the Word condition, participants matched the stimuli to an exact stimulus referent, or the exact stimulus presented two trials before. For the Word condition, the participants decided whether the current word rhymed with the target (e.g., “dream” would be a positive match with “steam”).

54Each N-Back run featured eight blocks (a 0-Back and a 2-Back for each of the 4 stimulus categories). Each block included a cue and nine trials. During the first cue stimulus, participants also saw the block type, either 2-Back or 0-Back. During 2-Back blocks participants looked for matches (identical images or rhyming words) with the stimulus two trials back, and during 0-Back blocks participants looked for matches to the cue. The background was black (matching the HCP format). All blocks included two target and two lure (repeated nontarget) trials. Targets and lures were equally likely to appear in each viable trial position within and across runs. Participants pressed a button for every trial, indicating match (right index finger) or no match (left index finger).

55Each run lasted 4 min 44 s (284 frames, with the first 12 frames removed for T1 equilibration). After 12 s of fixation overlapping the initial stabilization frames, an additional block of 12 s of fixation was followed by blocks of the N-Back task interspersed with 15-s fixation blocks (the fixation blocks came after two 25-s N-Back task blocks). Across runs, 0-Back and 2-Back blocks, categories, and their interactions were counterbalanced. Each trial was 2.5 s in duration (2 s of stimulus presentation followed by 0.5 s of fixation). The fixation crosshair was white for the extended fixation blocks and green during the N-Back task blocks. Within a run, all categories were seen before a category repeated. Eight runs were collected for each participant. Runs where participants missed responses in more than two trials were excluded from analysis.

56The Sentence Processing task was adapted from Fedorenko et al. to examine domain-specialized processing related to accessing word meaning and phrase-level meaning. The target task involved sentences presented one word at a time. The reference control task was presentation of nonword strings that were matched in length and visually similar. The goal of this task was to activate the LANG network [see Braga et al.].

57Participants passively read real sentences (“IN THE MORNING THE TAILOR WAS SHOWING DIFFERENT FABRICS TO THE CUSTOMER”) or pronounceable nonword strings (“SMOLE MUFRISONA VEDER SMOP FO BON FE PAME OMOSTREME GURY U FO”). The centered stimuli were presented one word (or nonword) at a time (0.45 s per word). After each word or nonword string, a cue appeared for 0.50 s, prompting the participants to make a right index finger button press. Stimuli, generously provided by the Fedorenko laboratory, never repeated. Word or nonword strings (6 s each) were presented in 18-s blocks of three strings. Extended fixation blocks (18 s each) appeared at the start of each run and after every fourth string block. The primary comparison of interest was the contrast between sentence and nonword blocks.

58Each run lasted 5 min 0 s (300 frames, with the first 12 frames removed for T1 equilibration). Six runs were collected for each participant. Runs were excluded if the participant did not read the stimuli (observed through eye monitoring) or missed responses.

59The Theory-of-Mind tasks were adopted from Saxe and colleagues to explore domain-specialized processing associated with representation of another person’s mental states. In the False Belief paradigm, participants viewed a brief story and then, on a separate screen, a question about that story. In the False Belief condition, the target stories described events surrounding a person’s perspective, followed by a question about the thoughts and beliefs of that person. In the control False Photo condition, stories described similar situations involving objects (e.g., in photos, on maps, and in descriptions). In the Emotional/Physical Pain Stories paradigm (subsequently abbreviated “Pain”), the target stories described a situation that evoked personal emotional pain (Emo Pain condition) and were contrasted with control stories of similar length and complexity involving physical pain (Phys Pain condition). Participants rated the level of pain from “None” to “A lot” during the question period. These two paradigms yield similar task activation maps. Here the task contrasts of False Belief versus False Photo and Emo Pain versus Phys Pain were combined, with the goal of activating DN-B (extending from Ref. ).

60Each run consisted of a series of stories and questions (15 s per individual story/question pairing). For both paradigms, each run included five target trials (False Belief or Emo Pain) and five control trials (False Photo or Phys Pain). Fifteen-second fixation periods occurred between trials. Stimuli never repeated.

61Each run lasted 5 min 18 s (318 frames, with the first 12 frames removed for T1 equilibration). Eight runs were collected for each participant: four of the False Belief paradigm and four of the Pain paradigm. We implemented an exclusion criterion to exclude any run with more than one missed trial. No runs met this criterion.

62The Episodic Projection task was adapted from Andrews-Hanna et al. and DiNicola et al. to encourage processes related to remembering the past and imagining the future (prospection). In the target task conditions, participants viewed a brief scenario that oriented to a situation in the past (Past Self) or future (Future Self) simultaneously with a question about the event that encouraged participants to imagine the specific scenario described. The similarly structured control condition asked the participants about a present situation (Present Self). The task contrasts of Past Self versus Present Self and Future Self versus Present Self were combined, with the goal of activating DN-A [extending from DiNicola et al.]. Of relevance, detailed behavioral analyses of these contrasts have suggested that the main component process tracking increased response in DN-A is the process of mentally constructing scenes (Refs.,; see also Ref. ). Thus, the task contrast used here taps into domain-specialized processing related to spatial/scene processing (see Ref. for discussion).

63Each run contained a series of scenarios with questions (10 s of scenario/question presentation followed by 10 s of fixation). Thirty questions appeared per run, with three for each condition of relevance (Past Self, Future Self, Present Self). Additional conditions were included toward goals distinct from those targeted here. For our analyses, we focus on the condition contrasts that have previously dissociated DN-A from DN-B in DiNicola et al.. All scenarios were unique.

64Each run lasted 10 min 17 s (617 frames, with the first 12 frames removed for T1 equilibration). Ten runs were collected for each participant that included 90 relevant trials across runs (30 of each of the 3 conditions). Runs with more than two missed trials were excluded.

65Functional task data were analyzed with the general linear model (GLM) as implemented by FSL’s first-level FEAT (FSL version 5.0.4; Ref. ). All conditions were included in each model design, even those not relevant to the contrasts of interest, except for the Oddball Effect task contrast, which coded the targets against the implicit baseline. The data were high-pass filtered with a cutoff of 100 s (0.01 Hz) to remove low-frequency noise within each run. GLM outputs included, for each contrast, β values for each vertex that were converted, within FEAT, to z values. Within each participant, z-value maps from all runs were averaged together by using fslmaths to create a single cross-session map for each contrast of interest. For the N-Back task, we ran both block-level and condition-level GLMs. Block-level GLM outputs included z-value maps for each trial block, which were averaged by condition across runs. Condition-level GLM outputs included z-value maps for 2-Back and 0-Back conditions, which were averaged across runs. In both cases, a single cross-session contrast map was then created by taking the difference between condition mean maps.

66Task contrasts were designed to functionally target specific networks and dissociate response properties between networks. Two convergent methods were used for visualization and quantification. First, z-value maps were compared visually by overlaying the borders of networks onto the task contrast maps on the same cortical surface (fsaverage6 cortical surface). This form of visualization allowed comprehensive assessment of task response patterns. Contrast z-value maps were manually thresholded to best demonstrate the task activation patterns for each participant. The PSYCH-FIXED lookup table within Connectome Workbench was used for the color scale.

67Second, a priori networks within individuals were used to formally quantify differences in response levels between networks, including direct tests for significant differences between networks and between task contrasts. For each task contrast, the average z value was calculated for all vertices within each selected network, combining across hemispheres. Mean z values were computed for each task run, and the cross-run mean z values for each network were then plotted in a bar graph, along with the standard error of the mean across participants. This analysis has the advantage of quantifying the magnitude and variance of the response in each a priori-defined network for each participant, without any subjective decisions.

68For both approaches to task response analysis, the networks were defined within the individuals before examination of the task maps, to avoid the possibility of bias.

69Functional connectivity between brain regions was calculated in MATLAB (version 2019a; http://www.mathworks.com; MathWorks, Natick, MA) using Pearson’s product moment correlations. FreeSurfer v6.0.0, FSL, and AFNI were used during data processing. The estimates of networks in volume space were visualized in FreeView v6.0.0. The estimates of networks on the cortical space were visualized in Connectome Workbench v1.3.2. Statistical analyses were performed with R v3.6.2. Model-free seed region confirmations were performed in Connectome Workbench v1.3.2. Network parcellation was performed using code from Kong et al. on GitHub (https://github.com/ThomasYeoLab/CBIG/tree/master/stable_projects/brain_parcellation/Kong2019_MSHBM).

70Networks were estimated for the refinement stage data with a 15-network MS-HBM model. Figure 1 and Figs. A1–A5 in appendix display the main results for S1 and S2 on the surface, and the Supplemental Materials display the comprehensive results and quality control visualizations on the surface and in the native-space volume.

71The first results pertain to data quality. The SNR maps are displayed on the cortical surface (Fig. 1 and Fig. A3). Most of the cortical mantle possessed high SNR. As expected, given signal dropout near the sinuses and the inner ear, there was variability in SNR across the cortical surface with orbitofrontal cortex (OFC) and adjacent ventrolateral prefrontal cortex (VLPFC), rostral inferior temporal cortex, and the temporal pole showing low SNR (see Supplemental Materials for additional visualizations). Network assignments in low-SNR regions should be interpreted cautiously.

72The primary result of our procedures was an estimated parcellation into distinct candidate networks. Figs. A1 and A4 display the 15-network estimates for S1 and S2. All networks, including local sensory and motor networks as well as distributed association networks, were identified in both participants. Although the general organization was shared between the two participants, the spatial boundaries were idiosyncratic. These patterns are elaborated upon in detail in the results of the novel 15 participants. For these first two individuals, we focused on validating the methods.

73The network estimates were based on a 15-network MS-HBM. In addition to assuming a specific number of networks, the method also employed group priors to constrain the estimates (see Supplemental Fig. S1). As such, it is possible that the resultant networks do not accurately reflect the underlying within-individual correlation patterns as one might expect. To explore this possibility and intuitively visualize the degree to which the model captures underlying correlation patterns, a model-free seed region-based correlation analysis was performed. A seed region was placed in an anterior position and separately a posterior position within each network within each individual. The resulting correlation maps are displayed in Figs. A2 and A5 in relation to the MS-HBM network boundaries.

74The estimated networks captured features of the correlation patterns remarkably well, including across small, distributed regions that might otherwise be overlooked. The alignments were not perfect. Specifically, the correlation patterns included most of the distributed regions in the MS-HBM solutions, and the patterns were largely selective to the estimated networks. Small deviations, in the form of extensions of the patterns beyond the network boundaries, were common, likely in part because the network estimates forced a winner-take-all assignment but also possibly because additional network details may be missed.5 The consistency between the general correlational structure and the network estimates in one sense is unsurprising because the underlying correlation matrix was employed by the network model. However, it is not obligated, and deviations could be seen if the model forced assignments or the model failed to capture the structure of the data.

75We next sought to explore what is gained by adopting the 15-network parcellation rather than the simpler 10-network parcellation. Figures A6 and A9 display the MS-HBM parcellation estimate for the 10-network and 15-network solutions for each participant. The first notable result is that, for most networks, there was little difference between the two models’ estimates. For example, the separation of DN-A and DN-B was well captured by both model solutions, with the distributed spatial patterns and idiosyncratic features quite similar between models. That is, if the goal were to study DN-A and DN-B, there is little gained by utilizing the more complex 15-network model. In both S1 and S2, many of the other major networks were also similar between the two parcellations, including FPN-A, FPN-B, SMOT-A, and SMOT-B. Thus, for networks well captured by the 10-network model, they appear to be roughly unchanged in the 15-network model. For other networks though, there were substantive differences.

76One motivation for investigating a 15-network model was that certain networks did not differentiate established distinctions at or around somatomotor cortex and visual cortex, as well as between multiple networks within or adjacent to the insula including separation of a cingulo-opercular network from a salience network (Ref.; see also Ref. for discussion). These features were captured in the 15-network MS-HBM. Specifically, the single visual network in the 10-network estimate was differentiated among dATN-B, VIS-C, and VIS-P in the 15-network solution (Fig. A7). The Salience (SAL) network in the 10-network estimate was differentiated into two separate networks, here labeled SAL/PMN and CG-OP (Fig. A11). The dATN in the 10-network estimate was differentiated into dATN-A and PM-PPr in the 15-network solution (Fig. A8), and a distinct AUD network emerged near LANG and SMOT-B (Fig. A10). Critically, seed region-based correlation patterns suggested that this expansion of networks from 10 to 15 captured clear features of the underlying correlation patterns (Figs. A7, A8, A10, and A11).

77One unexpected result was that our 15-network parcellation included a single network that has been variably described in the literature. What has been called the “parietal memory network”, with focus on the posterior midline, has often been discussed separately from the network referred to by Seeley and colleagues as the “salience network”. Here a single distributed network was identified that possessed the canonical features of both networks. The seed region-based correlation maps supported that the two networks discussed historically as distinct are likely a single network (Figs. A2, A5, and A11), a result that is further examined in the prospectively acquired and analyzed data.

78We next sought to address two related questions. First, are the network estimates described above reliable within individuals? Second, can they be obtained with a lesser amount of data? The resting-state fixation runs of S1 and S2 were divided into three datasets with roughly equal amounts of runs contributing to each data subset (20/20/22 runs of data for S1 and 20/20/21 runs of data for S2). The 15-network MS-HBM was estimated independently for each data subset. Results are displayed in Fig. 2.

79In S1, 82.7% on average of cortical vertices were assigned to the same networks across the independent datasets from within the individual. In S2, 86.9% of cortical vertices were assigned to the same networks. By contrast, overlap between the separate parcellations of S1 and S2 were 54.4%, 55.1%, and 55.4%, indicating that between-individual variability was substantially larger than within-individual variability.

80The reliability of network estimates derived from varying the amount of data, ranging from 2 to 20 runs of resting-state fixation data, is visually illustrated in Fig. 3 and Fig. 4. In each dataset, we independently applied the 15-network MS-HBM model. As shown in the figures, a few runs of data can capture the overall correlational patterns but fall short of revealing the idiosyncratic details of networks within the individual. As the quantity of data integrated into the model increased, both S1 and S2 demonstrated network estimates that displayed a noticeable trend toward greater consistency between independent datasets. This improvement is quantified in the Supplemental Materials. In general, robust networks emerged when we leveraged ∼20 runs of data for both S1 and S2. These findings suggest that cortical parcellations of the resolution and within-individual detail targeted here are replicable for models based on ∼20 runs of data. Notably, this is the amount of data collected for the 15 new participants in the implementation stage dataset analyzed throughout the rest of this article.

81Fifteen cerebral networks were estimated for all new participants. The 15 individuals were analyzed within subsamples (each n = 5) intended to replicate the MS-HBM’s network estimates in prospective participants, including novel discovery (P1–P5), replication (P6–P10), and triplication (P11–P15) datasets. Results were similar across all three subsamples, and the full parcellation for each individual is available in the Supplemental Materials on the surface and within the individual’s own native-space volume. Despite idiosyncratic spatial details of network organization, the broad properties were largely consistent. Three representative participants, one from each subsample, are displayed in Fig. 5, Fig. 6, and Fig. 7.

82VIS-C and VIS-P were identified in each participant, extending from the calcarine sulcus on the midline to the lateral surface. The extent of the two networks across the occipital lobe did not align them with individual visual areas but rather the expanded regions of the V1/V2/V3 retinotopic cluster, and likely adjacent retinotopic clusters. The multiple networks appeared to divide along the eccentricity gradient. The VIS-C network overlapped regions likely aligned to the central portions of the V1/V2/V3 retinotopic representations, whereas VIS-P overlapped the peripheral retinotopic representations (see Ref. ). The relation of VIS-C and VIS-P to task-elicited responses is directly explored in Somatomotor and Visual Networks Respond to Body Movements and Visual Stimulation in a Topographic Manner.

83Although the VIS-C and VIS-P networks contained vertices that were mostly contiguous, there were exceptions. Discontinuous islands were sometimes found in occipital-temporal cortex, possibly a reflection of separate extrastriate retinotopic clusters (e.g., at or near the MT/V5 hemifield representation). VIS-P also occasionally contained small, punctate representations near dorsolateral prefrontal cortex (DLPFC). These were the exceptions: the majority of the VIS-C and VIS-P network-included vertices were continuous and adjacent to one another, overlapping the expected location of early retinotopic visual cortex.

84Similarly, SMOT-A and SMOT-B were identified reliably as spatially continuous networks along the central sulcus, extending onto the midline and into the posterior insula. These two somatomotor networks also do not likely align to individual architectonic areas but rather extend across the pre- and postcentral gyrus, including primary motor as well as somatosensory areas. The extent along the midline and into the posterior insula further suggests that the networks span multiple body maps, not simply the dominant inverted body map along the central sulcus.

85A final predominantly local sensory network, AUD, was consistently identified near the superior temporal sulcus. This network extended across the full supratemporal plane including Heschl’s gyrus, and into adjacent regions.

86Multiple distributed networks were identified in each participant that were immediately adjacent to the local sensory and motor networks, with each network containing distributed regions that spanned multiple zones of cortex. dATN-A and dATN-B were adjacent to VIS-C and VIS-P but also with distant regions in the frontal cortex, likely at or near the frontal eye field (FEF; Refs., ). Similarly, CG-OP and PM-PPr radiated outward from the early somatomotor networks SMOT-A and SMOT-B. CG-OP and PM-PPr sometimes contained small islands indenting or even within the SMOT network boundaries that may relate to interspersed intereffector regions along the central sulcus. CG-OP and PM-PPr also included regions abutting and within the Sylvian fissure. The relations among the networks became even clearer in the flat map visualizations.

87A final network, SAL/PMN, displayed a spatial pattern that was adjacent to CG-OP in many locations but also with differences. Although SAL/PMN contained a prominent region in the anterior insula, the network’s positioning did not juxtapose the somatomotor networks. Rather, SAL/PMN was adjacent to a posterior midline cluster of association networks near regions of the canonical “default network” (e.g., Refs.,,, ). SAL/PMN consistently included a region within ACC anterior to CG-OP and a prominent set of regions along the midcingulate and the posterior midline. As noted for S1 and S2, the SAL/PMN network’s spatial pattern combined features described in prior work on the salience network [Ref.; see Seeley for discussion and Dosenbach et al. for related work] and the parietal memory network.

88The remaining regions of association cortex, which contain the majority of PFC, a large region of PPC extending into the temporoparietal junction (TPJ), and lateral temporal cortex (LTC), were populated by five distinct networks. With some exceptions, each of these five networks tended to possess regions in each of the distributed zones. The five networks were interwoven with local patterns of adjacencies that repeated across the cortex.

89Specifically, FPN-A and FPN-B were adjacent to one another throughout the cortical mantle. FPN-A and FPN-B displayed a distributed pattern consistent with the well-studied group-estimated frontoparietal control network, also referred to as the multiple-demand system. These two juxtaposed networks (FPN-A and FPN-B) consistently neighbored an additional clustered set of three networks: LANG, DN-B, and DN-A. These three additional networks were tightly juxtaposed among themselves on the lateral cortical surface including zones within PPC, LTC, and both DLPFC and VLPFC. DN-A and DN-B were interdigitated as well along the anterior and posterior midline, consistent with previous studies.

90Despite being adjacent in many locations, clear features distinguished the three networks. The LANG network surrounded the Sylvian fissure and included regions near the AUD network and in VLPFC at or near historically defined “Broca’s area.” The LANG network was generally larger in the left compared with the right hemisphere [see Braga et al.]. DN-A showed a strong correlation with the posterior parahippocampal cortex (PHC) [see also Reznik et al. for further details]. Additionally, DN-A occupied regions at or adjacent to the retrosplenial cortex (RSC) and ventral posterior cingulate cortex (PCC). DN-B prominently included anterior regions of the inferior parietal lobule extending into the TPJ [a region of particular focus, e.g., Saxe and Kanwisher, Jacoby et al.]. The posterior midline region of DN-B fell between regions of DN-A and specifically did not extend into RSC. DN-B also included a larger region of the LTC than DN-A; DN-A tended to include a small region or a few discontinuous regions in anterior LTC.

91Of importance, the spatial arrangements of the five networks (FPN-A, FPN-B, LANG, DN-B, and DN-A) repeated multiple times across the cortical mantle, a discovery that is expanded upon in the analyses of spatial juxtapositions on the flattened cortical surface.

92To enable clinical applications, particularly for presurgical planning and neuromodulation, an empirical method was developed to project the networks estimated on the surface into the native-space volume of individual participants. Figure 8, Figure 9, and Figure 10 display examples for three representative participants from the discovery (P1), replication (P6), and triplication (P11) datasets. The network estimates in the native-space volume were overlaid onto each individual’s T1w structural image and displayed in sagittal, coronal, and axial views. Similar maps for all available participants can be found in the Supplemental Materials.

93Several details become apparent that were lost in the surface visualizations, including that, within the native-space volume, regions that appear distant on the surface are often quite close in the volume. Furthermore, network assignments in one cortical gyrus were often similar to those in an adjacent gyrus, even those in separate lobes (and thus presumably quite distant from one another on the surface). For example, temporal pole regions, visualized particularly well in the sagittal plane, were often assigned to networks also present in the ventral prefrontal regions abutting the Sylvian fissure (e.g., regions linked to the DN-B and LANG networks). It remains an open question whether these between-gyri juxtapositions are a residual artifact of spatial blurring or an interesting feature of anatomy that minimizes absolute distances between functionally similar cortical regions.

94Certain aspects of the network estimates were impacted by signal loss. Low-SNR regions were observed in the OFC, ventral regions of VLPFC, and anterior regions of the temporal lobe (see Fig. 1, Fig. A3, and Supplemental Materials). When interpreting the network assignments, it is important to keep these spatially variable effects in mind. For example, a localized AUD network was observed across the supratemporal plane including Heschl’s gyrus. Inconsistent, discontinuous vertices were also labeled as part of the AUD network in the inferior temporal cortex and OFC, in the regions of the greatest signal dropout due to magnetic susceptibility differences. The network assignments in low-SNR regions should not be trusted.

95To demonstrate that the correlation properties of the within-individual data were captured in the network assignments, seed region-based correlation maps were visualized. Although there were minor differences between the MS-HBM network estimates and the seed region-based correlation maps,5 networks could be identified in all participants with both methods. Furthermore, the maps defined by anterior and posterior seed regions were similar, indicating that the seed region-based method was not dependent on a single vertex or one general region of cortex. Seed region-based confirmation for DN-A, DN-B, LANG, FPN-A, FPN-B, CG-OP, and SAL/PMN is displayed for three representative participants in Fig. 11, Fig. 12, and Fig. 13 and for all participants and additional networks in the Supplemental Materials.

96An overlap map of assignments for each network from the MS-HBM for the 15 participants is displayed in Fig. 14. Results revealed that the general organization of the networks was highly conserved across individuals, but with differences in the idiosyncratic spatial positioning and extents of the networks.

97A challenge in examining spatial overlap is that there is circularity in network definition since the process initiates with the same 15-network group prior, which could bias the networks to show more overlap than is truly in the data. To mitigate this concern, we also examined overlap using the network estimates derived from the seed region correlation maps. These maps are not constrained by the group prior and do not enforce a winner-take-all assumption, allowing deviations to emerge. Overlap maps of correlation patterns were obtained using anterior and posterior seed regions within each network for all 15 participants (Fig. A12).

98As another exploration of variability, the individual networks were plotted next to one another for all 15 participants, allowing another means to identify shared and idiosyncratic features of each participant’s estimate. The results are available in the Supplemental Materials.

99To provide atlases as open resources to the community for future use, the overlap of network estimates across the 15 participants was used to create 1) a network consensus atlas and 2) a series of thresholded network agreement atlases. These atlases can be used prospectively by researchers to construct regions that have high likelihood of being in one network or another without prior individual-level data.

100In the network “consensus” atlas, each vertex was assigned to the network that was most probable. For example, if a vertex was labeled DN-A in seven participants, DN-B in five participants, and LANG in three participants, it would be assigned to DN-A in the consensus atlas. Every vertex was assigned the label of its most likely network, including uncertain vertices. Supplemental Fig. S18 displays the consensus atlas.

101In the network “agreement” atlases, each vertex was assigned to its most common network but only if a certain number of participants showed agreement. For example, in the n ≥ 8 (53%) agreement atlas, displayed in Fig. 15, only vertices where eight or more participants agreed on the same network assignment were labeled. Higher-order networks exhibit a lower level of agreement among participants in the prefrontal, parietal, and temporal association cortices. Supplemental Fig. S12 displays agreement atlases for various count thresholds ranging from n ≥ 6 (40%) to n ≥ 12 (80%). These agreement atlases highlight cortical locations where many or most individuals would likely be assigned to the same network.

102To better reveal spatial relations among networks, a flattened cortical surface was constructed (Fig. 16). The 15 networks are displayed in representative participants in Fig. 17, Fig. 18, and Fig. 19 and for all participants in the Supplemental Materials. A broad observation was that higher-order networks nest outward from sensory and motor cortices.

103Specifically, the networks could be heuristically grouped into three levels beginning with first-order sensory and motor networks.6 The first-order networks were primarily locally organized, spatially arranged along the central sulcus (for SMOT-A and SMOT-B) and near the calcarine sulcus (for VIS-C and VIS-P). Surrounding these first-order networks were adjacent networks that radiated outward. We refer to these as second-order networks. CG-OP and PM-PPr surrounded SMOT-A and SMOT-B, and dATN-A and dATN-B were adjacent to VIS-C and VIS-P. In between these second-order networks were third-order networks (FPN-A, FPN-B, LANG, DN-B, and DN-A) that populated the large, expanded zones of higher-order association cortex. The flattened representation allowed further features to be appreciated.

104CG-OP and PM-PPr nearly fully surrounded both the anterior and posterior extents of the somatomotor networks, including the insular regions that are buried within the Sylvian fissure. Although CG-OP and PM-PPr were generally interdigitated around the first-order somatomotor networks, in several locations CG-OP fell distal to PM-PPr (meaning that PM-PPr directly juxtaposed SMOT-A and SMOT-B and CG-OP juxtaposed PM-PPr). Furthermore, whereas the PM-PPr network was adjacent to the somatomotor networks across its extent, CG-OP also involved distant regions in PFC and posterior association zones that were not adjacent to somatomotor networks. Thus, CG-OP possessed a partially distributed motif. Additional details of dATN-A and dATN-B were also evident. Of the two networks, dATN-B fell more proximal to the early visual networks and dATN-A more distal. dATN-A included distributed regions in frontal cortex at or near FEF, whereas dATN-B was more locally organized but not exclusively so.

105Comparing dATN-A and dATN-B with CG-OP and PM-PPr, as highlighted in panels A and B of Fig. 17, Fig. 18, and Fig. 19, revealed parallel features. The second-order networks were all spatially anchored near the early (first order) sensory and motor networks, appearing as if they grew out or formed from the earlier networks. CG-OP in several individuals extended into the precentral gyrus.5 And, despite being far apart in their major extents, both sets of networks had distributed components throughout association cortex, including adjacencies in frontal cortex. Thus, from the standpoint of a potential hierarchy among networks, these second-order networks possess a motif that anchors them to the early sensory and motor networks and simultaneously connects them to distributed zones of association cortex.

106We provisionally label SAL/PMN as a second-order network, but it has juxtapositions that differentiate it from the other second-order networks.6 Across much of its extent, SAL/PMN paralleled CG-OP with multiple juxtapositions. SAL/PMN differed in that it was not adjacent to early sensory and motor networks. Rather, SAL/PMN contained regions that were near the network labeled DN-A, especially along the posterior midline, where its regions could easily be confused with the large DN-A and DN-B regions that occupied much of the posterior midline.

107The zones farthest away from the sensory and motor regions were populated by five third-order association networks (FPN-A, FPN-B, LANG, DN-B, and DN-A). Each third-order network possessed regions distributed widely throughout association cortex. Moreover, regions of distinct third-order networks displayed side-by-side juxtapositions with a pattern that repeated similarly across multiple zones of cortex. We focus on these repeating clusters of five networks extensively below.

108The spatial extent of task-elicited responses to body movements and visual stimulation was mapped in direct relation to the network boundaries. The goal of these analyses was to explore whether within-individual network estimates predict task responses. We start here with descriptions of sensory and motor responses and explore responses to varied cognitive task demands below. In all cases, the network boundaries were established before examination of the task responses. Fig. A13 illustrates the body movement and visual stimulation mapping strategy, and Fig. 20 displays the detailed maps for one representative participant on the inflated and flattened surfaces. Additional participants are displayed in Fig. 21, and all participants with available task data are shown in the Supplemental Materials.

109Task activation maps revealed that body movement and visual stimulation task contrasts elicited responses that were aligned to, and generally filled in, the first-order network estimates (SMOT-A, SMOT-B, VIS-C, and VIS-P). That is, the idiosyncratic network estimates in each individual predicted the localization of the movement and visual stimulation responses. The visual responses often extended beyond the anterior boundaries of VIS-C and VIS-P, including portions of dATN-B, but generally followed the network boundaries. Several additional results are notable.

110First, the main body map along the central sulcus extended across networks (SMOT-A and SMOT-B), as did the retinotopic eccentricity gradient (VIS-C and VIS-P). Within the visual system, there was a clear correspondence between the two visual networks and eccentricity. VIS-C, as anticipated given its anatomical position, tracked the central representation. VIS-P covered the peripheral representation. A gap emerged for the most peripheral regions of the dorsal and ventral eccentricity portions of VIS-P, possibly because the visual stimulation did not extend fully to the periphery [see Park et al.]. Within the somatosensory and motor systems, there was a distinct gap between the representations of the hand and glutes, which may be an intereffector region, as discussed more below.

111Second, the response patterns did not align to expected boundaries of individual brain areas (i.e., V1, S1). The body movements activated regions of pre- and postcentral gyrus, spanning multiple motor and somatosensory areas. Examined in detail, the body movement responses suggest at least three distinct maps of the main body axis, labeled I, II, and III in Fig. 20. The largest distinct body map was found aligned to primary somatomotor cortex, exhibiting a medial-to-lateral progression from foot to hand to tongue (Fig. 20, labeled I). In the posterior insula, the body map was buried with a posterior to anterior orientation (Fig. 20, labeled II). On the medial wall, the body map progressed from anterior to posterior (Fig. 20, labeled III). Similarly, the visual responses spanned the extent of at least the V1/V2/V3 retinotopic cluster, with networks cutting across areas (verified through polar mapping as illustrated in Fig. A13E). Thus, the response patterns confirm that early somatomotor and visual networks group multiple areas together and split areas along topographic gradients (e.g., VIS-C vs. VIS-P; see Ref. ).

112Third, the gap in the body map between the representations of the hand and glutes was associated with a distinct distributed correlation pattern that created discontinuities within the primary motor body map (including breaks between effector regions along the precentral gyrus). Replicating Gordon et al., Fig. 22 illustrates the “intereffector” pattern. In almost every individual, placing a seed region within the gap revealed a robust distributed correlation pattern that typically contained at least two additional zones along the precentral gyrus: a region ventral to the face representation and a dorsomedial region near the foot representation. The pattern substantially overlapped the CG-OP network including regions surrounding the primary somatomotor networks and extending into the CG-OP regions distributed throughout the cortex. The intereffector seed region pattern also included prominent regions at and around the primary visual retinotopic cluster, most typically at the edges of the peripheral representation, consistent with the integration of the visual domain in addition to the motor domains noted by Gordon et al..

113The oddball task was designed to measure the transient response to uncommon visually salient targets that require participant response. The mapping strategy is illustrated in Fig. A14. On the flattened cortical surface, the within-individual a priori-defined networks CG-OP and SAL/PMN are displayed in relation to the Oddball Effect task contrast. The details of the Oddball Effect task contrast are shown for one representative participant in Fig. 23. Figure 24 illustrates that the features can be observed in additional participants, and in all participants as shown in the Supplemental Materials.

114The Oddball Effect task contrast response was widely distributed across the cortex. The response prominently involved the distributed regions of the CG-OP and SAL/PMN networks, including regions in the anterior insula as well as along the posterior midline. These collective regions have been the emphasis of prior studies separately focused on the salience network and the parietal memory network. Thus, as predicted by the hypothesis that SAL/PMN is a single network, the response pattern observed here extended across the full distributed extent of the network.

115In addition to the consistent responses across the distributed regions of CG-OP and SAL/PMN, additional responses were reliably observed: a response along the central sulcus in the left hemisphere near the estimated location of the hand representation and along the calcarine sulcus near the central representation of the visual field (contrast Fig. 23 with Fig. 20). The response in the hand region of the central sulcus was exclusively in the left hemisphere, consistent with the right-handed response.

116To quantify the selectivity of the task response, the mean z values for the Oddball Effect task contrast were calculated separately for each association network. The estimates were obtained within the bounds of each individual’s a priori-defined networks and then averaged across participants (N = 14). Results plotted in Fig. 25 revealed a strong, significant positive response to oddball targets in both the CG-OP [t(13) = 7.97, P < 0.001] and SAL/PMN [t(13) = 6.21, P < 0.001] networks. By contrast, for many of the third-order association networks the response was significantly negative [DN-A: t(13) = −11.76, P < 0.001, DN-B: t(13) = −8.81, P < 0.001, LANG: t(13) = −3.82, P < 0.01, FPN-B: t(13) = −3.02, P < 0.01], with FPN-A being the exception. FPN-A showed a weak, nonsignificant positive response [t(13) = 1.82, P = 0.09]. These observations suggest that the CG-OP and SAL/PMN networks are recruited during the Oddball Effect task contrast.

117Given the historical focus on the salience network and the parietal memory network as separate networks, and their proximity along the posterior midline to the historically defined default network, we replotted the Oddball Effect task contrast on the inflated surface (Fig. 26). For this visualization, the task map threshold was reduced to zero. Much of the full extent of the CG-OP and SAL/PMN networks was strongly activated. The positive response included the posterior midline regions that have been the focus of the parietal memory network as well as the anterior insula region that has been a focus of the salience network. An interesting feature is that islands of the CG-OP network that fell within the frontal midline showed positive responses in the within-individual maps (Fig. 26). These small responses, which were adjacent to large regions with an opposite response pattern, were absent in the group-averaged response (Fig. 26, bottom). The positive response was not selective to these two specific networks, with motor and visual responses as noted above. The positive response also extended into the region of the visual second-order networks (e.g., dATN-B).

118Critically, the networks at or near the historical default network, here estimated within individuals as encompassing at least DN-A and DN-B, were all strongly “deactivated,” meaning more active during the implicitly coded baseline reference than during the salient targets. That is, the contrast replicated the task deactivation pattern that originally generated interest in the default network in the presence of a robust positive response across the distributed extent of the SAL/PMN network. Thus, the separation of the effects along the posterior midline revealed a spatial dissociation between the second-order network SAL/PMN and the third-order networks DN-A and DN-B.

119Distributed throughout association cortex, in the zones roughly7 between the second-order networks, were the five association networks FPN-A, FPN-B, LANG, DN-B, and DN-A (Fig. 27). Among these networks, side-by-side juxtapositions repeated across multiple cortical zones (see I, II, III, and IV in Fig. 27). FPN-A and FPN-B were reliably positioned adjacent to one another and, as a pair, were adjacent to a repeating group of the three other networks, LANG, DN-B, and DN-A. We call these repeating clusters of five networks supra-areal association megaclusters or SAAMs. The reproducibility of the SAAMs across participants was striking and is illustrated for the posterior association zones in all 15 participants in Fig. 28. Although the idiosyncratic spatial details varied, multiple SAAMs were consistently observed. The remaining task analyses explored functional response properties of the association networks embedded within the SAAMs.

120The functional properties of the association networks comprising the SAAMs (FPN-A, FPN-B, LANG, DN-B, DN-A) were explored first in relation to domain-flexible demands on working memory and (in LANG, DN-B, and DN-A Respond Differentially to Distinct Cognitive Domains) in relation to domain-specialized processing functions. The hypothesis was that FPN-A and possibly FPN-B would modulate their response in relation to increasing working memory load across multiple verbal and nonverbal stimulus conditions (see Ref. ). The mapping strategy is illustrated in Fig. A15. On the flattened cortical surface, the within-individual a priori-defined networks FPN-A and FPN-B are displayed in relation to the N-Back Load Effect task contrast (collapsed across stimulus conditions). The N-Back Load Effect task contrast is shown in detail in Fig. 29 for one representative participant. Figure 30 illustrates that the features can be observed in additional participants, and in all participants with available task data in the Supplemental Materials.

121As hypothesized, the N-Back Load Effect task contrast increased activation within and near the boundaries of the FPN-A network and, to a lesser extent, the FPN-B network (Fig. 29 and Fig. A15). The widely distributed response included extensive regions of PFC as well as regions of PPC and the dorsal ACC, all canonical regions associated with domain-flexible cognitive control (e.g., Refs., –). As predicted by the network estimates, there was also a response in LTC and a small subregion of the anterior insula that is spatially distinct from that of other networks. Of equal importance was the consistent absence of response in the distributed association regions linked to the LANG, DN-B, and DN-A networks, including within the PPC and LTC. In essence, the N-Back Load Effect task contrast split the SAAMs and activated the portions linked to the FPN-A network, with minimal or no response in the juxtaposed portions associated with the LANG, DN-B, and DN-A networks.

122To quantify the selectivity of the task response, the mean z values output from the condition-level GLM for the N-Back Load Effect task contrast were calculated separately for each association network. The estimates were obtained within the bounds of each individual’s independent a priori-defined networks and then averaged (N = 15). Results plotted in Fig. 31 reveal a positive N-Back Load Effect response that was strongest in FPN-A [t(14) = 21.67, P < 0.001] and also present in FPN-B [t(14) = 6.45, P < 0.001]. SAL/PMN unexpectedly showed a significant positive response [t(14) = 7.91, P < 0.001] that was significantly weaker than either FPN-A [t(14) = −13.43, P < 0.001] or FPN-B [t(14) = −3.58, P < 0.01]. Thus, although SAL/PMN showed a response, the functional response was less relative to FPN-A and FPN-B, opposite to the pattern found above (contrast Fig. 31 with Fig. 25). The remaining networks, including the three additional networks that were adjacent within the SAAMs, showed a negative N-Back load effect. The effect was significantly negative for DN-A [t(14) = −4.85, P < 0.001] and DN-B [t(14) = −7.14, P < 0.001] but not LANG [t(14) = −0.81, P = 0.43]. These results provide evidence that the FPN-A network and possibly the FPN-B network are involved in processes enhanced by increasing working memory demands, whereas other juxtaposed networks, LANG, DN-B and DN-A, are functionally dissociated, consistent with the qualitative patterns visualized in the activation maps.

123To further investigate the domain flexibility of FPN-A and FPN-B, the mean z values output from the block-level GLM for each of the four stimulus conditions of the N-Back Load Effect task (Face, Letter, Word, and Scene) were separately plotted (Fig. 31). FPN-A [Face: t(14) = 11.74, Letter: t(14) = 16.03, Word: t(14) = 11.30, Scene: t(14) = 12.05; all P < 0.001] and to a lesser degree FPN-B [Face: t(14) = 5.13, Letter: t(14) = 5.60, Word: t(14) = 5.15, Scene: t(14) = 4.54; all P < 0.001] exhibited a significant response across all conditions of the N-Back Load Effect task contrast, supporting that their processing role generalizes across both verbal and nonverbal domains. That is, FPN-A and FPN-B responded to working memory demands more so than adjacent networks, and did so in a domain-flexible manner.

124Among the networks that populate the distributed zones of higher-order association cortex, FPN-A and possibly FPN-B responded in a domain-flexible manner to increasing working memory load. The adjacent trio of networks, LANG, DN-B, and DN-A, did not. In our final planned analyses, we explored the functional specialization of these additional three networks by examining Episodic Projection, Theory-of-Mind, and Sentence Processing task contrasts designed to emphasize distinct specialized domains of higher-order cognitive processing. This analysis framework and set of tests were envisioned prospectively [see Fig. 2 in DiNicola and Buckner].

125The mapping strategy is illustrated in Fig. A16. On the flattened cortical surface, the within-individual a priori-defined networks LANG, DN-B, and DN-A are displayed in relation to the three separate task contrasts simultaneously, to illustrate the adjacency of the responses in relation to each other and to the network boundaries. The details of one composite task contrast map are displayed for a representative participant in Fig. 32. Figure 33 illustrates additional participants, and all participants with available data are shown in the Supplemental Materials. Several results are notable.

126First, the composite activation patterns across the three task contrasts filled in the remaining zones of association cortex. Strikingly, the domain-specialized task responses are situated adjacent to, but separate from, the regions activated by domain-flexible working memory demands (contrast Fig. 32 with Fig. 29). This separation can be seen in many locations, with a particularly clear example visualized within the PPC, where the N-Back Load Effect task contrast showed a posterodorsal response relative to the three current task contrasts. The side-by-side juxtaposition of domain-specialized and domain-flexible regions was also observed within LTC and multiple locations throughout PFC.

127Second, within each juxtaposed cluster of domain-specialized regions, the region preferentially responding to the Sentence Processing task contrast abutted the region preferentially responding to the Theory-of-Mind task contrast, and these abutted the region preferentially responding to the Episodic Project task contrast. Although overlap and exceptions were found, the differential response patterns generally tracked the network separations between LANG, DN-B, and DN-A. The idiosyncratic positions and boundaries of the three networks in any given individual, LANG, DN-B, and DN-A, predicted the positions of the domain-specialized activation responses (Fig. 33).

128Thus, within each local zone the regions associated with the separate networks responded to their distinct specialized cognitive domains. Moreover, the spatially differentiated response patterns repeated across the multiple SAAMs (see I, II, III, and IV in Fig. 32). There were exceptions. For example, regions of task activation in VLPFC did not overlap well with the estimated networks in P12. The discrepancies tended to fall within anterior temporal regions and PFC regions where SNR is low, raising the possibility that technical variance played a role. To reveal the details of the task maps more fully, the Supplemental Materials include task maps for each task contrast separately in addition to the composite maps for all available participants.

129The response was quantified for each of the three task contrasts for each network to formally test for the hypothesized interaction. For each domain-specialized task contrast, the z values within the bounds of each individual’s three independent a priori-defined networks (LANG, DN-B, and DN-A) were obtained and then averaged (N = 13). The resulting mean z values are plotted in Fig. 34. A repeated-measures ANOVA on network-level task response revealed a significant 3 × 3 interaction between the effect of task contrast and network [F(4,48) = 77.82, P < 0.001]. Paired t tests then tested the individual contrasts, with the hypothesis that each network’s within-domain response would be significantly greater than either of the other two networks. All six of these planned comparisons were significant. The Episodic Projection task contrast recruited DN-A regions over those of DN-B [t(12) = 16.38, P < 0.001] and LANG [t(12) = 14.49, P < 0.001]. The Theory-of-Mind task recruited DN-B regions over those of DN-A [t(12) = 5.27, P < 0.001] and LANG [t(12) = 10.09, P < 0.001], and the Sentence Processing task contrast recruited the LANG regions over those of DN-A [t(12) = 6.55, P < 0.001] and DN-B [t(12) = 5.42, P < 0.001].

130Thus, in addition to the qualitative impressions (Fig. 32, Fig. 33, and Fig. A16), statistical tests revealed that the full interaction was significant, with all pairwise tests also significant in support of a triple functional dissociation across the three networks. These observations suggest that the parallel networks LANG, DN-B, and DN-A, with adjacent regions across multiple cortical zones, are specialized to support distinct higher-order cognitive domains.

131Our final post hoc analysis, as a reminder that our results and data patterns leave gaps, explored the possibility that we have not yet isolated the task demands that elicit an FPN-B network response. As depicted in Fig. 31, the response observed in FPN-B is notably weaker in comparison to that in FPN-A. Considering that FPN-B is adjacent to FPN-A, this raises the question of whether the observed response is located within the bounds of the FPN-B network or is due to spatial blur from the adjacent FPN-A network response. To weigh in on this question, we replotted the N-Back Load Effect, Episodic Projection, Theory-of-Mind, and Sentence Processing task contrasts simultaneously on a common surface in a representative participant, P6. The FPN-B network is outlined with a black border. As illustrated in Fig. 35, the regions of the FPN-B network do not display a distinct response in multiple zones including a clear gap in the large parietal association region despite adjacent responses in all of the other network regions. The response to the N-Back Load Effect within the FPN-B network might originate from signal blur from the adjacent FPN-A network. Our final result is thus the unmasking of an ambiguity. The present work does not disambiguate FPN-B from adjacent networks or provide clear evidence that we have found task demands that elicit a response.

132Detailed network estimates reveal a global organization that can be conceptualized as three levels of cortical hierarchy6: locally organized first-order sensory and motor networks, spatially adjacent second-order networks that link to distant regions, and third-order networks that populate and connect widely distributed zones of higher-order association cortex. Repeating side-by-side spatial juxtapositions among the third-order association networks form organized motifs that we call supra-areal association megaclusters or SAAMs. Within each SAAM, the regions linked to distinct association networks demonstrate differential task response properties. Certain networks contribute to domain-flexible cognitive control and others to domain-specialized processes involved in language, social, and spatial/episodic functions. We discuss the practical and conceptual implications of these findings, including how repeating organizational motifs might arise during development.

133In the present work, we explored the utility of a 15-network MS-HBM estimate of cerebral cortical organization that allowed the idiosyncratic details of each individual’s own anatomy to guide the solutions (e.g., Fig. A1 and Fig. A4). The method yielded robust, stable network estimates that were confirmed by analyses of seed region-based correlation (e.g., Fig. A2 and Fig. A5). All quantitative analyses and visual inspections of the data reinforced that the present 15-network estimate captured a great deal of the structured correlations present in the underlying data. From a methodological standpoint, the present results indicate that a MS-HBM can be used to estimate networks automatically and robustly within individuals. Several features of our network estimates revise or expand earlier ideas.

134First, the present network parcellation falls into a class of within-individual network estimates that refine group-based estimates. In group-based estimates, including multiple estimates from our laboratory, large monolithic networks have been identified that encompass extensive regions of association cortex (e.g., Refs., –). The present network estimates are broadly similar but separate the large group-based networks into multiple distinct parallel networks. For example, the network historically known as the default network overlaps four separate networks in the present parcellation including networks LANG, DN-B, DN-A, and SAL/PMN. Each of these four distinct networks can be identified in every individual in the present study. The multiple networks are not estimated to be “subnetworks” with shared regions or anatomical convergence but rather are distinct networks that are near one another and often blurred in group-averaged data (see also Refs.,,, –,, ). Thus, an advance of within-individual network estimates, including the present contribution, is to fully resolve adjacent networks that are difficult to separate through approaches that average over people.

135Second, among within-individual parcellation estimates, we settled on a 15-network solution because of our goal to separate nearby networks within the anterior insula (see Ref. ), as well as to better separate early sensory and adjacent networks. Our analyses confirmed that the newly proposed 15-network parcellation could capture correlational features absent in simpler network solutions, including our own 10-network solution previously estimated in Xue et al.; see Figs. A7, A8, A10, and A11. In addition to detecting distinctions among networks that have close juxtapositions in the anterior insula, the present 15-network parcellation also revealed clear separation of the estimated AUD network from the nearby LANG network. Multiple networks identified in the simpler network solutions remained in the 15-network estimates, indicating that the refinements did not come at the expense of the established networks.

136Third, the present parcellation identified a single distributed network, labeled SAL/PMN, that includes regions that have historically been studied separately as components of the salience network and the parietal memory network. Note that we do not say “joins” two previously described networks, as we suspect there have never been two separate networks. Rather, different research lineages may have focused on distinct regional components of what is ultimately the same network. This hypothesis will require further testing, but several lines of evidence lead to the present proposal that SAL/PMN is a single, coherent network. In every individual, the estimated SAL/PMN network included regions along the posterior midline and within the anterior insula (Fig. 14). Seed region-based correlation patterns recapitulated the automated network estimates: seed regions placed in PFC and posterior cortex revealed clear regional correlation in the anterior insula as well as multiple distinct posterior midline regions (Fig. A12).

137A recent study by Kwon et al. is consistent with the hypothesis of SAL/PMN as a single distributed network. Using both seed-based and clustering approaches on within-individual data, their study revealed a network consistent with the SAL/PMN as delineated in our investigation including anterior insula, the inferior parietal lobule, and the posteromedial regions canonical of the parietal memory network. Additionally, in another study by Lynch et al., a network exhibiting similar correlational structure as the SAL/PMN in our study was identified [see Extended Data Fig. 5 of Lynch et al.]. Furthermore, our independent task data focused on salience processing, via the oddball task, elicited robust responses in the distributed regions of the SAL/PMN network including the posterior midline and anterior insula (Fig. 23, Fig. 24, and Fig. 26). Meta-analyses of task-relevant oddball effects have previously noted responses in the anterior insula and portions of posterior midline. Thus, multiple studies and analyses raise the possibility that the SAL/PMN is a single coherent network.

138Despite these findings, it is important to acknowledge the results that have caused a distinction to be drawn between the parietal memory network and the salience network, including evidence from both anatomical and functional perspectives. The salience network has been emphasized to show strong connectivity to the ventral striatum and the parietal memory network to the posterior hippocampus. Functionally, the salience network has been implicated in detecting important environmental stimuli and orchestrating task switching, whereas the parietal memory network has been emphasized in paradigms involving stimulus repetitions. Given these complexities, further in-depth investigations are required. For example, given Zheng et al.’s findings, the present hypothesis that SAL/PMN is a single coherent network predicts that the posterior hippocampus might respond to low-level oddball events and task transitions. We are presently exploring this possibility (see Ref. ).

139Finally, it is important to note that the present estimates assume (and are optimized to detect) large-scale distributed networks. For this reason, our resultant parcellation is different from parcellations that are optimized to detect local gradients of change and/or directly estimate “area” boundaries [e.g., Refs.,,; for discussion see Refs., ]. Although there is some convergence between approaches, and it is possible to apply mutual constraints, our present parcellation is weighted to estimate networks based on long-range correlational properties, without weighting local gradients.

140A striking observation that is apparent in the flat map visualizations is the recurrent spatial grouping of the same five higher-order networks throughout association cortex (FPN-A, FPN-B, LANG, DN-B, DN-A). The clearest examples are found in PPC and LTC (Fig. 28), but the adjacencies are also present in multiple PFC zones (Fig. 27), as if a shared organizing force plays out repeatedly across different cortical territories. Each grouping of regions possesses similar spatial relations among the five networks: networks FPN-A and FPN-B are next to one another, and that pair of networks is adjacent to the trio of networks LANG, DN-B, and DN-A. These juxtaposed networks form SAAMs not only within the cerebral cortex but also within the cerebellum and the caudate.

141SAAMs possess several additional features. Although their global patterning, meaning spatial adjacencies between networks, is identifiable for multiple SAAMs within and across individuals, the orientations shift and the exact spatial positions vary. For example, within the PPC the axis that begins with the FPN-A/FPN-B pairing and ends with the LANG/DN-B/DN-A triad is oriented dorsal to ventral. Within the LTC, the axis is oriented ventral to dorsal (Fig. 27). Moreover, although the SAAMs are readily identified in every person in the PPC and LTC, usually with a discontinuity between the two SAAMs, the idiosyncratic spatial details vary from one person to the next. In some individuals the two zones appear fused (Fig. 28). It is thus unsurprising that group-averaged data, while revealing certain spatial features apparent in the SAAMs, blur over the fine spatial details that are robust and consistently apparent in the within-individual maps.

142The spatial juxtapositions that define the SAAMs in PPC and LTC are also present in multiple zones of the PFC. However, there is not always spatial separation. The boundaries of individual SAAMs in PFC are thus ambiguous. In Fig. 27 we note candidate SAAMs in VLPFC (labeled III) and DLPFC (labeled IV), recognizing that these are hypotheses. A future endeavor might explore how a repeating pattern could parsimoniously explain the juxtapositions in PFC with the assumption that multiple SAAMs are present like those observed in PPC and LTC, but with the additional complication that there are multiple adjacent SAAMs that collide into one another, perhaps as a consequence of their formation during development.

143A final detail regarding the SAAMs is subtle but potentially informative. Although the presence of five regions linked to the distinct networks is a consistent feature of PPC, LTC, VLPFC, and DLPFC, there are also partial sets of the network juxtapositions in other cortical zones. For example, along the midline there is clear representation of networks DN-A and DN-B in PMC and medial PFC (MPFC) but not consistently the other networks (Fig. 27). The partial SAAMs may provide an insight into the origins of the patterning. DN-A is a putative hippocampal-cortical network that has been extensively studied in humans (e.g., Refs.,,,,, ) and monkeys (e.g., Refs. –,, ). The hippocampal formation, via polysynaptic projections through entorhinal cortex and PHC, projects heavily to RSC and ventral PCC along the posterior midline, and also to MPFC. The exclusive assignment of PHC to DN-A and the predominance of DN-A along the midline may thus reflect connectivity to the hippocampal formation. The interdigitation of DN-A with other higher-order networks might emerge as the hippocampal-predominant projections intermix with other anatomical projection gradients in the apex association zones where the fully formed SAAMs are present.

144The default network, or default mode network, has received considerable attention among investigations of cerebral networks. In relation to estimating networks using resting-state functional connectivity, after the seminal description of the method, the default network was the first distributed association network to be characterized in humans and in monkeys. All group-based network estimates, even low-dimensional solutions that identify as few as seven networks, find a large, distributed network that has the spatial pattern of the default network (e.g., Refs. –). Thus, a critical issue to address, given the historical emphasis on the default network, is how the present network estimates relate to these earlier descriptions.

145Our present hypothesis is that the large monolithic (or core subnetwork) descriptions of the default network based on group-averaged data, including our own contributions (e.g., Refs., ), are inaccurate and reflect an artifact of spatial blurring. As noted above, the canonical group-averaged default network overlaps fully or partially four distinct networks: LANG, DN-B, DN-A, and SAL/PMN. The separation of these networks is anticipated in some prior group-based analyses. For example, Andrews-Hanna and colleagues noted that regions within PPC responding to social inference (theory of mind) tasks tended to activate an anterior region relative to tasks targeting remembering. This distinction likely captures the separation of DN-B and DN-A in PPC. Similarly, in a thorough analysis of functional connectivity in group data, the network identified here as LANG was separated from the canonical default network (Ref.; for discussion see Ref. ). However, the blurring induced by between-subject averaging, to date, has negated the ability to resolve the spatial details that fully distinguish the four nearby networks that are described here.

146A further observation emerges from our task-based results. In addition to the challenge of identifying the multiple, juxtaposed networks due to spatial blurring, there is a separate functional property that has anchored study of the default network. The default network was originally described based on task-induced deactivations, referring to the observation that the distributed association regions that comprise the default network are more active in passive tasks and fixation than active, externally-orientated tasks [Refs.,,; see Buckner and DiNicola for review]. When a contrast is made between active and passive tasks, a distributed pattern of “deactivations” emerges that is robust and overlaps with group-based estimates of the default network. What is surprising and interesting is that, even with the present high-resolution within-individual estimates, the task-induced pattern of deactivation remains broad and spans multiple networks.

147Specifically, the Oddball Effect task contrast reveals a broad task-induced deactivation pattern within individual participants (Fig. 26). That is, the regions deactivated by attending and responding to external stimuli span multiple association networks even when group averaging is not a factor. Figure 25 quantifies this effect: DN-A, DN-B, and LANG all show significant “deactivation,” with DN-A and DN-B being almost indistinguishable from one another, despite clear functional double dissociation during domain-relevant active tasks [e.g., Fig. 34; see also DiNicola et al.]. One possibility is that, although DN-A and DN-B are anatomically and functionally distinct networks, they may collectively be suppressed during certain externally oriented task events, perhaps as a result of a broad antagonistic process between externally versus internally oriented processing modes (Ref.; see also Refs. –). Thus, the phenomenon of task-induced deactivation, which is not selective to specific networks, may have reinforced an impression that there is a coherent monolithic function across large swaths of association cortex, a possibility refuted by a growing number of robust functional dissociations (e.g., Refs.,,,, ).

148Another relevant observation surrounds the relation between the default network and the present estimate of network SAL/PMN. The SAL/PMN network possesses regions distributed across the cortex, including multiple distinct regions along the posterior midline side by side with DN-A and DN-B network regions. The adjacencies make the regions easy to confuse. Despite their spatial proximity, Zheng et al. (Ref.; see their Fig. 6) noted that “deactivations” are restricted to the default network and separate from their estimate of SAL/PMN (labeled as the parietal memory network in their paper). The transient positive response in SAL/PMN observed here to salient oddballs is robust, including the regions along the posterior midline, separate from the juxtaposed DN-A and DN-B regions showing deactivation (Fig. 26). Moreover, SAL/PMN has small, focal regions of response in MPFC, which are also surrounded by DN-A and DN-B network regions. In the group-averaged map displayed in Fig. 26, bottom, there is no detectable positive response in MPFC. Each individual shows a response but in slightly different spatial positions from one person to the next. The positive task response in MPFC is likely lost in the process of spatial averaging.

149Our results thus converge with Zheng et al. to suggest that SAL/PMN is spatially and functionally distinct from the network historically described as the default network. The SAL/PMN network does not exhibit task-induced deactivation; rather, it displays an opposite functional response pattern: transiently activating to salient external task events, including in both posterior and anterior regions along the midline.

150Paul Flechsig contributed the powerful but simple idea that the cerebral cortex develops sequentially, radiating outward from motor and sensory cortex (see Refs. – for translations and discussion; see Ref. for further context). The basis of Flechsig’s hierarchy was the developmental timing of myelination of the fibers reaching the cortex. By his account “in the cerebral convolutions, as in all other parts of the central nervous system, the nerve-fibers do not develop everywhere simultaneously, but step by step in a definitive succession” [translated in Clarke and O’Malley, p. 548]. The cortical motor and sensory (and certain limbic) zones myelinate first. Next are the intermediate zones that surround the motor and sensory zones. The terminal zones myelinate in the final stage, beginning ∼4 mo after birth, and encompass prefrontal, temporal, and parietal regions thought of today as higher-order association cortex. The prescient lens of hierarchical cortical organization provides a framework to understand our findings.

151Specifically, the candidate assignments of first-, second-, and third-order networks6 are motivated by (and agree well with) Flechsig’s reference maps of sequential myelination (Fig. 36). In particular, the distributed regions late to myelinate (the terminal zones) are positionally similar to our estimated association zones containing the five higher-order networks that make up the SAAMs. These same general zones were emphasized more than a century ago as the regions distinguishing human and ape brains from those of smaller monkeys and have been supported, based on modern comparative anatomical approaches, to be disproportionately expanded in humans relative to monkeys. Taken together, the global spatial relations among networks (Fig. 17, Fig. 18, and Fig. 19) and the repeating fractionation of the higher-order association zones into five networks (Fig. 27 and Fig. 28) are consistent with processes that organize the cortex through distinct developmental stages.

152In a hypothesized first stage, cortical networks might progressively organize from the early sensory and motor areas that themselves are patterned through structured inputs. For example, retinotopic organization is imparted on early visual cortex via spontaneous retinal activity waves that are present before birth and carried to the cortex through the thalamic nuclei. These early organizing events may anchor the formation of the retinotopic clusters captured by our estimates of the VIS-C and VIS-P networks. The second-order networks may then organize tethered to these first-order networks, but with progressively more distributed regions, corresponding to Flechsig’s intermediate (or border) zones. Averbeck and colleagues have also proposed a similar nesting of networks outward from the primary motor and somatosensory areas (S1-M1) based on extensive analyses on anatomical connectivity patterns (Refs.,; see also Refs., ). The zones that generally fall between the regions of the second-order networks include large swaths of prefrontal, temporal, and posterior parietal association cortex that correspond to Flechsig’s terminal (or central) zones and are hypothesized to be the last to develop, forming our hypothesized third-order networks (Fig. 36). Thus, much of the cortical mantle may be patterned by a series of networks that nest outward from the primary cortical areas.

153In a second developmental stage, we hypothesize that, as the networks sequentially form, they may undergo a second process of fractionation and specialization. Our proposal of a distinct second process is specifically put forth to explain how juxtapositions might arise similarly across widely distributed (noncontiguous) zones of cortex, such as observed for the distinct SAAMs in LTC and PPC (Fig. 28). A specific prediction of this hypothesis is that, as development progresses, activity-dependent processes may eliminate and/or stabilize synapses that support specialization evident in the adult (for discussion see Refs. –). In the cortical mantle of humans, the expanded association zones may fractionate and specialize into the multiple juxtaposed networks that support higher-order cognition.

154By combining network estimation and task-based fMRI within the same individuals, the present results provide insight into the functional specialization of the networks. A broad observation was that the second-order distributed networks SAL/PMN and CG-OP were dissociated from the third-order association networks via their robust, transient response to oddballs (Fig. 25) consistent with prior studies. None of the third-order networks that populate the SAAMs displayed a robust transient positive response. In fact, four of the five networks within the SAAMs (FPN-B, LANG, DN-B, and DN-A) showed a significant negative response (network FPN-A was equivocal). By contrast, all of the third-order association networks responded robustly to ongoing task demands with distinct forms of functional specialization as described below.

155A first robust dissociation among the third-order networks came in their differential response to working memory demands. FPN-A in particular responded to high memory load in the N-Back Load Effect task contrast and did so similarly across verbal and nonverbal materials (Fig. 31). FPN-B’s response was quantitatively lower and may reflect spatial blur from FPN-A (Fig. 35).8 Furthermore, FPN-A displays the same general spatial pattern as the previously described multiple-demand network. Our data are thus convergent with the existing literature to suggest there is a distributed frontal-parietal network (or networks) that responds when tasks become more effortful, perhaps related to processing functions associated with cognitive control (e.g., Refs., ). The within-individual precision mapping allowed spatially precise network estimates to be made of FPN-A that predicted the idiosyncratic response patterns across participants.

156A few further details are of interest. First, across most individuals, FPN-A included a small region in the anterior insula (labeled in Fig. 29). This small region showed a N-Back load effect response surrounded by spatially distinct components of the CG-OP and SAL/PMN networks. Although separations are shown in some group-averaging parcellations, it can be easy to blur over or miss this buried insular region in group analyses. Our present estimates suggest that the anterior insula is a particularly challenging region of the cortex to study because multiple, distinct networks are spatially juxtaposed near where the cortex folds onto itself in the volume. Examinations of group data, especially data that average across participants within the volume, may be particularly vulnerable to distorting functional properties of this region.

157Second, the spatially circumscribed regions of each SAAM that aligned to FPN-A, with some exceptions in low-SNR regions, tended to show a robust N-Back load effect (Fig. 29 and Fig. 30). The adjacent regions of LANG, DN-B, and DN-A did not. Thus, the N-Back load effect functionally dissociated FPN-A from the LANG, DN-B, DN-A cluster multiple times across the distributed association zones including PPC, LTC, VLPFC, and DLPFC.

158Fedorenko and colleagues have previously noted that regions of the multiple-demand network lie side by side with functionally distinct domain-selective regions (specifically in the language domain). Our present results are consistent with their observations and reinforce that the functional dissociation is a general property of the association cortex including close spatial juxtapositions in temporal and parietal cortical association zones, not only within PFC. We interpret the repeating pattern across the cortical mantle to reflect that functional specialization is a property of the networks, including all their distributed regions (see also Ref. ). Furthermore, robust functional dissociations were present for higher-order cognitive domains beyond language (see also Ref. ). That is, although networks LANG, DN-B, and DN-A did not modulate in a domain-flexible manner to working memory demands, each network responded robustly and selectively to a distinct specialized domain of higher-order cognition.

159The most striking functional observation of the present study was the robust triple dissociation across networks LANG, DN-B, and DN-A, as predicted by DiNicola and Buckner. The LANG network responded when participants processed meaningful sentences; the DN-B network when participants engaged Theory-of-Mind tasks; and the DN-A network when participants remembered from their past or contemplated a personal future scenario. The triple dissociation was carried by a formal statistical interaction (Fig. 34) and could be visualized qualitatively on the flat maps of individual participants (Fig. 32 and Fig. 33; see also Refs., ). Considering that until recently we and others conceptualized these zones of association cortex as being deployed flexibly across a range of higher-order cognitive domains (e.g., Refs.,, ), this is a major revision to our understanding.

160Our composite results suggest that higher-order association cortex possesses at least three domain-specialized parallel networks supporting language, social behaviors, and remembering the past and imagining the future. These domain-specialized networks are themselves separate from domain-flexible networks that participate in cognitive control. We do not know how these networks interact or whether they remain functionally separate across multiple task classes, but the robust dissociations among juxtaposed regions demonstrated here suggest that there is more modularity in association cortex, including PFC, than has typically been considered.

161A key limitation of the present work is the reliance on correlational, indirect methods to infer network organization. The caveats surrounding interpreting such network estimates, and the empirical tests of their utility despite known limitations, are discussed elsewhere. Specific to the present work, it is notable that the boundaries in networks generally predicted task response patterns, bolstering confidence that valid organizational features are being described. However, there were exceptions and regions of mismatch, consistent with poor signal quality around the sinuses and inner ear (see Fig. 1 and Fig. A3). Network estimates in these poorly sampled regions of cortex may be distorted.

162There are also limitations to our modeling approach. In choosing the present parameters of the MS-HBM used to estimate networks, decisions were made that influence the estimates. Specifically, we choose to model 15 networks and initiated the model with a prior that arose from a group-averaged dataset. As the seed region analyses verified, the model captured within-individual correlational properties well, but not perfectly. Thus, a limitation in our present model is knowing whether one could do better and whether our specific decisions imparted bias. We assume the answer is yes to both questions. As our own work has evolved from a relatively crude 7-network estimate in average participant groups to a 10-network estimate within individuals, we expect the that the present network estimates will be refined further and eventually replaced. As a specific example, it is unclear that the present model fully captures the details of the recently described intereffector connectivity pattern. The structured correlations they observed, and we also find, are partially incorporated in our estimate of the CG-OP network (Fig. 22) but not entirely (e.g., see diamond in Fig. 20). Placing seed regions within the body map gap in the SMOT-B network reproduces the intereffector regions as well as the distributed pattern of the CG-OP network (Fig. 22), indicating that the spatial prior of the CG-OP network5 might benefit from revision in a future iteration. With the expectation of further improvements in mind, we are struck by how the present network parcellation captures multiple functional dissociations prospectively in task data, including idiosyncratic and small regions of response.

163Another limitation is that, although the task contrasts allowed for robust functional dissociations, the tasks were designed and implemented to differentiate networks, which is a different goal from interrogating in detail a hypothesized cognitive operation. That is, the limited task data we collected fall far short of systematically manipulating variables to clarify the component computations performed by each of the networks. As an example of a task domain needing further exploration, the Oddball Effect task contrast maps were fairly broad, spanning both the CG-OP and SAL/PMN networks and extending beyond them. One possibility is that we have not yet found the right task demands to dissociate the two networks. Another possibility is that the two networks do not respond differentially. The CG-OP and SAL/PMN networks may respond together to common task events linked to surprise and novelty. Within this possibility, rather than being distinguished by what they respond to, they may be distinguished by their anatomical interactions with other networks. CG-OP and SAL/PMN are both distributed networks that are near one another in many zones of cortex, but they are also adjacent to distinct networks in other portions of cortex. That is, CG-OP and SAL/PMN may both be networks that respond to transient orienting events but are ultimately functionally distinguished by interactions with their extended partner regions. The present work did not explore or test such possibilities.

164Another future direction pertains to the need to better understand the relation of traditional area estimates with the present network estimates. By “area” we mean the demarcation of regions of cortex as separate, defined zones using functional, architectonic, connectivity, and topography constraints. We previously noted discrepancies between functional connectivity patterns and areal boundaries (e.g., Refs.,, ), as have others (e.g., Ref. ). There are two topics to be considered.

165First, for regions of cortex that have well-recognized areas, our network borders do not align with the areal borders (e.g., Fig. A13E). For example, networks VIS-C and VIS-P group together V1, V2, and V3 and split them roughly along the 5° eccentricity line. The estimated networks likely reflect the dominant anatomical connectivity gradient that, within early visual areas, progresses along eccentricity. The V1/V2/V3 areal boundaries are distinguished by a local reversal in polar angle along the horizontal meridian. Thus, connectivity transitions between early visual areas are relatively subtle [for further discussion of this issue see Buckner and Yeo]. The somatomotor networks similarly group M1/S1 and multiple body maps that span architectonically distinct areas. One future direction is to understand the relation of the networks estimated here and the finer-scale anatomical differences that demarcate adjacent areas.

166The second related topic is the relation between the present network estimates and architectonic features in higher-order association cortex. This is a trickier topic. Varied perspectives have been put forth on whether association cortex possesses sharp areal boundaries that parallel those found in sensory systems (for discussion see Refs.,,, ). There is also an open question of whether, in practice, there are known stable features that can define areal borders in association cortex, especially when architectonics are considered in isolation (e.g., Refs., ). We will not resolve the debate here, but some of our observations are relevant to the discussion.

167Most critically, the extent and complexity of the network juxtapositions encompassed within the SAAMs are of such a spatial scale that they seem unlikely to align to traditional architectonic borders, at least those reflected in any of the commonly used atlases. In the spirit of supra-areal clusters reported in the visual system (see Refs.,; see also Refs., ), we refer to the repeated groupings of multiple association networks as supra-areal association megaclusters specifically to reinforce the possibility that they might span and split traditional architectonic patterns. One possibility is that future advances will find architectonic features that align to the transitions between SAAMs as well as between the multiple network regions within the SAAMs (perhaps via spatial transcriptomics; Ref. ). Alternatively, there may be broad patterning forces during development, such as captured in Flechsig’s maps of sequential myelination, that reflect processes that guide where SAAMs develop but that do not specify the details of the borders and regional specializations within the SAAMs. The local spatial arrangements might be carried by extrinsic anatomical connectivity differences that refine relatively late in development through activity-dependent processes, without rigid alignment to architectonic features.

168Data to inform these and other possibilities will emerge as the field charts development of association networks in nonhuman primates with direct anatomical techniques and in human infants with noninvasive approaches.

169The present study examined the organization of cerebral networks within intensively sampled individual participants. We provide the resulting network estimates and the raw data used to derive them as open resources for the community. Our initial explorations on the data uncovered a hierarchical organization that distinguishes three levels of cortical hierarchy: first-, second-, and third-order networks. Notably, regions of distinct third-order association networks consistently displayed side-by-side juxtapositions that repeated across multiple cortical zones, with clear and robust functional specialization among the embedded regions.

170The DU15NET atlases are available (https://freesurfer.net/fswiki/CorticalParcellation_DU15NET). Additional parcellation, SNR, and task maps related to this manuscript are available on Balsa (https://balsa.wustl.edu/study/zK166). Task descriptions, contrast descriptions, and code are provided on Harvard Dataverse (https://doi.org/10.7910/DVN/AVB4BW). Individual participant data are available in the NIH repository (https://nda.nih.gov).

171This work was supported by NIH Grant MH124004, NIH Shared Instrumentation Grant S10OD020039, and National Science Foundation (NSF) Grant 2024462. L.M.D. was supported by NSF Graduate Research Fellowship Program Grant DGE1745303. W.S. was supported by the Paul and Daisy Soros Foundation.

172Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

173No conflicts of interest, financial or otherwise, are declared by the authors.

174J.D., L.M.D., P.A.A., N.S.-G., W.S., M.C.E., and R.L.B. conceived and designed research; J.D., L.M.D., P.A.A., N.S.-G., and R.L.B. analyzed data; J.D., N.S.-G., W.S., S.K., J.L., and R.L.B. prepared figures; J.D., L.M.D., P.A.A., N.S.-G., and R.L.B. interpreted results of experiments; J.D. and R.L.B. drafted manuscript; J.D., L.M.D., P.A.A., N.S.-G., W.S., S.K., J.L., A.X., B.T.T.Y., M.C.E., and R.L.B. edited and revised manuscript; J.D., L.M.D., P.A.A., N.S.-G., W.S., S.K., J.L., A.X., B.T.T.Y., M.C.E., and R.L.B. approved the final version of the manuscript.