Time-resolved detection of stimulus/task-related networks, via clustering of transient intersubject synchronization.

Several methods are available for the identification of functional networks of brain areas using functional magnetic resonance imaging (fMRI) time-series. These typically assume a fixed relationship between the signal of the areas belonging to the same network during the entire time-series (e.g., positive correlation between the areas belonging to the same network), or require a priori information about when this relationship may change (task-dependent changes of connectivity). We present a fully data-driven method that identifies transient network configurations that are triggered by the external input and that, therefore, include only regions involved in stimulus/task processing. Intersubject synchronization with short sliding time-windows was used to identify if/when any area showed stimulus/task-related responses. Next, a first clustering step grouped together areas that became engaged concurrently and repetitively during the time-series (stimulus/task-related networks). Finally, for each network, a second clustering step grouped together all the time-windows with the same BOLD signal. The final output consists of a set of network configurations that show stimulus/task-related activity at specific time-points during the fMRI time-series. We label these configurations: "brain modes" (bModes). The method was validated using simulated datasets and a real fMRI experiment with multiple tasks and conditions. Future applications include the investigation of brain functions using complex and naturalistic stimuli.

INTRODUCTION

Over the last 20 years, functional magnetic resonance imaging (fMRI) has been extensively used to investigate sensory, motor, and cognitive functions in the human brain. The vast majority of these studies made use of experimental designs entailing the manipulation of one (or more) specific stimulus or task parameters to generate a set of relevant experimental “conditions.” Changes of brain activity associated with these conditions are assessed by fitting the BOLD signal with predictors based on the type and timing of conditions (GLM: general linear model) [Friston et al., 1995]. More recently, there has been an increasing interest in the study of brain functioning using more complex and realistic stimuli, such as videos and movies [e.g., Bartels and Zeki, 2004; Fiser et al., 2004; Hasson et al., 2008; Kayser et al., 2004; Nishimoto et al., 2011; Wolf et al., 2010; see also Peelen and Kastner, 2014]. Unlike traditional paradigms, in these cases there is little a priori knowledge about the timing/occurrence of the relevant events and, hence, the construction of “predictors” for data‐fitting is challenging [but see Bartels et al., 2008; Bordier et al., 2013; Lahnakoski et al., 2012a; Ogawa et al., 2013; Raz et al., 2012]. Indeed, most studies using naturalistic stimuli sought to identify patterns of activity based only on the data structure [i.e., data‐driven approaches, e.g., independent component analyses, Bartels and Zeki, 2005; Lahnakoski et al., 2012b; for review see Calhoun and Pearlson, 2012; intersubject correlation analyses, Hasson et al., 2004; cluster analyses, Heller et al., 2006]. Nonetheless, the available data‐driven methods entail several limitations that, here, we seek to overcome with a new approach that combines data‐driven and multivariate clustering techniques.

The data‐driven method that has been most commonly used for the analysis fMRI time‐series is the spatial Independent Component Analysis (ICA) [Calhoun et al., 2001; Comon, 1994]. Spatial ICA identifies functional networks that are spatially independent patterns [Mckeown and Sejnowski, 1998]. Using spatial ICA, a large number of functional networks have been identified even in the absence of any stimulus or task [i.e., ICA of resting‐state data, see Damoiseaux et al., 2006; Veer et al., 2010]. Most of these components/networks are present also in task‐related fMRI time‐series [see Smith et al., 2009]. The latter can make it difficult to decide whether any given IC is related to the task or, instead, reflects some task‐unspecific underling functional coupling (unless one utilizes additional information about the experimental design, e.g., by testing for correlations between the IC time‐courses and GLM predictors). Within the ICA framework, one method that may best identify dynamic and overlapping networks associated with the processing of stimuli/tasks is tensor independent component analysis (tensor‐ICA) [Beckmann and Smith, 2005]. The algorithm is based on parallel factor analysis [Harshman, 1970] and decomposes the data in independent components (ICs), each including a “temporal,” a “spatial,” and a “subjects” mode. The subjects‐mode indexes the consistency of the IC across subjects and has been primarily exploited for multi‐subjects IC analyses or to perform between‐groups comparisons [e.g., Rombouts et al., 2009]. However, it should be noticed that ICs with a high subjects‐mode will comprise areas that share similar time‐courses across subjects. The latter will occur when a network responds to some common external input or task, assuming that all subjects were presented with the same sequence of events (cf., intersubject synchronization, below). Accordingly, tensor‐ICA should enable identifying overlapping spatiotemporal networks that are specifically associated with stimulus/task processing and exclude any spontaneous fluctuations that will not be “synchronized” across different subjects.

A different approach for the identification of stimulus/task‐related brain activity in the absence of any a priori information was proposed by Hasson et al. (ISS: intersubject synchronization) [Hasson et al., 2004]. The ISS works by computing, voxel‐by‐voxel, the correlation between the BOLD time‐series acquired in different participants. Only brain regions that are involved in stimulus processing should show correlated activity across subjects, assuming that all subjects were presented with the same stimuli. The method is fully data‐driven and does not require any a priori knowledge about the stimuli [see also Pajula et al., 2012]. The ISS has been used to investigate brain activity associated with several different processes, including memory encoding [Hasson et al., 2008], auditory oddball task [Hejnar et al., 2007], or narrative speech comprehension [Wilson et al., 2008]; and it is often used to examine brain associated with complex, naturalistic stimuli like movies [Hasson et al., 2004; Golland et al., 2007; Jaaskelainen et al., 2008; see also Hasson et al., 2010, for review]. However, unlike tensor‐ICA, the ISS is a mass uni‐variate approach that considers every single voxel/region in isolation and does not provide us with any information about possible inter‐regional effects (functional networks).

A further issue with the available data‐driven methods, including both tensor‐ICA and ISS, is that they are generally applied to the entire fMRI time‐series. With this, there is an implicit assumption that the relationship between the stimuli and brain activity (ISS) or connectivity (tensor‐ICA) is stable over time. This point has been addressed by studies that used sliding windows to capture the dynamic nature functional connectivity, primarily using resting‐state datasets [see Hutchison et al., 2013a; for review]. Sliding windows have been combined with various methods, with the aim of assessing whether a given index of connectivity changes over the duration of the experiment. These indexes may be derived simply by correlating the BOLD signal in two brain areas; or can include more complex measures such a inter‐regional correlation matrices (tracking the connectivity between a seed‐region and an extensive set of other regions) [e.g., Allen et al., 2014] and networks derived from spatial ICA [e.g., Kiviniemi et al., 2011; Morton and Hutchison, 2014; see also Chang and Glover, 2010; Grigg and Grady, 2010; Majeed et al., 2011, for other indexes of inter‐regional connectivity]. Even at rest, these approaches revealed dynamic reconfigurations of brain networks, reflecting both the fluctuation of high‐level processes (e.g. arousal, attention, memory, etc.), as well as spontaneous ongoing activity [Hutchison et al., 2013b; see also Deco et al., 2011; Sadaghiani and Kleinschmidt, 2013, for reviews].

The relevance of these fluctuations for stimulus processing and/or task performance has been typically assessed by means of correlations between the chosen connectivity index and some behavioral measure. For example, using short time‐windows, Thompson et al. showed that the level of (anti‐)correlation between “the default mode network” and the “task positive network” in a 12.5 s window centered just before the target onset predicted the speed of target detection on a trial‐by‐trial basis [Thompson et al., 2013]. Recently, Raz et al. used 30 s sliding‐windows to extract dynamic patterns of inter‐regional connectivity during movie watching [Raz et al., 2012]. The connectivity indexes were correlated with behavioral and parasympathetic measures related to the subjects' emotional state, on a window‐by‐window basis. The results showed significant correlations between the connectivity of the limbic system with the prefrontal cortex and the sadness ratings, demonstrating the functional relevance of the dynamic connectivity changes [Raz et al., 2012; see also Baldassarre et al., 2012 and van den Heuvel and Sporns, 2013] (who demonstrated the task‐relevance of resting‐state networks using inter‐individual differences, rather than within‐scan time‐resolved measures).

Here, we propose combining ISS with sliding‐windows and cluster analyses with the aim of capturing dynamic inter‐regional effects that are relevant for stimulus/task processing. Cluster analyses have been extensively used to analyze resting‐state data for brain parcellization [e.g., Craddock et al., 2012; Thirion et al., 2014] and to index connectivity at rest [Cordes et al. 2002]. Moreover, time‐resolved applications of cluster analysis contributed to the identification of networks' dynamics at rest [Allen et al., 2014; Liu et al., 2013; Yaesoubi et al., 2014]. Using sliding‐windows and hierarchical clustering, Yang et al. [2014] reported that the duration of specific “network states” correlated with individual neuropsychological scores, again indicating a possible relevance of networks' fluctuations at rest (see also above). Using a task‐based paradigm, Ren et al. [2014] compared hypothesis‐driven GLM analyses and data‐driven clustering algorithms. Overall the results were consistent across methods (with a major advantage in terms of computation times, for the new clustering algorithm proposed by the authors), but the study examined a very simple On‐Off visual activation paradigm and did not include any time‐resolved analysis. To our knowledge, only one previous work combined clustering algorithms with ISS analysis. Kauppi et al. [2010] used a modified version of the k‐means clustering algorithm to analyze ISS coefficients obtained while subjects were presented with complex stimuli (movie watching). However, the aim of the clustering procedure was to assess the reliability of the ISS coefficients across pairs of subjects, rather than identifying dynamic inter‐regional effects (functional networks).

