Here we provide an outline of our experimental framework. The current study used the same dataset analyzed in our previous study [22], specifically designed for Granger causality and behavior joint analysis:

We preprocessed the images and carried out a GLM analysis to identify the regions of interest (ROIs).

i. We extracted and preprocessed the fMRI time series of each ROI using Granger causality analysis (GCA) to assess the ROIs’ directional connections.

The third and fourth steps aimed to intuitively depict the task-positive networks’ functional hierarchy during this visual-spatial attention task.

The detailed information of the experiment and analysis is provided below.

The current study used the same dataset analyzed in our previous study [22]. Specifically, twenty young, healthy right-handed human participants with normal or corrected-to-normal vision participated in the experiment. All participants had no history of taking medicine, psychiatric illness, or any brain surgery that might affect their central neural system. Because the attention task was very demanding and we want to observe the natural fluctuation of behavior mainly related to the varying attention level instead of other factors, the participants had to be well trained and well prepared for the six runs of task-scanning. Each participant’s training procedure included multiple out-scanner training sessions and one in-scanner training in 2-4 days before the MRI scanning day. Each of these training sessions lasted 0.5~1 hour. All participants underwent a screening and a short warm-up session before the MRI sessions. Only 13 of them finished the whole six runs. One of them was later excluded from the analysis because of severe image artifacts, which reflected shadow artifacts on both sides of the brain. The included participants were 24 $\pm$ 1.52 years old (8 females and 4 males). All participants signed a written informed consent beforehand, which abided by the Helsinki Declaration, and all research activities were authorized by the Brain Imaging Center at Beijing Normal University.

The current study adopted a mixed blocked/event-related design. The task protocol contained six runs, each with four blocks balanced in an ABBA or BAAB arrangement to counterbalance the temporal confound. The experimental timeline is schematically illustrated in Fig. 1A. Each block lasted one minute and was followed by a 20 second fixation period. The attention (A) blocks and passive-view (B) blocks shared the same timeline and stimuli except for the color of the crosshair at the fixation point (light red and light green, balanced across participants). In each trial of attention (A) blocks, participants were cued to direct and maintain covert attention to the left or right hemifield. Following a 2500 ms delay, a standard or a target stimulus of 100 ms in duration appeared either in the attended hemifield (valid trial) or the unattended hemifield (invalid trial). The standard stimulus was a circular checkerboard. The target stimulus was also a circular checkerboard but slightly smaller than the standard stimulus (10% smaller in radius). The standard stimulus appeared 80% of the time with 50% validity, and the target stimulus appeared 20% of the time with 50% validity. The trials were pseudo-randomly arranged so that the valid trials and the invalid trials in each block were evenly matched. Participants were required to make a speeded keypress response only to the valid targets. Since the proportion of valid target stimuli to all stimuli was small (10%), motor processing activation was weak at the block level. It would not affect the attention task’s activation results. On the other hand, it would avoid motor processing related component which might contaminate the time series for GC analysis at the block level. In passive-view (B) blocks, the stimulus presentation schedule remained the same, but neither attention nor response was required. The participant just maintained fixation. Each block lasted ~60 s (15 trials), with 20 s fixation periods inserted between successive blocks. The participants had to be well trained and well prepared to maintain their performance during the six runs of task-scanning, which lasted about an hour. The training procedure of each participant included multiple out-scanner training sessions and one in-scanner training in the 2-4 days before the scanning day. Each of these training sessions lasted 0.5~1 hour. The participants were fully instructed and went through the out-scanner training to get familiar with the task.

Fig. 1.

Experimental paradigm: region of interest (ROI) defined using task activation and mean ROI blood oxygen level-dependent (BOLD) changes. (A): An example of the stimulus and timeline for one of the trials in an attention block. (B): The task-positive networks activated in the current task (T $>$ 5.20, P 0.002, FDR corrected) are illustrated using different colors on a three-dimensional brain surface template (Caret 5 visualization software). (C): A constant elevation in BOLD change can be observed in the attention condition for all task-positive ROIs. The gray area denotes the block period.

