Introduction

The thalamus has been proposed to be important for communication among cortical regions because of its anatomical organization, neural response properties, synapse types, and prominent efferent and afferent connectivity with the cortex (Jones, 2007; Nakajima & Halassa, 2017; Ouhaz et al., 2018; Saalmann & Kastner, 2011; Sherman, 2016; Sherman & Guillery, 2006; Theyel et al., 2010). Some thalamic nuclei, such as the lateral geniculate nucleus, are well known for relaying sensory input to the cortex (Sherman, 2017; Sherman & Guillery, 2013, 2014). Other higher-order nuclei have been proposed to serve as relay stations for trans-thalamic corticocortical communication, forming pathways parallel to the direct corticocortical pathways (Sherman, 2012; Sherman & Guillery, 2011). Importantly, some recent studies have shown that thalamic function might be more complex than simply relaying information (Halassa & Sherman, 2019; Wolff & Vann, 2019). For example, Schmitt et al. (2017) found that the mouse thalamus controls functional connectivity in some higher-order regions, such as the prefrontal cortex. In line with the findings in mice and monkeys, computational models of control by the thalamocortical circuit also have been proposed (Kao et al., 2020; Logiaco et al., 2019). Some studies have proposed that control mechanisms might rely on modulatory signals from the thalamus (Casas-Torremocha et al., 2019; Saalmann et al., 2012; Viaene et al., 2011). In summary, recent progress indicates a more complex role of the thalamus in cognition than previously thought (Halassa & Sherman, 2019).

The aforementioned evidence is mainly derived from postmortem studies and invasive animal studies on particular trans-thalamic corticocortical pathways. Whole-brain neuroimaging and studies on mental disorders in humans provided more insight into the organization and function of the thalamocortical pathways in human cognition (Simon et al., 2018; Zhang & Li, 2017). For example, noninvasive diffusion tensor imaging (DTI) of the human brain demonstrated prominent and parceled structural connections between certain components of the thalamus and cortical regions (Behrens et al., 2003; Bohsali et al., 2015). DTI in patients with Alzheimer’s disease and mild cognitive impairment suggested that the attenuated deactivation of the default mode network (DMN) in cognitive decline might be the result of impaired white matter circuits connecting the thalamus and cortex (Alderson et al., 2017). Lesion studies have shown that thalamic stroke impairs the performance of patients in prospective learning tasks, indicating the importance of the thalamus in memory processes (Cheng et al., 2010). Noninvasive functional magnetic resonance imaging (fMRI) also has been used to investigate thalamocortical interactions (Ide & Li, 2010) and how they are compromised in mental disorders (Ide et al., 2018; Simon et al., 2018; Zhang et al., 2016).

The aforementioned studies have shown that the thalamus seems to be involved in a variety of advanced cognitive control-related processes, such as learning, memory, task switching, attention (Baxter, 2013; Bradfield et al., 2013; Fresno et al., 2019; Funahashi, 2013; Habib et al., 2013; Le et al., 2020; Marzinzik et al., 2008; McAlonan et al., 2008; Phillips et al., 2019; Radanovic et al., 2003; Rikhye, Gilra, & Halassa, 2018a; Saalmann, 2014; Wolff & Vann, 2019), and computing of contextual signals that rapidly reconfigure task-relevant cortical representations (Rikhye et al., 2018b), which implies that the thalamus plays an important role in cognitive control.

However, postmortem and DTI studies have provided limited information about the cognitive functions of thalamic connections, even though this has become increasingly compelling in recent years. Invasive animal studies provide useful insights into how thalamic nuclei regulate cortical-cortical communications through specific cognitive control-related pathways, but comprehensive investigations into the interaction between the thalamus and the multiple cognitive control-related cortical regions have rarely been undertaken.

It is worth noting that many functional imaging studies provide a noninvasive and whole-brain perspective to study the function of the thalamus. Many studies have reported that thalamic activation is associated with multiple cortical areas related to cognitive control. Using an automated meta-analysis which collected human fMRI results from 598 studies associated with the term “cognitive control” (www.neurosynth.org), we found that bilateral middle thalamic regions are highlighted together with several cortical systems related to cognitive control, including: (1) the core control network (Corbetta et al., 2008; Dosenbach et al., 2006; Wen et al., 2018), which comprises a set of prefrontal regions, such as the dorsal anterior cingulate cortex (dACC) and the bilateral anterior insula (AI), also frequently termed the cingulo-opercular network (CON); (2) the frontoparietal network (FPN) (Corbetta et al., 2008); and (3) the DMN, which is usually deactivated during cognitive tasks (Buckner et al., 2008). The CON has been proposed to function as a core system for initiating and maintaining task sets as well as regulating and monitoring the activities of downstream systems (Dosenbach et al., 2006; Wen et al., 2013a). The FPN is recognized as being important for attention control (Corbetta et al., 2008; Scolari et al., 2015; Wen, Yao, et al., 2012a) and working memory (Harding et al., 2015; Wen et al., 2018). The DMN, although often deactivated in demanding cognitive control-related tasks, also is important for cognitive control, and its deactivation and functional integrity are associated with performance in goal-directed tasks (Buckner et al., 2008; Brewer et al., 2011; Wen et al., 2013a; Hearne et al., 2015). In the current study, we refer to these networks collectively as cognitive control-related networks.

Although the thalamus is thought to form a complex large-scale network with the foregoing cortical systems (Dosenbach et al., 2006), its function in this large-scale network needs further investigation. The role of the thalamus in corticocortical communication has rarely been tested in large-scale human brain networks related to cognitive control. However, in the intact human brain, understanding whether and how the thalamus communicates with the aforementioned networks, which include a large set of regions prominent in cognitive control, could be important for determining how brain functional networks communicate with the thalamus.

To address this question, we conducted a comprehensive examination to characterize the thalamus in terms of its communication with multiple networks related to cognitive control. For this study, we used fMRI, which is a noninvasive technique that records whole-brain activity and is a suitable approach for the calculation of thalamic-cortical connections. First, for resting-state fMRI datasets, we examined the intrinsic connectivity pattern of the thalamus in the cognitive control-related network from three angles: 1) using a novel multi-seed Granger causality (GC) mapping approach, we directly examined the directed information flow pattern of the thalamus regarding the cognitive control-related network; 2) by combining dynamic functional connectivity (DFC) calculation with principle component analysis (PCA), we investigated whether dynamically co-varying patterns exist among the functional connections (FC) between the thalamus and cortex (see Methods); 3) by employing time-varying community analysis based on DFC, we examined whether the thalamus has high functional flexibility in forming modules with other nodes, making it adaptive for communication in complex networks. The above analyses were conducted on two independent resting-state fMRI datasets. Second, for the task fMRI dataset, we examined the task effective connectivity pattern of the thalamus by repeating the multi-seed GC mapping on an attention task dataset, mapping the directed information flow pattern of the thalamus during the task. For convenience and accuracy of expression, we have used the term thalamus-cortex (T-C) to describe the communication/connections between the thalamus and the cortex to distinguish them from those of the direct cortical cortex pathway.

Method and materials

Overview of the analyses

To clearly describe the main ideas of the various analyses used in this study, we have used this section briefly outlines the steps and materials used in the analysis (see Fig. 1 for a schematic diagram), which can be divided into two parts, the resting-state analysis and the task analysis. The following sections provide detailed information in this regard.

Fig. 1
figure 1

Schematic summary of the analysis

Full size image