Here, we use ISS with sliding windows to target brain areas where activity is stimulus/task‐dependent. Next, we use a first clustering step to identify sets of regions that engage concurrently and in a repeated manner, over the duration of the entire experiment (ISS‐clustering). This provides us with sequences of time‐points when specific groups of brain areas become engaged (network configurations), without making use of any a priori information about the stimuli. Because the first clustering is based on the significance of ISS, rather than the BOLD signal, the brain activity can differ between areas belonging to the same network and between successive time‐points. Accordingly, under the hypothesis that each network configuration may express different temporal patterns of BOLD activity during the experiment [see also Liu and Duyn, 2013], we apply a second cluster analysis to group time‐windows with the same signal (BOLD‐clustering). The output of the algorithm consists of stimulus/task‐related networks, each including: (1) a set of brain areas; (2) a sequence of time‐points when the network is engaged; and (3) a distinctive BOLD signal for each area belonging to the network. We label these network configurations: “brain Modes” (bModes). The bModes represent networks of brain areas (where) that engage transiently at specific time‐points (when) during the fMRI time‐series. It should be noted that, while the method was explicitly designed to identify “when and where” brain activity is stimulus/task related, the method does provide us any information about “what” in the stimuli or tasks triggers the ISS and generates the bModes, please see also Discussion section.

We first present a detailed description of the method, and then apply the method both to simulated fMRI time‐series and to a real whole‐brain fMRI experiment. It should be noted that the method was developed primarily to investigate brain activity using complex stimuli (e.g., during movie watching), but there are no well‐established methods for assessing task‐related inter‐regional dynamics using such complex stimuli. Therefore, we validated the method with a protocol that included “known” sequences of events and that could be analyzed also with a standard hypothesis‐driven GLM [i.e., data fitting with a priori defined predictors; e.g., see also Pajula et al., 2012, who resorted to an analogous approach for ISS validation; and Calhoun et al., 2001, for the evaluation of spatial and temporal ICA]. We selected particular stimuli and tasks, and combined these in a specific manner, with the aim of generating patterns of coactivation across brain regions (network configurations) that changed dynamically over time (see Fig. 2A1,B1). We used a linguistic task and a spatial judgment task, both entailing visual stimuli and manual responses. With this, we sought to activate task‐common regions in visual and motor cortices, as well as task‐specific regions in associative areas of the left and right hemisphere. Moreover, we presented irrelevant auditory stimuli that were sometime time‐locked to the spatial judgment task and sometime presented in the rest period between the task blocks (see Fig. 2A1). These auditory stimuli will engage the auditory cortex that, thus, should sometime coactivate with the “task‐common” visual/motor areas and with any “task‐specific” region associated with the spatial judgment task, and sometime not (i.e., when the sounds were presented in the rest periods). Accordingly, despite not using naturalistic stimuli, the validation protocol included relatively complex combinations of events, each expected to trigger a specific network configuration. The main goal of the method was to retrieve the different network configurations and their times of occurrence within the fMRI time‐series, in a fully‐data driven manner. For comparison, the same dataset was also analyzed with tensor‐ICA that should also identify task‐related networks, but without providing us with any direct information about the timing of the different stimuli/tasks in the time‐series.

Figure 2

A1. Experimental protocol used both for the simulations and for the main fMRI experiment. Two conditions were presented according to a block‐design and corresponded to the language task (T1) and the spatial task (T2) during fMRI, see Methods for details. The third condition comprised single events (T3) that in the fMRI experiment corresponded to the presentation of task‐irrelevant sounds. Note that the sounds were sometime presented in isolation while sometime they were coupled with the onset of a T2 block (T3 events marked with asterisks). B1. Areas and signals included in the simulations. The simulations comprised nonoverlapping (task‐specific) blocked responses in areas X and Y, “common” activation in area Z; transient responses in area W (sometime time‐locked with the block‐onsets in area Y, cf., asterisks) and the control area K, without any signal. B2. Results of one of the simulations (18 subjects, SNR = 15, time‐window = 16 s, P‐corr. = 0.05, see main text). The transient ISS analysis, with the binary matrix associated with the simulated dataset. For each area, the black lines represent “when” the transient‐ISS was significant. B3. Results of ISS‐clustering step that identified four clusters (A–D), each including one or more areas and engaging at specific time‐points (cf., black lines). C1. The results of the transient ISS (i.e., binary matrix) for the main fMRI experiment. This comprises 116 rows, corresponding to the AAL areas, and 622 columns corresponding to the concatenation of two fMRI runs (see also Methods). C2. The ISS‐clustering step of the main fMRI experiment identified four clusters, each repeating several times over the entire fMRI time‐series (see the black lines associated with each cluster).

Figure 1

Schematic representation of the main steps for the computation of the brain modes (bModes). This comprised: (1) Transient intersubject synchronization (transient‐ISS); (2) ISS‐clustering; and (3) BOLD‐clustering. In the transient‐ISS, separately for each area, sliding windows are used to compute a linear regression between the BOLD signal in each pair of subjects. Nonparametric methods are then applied on the resulting parameter estimates (betas) to determine whether the ISS in each area and each time‐window is statistically significant (see Methods). The result is a binary matrix of functionally relevant areas and time‐windows (the significant points are displayed in black). The binary matrix is then submitted to the ISS‐cluster analysis (k‐mean). The result of this is a set of clusters/networks, each with a distinctive set of time‐windows indexing when the cluster was engaged (the relevant time‐points are displayed in black). Separately for each cluster, the BOLD signal of each area and time‐window is submitted to the BOLD‐clustering step that allows us to group together time‐windows with similar BOLD signal. The final output consists of task‐related networks that are engaged at specific time‐points during the fMRI time‐series and that display a specific BOLD signal (i.e., the bModes).

MATERIAL AND METHODS

Computation of the bModes

The computation of the bModes comprises three main steps (see Fig. 1): (1) ISS using short sliding‐windows, with the aim of detecting stimulus/task‐related transients changes of the BOLD signal (transient intersubject synchronization [transient ISS]); (2) Clustering of the synchronized brain areas, with the aim of identifying network configurations comprising areas that became engaged at the same time (ISS‐clustering); and (3) Clustering according to the BOLD signal in the different areas of each cluster, with the aim of grouping time‐windows with similar BOLD signal (BOLD‐clustering).

Transient ISS

The aim of the first step of the analysis is to identify time‐points within the fMRI time‐series when there is a stimulus/task related change of the BOLD signal. To do this, we perform a variant of the intersubject correlation analysis proposed by Hasson et al. [2004] but using short time windows [see also Glerean et al., 2012].

Before computing the transient‐ISS, the fMRI time‐series were averaged within regions of interest, here using AAL defined regions and generating 116 time‐series for each subject [Tzourio‐Mazoyer et al., 2002] (please see also the Discussion section about the possibility of using other templates). Next, a multiple regression model was used to high‐pass filter the data and to remove signal related to head movements. For this, the model included a cosine basis‐set and the movement parameters estimated during realignment [Friston et al., 1995]. All the analyses presented below made use of a high pass filtering cut‐off = 0.008 Hz. These settings did not affect the sustained signal associated with the task‐blocks (cf., Simulations section) and were found not to produce any false positives in simulated datasets, irrespective of the presence/absence of low‐frequency drifts (not shown). However, higher cut‐off may be considered [i.e., related to the length of the sliding‐window: e.g., “1/window‐length,” see Leonardi and Van De Ville, 2015] particularly if the dataset is not expected to contain any regions with sustained periods of activity. At this stage, other covariates of no‐interest could be added, such as the global signal in each volume or covariates related to the individual behavior (e.g., eye‐movements recorded during fMRI scanning).

