The Brain’s Topographical Organization Shapes Dynamic Interaction Patterns That Support Flexible Behavior Based on Rules and Long-Term Knowledge

Spatial working memory task

Participants were required to maintain in memory four or eight sequentially presented locations in a 3 × 4 grid (Fedorenko et al., 2011), giving rise to easy and hard spatial working memory conditions. Stimuli were presented at the center of the screen across four steps. Each of these steps lasted for 1 s and highlighted one location on the grid in the easy condition and two locations in the hard condition. This was followed by a decision phase, which showed two grids side by side [i.e., two-alternative forced choice (2AFC) paradigm]. One grid contained the locations shown on the previous four steps, while the other contained one or two locations in the wrong place. Participants indicated their response via a button press and feedback was immediately provided within 2.5 s. Each run consisted of 12 experimental blocks (6 blocks per condition and 4 trials in a 32 s block) and 4 fixation blocks (each 16 s long), resulting in a total time of 448 s.

Math task

Participants were presented with an addition expression on the screen for 1.45 s and subsequently made a 2AFC decision, indicating their solution within 1 s. The easy condition used single-digit numbers, while the hard condition used two-digit numbers. Each trial ended with a blank screen lasting for 0.1 s. Each run consisted of 12 experimental blocks (with 4 trials per block) and 4 fixation blocks, resulting in a total time of 316 s.

Semantic feature matching task

Participants were required to make a yes/no decision, matching probe and target concepts (presented as words) according to a particular semantic feature (color or shape), specified at the top of the screen during each trial. The feature prompt, probe word, and target words were presented simultaneously. Half of the trials were matching trials in which participants would be expected to identify shared features, while the other half were non-matching trials in which participants would not be expected to identify shared features. For example, in a color matching trial, participants would answer “yes” to the word pair DALMATIANS–COWS, due to their color similarity, whereas they would answer “no” to COAL–TOOTH as they do not share a similar color.

We parametrically manipulated the degree of feature similarity between the probe and target concepts, using semantic feature similarity ratings taken from a separate group of 30 participants on a 5-point Likert scale. For instance, in color matching trials, the degree of color similarity between DALMATIANS and COWS was found to be very high (i.e., 4.8), while that between PUMA and LION was relatively low (i.e., 4.0), despite participants believing that the two trials had similar color. Conversely, in color non-matching trials, the degree of color similarity between CROW and HUMMINGBIRD was relatively high (i.e., 2.5), whereas that between COAL and TOOTH was very low (i.e., 1.2), even though participants perceived no similarity in color. For the matching trials, greater feature similarity facilitates the decision-making process, while for the non-matching trials, greater feature similarity makes the decision more difficult. This parametric design allowed us to model the effect of the difficulty of semantic decision-making in the neural data, and to test whether control subnetworks showed similar or opposite activation patterns.

This task included four runs and two conditions (two features, color and shape), presented in a mixed design. Each run consisted of four experimental blocks (two 2 min30 s blocks per feature), resulting in a total time of 10 min12 s. In each block, 20 trials were presented in a rapid event-related design. To maximize the statistical power of the rapid event-related fMRI data analysis, the stimuli were presented with temporal jitter randomized from trial to trial (Dale, 1999). The intertrial interval varied from 3 to 5 s. Each trial started with a fixation, followed by the feature, probe word, and target word presented centrally on the screen, triggering the onset of the decision-making period. The feature, probe word, and target word remained visible until the participant responded or for a maximum of 3 s. The condition order was counterbalanced across runs, and run order was counterbalanced across participants. Half of the participants pressed a button with their right index finger to indicate a matching trial and responded with their right middle finger to indicate a non-matching trial. The other half of the participants pressed the opposite buttons.

Semantic association task