Part I included the analyses of resting-state fMRI data recorded from healthy participants at two independent institutes. Resting-state fMRI dataset A included data recorded from 197 participants at the Beijing Normal University. Resting-state fMRI dataset B included data recorded from 204 participants at Southwest University. We also used two region of interest (ROI) sets (ROI set X and ROI set Y) to test the robustness of the results when the central coordinates of the ROIs changed but remained within the range of cognitive control-related regions.

ROI set X included ROIs representing the task-positive network (TPN1) and the default mode network (DMN1) previously defined using independent component analysis (ICA) (Wen et al., 2012a, 2012b; Wen et al., 2013a). ROI set Y included the resultant ROIs of a meta-analysis for task-positive networks related to the term “cognitive control” (TPN3) and for the default mode network regions related to the term “default mode” (DMN2). Three analyses were performed by employing each of the four combinations of the two datasets and the two ROI sets.

  • Analysis 1: Multi-seed GC mapping, which falls within the scope of effective connectivity analysis for directed network modeling, was performed using the cognitive control-related ROIs set. It was designed to directly and intuitively examine whether the thalamus has intense information inflow and outflow with respect to cognitive control-related networks.

  • Analysis 2: DFC estimation combined with PCA was employed to examine whether the subcomponents of the thalamocortical DFC showed dynamically co-varying patterns, which, if present, could be an indicator of information propagation among a set of pathways (connections).

  • Analysis 3: Based on DFC, we applied community analysis, which examined the community structures containing intrinsic functional modules in the brain, to evaluate the extent to which a node region flexibly changes its affiliation to the community structures. This step was designed to examine whether the thalamic ROI has high functional flexibility, which would make it adaptive for communication in complex networks.

The analysis in Part II involved repeating the multi-seed GC mapping analysis on a task dataset (dataset C) to investigate whether the thalamus also has intense information inflow and outflow with respect to cognitive control-related networks during an attention task. This analysis was performed on two ROI sets for validation.

Details of the datasets, ROI sets, and analysis procedure are described in the following sections.

Participants and experiments

Three independent datasets were examined in this study. In Part I, the resting-state fMRI datasets A and B were examined using multi-seed GC mapping, DFC-PCA, and DFC-community analysis. In Part II, task fMRI dataset C was examined using multi-seed GC mapping.

Resting-state fMRI dataset A: Participants were recruited from Beijing Normal University (after one left-handed female was omitted, n = 197; 122 females, mage ± SD = 21.17 ± 1.83; 75 male, mage ± SD = 21.14 ± 1.83). All participants were right-handed and had no history of neurological or psychiatric disorders. They were required to close their eyes and remain still without thinking about anything in particular. Resting-state data was collected using a 3.0-Tesla scanner (Siemens Trio TIM); TR = 2,000 ms, TE = 30 ms, flip angle = 90 degrees, field of view = 200 × 200 mm2, matrix size = 64 × 64, voxel size = 3.12 × 3.12 × 3.6 mm3, slices = 33, time points = 225. All data were acquired from the open resource website https://www.nitrc.org/frs/shownotes.php?release_id=819.

Resting-state fMRI dataset B: Participants were recruited from Southwest University (n = 204; 125 females, mage ± SD = 19.96 ± 1.22; 132 male, mage ± SD = 20.30 ± 1.40). All participants were right-handed and had no history of neurological or psychiatric disorders. Resting-state data was collected using a 3.0-Tesla scanner (Siemens Trio TIM); TR = 2000 ms, TE = 30 ms, flip angle = 90 degrees, matrix size = 64 × 64, field of view  =  220 × 220 mm2, voxel size = 3.4 × 3.4 × 3.6 mm3, slices = 32, time points = 242. All data were acquired from the open resource website http://fcon_1000.projects.nitrc.org/indi/retro/southwestuni_qiu_index.html

Task fMRI dataset C: Participants were recruited from Beijing Normal University (n = 300; 168 females, mage ± SD = 21.04 ± 1.15; 132 male, mage ± SD = 21.21 ± 1.16). None of the participants had a history of neurological or psychiatric disorders, and 280 of them were right-handed. Each participant performed the visual-spatial attention task developed by Wojciulik and Kanwisher (1999), in which the participants were required to covertly direct and maintain attention to a cued peripheral dot or a central dot and detect the slight size change of the attended dot at a random moment during each 4-second period. The task lasted for 5.5 minutes, with 18 condition blocks each lasting for 16 seconds; 12 of them were attend-to-peripheral blocks (Wojciulik & Kanwisher, 1999). Task data were collected using a 3.0-Tesla scanner (Siemens Trio TIM); TR = 2,000 ms, TE = 30 ms, flip angle = 90 degrees, field of view = 200 × 200 mm2, matrix size = 64 × 64, voxel size = 3.12 × 3.12 × 4.8 mm3, slices = 30, time points = 162.

Data preprocessing

The first five time points collected from each participant or each run were discarded to remove unstable data. All data were preprocessed through slice timing, realignment (with a bounding box [−90 −126 -72; 90 90 108]), normalization to the Montreal Neurological Institute (MNI) template (Friston et al., 1995), and spatial smoothing with an 8 mm FWHM Gaussian core using SPM8 (http://www.fil.ion.ucl.ac.uk/spm/). Because previous studies suggested removing the nuisance component in the fMRI data (Beckmann et al., 2005; Calhoun et al., 2001), we used the GIFT v4.0a toolbox (http://mialab.mrn.org/software/gift/index.html) to perform ICA for each participant to remove spurious signals. Specifically, each participant’s data were decomposed into several independent components (ICs). The ICs that were mostly distributed along the edge of the brain, in the white matter, in the cerebrospinal fluid (CSF), or had less power within the frequency band (0.01–0.08 Hz) were considered noncognitive components, and the data were reconstructed based on the ICs excluding these (Kaushik et al., 2013; Mckeown & Sejnowski, 1998). Scrubbing (Power et al., 2012) also was used to remove the possible effect of severe head motion, and time points with a framewise-displacement (FD) greater than 0.5 mm were replaced by interpolation (Han et al., 2016; Satterthwaite et al., 2012). Global scaling was conducted to remove global noise. Twenty-four head-motion parameters (Friston et al., 1996), white matter signal, and CSF signal were regressed out to further eliminate head motion effects and other confounding factors. For resting-state data, detrending and band-pass filtering (0.01–0.08 Hz) were used to acquire a stable low-fluctuating blood oxygenation level dependent (BOLD) signal. For task data, detrending was used, but band-pass filtering was not applied to avoid removing fast components related to the task.

ROI selection

Cognitive control-related networks included the CON, FPN, and DMN networks. The CON and FPN are consistently activated while the DMN is deactivated during high-cognitive control tasks, such as attention. We defined multiple ROI sets to represent the foregoing networks using different strategies, which provided a generally consistent spatial pattern of the cognitive control-related network but slightly different central coordinates of the major nodes. This setting aimed to test the robustness of the results against slight changes in ROI coordinates.

Three TPN ROIs were defined (TPN1, TPN2, and TPN3).

  1. 1)

    TPN1 ROIs were a set of task-positive ROIs defined in previous studies based on a spatial attention task (Wen et al., 2013a; Wen et al., 2012a, 2012b), in which the TPN1 regions were determined using GLM analysis by contrasting the attention condition and the passive-view condition (p < 0.002, T > 5.20, FDR-corrected; Fig. 2). ROIs were defined as spheres with a 5-mm radius centered at the voxels with local T value maxima. Please refer to Wen et al. (2012a, b, 2013a) for more details.

  2. 2)

    TPN2 ROIs were defined based on the current dataset by following the analysis procedure of Wojciulik and Kanwisher (1999). The task-positive regions were determined using GLM analysis by contrasting the peripheral and central conditions at the individual level and using a one-sample t-test to yield group activation (p < 0.00001, T > 22, FWE-corrected; Fig. 2). ROIs were defined as spheres with a 5-mm radius centered at the voxels with local T value maxima.

  3. 3)

    TPN3 ROIs were defined using meta-analysis mapping (www.neurosynth.org; Yarkoni et al., 2011) of the term "cognitive control" (uniformity test: z > 3.60, p < 0.01, FDR-corrected; Fig. 1). ROIs were defined as spheres with a 5-mm radius centered at the voxels with local z value maxima. The TPN3 ROIs were also used to demonstrate the validity of the TPN1/2 ROIs as representative nodes for task-positive networks related to cognitive control (Fig. 2).

