The goal of the described approach is to determine at what moments of the paradigm (temporal perspective), and between which regions (spatial perspective), significant reconfigurations in functional connectivity occur on functional magnetic resonance imaging recordings during which a time-locked stimulus is played.
Task-based functional magnetic resonance imaging bears great potential to understand how our brain reacts to various types of stimulation; however, this is often achieved without considering the dynamic facet of functional processing, and analytical outputs typically account for merged influences of task-driven effects and underlying spontaneous fluctuations of brain activity. Here, we introduce a novel methodological pipeline that can go beyond these limitations: the use of a sliding-window analytical scheme permits tracking of functional changes over time, and through cross-subject correlational measurements, the approach can isolate purely stimulus-related effects. Thanks to a rigorous thresholding process, significant changes in inter-subject functional correlation can be extracted and analyzed.
On a set of healthy subjects who underwent naturalistic audio-visual stimulation, we demonstrate the usefulness of the approach by tying the unraveled functional reconfigurations to particular cues of the movie. We show how, through our method, one can capture either a temporal profile of brain activity (the evolution of a given connection), or focus on a spatial snapshot at a key time point. We provide a publicly available version of the whole pipeline, and describe its use and the influence of its key parameters step by step.
Functional magnetic resonance imaging (fMRI) has become the tool of choice to non-invasively monitor the changes in brain activity resulting from external stimulation. More specifically, vivid interest has emerged about the understanding of statistical interdependence between regional activation time courses, known as functional connectivity (FC)1 and typically computed as Pearson's correlation coefficient. Functional interplays across the brain have been extensively shown to reconfigure as a function of the underlying task2,3,4.
Two analytical directions have separately been followed to go beyond this introductory characterization: on the one hand, the response induced in a given brain region by a time-locked stimulus was observed to strongly correlate across distinct subjects5. Quantifying this inter-subject correlation (ISC) showed potential to refine our understanding of cognition6,7,8,9 and brain disorders10,11. Further, this cross-subject correlational approach was also extended to the assessment of cross-regional synchronicity12, in what became known as the inter-subject functional correlation (ISFC) approach13.
On the other hand, the dynamic flavor of FC reconfigurations started to receive increased attention (see Hutchison et al.14, Preti, Bolton and Van De Ville15, Gonzales-Castillo and Bandettini16 for recent reviews on the resting-state and task-based sides of this question). In particular, whole-brain FC changes over time can be tracked through consecutive correlation measurements over a gradually shifted temporal sub-window17,18, revealing additional insight in the context of behavioral tasks19,20.
Here, we present a methodological framework that combines those two avenues. Indeed, we compute ISFC in sliding-window fashion to track the evolution of cross-regional synchronicity between the subjects exposed to a time-locked, naturalistic paradigm. Through the cross-subject aspect of the method, analyses are focused on stimulus-driven effects, while spontaneous fMRI changes (which are uncorrelated across subjects) are strongly damped. This is important because resting-state and task-evoked activity patterns are increasingly understood to be characterized by distinct properties21,22.
As for the dynamic component of the method, it enables a more complete and accurate characterization of task stimuli, particularly when probing a naturalistic paradigm in which a diverse set of cues (auditory, visual, social, etc.) are combined over time. Further, as the sound statistical evaluation of significant dynamic fluctuations has been hotly debated23,24, our approach takes particular care of this aspect of the analyses by isolating significant ISFC changes through comparison to appropriate null data.
We illustrate the method on a set of healthy subjects exposed to an audio-visual movie stimulus, for whom we show that the temporal and spatial ISFC change profiles arising from localized movie sub-intervals can be accurately extracted. In doing so, we also describe the influence of the main analytical parameters to be selected by the user. The presented findings are based on part of formerly published data25,26.
The following protocol has been approved by the local ethics committee (Biomedical Inserm 365 protocol C08-39).
1. Pre-imaging
2. Imaging
3. Data and Software Preparation
4. Data Preprocessing
Supplementary Figure 1: Example screenshot from the first preprocessing graphical user interface window. Voxel-wise time courses of gray matter voxels following the selected preprocessing options (top right plot), and covariates that may be used in the preprocessing (from top to bottom: cerebrospinal fluid/white matter average time courses, translational motion parameters, and rotational motion parameters. Please click here to view a larger version of this figure.
Supplementary Figure 2: Example screenshot from the second preprocessing graphical user interface window. Regional time courses following atlasing, before (top plot) and after (bottom plot) scrubbing and filtering according to selected parameters. Each curve depicts one regional time course randomly selected amongst all the available ones. Please click here to view a larger version of this figure.
5. Sliding-window ISFC Computations
Supplementary Figure 3: Example screenshot from the first ISFC-related graphical user interface window. (Top plot) Schematic representation of how often each considered session has its ISFC measurements computed (i.e., is not selected within the reference group). (Bottom plot) On an indicative subject, ISFC time courses computed for fifty example connections, selected as the ones exhibiting the largest summed absolute ISFC values across time. Please click here to view a larger version of this figure.
Supplementary Figure 4: Example screenshot from the second ISFC-related graphical user interface window. (Top left plot) On an indicative subject, ISFC time courses computed for three example connections, selected as the ones exhibiting the largest amount of significant ISFC excursions and displayed with their associated computed significance thresholds (horizontal lines). (Bottom left plot) For the same connections, associated excursion time courses averaged across subjects, with two-tailed 95% confidence intervals displayed as error measure. (Right plot) Spatial ISFC pattern (averaged ISFC excursions across subjects) for a selected time point indicated by a vertical black line on the ISFC and excursion plots. Positive ISFC excursions are shown in yellow, and negative ones in pink. The size and color code of the nodes are proportional to their degree. Please click here to view a larger version of this figure.
Here, we considered n = 15 typically developing (TD) subjects for whom we obtained written, informed consent. All were right-handed males (23.42 ± 7.8 years old). The chosen paradigm was an audio-visual scientific documentary for youngsters about the dangers of sun exposure. It contains a large array of visual, auditory and social stimuli, and can be watched at https://miplab.epfl.ch/index.php/miplife/research/supplement-asd-study.
We acquired two sessions per subject (RUN1 and RUN2) in which the assessed movie was displayed from 5 to 353 s (5.8 min duration). A resting-state segment also followed from 386 to 678 s (4.9 min duration). In addition, one solely resting-state session (RUN3) was acquired for each subject (excluding one who suffered from claustrophobia), lasting for 310 s (5.2 min). Example movie scenes and the timing of acquired data are summarized in Figure 1A. Importantly, the acquisition protocol was not optimal in the sense that resting-state recordings acquired just after movie exposure may be partly corrupted by spillover effects27; we make use of this data in the present findings to have a satisfying amount of samples for statistical thresholding, but this should be avoided whenever possible.
We excluded all sessions for which more than 10% of frames were scrubbed, at a threshold of 0.5 mm, and considered the parcellation from Craddock et al.34 (two-level temporal correlation algorithm) to generate regional time courses, for a total of 299 different brain regions.
ISFC was computed separately on (1) the movie-watching subparts of RUN1 and RUN2, (2) the resting-state subparts of RUN1 and RUN2, and (3) the resting-state RUN3 recordings. We used a window length W = 10 TR for the main presented results, and compare them to a lower value of W = 5 TR. Step size always remained equal to 1 TR. Bootstrapping was performed over 250 folds, including 6 session segments in each reference group.
Figure 1B displays ISFC time courses generated at W = 10 TR and W = 5 TR for three different representative connections: connection 1 involved a left inferior parietal region related to the expectation of moving objects (MNI coordinates: 41,9,32)35, and a right frontal opercular area linked to response inhibition (-34,-52,45)36. This latter region was also implicated in connections 2 and 3, respectively with an area implicated in sensory coordination (54,6,34)37, and one tied to the processing of the meaning of words (6,62,9)38.
A comparison across window lengths reveals that in the W = 5 TR setting, temporal variance in the subjects is overall larger in both the movie-watching and resting-state segment cases as compared to W = 10 TR, a known phenomenon in sliding-window analyses39. For connection 1, regardless of the window length, a localized subpart of the movie-watching recording (at around 55 s) shows a strong, synchronized ISFC increase across subjects, which largely exceeds the range of values taken in the resting-state case. Thus, we expect to capture this temporal subpart as a significant ISFC transient with our thresholding method.
For connection 2, we observe similar temporal dynamics, but for W = 5 TR, the increase becomes less easy to disentangle as compared to the resting-state time courses, due to the larger sliding-window methodology-related noise. As for connection 3, it reflects a case in which there is no clear response to the movie, and thus, the fluctuations from movie-watching and resting-state time courses are similar. The expected outcome at this analytical stage is a mix between connections that show clear stimulus-induced reconfigurations, and connections that do not respond.
Figure 1: Acquisition timing and example ISFC time courses. (A) The movie watched by the subjects involved a wide array of social situations (example images 1 and 4), scientific explanations with colorful panels (example images 2 and 5), and landscape sceneries (example image 3). Three sessions were acquired per subject: two (RUN1 and RUN2) included the movie stimulation (from 5 to 353 s, highlighted in green) followed by a resting-state period (from 386 to 678 s, shown in yellow), while one (RUN3) solely consisted in a resting-state recording (310 s duration, displayed in orange). (B) For three indicative connections (C1, C2 and C3, respectively dark green/red, light green/orange and turquoise/yellow traces), evolution of ISFC over time during movie-watching (cold colors) or resting-state (hot colors). For W = 10 TR (left panel), movie-watching ISFC changes more largely stand out as compared to W = 5 TR (right panel). Each trace reflects the ISFC time course of one session. This figure has partly been modified from Bolton et al.25. Please click here to view a larger version of this figure.
Figure 2A displays the results following statistical thresholding of ISFC time courses, for the same three connections as above. A time course value of 1 means that all subjects underwent the same ISFC increase at the same time point; a value of 0 means that no subject underwent a significant ISFC change; a value of -1 represents a synchronous ISFC decrease across all subjects. As before, we contrast W = 5 TR and W = 10 TR, and we also highlight two α-value cases: α = 0.01%, and α = 5%.
Fitting with the above observations, a lower window length reduces the amount of extracted significant ISFC changes. For connection 1, both W = 5 TR and W = 10 TR, however, extract the same particular moment (t = 55 s) as showing a strong ISFC increase. Taking a hemodynamic delay of roughly 5 s into account, this corresponds to a subpart of the movie when colored lines were extending towards a doll, and abruptly stopped just in front of it (46-49 s), fitting with the role of the involved regions in moving object expectation and response inhibition35,36.
When increasing α from 0.01% to 5%, one can observe a much lower specificity of the detected ISFC transients, likely including many false positives and expectedly showing much less temporal synchrony.
As another perspective that can be set on the data, Figure 2B shows the whole-brain spatial maps of significant ISFC changes at t = 55 s. It can be seen that the response to the movie scene extends far beyond the example connections described here.
Figure 2: Temporal and spatial snapshots of ISFC patterns. (A) ISFC transient time courses, averaged across subjects, for three indicative connections (C1, C2 and C3, respectively dark green, light green and turquoise traces). The movie scene that drove the ISFC changes is highlighted in light grey, and depicted by example images. For W = 10 TR (left column of plots), ISFC changes are more strongly detected than for W = 5 TR (right column of plots). For α = 0.01% (top row of plots), specificity to localized movie cues is larger than for α = 5% (bottom row of plots). Each trace reflects the ISFC transient time course of one session, and the two-tailed 95% confidence intervals are displayed as error measure. (B) For W = 10 TR and α = 0.01%, there is a neat, restricted spatial pattern of ISFC transients at t = 55 s (the peak ISFC transient value for C1); for W = 5 TR and α = 5%, connections undergoing a significant ISFC change at this time are much more numerous. Note that we assume a hemodynamic delay of around 5 s in the described temporality (i.e., a value of 55 s here relates to the movie stimulus at 50 s). This figure has partly been modified from Bolton et al.25. Please click here to view a larger version of this figure.
On a dataset of healthy subjects, we demonstrated how synchronous cross-subject increases and decreases in FC, the ISFC transients, would match temporally localized movie cues, providing information that goes beyond a static description. Although the use of cross-subject correlation measures enables to focus the analysis on stimulus-driven functional reconfigurations, one must also be aware that it limits the findings to the effects that are shared across the studied population: hence, low-level sensory processing is expected to be over-represented compared to frontal processing40. To bypass this limitation, new methods that also have the ability to extract the regions that most strongly vary across subjects are being developed41.
Another limitation from the introduced methodology arises from the sliding-window aspect, as temporal resolution of ISFC transient time courses is lowered compared to frame-wise approaches15. As we showed, a trade-off is needed between a sufficiently low window length to properly resolve dynamic ISFC reconfigurations, and a large enough size to obtain robust estimates. Two critical steps in our framework ensure that the extracted ISFC transients reflect truly occurring changes in connectivity: first, the high-pass filtering of regional time courses with the inverse of the window length32; second, the use of resting-state ISFC data for the generation of a relevant null distribution, with identical acquisition parameters as compared to the stimulus-related data. Of course, the latter also requires a lengthier global acquisition time, so that resting-state data can be gathered on top of stimulus-related sessions. As an alternative approach to avoid the additional resting-state recordings, we also offer the possibility to generate phase randomized data directly from the stimulus-related time courses, an approach often used in dynamic functional connectivity analyses23,24. Further evaluation on a subset of sessions revealed that although the resting-state null method is more conservative, and thus less prone to false positives, the global patterns of ISFC excursion detection were similar across both schemes (see Supplementary Figure 5).
Supplementary Figure 5: Detection of ISFC excursions across null data generation methods. For resting-state (left column, blue plots) or phase randomization (right column, red plots) null data generation methods, percentage of ISFC excursions extracted across connections. The bottom plots are an inset on the connections emanating from the first three considered brain regions. Error represents standard deviation across subjects. Please click here to view a larger version of this figure.
The duration of the resting-state acquisitions actually relates to a critical parameter of the analyses: the α-value. As exemplified above, a too lenient choice will lead to a large amount of false positives in the detected ISFC transients. The larger the amount of available resting-state data, the more stringent the achievable false positive rate, because thresholding can be based on more extreme values from the null distribution. As an indication, for n = 299 atlas regions as here and given our tally of 5,762 resting-state data points, we could at best achieve an α-value close to 0.01% (see step 5.2.6 for mathematical details).
Another key point pertaining to any fMRI analysis lies in the rigorous removal of possible motion-related artifacts from the analyzed data30,42. In particular, if one wishes to apply the introduced pipeline to diseased populations exhibiting marked motion in the scanner, we recommend that on top of including motion variables as covariates in the performed statistical analyses, additional preprocessing steps be run, such as wavelet denoising43 or ICA-AROMA44. Group comparison, for instance to compare ISFC transients between a healthy and a diseased group, can readily be performed by running the described approach in parallel on both groups of interest (see Bolton et al.25 for an example on a population diagnosed with autism spectrum disorders). However, a difference between the groups can then arise in two distinct settings: (1) an absent ISFC change in one group, or (2) a more heterogeneous evolution in that group. To disentangle those two factors, the pipeline should be run once more for the diseased group, using the healthy subject set as the reference group at the bootstrapping step. The former case would still result in an absent response, whereas the latter would not.
On top of what we described here, the introduced methodology also opens up promising future avenues: from an analytical side, ISFC transient maps could be viewed as brain graphs from which metrics quantifying brain connectivity could be derived45, or dynamic ISFC states could be extracted through clustering approaches and assessed in terms of their spatial and temporal characteristics17,46. In addition, one could also envisage the use of more sophisticated connectivity measurement tools than Pearson's correlation coefficient to reveal subtler sides of FC47,48.
From the experimental side, the application of our pipeline to a more extended set of paradigms is a promising perspective: for example, instead of a movie as studied here, one could envisage to use a piece of music49 or a narrative story13,50 as a time-locked stimulus. Alternatively, it could even be envisioned, through hyperscanning51, to probe naturalistic social communication52,53.
The authors have nothing to disclose.
This work was supported in part by each of the following: the Swiss National Science Foundation (grant number 205321_163376 to DVDV), the Bertarelli Foundation (to TB and DVDV), the Center for Biomedical Imaging (CIBM), and the National Agency for Research (tempofront grant number 04701 to ALG). The authors would like to thank Roberto Martuzzi and Giulia Preti for their contribution to the video content of this work as, respectively, the MRI operator and scanned volunteer.
Freesurfer version 6.0 | Laboratory for Computational Neuroimaging, Martinos Center for Biomedical Imaging, Boston (MA), USA | https://surfer.nmr.mgh.harvard.edu/fswiki/DownloadAndInstall | A MATLAB-compatible toolbox enabling to carry out various processing, visualisation and analytical steps on functional magnetic resonance imaging data |
MATLAB_R2017a | MathWorks | https://ch.mathworks.com/downloads/ | Working version of the MATLAB computational software (version 2014a or more recent should be used) |
Statistical Parametric Mapping version 12.0 (SPM12) | Wellcome Trust Center for Neuroimaging, University College London, London, UK | https://www.fil.ion.ucl.ac.uk/spm/software/spm12/ | A MATLAB-compatible toolbox enabling to perform statistical analyses on functional magnetic resonance imaging data |
Tim-Trio 3 T MRI scanner | Siemens | https://www.healthcare.siemens.ch/magnetic-resonance-imaging/for-installed-base-business-only-do-not-publish/magnetom-trio-tim | Magnetic resonance imaging scanner in which subjects have their functional brain activity recorded (at 3 T) |