All MRI data were acquired using a 3T Magnetom Trio whole-body MRI system (Siemens AG, Erlangen, Germany) at Beijing Normal University MRI Center. The functional scanning was performed using a T2*-weight echo-planar imaging sequence (echo time, 30 ms; repetition time, 2000 ms; flip angle, 90${}^{\circ }$) with 33 axial slices in each volume (field of view, 200 $\times$ 200 mm; matrix size, 64 $\times$ 64; slice thickness, 3.60 mm; voxel size, 3.13 mm $\times$ 3.13 mm $\times$ 3.60 mm). For each participant, 300 whole-brain resting-state volumes (10 minutes) were recorded but were not included in the present study. After the resting-state session, 180 whole-brain volumes were recorded in each of the six task runs. Participants were fully informed about the requirement of keeping head still in the scanner, and the head movement was strictly controlled using a set of memory foam adapters. Anatomic images were acquired with a T1-weighted 128 slice MPRAGE sequence (repetition time, 2530 ms; echo time, 3.39 ms; flip angle, 7${}^{\circ }$; inversion time, 1100 ms; voxel size, 1 mm $\times$ 1.33 mm $\times$ 1 mm).

The fMRI data were preprocessed using SPM8 software (http://www.fil.ion.ucl.ac.uk/spm). The preprocessing protocol included slice timing, motion correction, anatomical co-registration, normalizing to a Montreal Neurological Institute (MNI) space (voxel size, 3 mm $\times$ 3 mm $\times$ 3 mm), and spatial smoothing using an 8 mm FWHM Gaussian core. The scrubbing procedure, global scaling, and regression of white matter signals, cerebrospinal fluid (CSF), and the 24 head motion parameters were applied to reduce nuisances further. The 24 head motion parameters refer to Friston’s 24-parameter model of head motion [31, 32], which incorporates the 6 standard head motion parameters, the derivative of the standard motion parameters to account for a one-frame delay in the effect of motion on the blood oxygen level-dependent (BOLD) signal [33], and the 12 corresponding squared items [31].

The forgoing preprocessed data with spatial smoothing were fed to a GLM for activation analysis and ROI selection. The data were also preprocessed using the above steps except spatial smoothing. These none-spatial-smoothing versions of preprocessed data were fed to GCA for network analysis.

To define the ROIs for network analysis, we conducted a GLM analysis using SPM8 software. For first-level analysis, the regressor is generated by convolving the rectangular function representing the block sequence with a canonical hemodynamic response function (HRF). The individual activation maps were generated using the contrasts of attention condition against passive-view condition. For second-level random effect analyses, the individual contrast maps were fed to a one-sample t-test to yield a group-level activation map. False discovery rate (FDR) control was applied to correct for multiple comparisons (t $>$ 5.20, P 0.002, FDR-corrected). The ROIs were generated by intersecting the group-level activation map with spheres of 5 mm in radius centered at the voxels with maxima local t values.

To avoid false-positive results which may confound the ROI selection and to make sure the ROIs defined matched the well-proposed CON, DAN, and VAN network, we first compared the GLM activation results with the spatial pattern of the region associated with “spatial attention” according to the online meta-analysis (https://neurosynth.org/, 147 studies, uniformity test, P 0.01, FDR-corrected). Second, we compared the GLM activation results with the classical DAN, VAN, and CON spatial patterns reported in previous literature [6, 7, 17].

To elucidate the functional hierarchy of the task-positive networks, a directed network model would be meaningful. We chose GCA, which is widely used in examining the directed influence between time series to accomplish this goal [15, 22, 34, 35, 36]. The fundamental idea of GC is if the history of time series X facilitates the prediction of time series Y’s future, then we say there is a Granger causal influence from X to Y [37]. GC value of X$\to$Y indicates the strength of the information flow from X to Y.

One of the mathematical realizations of estimating the Granger causality is comparing the autoregressive (AR) prediction performance of the univariate prediction and multivariate regression (MVAR) performance. For example, the GC from X to Y can be defined as

(1)$$F_{X\to Y}=\ln \frac{{\mathrm{\Gamma }}_1}{{\mathrm{\Gamma }}_2},$$
where $\mathrm{\Gamma }1$ denotes the variance of the residue of the univariate AR model fitting and $\mathrm{\Gamma }2$ the covariance of the residue matrix of the bivariate (or multivariate) AR model fitting. More details of the mathematical realizations were provided in Wen et al., 2013 [38].

Assuming X is one ROI and Y being another, our GC calculation between X and Y has three major steps. 1) Extracting the time series: Time series of each voxel in X and Y were extracted from the preprocessed functional images without spatial smoothing. The time series were then converted to percentage BOLD signal changes by subtracting the mean signal value during the inter-block baseline fixation period and then dividing the difference by the mean value. 2) Making the time series stationary and zero-mean: For each voxel, each block’s percentage change signals were averaged within each condition (attend or passive-view) to yield the block-wise BOLD response. For the attention condition, the block-wise response was subtracted from the percentage BOLD change in each block and each voxel to yield residual BOLD time-series. For each block, the residual BOLD time-series’ first five-time points were discarded to eliminate the transient effects. The temporal mean of the remained time points of each block was removed to meet the zero-mean requirement assumed by autoregressive model estimation in GCA [34]. 3) Calculating GC values: each voxel in X were paired with a voxel in Y, and GC values were calculated for each voxel pair and averaged across all pairs to yield overall GC values between X and Y, including GC of X$\to$Y and GC of Y$\to$X. Based on the Bayesian information criterion, the order of the AR model was determined to be 1 [39, 40, 41, 42].