Fig. 2
figure 2

Task-positive network (TPN) regions and default mode network (DMN) regions related to cognitive control. A. Spatial overlap of the task-activated regions and results of the meta-analysis mapping regarding the term “cognitive control.” Red regions represent the previously defined TPN1 regions (Wen, Yao, et al., 2012a; Wen, Mo, & Ding, 2012b; Wen et al., 2013a); blue regions (TPN2) represent the group activation based on dataset C recorded in the current spatial attention task (n = 300, T > 22, p < 0.00001, FWE-corrected); green regions (TPN3) represent the results of the meta-analysis for the term “cognitive control” (598 studies, z > 3.60, p < 0.01, FDR-corrected); and white regions show the overlap of the three maps. B. Slice view of the local maximum activation points in the middle ROI displayed using BrianNavigator. C. DMN1 regions previously defined (Wen, Mo, & Ding, 2012b; Wen et al., 2013a). Brain surface view of the group deactivation (navy blue) by contrasting passive view blocks against attention blocks (T > 4.75, p < 0.005, FDR-corrected), the group resting-state IC map (T > 4.75, p < 0.005, FDR-corrected) referring to DMN (purple), and their overlap (green).

Full size image

The locations of thalamic ROIs were specifically tested by projecting Talairach space thalamic ROI centers using the Nonlinear Yale MNI to Talairach Converter (Lacadie, Fulbright, Rajeevan, Constable, & Papademetris, 2008) and superimposing the converted coordinates onto a detailed 3D brain atlas in the BrianNavigator software (http://www.thehumanbrain.info/brain/brain_navigator.php) (Fig. 1B).

The DMN ROIs were defined using two approaches (DMN1, DMN2).

  1. 1)

    The DMN1 ROIs (Buckner et al., 2008) are our previously defined DMN ROIs using ICA (Wen et al., 2013a; Wen, Mo, & Ding, 2012b). The DMN ROIs were defined as spheres with a 5-mm radius centered at the voxels with local T value maxima (p < 0.0002, FDR-corrected). These DMN ROIs also fell into the deactivation map on comparing the passive condition with the attention condition. Please refer to Wen, Mo, & Ding, 2012b; Wen et al., 2013a) for more details. The reason for not merely using the deactivation map to define the DMN regions is that the global scaling required by connectivity analysis might simultaneously exaggerate the deactivation results and falsely incorporate some deactivated regions that did not belong to the proposed DMN. ICA based on resting-state fMRI has been found to be stable in defining classical DMN regions (Buckner et al., 2008). Here, we adopted ICA to define the DMN ROIs, but used task deactivation for double-checking.

  2. 2)

    DMN2 ROIs were obtained by using meta-analysis mapping (www.neurosynth.org; Yarkoni et al., 2011) of the term "default mode" (association test: z > 3.60, p < 0.01, FDR). The DMN2 ROIs were also used to demonstrate the validity of the DMN1 ROIs as representative nodes for the DMN (Fig. 2).

Following the aforementioned procedure combining TPN and DMN, we defined 18–19 ROIs representing the major nodes of the CON, FPN, thalamus, and DMN. The ROIs included the dACC, right AI (rAI), and left AI (lAI) representing the CON; the bilateral frontal eye field (FEF) and intraparietal sulcus (IPS), right anterior and posterior medial frontal gyrus (raMFG and rpMFG), and right temporoparietal junction (rTPJ) representing the FPN; the dorsal and ventral medial prefrontal cortex (dMPFC and vMPFC), bilateral inferior parietal lobe (IPL), and posterior cingulate cortex (PCC) representing the DMN; and the bilateral middle thalamus. The dACC, bilateral AI, FEF, IPS, raMFG, rpMFG, dMPFC, vMPFC, and PCC were examined in our previous studies (Wen et al., 2013a; Wen, et al., 2012a, 2012b).

As mentioned in the overview, for cross-validation, we defined multiple ROI sets as follows:

  • ROI set X: TPN1 and DMN1 ROIs;

  • ROI set Y: TPN3 and DMN2 ROIs;

  • ROI set Z: TPN2 and DMN1 ROIs.

Note that because we did not obtain resting-state fMRI from the participants of task dataset C, we selected our previous resting-state ICA-based DMN ROIs (DMN1) to represent the DMN ROIs and combined them with TPN2 ROIs to form ROI set Z. Please refer to Fig. 1 for details of the combinations.

Measuring information flow using GC

GC analysis is an effective connectivity approach that is useful in directed network modeling. Let X and Y be two time series; the basic idea of GC is that if using both X’s and Y’s past, one can predict Y’s future better than merely using Y’s past, then we say that there is a Granger causal influence from X to Y (Granger, 1967). Specifically, this idea could be described using the multivariable autoregressive (MVAR) model (Ding et al., 2000; Chen, Bressler, & Ding, 2006; Wen, Rangarajan, & Ding, 2013c) to estimate GC.

Let the two stationary time series be denoted by Xt and Yt. Individually, Xt and Yt can each be represented by the following univariate autoregressive (AR) models:

$$ {\displaystyle \begin{array}{l}{X}_i=\sum \limits_{j=1}^{\infty }{a}_{1j}{X}_{t-j}+{\varepsilon}_{1j,}\operatorname{var}\left({\varepsilon}_{1t}\right)={\Sigma}_{\mathrm{l}}\\ {}{Y}_t=\sum \limits_{j=1}^{\infty }{d}_{1j}{Y}_{t-j}+{\eta}_{1t,}\operatorname{var}\left({\eta}_{1t}\right)={\Gamma}_{\mathrm{l}}\end{array}} $$
(1)

They can be jointly represented as the following bivariate AR model:

$$ {\displaystyle \begin{array}{c}\ {X}_t=\sum \limits_{j=1}^{\infty }{a}_{2j}{X}_{t-j}+\sum \limits_{j=1}^{\infty }{b}_{2j}{Y}_{t-j}+{\varepsilon}_{2t}\\ {}{Y}_t=\sum \limits_{j=1}^{\infty }{c}_{2j}{X}_{t-j}+\sum \limits_{j=1}^{\infty }{b}_{2j}{Y}_{t-j}+{\eta}_{2t}\end{array}} $$
(2)

If the bivariate AR model provides a better prediction of Y than the univariate AR model, then there is a Granger causal influence from X to Y, namely, an information flow from X to Y. More details on the calculation of GC are provided in our previous article (Wen, Rangarajan, & Ding, 2013b). GC measures directed information flow and has been widely used in directed network modeling based on multi-unit recordings, EEG, and fMRI, providing insightful results for complex interactions among large-scale networks.

