Published: July 19, 2017

Introduction {#sec1}
============

The cerebral cortex possesses a complex tapestry of networks that interact and compete in the service of information processing. Building on a history of assigning specialized functions to brain regions, early seminal work by Norman Geschwind, Marsel Mesulam, and others proposed ideas about how distributed regions might interact to perform high-level tasks (e.g., [@bib34], [@bib35], [@bib69], [@bib70]). A leap in progress occurred when networks began to be conceptualized within the framework of anatomical connectivity patterns in the macaque, following the availability of both retrograde ([@bib72]) and anterograde ([@bib18]) tracers. For example, by charting the laminar pattern of anatomical projections, a hierarchical pathway emerged from areas involved in visual perception through to areas enabling motor actions ([@bib65], [@bib101], [@bib3], [@bib11]; see [@bib92] for discussion). This canonical distributed network comprises primary and secondary visual areas, parietal association areas, and motor areas. These interconnected areas form a partially modular distributed network that interacts with, but is anatomically distinguishable from, other processing hierarchies ([@bib24], [@bib101], [@bib106]; see also [@bib64]). Notably, this canonical network involves areas distributed across temporal, parietal, and frontal cortices. As will be illustrated below, this distributed pattern is a general motif that is apparent across multiple large-scale networks (see [@bib38]).

Human neuroimaging studies are particularly useful for characterizing distributed networks because they survey the whole brain at once. [@bib17] highlighted a network currently known as the "dorsal attention network" (dATN) in the human neuroimaging literature, which is likely homologous to later-stage components of the sensory-motor hierarchy described above in the macaque ([@bib108], [@bib81]). Network analysis based on intrinsic functional connectivity (FC) consistently reveals the dATN (e.g., [@bib30], [@bib109]; see also [@bib7]). Detailed analysis in relation to retinotopic areas recapitulates, to first approximation, the full hierarchical pathway described in the macaque ([@bib113]). Thus, while there are caveats to interpreting networks observed by FC ([@bib15], [@bib77], [@bib97], [@bib83]), the results can generate hypotheses about the organization of multiple large-scale networks that populate human association cortex.

Group-based studies using FC suggest that association cortex comprises about five major distributed networks (e.g., [@bib113], [@bib82], [@bib20]). These networks have sufficiently modular properties to be consistently identified as isolated networks, even when constraints are relaxed to emphasize interactions between networks ([@bib114]). Among these, the "default network" (DN) has been extensively studied. The DN is anatomically separate from the dATN and is estimated to have expanded in hominin evolution ([@bib13], [@bib63]; see also [@bib46]). Regions across the DN are correlated to a level akin to the local sensory and motor networks despite being widely distributed ([@bib82]; see also [@bib41]). Tracer studies suggest that the DN discovered in the human may have an anatomically connected homolog in the macaque ([@bib14], [@bib8]; see also [@bib62], [@bib13]). Several other networks, near the dATN and DN, each have their own distributed organization ([@bib113], [@bib82], [@bib20]).

Several principles emerge through examining relations among multiple large-scale networks. First, the distributed networks all have an organization that is similar to the dATN and DN, with each network possessing frontal, temporal, posterior parietal, frontal midline, and posterior midline components. Depending on how the network is estimated, certain components can be underemphasized, e.g., the dATN is sometimes described without a frontal midline component, especially when "winner-take-all" network assignments are used (e.g., [@bib113]). Targeted analyses often reveal there is a midline component (e.g., Figure 5 in [@bib30]; Figure 32 in [@bib113]). The distributed networks follow a general motif that is *roughly* conserved even though each network contains spatially distinct regions. [@bib82] further noted that the networks have similar spatial arrangements at their interfaces in multiple zones of cortex (e.g., temporal, parietal, frontal). That is, if a network lies side by side with another network in parietal cortex, it is also likely to do so in frontal cortex. This echoes features proposed by [@bib38] (her Figure 4) and [@bib69] (his Figure 4) based on anatomical data. Second, while the broad architecture follows a general motif, there are differences that distinguish the networks. For example, while the dATN is coupled to visual and motor regions, the DN is instead coupled to the hippocampal formation, perhaps reflecting a mnemonic functional anchor ([@bib42], [@bib107], [@bib14], [@bib5]).

The discovery that cerebral association cortex possesses multiple distributed networks is a major milestone for the field. However, reliance on group-averaged estimates raises questions. A first open question concerns whether the present convergence on major networks reflects the correct level of description. Similarities between human neuroimaging and monkey anatomical findings suggest that the major networks described to date capture true organizational features. Nonetheless, several reports note distinctions that fractionate the larger major networks, either locally within a region or in complex ways across distributed regions. For example, [@bib97] presented a clustering analysis that illustrated both high-level structure that recapitulated the major networks and a substructure that included spatial distinctions between nearby regions of cortex (see also [@bib96]). In the study by [@bib113], they showed that the major networks broke down further when more fine-grained network structure was examined (their 17-network parcellation). [@bib82] noted at least two network subcomponents ("subgraphs" within their framework) that were not captured in the major network descriptions. In detailed analysis of the DN, there have been several proposed schemes to delineate subnetworks (e.g., [@bib5], [@bib6], [@bib58], [@bib59], [@bib12]). We do not attempt to integrate these various findings here but point out that analyses of group data have always possessed features that are not fully accommodated by assuming a small number of major networks.

A second open question is whether there are global properties of organization that span across the networks. [@bib63] recently proposed a macroscale organization of networks that moves outward from sensory-motor networks on one end to the DN on the other end. This echoes an idea advanced by [@bib71] that the cortex exhibits a hierarchical organization progressing from unimodal areas to integrative transmodal areas (see also [@bib51]). The repeating motif and spatial juxtaposition of multiple networks suggest that the macroscale organization of the distributed networks might be partially explained by developmental anchors or gradients ([@bib13], [@bib63], [@bib53]). Detailed analysis of cortical network architecture may provide evidence for or against macroscale gradients that span networks.

To make progress on these questions, we need to move from group-level description to the finer spatial scale that is accessible when studying organization within the individual. The majority of the literature on network architecture is based on averaging data from groups of spatially normalized individuals. Group averaging was necessary in positron emission tomography (PET) studies as radiation places severe restrictions on repeat imaging. A major technical breakthrough came with the ability to average brain volumes across individuals to boost signal to noise (e.g., [@bib29], [@bib21], [@bib33]). However, group-averaging approaches are limited due to blurring over anatomical and functional variability ([@bib99], [@bib94]).

Anatomical variability refers to the complex geometry of the brain's gross structural morphology (e.g., sulci and gyri) that must be normalized across subjects. Surface-based normalization has improved this technical hurdle for the cerebral cortex ([@bib26], [@bib105], [@bib86]), but misregistration is still a challenge. Functional variability refers to the spatial arrangement of functional zones on the cortical surface. Functional organization is likely derived from microscopic anatomical features, such as areal organization or anatomical subdomains that possess distinct connections. The point here is that these local organizational properties are not fully accounted for by gross anatomical differences between individuals. Histological data have illustrated that the border locations, size, and shape of architectonic areas vary on the surface between individuals (e.g., [@bib84], [@bib85], [@bib1], [@bib2], [@bib16], [@bib27], [@bib45]). Thus, even if individuals could be brought into complete anatomical alignment, spatial blurring across architectonically defined areas would persist. An additional consideration is that cortical organization may possess complex topography that has side-by-side juxtapositions and finely interleaved structure. Fine spatial arrangements may reflect the complex geometry of small cortical domains (e.g., [@bib39], [@bib37]), such as has been hypothesized for face patches ([@bib74]), or the fine-grained structure due to topographically arranged projections such as is present for eccentricity in the early visual system (e.g., [@bib65]). Normalization across subjects may blur these spatial features.

Recent human neuroimaging studies focusing on individuals have noted fine spatial details that are lost or attenuated by group averaging. [@bib23] demonstrated that a language-preferential region in prefrontal cortex (PFC) appears "as an island" between regions showing domain-general response properties. The exact positioning of the island moves from person to person. [@bib73] characterized interdigitated regions of PFC belonging to separate auditory and visual networks whose exact positions vary considerably across subjects (Figure S1 in [@bib73]). [@bib55] recently highlighted a network implication of closely juxtaposed small regions. Two distinct networks were found in an individual using nearby seed regions in left PFC; however, when the same two regions were applied to group-averaged data, a single distributed network emerged (Figure 7 in [@bib55]). In an analysis motivated by architectonic subdivisions of nearby prefrontal areas BA 44 and BA 45, [@bib49] (their Figure 11) observed a similar separation of networks in the individual. [@bib37] noted that whereas in the majority of cases a single contiguous region of lateral PFC was found to belong to a left-lateralized language network (Figure S8 in [@bib37]), a minority of individuals displayed two "split" non-contiguous regions. In a comprehensive analysis of individual-specific network features, [@bib40] identified numerous reliable features that were not captured in group estimates of network organization. These features form a distributed set of patches across the cortical mantle that were "too infrequently present and/or spatially variable relative to their size to emerge in group-average data." Thus, analysis of individual brains reveals regional and network features that are underemphasized in group-averaged analyses.

Motivated by these findings, we conducted an extensive set of analyses focused on the individual. We discovered that the DN could be reliably subdivided into parallel networks within the individual. Similar separations were made for other major networks. Regions of the separate networks lay side by side across distributed zones of cortex and were sufficiently variable between individuals to obscure their existence in group-averaged analyses. To make these observations, we analyzed data from four individuals scanned repeatedly over 24 sessions.

Results {#sec2}
=======

High Signal-to-Noise and Full-Brain Coverage Was Achieved in Each of Four Individuals {#sec2.1}
-------------------------------------------------------------------------------------

The present study acquired extensive data over many functional MRI (fMRI) sessions in the same individuals. Two estimates of data quality, slice-based temporal signal-to-noise ratio (tSNR) and fractional amplitude of low-frequency fluctuations (fALFF; [@bib115]), were calculated for each participant. Across each subject's 24 scans, slice-based tSNR ranged between 127.2 and 336.3 for S1, 264.0 and 461.2 for S2, 123.2 and 286.6 for S3, and 204.7 and 493.1 for S4. Absolute motion (i.e., the accumulated displacement in each run) ranged between 0.303 and 1.505 mm for S1, 0.209 and 1.181 mm for S2, 0.474 and 1.651 mm for S3, and 0.267 and 1.765 mm for S4. Spatial properties of tSNR and fALFF are affected by susceptibility artifacts (e.g., [@bib78]). [Figure S1](#mmc1){ref-type="supplementary-material"} displays voxel-based tSNR and fALFF maps projected onto the cortical surface. High voxel-based tSNR and fALFF were achieved across nearly the entire cortical mantle, including ventral PFC and portions of the anterior and ventral temporal lobe.

Distinct Distributed Networks Fractionate the Canonical Default Network within the Individual {#sec2.2}
---------------------------------------------------------------------------------------------

The goal of the present study was to identify spatially detailed features of network organization in the individual. For each participant, half the data were used for discovery (n = 12) of network features that were later tested in the remaining independently collected sessions (n = 12). The discovery datasets were analyzed blind to the hypothesis-testing datasets. An interactive seed selection and FC map viewing platform was established using a bespoke high-resolution cortical mesh so as to minimize spatial blurring during interpolation and allow networks to be defined with high precision by selecting individual vertices.

In the discovery data, the observation was made in one subject that two similar, but distinct, networks could be observed from two seed vertices selected from adjacent regions of left lateral PFC ([Figure 1](#fig1){ref-type="fig"}). Both networks were distributed across inferior parietal, lateral temporal, medial PFC, and posteromedial cortices and resembled the canonical DN but with distinct nodes. The two networks had components that were closely neighboring but separate in numerous cortical zones, including temporal and ventral medial PFC, suggesting that the DN may be comprised of distinct neighboring networks. We refer to these hypothesized networks as Network A and Network B ([Figure 1](#fig1){ref-type="fig"}). Investigating additional participants yielded a similar distinction.Figure 1Two Parallel Interdigitated Distributed Networks at or near the Canonical Default Network Are Revealed by Functional Connectivity within IndividualsEach row illustrates functional connectivity (FC) maps from a single subject (S1--S4). Two networks were observed in each individual. Subject-specific seed regions were placed in the left lateral PFC of the discovery dataset (white filled-in circle). The seed region labeled A^∗^ yielded Hypothesized Network A (left) and the seed region labeled B^∗^ yielded Hypothesized Network B (right). Note that throughout the cortex, Networks A and B are adjacent to one another with slightly varied positions from individual to individual. Hypothesized Network A includes posterior inferior parietal lobule, lateral temporal cortex, ventromedial PFC, retrosplenial/ventral posteromedial cortex, and parahippocampal cortex. Hypothesized Network B includes the temporoparietal junction, lateral temporal cortex, an inferior region of ventromedial PFC, a dorsal region of anteromedial PFC, and posterior cingulate cortex. Regions (hollow circles A and B) were selected to formally test the distinction between the two networks in independent data. The surfaces are rotated by 19° along the y plane to better show the ventromedial PFC and intraparietal sulcus. The same views are used in accompanying figures.

An important difference between the networks was that Network A showed correlation with a region in parahippocampal cortex (PHC), whereas no such evidence was found in Network B, even when additional and more sensitive analyses were carried out focused on this region ([Figure S2](#mmc1){ref-type="supplementary-material"}). The presence of correlation with medial temporal lobe structures is anticipated (e.g., [@bib42], [@bib107]). However, that the coupling was to one hypothesized network and not the other was unexpected.

The recurrence of the distributed pattern of regions across the four participants provided strong evidence for two dissociable networks ([Figure 1](#fig1){ref-type="fig"}).

Double Dissociation of the Two Networks within the Individual {#sec2.3}
-------------------------------------------------------------

The independent replication data were used to formally test the hypothesis that Network A was dissociable from Network B. Using only the discovery datasets, a priori regions (single vertices) were selected that maximized the separation of the networks from the two lateral PFC seed regions in the main regional zones of the cortex (temporal, inferior parietal, posteromedial, and medial PFC) and a region in PHC. These a priori regions were targeted within each subject to locations where contiguous vertices could be observed in Network A and Network B ([Figure 1](#fig1){ref-type="fig"}). Using the independent replication dataset (n = 12 in each subject), time courses were extracted from each subject's a priori regions. The r-to-z transformed Pearson's product moment correlation was computed between the two lateral PFC seed regions and each distributed test region ([Figure 2](#fig2){ref-type="fig"}). Two-way ANOVA was used to test the dissociation between the two PFC seed regions and the regions in each cortical zone. The critical test was whether there would be significant interactions between seed and target regions belonging to Networks A and B. The presence of interactions across all distributed zones of the cortex would provide strong evidence that there was a network dissociation across the cortex.Figure 2Parallel Distributed Networks Are Statistically Dissociated Using Independent DataFunctional correlation strength was computed between the two PFC seed regions yielding Network A and Network B and the pairs of adjacent seed regions in lateral temporal (Temporal), inferior parietal (Parietal), Posteromedial, and Medial PFC cortices (regions shown in [Figure 1](#fig1){ref-type="fig"}). This yielded a 2 × 2 contrast for each zone of cortex (e.g., Networks A and B's PFC regions against Networks A and B's Temporal regions). An additional seed region in PHC was grouped with Network B's posteromedial region. Correlations with Network A's PFC region are shown in yellow and Network B's in red. Bars represent the mean from the 12 sessions of the hypothesis-testing dataset with SEM. All 20 2 × 2 ANOVA tests were significant (^∗∗^p \< 0.01), with most showing a cross-over interaction.

Critically, the multiple tests across cortical zones and subjects targeted convergent evidence for network dissociation. Twenty 2 × 2 ANOVAs were carried out across four subjects and five cortical zones (temporal, inferior parietal, posteromedial, medial PFC, and PHC). The PHC region from Network A was paired with the posteromedial region from Network B. All 20 ANOVAs were individually significant (p \< 0.01) with most showing crossover interactions ([Figure 2](#fig2){ref-type="fig"}). One exception was the medial PFC interaction that, while significant in all four subjects, demonstrated a clear crossover in two subjects and minimal difference for the most ventral region in two subjects. Another exception was the PHC interaction, which did not show a clear crossover interaction in S1. Note that, in most cases, the crossover interactions are present even where regions are extremely close to one another.

While the above analyses formally test the double dissociation, another source of evidence is that the spatial patterns replicate within individuals across discovery and replication (hypothesis-testing) datasets ([Figure S3](#mmc1){ref-type="supplementary-material"}). Additional analyses examined alternative methods for identifying these separate networks, including using data-driven clustering and estimation of connectivity patterns in the volume in addition to the surface ([Figure S4](#mmc1){ref-type="supplementary-material"}). Visualization in the native volume is important because projection to the cortical surface can induce fractionations of single regions into multiple regions if they fall near sulcal boundaries. The combined results illustrate a robust double dissociation between two distributed networks that fractionate the canonical DN across numerous cortical zones ([Figure 3](#fig3){ref-type="fig"}).Figure 3Parallel Distributed Networks Contain Juxtaposed Regions in Numerous Cortical ZonesThe two dissociated networks near the canonical default network, Network A and Network B, are shown for two subjects (S1 and S4) in a schematic form on the same cortical surface representation. The dashed boxes highlight nine cortical zones where neighboring representations of the two networks were found including: (1) dorsolateral PFC, (2) inferior PFC, (3) lateral temporal cortex, (4) inferior parietal lobule extending into the temporoparietal junction, (5) posteromedial cortex, (6) midcingulate cortex, (7) dorsomedial PFC, (8) ventromedial PFC, and (9) anteromedial PFC. Some zones, including the dorsal region along the PFC (labeled 7), are subtle, but consistent, in all subjects, suggesting that there exists small, closely juxtaposed components of the two dissociated networks.

The Importance of Examining Network Organization within the Individual {#sec2.4}
----------------------------------------------------------------------

The organization of the dissociated networks suggests why they might evade group-averaged analyses. Spatial "hurdles" to group averaging have been previously reported (e.g., [@bib23]'s Figure 1; [@bib55]'s Figure 7). To quantify the effect of spatial misalignment between individuals, we took an approach using spatial yoking between the individuals. For each individual, her spatially optimized regions from the discovery dataset were applied to her independent replication data, yielding an unbiased correlation matrix ([Figure S5](#mmc1){ref-type="supplementary-material"}). Two clear clusters of strong within-network and minimal between-network correlations were observed in each individual, with S1 and S4 showing the strongest patterns. This illustrates that widely distributed regions can show strong correlation with one another, while spatially adjacent regions can be embedded in distinct correlated clusters. The same correlation matrices showed minimal structure when regions from one person were applied to another (e.g., Subject 1's regions were used to generate a matrix using Subject 2's fMRI data), highlighting the importance of respecting the exact spatial details present within an individual.

Topography of Multiple Distinct Networks within the Individual {#sec2.5}
--------------------------------------------------------------

The above analyses establish evidence for two distinct neighboring networks that are likely components of what has been studied as the DN. This unexpected observation prompted us to explore how these two networks relate to additional large-scale networks, the frontoparietal control network (FPN) and the dATN. Two questions drove these analyses. First, are the other large-scale networks themselves fractionated and, if so, how are the newly detected networks organized? Second, taken as a group, do the multiple large-scale networks possess consistent spatial relations between networks?

For these analyses, all 24 data sessions were combined to provide best estimate maps for each participant. Seed regions were placed in the left lateral PFC and the frontal eye fields to identify the FPN and dATN, respectively. For the FPN, two seed regions were selected that revealed networks that maximized the following features: (1) the networks contained strong representations in intraparietal, inferior temporal, and dorsomedial PFC, and (2) the networks occupied neighboring, but distinct, regions in each cortical location. For the dATN, similar criteria were applied, but the distribution instead included superior parietal, occipitoparietal, and occipitotemporal components.

As with the DN, both the FPN and dATN were fractionated into two distinct, parallel networks within the individual, identifying a total of six networks with close spatial arrangements ([Figure 4](#fig4){ref-type="fig"}). The effects were clearest in S1 and S4 and suggestive in the other subjects. In each case, the networks both resembled the canonical group-average network but inhabited separate subregions ([Figure 5](#fig5){ref-type="fig"}).Figure 4Multiple Parallel Interdigitated Distributed Networks at or near the Canonical Frontoparietal Control and Dorsal Attention Networks Estimated by Functional Connectivity within IndividualsBest estimate maps (using all 24 sessions in each individual) of Networks A and B that fractionate the default network are illustrated (top). Maps from two subjects (S1 and S4) are displayed. The canonical Frontoparietal Network (middle) and Dorsal Attention Network (bottom) also each fractionate into two juxtaposed networks. Seed regions are illustrated by filled white circles.Figure 5Relationship of Parallel Interdigitated Networks to Canonical Networks from Group-Averaged DataEach row illustrates how the networks identified in two individuals (S1 and S4) correspond to the well-characterized topography of group-derived networks. The black border represents the outline of the canonical default, frontoparietal control, and dorsal attention networks (top to bottom) calculated using data from 1,000 subjects that were parcellated into seven networks (from [@bib113]). The correlation maps from each seed (white filled circle) are shown in color. Broadly, the networks can be seen to occupy separate, closely juxtaposed regions that fall within the canonical network borders in most cases. Exceptions can also be found, such as in the inferior frontal cortex in Default Network A and in the parietal lobe in Frontoparietal Control Network B, where the individual's connectivity map strays outside the group network borders.

Important organizational details were evident when all six networks were visualized simultaneously. [Figures 6](#fig6){ref-type="fig"} and [7](#fig7){ref-type="fig"} focus on the anterior midline and parietal and temporal lobes to illustrate two organizational features. First, regions from each network are often located in distinct locations in each cortical zone. The white lines in [Figures 6](#fig6){ref-type="fig"} and [7](#fig7){ref-type="fig"} serve as landmarks to highlight non-overlapping features of the networks. Second, the networks display a fine-scale interdigitation. For example, in the anterior midline, a broad ventral to dorsal progression is observed ([Figure 6](#fig6){ref-type="fig"}); however, particularly for the DN, the representations from Network B were positioned in between representations from Network A. Similarly in the temporal lobe, the representation from Network A is largely surrounded by representations from Network B ([Figure 7](#fig7){ref-type="fig"}). In the parietal lobe, a broad posteroventral to anterodorsal progression is observed across the networks, with each of the six networks inhabiting a distinct region ([Figure 6](#fig6){ref-type="fig"}). The fractionated Networks A and B of the canonical dATN showed three or more separate regions, as visualized on the surface, which sequentially alternated along the intraparietal arc of the canonical dATN ([Figures 4](#fig4){ref-type="fig"}, [5](#fig5){ref-type="fig"}, and [6](#fig6){ref-type="fig"}).Figure 6Detailed Anatomy of Six Distinct Networks: Parietal and Medial Prefrontal CorticesThe fine-scale interdigitation of the six identified networks is highlighted in regions where the macroscale organization is evident. White lines serve as landmarks so that the relative position of each network can be appreciated across panels: Networks A and B of the Default Network (DN-A, DN-B), Frontoparietal Control Network (FPN-A, FPN-B), and Dorsal Attention Network (dATN-A, dATN-B). In each row, FC maps from an individual are displayed for either medial frontal cortex (top two rows) or lateral parietal cortex (bottom two rows).Figure 7Detailed Anatomy of Six Distinct Networks: Lateral Temporal CortexIn a similar format to [Figure 6](#fig6){ref-type="fig"}, each column displays FC maps from an individual to illustrate the fine-scale interdigitation in the lateral temporal cortex.

[Figures 8](#fig8){ref-type="fig"} and [S6](#mmc1){ref-type="supplementary-material"} show a diagrammatic representation of the six different networks in two subjects to highlight the close, parallel nature of their organization within frontal, parietal, and temporal lobes. For an additional analysis, which should be considered descriptive, a correlation matrix was constructed ([Figure S7](#mmc1){ref-type="supplementary-material"}) using regions chosen from all six networks ([Figure S8](#mmc1){ref-type="supplementary-material"}). Data from all 24 sessions were used to construct each matrix. Given that the regions were defined and tested on the same data, the specific quantitative values of the within-network correlations should not be interpreted; however, the between-network correlations reflect an unbiased estimate of interactions between networks. Of interest, certain networks showed hints of interactions with other networks. For example, the FPN-A showed slightly elevated correlation with DN-A; the FPN-B showed slightly elevated correlation with dATN-A. These suggestive interactions may be due to spatial blurring or to biologically meaningful factors and are presented purely for their ability to generate future hypotheses.Figure 8Diagrammatic Representation of Six Parallel Distributed Networks within One IndividualThe central figure shows an illustration of the six networks overlaid on the same cortical surface. The top panel shows the lateral view, and the lower panel shows the medial view. The different colors correspond to the canonical network that each network resembles (red, Default Network, DN-A and DN-B; blue, Frontoparietal Network, FPN-A and FPN-B; green, Dorsal Attention Network, dATN-A and dATN-B). The names of the networks are based on prior literature, recognizing that the novel organization identified here may lead to a reconsideration of the functional domains. Data shown are from S4 (see also [Figure S6](#mmc1){ref-type="supplementary-material"}).

Discussion {#sec3}
==========

The present study examined the organization of large-scale distributed networks within the individual. We discovered that the canonical DN fractionates into two parallel networks that have juxtaposed regions throughout the cerebral cortex. Motivated by this discovery, we examined the dATN and FPN and found that each of these canonical networks also fractionates into two parallel networks. The organization of the six identified networks was charted and found to have a spatial progression in multiple zones of cortex. These results are consistent with the ideas that the large-scale networks (1) share a conserved motif and (2) are embedded within a broad macroscale organization. As a technical point, the present results underscore a need to move from group-based analyses to examination of detailed anatomy within the individual.

Canonical Networks Fractionate into Distinct Networks within the Individual {#sec3.1}
---------------------------------------------------------------------------

The most pressing finding reported here is that three major networks (DN, FPN, and dATN) are each subdivided into parallel spatially juxtaposed networks ([Figures 1](#fig1){ref-type="fig"}, [3](#fig3){ref-type="fig"}, [4](#fig4){ref-type="fig"}, and [5](#fig5){ref-type="fig"}). Regions of the separate networks lay side by side to one another across several cortical zones ([Figure 3](#fig3){ref-type="fig"}) and exhibited complex topography, with one network's region sometimes surrounded by another's ([Figures 6](#fig6){ref-type="fig"} and [7](#fig7){ref-type="fig"}). It is thus unsurprising that detailed within-subject analysis is needed to visualize the separate networks. We focused on the DN first to show that the within-subject network fractionation is reproducible across ([Figure 1](#fig1){ref-type="fig"}) and within individuals ([Figure S3](#mmc1){ref-type="supplementary-material"}) and is statistically robust in an independent sample ([Figure 2](#fig2){ref-type="fig"}). We next characterized fractionations of the FPN and the dATN ([Figures 4](#fig4){ref-type="fig"}, [5](#fig5){ref-type="fig"}, [6](#fig6){ref-type="fig"}, and [7](#fig7){ref-type="fig"}).

These findings raise the prospect that the canonical networks studied in group-averaged data consist of distinct functional networks that are blurred together by spatial averaging. The six networks identified here appear to be fractionations of networks identified in group studies (e.g., [@bib113], [@bib82], [@bib20]). This point is important to emphasize: it is not the case that more detailed analysis carved network organization orthogonally to prior schemes, but rather fractionated the existing lower-resolution frameworks ([Figure 5](#fig5){ref-type="fig"}). It remains to be determined whether these different networks mediate separable functions across different task contexts. Such a finding could help explain the heterogeneity in cognitive functions and clinical conditions ascribed to the canonical networks (e.g., [@bib14], [@bib98], [@bib68], [@bib54], [@bib57], [@bib6]).

Relations to Prior Observations {#sec3.2}
-------------------------------

There is growing consensus that idiosyncratic features exist within individuals that are not captured (or are attenuated) by examining group central tendencies (e.g., [@bib22], [@bib23], [@bib75], [@bib55], [@bib37], [@bib48], [@bib49], [@bib39], [@bib40]). The DN was originally hypothesized based on the distributed pattern of regions that increase activity in passive relative to active non-self-referential tasks ([@bib4], [@bib93], [@bib66]; for review, see [@bib44]). The DN was later estimated using FC by placing a moderately sized seed region in the center of posteromedial cortex (often labeled posterior cingulate cortex or "PCC") and plotting the correlation pattern ([@bib41]). One immediate feature of the resulting network that raises the possibility of further subdivision is its large size. For example, the canonical group-averaged FC estimate of the DN contains an extensive correlation pattern extending from the dorsal extent of the frontal midline to ventromedial and orbitomedial PFC (e.g., [@bib28], [@bib14]).

A number of group-based studies, including our own, have sought to fractionate the DN into subnetworks (e.g., [@bib61], [@bib62], [@bib14], [@bib32], [@bib5], [@bib6], [@bib58], [@bib59]). While the present fractionation is not contained within these past efforts, prior observations have noted juxtaposed nodes in different zones of cortex. There is evidence for a dorsal to ventral separation in the posterior midline ([@bib62], [@bib58]; see also [@bib110]) that likely corresponds to the present DN-A and DN-B networks, although the topography is not fully captured by a simple linear axis in all individuals (e.g., see S2 in [Figure 1](#fig1){ref-type="fig"}). Prior reports have also noted that a subnetwork of the DN is coupled to the hippocampal formation (e.g., [@bib5]), consistent with anatomical connectivity in the macaque (see [@bib52] for review). The present hypothesized DN-A aligns well to the "hippocampal" subnetwork of [@bib5], while DN-B does not. This separation also informs our understanding of the IPL.

While many analyses of the DN, especially those based on increased response during passive tasks, reveal a large region covering much of IPL, social tasks involving mentalizing (theory-of-mind tasks) activate the more anterior temporoparietal junction (TPJ; [@bib88]) while tasks involving episodic remembering activate a caudal region (e.g., [@bib6]) near what might be a monkey homolog of Opt ([@bib80]; for discussion, see [@bib113]). Group-level FC analyses have noted that the posterior portion of IPL is preferentially coupled to the PHC, while the anterior portion is not (see [@bib113]'s Figure 30). Furthermore, monkey anatomical tracing studies consistently show that PHC projects to a circumscribed portion of area 7A within Opt (e.g., [@bib56], Case M-2-90; [@bib10], Cases 1 and 5). These collective findings are consistent with the distinction between DN-A, which is coupled to the PHC and posterior IPL, and DN-B, which involves a more anterior IPL region. What is novel in the present work is that this distinction is now shown to be one spatial component of a much broader separation of two parallel large-scale distributed networks.

We know of no precedent for one feature of our results. The ventral portion of the frontal midline is considered a projection zone of limbic structures, including the amygdala ([@bib79]) and hippocampal formation ([@bib87]). However, in the present study, the most ventral representation of DN-B is inferior to DN-A. This is unexpected because DN-B is differentiated from DN-A by its absence of coupling to the PHC and retrosplenial cortex ([Figure S2](#mmc1){ref-type="supplementary-material"}). This may indicate that a subregion of ventromedial PFC is tied to a large-scale distributed network that is functionally separated from a direct limbic influence. The fine interdigitation shown in [Figure 6](#fig6){ref-type="fig"} reveals why group averaging will likely blur the two networks. In each individual, there are multiple regions for each network that are interposed and whose positions spatially shift between individuals. Anatomical tracing in the monkey will be needed to substantiate whether there exists a ventral midline region that is minimally connected to limbic structures.

Parallel Large-Scale Distributed Networks Are an Organizing Principle of Association Cortex {#sec3.3}
-------------------------------------------------------------------------------------------

A notable feature of the identified networks is that each contains components in frontal, parietal, temporal, and frontal midline regions. This repeating pattern or motif has been discussed previously in group data (e.g., [@bib113], [@bib82]). What is striking is how well the same distributed pattern accounts for much of the newly fractionated networks.

[@bib38] suggested that this distributed motif was a general organizing principle of association circuits. Based on results from double-labeling tracer injections by [@bib89], [@bib90], she posited (1) that prefrontal and parietal areas are embedded within densely interconnected distributed circuits that include midline and temporal areas and (2) that this same motif repeats across nearby zones forming closely adjacent parallel networks. Our results are consistent with her ideas.

One interesting finding is that there can be clear distinctions even between networks with closely juxtaposed regions. As noted, DN-B does not couple to the hippocampal formation, while DN-A does. Similarly, dATN-B shows coupling to retinotopic visual regions along the midline, preferentially to the peripheral field representation, while dATN-A does not (see [Figure S8](#mmc1){ref-type="supplementary-material"}). The presence of such differences may shed insight into the functional role of the parallel networks and their development. One possibility is that there are broad constraints that establish the same motif but, seeded by competing inputs from limbic and sensory systems, activity-dependent processes differentiate the networks during development.

A second important finding is the fine spatial scale that differentiates neighboring networks. The spatial interdigitation of distinct regions was notable. In the lateral temporal lobe, one network's region could be surrounded by another's. This fine-scale interdigitation has implications for interpreting prior network estimations. For example, [@bib69] proposed a cortical network important to spatial attention. Extensive findings illustrate that certain areas near to the intraparietal sulcus form part of a sensory-motor hierarchy, sharing reciprocal projections with extrastriate visual cortex and the frontal eye fields ([@bib65], [@bib101], [@bib3], [@bib11]). However, Mesulam's ideas drew on injections within the IPL area 7a in or near Opt (e.g., [@bib72]). This specific parietal association area is near to macaque area LIP but has a projection fingerprint that spares distant extrastriate areas while including cingulate and PHC ([@bib3]). Given the present results, it is reasonable to suppose that past analyses of parietal association cortex may have lumped together injections in separate parallel, distributed networks.

Evidence for a Macroscale Organization of Association Cortex that Spans Networks {#sec3.4}
--------------------------------------------------------------------------------

An intriguing finding is revealed when the spatial relations across all the networks are considered together ([Figures 8](#fig8){ref-type="fig"} and [S6](#mmc1){ref-type="supplementary-material"}). While the networks have complex interdigitated relationships, there are also macroscale gradients that share the same general progression in multiple zones of cortex. In parietal association cortex, there is a caudal to rostral progression from DN-A through to dATN-B ([Figure 6](#fig6){ref-type="fig"}). In temporal association cortex, there is a rostral to caudal progression through the same networks ([Figure 7](#fig7){ref-type="fig"}). The axis is imperfect, and it is sometimes unclear which network should be ordered before the other, but the general ordering that repeats across distinct zones suggests a broad macroscale organizing principle.

[@bib63] recently argued that association cortex possesses a macroscale gradient of networks from sensory-motor networks on one end to the DN on the other. In agreement, the possibility was recently raised that association networks evolved from a prototypical distributed sensory-motor network followed by a period of cortical expansion, which freed up zones of association cortex from the constraints of primary sensory input ([@bib13], [@bib53]). The parallel and sequential nature of the presently defined networks adds further support to these ideas.

Limitations and Technical Considerations {#sec3.5}
----------------------------------------

An assumption behind our interpretation of the results is that FC across distributed regions is sufficiently constrained by direct and polysynaptic anatomical circuits to provide insight into the organization of distributed networks. Parallels with macaque anatomy reinforce this assumption. However, details of the results may be revised to the degree that factors beyond stable anatomical constraints contribute to the patterns. Exploration of monkey anatomy using multiple tracer injections from adjacent regions in the same animal would be a valuable complement to the present work. The anatomical origin of the fine spatial details are important to resolve because of their implications for clinical endeavors, including presurgical planning and targeted neuromodulation. The present results suggest that targeted neuromodulation of nearby cortical zones could have distinct effects because they are embedded within anatomically separate networks. While this general notion has been appreciated previously (e.g., [@bib31]), what is surprising is the fine spatial scale by which cortical zones participating in distinct networks are interdigitated.

Given that we were able to fractionate established large-scale networks by pushing the practical resolution of fMRI by targeting the individual, it seems likely that our present estimates might also be fractionated further if higher resolution was achieved. Despite efforts to minimize spatial blurring, we were unable to confidently delineate all networks in all individuals. Two subjects, S2 and S3, produced maps that were noticeably blurrier. This difference may be due to several factors, such as differences in head motion, SNR, and misregistration. The reported summary measures do not clearly indicate a cause. The fractionation of the DN into two parallel, distributed networks proved to be the most robust finding observed clearly in all subjects.

The present study is also limited in that we studied only cerebral cortical organization. Ongoing work is characterizing the topography of these networks in subcortical and cerebellar structures.

Conclusions {#sec3.6}
-----------

The present work reveals that there are parallel large-scale distributed networks that are spatially juxtaposed across the cerebral cortex. The spatial scale of these networks is such that they become evident only when analyzed within the individual. Discovery of the presence and description of details of these networks provide a foundation for future study of their functions.

STAR★Methods {#sec4}
============

Key Resources Table {#sec4.1}
-------------------

REAGENT or RESOURCESOURCEIDENTIFIER**Software and Algorithms**FreeSurfer[@bib25]<http://surfer.nmr.mgh.harvard.edu>MATLABMathWorks[www.mathworks.com](http://www.mathworks.com){#intref0015}FSL[@bib95]<https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/>Connectome Workbench[@bib60][www.humanconnectome.org](http://www.humanconnectome.org){#intref0025}

Contact for Reagent and Resource Sharing {#sec4.2}
----------------------------------------

Further information and requests for resources should be directed to and will be fulfilled by the Lead Contact, Rodrigo M. Braga (<rodrigo.braga@imperial.ac.uk> or <rbraga@fas.harvard.edu>).

Experimental Model and Subject Details {#sec4.3}
--------------------------------------

### Participants {#sec4.3.1}

Four healthy female right-handed young adults (ages 21 to 26) were recruited from the greater Boston community for a study that involved 24 separate MRI scanning sessions as well as extensive behavioral monitoring over a period of approximately 16 weeks. None of the participants were local (Harvard) students or institutional employees, and all were paid for participation with milestone payments at the end of the trial period (two weeks), after 12 MRI sessions, and after all 24 MRI sessions. Participants were screened to exclude a history of neurological or psychiatric illness, or ongoing use of psychoactive medications. Seven participants were enrolled and three chose not to continue within the first two weeks (which were described to all participants as a trial period). Four participants were enrolled for the full 24-session extended study and all of these individuals completed all intended sessions. The present paper concerns only the fMRI data, but enrolled individuals also participated in extensive behavioral testing and daily monitoring of behavior via smart phones (Beiwe; [@bib100]), sleep and activity monitoring via an Actigraph 2 wrist wearable (Philips Respironics, Murrysville, PA, USA), as well as hearing tests at several times during the study (Model 2500 Microprocessor Audiometer, AMBCO, Tustin, CA, USA). One subject required vision correction using MRI compatible glasses. Participants provided written informed consent in accordance with the guidelines set by the Institutional Review Board of Harvard University.

Method Details {#sec4.4}
--------------

### MRI Data Acquisition {#sec4.4.1}

Data were collected on a 3T Siemens Prisma-fit MRI scanner (Siemens Healthcare, Erlangen, Germany) using the vendor's 64-channel phased-array head-neck coil. Heads were immobilized with Siemens small foam head coil wedges. Each of the 24 MRI sessions included one 7 m 2 s run of resting state data (passive fixation) to estimate intrinsic functional connectivity ([@bib9]) as well as a number of other acquisitions (structural, ASL, and task-based functional runs). The resting state run was collected in the same fixed order during every session near to the beginning of the session to optimize compliance. Participants were instructed to remain still, stay awake and to fixate a centrally presented crosshair presented in black on a light gray background. The gray background color was used instead of white to reduce glare, eye fatigue and discomfort. The position of the screen was adjusted at the beginning of each session to ensure a comfortable, central viewing position to minimize muscle tension and head motion. Before each session participants were encouraged to spend a few minutes finding a comfortable lying position that they could maintain for the entire session. The scanner room lights were kept on to deter participants from becoming drowsy.

Eye closures and movements were monitored using the Eyelink 1000 Core Plus with Long-Range Mount. A video of the eye tracker output was recorded in order to quantify compliance and arousal. Additional in-scanner physiological monitoring (Biopac Systems Inc, Goleta, CA, USA) included a respirator belt around the chest to monitor breathing (Biopac, TSD221-MRI), electrodermal electrodes (Biopac, EL508) attached to participants' right sole to measure galvanic skin response, and a band pulse-oximeter (Biopac, OXY-MRI-SENSOR) attached to participants' right middle toe to measure oxygen saturation and pulse rate. See [@bib111] for details. The oximeter and electrodes were place on participants' feet to keep their hands free to make responses using button boxes during the in-scanner tasks. Immediately before each run, participants were asked to remain still for the entirety of the upcoming run. Following each run, participants were given feedback about noted movements and encouraged to stay still, in order to establish an expectation that their compliance and movement were being carefully watched.

Functional imaging data were acquired using a multi-band gradient-echo echo-planar pulse sequence ([@bib91]) with acquisition parameters: TR 1000 ms; TE 32.6 ms; flip angle 64°; 2.4 mm isotropic voxels; FOV 211 mm x 211 mm x 156 mm; 65 slices fully covering the cerebral cortex and cerebellum. Slice acquisition used interleaved simultaneous multi-slice 5x acceleration. The sequence was a custom sequence generously provided by the Center for Magnetic Resonance Research (CMRR) at University of Minnesota. Whole brain coverage and minimization of signal dropout due to magnetic susceptibility were achieved by aligning slices to a plane 25 degrees from the anterior commissure-posterior commissure plane toward the coronal plane ([@bib112], [@bib67]). This was implemented using an automated alignment procedure to ensure consistency across sessions ([@bib102]) and, in pilot acquisitions, was found to increase signal-to-noise in ventromedial prefrontal cortex (PFC). A rapid T1-weighted structural image was also acquired in each session using a multi-echo MPRAGE three-dimensional sequence ([@bib103]) with acquisition parameters: TR 2200 ms; TE 1.57, 3.39, 5.21, 7.03 ms; flip angle 7°; 1.2mm isotropic voxels; 144 slices; FOV 230 mm x 230 mm x 173 mm, in-plane GRAPPA acceleration 4 (see [@bib47] for empirical results and discussion of comparability of this brief sequence to traditional longer acquisitions).

### Data Preprocessing {#sec4.4.2}

Resting-state data were processed using methods similar to those previously described ([@bib104], [@bib113]): (1) 12 initial volumes from each run were discarded to allow for T1-equillibration, (2) head motion was corrected using rigid body translation and rotation (FSL, <https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/>; [@bib50], [@bib95]), (3) data were temporally low-pass filtered at a threshold of 0.08 Hz, and (4) nuisance variables (6 motion parameters, mean whole-brain signal, mean ventricular signal, mean deep cerebral white matter signal) and their temporal derivatives were regressed. Structural data were processed using the FreeSurfer version 4.5.0 software package (<http://surfer.nmr.mgh.harvard.edu>; [@bib25]). For each anatomical image (one per session), a surface mesh representation of the cortex was reconstructed and registered to a common spherical coordinate system by aligning the major sulcal patterns to the FreeSurfer average template ([@bib26]). The preprocessed functional images from each session were aligned to the cortical surface mesh reconstructed from that session's anatomical image using boundary-based registration ([@bib43]). Functional data were then propagated to the common spherical coordinate system via sampling (trilinear interpolation) from the middle of the cortical ribbon in a single interpolation step.

Functional data were sampled to the fsaverage6 surface mesh ([@bib26]) containing 40,962 vertices per hemisphere, and a 2mm full-width-at-half-maximum (FWHM) smoothing kernel was applied to the data in the surface space. A mesh resolution of 40,962 vertices was chosen to reduce blurring during the trilinear interpolation step and hence maximize the potential for observing network distinctions, while keeping the computational burden manageable. A bespoke cortical surface template containing 40,962 vertices per hemisphere was produced using the Connectome Workbench's command suite ([@bib36]). This was done so that the functional connectivity analyses could be performed and visualized interactively within the Workbench's flexible surface-based visualization software, wb_view ([@bib60]). The bespoke template was created by combining the left and right pial surfaces from the fsaverage6 freesurfer template into the CIFTI format using the Workbench commands. The pial surfaces were then selectively inflated (smoothing cycles: 3, smoothing strength 0.7, smoothing-iterations 13, inflation factor 1.02) using Workbench to allow visualization of the major cortical folds while maintaining the majority of the cortical surface visible.

Voxel-based tSNR maps were computed by taking the motion-corrected time series from each functional run and dividing the mean signal at each voxel by its standard deviation over time. The tSNR maps were then averaged across functional runs within the Discovery (n = 12), Replication (n = 12) and Full datasets (n = 24), and the resulting mean tSNR maps were projected to the cortical surface for visualization using FreeSurfer ([Figure S1](#mmc1){ref-type="supplementary-material"}). An additional metric of data quality, fractional Amplitude of Low Frequency Fluctuations (fALFF), was also computed ([Figure S1](#mmc1){ref-type="supplementary-material"}). fALFF maps were produced by normalizing the total power in the low (0.01 -- 0.08 Hz) frequency range by the total power across all frequencies ([@bib115]).

### Discovery of Networks within the Individual {#sec4.4.3}

For each participant, half of her data were used in a discovery manner to identify networks that would be later tested in the remaining independently collected data sessions. The odd-numbered sessions (i.e., 1^st^, 3^rd^, 5^th^, etc) formed the discovery dataset (n = 12) while the even-numbered sessions (i.e., 2^nd^, 4^th^, 6^th^, etc) were set aside as the hypothesis-testing dataset (next section). The discovery datasets in each of the four participants were analyzed blind to the hypothesis-testing datasets.

For the discovery analysis, Pearson's product moment correlations between the fMRI time series at each cortical surface vertex were computed. This resulted in an 81,924 × 81,924 element cross-correlation matrix (40,962 vertices per hemisphere) for each of the 12 fMRI runs from the discovery dataset. The matrices were then r-to-z transformed and averaged to yield a mean matrix with high stability. These discovery matrices were used to explore network organization. The mean cross correlation matrices were assigned to the bespoke cortical template so that individual seed vertices could be selected and their functional connectivity maps interactively viewed using wb_view ([@bib60]). Individual vertices were selected from the general vicinity of expected locations of target networks in PFC (as estimated from independent group-averaged data; n = 1000 from [@bib113]) and expectations from the literature (e.g., [@bib82], [@bib113]).

First, a seed vertex was selected from lateral PFC near regions which form part of the canonical default network, and the resulting maps visualized. If the functional connectivity map revealed a network of distributed regions showing a robust correlation with the seed region (z(r) ≈0.6), the seed vertex number was recorded for further analysis. If the map revealed no strong correlations, or a diffuse network of correlations that were observably lower (z(r) ≈0.4), this was taken as evidence that the seed region was sampling an area of signal dropout or containing a mixture of signals, respectively, and a different seed region was selected. The colorbar scale was set between 0.2 and 0.6, using the JET256 palette in wb_view, so that the correlation structure could be adequately represented and these subtle differences observed. After a robust network was detected and recorded, a different seed vertex was selected that satisfied 4 criteria: the second seed region i) was within the lateral PFC, ii) was in the vicinity of the previously identified network's seed, iii) was in a region of the cortical surface that showed low correlation (z(r)\<0.3) with the previous network's seed region, and iv) also showed robust correlation with a distributed set of regions. Thus the analysis converged on robust networks that had closely neighboring representations within the lateral PFC.

The goal of this discovery procedure was not to confirm expectations from the group data, but to allow the individual participant maps to be extensively interrogated, moving outward from properties expected from the group maps. Given prior analyses of within-subject data (e.g., Figure 4 from [@bib19]; Figure 3 from [@bib107]; Figure 9 from [@bib104]), it was unsurprising that many maps constructed within the individual participants resembled canonical networks discovered in group-analyzed data. The targets of our exploration were network features that are not fully captured by group analyses (e.g., [@bib55], [@bib39], [@bib40]; see also [@bib23]).

As the results will demonstrate, seed regions placed in nearby regions of lateral PFC revealed two important features that led to the regions selected for the hypothesis-testing phase of analysis. The first feature was that, like typical group-based analyses, the resulting distributed networks contained inferior parietal, temporal, medial prefrontal, and posteromedial cortical components near to what has been described as the 'default network' (DN). Second, nearby seed regions yielded distributed networks that were closely neighboring but separate throughout much of the distributed organization of the network. That is, two distinct networks were closely interdigitated throughout the topography of the canonical DN suggesting the hypothesis that the DN may be comprised of multiple neighboring networks. We refer to these hypothesized networks as Network A and Network B ([Figure 1](#fig1){ref-type="fig"}).

### Hypothesis Testing to Dissociate Networks within the Individual {#sec4.4.4}

The discovery phase of data analysis led to the hypothesis that two distinct interdigitated networks exist that each are within or near to the canonical DN ([Figure 1](#fig1){ref-type="fig"}). The two networks were present in each individual participant. The goal of the hypothesis-testing phase was to use independent data in each participant to support or refute the possibility of two dissociable networks. To conduct this analysis, two distinct neighboring lateral PFC seed regions were selected within each participant as well as pairs of regions throughout the cortex based only on the discovery datasets that maximized the separation of the distributed networks. A priori regions (single vertices) were selected in each of the main regional zones of the cortex (temporal, inferior parietal, posteromedial, and medial prefrontal). These a priori regions were targeted to locations where contiguous vertices could be observed in Network A and Network B ([Figure 1](#fig1){ref-type="fig"}). A seed region was also selected in parahippocampal cortex (PHC) to quantify the representation of Network A in this region, given a prior literature linking the hippocampal formation to the DN ([@bib42], [@bib107], [@bib52]). These regions were then statistically tested in the independent data to formally dissociate the two networks.

The critical test was whether there would be significant interactions between the two lateral PFC seed regions and Network A and Network B regions in each of the zones of cortex in the independent hypothesis-testing datasets. The presence of significant interactions would be evidence for regional dissociations. The presence of interactions across all distributed zones of the cortex would be strong evidence that there was a complete network dissociation across the cortex, in essence establishing dissociable but adjacent networks. For each of the 12 hypothesis-testing data sessions within each participant, values representing the r-to-z transformed Pearson's product moment correlation for each of the two lateral PFC seed regions were computed in relation to each of the distributed test regions. Statistical tests were performed as a two-way ANOVA. In each ANOVA, the independent (classification) variables were the PFC seed location (A and B) and the a priori test regions (A and B) within one of the separate zones of cortex. The dependent variable was the r-to-z transformed Pearson's product moment correlation between each seed and test region. Each cell contained a correlation measure from each of the 12 sessions in the hypothesis-testing dataset, producing a balanced 2 X 2 factorial design with 12 elements in each cell. This 2-factor analysis was repeated for each of the separate cortical zones (temporal, inferior parietal, posteromedial, and medial prefrontal). Statistical significance level was set at p \< 0.01. Critically, the multiple tests across the network and across subjects were not independent in the sense that they were targeting convergent evidence for network dissociation. Thus multiple, repeated significant results in the ANOVA across regions and participants would provide a high level of certainty for dissociation. Sporadic significance that occurred in 1 in 20 or 1 in 100 tests that showed no specific pattern might be indicative of false positives. As the results will reveal, the data patterns and hypothesis-directed statistical tests were quite clear in the weight of their evidence.

### Effects of Misalignment between Individuals {#sec4.4.5}

The details of network organization that are revealed by our analyses within individuals suggest a level of spatial and anatomical specificity that would likely be lost or underappreciated when central tendencies across participants are extracted in group-averaged data or when details of anatomy in one participant are assumed to apply to another. That intuition can be appreciated visually by examining [Figures 1](#fig1){ref-type="fig"}, [3](#fig3){ref-type="fig"}, [4](#fig4){ref-type="fig"}, [5](#fig5){ref-type="fig"}, [6](#fig6){ref-type="fig"}, and [7](#fig7){ref-type="fig"}, and also in results from other papers (e.g., Figure 7 from [@bib55]). To formally explore this issue, we asked to what degree network structure present using individually-tailored regions is lost when one person's anatomy is assumed to apply to another person.

For these analyses, in each of the four participants we extracted the hypothesized cortical regions that were components of DN-A and DN-B from the discovery datasets (lateral prefrontal, temporal, inferior parietal, posteromedial, and medial prefrontal, but not the parahippocampal region; [Figure 1](#fig1){ref-type="fig"}), yielding 10 individually-tailored regions per subject. These regions were then used to construct a 10 × 10 matrix of correlations in the hypothesis-testing dataset for each participant. To address the question of whether misalignment between individuals affects the results, we recalculated the matrices in each participant using the regions defined in the other three participants ([Figure S5](#mmc1){ref-type="supplementary-material"}). The goal of this analysis is to provide some level of visualization of what is maintained and what is lost when misalignment is present across individuals. Our focus here is on the specific dissociation between DN-A and DN-B as an example of an important feature of functional-anatomical organization (see [@bib39] for a conceptually similar analysis performed for algorithmically defined individually-tailored patches). As can be seen in [Figure S5](#mmc1){ref-type="supplementary-material"}, the clustering of individually tailored regions into the two distinct networks breaks down when seed region locations from a different subject are used. This result illustrates that the spatial variation between individuals is enough to obscure the clear and reproducible network dissociation that is uncovered within the individuals.

### Additional Explorations Targeting Adjacent Networks within the Individual {#sec4.4.6}

The main analyses explored and tested for dissociation of two distinct neighboring networks that fractionate the canonical DN. This was an unexpected result and encouraged further exploration to determine whether there were further network fractionations. For these additional, post hoc explorations the full dataset for each participant was used to dissociate networks. Results are reported for network properties that were evident in at least two separate participants.

### Confirmation of Interdigitated Network Representations in the Brain Volume {#sec4.4.7}

Surface projection provides a convenient means to visualize cortical organization and the juxtaposition of networks on the continuous surface. However, projection to the surface involves a spatially-nonlinear transformation that can project nearby voxels in the volume to distant positions on the surface (such as occurs when folds near the crowns of separate gyri abut one another). It is thus important to check that the observed interdigitation between the dissociated networks is present in both surface and volume representations and is not the result of the cortical sampling procedure used. Moreover, volume-based analysis will be required for many applied endeavors including presurgical planning and localization for neuromodulation (e.g., transcranial magnetic and direct-current stimulation).

To explore these issues, we reproduced the findings of two dissociated networks linked to the default network in the volume. The functional data from each subject were preprocessed as described for the surface-based pipeline, with the exceptions that smoothing was performed in the volume at 2mm FWHM prior to regression of nuisance variables and bandpass filtering, and that data were not projected to the surface. The first 8 resting state runs from each subject were concatenated in time and AFNI's InstaCorr (R. Cox and Z. Saad, 2010, International Conference on Resting-State Functional Brain Connectivity, conference; <https://afni.nimh.nih.gov/afni>) software was used to interactively select seed voxels and view the resulting FC maps. Two seed voxels were selected in left lateral PFC that maximized the separation between the two FC map representations in posteromedial cortex, namely a ventral representation in retrosplenial cortex for DN-A, and the dorsal posterior cingulate cortex representation for DN-B. When the two candidate seed regions were identified, the seed voxel locations were recorded and FC maps were produced in all 24 runs of the data for each subject.

Pearson's product moment correlations were calculated between the fMRI time series at each voxel and the seed voxel. The maps were then r-to-z transformed and averaged to produce a single FC map for each of the two seed regions. The question was whether the two parallel networks would retain the critical features observed in the cortical surface when viewed in the volume ([Figure S4](#mmc1){ref-type="supplementary-material"}). More broadly, while the surface-based visualization proved most effective for discovering organizational details, volume-based analysis is important to fully understand the underlying organizational features, including the possibilities that complex geometry of the sulcal patterns and also non-neuronal structures such as vessels may contribute to estimated features of cortical organization.

### Confirmation of Observed Parallel Networks Using Data-Driven Network Parcellation {#sec4.4.8}

The network dissociations reported were discovered and replicated using seed-based methods within-subjects that generalized across subjects. If the network dissociations are robust, they should also be observable using multiple approaches. That is, the path we initially employed to uncover them should not be the only way to reveal their presence and their specific spatial organization. To investigate whether the fractionation of the default network generalized across approaches, we explored network organization within the individual participants using data-driven clustering techniques.

We concatenated the surface-projected timeseries data from the Discovery dataset (n = 12), and used MATLAB's *kmeans* function (version R2015a; MathWorks, Natick, MA) to parcellate the vertex timeseries into clusters. As discussed extensively in our earlier analysis of group-based network estimates ([@bib113]), clustering approaches yield multiple solutions to data parcellations at many levels of clustering. Here we explored a relatively low dimensional fractionation (k = 12) that was able to yield the key dissociated networks including many of their distributed subcomponents.

### Additional Exploration of Network Organization along Regions of the Medial Temporal Lobe with Low Signal-To-Noise Ratio {#sec4.4.9}

As the results unfolded, the DN-A was found to possess a clear representation in the PHC, while no such representation was detected for DN-B. As portions of the medial temporal lobe are susceptible to signal drop-out in fMRI ([@bib78]), additional analyses were conducted to try to find evidence for a representation of DN-B in this region. Additional analyses were focused on the subject that showed the most robust network dissociations (S4). First, the FC maps for DN-A and DN-B from the PFC seeds were disattenuated by using the reliability of the functional connectivity maps across runs as an estimate of the signal dropout ([@bib76]; [Figure S2](#mmc1){ref-type="supplementary-material"}A). This revealed two representations belonging to DN-B in the inferior temporal lobe along zones of susceptibility artifact, but in regions located well outside the PHC. Second, seed vertices were selected in a posterior to anterior progression along the medial temporal lobe, and the resulting FC maps were observed ([Figure S2](#mmc1){ref-type="supplementary-material"}B). This analysis also did not reveal evidence for a representation of DN-B in the vicinity of the representation of DN-A. We also further explored whether other networks might be detected along the inferior and anterior temporal lobe, once signal dropout was partially accounted for. We describe these additional observations as potential avenues for further more detailed investigations as well as a reminder that, even with our numerous steps to increase signal-to-noise through signal averaging and our use of small voxels and acceleration during acquisition, signal loss due to susceptibility artifacts is still a problem in certain zones of cortex.

### Test-Retest Reliability of Parallel Interdigitated Networks {#sec4.4.10}

The FC maps for DN-A and DN-B were produced from the Discovery and Replication datasets ([Figure S3](#mmc1){ref-type="supplementary-material"}). Additionally, the mean connectivity matrix from the Discovery (n = 12) and Replication (n = 12) datasets were correlated to show the consistency of the connectivity patterns within each subject (S1: r = 0.88; S2: r = 0.95; S3: r = 0.87; S4: r = 0.94, all p \< 0.001).

Quantification and Statistical Analysis {#sec4.5}
---------------------------------------

This study includes n = 4 participants, each of which were scanned over 24 fMRI sessions. Each participant's imaging data were divided into discovery (n = 12; odd-numbered runs) and replication (n = 12; even-numbered runs) samples to allow for data exploration and statistical testing using independent data points. Functional connectivity between brain regions was calculated in MATLAB (version 7.4; <http://www.mathworks.com>; MathWorks, Natick, MA) using Pearson's product moment correlations which were r to z transformed prior to averaging or comparison. Statistical tests were performed as a two-way ANOVA using MATLAB's *anova2* function (version R2015a). Statistical significance was set to p \< 0.01. Network parcellation was performed using MATLAB's *kmeans* function (version R2015a).

Author Contributions {#sec5}
====================

R.M.B and R.L.B. designed the study, analyzed the data, interpreted the experiments, and wrote the paper. R.M.B. collected the data.

Supplemental Information {#app2}
========================

Document S1. Figures S1--S8Document S2. Article plus Supplemental Information

We thank the Harvard Center for Brain Science neuroimaging core and FAS Division of Research Computing. R. Mair and S. McMains assisted in optimizing acquisition. L. Farfel, M. Marotta, and R.M. Hutchison assisted in data acquisition. R.M.H. and T.M. O'Keefe assisted in data preprocessing. R.M.B. was supported by Wellcome Trust grant 103980/Z/14/Z. This work was supported by Kent and Liz Dauten, NIH grant P50MH106435, and Shared Instrumentation Grant S10OD020039. The multi-band EPI sequence was generously provided by the Center for Magnetic Resonance Research (CMRR) at University of Minnesota.

Supplemental Information includes eight figures and can be found with this article online at [http://dx.doi.org/10.1016/j.neuron.2017.06.038](10.1016/j.neuron.2017.06.038){#intref0060}.

[^1]: Lead Contact