Independently for each AAL region, the time‐series were divided into short segments using sliding windows (sliding step = 1 pt) [Gembris et al., 2000]. For the analyses of the real fMRI experiment, the size of the window was set to be equal to the duration of the hemodynamic response function (HRF), that is, approximately 16 s [cf., also Glerean et al., 2012]. The appropriateness of this choice was confirmed by our simulations that systematically assessed the effect of the window‐size (see below, and Fig. 4B). For each region and each segment, a linear regression was used to compute the “synchronization” between each pair of subjects, that is, the data of one subject acted as the predictor for the data of a different subject. The statistical assessment of the resulting (nonindependent) parameter estimates was performed using a nonparametric method. We sampled 1,000 bootstrap means [Efron and Tibshirani, 1998] and estimated the correspondent density function [kernel density estimation: Bowman and Azzalini, 1997]. The inverse of the probability function was used to assess whether the ISS was significantly larger than zero with a 5% probability threshold, corrected for the number of sliding windows (see also below and Fig. 4B, for simulations concerning the thresholds applied to the transient‐ISS analysis). This procedure was repeated for all regions and all sliding windows, generating a binary matrix that represents “where” (areas, rows) and “when” (windows, columns) the BOLD signal shows a transient burst of task‐related activity (see Fig. 1, panels on the left).

Figure 4

Results of the simulations using different parameters, see also main text. A. The Matthews correlation coefficient (Mcc) plotted against the signal‐to‐noise ratio, for different sample sizes (time‐window = 16 s; ISS‐threshold set to P‐corr = 0.05). B. The Mcc computed for sliding time‐windows of different length and with the ISS‐threshold set at different levels of significance (18 subjects, SNR = 15). Each simulation condition was run five times and the error bars represent the standard deviation of the Mcc across these five repetitions. The P‐values refer to the thresholds applied to the transient‐ISS analysis, with (P‐corr) or without (P‐unc) correcting for the number of sliding windows, see also Method section.

Clustering of the “synchronized” brain areas (ISS‐clustering)

The second step of the analysis aimed to group areas that show similar patterns of transient synchronization during the fMRI time‐series. For this, the ISS binary matrix (cf., Step 1, above) was submitted to a k‐mean cluster analysis [Seber, 1984]. The clustering procedure minimizes the distance between k centroids and the data. Because of the binary format of the data, the Hamming method—which quantifies the number of bits that differ between two sequences—was used to compute the distances [Hamming, 1950]. The number of centroids (k) was determined using the elbow method that takes into account the increase of the percentage of variance explained when the number of centroids is progressively increased [Thorndike, 1953]. This method determines the number of clusters, so that adding another cluster does not improve the modeling by more than 5%. Across all the simulations (290 in total, see below), k was found to range between 0 and 8; in the main fMRI experiment k was found to be equal to 4. To minimize the effect of the initial seed, the clustering was repeated n times using different seeds (n = number of windows divided by k), and the clustering with the lowest sum of square distance was retained as the final result. For instance, the final results for the fMRI experiment were obtained by running the k‐mean cluster analysis 153 times (i.e., number of windows = 615; k = 4; n = 153).

This second step of the analysis provided us with transient network configurations, comprising areas that became engaged at the same time and in a repeated manner. Each cluster comprises of a set of areas and a sequence of relevant time windows (see Fig. 1, center panel). The latter represents “when” the network was engaged in the fMRI time‐series.

Clustering according to the BOLD signal (BOLD‐clustering)

The ISS‐clustering procedure generated groups of areas and windows (clusters), irrespective of the BOLD signal in each area and each window: only the timing and the consistency of the synchronization across subjects contributed to the definition of these inter‐regional patterns. This may include the simultaneous coactivation of several regions; but also anti‐correlations between regions, with the BOLD signal increasing in one area and decreasing in a different area (see also Discussion section). The aim of the last step of the analysis was to further partition each cluster according to the BOLD signal in the relevant windows.

For each cluster, we extracted the BOLD signal of each area (mean across subjects) and each time window. The signals of all the areas belonging to the same cluster were concatenated, separately for each window (see Fig. 1, panels on the right). For each cluster, this generated a W×T matrix, where W is the number of windows and T is the number of areas by the number of data points per window. Accordingly, each row of a W×T matrix can include data belonging to different areas. A new k‐mean cluster analysis was then performed on each of these matrices, with the aim of grouping windows with similar patterns of activity. Because the W×T matrix is not binary, the data‐to‐centroids distances were now computed using the square Euclidean method. The number of centroids was estimated using the Dunn's index [Dunn, 1974; see also Pal and Biswas, 1997] that aims to identify dense and well‐separated clusters. This index can range between 0 and infinite and is defined as the maximization of the ratio between the smallest inter‐cluster distance to the largest intra‐cluster distance. The Matlab‐code to compute the brain Modes is available at: http://www.slneuroimaginglab.com/mt-tools.

Simulations

The method was evaluated using simulated fMRI time‐series and varying several parameters, see below. The signal/paradigm was identical in all the simulations and comprised five areas (X, Y, Z, W, K) with specific patterns of activation (see Fig. 2B1). Area X and Y included 10 blocks of nonoverlapping sustained activation, while area Z included all 20 blocks simulated in X and Y. Thus, areas X and Y simulated regions responding selectively to one or another stimulus/task, while area Z simulated a “common” region responding to both tasks. Area W included 16 transient events, 6 of which were time‐locked with the onset of the activation blocks in area Y (cf., events marked with an asterisk, in Fig. 2A1). Area K served as a control region and did not include any signal. The time‐series comprised 622 data points (volumes) that would correspond to two fMRI‐runs of approximately 10 min each, with a repetition time of 2 s; cf., also main fMRI experiment below.

We used the R package “neuRosim” to simulate the BOLD signal in each area [Welvaert et al., 2011]. We manipulated the number of subjects (sample size), the signal‐to‐noise ratio (SNR), the size of the sliding‐windows, and statistical threshold applied to the ISS analysis. For all the simulations, the baseline value of the time‐series was set to 100 and the response amplitude was set to 1% [see Welvaert and Rosseel, 2013]. The signal in each region and each subject was convolved with the HRF. Rician noise was added to the simulated time‐series, independently for each subject and area. The simulated data neither included any slow drift component nor simulated physiological noise. The average SNR was defined as:

With  = the average magnitude of the signal; and  = the standard deviation of the noise [Kruger and Glover, 2001].

In the first set of simulations, we varied the number of subjects and the SNR. The sample size was set equal to 8, 14, 18, or 25 subjects and the SNR ranged between 1 and 30; see Figure 4A. For this set of simulations, the size of the sliding windows was equal to 16 s and the results were evaluated using an ISS‐threshold of P‐corr. = 0.05. A second set of simulations varied the size of the sliding windows between 12 and 40 s (i.e., 6–20 TRs) and the ISS‐threshold, see Figure 4B. The ISS‐threshold was set to P‐corr. = 0.05, P‐corr. = 10−4, or P‐uncorr. = 0.05. Please note that here setting the threshold to P‐corr. = 10−4 corresponds to applying a Bonferroni‐type correction for multiple comparisons considering the number of sliding windows, as well as 200 areas (i.e., more than the number of ALL regions used for the analysis the fMRI experiment). For this second set of simulations, the sample size was fixed to 18 subjects and the SNR was set to 15.

Each simulation condition was repeated five times, thus 200 different datasets were constructed for the first set of simulations (4 “simple sizes” x 10 “SNR” x 5 “repetitions,” Fig. 4A) and 90 for the second set (6 “windows” × 3 “thresholds” × 5 “repetitions,” Fig. 4B). For each single simulation, we quantified the accuracy of the output using the Matthews correlation coefficient (Mcc) [Matthews, 1975]. This quantifies the performance of a binary classification by comparing the “true” versus the “observed” classification. Importantly, this method takes into account both correct detections as well as any false positive. Specifically, the Mcc is defined as:

With TP corresponding to the number of true positive (hits), TN of true negative (correct rejections), FP of false positive (false alarms), and FN of false negative (misses).The Mcc‐values range between +1 and −1, indexing a perfect match versus a total disagreement between the true and the observed classification. Values around zeros indicate random performance of the classification procedure.

Here, for each simulation we obtained a set of bModes, each including a number of areas and a sequence of time‐points indicating when the bMode was active. We scored as a TP each time‐point of a bMode, when both the timing and the combination of areas were correct. The timing was considered “correct,” if the true signal happened within the size of the sliding‐window used for the simulation. Any bMode occurrence without corresponding true signal was scored as FP. If the method failed to detect a simulated combination of areas at a given time‐point, this was counted as a FN. Finally, the number of TN was obtained by subtracting the total number TP + FP + FN from the total number of sliding windows (622 in our simulations). It should be noticed that at low levels of SNR and with small sample sizes, there was no significant ISS and the Mcc could not be computed (TP = 0, FP = 0, with the Mcc‐denominator = 0). In these cases, the Mcc was set equal to 0. Table 1 shows the results of one of the simulations, with the “true” simulated signal and the corresponding counts of TP, TN, FP, and FN.

Table 1

Summary of the results obtained for one of the simulated dataset (number of subjects = 18; SNR = 15; size of the sliding time‐windows = 16 s)