Participants were asked to decide if pairs of words were semantically associated or not (i.e., yes/no decision as above) based on their own experience. Overall, there were roughly equal numbers of “related” and “unrelated” responses across participants. The same stimuli were used in both the semantic feature matching task and semantic association task. For example, DALMATIANS and COWS are semantically related; COAL and TOOTH are not. The feature and association tasks were separated by 1 week. Similarly, we parametrically manipulated the semantic association strength between the probe and target concepts, using semantic association strength ratings taken from a separate group of 30 participants on a 5-point Likert scale. For example, in related trials, the association strength between PUMA and LION is very strong, while for COWS and WHALE, it is relatively weak (although they are still both animals). In unrelated trials, the association strength between KINGFISHER and SCORPION is relatively high, while that between BANANA and BRICK is very low, although participants thought neither were related. For the related trials, stronger associations would facilitate decision-making, while for unrelated trials, stronger associations interfere with the decision-making. This parametric design allowed us to model the effect of the difficulty of semantic decision-making in the neural data and to test whether control subnetworks showed similar or opposite activation patterns.

This task included four runs, presented in a rapid event-related design. Each run consisted of 80 trials, with approximately half being related and half being unrelated trials. The procedure was the same as the feature matching task, except only two words were presented on the screen.

Image acquisition of HCP dataset

MRI acquisition protocols of the HCP dataset have been described (Barch et al., 2013; Glasser et al., 2013). Images were acquired using a customized 3 T Siemens Connectome scanner having a 100 mT/m SC72 gradient set and using a standard Siemens 32-channel radiofrequency receive head coil. Participants underwent the following scans: structural [at least one T1-weighted (T1w) magnetization-prepared rapid gradient echo (MPRAGE) and one 3D T2-weighted (T2w) sampling perfection with application optimized contrast (SPACE) scan at 0.7 mm isotropic resolution] and rsfMRI (4 runs × 14 min and 33 s). Since not all participants completed all scans, we only included 339 unrelated participants from the S900 release. Whole-brain rsfMRI and task fMRI data were acquired using identical multiband echoplanar imaging (EPI) sequence parameters of 2 mm isotropic resolution with a repetition time (TR) of 720 ms.

Subjects were considered for data exclusion based on the mean and mean absolute deviation of the relative RMS motion across either four rsfMRI scans or one diffusion MRI scan, resulting in four summary motion measures. If a subject exceeded 1.5 times the interquartile range (in the adverse direction) of the measurement distribution in two or more of these measures, the subject was excluded. In addition, functional runs were flagged for exclusion if >25% of frames exceeded 0.2 mm FD (FD_power). These above exclusion criteria were established before performing the analysis (Faskowitz et al., 2020; Sporns et al., 2021). The data of 91 participants were excluded because of excessive head motion and the data of another 3 participants were excluded because their resting data did not have all the time points. In total, the data of 245 participants were analyzed after exclusions.

Image acquisition of York nonsemantic dataset

MRI acquisition protocols have been described previously (Wang et al., 2020, 2021). Structural and functional data were collected on a Siemens Prisma 3 T MRI scanner at the York Neuroimaging Centre. The scanning protocols included a T1w MPRAGE sequence with whole-brain coverage. The structural scan used the following: acquisition matrix of 176 × 256 × 256 and voxel size 1 × 1 × 1 mm³; TR, 2,300 ms; and echo time (TE) = 2.26 ms. Functional data were acquired using an EPI sequence with an 80º flip angle and using generalized autocalibrating partially parallel acquisitions (GRAPPA) with an acceleration factor of 2 in 3 × 3 × 4 mm voxels in 64-axial slices. The functional scan used the following: 55 3-mm-thick slices acquired in an interleaved order (with 33% distance factor); TR, 3,000 ms; TE, 15 ms; and field of view (FOV), 192 mm.

Image acquisition of York semantic dataset