For practical GC analysis based on fMRI data, we applied an additional preprocessing step, following a previous study, which provides an effective deconvolution method based on hemodynamic response function (HRF) blind detection (Wu Liao, Stramaglia, Ding, Chen, & Marinazzo, 2013); we used the HRF toolbox (http://users.ugent.be/∼dmarinaz/code.html) to deconvolve the preprocessed data to retrieve the signal change at the neural level from BOLD and alleviate the possible confounds caused by HRF inhomogeneity. GC estimation was performed based on the deconvolved fMRI time series. According to the BIC and AIC model order estimations and our previous studies (Harrison et al., 2003; Rogers, Katwal, Morgan, Asplund, & Gore, 2010; Wen, Mo, & Ding, 2012b), the order of the AR model was set to 2 in this study.

Mapping the inflow-outflow stations using multi-seed GC mapping

We introduced a multi-seed GC mapping analysis to measure the outflow and inflow degrees of each voxel when constructing a network with the ROI set which includes m seed ROIs representing the network of interest. That is, the causal influence each in-brain voxel exerted upon and received from the m ROIs was assessed by calculating the GC between the time series of the voxel and that of the m ROIs. The in-/out-degree of the current voxel was calculated by averaging the m in-/out-GC values. Repeating this procedure for each voxel, we generated the in- and out-degree maps for the network represented by the m ROIs.

Compared with FC based on correlation analysis, GC estimation is more vulnerable to noise, in addition to the lower signal-to-noise ratio (SNR) found in the thalamus (Demetriou et al., 2018), the thalamic results might be obscured by noise in the individual raw GC map. Therefore, we spatially smoothed the individual maps using an 8 mm FWHM Gaussian core to improve the SNR of the maps and then z-scored them. A one-sample t-test was then performed to generate the group in- and out-degree maps using SPM8. The overlap between the in- and out-degree maps was considered to highlight the region having both high in- and out-degree, which indicates the important role of receiving and sending information regarding the particular network of interest.

This multi-seed GC mapping analysis was applied to the combinations of resting-state and task datasets and ROI sets, as illustrated in Fig. 1.

Extraction of the dynamically co-varying pattern of dynamic functional connectivity based on principal component analysis

The DFC-PCA aimed to investigate the inherent dynamic patterns that exist among multiple connections. Its basic idea is to consider the functional connectivity between any of the two given regions as a time-varying variable. At each moment, the instantaneous functional connectivity is considered as a linear combination of multiple inherent components. Each component can be represented as a vector in the N-dimensional space with each dimension corresponding to an ROI pair (a connection), where N denotes the number of connections. Each component corresponds to a time-varying weight value (or contribution score). Some components may contribute to the overall FC prominently, which have been vividly and repeatedly described as (connectivity) patterns in the literature. These components are the so-called principal components (PCs). If a PC has large absolute values (projections) on a certain set of connections (subspace), we say that these connections comprise components that fluctuate simultaneously.

Our DFC-PCA was derived from the standard pipeline proposed by Leonardi et al. (2013), except that a bootstrapping procedure was added (10,000 iterations) to overcome sampling errors and obtain stable PCs. The analysis was performed using the following steps.

  • Step 1: For each participant, a sliding window with a length of 25 TRs (50 seconds) was set to move along the time axis with a step of 1 TR. As a previous study (Preti et al., 2017) suggested that a window within 30–60 seconds is optimal for capturing DFC characteristics, we set the window length to 20 TRs (40 seconds) to validate the robustness of the DFC results (please refer to the Supplementary material). In each sliding window (each time point), the correlation between the fMRI time series of any of the two ROIs was calculated to yield an m × m FC matrix for each time point.

  • Step 2: The FC matrices were then vectorized and stacked together, their absolute values were Fisher-transformed and z-scored to remove the influence of individual differences, and the temporal mean of each ROI pair was removed. The column number N of the individual stacked FC matrix was the number of the connections (the N ROI pairs among an ROI set containing m ROIs, N = m × [m – 1] / 2) and the row number of the stacked FC was the number of temporal windows.

  • Step 3: In each iteration, we randomly selected 100 participants (about half of the total participants) and concatenated their individual stacked FC matrices along the time axis, with each row indicating an observation of connectivity strength in each sliding window and each column a group time series of the DFC of an ROI pair (connection). For convenience, we termed this concatenated matrix the DFC time series. The DFC time series was then fed into the PCA to generate N PC vectors, each of which comprised N values representing an underlying temporal co-varying pattern of the DFC of the m ROIs in the resting-state.

The PCA regarded the concatenated matrices as a weighted sum of several PCs (or connectivity states).

$$ FC(t)=W1(t)\times PC1+W2(t)\times PC2+W3(t)\times PC3+ $$
(3)

Here, FC(t) is a time-varying connectivity vector, and PCi are PCs that are regarded as inherent connectivity patterns (or states) that contribute to the overall FC. The contribution of PCi is described using Wi(t), which varies from time to time. The higher the variance of Wi(t), the higher the contribution of PCi to the variance of FC(t), and PCi is ranked in a more frontal position in the PC arrays. For each PC rank, the results were averaged across all iterations to yield a stable mean PC for that rank.

Step 4: Based on previous research, we used Kendall's coefficient of concordance (w) to measure the consistency of the PCA results across the iterations in each PC rank (Baumgartner et al., 1999; Liu et al., 2010). The mean PC of the ranks with consistency higher than 0.5 were reserved as prominent dynamically co-varying patterns, while other PCs were considered too noisy and were excluded from the analysis (Gisev et al., 2010; Levitan et al., 2008).

Step 5: The foregoing prominent PC vectors were pooled together, and their standard deviations σ and μ were calculated. The values of μ + 1.96 σ and μ - 1.96 σ were used as the upper and lower thresholds, respectively, to highlight the significant dynamically co-varying pattern of the FC. Each thresholded PC vector was reconstructed into an m × m DFC pattern matrix, which denoted an inherent, principal dynamically co-varying pattern of the connectivity along with the time-varying contribution to the overall DFC. In other words, the reconstructed matrix depicted some basic connectivity patterns that were temporally stable and contributed to the overall DFC with prominent weights.

Examination of node flexibility using DFC-based community analysis

In recent years, community analysis has been used as an influential tool for detecting functional modules of human brain networks (Allen et al., 2011; Bassett et al., 2013; Chang & Glover, 2010; Mucha et al., 2010). Groups of nodes that are more tightly connected may assemble into structures called communities (Bassett et al., 2013). These community structures contain functional modules in the brain, and changes in community structure could embody meaningful dynamic patterns of functional connectivity that relate to changes in cognitive function. Based on the DFC matrix, we computed the flexibility to quantify the frequency at which the allegiance of each ROI to whole network communities changed over time. Specifically, for the sliding window method with n windows, the maximum time that a brain region could be affiliated with community changes is n - 1, and m is the actual number of times that the affiliated community of a brain region changes. Thus, the flexibility of an ROI is the ratio of m to n - 1. This is an important measurement of the dynamics and activation of brain regions (Bassett et al., 2013). A high flexibility value indicates that the community affiliation of that particular region changes often, indicating that the ROI dynamically cooperates with multiple modules/systems, which could be seen as one of the necessary properties of a node that is important for network communication.

It is worth noting that PCA-based DFC analysis has shown that minor PCs might contribute to connectivity variability and are more vulnerable to noise. Therefore, we used the aforementioned PCA before community analysis to extract the major components of DFC and calculated the flexibility based on the major PC-reconstructed DFC. Specifically, PCs that contained the top 80% of the variation were reserved for community analysis to determine the flexibility of each ROI in the cognitive control-related networks.

Results

ROIs related to attention and cognitive control

Both TPN1 and TPN2 highlighted the right and left middle thalamic ROIs, which were co-activated with the CON and FPN. The regions centered by the activated ROIs were highly consistent with the CON and FPN regions outlined in the meta-analysis mapping of the term "cognitive control" (Fig. 1). They are also consistent with many classical findings that report co-activation of the CON and FPN along with simultaneous deactivation of the DMN during many goal-directed tasks (Brewer et al., 2011; Dosenbach et al., 2006; Wen et al., 2018), as well as those reported to have coupled spontaneous BOLD activity during the resting-state (Buckner et al., 2008; Fox et al., 2005; Fox et al., 2006; Fransson, 2010). Therefore, we validated the selection of the 18~19 ROIs as the major nodes representing the cognitive control-related networks (Tables 1, 2, 3, 4 and 5).

Table 1 Selected task-positive ROIs (TPN1) previously defined using the group activation of a visual-spatial attention task (Wen, Yao, et al., 2012a; Wen, Mo, & Ding, 2012b; Wen et al., 2013a)
Full size table
Table 2 Selected task-positive ROIs (TPN2) based on group activation of the current spatial attention task
Full size table
Table 3 Selected task-positive ROIs (TPN3) based on the meta-analysis (“cognitive control,” 598 studies)
Full size table
Table 4 Selected DMN ROIs (DMN1) previously defined using ICA (Wen, Mo, & Ding, 2012b; Wen et al., 2013a)
Full size table
Table 5 Selected DMN ROIs (DMN2) based on the meta-analysis (“default mode,” 777 studies)
Full size table

Part I: Regions showing inherent intense causal interactions with cognitive control-related networks

According to the results of the multi-seed GC mapping analysis (Fig. 3), several regions showed high in- or/and out-degree when associated with cognitive control-related networks (Dataset A + ROI set X/Y: p < 0.0001 FDR-corrected; Dataset B + ROI set X/Y: p < 0.00001, FDR-corrected). Specifically, the CON regions, for example, the dACC and AI, and posterior part of the PCC showed high out-degree, the anterior MFG and anterior PCC showed high in-degree, and the rAI and right middle thalamus showed both high in- and out-degree (overlapped regions in yellow, Fig. 4), indicating the functional role of both receiving and sending information in the cognitive control-related networks.

Fig. 3
figure 3

Multi-seed Granger causality mapping demonstrates consistent overlap (yellow) of both high-input (red) and high-output (green) regions during resting in middle thalamus across combinations of different ROI sets and resting-state datasets. A. High in- and out-degree maps regarding ROI set X and Dataset A (n = 197, p < 0.0001, FDR-corrected). B. High in- and out-degree maps regarding ROI set X and Dataset B (n = 204, p < 0.00001, FDR-corrected). C. High in- and out-degree maps regarding ROI set Y and Dataset A (n = 197, p < 0.0001, FDR-corrected). D. High in- and out-degree maps regarding ROI set Y and Dataset B (n = 204, p < 0.00001, FDR-corrected). The yellow area denotes the regions showing both high in- and out-degree Table 5.

Full size image
Fig. 4
figure 4

PCA results indicated the dynamically co-varying pattern among the dynamic functional connections. Panels A, B, C, and D show the first six principal components (PCs), reflecting the dynamically co-varying patterns of the dynamic functional connections in the two independent datasets and using ROI set X and ROI set Y. For each PC, the color of the cells indicates the extent to which the connections are coupled with each other. The more intense red connections were more strongly coupled with each other, while the more intense blue connections were more strongly coupled with each other and had dynamics inverse to the red connections. The green connections were decoupled with either of the red or blue connections. A. connectivity coupling patterns (PCs) based on dataset A and ROI set X. B. connectivity coupling patterns (PCs) based on dataset B and ROI set X. C. connectivity coupling patterns (PCs) based on dataset A and ROI set Y. D. connectivity coupling patterns (PCs) based on dataset B and ROI set Y. E. An example of the dynamically co-varying functional connections between the thalamus and other nodes in PC3 of panel B (the green nodes in panel E are the bilateral middle thalamus).

Full size image

We also compared the multi-seed mapping results with the thalamus and those without the thalamus to show the percentage drop in in- and out-degree, which reflected the contribution of thalamic seeds. The result (Figure S3) showed similar trends for both resting-state datasets: when removing the thalamus from the model, the drop in the in-degree was approximately 10–12% for both the thalamus (defined by the AAL template in MRIcron) and cortical regions significantly receiving input (red regions in Fig. 3, excluding the thalamic parts). On the other hand, the drop in the out-degree was approximately 10% for thalamic regions, and approximately 5% for cortical regions significantly exerting output (green regions in Fig. 3, excluding the thalamic parts). More details regarding the comparison are provided in the Supplementary material.

Part I: the intrinsic dynamically co-varying patterns of DFC

Using DFC-PCA, we acquired the PCs that indicate the dynamically co-varying patterns of the networks related to cognitive control and ordered them according to their contribution to the overall DFC. The first six prominent PCs were preserved as representing stable patterns according to their Kendall's coefficients in the bootstrap test (see Methods for more details). Consistent across all four combinations between the resting-state datasets A/B and the ROI sets X/Y, the results showed that PC1 reflected a dynamically co-varying pattern of the direct corticocortical connections, while some of the other PCs reflected dynamically co-varying patterns in almost all corticothalamic connections (please see the “hot” cross shape in PC3 of dataset A/B + ROI set X, dataset B + ROI set Y, and PC4 of dataset A/B + ROI set Y). In addition, PC2 showed dynamically co-varying patterns in FPN-related connections, including some corticothalamic and CON-related connections. The patterns depict the co-variation of the dynamic components of the connections between the thalamus and cortex. The results intuitively demonstrated that communication between the thalamus and cortex and direct corticocortical communication coexist among the cortical sites related to cognitive control. These results were preserved when setting the window length to 20 TRs (40 seconds) for DFC estimation (Supplementary material, Figure S1).

Part I: Flexibility of cognitive control-related ROIs

The results of community analysis for all four combinations between the resting-state dataset A/B and the ROI set X/Y consistently showed that the left and right thalamus, right temporoparietal junction (rTPJ), right middle frontal gyrus (rMFG), and left frontal eye field (lFEF) were the most prominent nodes with high flexibility (Fig. 5). Among these highly flexible nodes, the rTPJ and rMFG have been proposed as important cortical stations for the interaction of attention control signals (Corbetta et al., 2008; Fox et al., 2006). All CON and DMN nodes showed low flexibility for all four combinations. The thalamic flexibility patterns were preserved when setting the window length to 20 TRs (40 seconds) for DFC estimation (Supplementary material, Figure S2).

Fig. 5
figure 5

Z-scored flexibility of bilateral middle thalamus (rmTha/lmTha in orange) and the 17 other ROIs based on community analysis. A. The flexibility of the nodes of ROI set X for dataset A and dataset B, respectively; the flexibility of mid-thalamic ROIs is significantly higher than those of raMFG, rFEF, all CON ROIs, and all DMN ROIs (p < 0.01, FDR-corrected). B. The flexibility of the nodes of ROI set Y for dataset A and dataset B, respectively; the flexibility of mid-thalamic ROIs is significantly higher than those of most FPN ROIs (except rTPJ), all CON ROIs, and all DMN ROIs (p < 0.01, FDR-corrected). C. Averaged flexibility at the network level for all combinations between dataset A/B and ROI set X/Y; the mean flexibility of the thalamus is significantly higher than the mean flexibility of the ROIs in all other networks (p < 0.001, FDR-corrected).

Full size image

Part II: Regions showing intense causal interactions with cognitive control-related networks during tasks

The results of the multi-seed GC mapping analysis using task dataset C (Fig. 6) highlighted several regions that showed high in- or/and out-degree when associated with the cognitive control-related networks represented by either ROI set X or ROI set Z (p < 0.0001 FDR-corrected; p < 0.0001, FDR-corrected). Specifically, the right middle thalamus, together with dACC, rFEF, and rpMFG, which have been proposed as important prefrontal nodes for top-down control, showed both high in- and out-degree (overlapping regions in yellow, Fig. 6). Other cortical regions, such as the PCC, vMPFC, bilateral IPL, rTPJ, raMFG, and bilateral lateral occipital cortex, mainly showed high in-degree.

Fig. 6
figure 6

Multi-seed Granger causality mapping base on the task dataset C demonstrates consistent overlap (yellow) of both high-input (red) and high-output (green) regions in the right middle thalamus during the visual-spatial attention task. A. High in- and out-degree maps (n = 300, p < 0.0001, FDR-corrected,) regarding ROI set X. B. High in- and out-degree maps (n = 300, p < 0.0001, FDR-corrected) regarding ROI set Z

Full size image

Finally, we compared the task-state multi-seed mapping results with the thalamus and those without the thalamus to show the in- and out-degree percentage drops. The result (Figure S3) showed similar trends for dataset C: when removing the thalamus from the model, the drop in the in-degree was approximately 11% for both the thalamus (defined by the AAL template in MRIcron) and cortical regions significantly receiving input (red regions in Fig. 4, excluding the thalamic parts). On the other hand, the drop in the out-degree was approximately 10% for thalamic regions, and approximately 8% for cortical regions significantly exerting output (green regions in Fig. 4, excluding the thalamic parts). Further details regarding the comparison are provided in the Supplementary material.

Discussion

Using several types of analysis, this study tested the functional role of multiple regions that form a set of complex networks, both task-positive and task-negative, related to cognitive control. The results showed that thalamic ROIs tended to exhibit important characteristics in every aspect examined. First, for both resting-state and the attention task, the multi-seed GC mapping showed that the right middle thalamic region had both high in- and out-degrees for the directed network formed by the ROIs, indicating its function of both sending information to and receiving information from cognitive control-related networks. Second, in the resting-state, parallel to the direct corticocortical connections, all the dynamic T-C functional connections showed a dynamically co-varying pattern that predominantly contributed to the overall functional connectivity. Second, when considering the flexibility of the nodes that are dynamically organized into communities, the thalamic ROIs were included in the list of those that had high flexibility, allowing for adaptability in communication within the complex network.

Information flow between the thalamus and cortex

The challenge in identifying the role of nodes in a non-directed complex network is that the FC measurement of a connection does not describe how the information propagates from one node to another. Classic FC parameters that describe the importance of a node, such as node degrees, can highlight important nodes called network hubs. However, important nodes with a high degree may contribute differently to the network. That is, nodes that are prominent in exerting influences (sources) on others, those that are prominent in receiving influence (sinks) from others, and those that are prominent in both receiving the influence of other nodes and exerting an influence on other nodes, all show a high degree in classical functional connectivity analysis (non-directed view), although their functions are very different in a directed network. Using only non-directed network models is inadequate to provide a clear functional distinction between such stations and other important nodes. It is worth noting that in graph theory, being prominent in having both in-connections and out-connections in a directed network may be associated with the function of “relay,” but this could be inaccurate in practice. For example, if we observe information flow from X→Y and from Y→Z, we cannot conclude that Y is relaying information from X→Z, because X→Y and Y→Z could be independent. For a convenient and accurate expression, we call the nodes that show both high in-degree and high out-degree “complex” nodes, which are in contrast to the simple “source” nodes and “sink” nodes.

The resting-state multi-seed GC mapping results are useful for distinguishing functional regions that act as sources, sinks, or complex stations. For both resting-state datasets, the main regions of the dACC and bilateral AI significantly influenced other regions, while the raMFG primarily receives influences, which is consistent with their proposed roles in cognitive control (Fig. 4). Currently, the dACC and bilateral AI are thought to comprise the core network for top-down control (Dosenbach et al., 2006; Menon & Uddin, 2010; Wen et al., 2013b; Ham et al., 2013; Parro et al., 2018), whereas the raMFG is considered a part of the ventral attention network that receives top-down control signals and bottom-up sensory signals (Corbetta et al., 2008; Farrant & Uddin, 2015; Fox et al., 2006; Wen, Yao, et al., 2012a). However, in the current analysis, the right middle thalamic regions showed both high in- and out-degree, indicating a specific role in both receiving and sending information among the cortical sites.

It is worth noting that the rAI, though majorly highlighted as having a high out-degree (Fig. 4, green area in the rAI), also shows overlapping regions with high in-degree. This result is consistent with previous findings showing that the rAI is an important cortical site for mediating control signals to other systems (Cai et al., 2016; Hilger et al., 2017; Menon & Uddin, 2010). On the other hand, many previous studies have suggested that the PCC is an important hub of the DMN (Buckner et al., 2008). The current results show that the anterior part of the PCC mainly receives influence, the posterior part of the PCC mainly exerts influence, and that these two regions had very little overlap. Taken together, the multi-seed GC mapping results of known cortical regions, such as the dACC, AI, and PCC, were consistent with prior knowledge of their functional roles, indicating the efficacy of this method in detecting and distinguishing relay stations from other important nodes. Therefore, the cortical results can be regarded as a methodological validation of the thalamic results.

The task-state multi-seed GC mapping results confirmed the in-and-out function of the right middle thalamus as a complex node. On the other hand, the spatial pattern of the regions of high in-degree and high out-degree is different from that during the resting-state. The DMN regions become prominent only when receiving information. This is consistent with the idea that the DMN may be a source of task-irrelevant internal interference (see also the result shown in Fig. 4), which may impair the performance of the ongoing external task and therefore need to be suppressed during it (Buckner et al., 2008; Wen et al., 2013a). The cortical regions showing high out-degree included the dACC, rFEF, and rpMFG, demonstrating a spatial pattern in line with the previously proposed PFC regions, which are thought to exert top-down control signals (Corbetta et al., 2008; Dosenbach et al., 2006; Wen, Yao, et al., 2012a; Wen et al., 2013a). The results of the cortical sites that were consistent with the previous proposal provided support to the thalamic results, which indicate the important role of receiving and sending information in a large-scale network related to cognitive control.

Additionally, the systematic comparison of multi-seed mapping results with the thalamus and those without the thalamus was helpful in understanding the contribution of the thalamus. Given that in multi-seed mapping, the degree of each seed is equivalent to summing the connectivity strength between the voxel and each seed ROI, the 10–12% drop in the input/output degree in thalamic regions (Figure S3) indicates that removing the specific thalamic seeds from the multi-seed set only causes a limited degree drop, which means that the connection between the thalamic voxels and the highlighted cortical regions contributed significantly to the degree of the thalamic regions. The 6–12% drop in the input/output degree within highlighted cortical regions (Figure S3) reminds us that the thalamus is not the only region that has strong effective connectivity with the multi-seed ROIs, which is consistent with the finding that overlapping regions (Figs. 3 and 6, marked in yellow) were found in multiple cortical sites.

Contribution to the network interaction model of attention

The task results also provided a novel perspective on how large-scale networks interact under attention-demanding conditions. For large-scale networks activated during attention, researchers have proposed certain models that interpret the interaction between well-defined functional networks. For example, Corbetta et al. (2008) proposed an interactive model between the dorsal attention network and the ventral attention network to explain how control signals propagate between the two frontoparietal cortical systems interacting during attention. This model is being accepted by an increasing number of researchers (Wen et al., 2012b; Vossel et al., 2014; Shulman et al., 2017; Ebrahimpour et al., 2019). We have attempted to expand the exploration to the network interaction between the attention-positive and attention-negative networks, such as between the CON and DMN, during attention (Wen et al., 2013a). However, the foregoing model did not consider the role of the thalamus in communicating with attention-related networks. On the other hand, several important animal studies aimed at investigating specific thalamic nuclei have suggested that the thalamus might relay corticocortical communication or even play a more complex role than that of a simple relay deploying attention (Halassa & Sherman, 2019). Nelson, Dosenbach, Cohen, Wheeler, Schlaggar, & Petersen, (2010) proposed a multi-network task-control architecture, mainly based on functional connectivity analysis during the resting-state. This architecture speculates that the thalamus functions as a component of the cingulo-opercular system based on functional connectivity strength, but does not delineate its role in communicating with multiple networks. From the perspective of the large-scale network of the human brain and using GC as a tool for measuring information flow, our task results intuitively demonstrated that the right middle thalamus has a simultaneous input and output pattern similar to that of the dACC, FEF, and rpMFG, which are already known to be important for attention-related control. Specifically, the dACC is considered the core node of the cingulo-opercular system, and the right FEF is considered the core node of the frontoparietal system. Both systems are considered essential for attention control. This pattern distinguishes the middle thalamus, dACC, and FEF from other ROIs related to cognitive control. It indicates that communication between the middle thalamus and the network is also important for attention control.

The dynamically co-varying pattern of functional connections between the thalamus and cortex

The PCA-based DFC analysis was used to detect the underlying patterns of functional connectivity that repeatedly appeared across time and participants. These patterns were identified as connectivity matrices in the PC space, indicating the coupling strength between the connections. Given that each connection was calculated according to the Pearson correlation between the BOLD activities in a short moving temporal window, the index of each element in the matrix can be thought of as “the coupling of the coupling strength.” In other words, the strengths of connections A and B are two variables that vary with time, and the extent to which they co-vary could be important in estimating the building blocks of dynamic connectivity (Leonardi et al., 2013). As discussed above, this could be important to assess whether some pathways comprise dependent components.

The three or four most prominent PCs embodied patterns that were consistent across the two resting-state datasets. The first PC reflected the general coupling between most corticocortical connections. The second, third, or/and fourth PCs demonstrated that the T-C connections (rows 11–12 and columns 11–12 in the connectivity matrix) showed strong covariant patterns (Fig. 2). Specifically, for the second PC, the T-C connections also co-varied with many trans-FPN connections (including CON-FPN, FPN-FPN, and FPN-DMN connections), which have been proposed to be essential in cognitive control (Fox et al., 2005; Stern et al., 2012; Vincent et al., 2008). Thus, the PCA results suggested that the flow of information between cognitive control-related networks might include three major and stable components. The first component reflects direct communication between cortical ROIs. The second component reflects communication through trans-FPN and T-C connections. The third component corresponds to communication that mainly occurs through T-C connections.

Flexibility of the thalamic ROIs in the cognitive control-related network

The PCA-based DFC analysis results provided insight into the components of the T-C connections that couple with different systems. The flexibility calculation, which was based on community analysis, examined the capability of nodes to switch membership among different communities. Nodes with high flexibility are more likely to assume responsibility for relaying to multiple networks.

The results of the current study indicate that thalamic ROIs can flexibly mediate different signals associated with different cognitive functions. The results of the flexibility analysis were consistent with the results of the PCA. Specifically, the community results showed that the thalamic ROIs were highly flexible, which allows them to be well adapted for communication in complex networks. Other regions showed a high degree of flexibility, chiefly including a group of FPN regions that are thought to play important roles in a variety of cognitive functions, such as selective attention, attention reorientation, and working memory (Corbetta et al., 2008; Dosenbach et al., 2008; Harding et al., 2015). For example, previous resting-state and task fMRI studies have suggested that the ventral frontoparietal nodes, including the rTPJ and raMFG, are stations where bottom-up and top-down control signals interact (Corbetta et al., 2008; Fox et al., 2006; Wen, Mo, & Ding, 2012b). As mentioned above, in the results of the PCA-DFC analysis, we observed strong coupling between the FPN connections, trans-FPN connections, and T-C connections. The current PCA-DFC and the flexibility analysis not only confirmed a set of previously defined cortical relays in the FPN but also suggested that thalamic ROIs may also play an important role in control-related network communication.

Issues in DFC analysis

Theoretically, the use of DFC analysis for estimating the functional connectivity within a short sliding window may suffer from an issue similar to that encountered in time-frequency decomposition. One can increase the window length to obtain a better estimation of FC in a single-window, but at the cost of sacrificing the temporal details of FC changes. One can shorten the window to capture more temporal details, but lose the precision of FC estimation, making the analysis more vulnerable to sampling errors or interference due to factors such as head motion, scanner noise, and physiological noise. This motivated us to adopt an improved DFC analysis strategy involving bootstrapping and stable PC detection, and to perform the analysis on multiple datasets with a large sample size.

The issue of artifacts (e.g., head motion) contaminating DFC analyses has been the subject of much recent debate (Laumann et al., 2017). In addition to using multiple approaches to alleviate the impact of head motion, we also tested whether head motion could affect our results. We extracted the DFC time series data of all ROI pairs (including the thalamocortical pairs) and correlated them with each head-motion parameter; this test was performed for each participant. Only a small number of significant correlations (p < 0.05, uncorrected) were found in the random ROI pairs of random participants, which were not significant after correction for multiple testing. Namely, we found no systematic correlation patterns between our DFC time series and head motion.

GC

The GC method directly calculates the directed influence between time series, intuitively showing how a node exerts influence upon and is influenced by other nodes. Unlike the static correlation method, which only considers the point-to-point relationship between time series and ignores temporal mutual information, GC relies on the temporal dynamics of time series (Wen, Mo, & Ding, 2012b). It has been successfully used to explore the directed influence between neuroelectronic signals and has provided many insightful results in fMRI studies (Shirer et al., 2012; Wen, Yao, et al., 2012a; Wen et al., 2013a; Anderson et al., 2015; Solo, 2016; Wang et al., 2016). However, when applied to fMRI, the results of GC analysis may be confounded because of the inhomogeneity of the HRF (Smith et al., 2011). Although there is an ongoing debate on whether GC is a viable method for analyzing fMRI time series (Smith et al., 2011; Friston et al., 2013; Wen et al., 2013b; Seth et al., 2015; Barnett et al., 2017), many researchers have reached the consensus that 1) GC at the fMRI level and underlying neural causal influences are correlated; 2) caution needs to be used in the interpretation of results when considering noise and HRF inhomogeneity; and 3) HRF inhomogeneity detection and deconvolution may help increase the credibility of fMRI-GC. In this study, a deconvolution method based on HRF-blind detection (Wu et al., 2012) was used to avoid the possible confounding effect of HRF inhomogeneity and render the results more interpretable at the neural level.