bMode results		Simulated signal					Evaluation
bModes		Areas				Nr.	True		False
Label	Areas	Z	X	Y	W		Positive	Negative	Positive	Negative
A1‐6	W	—	—	—	Ev.	10	10	—	0	0
B1‐2‐3	YZW	B‐on	—	B‐on	Ev.	6	6	—	1	0
C2	XZ	B‐on	B‐on	—	—	10	10	—	0	0
C1	XZ	B‐off	B‐off	—	—	10	10	—	0	0
D2	YZ	B‐on	—	B‐on	—	4	4	—	0	0
D1	YZ	B‐off	—	B‐off	—	10	9	—	0	1
—	—	—	—	—	—	—	—	571	—	—
						Total	49	571	1	1
							Mcc = 0.978

On the left, the results/output of the new method (see also Fig. 3). In the center, the “true” simulated signal with the pattern of activation in the different areas (B‐on: block‐onset; B‐off: block‐offset; Ev. event, see also Fig. 2A1). On the right, the evaluation the results that took into account the True positives/negatives (hits and correct rejections), and the False positives/negatives (false alarms and misses). The count of the true negative corresponds to the number of sliding windows (622) minus the sum of true positive, false positive, and false negative. With these counts we computed the Matthews correlation coefficient (Mcc) that provided us with a quantitative measure of the accuracy of the results (see also main text).

Figure 3

Final bModes for one of the simulations (see Fig. 2B). The BOLD‐clustering allowed us to group together time‐windows with the same BOLD signal, thus generating the final bModes. Each bMode comprises a set of areas, a set of time‐windows and a specific BOLD signal for each area belonging to the bMode. We also show the “true” simulated signal (cf., Fig. 2B1) that here was averaged using the time‐windows associated with each bMode, see plots in the insets. The latter demonstrate that the method successfully identified and grouped together time‐points corresponding to specific events (or combination of events) in our simulated dataset.

Main fMRI Experiment

Participants

Nineteen subjects (aged 20–40, mean = 28.1 years, 11 females and 8 males) with no history of neurological or psychiatric illness participated in this study. They had normal or corrected‐to‐normal visual acuity and reported no difficulty of hearing. One subject was discarded from the analysis because of large head movements during the data acquisition, leaving 18 participants for the analyses. The Ethical Committee of Santa Lucia Foundation approved this study. All subjects gave written informed consents prior to the scanning session.

Paradigm

The aim of the current method is to detect changes of network configurations that are related stimulus and task processing, and to do this in a fully data‐driven manner. To induce these changing patterns inter‐regional of activation (network configurations), we made use of three tasks/stimuli with specific characteristics and presented with specific temporal relationships with each other. First, we selected two tasks that should activate both “common” brain regions, as well as “task‐specific” areas. We made use of a language task and a spatial perception tasks (see section below, for details). Both tasks involved visual stimuli and motor responses, and should therefore activate common regions in visual and motor cortex. In addition, the language task should activate frontotemporal regions in the left hemisphere [Binder et al., 1997]; while the spatial task should activate frontoparietal regions in the right hemisphere [Fink et al., 2001; see also Pugh et al., 1996]. The two tasks were presented in different blocks. We expected to identify network configurations including both task‐common and task‐specific areas, with sequences of time‐windows consistent with the timing of corresponding tasks: for example, visual and motor areas together with language regions in the left hemisphere, when the subjects performed the language task; and different network configurations including the same visual and motor regions, but now together with areas in the right hemisphere, during the visuospatial judgment task.

Moreover, the experimental paradigm included task‐irrelevant auditory events. These will activate regions in the auditory cortex, not involved in the two other tasks, and should not activate any visual or motor region. Accordingly, these stimuli will engage a third, independent set of areas (auditory‐network). However, we induced the coactivation of this auditory‐network and the spatial‐task network by controlling the temporal relationship between the onset of the spatial‐task and auditory events: sometime the sounds were presented in isolation (i.e., during the rest period between the task‐blocks), while other times they were time‐locked with the onset of the spatial‐task (see also Fig. 2A1, events marked with an asterisk). Thus, we expected to identify a network configuration comprising auditory regions alone, with a sequence of time‐windows corresponding to the events presented in isolation; and a different configuration comprising both auditory areas and the areas involved in the spatial judgment task, now with time‐windows corresponding to the sounds co‐occurring with the spatial‐task.

In sum, our fMRI paradigm was designed to trigger several different patterns of coactivation at specific time‐points during the experiment. These included: (1)Visual and motor regions with left‐hemisphere areas (language task); (2) Visual and motor regions with right‐hemisphere areas (spatial‐task, without sound); (3) A similar network, but also including auditory areas (spatial task with co‐occurring sounds); and (4) Auditory areas alone (sound only). With the new method, we sought to isolate these different inter‐regional patterns and to recover the corresponding timings in a fully data‐driven manner.

Stimuli and tasks

The fMRI protocol comprised two active tasks that were presented in blocks of 30 to 36 s, plus task‐irrelevant auditory events. The first task required the semantic categorization of written words [Leube et al., 2001]. On each trial, a word was presented on the screen for 500 ms and the participants were asked to report whether this belonged to a “living” (animals) or a “non‐living” (vehicles, tools or furniture) category. Participants responded by pressing one of two buttons with the index/middle finger of the right hand. Each block included 16–19 trials (50% living and 50% non‐living), with a fix inter‐trial interval of 1500 ms. The same pool of words was used for all participant (87 for each category), but the order of presentation was randomized across participants. The second task was a spatial bisection judgment [Fink et al., 2001]. Each trial consisted of an horizontal line (length =6°–7° visual degrees) that was pre‐bisected by a shorter vertical line (2.5°). On half of the trials, the bisection was “symmetrical” with respect to the horizontal line; while in the other half of the trials the bisection was “asymmetrical” with the vertical line intersecting the horizontal line either to the left or to the right of the center. The stimulus display was presented for 500 ms, followed by a fix intertrial interval of 1500 ms. The task of the participants was to report whether the bisection was symmetrical or asymmetrical. Participants responded by pressing one of two buttons with the index/middle finger of the right hand. Each block lasted for 30 to 36 s and included 16–19 trials.

Each fMRI‐run included five blocks of the words‐task and five blocks of the spatial‐task, presented in a pseudo random order. Together with the two active tasks, during each fMRI‐run we presented eight task‐irrelevant auditory events, consisting of bursts of white‐noise (duration = 1 s). Three of these auditory events were time‐locked to the onset of a spatial‐task block, while the other five were presented during the interblock interval. Each participant underwent two fMRI‐runs lasting approximately 10 min each. The data analyses considered both runs together, thus including 10 blocks of the words‐task, 10 blocks of the spatial‐task, and 16 auditory events (10 in isolation and 6 coupled with the spatial‐task, see Fig. 2A1). The most important constraint to perform any ISS analysis is that all subjects are presented with the same stimuli [Hasson et al., 2004]. Accordingly, the order and the timing of the three conditions (i.e., the blocks of the two tasks and the auditory events) were identical for all 18 subjects.

FMRI acquisition and preprocessing

A Siemens Allegra (Siemens Medical Systems, Erlangen, Germany) 3T scanner equipped for echo‐planar imaging (EPI) was used to acquire the functional resonance images. A head‐sized quadrature volume coil was used for radio frequency transmission and reception. Mild cushioning minimized head movement. Thirty‐two slices volumes were acquired using blood oxygenation level dependent contrast (192 mm × 192 mm × 120 mm, in‐plane resolution = 64 × 64, pixel size = 3 mm × 3 mm, thickness = 2.5 mm, 50% distance factor, TR = 2.08 s, TE = 30 ms), covering the entire cortex. The total duration of each functional run was 10 min and 55 s, with the acquisition of 315 volumes.

The fMRI data were preprocessed using SPM8 (Wellcome Trust Centre for Neuroimaging, London, UK). After discarding the initial four volumes of each run, the remaining volumes were slice‐time corrected, head‐motion realigned and normalized to the standard MNI EPI template space (voxel‐size resampled to 3 × 3 × 3 mm3). Finally, the data were spatially smoothed with a 8 × 8 × 8 full‐width at half‐maximum Gaussian kernel.

GLM analysis

First, the imaging data were analyzed using a standard hypothesis‐based approach in SPM8. For each subject, a GLM was used to fit the time‐series at each and every voxel. The regression model included four predictors of interest corresponding to the two active tasks, the sound alone and the sounds coupled with the onsets of the spatial‐task. Each model included also subject‐specific realignment‐parameters and the session constants, as effects of no interest. The task‐blocks (words and spatial) were modeled with a boxcar function time‐locked to the onset of the task with a duration corresponding to the block duration (15–18 TR). The auditory events were modeled with delta functions (duration = 0). The predictors were convolved with the canonical HRF and the time‐series were high‐pass filtered at 0.008 Hz. For each subject, we computed five contrast‐images corresponding to the parameter estimates for the four conditions of interest averaged across the two fMRI‐runs, and the difference between the “words” and the “spatial” tasks. Group‐level analyses were performed with five separate t‐tests. One model was used to directly compare the two active tasks (see Fig. 5A,B), while the other four tests identified activations and deactivations for each condition compared to rest (Fig. 5C–F). The statistical thresholds were set to voxel‐level P‐FWEcorr. = 0.05, corrected for multiple comparisons considering the whole brain as the volume of interest. In addition, the contrast‐images associated with the four conditions were also entered in an omnibus one‐way ANOVA. This was used to compute an overall F‐map highlighting all the areas where the GLM fitted the data, irrespective of condition and activation vs. deactivation (see Fig. 7A).