Investigating the change of GC values across different conditions on the group level is more meaningful than merely observing the raw GC values at the individual level. The former may reveal the cognitive significance of the directional connections and to mitigate the confounds caused by individual differences and noise background [15, 29, 39]. Accordingly, we employed a framework to correlate behavioral performance with the causal influence between brain regions. The framework formed the foundation of the network construction based on behavior-correlated inter-ROI connections. Specifically, for each subject, the GC values and behavioral performance (either accuracy or response time (RT)) for each attention block were converted into z-scores. For the convenience of combining the two behavioral measures, we multiplied the RT z-score by -1 so that larger scores for both measures indicated better performance. The attention blocks were then sorted according to z-scores and assigned to 10 levels, each containing three neighboring blocks. The sorting assured that the first level denoted the worst performance (lowest accuracy or longest RT), and the last level, the best performance (highest accuracy or shortest RT). For each voxel pair and level, the three blocks’ GC values were averaged to represent the GC strength corresponding to the performance at that level. Spearman’s rank correlation analysis was then performed to examine the relationship between level-GC and level-performance, assessed using either accuracy or RT.

Both activation analysis and GCA were calculated at the block level. Therefore, all events were included in GLM and GC estimation. When calculating block-level mean RT, the responses to standard stimuli (false alarm) were excluded. However, the RT of the trials of the missing target was set at 2 times the participant’s average hit RT because previous behavioral studies have shown that in visual-motor experiments, RT in target-absent trials is usually twice as long as the average RT of hit trials [43, 44].

To identify whether a directional connection, for example, A$\to$B, was behaviorally significant, we combined the GC-accuracy and GC-RT correlation results. Specifically, if the GC value of A$\to$B and accuracy or RT (or both) were significantly positive at P 0.05, A$\to$B was behavior-enhancing; in other words, a stronger A$\to$B was associated with better performance. By contrast, if the GC value and accuracy or RT (or both) were significantly negative at P 0.05, A$\to$B was behavior-degrading; in other words, a stronger A$\to$B was associated with worse performance. If a correlation was significantly positive for one behavioral measure but significantly negative for the other, the directional connection’s role was considered ambiguous. We observed no ambiguous connections in our study.

The GC values of all cross-network ROI-pairs were calculated and averaged to yield the network-pairs’ GC values. For example, let A and B be two separate networks and a${}_i$ be the ith ROI in A, and b${}_j$, the jth ROI in B. We calculated the GC values of a${}_i\to$ b${}_j$ for all i-j combinations and those in the opposite direction to yield the inter-network GC of A$\to$B and B$\to$A, respectively. The inter-network GC values were then correlated with the behavioral performance using the same conventions introduced above to identify behavior-significant interactions on the inter-network level.