The calculation of GC followed a novel strategy known as multi-seed mapping. Traditional seed mapping methods consider the relationship between whole-brain voxels and only one seed region at a time. The multi-seed mapping method examines the relationship between a voxel and multiple seed ROIs. Strictly speaking, it examines the role of a voxel after it joins a predefined network of interest with multiple seeds. This role can be measured using parameters such as degree, centrality, etc. The traditional seed mapping method can be considered as a special type of multi-seed mapping of weighted node degree with only one seed and using correlation as the metric. The novel method expands the concept of seed mapping, in which the “seed” becomes a network instead of a single region. Multi-seed mapping can be regarded as a combination of the model- and data-driven strategies, especially suitable for searching for nodes that are important to a specific network of interest, which in the current case is the cognitive control-related network.

In the fMRI scenario, unlike the traditional functional connectivity strength measured using cross-correlation, which measures the similarity between the major fluctuations of BOLD activities, GC considers the subtle temporal dependence, which is only a small portion of the total interdependence between the BOLD time series recorded by fMRI (Wen et al., 2012b, Wen et al., 2013a). The GC method is therefore more vulnerable to noise (Wen et al., 2013a). This vulnerability is reflected in the mapping result, which is noisier than the CC mapping results, and could impair subsequent statistical testing.

Theoretically, spatial smoothing might result in a loss of spatial resolution, but it actually has advantages that outweigh the drawbacks. When conducting the z-normalization on the GC values, we subtracted the mean (m) of the map from the raw value and divided the remainder by the map's standard deviation(s). However, when feeding into a statistical test, such as the t-test, the standard deviation(s) becomes larger if the map is noisier. Thus, the chance of detecting z-scored values significantly deviant from 0 would be smaller, inducing a higher false-negative rate. Lowering the statistical threshold to highlight the "strong" values actually impairs the spatial resolution of the results by introducing more false-positive results. On the contrary, using spatial smoothing on the raw GC map will significantly reduce the standard deviation components caused by noise and better reflect the spatial pattern.