Figure 5

Results of the standard, hypothesis‐driven analyses of the main fMRI experiment. A. Direct comparison of “Words > Spatial.” The threshold was set to P‐FWE‐corr. = 0.05. B. Direct comparison of “Spatial > Words.” The threshold was set to P‐unc. = 0.001, for illustrative purposes (see Table II, with the statistics corrected for multiple comparisons). C. “Words vs. Rest”; D. “Spatial vs. Rest”; E. “Sound (alone) vs. Rest”; and F. Sounds time‐locked with spatial block‐onset “(Sound & Spatial) vs Rest.” Activations are shown in red, deactivations (i.e., Rest > stimuli/task) in blue. For C–F, the thresholds were set to P‐FWE‐corr. = 0.05. Activations are rendered on the MNI template in SPM8.

Figure 7

A. Results of the omnibus F‐test that assessed the overall fitting of the GLM model, including all tasks and stimuli. The threshold was set to P‐FWE‐corr. = 0.05. B. Projection of all the bModes (clusters A–D) that were detected in the main, task‐related fMRI experiment. C. Projection of the 36 IC‐maps with temporal modes that correlated at least with one of the four GLM predictors (cf., also Supporting Information Fig. S2). D. Projection of the 24 IC‐maps with temporal modes that did not correlate with any of the GLM predictors (see also Supporting Information Fig. S2). Activations, bModes, and IC‐maps are rendered on the MNI template in SPM8.

Tensor‐ICA