The GLM analysis yielded an activation map highlighting three major task-positive networks: CON, DAN, and VAN (t $>$ 5.20, P 0.002, FDR-corrected). The activation of the CON comprised dACC and bilateral AI [7, 26]; the activation of the DAN comprised the bilateral FEF and IPS; and the activation of the VAN comprised the right anterior MFG (raMFG), the right posterior MFG (rpMFG), and TPJ (Fig. 1B) [6], were all observed. The coordinates of the center voxels of the ROIs are listed in Table 1. The mean percentages of BOLD signal changes extracted from these regions showed a constant elevation in attention condition (Fig. 1C). Among the 10 ROIs, the dACC and right AI (rAI) had the most robust signal increases in the attention condition. Besides, comparing the activation results with the meta-analysis of spatial attention tasks, we found that the activation pattern was consistent with the classical networks significantly activated in previous spatial attention tasks (Fig. 2). The group activation results were also consistent with the classical DAN, VAN, and CON spatial patterns reported in previous literature [6, 7, 17] (Fig. 1B).

Fig. 2.

Comparing the group activation results of the current spatial-attention task with the spatial pattern of the regions associated with “spatial attention” according to the online meta-analysis (https://neurosynth.org/, 147 studies, uniformity test, P 0.01, FDR-corrected).

To assess the causal connections between different networks’ ROIs and the behavioral significance of those connections, we calculated the GC values between the ROIs and correlated these values with behavioral performance (see Methods section for details). Because we previously investigated the causal interactions between DAN and VAN, we did not repeat these calculations [22]. The current work mainly focused on the connections between CON and DAN/VAN, which were not assessed in our previous studies and were less often discussed in other studies. The causal connections significantly correlated with behavioral performance are shown in Fig. 3A. Generally, most of the behaviorally significant causal connections from the CON ROIs to the DAN and VAN ROIs were behavior-enhancing. In the opposite direction, the behavioral significance became inconsistent across the ROI-pairs. The connections from the bilateral FEF or raMFG to the dACC and those from the left IPS (lIPS) or the right temporoparietal junction (rTPJ) to the rAI were behavior-degrading, while most of the connections from the bilateral FEF or rIPS to the AI were behavior-enhancing. Fig. 3B illustrated an example of identifying behaviorally significant connections between dACC and rFEF in which the causal influence of dACC$\to$rFEF was significantly and positively correlated with the RT score, while that in the opposite direction was significantly and negatively correlated with the RT score.

Fig. 3.

Behaviorally significant causal connections on the both region of interest (ROI) and network levels. (A): A schematic summary of the behaviorally significant connections between CON and DAN/VAN ROIs. Each cell refers to a directional connection. The abbreviations in the cell indicate the behavior measurements significantly correlated with the causal influence. AC, accuracy; RT, response time. Red color denotes positive correlation (behavior enhancement), while blue denotes negative correlation (behavior degradation). (B): Inter-network causal influences as functions of accuracy. The conventions are the same as above. (C): A schematic depiction of the behaviorally significant causal interactions between the CON and the other two task-positive networks. Red color denotes behavior-enhancement, and blue denotes behavior-degradation. The previously reported [22] causal interactions between DAN and VAN and their behavioral significance were also displayed.

The inter-ROI causal influences were averaged within four categories CON$\to$DAN, DAN$\to$CON, CON$\to$VAN, and VAN$\to$CON-to assess network-level interactions. Analysis of the correlations of the inter-network interactions with the behavioral performance showed that stronger CON$\to$DAN and CON$\to$VAN were associated with higher accuracy (dACC$\to$AI, R = 0.92, P = 0.00047, uncorrected; dACC$\to$VAN, R = 0.66, P = 0.044, uncorrected; scatter plot in Fig. 2B, left and middle panels). By contrast, a stronger VAN$\to$CON was associated with lower accuracy (R = -0.82, P = 0.0068 uncorrected; see Fig. 2B, right panel). The results showed different behavior consequences between DAN$\to$CON and VAN$\to$CON. The pattern of behaviorally significant network interactions between CON and DAN/VAN was intuitively depicted in Fig. 3C.

During attention blocks, the participants frequently directed attention according to cue and maintained attention for seconds (top-down attention control) and frequently be affected by invalid stimulus (bottom-up interference), which cause interference even reorienting. The three processes were proposed to activate CON, DAN, and VAN. By contrasting the attention blocks against the passive view blocks, our primary purpose is to validate the three networks’ activation and provide reliable ROIs for GC analysis rather than inspecting complex visual-spatial components at the trial level, which is not the focus of the current study and had been done by many previous studies. To avoid defining ROIs based on false-positive activation results, we carefully compared the current study’s activation map with those independently reported in previous literature, especially those dedicated to defining the classical DAN, VAN, and CON according [6, 7, 17, 18, 45]. Our results of the activation were consistent with the previous studies. Further comparison with the meta-analysis could be regarded as a double-check of the ROI selection efficacy (Fig. 2).

Visual attention is considered to be controlled by the frontoparietal attention system, including the DAN and VAN. The CON is thought to integrate internal and extra personal information to regulate other brain areas and guide behavior, which is crucial for task-control [8, 46, 47]. On the other hand, the DAN initiates and maintains goal-directed top-down signals important for attention control [17, 48, 49]. The findings that stronger CON$\to$DAN leads to better behavioral performance suggests a network-level mechanism in which top-down control signals pass from the CON to the DAN. Consistent with previous fMRI studies, the top-down influences of the CON contribute to goal-oriented behavior in task-relevant information [25, 50].

The finding that stronger CON$\to$VAN was associated with behavior enhancement implies the presence of top-down control signals from the CON to VAN to filter behaviorally irrelevant input and enhance the efficacy of attention processing in the attended domain [17, 45]. On the pathological level, attention deficit hyperactivity disorder (ADHD) patients showed an abnormally heightened processing of irrelevant information, while interference control was continuously compromised. Moreover, ADHD patients have shown reduced activation in dACC than in healthy controls [51, 52, 53].

Our results also showed bidirectional interactions between VAN and CON. Emerging evidence has shown that attention reoriented to a new source output from the VAN interrupts ongoing selection in the DAN, which, in turn, shifts the attention toward the novel object of interest [15, 22]. Our results showed that a stronger VAN$\to$CON was associated with lower accuracy, reflecting the interference effect of VAN. The result, together with our previous finding [22], indicates that VAN not only sends bottom-up interference to DAN but also to CON, disrupts the processing conducted by CON in attention-demanding tasks, thus potentially acting as a “circuit breaker”.

Our results demonstrated the relationship between CON and the frontoparietal attention system, and we further proposed a functional hierarchy model in the visual attention task. Our previous study emphasized the relationship between DAN and VAN [22]. Combined with the present results, the CON appeared to be the highest-ranking network in the hierarchy, suggesting that it may regulate the frontoparietal attention network by transmitting top-down signals to accomplish attention goals. Generally, the top-down processing from the CON to the DAN specializes in selecting and linking stimuli and responses to guarantee attention performance. In contrast, the CON regulation of top-down signals prevents interference by VAN to ignore distraction. VAN may also break the attentional set maintained by CON to enable attentional reorienting. These findings add to our understanding of the brain’s functional hierarchy from the perspective of network connectivity.

Granger causality analysis was employed in the present study to construct a functional hierarchy network employing BOLD fMRI data during the attention task. Granger causality analysis is an exploratory approach, which is not restricted to the preselection of interacting regions and assumptions about the structure and the direction between the brain regions. Therefore, GCA (unlike SEM and dynamic causal modeling (DCM)) does not appear to have the issues of the model’s misspecification and inaccurate results [39]. While the applicability of Granger causality to fMRI data is debated [54], it is noteworthy that from a statistical point of view, if a time series is analyzable by functional connectivity measures such as temporal correlation and coherence [55], it is analyzable by Granger causality. Recent work shows that both resting-state and task-state fMRI data are well described by autoregressive models, which are the basis for deriving Granger causality [9, 40, 41, 56].

Admittedly, GCA is suggested vulnerable to noise [54] and sometimes limited by the AR model order. However, by applying GCA, the information streams between ROIs would be more directly measured and intuitively displayed. When ROI A is influencing ROI B, which means there is a steady information stream from A to B. Thus, when predicting ROI B’s future activity, the activity of ROI A would contribute to the prediction in addition to using the history of ROI B merely. By applying GCA, such an information stream could be quantified by the coefficient of the history of A on future B, which reveals the synchrony of activation and reflects dependencies between ROIs. In the present study, GCA is used to build the hierarchy of task-positive networks, which not only depicted the organization of the structure but also demonstrated the information streams and provided the bases for further analyses on the functions of the ROIs (i.e., target brain regions) in a visual-spatial attention task.

It is worth noting that avoiding false-positive results in the first step of ROI definition using GLM activation result is crucial for the subsequent analysis [57]. Therefore, we did not define the ROI only by considering the task GLM activation results based on the current dataset. In practice, we carefully compare the current study’s activation map with those independently reported in previous literature and the activation map by meta-analysis (see Method). As mentioned, our activation results were consistent with the previous studies in the literature and the meta-analysis, which excluded false-positive confound.

Besides, our GC analysis did not rely on the same information that the GLM analysis relied on. Traditionally, we represent the regional BOLD activity as y(t) = x(t) + $ϵ$(t), where x(t) denotes the task-evoked amplitude change of the BOLD signal, which is mainly considered in GLM analysis, and $ϵ$(t) is considered as noise. Further, according to previous studies [29, 58], $ϵ$(t) includes the so-called ongoing activity s(t) and the true noise n(t), and s(t) comprises subtle temporal information that may infer the information transmitting between the brain regions. Therefore, the representation of BOLD activity can be reframed as y(t) = x(t) + s(t) + n(t). The GC analysis only focused on $ϵ$(t), independent of the task-evoked component x(t). The results of the GLM, which is a univariate analysis, does not affect the strength of GC, which examines the temporal interdependence between ongoing activities of two regions and falls into the scope of multivariate analysis.

Although the number of participants was a limitation of our study, we adopted a mixed blocked/event-related design, which allowed us to reveal the target networks’ hierarchy. The experiment was designed specifically for GC-behavior joint analysis, which is suitable for evaluating the behavior significance of effective connectivities. First, the experiment had to contain a sufficient number of blocks for a single condition (12 attention blocks in six runs in our case); the more blocks, the merrier. Second, each block needed to be sufficiently long for GC estimation; the longer, the better. Third, because the attention task was very demanding and we wanted to observe the natural fluctuation of behavior mainly related to the varying strength of GC rather than other factors, the participants had to be well trained and well prepared to maintain their performance during the six runs of task-scanning, which lasted for approximately 1 hour. The first and the second methodological demands linked to the third practical issue greatly limited our ability to employ this kind of design in massive recoding studies.

The second limitation is that we could only perform GC-behavior joint analysis on the block level signals, which contain frequent trial events. By itself, the current GC-behavior joint analysis cannot distinguish whether the control signal is related to the task level control or the trial level control. Although CON is considered to perform at the task-level [7], while the interaction between DAN and VAN is discussed at the trial-level [6, 17], and it is possible that the causal influence exerted from CON may reflect task-level control signals, while those from DAN and VAN may reflect trial-level control signals, a more sophisticated framework that could elaborate the control signals on trial-level is still needed in future studies.

It is worth noting that the current study only considered a specific visual spatial-attention task paradigm. Whether this functional hierarchy survives other attention paradigms such as visual feature attention paradigm, auditory attention paradigm, or even other attention-demanding paradigms remains unclear. Therefore, studies using more sophisticated designs, a larger sample size, and more attention paradigms should be carried out to elucidate further the generality of the functional hierarchy of CON, DAN, and VAN in the future.