Whole-brain structural and functional MRI data were acquired using a 3 T Siemens MRI scanner utilizing a 64-channel head coil, tuned to 123 MHz at York Neuroimaging Centre, University of York. The functional runs were acquired using a multiband multiecho (MBME) EPI sequence, each 11.45 min long [TR, 1.5 s; TE, 12, 24.83, 37.66 ms; 48 interleaved slices per volume with slice thickness of 3 mm (no slice gap); FOV, 24 cm (resolution matrix = 3 × 3 × 3; 80 × 80); 75° flip angle; 455 volumes per run; 7/8 partial Fourier encoding and GRAPPA (acceleration factor, 3, 36 ref. lines); multiband acceleration factor, 2]. Structural T1w images were acquired using an MPRAGE sequence (TR, 2.3 s; TE, 2.3 s; voxel size, 1 × 1 × 1 isotropic; 176 slices; flip angle, 8°; FOV, 256 mm; interleaved slice ordering). We also collected a high-resolution T2w scan using an EPI sequence (TR, 3.2 s; TE, 56 ms; flip angle, 120°; 176 slices; voxel size, 1 × 1 × 1 isotropic; FOV, 256 mm).

Image preprocessing of HCP dataset

We used HCP's minimal preprocessing pipelines (Glasser et al., 2013). Briefly, for each subject, structural images (T1w and T2w) were corrected for spatial distortions. FreeSurfer v5.3 was used for accurate extraction of cortical surfaces and segmentation of subcortical structures (Dale et al., 1999; Fischl et al., 1999). To align subcortical structures across subjects, structural images were registered using nonlinear volume registration to the Montreal Neurological Institute (MNI152) space. Functional images (rest and task) were corrected for spatial distortions, head motion, and mapped from volume to surface space using ribbon-constrained volume to surface mapping.

Subcortical data were also projected to the set of extracted subcortical structure voxels and combined with the surface data to form the standard Connectivity Informatics Technology Initiative (CIFTI) grayordinate space. Data were smoothed by a 2 mm full width at half maximum (FWHM) kernel in the grayordinates space. Rest data were additionally identically cleaned for spatially specific noise using spatial independent component analysis (ICA) + FMRIB’s ICA-based X-noiseiifier (FIX) (Salimi-Khorshidi et al., 2014) and global structured noise using temporal ICA (Glasser et al., 2018). For accurate cross-subject registration of cortical surfaces, a multimodal surface matching (MSM) algorithm (Robinson et al., 2014) was used to optimize the alignment of cortical areas based on features from different modalities. MSMSulc (“sulc”: cortical folds average convexity) was used to initialize MSMAll, which then utilized myelin, resting-state network, and rfMRI visuotopic maps. Myelin maps were computed using the ratio of T1w/T2w images (Salimi-Khorshidi et al., 2014). The HCP's minimally preprocessed data include cortical thickness maps (generated based on the standardized FreeSurfer pipeline with combined T1/T2 reconstruction). For this study, the standard-resolution cortical thickness maps (32,000 mesh) were used.

Image preprocessing of York non-semantic and semantic dataset

The York datasets were preprocessed using fMRIPrep 20.2.1 (Esteban et al., 2018; RRID:SCR_016216), which is based on Nipype 1.5.1 (Gorgolewski et al., 2011; RRID:SCR_002502). The details about preprocessing using fMRIPrep are provided in the subsequent subsections.

Individual-specific parcellation

