Explore resting state BOLD fMRI data (work in progress ...)¶

BMED360-2021 03-Resting-state-fmri-explore.ipynb

(using selected datasets from the Brain-Gut BIDS collection and open data from Nilearn and OpenNeuro)

Learning objectives:¶

Three major topics

• The nature of 4D (3D+time) resting state BOLD fMRI data¶

See also: What is meant by resting state fMRI? How is it used? in [mriquestions.com]
and Gili Karni's Identifying Resting-State Networks from fMRI Data Using ICAs - A hands-on tutorial: from data collection to feature extraction [here] and [here]
List of functional connectivity software on Wikipedia is [here].

• BIDS: Brain Imaging Data Structure¶

The Brain Imaging Data Structure (BIDS) is a standard prescribing a formal way to name and organize MRI data and metadata in a file system that simplifies communication and collaboration between users and enables easier data validation and software development through using consistent paths and naming for data files. (See: https://bids.neuroimaging.io and https://training.incf.org/lesson/introduction-brain-imaging-data-structure-bids).
We will be using pybids (https://github.com/bids-standard/pybids) (included by pip install pybids in the environment.ymlfile) to access the ./data/bids_bg_bmed360 information (e.g. participants.tsv)

• Open science with Open data (BIDS validated)¶

E.g. Explore the OpenNeuro collection: https://openneuro.org/datasets/ds000030/versions/1.0.0 + https://f1000research.com/articles/6-1262/v2
How to download such a large dataset (Files: 9010, Size: 79.28GB, Subjects: 272) using the AWS Command Line Interface (https://aws.amazon.com/cli). The derivatives are the preprocessed data using fmriprep.

aws s3 sync --no-sign-request s3://openneuro.org/ds000030 ds000030-download/
aws s3 sync --no-sign-request s3://openneuro/ds000030/ds000030_R1.0.5/uncompressed/derivatives/ ds000030-derivatives-download/

Bilder, R and Poldrack, R and Cannon, T and London, E and Freimer, N and Congdon, E and Karlsgodt, K and Sabb, F (2020). UCLA Consortium for Neuropsychiatric Phenomics LA5c Study. OpenNeuro. [Dataset] doi: 10.18112/openneuro.ds000030.v1.0.0

Some instructive videos on resting state fMRI and the BIDS - Brain Imaging Data Structure¶

Setup¶

For using Google Colab¶

--> (some of) the following libraries must be pip installed):

Set up libraries and folders¶

Check your platform for running this notebook

Download the data¶

Download a files (not on GitHub repository) from Google Drive using gdown (https://github.com/wkentaro/gdown)

BIDS: Brain Imaging Dataset Specification¶

Recently, there has been growing interest to share datasets across labs and even on public repositories such as openneuro. In order to make this a succesful enterprise, it is necessary to have some standards in how the data are named and organized. Historically, each lab has used their own idiosyncratic conventions, which can make it difficult for outsiders to analyze. In the past few years, there have been heroic efforts by the neuroimaging community to create a standardized file organization and naming practices. This specification is called BIDS for Brain Imaging Dataset Specification.

As you can imagine, individuals have their own distinct method of organizing their files. Think about how you keep track of your files on your personal laptop (versus your friend). This may be okay in the personal realm, but in science, it's best if anyone (especially yourself 6 months from now!) can follow your work and know which files mean what by their titles.

Here's an example of non-Bids versus BIDS dataset found in this paper:

Here are a few major differences between the two datasets:

1. In BIDS, files are in nifti format (not dicoms).
2. In BIDS, scans are broken up into separate folders by type of scan(functional versus anatomical versus diffusion weighted) for each subject.
3. In BIDS, JSON files are included that contain descriptive information about the scans (e.g., acquisition parameters)

Not only can using this specification be useful within labs to have a set way of structuring data, but it can also be useful when collaborating across labs, developing and utilizing software, and publishing data.

In addition, because this is a consistent format, it is possible to have a python package to make it easy to query a dataset. We recommend using pybids.

The dataset we will be working with has already been converted to the BIDS format (see download localizer tutorial).

You may need to install pybids to query the BIDS datasets using following command !pip install pybids (in colab), or
conda install -c aramislab pybids (old version 0.5.1 !) in the bmed360 conda environment.

The BIDSLayout¶

Pybids is a package to help query and navigate a neurogimaging dataset that is in the BIDs format. At the core of pybids is the BIDSLayout object. A BIDSLayout is a lightweight Python class that represents a BIDS project file tree and provides a variety of helpful methods for querying and manipulating BIDS files. While the BIDSLayout initializer has a large number of arguments you can use to control the way files are indexed and accessed, you will most commonly initialize a BIDSLayout by passing in the BIDS dataset root location as a single argument.

Notice we are setting derivatives=True. This means the layout will also index the derivatives sub folder, which might contain preprocessed data, analyses, or other user generated files.

When we initialize a BIDSLayout, all of the files and metadata found under the specified root folder are indexed. This can take a few seconds (or, for very large datasets, a minute or two). Once initialization is complete, we can start querying the BIDSLayout in various ways. The main query method is .get(). If we call .get() with no additional arguments, we get back a list of all the BIDS files in our dataset.

Let's return the first 10 files

When you call .get() on a BIDSLayout, the default returned values are objects of class BIDSFile. A BIDSFile is a lightweight container for individual files in a BIDS dataset.

Here are some of the attributes and methods available to us in a BIDSFile (note that some of these are only available for certain subclasses of BIDSFile; e.g., you can't call get_image() on a BIDSFile that doesn't correspond to an image file!):

• .path: The full path of the associated file
• .filename: The associated file's filename (without directory)
• .dirname: The directory containing the file
• .get_entities(): Returns information about entities associated with this BIDSFile (optionally including metadata)
• .get_image(): Returns the file contents as a nibabel image (only works for image files)
• .get_df(): Get file contents as a pandas DataFrame (only works for TSV files)
• .get_metadata(): Returns a dictionary of all metadata found in associated JSON files
• .get_associations(): Returns a list of all files associated with this one in some way

Let's explore the first file in a little more detail.

If you want a summary of all the files in your BIDSLayout, but don't want to have to iterate BIDSFile objects and extract their entities, you can get a nice bird's-eye view of your dataset using the to_df() method.

Let's have a look at our fMRI data¶

First, a few handy functions:

Plot an image representation of voxel intensities across time

Power, J. D. A simple but useful way to assess fMRI scan qualities. Neuroimage 2017;154:150-158 [link]

FSL: Brain Extraction Tool (BET2)¶

CPU times: user 91.3 ms, sys: 38.7 ms, total: 130 ms
Wall time: 8.16 s

--> sub-102_ses-1_T1w_brain.nii.gz and sub-102_ses-1_T1w_brain_mask.nii.gz

FSL: FEAT¶

(FEAT = FMRI Expert Analysis Tool - for preprocessing and statistical analysis of FMRI data)

First-level analysis | Preprocessing

Data (4D): ~/prj/BMED360/Lab6-Networks-Graphs/sub-103_ses-1_task-rest_bold.nii.gz
Total volumes: 240; Delete volumes:0; TR: 2.0 s; High pass filter cutoff: 60 s
Pre-stats:
Motion correction (MCFLIRT),
Slice timing correction: None
BET brain extraction: Yes
Spatial smoothing: FWHM 5 mm
Temporal filtering Highpass, MELODIC ICA data exploration
Registration:
Man structural image: ~/prj/BMED360/Lab6-Networks-Graphs/sub-102_ses-1_T1w_brain.nii.gz; Linear Normal search, BBR
Standard space: /usr/local/fsl/data/standard/MNI151_T1_2mm_brain.nii.gz; Linear: Normal search, 12 DOF
Misc: Brain/background threshold: 10%, Noise level: 0.66; Temporal smoothness: 0.34; Z-threshold: 5.3

Preprocessing:Stage 1) --> sub-102_ses-1_task-rest_bold_filtered_func_data.nii.gz¶

Filtered function ICA report --> sub-102_ses-1_task-rest_bold.feat/filtered_func_data.ica¶

FSL: MELODIC¶

(MELODIC = Multivariate Exploratory Linear Optimized Decomposition into Independent Components)

Uses Independent Component Analysis (ICA) to decompose a single or multiple 4D data sets into different spatial and temporal components. For ICA group analysis, MELODIC uses either Tensorial Independent Component Analysis (TICA, where data is decomposed into spatial maps, time courses and subject/session modes) or a simpler temporal concatenation approach. MELODIC can pick out different activation and artefactual components without any explicit time series model being specified.

A paper investigating resting-state connectivity using independent component analysis has been published in Philosophical Transactions of the Royal Society. For detail, see a technical report PDF.

Automatically generated FSL reports¶

Registration¶

Analysis methods:
FMRI data processing was carried out using FEAT (FMRI Expert Analysis Tool) Version 6.00, part of FSL (FMRIB's Software Library, www.fmrib.ox.ac.uk/fsl). Registration to high resolution structural and/or standard space images was carried out using FLIRT [Jenkinson 2001, 2002].

References:
[Jenkinson 2001] M. Jenkinson and S.M. Smith. A Global Optimisation Method for Robust Affine Registration of Brain Images. Medical Image Analysis 5:2(143-156) 2001.
[Jenkinson 2002] M. Jenkinson, P. Bannister, M. Brady and S. Smith. Improved Optimisation for the Robust and Accurate Linear Registration and Motion Correction of Brain Images. NeuroImage 17:2(825-841) 2002.

Pres-stats¶

Analysis methods:
FMRI data processing was carried out using FEAT (FMRI Expert Analysis Tool) Version 6.00, part of FSL (FMRIB's Software Library, www.fmrib.ox.ac.uk/fsl). The following pre-statistics processing was applied; motion correction using MCFLIRT [Jenkinson 2002]; non-brain removal using BET [Smith 2002]; spatial smoothing using a Gaussian kernel of FWHM 5mm; grand-mean intensity normalisation of the entire 4D dataset by a single multiplicative factor; highpass temporal filtering (Gaussian-weighted least-squares straight line fitting, with sigma=30.0s). ICA-based exploratory data analysis was carried out using MELODIC [Beckmann 2004], in order to investigate the possible presence of unexpected artefacts or activation.

References:
[Jenkinson 2002] M. Jenkinson and P. Bannister and M. Brady and S. Smith. Improved optimisation for the robust and accurate linear registration and motion correction of brain images. NeuroImage 17:2(825-841) 2002.
[Smith 2002] S. Smith. Fast Robust Automated Brain Extraction. Human Brain Mapping 17:3(143-155) 2002.
[Beckmann 2004] C.F. Beckmann and S.M. Smith. Probabilistic Independent Component Analysis for Functional Magnetic Resonance Imaging. IEEE Trans. on Medical Imaging 23:2(137-152) 2004.

MELODIC data exploration report¶

Analysis methods:
Analysis was carried out using Probabilistic Independent Component Analysis [Beckmann 2004] as implemented in MELODIC (Multivariate Exploratory Linear Decomposition into Independent Components) Version 3.15, part of FSL (FMRIB's Software Library, www.fmrib.ox.ac.uk/fsl). The following data pre-processing was applied to the input data: masking of non-brain voxels; voxel-wise de-meaning of the data; normalisation of the voxel-wise variance; Pre-processed data were whitened and projected into a 60-dimensional subspace using probabilistic Principal Component Analysis where the number of dimensions was estimated using the Laplace approximation to the Bayesian evidence of the model order [Minka 2000, Beckmann 2004]. The whitened observations were decomposed into sets of vectors which describe signal variation across the temporal domain (time-courses) and across the spatial domain (maps) by optimising for non-Gaussian spatial source distributions using a fixed-point iteration technique [Hyvärinen 1999]. Estimated Component maps were divided by the standard deviation of the residual noise and thresholded by fitting a mixture model to the histogram of intensity values [Beckmann 2004].

References:
[Hyvärinen 1999] A. Hyvärinen. Fast and Robust Fixed-Point Algorithms for Independent Component Analysis. IEEE Transactions on Neural Networks 10(3):626-634, 1999.
[Beckmann 2004] C.F. Beckmann and S.M. Smith. Probabilistic Independent Component Analysis for Functional Magnetic Resonance Imaging. IEEE Transactions on Medical Imaging 23(2):137-152 2004.
[Everson 2000] R. Everson and S. Roberts. Inferring the eigenvalues of covariance matrices from limited, noisy data. IEEE Trans Signal Processing, 48(7):2083-2091, 2000
[Tipping 1999] M.E. Tipping and C.M.Bishop. Probabilistic Principal component analysis. J Royal Statistical Society B, 61(3), 1999.
[Beckmann 2001] C.F. Beckmann, J.A. Noble and S.M. Smith. Investigating the intrinsic dimensionality of FMRI data for ICA. In Seventh Int. Conf. on Functional Mapping of the Human Brain, 2001.
[Minka 2000] T. Minka. Automatic choice of dimensionality for PCA. Technical Report 514, MIT Media Lab Vision and Modeling Group, 2000.

• Single-session ICA: This will perform standard 2D ICA on each of the input files. The input data will each be represented as a 2D time x space matrix. MELODIC then de-composes each matrix separately into pairs of time courses and spatial maps. The original data is assumed to be the sum of outer products of time courses and spatial maps. All the different time courses (one per component) will be saved in the mixing matrix melodic_mix and all the spatial maps (one per component) will be saved in the 4D file melodic_IC.
  When using separate analyses, MELODIC will attempt to find components which are relevant and non-Gaussian relative to the residual fixed-effects within session/subject variation. It is recommended to use this option in order to check for session-specific effects (such as MR artefacts). You will need to use this option if you want to perform MELODIC denoising using fsl_regfilt. When using single-session ICA the component are ordered in order of decreasing amounts of uniquely explained variance.

Melodic results will be in filtered_func_data.ica

To view the output report point your web browser at filtered_func_data.ica/report/00index.html

CPU times: user 8.85 s, sys: 2.66 s, total: 11.5 s
Wall time: 14min 30s

Selected components¶

(Thresholded IC map, alternative hypothesis test at p > 0.5: IC_5_thresh.png / Timecourse: t5.png / Powerspectrum of timecourse: f5.png)

MELODIC Component 5 | 2.33 % of explained variance; 1.83 % of total variance
MELODIC Component 13 | 1.90 % of explained variance; 1.50 % of total variance
MELODIC Component 18 | 1.84 % of explained variance; 1.45 % of total variance

MELODIC Component 5¶

Thresholded IC map (IC_5_thresh.png), alternative hypothesis test at p > 0.5

Timecourse (t5.png)

Powerspectrum of timecourse (f5.png)

Plot the power spectral density and the magnitude spectrum.¶

Compute the magnitude spectrum of x

The power spectral density 𝑃𝑥𝑥 by Welch's average periodogram method, and then compute the magnitude spectrum of x.

Connectome example from Nilearn¶

This example is adopted form https://nilearn.github.io/auto_examples/03_connectivity/plot_sphere_based_connectome.html

Extract signals on spheres and plot a connectome¶

This example shows how to extract signals from spherical regions. We show how to build spheres around user-defined coordinates, as well as centered on coordinates from the Power-264 atlas [1], and the Dosenbach-160 atlas [2].

References

[1] Power, Jonathan D., et al. “Functional network organization of the human brain.” Neuron 72.4 (2011): 665-678.

[2] Dosenbach N.U., Nardos B., et al. “Prediction of individual brain maturity using fMRI.”, 2010, Science 329, 1358-1361.

We estimate connectomes using two different methods: sparse inverse covariance (not shown) and partial_correlation, to recover the functional brain networks structure.

We’ll start by extracting signals from Default Mode Network regions and computing a connectome from them.

Coordinates of Default Mode Network¶¶

Extracts signal from sphere around DMN seeds¶

We can compute the mean signal within spheres of a fixed radius around a sequence of (x, y, z) coordinates with the object nilearn.input_data.NiftiSpheresMasker. The resulting signal is then prepared by the masker object: Detrended, band-pass filtered and standardized to 1 variance.

Display time series¶

Compute partial correlation matrix¶

Using object nilearn.connectome.ConnectivityMeasure: Its default covariance estimator is Ledoit-Wolf, allowing to obtain accurate partial correlations.

Display connectome¶

We display the graph of connections with nilearn.plotting.plot_connectome

Display connectome with hemispheric projections. Notice (0, -52, 18) is included in both hemispheres since x == 0.

3D visualization in a web browser¶

An alternative to nilearn.plotting.plot_connectome is to use nilearn.plotting.view_connectome, which gives more interactive visualizations in a web browser. See 3D Plots of connectomes for more details.

We have represented the functional connectivity of DMN as a weighted undirected graph:¶

A final example with an OpenNeuro dataset¶

Feature extraction using ICA from resting state networks¶

(Using preprocessed data with fmriprep from HC subject sub-10228 in the OpenNeuro collection ds000030)

Adopted from Gili Karni, Princeton University:
https://github.com/karnigili/Computational-Neuroscience/blob/master/rs_fmri_networks_ICA.ipynb (MIT License)

Feature extraction using ICA from resting state networks¶

Preparing for atlas-based analysis¶

We 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 atals is based on MNI defined regions and the Smith's atlas reflecs regions resulting from an independent analyses of resting and activation brain dynamics.

Smith's includes 10 resting-state networks (rsns)

ICA analysis (Varoquaux et al., 2010)¶

Varoquaux G. et al. “A group model for stable multi-subject ICA on fMRI datasets”, NeuroImage Vol 51 (2010), p. 288-299

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

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

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.

Goldstein JM, Seidman LJ, Makris N, Ahern T, O'Brien LM, Caviness VS Jr, Kennedy DN, Faraone SV, Tsuang MT. Hypothalamic abnormalities in schizophrenia: sex effects and genetic vulnerability. Biol Psychiatry. 2007 Apr 15;61(8):935-45

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. [link]