Patient lesion and neuroimaging studies have identified a rostral-to-caudal functional gradient in the lateral frontal cortex (LFC) corresponding to higher-order (complex or abstract) to lower-order (simple or concrete) cognitive control. At the same time, monkey anatomical and human functional connectivity studies show that frontal regions are reciprocally connected with parietal and temporal regions, forming parallel and distributed association networks. Here, we investigated the link between the functional gradient of LFC regions observed during control tasks and the parallel, distributed organization of association networks. Whole-brain fMRI task activity corresponding to four orders of hierarchical control [Badre, D., & D'Esposito, M. Functional magnetic resonance imaging evidence for a hierarchical organization of the prefrontal cortex. Journal of Cognitive Neuroscience, 19, 2082–2099, 2007] was compared with a resting-state functional connectivity MRI estimate of cortical networks [Yeo, B. T., Krienen, F. M., Sepulcre, J., Sabuncu, M. R., Lashkari, D., Hollinshead, M., et al. The organization of the human cerebral cortex estimated by intrinsic functional connectivity. Journal of Neurophysiology, 106, 1125–1165, 2011]. Critically, at each order of control, activity in the LFC and parietal cortex overlapped onto a common association network that differed between orders. These results are consistent with a functional organization based on separable association networks that are recruited during hierarchical control. Furthermore, corticostriatal functional connectivity MRI showed that, consistent with their participation in functional networks, rostral-to-caudal LFC and caudal-to-rostral parietal regions had similar, order-specific corticostriatal connectivity that agreed with a striatal gating model of hierarchical rule use. Our results indicate that hierarchical cognitive control is subserved by parallel and distributed association networks, together forming multiple localized functional gradients in different parts of association cortex. As such, association networks, while connectionally organized in parallel, may be functionally organized in a hierarchy via dynamic interaction with the striatum.
The lateral frontal cortex (LFC) has an established role in cognitive control function. However, the functional organization of the frontal cortex remains widely debated (Fedorenko, Duncan, & Kanwisher, 2013; Curtis & D'Esposito, 2004; Petrides & Pandya, 1994; Goldman-Rakic, 1988). One hypothesis proposes that the LFC is organized hierarchically to support behavior according to abstract, complex, or multistep rules (see Badre & D'Esposito, 2009; Badre, 2008). This hypothesis has been motivated in part by the frequent observation of an apparent functional gradient from caudal to rostral LFC in relation to progressively more abstract and complex tasks (Badre & D'Esposito, 2007, 2009; Christoff, Keramatian, Gordon, Smith, & Mädler, 2009; Badre, 2008; Buckner, 2003; Koechlin, Ody, & Kouneiher, 2003). Though there have been some notable exceptions to this pattern observed in some experiments (Crittenden & Duncan, 2014; Reynolds, O'Reilly, Cohen, & Braver, 2012), rostrocaudal functional distinctions have been shown to be replicable across a variety of experiments that have controlled for variables like overall difficulty and individual differences in the control structure of the task (Badre & Nee, 2018; Nee & D'Esposito, 2016; Bahlmann, Blumenfeld, & D'Esposito, 2015; Dixon, Fox, & Christoff, 2014; Nee, Jahn, & Brown, 2014; Nee & Brown, 2012).
Beyond the observation of a functional gradient, there has been evidence from fMRI effective connectivity (Nee & D'Esposito, 2016; Koechlin et al., 2003) and studies of patients with focal neurological damage (Azuar et al., 2014; Badre & D'Esposito, 2009) that separate regions of LFC along this gradient form a functional hierarchy whereby mostly rostral, higher-order regions differentially influence the activity of mostly caudal, lower-order regions. Though there is some debate as to which region is the ultimate “top” of this functional hierarchy (Badre & Nee, 2018; Nee & D'Esposito, 2016; Parkin, Hellyer, Leech, & Hampshire, 2015; Goulas, Uylings, & Stiers, 2014; Koechlin et al., 2003), there is strong evidence that interregional interactions increase along the rostrocaudal axis of LFC as rules become more complex, and they do so based on a priority relationship among regions depending on the task. For example, in a recent electrocorticography study, Voytek et al. (2015) showed that, during rule-guided behavior, rostral regions of LFC asymmetrically influenced gamma amplitude changes in caudal LFC via directional theta phase-to-gamma amplitude coupling. Importantly, this hierarchical, interregional theta phase-to-gamma amplitude coupling only occurred when the neurosurgery patients performed higher-order, abstract rules. Thus, there is evidence that regions along the rostrocaudal axis are not only functionally distinguishable but also, under certain task conditions, can interact with each other hierarchically to support higher-order rule use.
What might support this functional hierarchy? One class of hypotheses focuses on corticocortical connectivity (reviewed in Badre & D'Esposito, 2009) with putative asymmetry in connectivity among caudal-to-rostral frontal regions supporting increasingly higher-order control representations. Though early theory focused on corticocortical connections, recent hypotheses propose that these hierarchical relationships are established functionally via dynamic interaction with other mediating regions, such as the striatum. One such hypothesis proposes that the striatum supports hierarchical relationships among separate LFC regions via a cortico-basal ganglia-thalamic gating function (Chatham & Badre, 2015; Frank & Badre, 2012). Selection of a task-relevant LFC control representation may be accomplished by a striatal gating mechanism in which the LFC transmits candidate rule information to a local region of the striatum, which in turn selectively amplifies cortical processing of adaptive rules in the LFC region through its disinhibitory pallido-thalamo-cortical loop (Hazy, Frank, & O'Reilly, 2007; O'Reilly & Frank, 2006; Frank, Loughry, & O'Reilly, 2001; Braver & Cohen, 2000). Frank and Badre proposed that this striatal gating mechanism could be elaborated to support hierarchical rule use. Specifically, rostral, higher-order LFC regions could send top–down contextual information to the striatum that, through cortico-basal ganglia-thalamic loops, gate cortical processing in caudal, lower-order LFC regions, thereby allowing contextual representations maintained in rostral LFC to influence the selection of contexts or actions by caudal LFC (Badre & Frank, 2012; Frank & Badre, 2012). These properties may be reflected in the architecture of corticostriatal connections in two predictions. First, there might be an ordered relationship of corticostriatal connectivity reflecting candidate rule maintenance loops between the cortex and subregions of the striatum that maintain level-specific rule information. Second, there might be an asymmetry in the corticostriatal connectivity (Badre & Frank, 2012; Frank & Badre, 2012) such that higher-order cortical (e.g., more rostral LFC) regions connect with a broader caudal extent than caudal LFC regions do rostrally in the striatum, reflecting the influence of higher-order LFC regions on lower-order cortico-basal ganglia-thalamic loops. Thus, there may be an ordered but increasingly symmetric corticostriatal connectivity from higher- to lower-order cortical regions. Some initial evidence for these asymmetric connectional relationships between LFC and striatum has been observed using high angular resolution diffusion tractography (Verstynen, Badre, Jarbo, & Schneider, 2012).
Although much work has focused on the functional gradient within the LFC and its relationship to the striatum, these structures also interact with other cortical regions during cognitive control. Monkey anatomical tract-tracing studies show that the LFC has widespread connections to regions in frontal, parietal, and temporal association cortices (Saleem, Miller, & Price, 2014; Morecraft et al., 2012; Petrides & Pandya, 1999, 2006; Morecraft, Geula, & Mesulam, 1992; Cavada & Goldman-Rakic, 1989a, 1989b; Vogt & Pandya, 1987). Based on anatomical connectivity patterns, Goldman-Rakic (1988) suggested that there are multiple parallel, distributed association networks, each consisting of a complementary set of regions in frontal, parietal, and temporal association cortices and dedicated to a distinct function. Consistent with this model, recent network parcellations of the human cerebral cortex using resting-state functional connectivity MRI (fcMRI) also yield parallel, distributed association networks, recognizably associated with distinct functions (Lee et al., 2012; Power et al., 2011; Yeo et al., 2011) and functionally coupled to distinct regions of the striatum (Choi, Yeo, & Buckner, 2012). This parallel, distributed organization of association networks suggests that cognitive control arises from the interactions between the distributed regions of a network, rather than a single region acting alone (Cole, Bassett, Power, Braver, & Petersen, 2014; Mesulam, 1981, 1990; Goldman-Rakic, 1988). Indeed, numerous fMRI studies have reported activity in parietal and temporal cortices in addition to frontal cortex during cognitive control tasks, including those involving complex hierarchical cognitive control (Farooqui, Mitchell, Thompson, & Duncan, 2012; Kouneiher, Charron, & Koechlin, 2009; Badre & D'Esposito, 2007; Jubault, Ody, & Koechlin, 2007; Koechlin & Jubault, 2006; Koechlin et al., 2003). Importantly, however, frontal regions in these network parcellations from fcMRI are not arrayed in an obvious rostrocaudal manner (Power et al., 2011), nor do they show a clear hierarchical ordering in their resting correlation structure (Yeo et al., 2011). Thus, this raises a fundamental question: How do the rostral-to-caudal functional distinctions observed in LFC relate to the widely observed parallel, distributed organization of association networks?
In this study, we provide some initial answers to this question by contrasting task-based fMRI activations found during a hierarchical control study with parcellations from fcMRI. First, we assessed the overlap of fMRI task activity from a hierarchical cognitive control study (from Badre & D'Esposito, 2007) with a comprehensive estimate of human cortical networks in 1000 healthy, young adult participants using fcMRI (from Yeo et al., 2011) to directly determine their relationship. This comparison revealed that the separate regional activations identified with fMRI along the rostrocaudal axis of LFC were systematically associated with distinct functional networks, both in LFC and in parietal cortex. Next, given the striatal gating hypothesis and its importance for establishing hierarchical functional relationships, we investigated the functional connectivity between the LFC and parietal regions of these association networks and the striatum for agreement with the first of the predictions of the corticostriatal model. We note that we were unable to test the second prediction of asymmetric skew in the corticostriatal connectivity here with fcMRI because of confounds arising from the physical constraints of the striatum (see Discussion).
Hierarchical Rule Study
We reexamined fMRI task activity from a cognitively demanding hierarchical rule use study conducted by Badre and D'Esposito (2007). The relevant methods are described briefly here; see their work for further details.
The study involved 19 paid participants (10 men) who were 18–31 years, clinically normal, native English speakers, and right-handed. Written informed consent was obtained in accordance with procedures set by the Committee for Protection of Human Subjects at the University of California, Berkeley.
Task and Procedure
Conceptually, hierarchical rules consist of higher-order rules that specify the use of differing lower-order rules in different contexts. In Badre and D'Esposito's (2007) study, participants completed four independent hierarchical rule tasks that ranged from an easy first-order task to a highly demanding fourth-order task. Tasks had one, two, three, or four hierarchical levels of rule complexity (see Badre & D'Esposito, 2007, their Figure 2; also termed policy abstraction, see Badre, Kayser, & D'Esposito, 2010; Badre, 2008; Botvinick, 2008). In a single trial of each task, participants were shown a visual contextual cue and made a button press response according to the relevant rule. Trials had one, two, or four available rules at the highest order of each task (termed rule competition, see Badre & D'Esposito, 2007, their Figure 1), which served as a parametric variable approximating cognitive demand in the data analysis to identify regions specifically involved in cognitive control at each rule order.
MRI Data Acquisition
MRI data were collected on a 4-T Varian/Inova (Palo Alto, CA) MRI system. Structural data included anatomical images acquired by a T1-weighted gradient-echo multislice sequence and a high-resolution T1-weighted MP-FLASH 3-D sequence. The functional imaging data were acquired using a two-shot gradient-echo EPI sequence sensitive to BOLD contrast. Functional imaging parameters were as follows: repetition time (TR) = 2000 msec, echo time (TE) = 28 msec, 3.5 × 3.5 × 5 mm, 18 axial slices (whole-brain coverage), and 0.5 mm interslice gap voxels.
fMRI Data Preprocessing
fMRI data were corrected for slice acquisition-dependent time shifts and interpolated to a 1-sec resolution. Volumes that were large outliers in global signal were removed by replacing them with the global mean signal. Subsequent preprocessing using SPM2 (Wellcome Department of Cognitive Neurology, London, United Kingdom) included motion correction, normalization to Montreal Neurological Institute (MNI) space using a 12-parameter affine transformation and a nonlinear transformation, and resampling into 3-mm3 voxels. Data were then spatially smoothed with an 8-mm FWHM Gaussian kernel.
Data for each task were independently analyzed in SPM2 at the participant level using a general linear model, which included a regressor for rule competition. The resulting beta maps reflected the degree to which signal at each voxel parametrically correlated with rule competition at each order.
Cortical Network Parcellation and Resting-state Functional Connectivity Preprocessing
We used the 17-network cortical parcellation created by Yeo et al. (2011) and fully preprocessed fcMRI data from Choi et al. (2012) for the corticostriatal fcMRI analyses. Their relevant methods for creating the parcellation and preprocessing the fcMRI data are described briefly here; see their works for further details.
The study involved 1000 paid participants (427 men) who were aged 18–35 years (mean age = 21.3 years), clinically normal, native English speakers, and majority right-handed (91%). A subset (Discovery sample set, see Choi et al., 2012; Yeo et al., 2011) was used for the fcMRI analysis: 500 paid participants (213 male) were ages 18–35 years (mean age = 21.3 years), clinically normal, native English speakers, majority right-handed (91%), and with normal or corrected-to-normal vision. Participants were excluded from the study if their slice-based fMRI signal-to-noise ratio was low (<100; Van Dijk, Sabuncu, & Buckner, 2012) or artifacts were detected in the MRI data. Participants were prescreened for prior neurological or psychiatric illness, psychoactive medications, and MRI contraindications. Written informed consent was obtained in accordance with procedures set by institutional review boards of Harvard University or Partners Healthcare.
MRI Data Acquisition
MRI data were collected on matched 3-T Tim Trio scanners (Siemens, Erlangen, Germany) using the vendor-supplied 12-channel phased-array head coil. Structural data included a high-resolution multiecho T1-weighted magnetization-prepared gradient-echo image (multiecho MPRAGE). Structural scan (multiecho MPRAGE) parameters were as follows: TR = 2200 msec, inversion time = 1100 msec, TE = 1.54 msec for Image 1 to 7.01 msec for Image 4, flip angle = 7°, 1.2 × 1.2 × 1.2 mm, and field of view = 230. The functional imaging data were acquired using a gradient-echo EPI sequence sensitive to BOLD contrast. Functional imaging parameters were as follows: TR = 3000 msec, TE = 30 msec, flip angle = 85°, 3 × 3 × 3 mm voxels, field of view = 216, and 47 slices aligned to the AC–PC plane (whole-brain coverage). During the functional scans, participants were instructed to stay still, stay awake, and keep their eyes open.
fMRI Data Preprocessing
First four volumes of each run were discarded, slice acquisition-dependent time shifts were compensated per volume using SPM2, and head motion was corrected using rigid body translation and rotation using FMRIB Software Library (FSL). The data underwent further preprocessing specific to fcMRI analysis, including low-pass temporal filtering and regression for head motion, whole-brain signal, and ventricular and white matter signal.
Structural MRI Data Preprocessing and Surface Functional–Structural Data Alignment
For each participant, the cerebral cortex was extracted from the structural MRI image, modeled as a 2-D flat surface, and registered to a common spherical coordinate system using the FreeSurfer version 4.5.0 (surfer.nmr.mgh.harvard.edu). Functional images were aligned to the structural image using FsFast (surfer.nmr.mgh.harvard.edu/fswiki/FsFast). A 6-mm FWHM Gaussian smoothing kernel was applied to the fMRI data in the surface space, and data were downsampled to a 4-mm mesh.
Hybrid Surface- and Volume-based Alignment, Mapping between Surface- and Volume-based Coordinates, and Visualization
To perform the fcMRI analyses between the FreeSurfer cortical surface and volume, the surface and volume were aligned using a nonlinear volumetric registration algorithm (see Choi et al., 2012). Briefly, the striatum was modeled as a volume and aligned to the cortical surface. The participant's fMRI data were transformed into the common FreeSurfer nonlinear volumetric space and smoothed with a 6-mm FWHM smoothing kernel. Spatial correspondence between the FreeSurfer surface and volumetric coordinate systems was established by averaging over 1000 participants the transformation from each participant's native space to the FreeSurfer surface space and the transformation from the FreeSurfer nonlinear volumetric space. For figure images, the normalized FreeSurfer nonlinear volumetric data were then transformed into FSL MNI space. The striatum was defined using a FreeSurfer template mask of the striatum (Note that the mask does not include the tail of the caudate because of resolution limitations). All analyses were performed in FreeSurfer surface and volumetric spaces (fsaverage5 2 mm isotropic resolution). Results were displayed in MNI atlas space for the volume and the left and right inflated PALS cortical surfaces using Caret software (Van Essen, 2005) for the surface.
Regression of Cerebral Cortex Signal from Striatum
As described by Choi et al. (2012), the mean signal of cortical voxels adjacent to the left (within 4 voxels or 8 mm) or right (within 4.5 voxels or 9 mm) putamen was regressed from the entire left or right striatum at the participant level. This signal was determined to likely be spurious bleeding because (1) it was part of a continuous region of correlations primarily in the insula, but also extending into the claustrum and putamen; (2) the correlations in the striatum did not agree with the pattern of monkey anatomical projections from the insula to the striatum; and (3) regression of this signal revealed topographic motor correlations with the putamen as shown by monkey anatomy (see Choi et al., 2012, their Figure 13).
We noted another likely region of spurious signal covering the rostral ventral edge of the caudate and adjacent ventricles and spread into the septal wall in the fcMRI map for premotor dorsal (PMd). This signal was present even after regression of the mean ventricular signal during preprocessing. These correlations are contradictory with monkey anatomical projections from premotor cortex, which target the dorsal most part of the putamen and caudate, but not the rostral ventral caudate. Seeding the rostral ventral edge of the caudate, adjacent ventricular region, or septal wall yielded maps with broad, weak correlations in these regions, much of the entire medial wall of the cortex, and the lateral motor and premotor cortices. Because regression of this spurious signal from the data did not show qualitative changes to the resulting fcMRI maps, we decided not to perform this step. We note their presence with asterisks in Figures 4 and 5.
Creation of Network Parcellation
Yeo et al. (2011) divided the cerebral cortex into 18,715 uniformly distributed regions. At the individual participant level, for each region, a connectivity profile was computed of its Pearson's product moment correlation to 1175 uniform ROI vertices (∼16 mm apart). The resulting 18,715 × 1175 matrix was thresholded to the top 10% of correlations and binarized. Group matrix was obtained by averaging over 1000 participants. A k-means clustering method was used to empirically group regions into a defined number of networks based on their connectivity profiles. Clustering was repeatedly run with random splits of participants into two groups of n = 500 to assess the stability of different network number solutions. The 17-network parcellation was selected as a relatively stable solution with fine-grained detail.
Quantification of Task Activity Overlap on Network Parcellation
For each participant, SPM beta maps of activity parametrically modulating with rule competition were obtained for each of the four tasks. SPM beta maps were nonlinearly transformed from standard MNI space to FreeSurfer volumetric space and then projected to the FreeSurfer cortical surface for overlay with the cortical parcellation. Activity was thresholded at a very lenient threshold (0.05) to maximize sensitivity to overlap at the acceptable cost of biasing toward the null hypothesis (i.e., failing to detect differences in overlap). Activity was quantified by finding the mean beta value within the top three networks in the LFC that overlapped the most with the activity over the four orders (Figure 2) and contained regions that were consistent with prior knowledge of the location of the LFC functional gradient. These were a somatomotor-related network (green; Figure 1, top) and two association networks covering the middle (orange) or rostral lateral (maroon) LFC. Mean beta values were also obtained for parietal regions of the same networks. Statistical significance was determined by an Order × ROI (4 × 3) repeated-measures ANOVA. Mauchley's test of sphericity (i.e., the variances in the differences of all possible pairs of groups are the same) showed sphericity in the data for the parietal cortex but not the LFC. To address this, the Greenhouse–Geisser correction was applied to the LFC ANOVA to adjust the degrees of freedom and obtain a corrected F ratio.
To determine whether network or rostrocaudal position (estimated by y coordinate) significantly predicts the order of task activity, multiple and single linear regressions of these factors were performed on group mean task activity within the three networks in the left frontal or parietal cortex that were normalized by network and task and thresholded at 0.05. The statistical significance of the difference of model R2 values (R2network − R2y coordinate) was assessed using statistical bootstrapping to obtain its distribution. This difference distribution was created by resampling 1000 times from the set of single voxel values to obtain R2network, R2y coordinate, and their difference R2network − R2y coordinate. The 95% confidence interval was obtained from the difference distribution, and statistical significance was determined by examining if the confidence interval contained 0.
Resting-state Functional Connectivity Analysis
Analyses were conducted on fully preprocessed fcMRI data from the n = 500 subset (Discovery set) described above.
Cortical Seed Region Selection
In the LFC, fcMRI analyses were performed on two sets of cortical vertices. The first set consisted of the closest cortical surface vertices to the regions that were most parametrically correlated to rule competition in the first-, second-, third-, and fourth-order hierarchical rule tasks. The MNI coordinates for these cortical regions are as follows: −30, −10, 68 (first order; PMd); −38, 10, 34 (second order; pre-premotor dorsal [prePMd]); −50, 26, 24 (third order; inferior frontal sulcus [IFS]); and −36, 50, 6 (fourth order; frontal polar cortex [FPC]). The second set of cortical vertices was based on the first set but moved by hand to be approximately equidistant from nearby cortical network boundaries to avoid signal mixing with neighboring networks, which is a potential confound in the results. Seeding adjacent cortical vertices showed qualitatively similar fcMRI results. The MNI coordinates for these adjusted LFC regions are as follows: −23, −10, 56 (first order); −38, 4, 29 (second order); −54, 19, 22 (third order); and −37, 59, −2 (fourth order). In the parietal cortex, seed regions were selected to be located away from network boundaries and have representative fcMRI maps of their respective parietal parcellations. The MNI coordinates for parietal seed regions are as follows: −53, −26, 37 (somatomotor-related network); −52, −38, 51 (Association Network 1); and −52, −50, 47 (Association Network 2).
Resting-state Functional Connectivity Analysis
Striatal fcMRI maps for specific cerebral seed regions were obtained by computing the Pearson's product–moment correlation between the cortical region's fcMRI time course and the time courses of striatal voxels. Each cortical region consisted of a single surface vertex (∼4 × 4 mm) but should be considered spatially more extensive due to spatial smoothing. Group mean correlation z maps were obtained by converting the correlation (r) maps of individual participants to z maps using Fisher's r-to-z transformation and averaging across all participants in the group. For participants with multiple runs, the individual participant z maps were first averaged within participant and submitted to the group average. An inverse Fisher's r-to-z transformation was then applied to the group-averaged correlation z map, yielding a group mean correlation (r) map.
Corticostriatal Correlation Plots
To comprehensively evaluate the fcMRI correlations throughout the striatum, for each cortical region, the strongest correlation value at each coronal slice spanning the striatum in the group-level maps was plotted to create corticostriatal correlation curves. There were 21 coronal slices intersecting the striatum in FreeSurfer fsaverage5 2 mm space (spanning y = 24 to −20 in MNI space). Visual inspection showed that the correlation curves accurately reflected the overall pattern of correlations through the striatum. Each correlation curve was assessed for its mean by considering the corticostriatal correlation curves as a function r(y), where y is the y coordinate of the peak voxel in each of the 21 coronal slices through the striatum. We then calculated the means of the group mean correlation curves, defined as mean = . However, we could not assess the statistical significance using these mean values because some individual participant correlation curves contained negative correlation values, which lead to nonsensical statistics. Specifically, with σ = , a negative correlation value will lead to σ = , that is, an imaginary number, which cannot be used in statistical tests. Because we were only interested in assessing the relative statistical differences of these curves, not an inherent statistical property, we performed hypothesis tests on the nonnormalized means, defined as nonnormalized mean = first moment = Σir(yi) × yi, from the participants. Hypothesis tests consisted of Order × Mean (4 × 1) and Order × Skewness (4 × 1) repeated-measures ANOVAs for the LFC and Network × Mean (3 × 1) and Network × Skewness (3 × 1) repeated-measures ANOVAs for the parietal cortex. The relationship between Order and the x, y, and z coordinates of the peak corticostriatal correlation was estimated with single and multiple linear regressions.
Cognitive Control Recruits Gradients of Activity in the LFC and Parietal Cortex
We reexamined whole-brain fMRI maps from Badre and D'Esposito (2007) of activity parametrically covarying with rule competition (estimating cognitive control) specific to four orders of hierarchical rule use. The threshold of activity in these maps was lowered to compare in the most liberal way possible regional activation with the fcMRI cortical network parcellation. Given that our objective was to look for differences in overlap with separate functional networks, using a liberal threshold for the fMRI activation map biases in favor of the null hypothesis, as randomly spatially distributed spurious or noisy activation would reduce our sensitivity to locate differences.
Even at this highly liberal threshold, the rostral-to-caudal gradient of activity from the fourth to first orders was still distinct in the LFC, as originally reported (Badre & D'Esposito, 2007; Figure 1, top row). Activity was also present in distributed regions in the parietal cortex (Figure 1, middle row), medial frontal cortex (Figure 1, bottom row), and temporal cortex (Figure 1, top row). In particular, the activity in the parietal cortex was arrayed in a smaller and more overlapping but nonetheless distinct caudal-to-rostral functional gradient (Figure 1, middle row). Overall, the global, whole-brain pattern is one of localized activity in somatomotor cortex at lower orders of cognitive control that moves centrifugally in frontal and parietal cortices as higher orders of control are recruited. However, given the low statistical threshold used here, one cannot draw inferences about these patterns or distinguish any individual voxel activation from chance. Nevertheless, we report them, as they are the maps that are compared with the network parcellation in the subsequent analyses.
Cognitive Control Recruits Higher- and Lower-order Cortical Networks
We asked whether this broad patterning is random or if it has a systematic relationship with the independently defined cortical networks. To address this, the group mean fMRI map of activity parametrically covarying with rule competition at each rule order was overlaid with a previously reported 17-network fcMRI estimate of human cortical networks derived from 1000 healthy adult participants (Yeo et al., 2011; Figure 1, top). In the LFC, parametric activation related to the first-order rule task overlapped primarily with regions belonging to a somatomotor-related network (green), for the second and third orders with an association network (orange; Association Network 1), and the fourth order with a second, separate association network (maroon; Association Network 2; Figure 2, left column). Together, these three networks overlapped with 64% or more of the parametric activation for each order in the LFC.
Figure 3A shows a quantification of this activity across the four orders and three networks. Black bars corresponding to the first rule order show the greatest parametric activation in the somatomotor-related network and progressively less activation for Association Networks 1 and 2 (Figure 3A, top). The opposite pattern holds for the fourth-order parametric activation shown in light blue bars, which is least in the somatomotor-related network and progressively increases in Association Networks 1 and 2. The second and third rule orders show the greatest parametric activation in Association Network 1 and less activation in the somatomotor-related network and Association Network 2. A repeated-measures ANOVA (Greenhouse–Geisser corrected) showed a significant Order × ROI interaction effect, F(2, 42) = 8.99, p < .001.
To determine whether network or rostrocaudal position was a significant predictor of the order of task activity and, if so, which factor better explained the variance in order, multiple and single linear regressions of these factors were performed on the normalized group mean task activity in the three networks (Table 1). We sought to test whether the frontal BOLD activity observed at each order was better described as progressively following a rostrocaudal gradient, such as along a single dimension of abstraction, or, rather, followed the distributed network of frontal regions. These networks are rostrocaudally arrayed in one part of LFC, but there are exceptions that allow us to distinguish this pattern from a basic rostrocaudal gradient. Thus, we took the BOLD activity in the three networks of interest in the frontal cortex from the four task orders and tested the hypotheses that the order of the BOLD activity is better explained by a given voxel's network membership or, alternatively, purely by rostrocaudal location (approximated by the voxel's y coordinate). The multiple regression showed that both factors are significant predictors of order (βNetwork = 0.6867, p < .001; βy = 0.0201, p < .001) and explained 31.6% of the variance in order, F(2, 1805) = 417.9, p < .001. However, single linear regressions showed that network explained more of the variance in order (R2 = 0.301, F(1, 1806) = 776.9, p < .001) than y coordinate did (R2 = 0.213, F(1, 1806) = 489.2, p < .001). The 95% confidence interval for the distribution of the difference between these two model R2 values is [0.03, 0.14], which does not contain 0. These results indicate that network membership, rather than rostrocaudal position, is a stronger predictor of the order of a given voxel's BOLD activity in the LFC during hierarchical cognitive control.
|F Statistic/p Value||t Statistic/p Value||t Statistic/p Value||t Statistic/p Value|
|Left frontal cortex||2||1805||0.316a||0.6867||0.0201||−0.1769|
|Left frontal cortex||1||1806||0.301||0.8668||—||0.8865|
|Left frontal cortex||1||1806||0.213||—||0.0547||−1.2613|
|Left parietal cortex||2||1775||0.434a||0.5738||−0.0819||4.5647|
|Left parietal cortex||1||1776||0.361||0.9752||—||0.6143|
|Left parietal cortex||1||1776||0.366||—||−0.1350||7.6950|
|F Statistic/p Value||t Statistic/p Value||t Statistic/p Value||t Statistic/p Value|
|Left frontal cortex||2||1805||0.316a||0.6867||0.0201||−0.1769|
|Left frontal cortex||1||1806||0.301||0.8668||—||0.8865|
|Left frontal cortex||1||1806||0.213||—||0.0547||−1.2613|
|Left parietal cortex||2||1775||0.434a||0.5738||−0.0819||4.5647|
|Left parietal cortex||1||1776||0.361||0.9752||—||0.6143|
|Left parietal cortex||1||1776||0.366||—||−0.1350||7.6950|
Adjusted R2 for two factors.
We then tested how the caudal-to-rostral functional gradient in the parietal cortex overlapped with the parietal regions of the same three cortical networks overlapping with the LFC functional gradient (Figure 2B, bottom). In parietal cortex, the parametric activation related to each rule level overlapped with the corresponding network that overlapped that level in LFC. Specifically, the greatest activation in the parietal cortex was located in the somatomotor-related network (green) for the first-order task, Association Network 1 (orange) for the second- and third-order tasks, and Association Network 2 (maroon) for the fourth-order task (Figure 2, right column). Together, these three networks overlapped with 67% or more of the parametric activation for each order in parietal cortex.
Figure 3B shows a quantification of this parametric activation across the four orders and three networks. Again, first-order activity showed the greatest activation with the somatomotor-related network and progressively decreases in Association Networks 1 (orange) and 2 (maroon), the second- and third-order activities show the greatest activation with Association Network 1, and the fourth-order activity shows the least activation with the somatomotor-related network and progressively increases in Association Networks 1 and 2 (Figure 3B). A repeated-measures ANOVA showed a significant Order × ROI interaction effect, F(6, 108) = 9.74, p < .001.
To examine whether network or rostrocaudal position of task activity better explains order, a multiple linear regression was performed that showed that both network and rostrocaudal position (y coordinate) are significant predictors of order (βNetwork = 0.5738, p < .001; βy = −0.0819, p < .001) and explained 43.4% of the variance in order, F(2, 1775) = 682.3, p < .001 (Table 1). However, unlike the LFC, single linear regressions showed a similar amount of variance in order explained by network (R2 = 0.361, F(1, 1776) = 1003, p < .001) and y coordinate (R2 = 0.366, F(1, 1776) = 1027, p < .001). The 95% confidence interval for the distribution of the difference between these two model R2 values is [−0.06, 0.04], which does contain 0. These results indicate that network membership and rostrocaudal position are similarly strong predictors of the order of a given voxel's BOLD activity in the parietal cortex during hierarchical cognitive control. This is consistent with the single, rostrocaudally arrayed set of network regions in the parietal cortex, in contrast to the multiple, distributed nature of regions in the LFC.
Functional Connectivity between Higher- and Lower-order Cortical Networks and the Striatum
We next investigated whether the functional connectivity between those networks and the striatum is consistent with the corticostriatal gating model of hierarchical rule use (Badre & Frank, 2012; Frank & Badre, 2012). Based on this model, each LFC region should feature ordered, spatially distinct, and progressively asymmetric corticostriatal fcMRI connectivity for lower- versus higher-order regions. Specifically, we predicted that LFC regions would have correlations through the caudate of the striatum with the peak correlation at a distinct rostrocaudal site corresponding to the rostrocaudal location of the LFC region based on prior monkey and human corticostriatal connectivity results (Verstynen et al., 2012; Kemp & Powell, 1970). Importantly, we expected to see a similar pattern for parietal regions that are members of the same functional networks.
We first examined striatal functional connectivity with the four LFC regions most correlated with rule competition at each task order (PMd, prePMd, IFS, and FPC from Badre & D'Esposito, 2007) located in the somatomotor-related network and Association Networks 1 and 2. (Note that the exception was the third-order region [IFS], which was not located in Association Network 1, but rather in an adjacent network. However, because the majority of the covarying fMRI activity at the third order was located in Association Network 1, we nonetheless believe that this is the involved network.) However, because these seed regions were close to the borders of the association networks where there is a potential blurring of signals from neighboring networks, we created a second set of seed regions placed away from the network borders (PMd′, prePMd′, IFS′, and FPC′; Figure 4B, top). The first-order seed region (PMd′) was primarily correlated with the caudal putamen, consistent with its close functional link to the motor cortex, and weakly correlated with the caudal body and tail of the caudate (Figure 4A; note that correlations in the rostral ventral caudate are likely due to signal bleeding from adjacent ventricles and gray and white matter; see asterisks in Figure 4B–D and Methods). The second (prePMd′), third (IFS′), and fourth (FPC′) order regions showed longitudinal correlations throughout the rostrocaudal extent of the caudate and parts of the medial putamen (Figure 4A).
The strongest correlations for 21 evenly spaced coronal slices of the caudate were plotted to create a corticostriatal curve of correlation for each LFC region (Figure 4B). These curves revealed that, consistent with our first prediction, the fourth-order map had the most rostral peak correlation, followed by successively caudal peaks for the third- and second-order maps. A linear regression showed a dependence of peak correlation location in the y axis on order (βy = 3.26, p < .001). This was consistent with the rostral-to-caudal locations of the means of each curve, which were located at fourth order (z = 6.92), third order (z = 3.53), second order (z = −1.55), and first order (z = −14.99). A repeated-measures ANOVA showed a significant Order × Mean interaction effect, F(3, 1996) = 101.43, p < .001. Follow-up paired t tests were significant between all orders (p < .005).
In addition to the rostral–caudal (y) axis, the medial–lateral (x), and dorsal–ventral (z) axes were examined and found to have smaller shifts in mean peak correlation location from lower to higher order. From first to fourth order, the mean peak correlation shifted from medial to laterally for the x axis (Figure 4C) and from ventral to dorsally for the z axis (Figure 4D). A single linear regression showed a dependence of the peak correlation location in the x or z axis on order (βx = −0.50 and βz = 1.91, both ps < .001). A multiple linear regression of all three x, y, z axes on order indicated a statistically significant rostrolateral shift in peak correlation from lower to higher order (βx = −0.11, βy = 0.22, both ps < .001; βz = −0.06, p = .07).
We next examined the striatal functional connectivity for parietal seed regions in the somatomotor-related network and Association Networks 1 and 2 (Figure 5B, top). If parietal regions are recruited with LFC regions as part of cortical networks, we predicted a similar rostral-to-caudal differential of striatal connectivity from caudal-to-rostral (higher- to lower-order) parietal regions. Indeed, similar to the first-order LFC region (PMd′), the seed region in the somatomotor-related network was most strongly correlated with the putamen (Figure 5A) and weakly correlated with the tail of the caudate (Figure 5B–D; note that correlations in the rostral ventral caudate are likely due to signal bleeding from adjacent ventricles and gray and white matter; see asterisks in Figure 5B–D and Methods). Like the second-, third-, and fourth-order LFC regions (prePMd′, IFS′, FPC′), Association Networks 1 and 2 seed regions in the parietal cortex also showed longitudinal correlations in the caudate and medial putamen (Figure 5A).
A plot of peak correlations from 21 coronal slices through the caudate showed a rostrocaudal difference in the location of the peak correlation for each order (Figure 5B, bottom), consistent with the first prediction. A linear regression showed a dependence of peak correlation location in the y axis on order (βy = 4.58, p < .001). This was consistent with the means of the curves, which were located at Association Network 2 (z = 5.45), Association Network 1 (z = −2.18), and somatomotor-related network (z = −14.22). A repeated-measures ANOVA showed a significant Network × Mean interaction effect, F(2, 1497) = 132.97, p < .001. Follow-up paired t tests were significant between all networks (p < .001).
As with the LFC, the x (Figure 5C) and z (Figure 5D) axes showed smaller medial to lateral and dorsal to ventral shifts, respectively, in the positions of the peak correlations from lowest to highest order. A linear regression showed a dependence of the peak correlation location in the x or z axis on order (βx = −0.69, βz = 2.50, both ps < .001). A multiple linear regression of all three x, y, z axes on order did not yield any statistically significant regression coefficients (βx = −0.10, p = .54; βy = −0.14, p = .71; βz = 0.31, p = .57). Altogether, we found that parietal regions correlate with the striatum with similar connectivity as the corresponding LFC regions. These regions correlate with primarily rostrocaudal differences in peak correlation, consistent with the first model prediction.
Hierarchical cognitive control has been associated with a gradient of higher- to lower-order rostral-to-caudal regions in the LFC. At the same time, monkey anatomical tract-tracing and human resting-state fcMRI shows that the LFC is part of parallel and distributed association networks. How do these regional versus network views of LFC functional organization relate to one another? We report the following observations addressing this question. (1) Distributed activity at each order of cognitive control is not restricted to LFC but also includes regions of medial frontal, parietal, and temporal cortices. In particular, in addition to the previously observed rostral-to-caudal gradient of activity in the LFC, an approximate caudal-to-rostral gradient is present in the parietal cortex from higher- to lower-order cognitive control. (2) Activity related to each order of cognitive control preferentially overlaps with a particular network. The same network preferences are observed in both frontal and parietal cortex. Crucially, for at least LFC, network membership better accounted for hierarchical order than did absolute rostrocaudal position. (3) The rostral-to-caudal LFC and caudal-to-rostral parietal regions of the same cortical network both have similar patterns of striatal connectivity such that higher-order networks show peak correlation with more rostral parts of the striatum. Thus, to answer the question above, our results suggest that there are specific higher- and lower-order association networks that are differentially recruited during hierarchical cognitive control tasks. The distributed regions of different networks form apparent local functional gradients located throughout association cortex. The largest, most robust one of these is the LFC gradient, and at least one other functional gradient located in the parietal cortex is also involved. These results suggest that it is network membership rather than rostrocaudal locus that determines rank in a hierarchy. Furthermore, the respective functional influence of these networks on each other may partly arise from the differing influence each network has over the striatum. Thus, these data are consistent with a model whereby these association networks are connectionally organized mostly in parallel and broadly support cognitive control but can be functionally organized hierarchically via dynamic interaction with the striatum.
Several lines of study have suggested a hierarchical relationship among distributed association networks. Throughout the cortex, there are broad gradients of cytoarchitectonic differentiation (Pandya & Yeterian, 1990) and anatomical connectivity (Mesulam, 1998; Pandya & Yeterian, 1990; Jones & Powell, 1970) from sensory to sensory association to high-order association cortex, including in the LFC and parietal cortex. Regions of higher or lower hierarchical level within these distributed gradients are connected by long-range anatomical connections. For example, there are connections between motor-related regions in the LFC and parietal cortex and centrifugally farther pairs of LFC and parietal regions in association cortex (Margulies & Petrides, 2013; Pandya & Yeterian, 1990; Chavis & Pandya, 1976). Based on such observations, Fuster's model proposes a hierarchical organization of cognitive networks involving the LFC and parietal cortex (Fuster, 2009). Here, we identify three networks in particular that overlap primarily with hierarchical functional activation. Across prior functional studies of hierarchical cognitive control, LFC regions similar to the PMd, IFS, and FPC are readily observed to dissociate and approximately map onto this tripartite network distinction (Nee & D'Esposito, 2016; Bahlmann et al., 2015; Nee & Brown, 2013; Barbalat et al., 2011; Koechlin et al., 2003). We note, however, that there are cases in which regions belong to the same network but do not always coactivate during tasks. Notably, prePMd and the more rostral IFS both belong to Association Network 1 but have been observed in some studies to show separate response patterns (Nee & D'Esposito, 2016; Chatham, Frank, & Badre, 2014; Badre & D'Esposito, 2007; Koechlin et al., 2003). This suggests that the network parcellations delineate regions that are overall more functionally associated than others but do not necessarily coactivate within network nor dissociate across networks for all tasks.
A second key point is that, though there was clear preference for a given network for each level of control, the overlap also was not complete, instead reflecting a differential mixing of the networks. For example, as can be seen in Figure 2, activity related to the highest level of control is most associated with Association Network 2. However, it also has some overlap with Association Network 1. Conversely, the third-order control activity overlaps primarily with Association Network 1 but includes some overlap with Association Network 2, as well. What accounts for this mixing? It is possible that this gradient of network involvement reflects noise. Another possibility is that the network boundaries themselves are estimated from fcMRI data in separate participants. As both the task and individuals differ, slight differences in these network boundaries might blur sharper functional bounds. Finally, each network may not correspond to a particular single level of control per se, but solving each control problem requires interaction among the networks to differing degrees.
Consistent with the hypothesis that hierarchical relationships are established, at least in part, through interaction with the striatum, we found for both LFC and parietal components of each network that the overall pattern of striatal connectivity was primarily a higher- to lower-order rostrocaudal differential in peak correlation locus. This pattern is consistent with the first prediction of the corticostriatal gating model that these separate networks are regulated by distinct, local cortico-basal ganglia-thalamic loops. These rostral-to-caudal LFC and caudal-to-rostral striatal correlations agree with diffusion MRI studies showing a rostral-to-caudal correspondence in fibers from the LFC to the striatum (Verstynen et al., 2012; Draganski et al., 2008). This is also consistent with anatomical tract-tracing observations that cortical regions project longitudinally through the striatum (Selemon & Goldman-Rakic, 1985; Goldman & Nauta, 1977; Künzle, 1975) with, it is thought, the densest projections from a cortical region terminating at its most proximal striatal region (Kemp & Powell, 1970). Interestingly, however, the reverse organization is true for parietal cortex: caudal-to-rostral parietal regions have a rostral-to-caudal differential in striatal connectivity. A direct comparison of specific prefrontal and parietal monkey anatomical tract-tracing cases from the literature shows evidence for the organization seen here. Figure 6 shows striatal projection patterns for three concentric pairs of anatomically connected LFC and parietal regions, each of which has a unique pattern of striatal connectivity that is shared by the LFC and parietal region of each pair. Furthermore, these projection patterns resemble the lower- to higher-order human fcMRI correlation patterns for the LFC and parietal cortex (Figures 4A and 5A). Although it is difficult to discern where the densest projection sites are without directly examining the anatomical tissue, these cases show that prefrontal and parietal regions share similar patterns of anatomical projections to the striatum and that these projection patterns bear close resemblance to those seen with fcMRI here. In addition, these fcMRI patterns are consistent with overlapping prefrontal and parietal connection zones observed in the striatum from monkey and human connectivity studies (Jarbo & Verstynen, 2015; Selemon & Goldman-Rakic, 1985), including evidence that these cortical regions form interconnected networks with distinct patterns of associated striatal connectivity (Choi, Tanimura, Vage, Yates, & Haber, 2017; Choi et al., 2012; Yeterian & Van Hoesen, 1978; but see Selemon & Goldman-Rakic, 1985).
A second prediction of the corticostriatal gating model is that there may be an asymmetric skew of the corticostriatal correlation curves that is reflective of the influence of higher-order regions on caudal regions. In the present study, we were unable to test this prediction with fcMRI due to confounds arising from the physical constraints of the striatum, which create a ceiling-floor effect that artificially affects assessments of asymmetric connectivity across the striatum. Any curves that peak at the ends of the striatum (e.g., fourth-order curve) artificially appear to have greater asymmetry simply from being “cut off” before it resolves to 0. Thus, for technical reasons, we are unable to test for an asymmetric, hierarchical relationship between the curves with fcMRI. We note that evidence for an asymmetric structure of LFC to striatal connections was found previously with diffusion tractography (Verstynen et al., 2012). Future study using neurophysiological methods that can observe dynamic, real-time activity would be a next step to answering this question.
What is the role of the parietal cortex within these hierarchical networks? By virtue of the parietal cortex's high-order, multimodal role in processing and integrating complex visual and somatomotor properties, the parietal cortex likely performs this processing in an order-specific manner during hierarchical cognitive control. Collins and Frank (2013) have suggested and provided evidence for a more involved role of the parietal cortex that represents the visual properties of the stimulus, combines this with context-related information from the pFC, and sends this multiplexed context–stimulus information to the premotor corticostriatal loop for action selection. Within the multiple, hierarchical networks schema reported here, a similar process may occur at each order, where a parietal region processes the visual properties of the relevant, order-specific features of the context and stimulus and creates a multiplex of context–stimulus information that is sent to the pFC, parietal, and striatal regions of the downstream, lower-order corticostriatal network.
Altogether, this study broadens a regional view of the functional organization of the LFC by demonstrating the relationship between the functional gradient observed in the LFC during hierarchical cognitive control and an fcMRI-based estimate of distributed cortical networks. We observe evidence that higher- and lower-order association networks are differentially recruited across levels of hierarchical cognitive control and that the distributed regions of these networks form local functional gradients that have been observed throughout association cortex. Furthermore, the LFC and parietal regions of these networks have similar correlation patterns within the striatum that are consistent with cases reported in the anatomical literature. The present identification of specific hierarchically organized networks serves as a basis for further functional, anatomical, and neurocomputational investigations of between-network organization and interactions.
The authors thank Drs. Randy Buckner, Christopher Chatham, and Derek Nee for valuable discussion and Drs. Thomas Yeo, Avram Holmes, Fenna Krienen, and Julia Lehman for technical assistance. This work was supported by a National Science Foundation Graduate Research Fellowship (E. Y. C.). The authors declare no competing financial interests.
Reprint requests should be sent to Eun Young Choi, Department of Neurosurgery, Stanford University, 318 Campus Drive, W1 103, Stanford, CA 94305, or via e-mail: email@example.com.