Interrelationship between multiple indicators and method considerations

In summary, four indicators were used to characterize the role of the nodes in the cognitive control-related network, which includes thalamic and cortical ROIs (Table 6). Multi-seed GC mapping directly examined regions that intuitively show high in- and out-degree during the resting-state and task, respectively. The PCA-based DFC analysis was used to describe the dynamically co-varying pattern (PC) of the dynamic FC. The community analysis also examined the flexibility of the node ROIs from the perspective of the DFC.

Table 6 Abbreviated summary of the indicators and major results
Full size table

When using GC as the metric of the edges between the nodes, by constructing a directed network, the regions that show both high in- and out-degree can be regarded as prominent in both sending information to and receiving information from the remaining part of the predefined network.

It is worth noting that, although having both high in- and out-degree may appear to be necessary for relay, it is insufficient to define the relay station. For example, X receives a message from Y1 and sends a message to Y2, but the two messages are not the same. This is a very important point that motivated the DFC-based analysis. If a set of connections (or node pairs) share co-varying components and have common nodes connected with all the other nodes in the network, for example, the FPN and thalamic nodes for PC2 and the thalamic nodes for PC3, it is highly possible that these “common nodes” play a relay role in the communication between all nodes.

On the other hand, community analysis characterizes the flexibility of the nodes, providing information about how a node is dynamically connected with the sub-groups in the network. For example, let N1 and N2 be two sub-networks and R be a node important for the communication of both N1 and N2; one may observe that R is connected to N1 at one moment and N2 at another, exhibiting high flexibility. However, it is difficult to exclude the possibility that a node is stably connected to certain members of the network and shows low flexibility. Thus, one can still argue that highly flexible nodes may include some sources/sinks that are dynamically connected to different sub-groups.

