Dynamic Inter-subject Functional Connectivity Reveals Moment-to-Moment Brain Network Configurations Driven by Continuous or Communication Paradigms

Summary

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.

Abstract

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.

Introduction

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)¹ 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 task^2,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 subjects⁵. Quantifying this inter-subject correlation (ISC) showed potential to refine our understanding of cognition^6,7,8,9 and brain disorders^10,11. Further, this cross-subject correlational approach was also extended to the assessment of cross-regional synchronicity¹², in what became known as the inter-subject functional correlation (ISFC) approach¹³.

On the other hand, the dynamic flavor of FC reconfigurations started to receive increased attention (see Hutchison et al.¹⁴, Preti, Bolton and Van De Ville¹⁵, Gonzales-Castillo and Bandettini¹⁶ 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-window^17,18, revealing additional insight in the context of behavioral tasks^19,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 properties^21,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 debated^23,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 data^25,26.

Protocol

The following protocol has been approved by the local ethics committee (Biomedical Inserm 365 protocol C08-39).

1. Pre-imaging

1. Enroll a study population of subjects, obtaining written, informed consent for all of them. Seek approval from the local ethics committee.
2. Select a paradigm to investigate that can be applied to all subjects in time-locked manner.
   NOTE: Here, we used an audio-visual scientific documentary for youngsters (https://miplab.epfl.ch/index.php/miplife/research/supplement-asd-study).

2. Imaging

1. For each subject to consider in the analyses, perform at least one functional imaging session in which the scanned volunteer is subjected to the time-locked paradigm of interest.
   1. Use a 3 Tesla MRI scanner to acquire transverse slices through an echoplanar imaging sequence.
   2. Employ the following imaging parameters: voxel size = 3 mm x 3 mm x 3 mm, repetition time (TR) = 2 s, echo time = 50 ms, field of view = 192, 40 slices.
      NOTE: Faster TR values are encouraged within the scope of feasibility. The protocol can also be applied with a more restricted field of view (e.g., for analyses restricted to a specific brain sub-structure), which would enable either a better temporal resolution (lower TR), or a spatially more precise analysis.
   3. Leave a few seconds of recording (≥ 2 TR) before and after the presentation of the stimulus.
2. Perform at least one separate functional imaging session in which the scanned volunteer lies at rest in the scanner, eyes closed and instructed not to fall asleep.
   NOTE: Separate stimulus-related and resting-state acquisitions prevent otherwise possible interplays between the conditions (e.g., having watched the movie beforehand may leave a lasting trace to a subsequently acquired resting-state recording)²⁷. If it is not desired to go through the aforementioned additional resting-state acquisition, an alternative (albeit more prone to the detection of false positives; see Discussion) computational option in the pipeline replaces this data by surrogate time courses computed from the paradigm-related signals (see step 5.1.2).
3. Perform structural imaging.
   1. Use a 3 Tesla MRI scanner and a T1-weighted magnetization-prepared rapid acquisition gradient echo sequence.
   2. Employ the following imaging parameters: voxel size = 1 mm x 1 mm x 1 mm, field of view = 256, 176 slices.

3. Data and Software Preparation

1. For each session to analyze, ensure the existence of the following data files:
   1. A set of functional MRI volumes, present as separate 3D NIFTI or HDR/IMG files, with a consistent numbering scheme (e.g., "fMRI_0001", "fMRI_0002", etc.).
   2. A T1 structural MRI image, in NIFTI or HDR/IMG format.
   3. An atlas of interest in Montreal Neurological Institute (MNI) space, in NIFTI format.
      NOTE: An example of required input files is provided for a representative subject ("S17"), along with the full pipeline code, at https://c4science.ch/source/Intersubj_pipeline.git 
2. Download the latest version of the publically available Freesurfer software²⁸ (https://surfer.nmr.mgh.harvard.edu/fswiki/DownloadAndInstall).
3. Download the latest version of the publically available Statistical Parametric Mapping (SPM) MATLAB toolbox from https://www.fil.ion.ucl.ac.uk/spm/software/spm12/.
4. Open MATLAB (version 2017a or more recent) and locate the newly downloaded "freesurfer" and "spm12" folders. For each, right click on it and select the Add to Path > Selected Folders and Subfolders option.

4. Data Preprocessing

1. In the MATLAB terminal, type spm to launch the SPM12 main menu, and click on the fMRI button to access the preprocessing options devoted to fMRI data. Perform the following steps separately for each fMRI session to preprocess.
   1. Click on Realign (Est & Res), and in the newly open Batch Editor window, double click on Data > Session. In the newly open Session window, select all functional images to process. Then, click on the Done button, and afterwards, on the Run Batch icon from the Batch Editor window (green triangle). Wait until the realignment step finishes, as indicated in the MATLAB terminal window.
   2. Click on Coregister (Est & Res), and in the newly open Batch Editor window, double click on Reference image. In the newly open Reference Image window, select the average functional volume created in the following step, prefixed with "mean", and click on the Done button. Then, double click on Source image, and in the newly open Source Image window, select the T1 image. Click on the Done button, and afterwards, on the Run Batch icon from the Batch Editor window (green triangle). Wait until the coregistration step finishes, as indicated in the MATLAB terminal window.
      NOTE: The T1 image is overwritten at this step, so that the updated one lies in the same space as the functional volumes.
   3. Click on Segment, and in the newly open Batch Editor window, double click on Volumes. In the newly open Volumes window, select the T1 image, and click on the Done button. Then, in the Batch Editor window, double click on Deformation Fields and select the Inverse option. Click on the Run Batch icon (green triangle), and wait until the segmentation step finishes, as indicated in the MATLAB terminal window.
2. Type JOVE_GUI1 in the MATLAB terminal to open the first preprocessing graphical user interface window. Perform the following steps for each fMRI session to analyze.
   1. Click on Enter fMRI data, and select all realigned functional volumes created in step 4.1.1 (prefixed with "r"). For IMG/HDR files, select both the IMG and HDR volumes.
   2. Enter the TR of the data (in seconds) in the dedicated editable text window.
   3. Click on Enter T1 data and select the three probabilistic tissue type volumes created in step 4.1.3 (prefixed with "c1", "c2" and "c3").
   4. Click on Enter motion file, and select the text file containing motion parameters from the session of interest, created in step 4.1.1 and prefixed with "rp".
   5. Select the desired type of preprocessing, that is, whether the data should be detrended or not (respectively setting the dedicated radio button on or off), and what covariates should be regressed out (by selecting the appropriate option from the dedicated list).
      NOTE: The regression step is inspired from a function originating from the DPARSF toolbox²⁹. The white matter and cerebrospinal fluid signals from individual subjects are averaged over the voxels for which the respective template DPARSF probabilistic tissue map showed a signal larger than 0.99. In our analyses, we detrended the data, and regressed out white matter/cerebrospinal fluid time courses as well as constant, linear and quadratic trends.
   6. To preprocess the data, click on Preprocess, and wait for the display to appear in the window. The data can be re-preprocessed differently by modifying the options, and clicking again on the Preprocess button.
      NOTE: The gray matter plot is inspired from the representation suggested by Power et al.³⁰.
   7. To save the output for following steps, click on the Save button. To clear the content of the window, click on the Clear button.

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.

3. Type JOVE_GUI2 in the MATLAB terminal to open the second preprocessing graphical user interface window. Perform the following steps for each fMRI session to analyze.
   1. Click on Select data, and select the data file saved in step 4.2.7 (named "ISFC_VX.mat").
   2. Click on Select motion, and select the text file containing motion parameters from the session of interest, created in step 4.1.1 and prefixed with "rp".
   3. Click on Select atlas, and select the NIFTI file representing the atlas to use for parcellation.
   4. Click on Select inverse warp, and select the NIFTI file representing the deformation field from MNI to native space, created in step 4.1.3 and prefixed with "iy".
   5. Click on Select fMRI volume, and select any of the fMRI data volumes.
      NOTE: This step enables to access the header information of the functional data, hence why the actual chosen volume is not important.
   6. Enter the TR of the data (in seconds) in the dedicated editable text window.
   7. Enter scrubbing-related information: the type of scrubbing to perform (i.e., how many frames to scrub out before and after the tagged ones) in the "Scrubbing type" list, and the framewise displacement threshold value (Power's criterion³¹) above which an fMRI volume should be scrubbed in the "Scrubbing threshold" editable text window (in mm).
      NOTE: Cubic spline interpolation is performed on the scrubbed data points to replace them with estimated values from neighboring samples. In our analyses, we scrubbed one frame after the tagged volumes, and used a 0.5 mm threshold for scrubbing.
   8. Enter the size of the sliding-window W to use for ISFC computations (see step 5), in TRs.
      NOTE: This piece of information will enable filtering of the time courses through a function originating from the DPARSF toolbox²⁹, at f = 1/W Hz³². In our analyses, we used W = 10 TR as a trade-off value to capture dynamic fluctuations while conserving enough samples for robust estimations.
   9. Click on the Plot button to display indicative atlased time courses before (top plot) and after (bottom plot) the scrubbing and filtering steps. Verify, by visual inspection, that following the selected preprocessing steps, those output signals do not incorporate salient artifactual components.
   10. To save the outputs for following steps, enter a save name in the dedicated editable text window, and click on the Save button. To clear the content of the window, click on the Clear button.

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

1. Type JOVE_GUI3 in the MATLAB terminal to open the first ISFC-related graphical user interface window. Perform the following steps separately for each type of acquired fMRI session segment (stimulus-related segments, resting-state segments of stimulus-related sessions, and purely resting-state segments).
   1. Click on Load data, and select all the appropriate data files created through step 4.3.
   2. Select whether the selected session segments should undergo phase randomization.
      NOTE: Phase randomization can be used as an alternative option for the generation of null data from stimulus-related signals, if no resting-state recordings are available.
   3. Enter the TR of the data (in seconds) in the dedicated editable text window.
   4. Enter sliding-window parameters to use for the analysis in the dedicated editable text windows: window size (in TRs) over which connectivity measurements should be computed, and step size (in TRs) by which successive windows should be shifted.
      NOTE: In our analyses, we used a window size of 10 TR and a step size of 1 TR.
   5. Modify the "Session types" table to specify which of the loaded session segments were acquired upon the same experimental condition. Use increasing integer numbers from 1 onwards to tag different types of segments (e.g., if the stimulus was displayed for the first or for the second time in a given recording). Leave the table untouched if only one type of session segment was acquired.
      NOTE: A session in the present work may refer either to a combined movie/resting-state recording (termed RUN₁ and RUN₂ in Figure 1A), or to a purely resting-state recording (RUN₃). A session segment refers to a sub-portion of a session recording, either when the movie was watched, or when the subjects lied at rest. The above information is used in the subsequently described ISFC computations (see step 5.1.8) to limit the confounding influence of different session segment types.
   6. Enter bootstrapping-related parameters in the dedicated editable text windows: the number of bootstrapping folds over which to perform ISFC computations, and the number of subjects that should constitute the reference group for each fold of ISFC computations.
      NOTE: In our analyses, we used 250 bootstrapping folds, and 6 subjects into the reference group.
   7. Enter specifications about which sub-portion of time courses should be analyzed in the Timing parameters section, in the dedicated editable text windows. A start index and an end index (in TRs) should be provided. To analyze the whole recording duration, use 1 as start index and the number of samples as end index.
   8. Click on the Plot button to perform ISFC computations. Displays are gradually updated over time, along with the amount of elapsed bootstrapping folds. For a region pair (i,j) and a sliding window index τ, ISFC is computed as the average of cross-correlations between session segment s and all session segments from the reference group, within a sliding-window of length W; denote this reference group by Ψ, its number of subjects by N_Ψ, and let x_i^[s](t) the time course of region i for session segment s at time t; an ISFC estimate is then given by:
      
      ISFC measurements are computed over the specified amount of bootstrapping folds, and with the selected number of session segments used as a reference group at each fold (see step 5.1.6). If several session segment subtypes are included, a mixture of subtype samples always composes the reference group. The final output for each session segment is the average ISFC across all folds in which it was not included as a reference measurement.
      ​NOTE: The reference group is the set of session segments to which the functional time courses of session segment s are compared at each fold of the bootstrapping process. For the results to be more robust to outlier data points, ISFC is computed multiple times on a different reference group (that is, a different subset of session segments). Importantly, the acquisition time t does not match the sliding window index τ, as the latter is computed over a set of W data points, and relies on the window step size for successive estimates. The bootstrapping process was inspired from a former study by Byrge et al.³³.
   9. To save the outputs for following steps, enter a save name in the dedicated editable text window, and click on the Save button. To clear the content of the window, click on the Clear button.

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.

2. Type JOVE_GUI4 in the MATLAB terminal to open the second ISFC-related graphical user interface window.
   1. Click on Load ISFC data and select the stimulus-related ISFC output file(s) created in step 5.1.
   2. Click on Load null data and select, depending on the used null data generation scheme, either the resting-state ISFC, or the phase randomized stimulus-related ISFC output file(s) created in step 5.1.
   3. Click on Load codebook and select the codebook file created in step 4.3.
   4. Enter the TR of the data (in seconds) in the dedicated editable text window.
   5. Enter the sliding-window parameters used in the computations of step 5.1 (window size and step size, in TRs) in the dedicated editable text windows.
   6. Enter (in percent) the α-value at which the ISFC time courses should be thresholded to highlight significant changes in the dedicated editable text window.
      NOTE: Here and elsewhere, when referring to an α-value of 2.5%, it means that significance is achieved when a value is lower than the 2.5^th percentile, or larger than the 97.5^th percentile, of the null data. In our analyses, we had 5,762 resting-state data points to our disposal, and selected an α-value of 10^-4. This means that we wanted 0.01% of data samples to be larger or equal to the selected thresholds past which an ISFC excursion would be deemed significant. For comparison purposes, the α-level demanded by Bonferroni correction would be 0.05/44,551 = 1.12 x 10^-6, and the most stringent possible α-level enabled with our amount of data (n samples) would be .
   7. Click on the Plot button to perform the ISFC thresholding process, in which all available null ISFC measurements are aggregated, for a given connection, to construct a null distribution, following which stimulus-related ISFC measurements are thresholded according to the selected α-value. Time points at which a stimulus-related ISFC value statistically significantly exceeds the null distribution are tagged as -1/+1 for significant ISFC decreases and increases, respectively.
      NOTE: The thresholding process draws inspiration from the resting-state dynamic FC work of Betzel et al.²³.
   8. To visualize the ISFC spatial patterns at different time points, drag the slider below the ISFC excursion plot.

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.

Discussion

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 processing⁴⁰. To bypass this limitation, new methods that also have the ability to extract the regions that most strongly vary across subjects are being developed⁴¹.

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 approaches¹⁵. 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 length³²; 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 analyses^23,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 data^30,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 denoising⁴³ or ICA-AROMA⁴⁴. 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.²⁵ 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 derived⁴⁵, or dynamic ISFC states could be extracted through clustering approaches and assessed in terms of their spatial and temporal characteristics^17,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 FC^47,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 music⁴⁹ or a narrative story^13,50 as a time-locked stimulus. Alternatively, it could even be envisioned, through hyperscanning⁵¹, to probe naturalistic social communication^52,53.

Materials

Name	Company	Catalog Number	Comments
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)