Identifying Resting-State Networks from fMRI Data Using ICAs
A hands-on tutorial: from data collection to feature extraction
This tutorial offers an insight into reasonable practices of neuroimaging research via a basic hands-on tutorial for analyzing fMRI data. I discuss (1) posing a suitable research question, (2) data collection (3) quality control, and preprocessing (4) feature extraction.
The next post- Classifying ADHD from healthy controls using LSTMs with rs-fMRI data focuses on data analysis.
The code for this tutorial is available here. To be prepared, please download freesurfer, docker, and python.
I chose to examine resting-state networks among schizophrenia (SZ) patients and compare them to these healthy people.
I am curious whether schizophrenia displays functional differences, and mainly whether such differences could be found in the resting-state network.
I believe this is an interesting question (and a good case-study for a tutorial about fMRI analysis) since it offers an insight into functional connectivity analysis via the exploration of a disorder known to be associated with abnormal neural circuitry (Yu et al., 2012).
Data collection
No, this section does not concern with the operation of MRI machines, but rather it focuses on various online resources offering existing scans for public use. Among them, openfMRI, Oasis, openNeuro, and 1000 Functional Connectomes project.
Most datasets present raw, de-identified, data — while some include some preprocessing scripts. A useful dataset will be organized following some standard; a common one is the Brain Imaging Data Structure (BIDS) format. In a BIDS formatted dataset, you can expect to see the files organized in a hierarchy: first the subjects id, then the type of scan, and later (if applicable) different tasks or runs relevant for every kind of scan. The files themselves are most likely to be NIFTI files (Neuroimaging Informatics Technology Initiative) or compressed NIFTI files (.nii or .nii.gz).
I use this dataset from openneuro. It comprises scans of 290 individuals diagnosed with various physio-psychiatric disorders (schizophrenia, bipolar, ADHD, depressive disorder, PTSD, and more) as well as healthy controls. Additional to the variety of conditions, this dataset includes several types of scans capturing different tasks. To explore the research question posed above, I chose to focus on the resting-state fMRI scans of SZ patients.
fMRI data were acquired 3T Siemens Trio scanners. Scan parameters are as follows- slice thickness=4 mm, 34 slices, TR=2 s, TE=30 ms, flip angle=90°, matrix 64×64, FOV=192 mm, oblique slice orientation. During resting-state scans, participants were asked to remain relaxed and keep their eyes open; they have not been presented any stimuli or asked to respond during the scan (Poldrack, 2016A).
From the dataset description, we can see that the T1w (anatomical) scans have been defaced (anonymized via facial features removal) using freesurfer mri_deface. Defaced data, by definition, underwent some preprocessing before. Ideally, working with raw data ensures more uniform preprocessing, but if that is not the case, it is essential to know which steps have been taken already.
Since this post focuses on data preparation, I demonstrate the tutorial using a subset of the data. One SZ (sub-50006) and one control (sub-10228) subjects.The data directory (organized according to BIDS) should look like that:
directories, this includes the dataset description, a table with subject demographic information, and a task description file.
Let’s have a look at our data.
NIFTI images contain a header that can nicely provide you with a lot of the metadata about the image. This information will become handy when we need to correct the image for noise and other distortions. (more details here)
We can also find the number of slices (dim of the Z-axis) and the TR for the scan via the metadata. TR: 2.0 and the # of slices: 34. (hooray, these numbers seem to match those reported by the authors).
Additionally, we can easily visualize them using freesurfer. Inspecting the image allows us to identify obvious problems with the scan (weird orientation, artifacts, poor contrast, or even the lack of image all together). There are tools allowing to perform such quality control automatically (we will discuss them soon).
Quality Control and Preprocessing
Before running the preprocessing pipeline, let’s have a look at the quality of the data in hand. Quality control (QC) helps detect images that suffer severe problems- such as artifacts or heavy distortion.
QC for data used to analyze large-scale networks focuses particularly on temporal consistency. Power et al. (2011) detail several image quality metrics (IQM) focused on identifying non-spurious correlations. I use a software named mriQC, since it presents IQMs following peer-reviewed publications (among them power et al., 2011; Krüger et al., 2011). The list below describes the IQMs relevant to ensuring temporal consistency and avoiding spurious correlations. Images containing bad artifacts or very low SNRs would be excluded while the rest enter the preprocessing pipeline. DVARS and FD will be used during preprocessing to form a temporal mask (i.e. highlight specific slices to be removed).
IQMs
SNR (Signal-to-Noise Ratio)
- It is the primary spatial metric. Uses the structural data to extract the ratio of signal power to noise power. SNR is important to estimate the statistical power.
- It calculates the signal amplitudes over the RMS of the noise. A higher value indicates a more stable signal over repeated measurements, thus higher reliability of the measure.
tSNR (temporal Signal-to-Noise Ratio)
- Reports the median of a time-course SNR for a time series.
- It calculates the ratio between the average BOLD signal (across time), and its corresponding temporal standard-deviation map.
- tSNR is important when using the data later to compare signal changes (low tSNR may bias the data)
DVARS (Derivative of the VARiance)
- Measures of how much the intensity of a brain image changes in comparison to the previous time point. Thus, it is a temporal metric.
- It indexes the rate of change of BOLD via the derivative RMS variance over voxels. Slices with abnormally high values are marked for exclusion.
FD (Framewise Displacement)
- Measures how much the head changed position between sequential frames.
- It calculates the sum of absolute values of the derivatives of the realignment params. Slices with abnormally high values are marked for exclusion.
This command runs the QC for all participants in your data directory- this is a long and heavy calculation. Please make sure to allocate sufficient resources (or run in one by one).
Before looking at the metrics, I take a look at the structural visual report. The key takeaway from the report is the defacing process (seen here as a yellow mark over the BOLD noise visualization).
The functional report reports DVARS, FD, and an outliers measurement over time. These graphs indicate that sub-10228 provides a much cleaner scan than sub-50006.
The group functional report reveals high values for both tSNR and SNR.
After inspecting the data and making sure we are happy with its quality, let us perform preprocessing.
I use fMRIPrep, since it attempts to provide a standardized preprocessing pipeline and to address the challenge of robust and reproducible preprocessing for fMRI data. fMRIPrep is an analysis-agnostic tool that offers a flexible, minimal, and robust pipeline for fMRI Data. Additionally, it is easy to use, open-source software that keeps a transparent discussion about their limitations and development.
Pay attention to the flags I used- ‘- -fs-license-file’ indicates to fMRIPrep where is your freesurfer license file, ‘- -fs-no-reconall’ skips freesurfer reconstruction step (it is somewhat buggy and not essential here), and ‘- -output-spaces’ defines the anatomical spaces from which to resample the data (T1w refers back to the image itself, fsaverage5 refers to freesurfers template and template refers to the MNI map).
Similarly, this command runs fMRIPrep for all participants in your data directory. Be aware. fMRIPrep is heavy. More info here.
fMRIPrep uses both structural and functional images in its pipeline. It outputs the images post head-motion correction, slice-timing correction, and aligned into the same-subject’s structural space and in accordance with the Montreal Neurological Institute (MNI) space. MNI is a standardized space allowing normalize different brains into a single, comparable map.
fMRIPrep pipeline in a nutshell. For more info please visit the documentation.
- T1w images are analyzed to detect brain mask and brain tissue segmentations
- The T1w Reference scans are used for spatial normalization. Images are being mapped onto their native space as well as MNI (and any others you have specified when via the flag ‘-output-spaces’)
- The first functional realignment targets head motion by using FSL’s mcflirt algorithm to establish motion correction.
- Slice time correction using AFNI’s 3dTShift function to realign in time all slices to the middle of each TR.
- The Susceptibility Distortion Correction (SDC) focused on fixing spatial distortion caused by the inhomogeneity of the field inside the scanner.
- Frames that exceeded a threshold of 0.5 mm FD or 1.5 standardized DVARS were annotated as motion outliers.
The order of these steps supports a more robust correction (Power et al., 2012).
Compare the DVARS, FD, and outliers measurements of the raw data and the preprocessed one, we see that the preprocessing reduced the noise in the data.
Feature extraction: Resting-state analysis methods
Once the data is clean, denoised, and preprocessed- we can (finally) identify connectivity features. Remember, resting-state networks exhibit correlated temporal patterns across spatially independent regions (Jandric et al., 2018). Two common ways to examine such correlation patterns are via seed-based correlation analysis (SCA) and independent component analysis (ICA). SCA (Biswal et al., 1995) is a model-based method; thus, it requires an a-priori division of the brain into regions of interests (ROIs) (Lee et al., 2013). SCA is especially useful since it is effortless to implement and straightforward to interpret (Soares et al., 2016). Basically, it uses the time series of an a-priori selected seed region to identify whole-brain functional connectivity maps (based on the specified ROIs). ICA, on the other hand, is a model-free method and decomposes data from the whole brain into the time courses and spatial maps of recognized signals. These signals, the independent components, are a collection of regions that exhibits maximal spatial independence but co-varying temporally (Cole et al., 2012; Smith et al., 2013).
Nilearn is a python library that supports both these (and many other) neuroimaging analysis methods.
First, we need to upload the preprocessed images:
These tools are the ones we need in order to get an atlas of regions and produce a correlation matrix. As discussed above, in SCA- we use an a-priori atlas and in ICA- we create one based on the dataset.
The first function creates a masker. I like the Nilearn definition of a masker. They as us to “think of masker objects as swiss-army knives for shaping the raw neuroimaging data in 3D space into the units of observation relevant for the research questions at hand”. The masker will filter in only the parts of the data we are interested in. The masker uses the atlas to transform the fMRI scan.
Actually this is it. A transformed fMRI scan reflects a 2D matrix of regions to time steps (i.e. a [10,160] would reflect the values of ten regions across 160 time-steps).
But to have a nice, simple, visualization- we calculate the correlation matrix of the extracted regions and plot it (second and third functions).
SCA Analysis
I include two types of atlases. The first is the Harvard-Oxford probabilistic atlas (Makris et al., 2006; Frazier et al., 2005; Desikan et al., 2006; Goldstein et al., 2007) and the second is Smith’s ICA functional map (Smith et al., 2009).
The Harvard-Oxford Atlas is based on MNI defined regions thus demonstrates a well known a-priori segmentation of the brain to regions. On the other hand, the Smith’s atlas reflects regions resulting from an independent analysis of resting-state fMRIs and their demonstrated activation brain dynamics. Thus it is somewhat in between an SCA to an ICA. It does present a mask independent of the current data, but its origins are in an ICA of a similar dataset.
We use both these atlases to create maskers based on their extracted regions. Calculate their corresponding correlation matrixes and plot.
ICA Analysis
In ICA, we are not importing atlases but calculating them using the dataset itself.
Once the regions are generated, the process is as we did with the SCA. We create a mask and calculate the correlations.
In the next part of this tutorial — How do we use these data in statistical analysis. Are there differences in the resting-state networks between SZ subjects to the controls?
As a little teaser, this figure shows the correlation matrix of one SZ patient. Can you spot any obvious differences?
References
Biswal, B., Yetkin, F. Z., Haughton, V. M., and Hyde, J. S. (1995). Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn. Reson. Med. 34, 537–541. doi: 10.1002/mrm.1910340409
Cole, D. M., Smith, S. M., & Beckmann, C. F. (2010). Advances and pitfalls in the analysis and interpretation of resting-state FMRI data. Frontiers in systems neuroscience, 4, 8.
Desikan RS, Ségonne F, Fischl B, Quinn BT, Dickerson BC, Blacker D, Buckner RL, Dale AM, Maguire RP, Hyman BT, Albert MS, Killiany RJ. An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. Neuroimage. 2006 Jul 1;31(3):968–80.
Esteban O, Birman D, Schaer M, Koyejo OO, Poldrack RA, Gorgolewski KJ; MRIQC: Advancing the Automatic Prediction of Image Quality in MRI from Unseen Sites; PLOS ONE 12(9):e0184661; doi:10.1371/journal.pone.0184661.
Frazier JA, Chiu S, Breeze JL, Makris N, Lange N, Kennedy DN, Herbert MR, Bent EK, Koneru VK, Dieterich ME, Hodge SM, Rauch SL, Grant PE, Cohen BM, Seidman LJ, Caviness VS, Biederman J. Structural brain magnetic resonance imaging of limbic and thalamic volumes in pediatric bipolar disorder. Am J Psychiatry. 2005 Jul;162(7):1256–65
Jandric, D., Lipp, I., & Mumford, J. (2018). OHBM OnDemand How-to: Resting State fMRI Analysis. Retrieved 26 October 2019, from https://www.ohbmbrainmappingblog.com/blog/ohbm-ondemand-how-to-resting-state-fmri-analysis
Krüger et al., Physiological noise in oxygenation-sensitive magnetic resonance imaging, Magn. Reson. Med. 46(4):631–637, 2001. doi:10.1002/mrm.1240.
Lang, E. W., Tomé, A. M., Keck, I. R., Górriz-Sáez, J. M., & Puntonet, C. G. (2012). Brain connectivity analysis: a short survey. Computational intelligence and neuroscience, 2012, 8.
Lee, M. H., Smyser, C. D., and Shimony, J. S. (2013). Resting-state fMRI: a review of methods and clinical applications. AJNR Am. J. Neuroradiol. 34, 1866–1872. doi: 10.3174/ajnr.A3263
Lindquist, M. A. (2008). The statistical analysis of fMRI data. Statistical science, 23(4), 439–464.
Logothetis, N. K. (2008). What we can do and what we cannot do with fMRI. Nature, 453(7197), 869.
Lyon, L. (2017). Dead salmon and voodoo correlations: should we be sceptical about functional MRI?. Brain, 140(8), e53-e53.
Maclaren J, Herbst M, Speck O, Zaitsev M. Prospective motion correction in brain imaging: a review. Magn Reson Med 2013; 69: 621–636.
Makris N, Goldstein JM, Kennedy D, Hodge SM, Caviness VS, Faraone SV, Tsuang MT, Seidman LJ. Decreased volume of left and total anterior insular lobule in schizophrenia. Schizophr Res. 2006 Apr;83(2–3):155–71
Mueser, K. T., and Jeste, D. V. (Eds.). (2011). Clinical handbook of schizophrenia. Guilford Press
Poldrack, R. A., Congdon, E., Triplett, W., Gorgolewski, K. J., Karlsgodt, K. H., Mumford, J. A., … & Bilder, R. M. (2016A). A phenome-wide examination of neural and cognitive function. Scientific data, 3, 160110.
Poldrack, R. A., Baker, C. I., Durnez, J., Gorgolewski, K. J., Matthews, P. M., Munafo, M., … & Yarkoni, T. (2016B). Scanning the Horizon: Future challenges for neuroimaging research. bioRxiv, 059188.
Power, J. D., Cohen, A. L., Nelson, S. M., Wig, G. S., Barnes, K. A., Church, J. A., … & Petersen, S. E. (2011). Functional network organization of the human brain. Neuron, 72(4), 665–678.
Power JD, Plitt M, Kundu P, Bandettini PA, Martin A (2017) Temporal interpolation alters motion in fMRI scans: Magnitudes and consequences for artifact detection. PLOS ONE 12(9): e0182939. doi:10.1371/journal.pone.0182939.
Raichle, M. E., MacLeod, A. M., Snyder, A. Z., Powers, W. J., Gusnard, D. A., & Shulman, G. L. (2001). A default mode of brain function. Proceedings of the National Academy of Sciences, 98(2), 676–682.
Smith SM, Fox PT, Miller KL, Glahn DC, Fox PM, Mackay CE,
Filippini N, Watkins KE, Toro R, Laird AR, and Beckmann CF. 2009.
Correspondence of the brain’s functional architecture during activation and rest. Proc Natl Acad Sci USA (PNAS), 106(31):13040–13045.
Smith, S. M., Vidaurre, D., Beckmann, C. F., Glasser, M. F., Jenkinson, M., Miller, K. L., … & Barch, D. M. (2013). Functional connectomics from resting-state fMRI. Trends in cognitive sciences, 17(12), 666–682.
Soares, J. M., Magalhães, R., Moreira, P. S., Sousa, A., Ganz, E., Sampaio, A., … & Sousa, N. (2016). A hitchhiker’s guide to functional magnetic resonance imaging. Frontiers in neuroscience, 10, 515.
Yu, Q., A Allen, E., Sui, J., R Arbabshirani, M., Pearlson, G., & D Calhoun, V. (2012). Brain connectivity networks in schizophrenia underlying resting state functional magnetic resonance imaging. Current topics in medicinal chemistry, 12(21), 2415–2425.
