{
"cells": [
{
"cell_type": "markdown",
"source": [
"## Diffusion Tensor Imaging (DTI)\n",
"\n",
"Diffusion tensor imaging or \"DTI\" refers to images describing diffusion with a tensor model. DTI is derived from preprocessed diffusion weighted imaging (DWI) data. First proposed by Basser and colleagues ([Basser, 1994](https://www.ncbi.nlm.nih.gov/pubmed/8130344)), the diffusion tensor model describes diffusion characteristics within an imaging voxel. This model has been very influential in demonstrating the utility of the diffusion MRI in characterizing the microstructure of white matter and the biophysical properties (inferred from local diffusion properties). The DTI model is still a commonly used model to investigate white matter.\n",
"\n",
"The tensor models the diffusion signal mathematically as:\n",
"\n",
"\n",
"\n",
"Where  is a unit vector in 3D space indicating the direction of measurement and b are the parameters of the measurement, such as the strength and duration of diffusion-weighting gradient.  is the diffusion-weighted signal measured and  is the signal conducted in a measurement with no diffusion weighting.  is a positive-definite quadratic form, which contains six free parameters to be fit. These six parameters are:\n",
"\n",
"\n",
"\n",
"The diffusion matrix is a variance-covariance matrix of the diffusivity along the three spatial dimensions. Note that we can assume that the diffusivity has antipodal symmetry, so elements across the diagonal of the matrix are equal. For example: . This is why there are only 6 free parameters to estimate here.\n",
"\n",
"Tensors are represented by ellipsoids characterized by calculated eigenvalues () and eigenvectors () from the previously described matrix. The computed eigenvalues and eigenvectors are normally sorted in descending magnitude (i.e. ). Eigenvalues are always strictly positive in the context of dMRI and are measured in mm^2/s. In the DTI model, the largest eigenvalue gives the principal direction of the diffusion tensor, and the other two eigenvectors span the orthogonal plane to the former direction.\n",
"\n",
"\n",
"_Adapted from Jelison et al., 2004_\n",
"\n",
"In the following example, we will walk through how to model a diffusion dataset. While there are a number of diffusion models, many of which are implemented in `DIPY`. However, for the purposes of this lesson, we will focus on the tensor model described above.\n",
"\n",
"### Reconstruction with the `dipy.reconst` module\n",
"\n",
"The `reconst` module contains implementations of the following models:\n",
"\n",
"* Tensor (Basser et al., 1994)\n",
"* Constrained Spherical Deconvolution (Tournier et al. 2007)\n",
"* Diffusion Kurtosis (Jensen et al. 2005)\n",
"* DSI (Wedeen et al. 2008)\n",
"* DSI with deconvolution (Canales-Rodriguez et al. 2010)\n",
"* Generalized Q Imaging (Yeh et al. 2010)\n",
"* MAPMRI (Özarslan et al. 2013)\n",
"* SHORE (Özarslan et al. 2008)\n",
"* CSA (Aganj et al. 2009)\n",
"* Q ball (Descoteaux et al. 2007)\n",
"* OPDT (Tristan-Vega et al. 2010)\n",
"* Sparse Fascicle Model (Rokem et al. 2015)\n",
"\n",
"The different algorithms implemented in the module all share a similar conceptual structure:\n",
"\n",
"* `ReconstModel` objects (e.g. `TensorModel`) carry the parameters that are required in order to fit a model. For example, the directions and magnitudes of the gradients that were applied in the experiment. `TensorModel` objects have a `fit` method, which takes in data, and returns a `ReconstFit` object. This is where a lot of the heavy lifting of the processing will take place.\n",
"* `ReconstFit` objects carry the model that was used to generate the object. They also include the parameters that were estimated during fitting of the data. They have methods to calculate derived statistics, which can differ from model to model. All objects also have an orientation distribution function (`odf`), and most (but not all) contain a `predict` method, which enables the prediction of another dataset based on the current gradient table.\n"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"### Reconstruction with the DTI model\n",
"\n",
"Let's get started! First, we will need to grab **preprocessed** DWI files and load them! We will also load in the anatomical image to use as a reference later on!"
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"import bids\n",
"from bids.layout import BIDSLayout\n",
"from dipy.io.gradients import read_bvals_bvecs\n",
"from dipy.core.gradients import gradient_table\n",
"from nilearn import image as img\n",
"\n",
"bids.config.set_option('extension_initial_dot', True)\n",
"\n",
"deriv_layout = BIDSLayout(\"../data/ds000221/derivatives\", validate=False)\n",
"subj = \"010006\"\n",
"\n",
"# Grab the transformed t1 file for reference\n",
"t1 = deriv_layout.get(subject=subj, space=\"dwi\",\n",
" extension='nii.gz', return_type='file')[0]\n",
"\n",
"# Recall the preprocessed data is no longer in BIDS - we will directly grab these files\n",
"dwi = \"../data/ds000221/derivatives/uncorrected_topup_eddy/sub-%s/ses-01/dwi/dwi.nii.gz\" % subj\n",
"bval = \"../data/ds000221/sub-%s/ses-01/dwi/sub-%s_ses-01_dwi.bval\" % (\n",
" subj, subj)\n",
"bvec = \"../data/ds000221/derivatives/uncorrected_topup_eddy/sub-%s/ses-01/dwi/dwi.eddy_rotated_bvecs\" % subj\n",
"\n",
"t1_data = img.load_img(t1)\n",
"dwi_data = img.load_img(dwi)\n",
"\n",
"gt_bvals, gt_bvecs = read_bvals_bvecs(bval, bvec)\n",
"gtab = gradient_table(gt_bvals, gt_bvecs)"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Next, we need to create the tensor model using our gradient table and then fit the model using our data! We will start by creating a mask from our data and apply it to avoid calculating tensors on the background! This can be done using `DIPY`'s mask module. Then, we will our data!"
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"import dipy.reconst.dti as dti\n",
"from dipy.segment.mask import median_otsu\n",
"\n",
"dwi_data = dwi_data.get_fdata() # We re-use the variable for memory purposes\n",
"dwi_data, dwi_mask = median_otsu(dwi_data, vol_idx=[0], numpass=1) # Specify the volume index to the b0 volumes\n",
"\n",
"dti_model = dti.TensorModel(gtab)\n",
"dti_fit = dti_model.fit(dwi_data, mask=dwi_mask) # This step may take a while"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"The fit method creates a TensorFit object which contains the fitting parameters and other attributes of the model. A number of quantitative scalar metrics can be derived from the eigenvalues! In this tutorial, we will cover fractional anisotropy, mean diffusivity, axial diffusivity, and radial diffusivity. Each of these scalar, rotationally invariant metrics were calculated in the previous fitting step!\n",
"\n",
"### Fractional anisotropy (FA)\n",
"\n",
"Fractional anisotropy (FA) characterizes the degree to which the distribution of diffusion in an imaging voxel is directional. That is, whether there is relatively unrestricted diffusion in a particular direction.\n",
"\n",
"Mathematically, FA is defined as the normalized variance of the eigenvalues of the tensor:\n",
"\n",
"\n",
"\n",
"Values of FA vary between 0 and 1 (unitless). In the cases of perfect, isotropic diffusion, , the diffusion tensor is a sphere and FA = 0. If the first two eigenvalues are equal the tensor will be oblate or planar, whereas if the first eigenvalue is larger than the other two, it will have the mentioned ellipsoid shape: as diffusion progressively becomes more anisotropic, eigenvalues become more unequal, causing the tensor to be elongated, with FA approaching 1. Note that FA should be interpreted carefully. It may be an indication of the density of packing fibers in a voxel and the amount of myelin wrapped around those axons, but it is not always a measure of \"tissue integrity\".\n",
"\n",
"Let's take a look at what the FA map looks like! An FA map is a gray-scale image, where higher intensities reflect more anisotropic diffuse regions.\n",
"\n",
"_Note: we will have to first create the image from the array, making use of the reference anatomical_"
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"from nilearn import plotting as plot\n",
"import matplotlib.pyplot as plt # To enable plotting\n",
"%matplotlib inline\n",
"\n",
"fa_img = img.new_img_like(ref_niimg=t1_data, data=dti_fit.fa)\n",
"plot.plot_anat(fa_img, cut_coords=(0, -29, 20))"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Derived from partial volume effects in imaging voxels due to the presence of different tissues, noise in the measurements and numerical errors, the DTI model estimation may yield negative eigenvalues. Such *degenerate* case is not physically meaningful. These values are usually revealed as black or 0-valued pixels in FA maps.\n",
"\n",
"FA is a central value in dMRI: large FA values imply that the underlying fiber populations have a very coherent orientation, whereas lower FA values point to voxels containing multiple fiber crossings. Lowest FA values are indicative of non-white matter tissue in healthy brains (see, for example, Alexander et al.'s \"Diffusion Tensor Imaging of the Brain\". Neurotherapeutics 4, 316-329 (2007), and Jeurissen et al.'s \"Investigating the Prevalence of Complex Fiber Configurations in White Matter Tissue with Diffusion Magnetic Resonance Imaging\". Hum. Brain Mapp. 2012, 34(11) pp. 2747-2766)."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"### Mean diffusivity (MD)\n",
"\n",
"An often used complimentary measure to FA is mean diffusivity (MD). MD is a measure of the degree of diffusion, independent of direction. This is sometimes known as the apparent diffusion coefficient (ADC). Mathematically, MD is computed as the mean eigenvalues of the tensor and is measured in mm^2/s.\n",
"\n",
"\n",
"\n",
"Similar to the previous FA image, let's take a look at what the MD map looks like. Again, higher intensities reflect higher mean diffusivity!\n"
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"%matplotlib inline\n",
"\n",
"md_img = img.new_img_like(ref_niimg=t1_data, data=dti_fit.md)\n",
"# Arbitrarily set min and max of color bar\n",
"plot.plot_anat(md_img, cut_coords=(0, -29, 20), vmin=0, vmax=0.01)"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"### Axial and radial diffusivity (AD & RD)\n",
"\n",
"The final two metrics we will discuss are axial diffusivity (AD) and radial diffusivity (RD). Two tensors with different shapes may yield the same FA values, and additional measures such as AD and RD are required to further characterize the tensor. AD describes the diffusion rate along the primary axis of diffusion, along , or parallel to the axon (and hence, some works refer to it as the *parallel diffusivity*). On the other hand, RD reflects the average diffusivity along the other two minor axes (being named as *perpendicular diffusivity* in some works) (). Both are measured in mm^2/s.\n",
"\n",
""
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"### Tensor visualizations\n",
"\n",
"There are several ways of visualizing tensors. One way is using an RGB map, which overlays the primary diffusion orientation on an FA map. The colours of this map encodes the diffusion orientation. Note that this map provides no directional information (e.g. whether the diffusion flows from right-to-left or vice-versa). To do this with DIPY, we can use the color_fa function. The colours map to the following orientations:\n",
"\n",
"* Red = Left / Right\n",
"* Green = Anterior / Posterior\n",
"* Blue = Superior / Inferior\n",
"\n",
"_Note: The plotting functions in nilearn are unable to visualize these RGB maps. However, we can use the matplotlib library to view these images._"
],
"metadata": {}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"from scipy import ndimage # To rotate image for visualization purposes\n",
"from dipy.reconst.dti import color_fa\n",
"%matplotlib inline\n",
"\n",
"RGB_map = color_fa(dti_fit.fa, dti_fit.evecs)\n",
"\n",
"\n",
"fig, ax = plt.subplots(1, 3, figsize=(10, 10))\n",
"ax[0].imshow(ndimage.rotate(\n",
" RGB_map[:, RGB_map.shape[1]//2, :, :], 90, reshape=False))\n",
"ax[1].imshow(ndimage.rotate(\n",
" RGB_map[RGB_map.shape[0]//2, :, :, :], 90, reshape=False))\n",
"ax[2].imshow(ndimage.rotate(\n",
" RGB_map[:, :, RGB_map.shape[2]//2, :], 90, reshape=False))"
],
"outputs": [],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Another way of viewing the tensors is to visualize the diffusion tensor in each imaging voxel with colour encoding (we will refer you to the [`Dipy` documentation](https://dipy.org/tutorials/) for the steps to perform this type of visualization as it can be memory intensive). Below is an example image of such tensor visualization.\n",
"\n",
""
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"### Some notes on DTI\n",
"\n",
"DTI is only one of many models and is one of the simplest models available for modelling diffusion. While it is used for many studies, there are also some drawbacks (e.g. ability to distinguish multiple fibre orientations in an imaging voxel). Examples of this can be seen below!\n",
"\n",
"\n",
"\n",
"_Sourced from Sotiropoulos and Zalesky (2017). Building connectomes using diffusion MRI: why, how, and but. NMR in Biomedicine. 4(32). e3752. doi:10.1002/nbm.3752._\n",
"\n",
"Though other models are outside the scope of this lesson, we recommend looking into some of the pros and cons of each model (listed previously) to choose one best suited for your data!"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Exercise 1\n",
"\n",
"Plot the axial and radial diffusivity maps of the example given. Start from fitting the preprocessed diffusion image."
],
"metadata": {
"solution2": "hidden",
"solution2_first": true
}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"from bids.layout import BIDSLayout\n",
"from dipy.io.gradients import read_bvals_bvecs\n",
"from dipy.core.gradients import gradient_table\n",
"import dipy.reconst.dti as dti\n",
"from dipy.segment.mask import median_otsu\n",
"from nilearn import image as img\n",
"import nibabel as nib\n",
"\n",
"deriv_layout = BIDSLayout(\"../data/ds000221/derivatives\", validate=False)\n",
"subj = \"010006\"\n",
"\n",
"t1 = deriv_layout.get(subject=subj, space=\"dwi\",\n",
" extension='nii.gz', return_type='file')[0]\n",
"dwi = \"../data/ds000221/derivatives/uncorrected_topup_eddy/sub-%s/ses-01/dwi/dwi.nii.gz\" % subj\n",
"bval = \"../data/ds000221/sub-%s/ses-01/dwi/sub-%s_ses-01_dwi.bval\" % (\n",
" subj, subj)\n",
"bvec = \"../data/ds000221/derivatives/uncorrected_topup_eddy/sub-%s/ses-01/dwi/dwi.eddy_rotated_bvecs\" % subj\n",
"\n",
"t1_data = img.load_img(t1)\n",
"dwi_data = img.load_img(dwi)\n",
"\n",
"gt_bvals, gt_bvecs = read_bvals_bvecs(bval, bvec)\n",
"gtab = gradient_table(gt_bvals, gt_bvecs)\n",
"\n",
"dwi_data = dwi_data.get_fdata()\n",
"dwi_data, dwi_mask = median_otsu(dwi_data, vol_idx=[0], numpass=1)\n",
"\n",
"# Fit dti model\n",
"dti_model = dti.TensorModel(gtab)\n",
"dti_fit = dti_model.fit(dwi_data, mask=dwi_mask) # This step may take a while\n",
"\n",
"# Plot axial diffusivity map\n",
"ad_img = img.new_img_like(ref_niimg=t1_data, data=dti_fit.ad)\n",
"plot.plot_anat(ad_img, cut_coords=(0, -29, 20), vmin=0, vmax=0.01)\n",
"\n",
"# Plot radial diffusivity map\n",
"rd_img = img.new_img_like(ref_niimg=t1_data, data=dti_fit.rd)\n",
"plot.plot_anat(rd_img, cut_coords=(0, -29, 20), vmin=0, vmax=0.01)"
],
"outputs": [],
"metadata": {
"solution2": "hidden"
}
},
{
"cell_type": "code",
"execution_count": null,
"source": [
"# –––––––––––––--––– #\n",
"# –––– Solution –––– #\n",
"# ––––––––––––––--–– #\n"
],
"outputs": [],
"metadata": {
"tags": [
"nbval-skip"
]
}
}
],
"metadata": {
"celltoolbar": "Tags",
"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.7.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}