Considering the anatomical and functional variability across individuals (Mueller et al., 2013; Laumann et al., 2015; Braga and Buckner, 2017; Gordon et al., 2017), we estimated individual-specific areal-level parcellation using a multisession hierarchical Bayesian model (MS-HBM; Kong et al., 2019, 2021). To estimate individual-specific parcellation, we acquired “pseudo-resting state” time series, in which the task activation model was regressed from spatial working memory, math, feature matching, and semantic association fMRI data (Fair et al., 2007) using xcp_d (https://github.com/PennLINC/xcp_d). The task activation model and nuisance matrix were regressed out using AFNI's 3dTproject (for similar implementation, see Cui et al., 2020).

Using a group atlas, this method calculates intersubject resting-state functional connectivity variability and intrasubject resting-state functional connectivity variability and finally parcellates each single subject based on this prior information. As in Kong et al. (2019, 2021), we used MS-HBM to define 400 individualized parcels belonging to 17 discrete individualized networks for each participant, where the FPCN was divided into 3 subnetworks. This allowed us to explore the heterogeneity of FPCN. Specifically, we calculated all participants' connectivity profiles, created the group parcellation using the average connectivity profile of all participants, estimated the intersubject and intrasubject connectivity variability, and finally calculated each participant's individualized parcellation. This parcellation imposed the Markov random field spatial prior. We used a well-known areal-level parcellation approach, that is, the local gradient approach (gMS-HBM), which detects local abrupt changes (i.e., gradients) in resting-state functional connectivity across the cortex (Cohen et al., 2008). A previous study (Schaefer et al., 2018) has suggested combining local gradient (Cohen et al., 2008; Gordon et al., 2016) and global clustering (Yeo et al., 2011) approaches for estimating areal-level parcellations. Therefore, we complemented the spatial contiguity prior in contiguous MS-HBM (cMS-HBM) with a prior based on local gradients in resting-state functional connectivity. This encouraged adjacent brain locations with gentle changes in functional connectivity to be grouped into the same parcel. We used the pair of parameters (i.e., β value = 50; w = 30; and c = 30), which was optimized using our own dataset. The same parameters were also used in Kong et al. (2021). Vertices were parcellated into 400 cortical regions (200 per hemisphere). To parcellate each of these parcels, we calculated the average time series of enclosed vertices to improve the signal-to-noise ratio using Connectome Workbench software. This parcel-based time series was used for all the subsequent analyses. The same method and parameters were used to generate the individual-specific parcellation for participants in the HCP dataset using the resting-state time series except that the task regression was not performed.

Cortical geometry—global minimum distance to primary sensory-motor landmarks

To reveal how physical proximity to structural landmarks corresponding to primary sensory-motor regions influences the function of regions, we calculated the geodesic distance between each parcel and key landmarks associated with the primary visual, auditory, and somatomotor cortex. These values were used to identify the minimum geodesic distance to primary sensory-motor regions for each parcel. Three topographical landmarks were used: the central sulcus corresponding to the primary somatosensory/motor cortex; temporal transverse sulcus indicating primary auditory cortex; and calcarine sulcus, marking the location of primary visual cortex. Since the cortical folding patterns vary across participants, and the individual variability in cortical folding increases with cortical surface area (Van Essen et al., 2019), both the shapes of these landmarks and the number of vertices within each landmark might show individual differences. We used participant-specific landmark label files to locate the participant-specific vertices belonging to each landmark and participant-specific parcellation to locate the vertices within each parcel.

Geodesic distance along the “midthickness” of the cortical surface (halfway between the pial and white matter) was calculated using the Connectome Workbench software with an algorithm that measures the shortest path between two vertices on a triangular surface mesh (Mitchell et al., 1987; O’Rourke, 1999). This method returns distance values independent of mesh density. Geodesic distance was extracted from surface geometry (GIFTI) files, following surface-based registration (Robinson et al., 2014). To ensure that the shortest paths would only pass through the cortex, vertices representing the medial wall were removed from the triangular mesh for this analysis.

We calculated the minimum geodesic distance between each vertex and each landmark. Specifically, for the central sulcus, we calculated the geodesic distance between vertex i outside the central sulcus and each vertex within it (defined for each individual). We then identified vertex j within the central sulcus closest to vertex i and extracted this value as the minimum geodesic distance for vertex i to this landmark. To compute the minimum geodesic distance for parcel k to the central sulcus, we averaged the minimum distance across all the grayordinate vertices in parcel k to the vertices within the central sulcus. The same procedure was applied to calculate minimum geodesic distance between each parcel and all three sensory-motor landmarks (central sulcus, temporal transverse sulci, and calcarine sulcus). From these three minimum geodesic distances, we selected the lowest distance value (i.e., the closest landmark to parcel k) as the global minimum distance to sensory-motor regions for parcel k. We then averaged the mean minimum distance of all parcels within each network for each participant and then sorted the networks by the mean minimum distance across participants. Finally, we examined whether mean minimum distance of FPCN-A and FPCN-B were different by performing a paired t test.

Anatomical hierarchy—myelin content and cortical thickness

We measured myelin, which is a noninvasive and valid proxy for anatomical hierarchy, and captures this hierarchy more effectively than cortical thickness (Burt et al., 2018). Gray matter myelin content can be measured via the cortical T1w/T2w map, a structural neuroimaging map defined by the contrast ratio of T1- to T2-weighted (T1w/T2w) magnetic resonance images (Glasser and Van Essen, 2011; Glasser et al., 2013, 2016; Salimi-Khorshidi et al., 2014).

Human T1w/T2w maps were obtained from the HCP in the surface-based CIFTI file format. To generate these maps, high-resolution T1- and T2-weighted images were first registered to a standard reference space using a state-of-the-art areal-feature-based technique (Glasser and Van Essen, 2011; Glasser et al., 2013, 2016).

Each participant has a myelin map, with each vertex having a myelin value (i.e., T1w/T2w ratio). We calculated the myelin value of each parcel by averaging the myelin values of all vertices within the parcel for each participant, using the individual-specific parcellation. Similarly, we calculated the myelin value of each network by averaging the myelin values of all parcels within each network for each participant and then sorted the networks by the mean myelin value across participants. Finally, we examined whether the myelin values of the FPCN-A and FPCN-B differed by performing paired t tests.

Cortical thickness, which coarsely tracks changes in cytoarchitecture and myelin content, can be seen as a pragmatic surrogate for cortical microstructure. Therefore, we also examined cortical thickness of each network which measures the width of gray matter across networks. The procedure was similar to that for myelin mapping, except we used the cortical thickness maps and extracted the cortical thickness of each vertex for each participant.

Function hierarchy—principal connectivity gradient analysis

To examine the relative position of networks on the principal gradient axis of intrinsic connectivity, we performed dimension reduction analysis on the resting-state functional connectivity matrix of the HCP dataset. First, we calculated the resting-state functional connectivity for each run of each participant using the method mentioned below, Constructing fMRI functional connectivity matrices. We then averaged these individual connectivity matrices to generate a group-averaged connectivity matrix. We used the BrainSpace Toolbox (Vos de Wael et al., 2020) to extract 10 group-level gradients from the group-averaged connectivity matrix (dimension reduction technique, diffusion embedding; kernel, none; sparsity, 0.9), following the methodology of previous studies (Mckeown et al., 2020; Wang et al., 2020). Using identical parameters, gradients were then calculated for each individual using their average 400 × 400 resting-state functional connectivity matrix across four runs. These individual-level gradient maps were aligned to the group-level gradient maps using Procrustes rotation, facilitating the comparison between the group- and individual-level gradients (N iterations = 10). This analysis produced 10 group-level gradients and 10 individual-level gradients for each participant, explaining maximal whole-brain connectivity variance in descending order. For 238 out of 245 participants, the first gradient, which explained the maximal variance, was the principal gradient which captures the separation between unimodal and transmodal regions (Margulies et al., 2016). Next, we averaged the first gradient values of all parcels within each network for each participant for whom the first gradient was the principal gradient and then sorted the networks by the mean gradient values across participants. Finally, we examined whether the gradient values of FPCN-A and FPCN-B differed by performing paired t tests.

Comparing the evolutionary expansion and cross-species similarity between FPCN-A and FPCN-B

We used the evolutionary expansion map and cross-species similarity map provided by Xu et al. (2020), available at https://github.com/TingsterX/alignment_macaque-human. The surface areal expansion map was calculated as the ratio of the human area to the macaque area at each corresponding vertex on human and macaque surfaces (Xu et al., 2020). Cross-species similarity was calculated by comparing the whole-brain patterns of functional connectivity in macaques and humans (Xu et al., 2020). We extracted the value of each vertex and calculated the mean value for each parcel within the group parcellation. Finally, we compared the values between FPCN-A and FPCN-B by conducting independent t tests.

Feature extraction of the time-series data

We used the highly comparative time-series analysis toolbox, hctsa (Fulcher et al., 2013; Fulcher and Jones, 2017), to extract massive features from the time-series data. Using a time-series dataset, hctsa allowed us to transform each time series into a set of over 7,700 features, with each feature encoding a different scientific analysis method (Fulcher et al., 2013; Fulcher and Jones, 2017). The extracted features include, but are not limited to, distributional properties, entropy and variability, autocorrelation, time-delay embeddings, and nonlinear properties of a given time series (Fulcher et al., 2013; Fulcher, 2018). The hctsa feature extraction analysis was performed on the parcellated fMRI time series for each participant, each task, and each run separately. Following this procedure, we removed the outputs of the operations that produced errors and normalized the remaining features (∼6,900 features) across parcels using an outlier-robust sigmoidal transform. The normalized feature matrix (400 parcels × ∼7,000 features × 4 runs) was used for subsequent classification analysis and feature similarity analysis.

Classification analysis

To reveal network similarity in a data-driven approach, we performed classification analysis using the normalized feature matrix. This analysis was motivated by the observation that specific brain regions possess distinct features, equipping them to support different functions, while regions with similar features are likely to support similar functions, irrespective of whether they exhibit positive or negative functional connectivity. For instance, early visual areas and DMN are situated at opposing ends of the timescale hierarchy. The former is characterized by the shortest timescale, marked by rapid temporal autocorrelation decay, while the latter has the longest, characterized by gradual autocorrelation decay (Honey et al., 2012; Simony et al., 2016; Blank and Fedorenko, 2020; Ito et al., 2020; Raut et al., 2020). Consequently, neural representations in early visual areas are minimally influenced by prior knowledge, whereas DMN regions are significantly swayed by it (González-García et al., 2018). In addition, regions with more similar features often support parallel functions. For example, while the two large-scale networks in the transmodal regions, FPCN and DMN, typically exhibit negative functional connectivity, they share attributes like long timescales, indicating they can process inputs over extended periods. Thus, their neural representations are influenced by prior knowledge (González-García et al., 2018) and goal states maintained over time (Wang et al., 2021). These observations suggest that feature similarity can provide valuable information about the functional similarity of networks beyond functional connectivity. Furthermore, features can vary depending on tasks and brain states (Golesorkhi et al., 2021; Wolff et al., 2022), indicating their potential to reflect variations in network interaction across tasks.

Constructing feature similarity matrices

To investigate how subnetworks of FPCN interact with DAN and DMN, we calculated both task and resting-state feature similarity. For each normalized feature matrix, we calculated Pearson’s correlation coefficients and transformed them to Fisher z values to represent the pairwise feature similarity between the time-series features of all possible combinations of brain parcels. This process resulted in a 400 × 400 feature similarity matrix for each individual and each run, representing the strength of the similarity in the local temporal fingerprints of brain areas. Finally, we averaged the estimates of feature similarity both within networks and between pairs of networks to construct a network-by-network feature similarity matrix. The same method was used to calculate the resting-state feature similarity using the HCP dataset and construct a network-by-network feature similarity matrix.

Constructing fMRI redundancy matrices

To investigate how subnetworks of FPCN interact with DAN and DMN, we employed the concept of time-delayed mutual information (TDMI), a specialized version for analyzing temporal dependencies in time-series data. TDMI extends mutual information by introducing time delays, assessing the mutual information between a variable's current state and its past states with a time lag. This extension enables us to understand how past information predicts the current state over time. To calculate redundancy, we utilized the recently developed integrated information decomposition approach, as introduced by Mediano et al. (2021) and Luppi et al. (2022). Integrated information decomposition expands on the partial information decomposition theory, breaking down TDMI into redundant (information provided by both sources), unique (by one source but not the other), and synergistic information (jointly by their combination) concerning both past and present states of the involved regions. We specifically focused on temporally persistent redundancy, aiming to quantify the extent to which distinct brain regions redundantly carry information about the brain's future trajectory, serving as an indicator of functional interaction. This measure was calculated using the Gaussian solver implemented in the JIDT toolbox, based on their hemodynamic response function-deconvolved BOLD signal time series. We calculated the redundant interaction for each pair of brain regions, resulting in a 400 × 400 redundancy matrix for each participant, per task and per run. Finally, we averaged the estimates of redundancy within networks, and between pairs of networks, to construct a network-by-network redundancy matrix, similar to the one constructed from the feature similarity data.

While the organization of redundancy shares similarities with traditional functional connectivity, as pairs of regions with more correlated time courses are more likely to provide redundant information, they are not identical. This is evident by the moderate correlation between the functional connectivity matrix and the redundancy matrix (mean = 0.62; SD = 0.24; Luppi et al., 2022). Notably, different ensembles of regions can transiently shift from being redundancy dominated to synergy dominated (Varley et al., 2023). This shift suggests that two regions may have strong functional connectivity while not consistently sharing information throughout the entire observation period. This distinction between measures also indicates that persistent redundancy offers valuable insights into the dynamics of information sharing between brain regions that are not fully captured by functional connectivity.

Constructing fMRI functional connectivity matrices

To investigate how FPCN subnetworks interact with DAN and DMN, we further calculated both task and resting-state functional connectivity. We opted not to use the traditional psychophysiological interaction method for measuring task-state functional connectivity due to its potential to inflate activation-induced task-state functional connectivity, which may identify regions that are active rather than interacting during the task (Cole et al., 2019). Since task activations can spuriously inflate task-based functional connectivity estimates, it is necessary to correct for task-timing confounds by removing the first-order effect of task-evoked activations (i.e., mean evoked task-related activity, likely active during the task) prior to estimating task-state functional connectivity (likely interacting during the task; Cole et al., 2019). Specifically, we fitted the task timing for each task using a finite impulse response model, a method that has been widely used and shown to reduce both false positives and false negatives in functional connectivity estimation (Norman-Haignere et al., 2012; Cocuzza et al., 2020; McCormick et al., 2022). In semantic tasks, approximately five time points were modeled for each trial. For the nonsemantic task, roughly 2.5 time points were modeled for each trial.

Following task regression, we demeaned the residual time series for each parcel and quantified the functional connectivity using Pearson’s correlation for each participant, task, and run. Pearson’s correlation coefficients might be inflated due to the temporal autocorrelation in task fMRI time-series data (Allefeld et al., 2008). To address this, we corrected Pearson’s correlation using a novel correction approach, xDF, which accounts for both autocorrelation within each time series as well as instantaneous and lagged cross-correlations between the time series (Afyouni et al., 2019). This method provides an effective degrees of freedom estimator that addresses cross-correlations, thereby preventing inflation of Pearson’s correlation coefficients. Our goal was not to remove temporal autocorrelation but to enhance the precision and reliability of correlation estimates, reducing false functional connectivity between regions. We calculated xDF-adjusted z-scored correlation coefficients to assess the interregional relationships in BOLD time series, resulting in a 400 × 400 functional connectivity matrix for each participant, task, and run. Finally, we averaged these functional connectivity estimates within networks, and between pairs of networks, to construct a network-by-network functional connectivity matrix. The same method was used to calculate the resting-state functional connectivity of the HCP dataset and construct a corresponding network-by-network functional connectivity matrix, except without the task regression step.

Comparing feature similarity, redundancy, and functional connectivity difference between networks

To investigate whether FPCN-A showed greater feature similarity with DAN than FPCN-B did in each task, we calculated the average feature similarity between FPCN-A and DAN, and the feature similarity between FPCN-B and DAN, respectively, across all runs per participant per task, and then conducted paired t tests to compare the levels of similarity between these network pairs for each task. Similarly, we examined the feature similarity between FPCN-B and DMN versus FPCN-A and DMN. We conducted FDR correction at p = 0.05 to control for multiple comparisons. We investigate whether FPCN-A showed stronger redundancy and functional connectivity with DAN than FPCN-B did. We also examined whether FPCN-B showed greater redundancy and functional connectivity with DMN than FPCN-A did. The procedure was as above, except we extracted the redundancy and functional connectivity matrix, respectively.

We further investigated whether the task context influenced the feature similarity difference between network pairs using the maximum/minimum permutation test. We might expect that feature similarity between FPCN-A and DAN-A versus FPCN-A and FPCN-B would be greater in spatial working memory task than that in the association task because FPCN-A might shift more toward DAN in the spatial working memory task that requires an external goal but no memory retrieval. To examine this possibility, we calculated the mean feature similarity difference between FPCN-A and DAN-A versus FPCN-A and FPCN-B for each task and calculated the mean feature similarity difference between each task pair. To test for statistical significance, we permuted the task label 10,000 times; we then calculated the mean feature similarity difference between these two tasks to build a null distribution for each task pair. Since we included multiple task pairs, we used the permutation-based maximum mean feature similarity difference and minimum mean feature similarity difference values in the null distribution for each task pair to control the FWE rate (p = 0.05, FWE corrected). To evaluate significance, if the observed mean difference value was positive, we counted the percentage of times that mean difference values in the maximum null distribution were greater than the observed “true” mean difference values; in contrast, if the observed mean difference value was negative, we counted the percentage of times of mean difference values in the minimum null distribution were less than the observed “true” mean difference values.

Similarly, we investigated whether the task influenced the redundancy and functional connectivity difference between network pairs. The procedure was as above, except we extracted the redundancy and functional connectivity from the network-by-network redundancy and functional connectivity matrix, respectively. We conducted FDR correction at p = 0.05 to control for multiple comparisons.

Task fMRI univariate analysis

To reveal the functional differentiation between FPCN-A and FPCN-B, we examined whether they showed similar or opposite activation patterns. We identified regions that were activated or deactivated in the tasks by building a general linear model. We also examined regions where neural responses were modulated by task difficulty. For semantic tasks, we included one task mean regressor and one demeaned parametric regressor of semantic rating. We examined how the neural responses were negatively modulated by feature similarity rating in the matching trials and positively modulated in the nonmatching trials. For the semantic association task, we examined how neural responses were negatively modulated by semantic association strength rating in related trials and positively modulated in nonrelated trials. For the nonsemantic tasks, we included two regressors—easy and hard conditions—to reveal the regions showing greater activation in the hard than easy conditions. We also modeled incorrect trials as regressors of no interest in all tasks.

We extracted the β value for each parcel in these task conditions and tested whether they were significantly activated (i.e., above zero) or deactivated (i.e., below zero) relative to implicit baseline (i.e., fixation period) to explore the task mean effect. Then we considered parcels that showed stronger or weaker activation in the hard than that in the easy condition in the spatial working memory and math tasks and that the activations were positively modulated by semantic difficulty. These parcels were thought to support general executive control. Fixed-effects analyses were conducted using nilearn (Abraham et al., 2014) to estimate average effects across runs within each subject for each parcel. Then we conducted one-sample t tests to assess whether the estimated effect size (i.e., contrast) was significantly different from zero across all subjects. We applied FDR correction at p = 0.05 to control for multiple comparisons. Finally, we identified the network to which each parcel belonged according to Kong et al. (2021).