Both methods, based on the examination of dynamic characteristics of the functional connectivity, bear the same limitation in that functional connectivity provides no direction of information transmission, which is critical for delineating the complex stations in a network.

In summary, we explored the role of the thalamus from three perspectives: 1) multi-seed GC mapping, 2) DFC-PCA, and 3) DFC-flexibility testing. As discussed above, analyses 1) and 2) together may provide some insight regarding the relay function of the thalamus. However, analysis 3) by itself did not aim to address the exact relay function. The flexibility analysis served as a test to demonstrate the dynamic characteristics of the thalamus in adapting to communication in complex networks, which could be an important property. Our analysis aimed to provide a comprehensive picture of the characteristics of thalamic-cortical communication among cognitive control-related networks from the perspective of large-scale network analysis based on human fMRI, rather than concluding the thalamus’s function as a simple relay, which is insufficient according to the recently proposed broad and complex functions of the thalamus (Halassa & Sherman, 2019).

As mentioned above, using a single indicator may be inadequate for assessing the significance of the nodes in the network. Collectively, the middle thalamic regions were highlighted by all four indicators, strongly suggesting an essential role in mediating the signals from cortical sites to thalamic regions and vice versa. The results also showed that some cortical ROIs may play essential roles in the network. Some of these, such as the rAI, rTPJ, and rpMFG, have been proposed to be important hubs of the cognitive control-related network (Corbetta et al., 2008; Menon & Uddin, 2010).