Together with the validation of the bModes results using a standard GLM approach, which requires a priori knowledge about the timings and conditions (cf., above), we also sought to compare the current method with another fully data‐driven approach. For this, we submitted our task‐related fMRI data to the “tensor ICA” algorithm [Beckmann and Smith, 2005] implemented in MELODIC v3.14 tool (Multivariate Exploratory Linear Decomposition into Independent Components), part of FSL v.5.0.8 (FMRIB's Software Library, http://www.fmrib.ox.ac.uk/fsl). This implementation of ICA makes the assumption that the temporal response pattern of each IC is the same across the population. The algorithm extracts a set of ICs that characterize the signal variation in the temporal (time‐series), spatial (spatial maps) and subject domains. Here, we concatenated the two fMRI‐runs and submitted the data of all 18 subjects to the algorithm, which initially extracted 198 ICs. The number of components was determined with the “automatic dimensionality estimation” in Melodic, which uses the Laplace approximation to the Bayesian evidence of the model order [Beckmann and Smith, 2004; Minka, 2000]. We retained 60 ICs that showed reliable spatiotemporal components over subjects (i.e., mean of the subject‐mode larger than zero, at P < 0.001) and examined the relationship between these ICs and the experimental design by correlating the temporal mode of each IC with the four GLM predictors, using four separate regressions (see Supporting Information Fig. S2).

RESULTS

Qualitative description of the results for one simulated dataset

First, we describe in details the output of one of the simulations. The simulated dataset included 18 subjects (corresponding to the sample size of the fMRI experiment), a SNR equal to 15 (analogous to the SNR of the real fMRI data; SNR =16.6) [computed according to Constantinides et al., 1997], sliding windows of 16 s and an ISS‐threshold of P‐corr. = 0.05, that is, the same parameters used for the analysis of the main fMRI experiment. The timings of the simulated signals corresponded exactly to the timings of the stimuli/tasks used during fMRI (see Methods and Fig. 2A1).

Figure 2B.2 shows the output of the transient‐ISS analysis. This is a binary matrix with five lines corresponding to the five simulated “areas” (X, Y, Z, W, K) and 622 columns corresponding to the number of the sliding windows applied to the time‐series. The black lines (value in the matrix = 1) indicate “when” and “where” the ISS was significant, thus identifying stimulus/task related areas and the associated relevant time‐points (cf., also Fig. 1). This included the onset and the offset of the task‐blocks (areas X, Y, Z) and the transient events (area W); compare Figure 2B1,B2. It should be noticed that the onset of the significant ISS with respect of the simulated signal, as well as number of successive windows showing a significant ISS (i.e., width of the black lines), was somewhat variable. This was accounted for by the last step of the analysis (BOLD‐clustering), please see below.

Next, the ISS binary matrix was submitted to k‐mean clustering (cf., Fig. 1). This isolated four separate clusters, each including areas with the same pattern of significant ISS over the time‐series. Figure 2B3 shows “when” the clusters were active and the areas belonging to each cluster (where). Cluster A included area W only; cluster B also included area W, but now together with areas Y and Z; cluster C included area X and Z; and cluster D included area Y and Z. The areas belonging to each cluster and the timings of engagement of the clusters were consistent with the simulated signal (see below, Fig. 3 and Table 1).

As noted in the Methods section, the “transient ISS” and the “ISS‐clustering” analysis steps work on the consistency and the timing of the responses, not on the specific shape of the BOLD signal within each window. Thus, a single cluster can include windows with different signals: for example, cluster C included both the onsets and the offsets of the words‐task blocks. The second cluster analysis (BOLD‐clustering, see Fig. 1) grouped together time‐windows with similar signals. This then enabled us to average signals across time‐windows and to obtain the final bModes, each with a specific set of areas and time‐windows, but also with a specific BOLD signal at the times of occurrence. Figure 3 shows the results of this last step of the analysis, with the final set of bModes for this specific simulation. The cluster A was separated in A1 to A6 that differed in time of onset of the BOLD response within the window, with some bMode showing an earlier onset than others (Fig. 3A, for bModes A1–3). This reflects the variability of “when” the transient‐ISS became statistically significant compared to the actual/simulated response onset. When the simulated signal was averaged and plotted using the time‐points associated with the bModes, we observed the very same shift in time (cf., inset, in each panel). This demonstrates that the method identified the correct timing of the stimulus onsets, despite the variability of exactly when the transient‐ISS became significant. The BOLD‐clustering of cluster B showed the same effect (see Fig. 3B). Most importantly, the comparison of Figure 3A,B highlights that the method was able to correctly track the changes of coupling of area W, which sometime activated alone (bModes A1–6) and sometime activated together with the XZ network (bModes B1–3). For clusters C and D, the BOLD‐clustering separated correctly the block‐onsets and the block‐offsets, both for the XZ network (C1 vs. C2, in Fig. 3C) and the YZ network (D1 vs. D2, Fig. 3D). Figure 3C,D show that the method was successfully in tracking the changes of network membership of the “common” area Z, which became part of the XZ network or the YZ network at different moments in the time‐series.

Table 1 presents the quantification of the performance accuracy for this specific simulation. The method successfully detected and classified almost all the expected network configurations, at the correct time‐points of the simulated time‐series. Specifically, all the onsets and the offsets of the XZ network were detected successfully (bModes C1–2). The method discriminated well the presentation of isolated events in area W (see bModes A1–6) from when these events were coupled with the activation of the YZ network (YZW‐network: bModes B1–3). We observed a single case of wrong classification with bMode B3 (areas YZ, plus W) occurring at a time‐point actually corresponding to an onset of the YZ network without any co‐occurring event in area W (see also inset in Fig. 3B). This was scored as a false positive for bMode B3, as well as a false negative for bMode D1; that is the bMode corresponding to the correct YZ network configuration at that time‐point. The resulting Mcc‐value for this specific simulation was 0.978.

In sum, this simulation highlighted that the method was able to retrieve several different network configurations in fully data‐driven manner, and to do this with a high degree of accuracy (Table 1, and see section below). The method identified areas that switched between two different networks (“common” area Z switching between “task‐specific” areas X and Y; see bModes C and D); and areas that dynamically coupled or uncoupled with another network (area W activating alone or coupled with YZ‐network, see bModes A and B). The BOLD‐clustering step separated instances when the same set of regions were coupled, but showed a different BOLD signal (i.e., block‐onset vs. block‐offset, cf., bModes C1–2 and D1–2).

Quantitative evaluation of the method's performance

To quantify the performance of the new method, we manipulated several parameters and evaluated the accuracy of the method's output. The performance was assessed by comparing the simulated patterns of inter‐regional coactivation (“true signal,” see Fig 2B1) with the observed bModes, each including a set of areas and a sequence of time‐points (e.g., see Fig. 3 and Table 1). The accuracy of each simulation was quantified using the Mcc that takes into account the performance both on true positives (Hits and Misses), as well as on true negatives (correct rejections and false alarms), see Methods section.

A first set of simulations varied the sample size (number of simulated “subjects”) and the SNR. The size of the sliding‐windows was set to 16 s and the threshold of the ISS analysis was set to P‐corr. = 0.05. Figure 4A shows the method performance as a function of the sample size and SNR. As expected we found that, across sample sizes, the performance increased with increasing SNRs. The performance reached a plateau of approximately 0.9 with SNR = 8–10 for the largest sample size (n = 25); while for the small sample size (n = 8) this required an SNR = 20–25. With an SNR = 15 (corresponding to the value estimated for our main fMRI experiment, cf., above), a sample size between 15 and 20 subjects appeared suitable to obtain a good performance (Mcc approx. 0.9).

A second set of simulations varied the size of the sliding‐windows and the statistical threshold applied to the transient‐ISS analysis (Fig. 4B). With the current simulated signal, which included both blocks and events (see Fig. 2A1), we found that the best performance was obtained with sliding windows between 8 and 12 TRs (i.e., 16 to 24 s). Possibly, this reflects the best balance between having sufficient data points to compute the ISS reliably and, at the same time, including segments of the time‐series with a relevant change of BOLD signal (i.e., increases/decreases of the signal associated with the block‐onsets/offset or the events). Concerning the ISS‐threshold used to generate the binary matrix, the simulations showed a performance decrease when this was not corrected for multiple comparisons (P‐unc., see black line in Fig. 4B). This was primarily due to an increase of the number of “synchronized” time‐windows appearing in the binary matrix, despite the lack of any true signal in the corresponding areas or time‐points (i.e., Type I errors in the transient‐ISS analysis). This, in turns, generated an high number of false positives that included both the detection of bModes at time‐points without any true signal, as well as the detection of wrong network configurations. Conversely, increasing the corrected threshold from P‐corr. = 0.05 to P‐corr. = 10−4 did not affect the method's performance to any large extent. This may be relevant if one considers correcting the ISS‐analyses not only for the number of windows but also for the number of areas included in the analysis. In this context, a P‐corr. = 10−4 would correspond to a further Bonferroni‐type correction for 200 areas, that is, more than those considered in the main fMRI experiment (116 AAL‐defined regions).

Main fMRI Experiment: GLM Results

The same protocol used for the simulations was applied during whole‐brain fMRI scanning in healthy participants (cf., Fig. 2A1). One task was a category‐judgment language task (“living vs. non‐living” words) and the other task was a spatial judgment task (“symmetrical vs. asymmetrical” bisection lines). The transient events were auditory bursts of white noise. First, the fMRI data were analyzed with a standard GLM that made use of a priori knowledge of “what” stimuli/tasks were presented and “when” (hypothesis‐based approach, see Fig. 2B1).

The direct comparison between the two tasks revealed the expected dissociation, with activation primarily in the left hemisphere, including the left inferior frontal operculum and regions in the left occipital cortex, for the words‐task and in the right hemisphere, including the right superior parietal gyrus, plus right precuneus and right inferior temporal gyrus for the spatial‐task (Fig. 5A,B, and Table 2). When both tasks were compared with rest, we found the expected activation of “common” regions in visual occipital cortex and in the left motor cortex contra‐lateral to the hand used for responding (see Fig. 5C,D, activations shown in red; and Supporting Information, Tab. S1). Moreover, we found significant activation of the left frontal operculum for the words‐task and of the right supramarginal gyrus for the spatial‐task, see Supporting Information Tab. S1. The auditory events activated only the auditory cortex when the sounds were presented in isolation (Fig. 5E); while the sounds time‐locked with the onset of spatial‐task activated the auditory cortex, “task common” regions in visual and motor cortex, plus areas specific for the spatial‐task in the right premotor cortex (see Fig. 5F and Supporting Information Tab. S1).

Table 2

Direct comparisons between the two active tasks in the main fMRI experiment

	x	y	z	z‐sc.	P‐corr.	Vox.
Words > Spatial
L. Inf. Occip.	−15	−102	−7	5.90	< 0.001	18
L. Lingual	−36	−90	−13	5.76	< 0.001	54
L. Inf. Frontal Oper.	−42	15	23	5.22	= 0.005	21
L. Cereb. Crus 1	−45	−66	−28	5.39	= 0.002	22
Spatial > Words
R. Sup. Parietal	27	−72	56	4.82	= 0.035	4
R. Precuneus	3	−51	38	4.78	= 0.040	1
R. Inf. Temporal	57	−66	−10	4.74	= 0.048	1

The contrasts have been assessed with a standard, hypothesis‐driven GLM analysis. Results are reported at P‐FWE‐corr.= 0.05, corrected for multiple comparisons at the voxel level. The labeling of the areas correspond to the AAL template. L/R: left/right hemisphere. Coordinates are in MNI space (x, y, z; in mm). Z‐sc: Z scores. Vox: number of voxels in the activated cluster.

Overall, the GLM analyses confirmed that our experimental paradigm was successful in generating: “common” activation for the two active tasks (visual and motor areas); “task‐specific” activation for the words‐task (inferior fontal regions in the left hemisphere) and for the spatial‐task (premotor and parietal regions in the right hemisphere); and activations associated with the auditory events, including either only the auditory cortex (isolated sounds) or the auditory cortex plus visuo‐motor and right premotor cortex when the sounds were time‐locked to the onset of the spatial‐task. In addition, the GLM analyses revealed that during the two active tasks the BOLD signal decreased below rest‐baseline in several regions (see Fig. 5C,D, de‐activations shown in blue). For both tasks, this included the medial posterior parietal cortex (cf., default mode network) [Raichle et al., 2001].

Main fMRI Experiment: bModes

The time‐series of the main fMRI experiment were averaged within the 116 AAL regions and submitted to the new method. Visual inspection of the transient‐ISS binary matrix already indicated some regular pattern of ISS over both areas and time‐windows (Fig. 2C1). This was formally assessed with the first clustering step. The ISS‐clustering identified four clusters (A, B, C, and D; see Fig. 2C2), plus one additional cluster that occurred only once during the entire experiment and was not considered further. Cluster A comprised four areas corresponding to the auditory cortex bilaterally (see also Fig. 6A). Cluster B comprised 13 areas, including the auditory cortex but also the visual and motor areas, and frontoparietal regions in the right hemisphere (Fig. 6B). Cluster C comprised 11 areas, including visual and motor regions, the left premotor cortex, plus medial parietal cortex, and dorsal occipital regions (Fig. 6C). Cluster D included only the visual and motor cortices (Supporting Information Fig. S1; and Table 3, with the list of the AAL regions included in each cluster). Overall, there was a good correspondence between the areas identified with our data‐driven approach and those highlighted with standard data‐fitting GLM (compare panels A and B, in Fig. 7).

Figure 6

Final bModes obtained for the main fMRI experiment. The figure shows three out of the four clusters identified in main fMRI dataset, the fourth cluster is reported in the Supporting Information. For each cluster, we show the relevant areas rendered on the MNI brain template; and the further subdivision in the final bModes, each engaging at specific time‐points (see plots below the 3D renderings) and with a specific BOLD signal (see plots on the left). The areas belonging to each cluster are listed in Table III. Above each signal plot, we report how many times the bMode occurred during the time‐series (“n” time‐windows). The BOLD signal is an average of the signal in these time‐windows, separately for each area belonging to the bMode. To highlight that areas belonging to the same bMode can display substantially different BOLD signal, for bModes B1‐3 and C1‐2 the signal of different areas are plotted in different colors (cf., also Table III). The insets show the expected signal for the different tasks/stimuli (i.e., GLM design matrix, and cf., Fig. 2B1) averaged using the time‐windows associated with the corresponding bMode. These highlight that we successfully identified and grouped together time‐points corresponding to specific events in the experimental paradigm, see also main text.

Table 3

The areas belonging to the four clusters identified in the main fMRI experiment

	Cluster A		Cluster B		Cluster C		Cluster D
Visual			L. Sup. Occip.	(r)	L. Sup. Occip.	(g)
					L. Inf. Occip.	(y)	L. Inf. Occip.	(SI)
					R. Inf. Occip.	(y)	R. Inf. Occip.	(SI)
Motor			L. Precentral	(r)	L. Precentral	(y)	L. Precentral	(SI)
					R. Cereb. 6	(y)	R. Cereb. 6	(SI)
					L.SMA	(y)
Auditory	L. Heschl	(b)	L. Heschl	(b)
	R. Heschl	(b)	R. Heschl	(b)
	L. Sup. Temporal	(b)	L. Sup. Temporal	(b)
	R. Sup. Temporal	(b)	R. Sup. Temporal	(b)
L‐language			L. Inf. Frontal Op.	(r)	L. Inf. Frontal Op.	(y)
R‐attention			R. Inf. Parietal	(r)
			R. Rolandic Op.	(r)
Other			L. Insula	(r)	L. Precuneus	(g)
			R. Insula	(r)	R. Precuneus	(g)
			L. Middle Cingulate	(r)	L. Cuneus	(g)
			L. Rolandic Op.	(r)	R. Cuneus	(g)

The AAL areas belonging to each cluster. Areas are grouped according to whether they belonged to visual, motor or auditory cortex, the left‐hemisphere language network or right‐hemisphere attention network. The letters in the brackets refer to the color used in Figure 5 to plot the corresponding area (cf., 3D rendering and signal plots; b: blue; g: green; r: red; y: yellow; SI: Supporting Information, Fig. S1).

The BOLD‐clustering of Cluster A identified three bModes (see Fig. 6A). The BOLD signal in bModes A1 and A2 was consistent with an event‐related response to short auditory stimuli. We evaluated the specificity of these results by averaging the predicted signal used for the GLM analyses, but time‐locking the predictors based on the time‐points associated with the bModes (cf., insets, in the plots of Fig. 6). If the time‐windows associated with the bModes are unrelated to the experimental paradigm, these plots would simply show some noise. By contrast, if the method correctly detects the occurrence of specific events in time‐series, these plots will show the corresponding effect coded in the stimulus/task predictors (i.e., design matrix of the GLM). This procedure confirmed that bModes A1–A2 corresponded to the isolated auditory events (light gray, in the insets), without any co‐occurring active task (continuous and dotted‐lines). As we already observed in the simulated data, the separation of A1 and A2 was due to a temporal shift by one data‐point/TR, which can be attributed to the timing of “when” the transient‐ISS became significant. By contrast, the single occurrence of bMode A3 appeared to be time‐locked just after the onset of a spatial‐task block (cf., dotted‐line in the inset of the corresponding plot) that we interpret as an incorrect classification.

The BOLD‐clustering of Cluster B (13 areas) identified four bModes, see Fig. 6B. For illustrative purposes, we plotted the signal of the four “auditory” regions in blue and the other nine regions in red. For bModes B1–3, the average signal within the windows showed a BOLD increase all the areas (see signal plots in Fig. 6B). By time‐locking the predicted responses (GLM) to the time‐windows of the bModes, we confirmed that these bModes corresponded to the auditory events co‐occurring with the onset of the spatial‐task blocks (see insets, in each plot). The brain regions belonging to these bModes were also consistent with this, and included auditory areas, as well as visual and motor areas, plus frontoparietal regions in the right hemisphere (see Table 3). Again the separation of bModes B1, B2, and B3 reflected the variability of when the transient‐ISS became significant. The single occurrence of bMode B4 corresponded to an onset of the spatial‐task without any associated auditory event (see inset plot; and also note the lack of any auditory response in bMode B4, cf., blue lines) and was interpreted as an incorrect classification.

These results indicate that the algorithm successfully detected the activation of auditory cortex for auditory events presented in isolation (bModes A1‐2), and segregated these from responses to auditory events that co‐occurred with the onset of the spatial‐task (bModes B1–3). Note that these findings were obtained without any a priori information about neither the areas involved, nor the type or timing of the stimuli presented to the participants.

The BOLD‐clustering of Cluster C (11 areas) dissociated block‐onset and block‐offsets for the language task. For illustration purposes, in Fig. 6C the signal of the visual, motor, and left premotor cortex is plotted in yellow, while the signal in the medial parietal cortex is plotted in green. The bMode C1 identified the onsets of the language blocks (cf., inset in the corresponding panel) and was associated with an increase of activation in visual‐motor areas and in the left inferior pre‐motor cortex (in yellow) and a concurrent decrease of the BOLD signal in the medial parietal cortex (in green). The bMode C2 identified the block‐offsets, with a signal decrease in task‐related areas and signal increase in the medial parietal cortex.

The BOLD‐clustering of Cluster D (four areas) revealed a more complex collection of bModes. Inspection of the signal associated with the different bModes and comparisons with the corresponding predicted signal (GLM) indicated that most of these bModes corresponded with the onset and offsets of the task‐blocks, including both the words and spatial conditions (see Supporting Information Fig. S1).

Figure 7A,B displays an overview of the results that enables visually comparing the bModes and the GLM results. All the bModes and the results of the omnibus F‐test (GLM, see also method section) are projected side‐by‐side on a 3D rendering of the brain. This highlights a good agreement, despite the fact that the F‐map was computed at the voxel‐level (approx. 65.000 voxels), while the projected bModes considered only the 116 AAL regions. This correspondence confirms that our fully data‐driven method is able to specifically identify brain areas associated with stimuli/tasks processing.

Main fMRI Experiment: Tensor ICA

Finally, we submitted our main fMRI dataset to tensor‐ICA, which also makes use of the consistency of the BOLD response across subjects to extract functional networks, in a fully data‐driven manner. The tensor‐ICA identified 60 ICs with reliable spatiotemporal components over subjects (P < 0.001, see Methods). The relationship between these 60 ICs and the experimental design was formally assessed by correlating the temporal mode of each IC with the four stimulus/task predictors of the GLM model. This revealed that 36 ICs showed some correlation with the at least one of the predictors (see Fig. 7C). Specifically, a total of 32 ICs correlated with the predictor associated with the words‐task, 21 with the spatial‐task, 16 with the isolated sound‐events, and 6 with the sound‐events coupled with the spatial‐task (see Supporting Information Figs. S2, for details). Accordingly, the tensor‐ICA was successful in extracting several components related to stimulus/task processing; for example, see IC‐60 in Supporting Information Fig. S3A, which could be linked specifically with the words‐task. Nonetheless, 24 ICs included temporal modes that did not correlate with any of the GLM predictors (cf., Fig. 7D). Moreover, in most cases the link between a single IC and the stimuli/tasks was not straightforward, even for components with a temporal mode highly correlated with the GLM. For example, the temporal mode of IC‐31 showed a positive correlation with the isolated‐sound GLM predictor (see Supporting Information Fig. S2), but the corresponding thresholded IC‐map did not include the auditory cortex (cf., Supporting Information Fig. S3B). In sum, tensor‐ICA successfully identified several stimulus/task‐related effects in our fMRI dataset, but the high number of ICs and the complexity of the temporal modes made it difficult to relate the output of the tensor‐ICA with the experimental stimuli and tasks, see also Discussion below.

DISCUSSION

We present a new method for the analysis of fMRI time‐series that was designed to identify transient network configurations that show stimulus or task‐related responses, and to do this without any a priori information (fully data‐driven approach). We applied the method both to simulated time‐series and to a real fMRI experiment with multiple tasks/conditions. The results of the simulations demonstrated that the method can effectively detect transient changes of network configurations. The systematic manipulation of sample size, SNR, size of the sliding‐windows and of the statistical threshold applied to the transient‐ISS analysis indicated that the method should perform well in standard fMRI settings. Indeed, the results of the main fMRI experiment confirmed that method can work in real neurophysiological conditions. We were able to track “when” (time‐windows) and “where” (networks) different stimuli and tasks triggered specific network configurations, without using any a priori information about the stimuli and the tasks. The results of a standard hypothesis‐based GLM analysis of the same dataset were in good agreement with our data‐driven results.

The four main features of the method presented here are: (1) it selectively targets areas engaged in stimulus/task processing; (2) it identifies inter‐regional effects; (3) it tracks the dynamic changes of these inter‐regional effects; and (4) it is fully data‐driven. These characteristics makes it ideal to investigate brain functioning in complex and naturalistic experimental setting (e.g., during movie watching), when there is little a priori knowledge about the underling processes and the construction of predictors for hypothesis‐driven analyses (e.g., GLM) is challenging. However, it should be noticed that the aim of the current work was to validate the method and, therefore, we designed a protocol that allowed us to have some control of the patterns of activity across brain regions and any change of these over time. This enabled us to quantify the method performance in our simulations and, for the main fMRI experiment, to compare the results of the new method with a standard GLM analysis [see also Pajula et al., 2012]. While representing only a simplified version of the network dynamics likely to occur with naturalistic stimuli, the current experimental protocol included a relatively complex set of inter‐regional effects. We generated inter‐regional dynamics that included areas switching network “membership,” as well as areas coupling/uncoupling over time [see Hutchison et al., 2013a and Smith et al., 2012]. The first constraint was implemented by choosing two tasks that would engage both “task‐specific” and “task‐common” areas (i.e., language‐task and spatial‐task) and by including “task‐specific” and “common” areas in the simulations (i.e., areas X and Y, vs. area Z; see Fig. 2B1). For the second constraint, we introduced the task‐irrelevant auditory stimuli that could either co‐occur with the spatial‐task, thus triggering coactivation of the auditory cortex with the spatial‐task network, or be presented in isolation (cf., also area W in the simulations, and see Fig. 2B1). The simulations and the results of the main fMRI experiment showed that the current method was able to identify and track these different types of dynamic network re‐configurations.

The main value of the current method is that it can identify transient network configurations that are relevant for stimulus and task processing, without making use of any a priori information. Most of the available methods used to study stimulus/task‐dependent changes of connectivity require knowing when such changes occur [e.g., see psychophysiological interactions, Friston et al., 1997; dynamic causal modeling, Friston, 1994; Granger causality, Goebel et al., 2003; structural equation modeling, McIntosh and Gonzalez‐Lima, 1994], which is typically not the case for experiments that make use of complex and naturalistic stimuli. Conversely, data‐driven methods enable identifying patterns of activation without any a priori knowledge, but the link between the patterns of activity/connectivity and the stimuli/tasks might be difficult to establish. Most commonly, the latter has been addressed by looking for correlations between the extracted indexes of connectivity and some behavioral measure ([e.g. Raz et al., 2012; Thompson et al., 2013; see also Krishnan et al., 2011] for a review of Partial Least Squares [PLS] methods, including behavior‐PLS and task‐PLS that can test for relationships between inter‐regional patterns of activity and behavior/task).

Here, we adopted a conceptually different approach, seeking to identify regions selectively engaged in stimulus/task processing using ISS [Hasson et al., 2004]. Moreover, we computed the ISS using sliding windows, thus obtaining time‐resolved information about these stimulus/task‐relevant effects. Recently, Glerean et al. [2012] also combined intersubject correlation analysis, and seed‐based inter‐regional correlations, with sliding windows. The authors were able to track dynamic changes of connectivity during movie watching and demonstrated that window sizes of 8–10 TR are optimal both for intersubjects synchronization and for the analysis of inter‐regional connectivity. This window‐size range is consistent with the results of our simulations (see Fig. 4B). However, because of the uncontrolled nature of the stimuli, Glerean's study could not validate the time‐resolved output against some “known” sequence of events and the associated patterns of brain activation. By contrast, here we were able to compare the output of our data‐driven method with simulated sequences of activations and coactivation patterns; as well as with the results of a standard GLM analysis, in our main fMRI experiment.

Another data‐driven approach that can make use of the consistency of the responses across subjects to extract information from fMRI time‐series, is tensor‐ICA. When applied to the current fMRI dataset, tensor‐ICA identified several components (n = 36) with temporal modes correlating with the stimulus/task design (see Fig. 7C and Supporting Information Fig. S2). However, despite these correlations, the single ICs could not be easily linked to the different stimuli/tasks presented to the subjects (e.g., see IC‐31, in Supporting Information Fig. S3B). Part of the issue may relate to the statistical thresholding of the IC maps. For instance, the IC‐31 was found to correlate with the auditory‐events, but did not show any effect in the auditory cortex at the default threshold (Z‐score = 3.3; corresponding to P > 0.5, in the alternative hypothesis testing framework of Melodic in FSL). However, lowering the threshold to Z‐score = 2.0 started reveling some voxels also in the superior temporal cortex. Furthermore, we should note that here we performed a high‐order decomposition (189 ICs were extracted initially, see Methods) and some stimuli/task‐related effects may have been “parcellised” across several ICs, some of which may have not survived the thesholding procedure.

One difference between the tensor‐ICA and the current method is that tensor‐ICA provides us with a continuous decomposition of the time‐series (i.e., each ICs has a single temporal mode with a duration corresponding to entire time‐series) while the bModes are associated with specific time‐points within the fMRI time‐series. The combination of tensor‐ICA and sliding windows may also enable capturing such dynamical aspects, but this would require in‐depth testing; for example, to determine the size of the sliding windows required to obtain reliable decompositions. Most importantly, a key feature of our current method is that it permits characterizing time‐varying and overlapping networks related to the different stimuli/tasks that share common processes (e.g., auditory responses in superior temporal areas of bModes A and B, associated with the auditory events in isolation vs. sounds occurring together with the spatial‐task; see Fig 6). Thus, the bModes can be used to detect and categorize events that occur in a repeated manner during the experiment, under the assumption that each occurrence of the event generates a specific network configuration.

The notion of using network configurations to identify and classify “brain‐states” in a data‐driven manner has been used in several recent studies. For example, Leonardi et al. [2014] used functional connectivity dynamics to classify three different cognitive states: counting backward versus episodic retrieval versus silent singing. However, their approach focused on distinguishing cognitive states across fMRI‐sessions [see also Shirer et al., 2012], rather than detecting state‐changes within a single session. Crossley et al. [Crossley et al., 2013] used graph analysis to examine a large dataset of task‐related PET and fMRI experiments, which revealed the recurrent coactivation of distinct areas across experiments/tasks. These results confirmed the possibility of using network configurations to track cognitive function in a data‐driven manner, but again without providing us any time‐resolved information for event‐categorization within a fMRI time‐series/experiment. Conversely, Liu et al. proposed to generate temporally resolved coactivation patterns (CAPs) based on the recurrent formation of specific spatial patterns of activity at the level of the single MR volume [Liu et al., 2013, Liu and Duyn, 2013]. Similar to the bModes, the CAPs are computed with k‐means clustering and each CAP can include both regions with a positive BOLD signal and regions with a negative signal. However, we should note that CAPs have been developed in the context of resting‐state data analysis and were designed to extract bursts of spontaneous of activity. By contrast, a key feature of the current algorithm is that the bModes include only areas that involved in stimulus/task processing. Indeed, the analysis of a resting‐state dataset with the current method did not detect any bMode (data not shown).

Because the identification of network configurations via coactivation does not rely directly on the BOLD signal, the areas belonging to the same network can have markedly different signal [cf., also CAPs, above; and “competitive interactions,” in Crossley et al., 2013]. Here, the characteristics of the BOLD signal within each area were examined further with the second clustering step (BOLD‐clustering). The analysis of the main fMRI experiment highlighted that the coupling between “task‐related” regions (sensory and motor areas) and the medial parietal cortex reversed completely at the beginning versus the end of a task‐block (see cluster C, in Fig. 6). Thus, these bModes identified couplings between areas with anti‐correlated BOLD signal, and separated time‐points when the coupling involved activation in a set of regions and concurrent de‐activation in another set of regions (i.e., at the beginning of a task‐block: bMode C1) versus the opposite pattern of activation/de‐activation (i.e., the end of a task‐block: bMode C2). These patterns were not included in our simulated dataset, but are consistent with the existing literature showing that regions in medial parietal cortex de‐activate when participants start a cognitive task and activate again as soon as the participant goes back to rest [Fox et al., 2005, McKiernan et al., 2003]. The finding of bModes reflecting these complex inter‐regional interactions emphasize the potential of the current method to highlight inter‐regional patterns of activation that, for example, might not be considered/tested in a classical hypothesis‐driven GLM analysis.

Future optimization of the method should consider improving the spatial resolution. The analyses presented here considered 116 regions, averaging the BOLD signal over many voxels (up to 1633, for the AAL right middle frontal region). This initial averaging reduces the possibility of detecting any effect that is highly localized in the brain. Alternative applications of the method may consider different templates, for example using cytoarchitectonic‐based atlases [Eickhoff et al., 2005] or defining ad hoc templates based on independent functional data. The latter could include the investigation of network re‐configurations within the occipital visual cortex, by first constructing a specific template based on individual retinotopic maps. A second point that could be optimized concerns the automatic detection of bModes that differ only because of small temporal shifts the BOLD signal (e.g., bModes B1–3, in Fig. 6). These should be detected and averaged together by shifting the time‐windows and the BOLD signal by the correct number of TRs. Finally, the main purpose/application of the method is investigate network dynamics involved in the processing of complex and naturalistic stimuli; that is, when the construction of predictors for GLM analyses is particularly challenging. Because each bMode includes a set of relevant time‐windows (when), future work will focus on re‐examining the stimuli at the relevant time‐points with the aim of associating each bMode with specific aspects of the stimuli [cf., Hasson et al., 2004; see also Crossley et al., 2013, who combined data‐driven network identification with database‐derived “functional labels”]. This will add an additional dimension to the output of the algorithm, which would then include a specific “functionality” for each bMode. Additional applications may involve using bModes to detect instances of cooperative versus noncooperative behavior during hyper‐scanning [Krill and Platek, 2012, Scholkmann et al., 2013], thus seeking to exploit ISS related to common high‐level goals, over and above the presence of a common sensory input as during movie‐watching.

In conclusion, we present and validate a fully data‐driven technique that detects inter‐regional configurations of stimulus/task‐related areas and tracks the dynamics of these networks over time. The method does not rely on the temporal relationship between the BOLD signal in the areas belonging to the same network but rather identifies network configurations based on the repeated engagement of the same set of task‐related regions over the fMRI time‐series. This new approach will contribute to the growing field of neuroimaging studies seeking to understand cognitive functions in naturalistic, real‐world‐like conditions [Peelen and Kastner, 2014].

Supporting Information

Click here for additional data file.