MRI Data Basics and the MPI-Leipzig Dataset.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# MRI data basics\n",
    "\n",
    "From [Nipype Beginner's Guide](https://miykael.github.io/nipype-beginner-s-guide/neuroimaging.html) (Some parts actually word-for-word copied; I suggest going there if you want to learn rather than have your memory jogged).\n",
    "\n",
    "The standard anatomical volume (1mm voxels) is 256x256x256 voxels.\n",
    "\n",
    "The scanner measures one plane (slice) of the brain after another (generally a horizontal plane). The resolution of the measured volume depends on the in-plane resolution, the number of slices/the thickness of the slices, and any possible gaps between them.\n",
    "\n",
    "The quality of the MRI depends on:\n",
    "> - The **resolution**\n",
    "> - The **repetition time (TR)**: time required to scan one volume\n",
    "> - The **acquisition time (TA)**: time required to scan one slice $TA = TR - (TR/n_{slices})$\n",
    "> - The **field of view (FOV)**: the extent of the slice (ex. 256mm x 256mm)\n",
    "\n",
    "The raw data format output by MRI scanners varies with the type of scanner, and is often saved in [k-space](https://en.wikipedia.org/wiki/K-space_%28magnetic_resonance_imaging%29) format. Common formats are **DICOM** and **PAR/REC**. I think DICOM and PAR/REC are scanner-specific formats encoded in k-space format, but I'm not sure.\n",
    "\n",
    "I think that the k-space formatted files are converted into **NIfTI** or **Analyze** formats. These contain an **image** and **header**. For **NIfTI** files (extension .nii-file) the image and header are in the same file. **Analyze**, which is older, produces 2 files (.img for the image and .hdr-file for the header)\n",
    "\n",
    " - The **image** is the data in a 3D matrix that contains a value for each voxel.\n",
    " - The **header** contains metadata\n",
    "   - Voxel dimension\n",
    "   - Voxel extend in each dimension\n",
    "   - The number of measured time points\n",
    "   - A transformation matrix that places the 3D matrix in the image in a 3D coordinate system\n",
    "   - etc\n",
    "\n",
    "#### sMRI data\n",
    "\n",
    "> - High-resolution images used as reference images for\n",
    ">   - corregistration\n",
    ">   - normalization\n",
    ">   - segmentation\n",
    ">   - surface reconstruction\n",
    "> - Voxel resolution ranges from 0.2-1.5mm depending on the magnetic field of the scanner (1.5T, 3T, or 7T)\n",
    "> - Grey matter structures are dark\n",
    "> - White matter structures are light\n",
    "\n",
    "#### fMRI data\n",
    "\n",
    "> Neural activity results in a characterized curve of blood oxygen level. First an initial dip in oxygen; then an increase in oxygen (peak reached after 4-6 seconds); an undershoot in oxygen levels (after 10s in figure); then eventual stabilization.\n",
    ">\n",
    "> - The signal is the Blood Oxygen Level Dependent (BOLD) response\n",
    "> - Voxel resolution is 2-4mm depending on the magnetic field of the scanner\n",
    "> - Grey matter structures are light\n",
    "> - White matter structures are dark\n",
    ">\n",
    "> Types of experimental designs:\n",
    ">\n",
    "> - **event-related design:** Stimuli are administered for a short period. BOLD responses will be short and will manifest as peaks\n",
    "> - **block design:** Multiple stimuli of similar nature are shone in a \"block\" or \"phase\" of 10-30s. The peaks will be elevated for a longer time, creating a plateau; thus the underlying activation increase should be easier to detect.\n",
    "> - **resting-state design:** Absence of stimulation. Sometimes done to analyze functional connectivity of the brain.\n",
    "\n",
    "#### dMRI data\n",
    "\n",
    "> Done to obtain information about white matter connections. Measures the diffusion of water via mean diffusivity (MD), fractional anisotropy (FA), and Tractography.\n",
    ">\n",
    "> Includes\n",
    ">\n",
    "> - diffusion tensor imaging (DTI)\n",
    "> - diffusion spectrum imaging (DSI)\n",
    "> - diffusion weighted imaging (DWI)\n",
    "> - diffusion functional MRI (dfMRI)\n",
    "\n",
    "> It is a new field in MRI and has problems with its sensitivity to correctly detect fiber tracts and their underlying orientation. Standard DTI has almost no chance of reliably detecting kissing (touching) or crossing fiber tracts.\n",
    ">\n",
    "> High-angular-resolution diffusion imaging (HARDI) and Q-ball vector analyses were developed to account for this disadvantage.\n",
    "\n",
    "------------------------------\n",
    "\n",
    "# MRI Data Analysis Steps\n",
    "\n",
    "Still from [Nipype Beginner's Guide](https://miykael.github.io/nipype-beginner-s-guide/neuroimaging.html)\n",
    "\n",
    "1. Preprocessing: spatial and temporal preprocessing of the data to prepare it for the 1st and 2nd level inferential analysis\n",
    "2. Model Specification and Estimation: specifying and estimating parameters of the statistical model\n",
    "3. Statistical Inference: Making inferences about the estimated parameters using appropriate statistical methods\n",
    "\n",
    "## Preprocessing\n",
    "\n",
    "### Slice Timing Correction (fMRI only)\n",
    "\n",
    " - Requires knowing the order of slice acquisition (top-down, bottom-up, or interleaved)\n",
    "\n",
    "### Motion Correction (fMRI only)\n",
    "\n",
    " - aka **Realignment**\n",
    " - Corrects for head movement. Aligns data to a reference time volume (usually the mean image of all timepoints, but can be the first or another time point)\n",
    " - Head movement has 6 parameters (translation along X, Y, and Z and rotation around X, Y, and Z)\n",
    " - Realignment uses an affine rigid body transformation to manipulate the data in these parameters.\n",
    " - A \"good\" subject will not translate over +/- 0.6mm and will not rotate over +/- 1 degrees (not sure, but seems to be good?)\n",
    "\n",
    "### Artifact Detection (fMRI only)\n",
    "\n",
    " - Identify and label images acquired during extreme rapid movement that should be excluded from further analysis\n",
    " - Check the translation and rotation graphs for sudden movement greater than 2 SD from the mean (or for movement greater than 1mm)\n",
    "\n",
    "### Corregistration\n",
    "\n",
    " - Aligns the functional image with the reference structural image.\n",
    " - Allows further transformations on the anatomical image, like normalization, to be directly applied to the functional one\n",
    "\n",
    "### Normalization\n",
    "\n",
    " - Used to compare the images of different subjects\n",
    " - Translates the images onto a common shape and size (maps to a reference-space)\n",
    " - Always includes a template and source images\n",
    "\n",
    "> - Template image is the standard brain in reference-space. Can be a Talairach-, MNI-, or SPM-template, or another reference image.\n",
    "> - Source image (normally a higher-resolution structural image) is used to calculate the transformation matrix, which is used for the rest of your images.\n",
    "\n",
    "### Smoothing\n",
    "\n",
    "Used on both structural and functional data\n",
    "\n",
    " - Applies a filter to the image\n",
    " - Increases the signal to noise ratio (filter highest frequencies)\n",
    " - Makes the larger scale changes more apparent\n",
    " - Reduces spatial differences between subjects, so easier to compare\n",
    " - Lose resolution!\n",
    " - Can cause functionally different regions to combine - surface based analysis with smoothing on the surface may be better choice\n",
    " - Applies a 3D Gaussian kernel to the image; amount of smoothing is determined by its full width at half maximum (FWHM) parameter\n",
    " - If you're studying a small region, a large kernel might smooth the data too much; should be smaller than/equal to the activation you're trying to detect\n",
    " - Some authors suggest using 2x{Voxel dimensions} as a reasonable starting point\n",
    "\n",
    "### Segmentation (sMRI only)\n",
    "\n",
    " - Divides the brain into neurological sections according to a template specification\n",
    "   - GM, WM, and CSF done with SPM's Segmentation\n",
    "   - Segmenting into specific functional regions and their subregions done with FreeSurfer's recon-all\n",
    "   - Can aid normalization, use specific segmentation as a mask, use segmentation as definition of ROI\n",
    "\n",
    "## Model Specification and Estimation\n",
    "\n",
    "To test a hypothesis, we need to specify a model that incorporates the hypothesis and accounts for the expected function of the BOLD signal, the movement during measurement, experiment specify parameters and other regressors and co-variants. This is usually represented by a Generalized Linear Model (**GLM**).\n",
    "\n",
    "### Making a GLM\n",
    "\n",
    "Describes a response $y$ (ex: the BOLD response in a voxel), at the time points with data in terms of all of its contributing factors ($X\\beta$) and error ($\\epsilon$).\n",
    "\n",
    "$$\n",
    "y = X\\beta + \\epsilon\n",
    "$$\n",
    "\n",
    " - $y$ (dependent variable)\n",
    " - $X$ (independent variable; predictor; experimental conditions, stimulus information, expected shape of BOLD response) These are stored in a design matrix\n",
    " - $\\beta$ (parameters; regression coefficient/beta weights; quantifies how much each predictor independently influences the dependent variable)\n",
    " - $\\epsilon$ (error; assumed to be normally distributed)\n",
    "\n",
    "\n",
    "> **Caveats for the GLM approach:**\n",
    ">\n",
    "> - We need to take the time delay and the HRF (hemodynamic response function) shape of the BOLD response into account when creating the design matrix\n",
    "> - We need to high pass filter the data and add time regressors of 1st, 2nd, ... order to correct for low-frequency drifts in the measured data (the low-frequency noise are caused by scanner drift and other non-experimental effects)\n",
    "> - The high pass filter is established by setting up discrete cosine functions over the time period of data acquisition.\n",
    "> - Each regressors in the model decreases the degrees of freedom in statistical tests\n",
    "\n",
    "### Estimating the Model\n",
    "\n",
    "Applying the model on the time course of each and every voxel. SPM will create images every time an analysis is performed.\n",
    "\n",
    "> - **Beta images:** images of estimated regression coefficients (parameter estimate). These contain information about the size of the effect of interest. A given voxel in each beta image will have value related to the size of effect for the explanatory variables.\n",
    "> - **Error image (ResMS-image):** Residual sum of squares/variance image. Measurement of within-subject error at the 1st level or between-subject error at the 2nd level analysis. Used to produce spmT images (see below).\n",
    "> - **con images (con-images):** Produced during contrast estimation from linearly combining beta images\n",
    "> - **T images (spmT-images):** Produced during contrast estimation from combining the beta values of the con-images with error values of the ResMS-image to calculate the t-value at each voxel.\n",
    "\n",
    "## Statistical inferences\n",
    "\n",
    "**1st level analysis (within-subject analysis):** The data does not need to be normalized. The design matrix for this level controls for movement, respiration, heart beat, etc.\n",
    "\n",
    "**2nd level analysis (between-subject analysis):** Requires normalization and transformation onto reference-space of subject specific data. Requires the contrasts of the 1st level analysis. The design matrix controls for age, gender, socio-economic parameters. Here we specify the group assignment of each subject.\n",
    "\n",
    "### Contrast Estimation (for both 1st and 2nd level analysis)\n",
    "\n",
    "Specify how to weight the different regressors of your design matrix and combine them into one image. Essentially, how can you numerically weight what you are interested in such that it's on the same scale.\n",
    "\n",
    "> Examples:  \n",
    "> faces (+1) vs resting (-1)  \n",
    "> session 1 (+1) vs session 2 (-1))\n",
    "\n",
    "### Thresholding\n",
    "\n",
    "Specify the level of significance you test your data on. Correct for multiple comparison and specify the parameters of the result you are looking for.\n",
    "\n",
    "> Examples:\n",
    ">\n",
    "> - FWE-correction: the family-wise error correction can correct for multiple comparisons\n",
    "> - p-value\n",
    "> - voxel extend: specify the minimum size of a significant cluster with the lowest number of voxels it needs to contain (ex 100)\n",
    "\n",
    "------------------------------------\n",
    "\n",
    "# Brain screening protocol (MRI)\n",
    "\n",
    "Simple non-contrast MRI protocol - mainly for use in pain, not studies.\n",
    "\n",
    "T1 weighted\n",
    " - plane: saggital (or volumetric 3D)\n",
    " - sequence: fast-spin echo (T1 FSE) or gradient (T1 MPRAGE)\n",
    " - purpose: anatomical overview\n",
    " \n",
    "T2 weighted\n",
    " - plane: axial\n",
    " - sequence: T2 FSE\n",
    " - purpose: evaluations of basal cisterns, ventricular system and subdural spaces, and good visualization of flow voids in vessels\n",
    " \n",
    "FLAIR\n",
    " - plane: axial\n",
    " - sequence: FLAIR\n",
    " - purpose: assessment of white matter disorders\n",
    "\n",
    "Diffusion weighted imaging (DWI)\n",
    " - plane: axial\n",
    " - sequence: DWI B=0, B=1000 and ADC\n",
    " - purpose: multiple (identification of stroke to demyelination)\n",
    "\n",
    "Susceptibility weighted imaging (SWI)\n",
    " - plane: axial\n",
    " - sequence: susceptibility weighted imaging (ideal) or T2\\*\n",
    " - purpose: identify blood products or calcification\n",
    " \n",
    "\n",
    "------------------------------------\n",
    "\n",
    "# Structure of data in the MPI-Leipzig dataset (ds000221)\n",
    "\n",
    "The Brain Imaging Data Structure ([BIDS](http://bids.neuroimaging.io/bids_spec1.0.1.pdf) defines an organizational system relying on folder and file names.\n",
    "\n",
    "The DICOM files are organized:\n",
    "> dicomdir/\n",
    "> > [timestamp]/\n",
    "> > > .dcm files\n",
    "\n",
    "The other files follow\n",
    "> participants.tsv\n",
    "> sub-01/\n",
    "> > anat/\n",
    "> > > .nii.gz\n",
    "> >\n",
    "> > func/\n",
    "> > > .nii.gz  \n",
    "> > > .json\n",
    "> >\n",
    "> > dwi/\n",
    "> > > .nii.gz  \n",
    "> > > .json  \n",
    "> > > .bval  \n",
    "> > > .bvec\n",
    "\n",
    "## Data types\n",
    "\n",
    "Functional group of different types of MRI data\n",
    " - func: task baased and resting state functional MRI\n",
    "     - Follows the template:\n",
    "       `sub-<participant_label>[_ses-<session_label>]_task-<task_label>[_aqc_<label>][_rec-<label>][_run-<index>]_bold.nii[.gz]`\n",
    " - dwi: diffusion weighted imaging\n",
    "     - Follows the templates:\n",
    "         - `sub-<participant_label>[_ses-<session_label>][_aqc_<label>][_run-<index>]_dwi.nii[.gz]`\n",
    "         - `sub-<participant_label>[_ses-<session_label>][_aqc_<label>][_run-<index>]_dwi.bval`\n",
    "         - `sub-<participant_label>[_ses-<session_label>][_aqc_<label>][_run-<index>]_dwi.bvec`\n",
    "         - `sub-<participant_label>[_ses-<session_label>][_aqc_<label>][_run-<index>]_dwi.json`\n",
    "     - The bvec and bval files are in the FSL format.     \n",
    " - fmap: field inhomogeneity mapping data (field maps)\n",
    " - anat: structural imaging\n",
    "     - Follows the template:\n",
    "       `sub-<participant_label>[_ses-<session_label>][_acq-<label>][_rec-<label>][_run-<index>]_<modality_label>.nii[.gz]`\n",
    "     \n",
    "\n",
    "> **NOTES** \n",
    "> - Derivatives of the original data should be kept in a derivatives/ folder\n",
    "> - Elapsed time is in seconds\n",
    "> - Frequency is in Hz\n",
    "> - YYYY-MM-DDThh:mm:ss\n",
    "\n",
    "##  Modality labels in openfmri/ds000221\n",
    "\n",
    "### [FLAIR](https://en.wikipedia.org/wiki/Fluid-attenuated_inversion_recovery)\n",
    "\n",
    "FLAIR files.\n",
    "\n",
    "It's a pulse sequence used in MRI; can be 2D or 3D. Can be used to suppress CSF effects. Helpful for studying lacunar infarction, MS plaques, Subarachnoid heamorrhage, head trauma, and meningitis. The inversion time (TI) chooses which signal is nullified.\n",
    "\n",
    "Also in [Radiopaedia](https://radiopaedia.org/articles/fluid-attenuation-inversion-recovery). Removes signal from CSF. It is part of almost all protocols for imaging the brain. Useful in detection of subtle changes at the periphery of the hemispheres and in the periventricular region.\n",
    "\n",
    "### [T1map](https://mandymejia.wordpress.com/2014/02/04/what-is-quantitative-mri-and-why-does-it-matter/)\n",
    "\n",
    "Quantitative T1 map files. \n",
    "\n",
    "T1 relaxation time is the time it takes for the protons to \"relax\" back to their original position. It can vary for different tissue types: MS patients have widespread increases in T1 relaxation times. The map of the T1 relaxation times is the \"T1map\".\n",
    "\n",
    "T1 map can have a high level of noise and poor signal to noise ratio when compared to T1-weighted image.\n",
    "\n",
    "T1 mapping may be used to\n",
    " - optimize parameters for a sequence (Ernst angle is $acos(e^{-TR/T1})$) ([jbarral](http://www-mrsrl.stanford.edu/~jbarral/t1map.html))\n",
    " - monitor diseased tissue\n",
    " - measure Ktrans in DCE-MRI\n",
    " - derive other quantitative parameters (bound pool fractions)\n",
    " - determination of local white matter and grey matter ([Deichmann 2005](https://www.ncbi.nlm.nih.gov/pubmed/15968665))\n",
    "\n",
    "### [T1w](https://radiopaedia.org/articles/t1-weighted-image)\n",
    "\n",
    "T1 weighted image.\n",
    "\n",
    "T1 weighted image is one of the basic pulse sequences in MRI and demonstrates differences in the T1 relaxation times of vectors.\n",
    "\n",
    "It relies on the longitudinal relaxation of a tissue's net magnetisation vector (NMV). (How long the protons' spins realign with the main magnetic field). \n",
    "\n",
    "It tends to have short **TE** (echo time - the time between the application of the radiofrequency excitation pulse and the peak of the signal induced in the coil; ms) and **TR** (repitition time - the time from the application of an excitation pulse to the application of the next pulse; ms) times. You need to pick a **TR** shorter than the tissues' recovery time to allow for contrast between the tissues.\n",
    "\n",
    "Fat is bright; water is dark.\n",
    "\n",
    "T1-weighted sequences include:\n",
    " - T1W spin echo (SE)\n",
    " - T1W gradient echo (GRE)\n",
    " - gadolinium postcontrast sequences (gradient echo sequences)\n",
    " - time of flight 2D or 3D MR angiography sequences\n",
    " - contrast-enhanced sequence (in-phase and out-of-phase)\n",
    "\n",
    "\n",
    "### defacemask \n",
    "\n",
    "A defacing mask. \"Mask used for defacing\"\n",
    "\n",
    "### [bold](https://radiopaedia.org/articles/bold-imaging)\n",
    "\n",
    "Imaging data acquired during BOLD imaging. If it is a task based fMRI, a corresponding task events file will be provided: `<matches>_events.tsv`. Look for a similiarly named `.json` file for metadata.\n",
    "\n",
    "Deoxygenated hemoglobin is paramagnetic whereas oxygenated hemoglobin is not, and therefore the former will cause local dephasing of protons, and thus reduce the returned signal from the tissues in the immediate vicinity. Heavily T2\\* weighted sequences are used to detect this change, which is in the order of 1-5%. \n",
    "\n",
    "Limitations:\n",
    " - cerebral blood flow (CBF) is an indirect marker of activity\n",
    " - the smallest unit of the brain that is able to have its blood flow individually regulated is ~mm in diameter\n",
    " - CBF increases in response to increased activity, so there is a 2-6s lag\n",
    " - T2\\* sequences are susceptible to field inhomogeneity due to bone-air interface, blood products, rapid flow in large veins, and metal\n",
    " - The change detected is small, so small movement artefacts can lead to poor images\n",
    "\n",
    "### [T2w](https://radiopaedia.org/articles/t2-weighted-image)\n",
    "\n",
    "T2 weighted image.\n",
    "\n",
    "Highlights differences in the T2 relaxation time (spins decaying from the aligned precession in the transverse plane) of tissues. Requires long **TE** and **TR** times.\n",
    "\n",
    "Each tissue has an inherent T2 value (external factors such as magnetic field inhomogeneity can increase the T2 relaxation time - these are captured in T2\\* weighted imaging).\n",
    "\n",
    "Does not work well with paramagnetic contrast agents, unlike T1.\n",
    "\n",
    "Fat is intermediate bright. Fluid is bright.\n",
    "\n",
    "### [dwi](https://radiopaedia.org/articles/diffusion-weighted-imaging-1)\n",
    "\n",
    "Possibly an isotropic diffusion map.\n",
    "\n",
    "Water diffusion in the brain is hindered by cell structures, so water diffusion can occur:\n",
    " - within the intracellular fluid\n",
    "     - within the cytoplasm\n",
    "     - within organelles\n",
    " - within extracellular fluid\n",
    "     - interstitial fluid\n",
    "     - intravascular\n",
    "     - lymphatic\n",
    "     - in cavities, eg. ventricles\n",
    " - between intra- and extracellular compartments\n",
    "\n",
    "The **Apparent diffusion coefficient (ADC)** ($mm^2/s$) is a measure of the magnitude of diffusion of water within tissue. They are calculated automatically by \"the software\" and then displayed as a parametric map. ADC measurements are recorded for a given region by drawing regions of interest on the ADC map.\n",
    "\n",
    "$ ADC < 1.0*10^{-3} $ indicates restriction, but depends on the organ being imaged and the pathology.\n",
    "\n",
    "Useful values ($10^{-6} mm^2/s$):\n",
    " - white matter: 670-800\n",
    " - cortical grey matter: 800-1000\n",
    " - deep grey matter: 700-850\n",
    " - CSF: 3000-3400\n",
    "\n",
    "**B value** measures the degree of diffusion weighting applied, indicating the amplitude ($G$), time of applied gradients ($\\delta$), and duration between the paired gradients ($\\Delta$).\n",
    "\n",
    "$$\n",
    " b = \\gamma^2G^2\\delta^2(\\Delta-\\delta/3)\n",
    "$$\n",
    "\n",
    "The ADC is calculated using $b$ values. A good rule of thumb is to choose the $b$ value such that $b*\\text{ADC} \\approx 1$\n",
    "\n",
    "Diffusion tensor imaging is an extension of diffusion weighted imaging. I'm not sure how one gets from A to B\n",
    "\n",
    "### [fieldmap](https://en.wikibooks.org/wiki/Neuroimaging_Data_Processing/Field_map_correction)\n",
    "\n",
    "Probably a field map.\n",
    "\n",
    "A field map is an image of the intensity of the magnetic field across space. The magnetic field may change around air-filled sinuses and bone, which can lead to distortion of the image.\n",
    "\n",
    "SPM provides geometric distortion correction:\n",
    " - if the distortion is along the anterior-posterior axis, not affecting the symmetry of the brain, use `Realign & Unwarp`.\n",
    " - if not, use the `FieldMap toolbox` and `VDM utility`.\n",
    "\n",
    "If included with a BOLD file, according to BIDS, the BOLD file will have the metadata `PhaseEncodingDirection`, `EffectiveEchoSpacing`, `TotalReadoutTime`, and `EchoTime`.\n",
    "\n",
    "If included with a dwi file, according to BIDS, the dwi file will have the metadata `PhaseEncodingDirection`, `EffectiveEchoSpacing`, and `EchoTime`.\n",
    "\n",
    "Optionally includes an `IntendedFor` field with one or more relative paths stripped from `sub-<label>`.\n",
    "\n",
    "### [GRE](https://radiopaedia.org/articles/gradient-echo-sequences-1)\n",
    "\n",
    "A gradient echo sequence file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}