Compared with previous single-method and single-variable studies, this study benefits from the perspective of multiple measurements. In addition, to illustrate the reproducibility of the results, this study verified the results in different fMRI datasets (both resting-state and task).

Limitations

GC analysis is regarded as a method for depicting information flow, and specifically, whether the past of one time series can facilitate the prediction of the future of another time series. Its temporal information-based nature as well as its practical realization using the MVAR model makes GC analysis fall within the scope of multivariate analysis, which is different from the traditional co-activation analysis. Although GC has been considered as a convenient tool for measuring the information flow between brain regions, we need to realize that a further encoding-decoding model for task-relevant variables would be necessary to better address whether and how the thalamus propagates information regarding task-relevant variables, which is beyond the capability of the current design and should be considered in future studies.

Although this study focused on examining the interactions between the thalamus and cortex among the cognitive control-related networks, we need to be very careful when drawing the exclusive conclusion that the middle thalamic ROIs were specifically geared toward coordinating cognitive control-related regions. This study was not designed to exclude the possibility that the middle thalamic ROIs may play a role in sending or receiving information to or from other cortical systems. Answering whether the role is specific to the cognitive control-related networks or general for all kinds of networks may also require further comprehensive investigation across different systems in the future.

It is worth noting that the current study did not directly examine T-C connectivity and behavioral performance; therefore, it is lacking in terms of behavioral specificity for further interpreting the function of the thalamus. To fully understand the role of the thalamus in cortical network interactions and the behavioral significance of these interactions, more sophisticated designs and methods that enable direct investigation of the relationship between task GC and behavioral performance during cognitive control are required.

Conclusions

The current study conducted a comprehensive exploration of the role of the thalamus in the communication between multiple cognitive control-related networks. From the perspective of effective connectivity, the thalamic regions exhibit both high in- and out-degree, functioning as nodes of a large-scale complex network that comprises multiple intrinsic networks related to cognitive control. From the DFC perspective, the middle thalamic regions show dynamically co-varying patterns in their intrinsic communications, and show high flexibility in the dynamic community analysis. These results indicate that the mid-thalamic region plays a complex and important role in the communication between the nodes of cognitive control-related networks.