Using the free water elimination model to remove DTI free water contamination

As shown previously (see example_reconst_dti), the diffusion tensor model is a simple way to characterize diffusion anisotropy. However, in regions near the ventricles and parenchyma, anisotropy can be underestimated by partial volume effects of the cerebral spinal fluid (CSF). This free water contamination can particularly corrupt Diffusion Tensor Imaging analysis of microstructural changes when different groups of subjects show different brain morphology (e.g. brain ventricle enlargement associated with brain tissue atrophy that occurs in several brain pathologies and aging).

A way to remove this free water influences is to expand the DTI model to take into account an extra compartment representing the contributions of free water diffusion [Pasternak2009]. The expression of the expanded DTI model is shown below:

where \(\mathbf{g}\) and \(b\) are diffusion gradient direction and weighted (more information see example_reconst_dti), \(S(\mathbf{g}, b)\) is the diffusion-weighted signal measured, \(S_0\) is the signal in a measurement with no diffusion weighting, \(\mathbf{D}\) is the diffusion tensor, \(f\) the volume fraction of the free water component, and \(D_{iso}\) is the isotropic value of the free water diffusion (normally set to \(3.0 imes 10^{-3} mm^{2}s^{-1}\)).

In this example, we show how to process a diffusion weighting dataset using an adapted version of the free water elimination proposed by [Hoy2014].

The full details of Dipy’s free water DTI implementation was published in [Henriques2017]. Please cite this work if you use this algorithm.

Let’s start by importing the relevant modules:

import numpy as np
import dipy.reconst.fwdti as fwdti
import dipy.reconst.dti as dti
import matplotlib.pyplot as plt
from dipy.data import fetch_hbn
import os.path as op
import nibabel as nib
from dipy.core.gradients import gradient_table

Without spatial constrains the free water elimination model cannot be solved in data acquired from one non-zero b-value [Hoy2014]. Therefore, here we download a dataset that was acquired with multiple b-values.

data_path = fetch_hbn(["NDARAA948VFH"])[1]
dwi_path = op.join(
       data_path, "derivatives", "qsiprep", "sub-NDARAA948VFH",
       "ses-HBNsiteRU", "dwi")

img = nib.load(op.join(
       dwi_path,
       "sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi.nii.gz"))

gtab = gradient_table(
       op.join(dwi_path,
"sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi.bval"),
       op.join(dwi_path,
"sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi.bvec"))

data = np.asarray(img.dataobj)

The free water DTI model can take some minutes to process the full data set. Thus, we use a brain mask that was calculated during pre-processing, to remove the background of the image and avoid unnecessary calculations.

mask_img = nib.load(
       op.join(dwi_path,
"sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-brain_mask.nii.gz"))

Moreover, for illustration purposes we process only one slice of the data.

mask = mask_img.get_fdata()

data_small = data[:, :, 50:51]
mask_small = mask[:, :, 50:51]

The free water elimination model fit can then be initialized by instantiating a FreeWaterTensorModel class object:

fwdtimodel = fwdti.FreeWaterTensorModel(gtab)

The data can then be fitted using the fit function of the defined model object:

fwdtifit = fwdtimodel.fit(data_small, mask=mask_small)

  0%|          | 0/5148.0 [00:00<?, ?it/s]
  1%|1         | 76/5148.0 [00:00<00:06, 752.96it/s]
  3%|3         | 161/5148.0 [00:00<00:06, 805.77it/s]
  5%|4         | 242/5148.0 [00:00<00:06, 790.45it/s]
  6%|6         | 322/5148.0 [00:00<00:06, 775.86it/s]
  8%|7         | 405/5148.0 [00:00<00:05, 794.43it/s]
 10%|9         | 493/5148.0 [00:00<00:05, 821.23it/s]
 11%|#1        | 579/5148.0 [00:00<00:05, 831.67it/s]
 13%|#2        | 666/5148.0 [00:00<00:05, 843.53it/s]
 15%|#4        | 751/5148.0 [00:00<00:05, 842.83it/s]
 16%|#6        | 839/5148.0 [00:01<00:05, 852.91it/s]
 18%|#7        | 925/5148.0 [00:01<00:04, 853.56it/s]
 20%|#9        | 1011/5148.0 [00:01<00:04, 853.76it/s]
 21%|##1       | 1097/5148.0 [00:01<00:04, 850.94it/s]
 23%|##2       | 1183/5148.0 [00:01<00:04, 842.55it/s]/opt/homebrew/Caskroom/miniforge/base/envs/dipy-39-x86/lib/python3.9/site-packages/scipy/optimize/_minpack_py.py:492: RuntimeWarning: Number of calls to function has reached maxfev = 1800.
  warnings.warn(errors[info][0], RuntimeWarning)

 25%|##4       | 1268/5148.0 [00:01<00:05, 759.44it/s]
 26%|##6       | 1346/5148.0 [00:01<00:05, 711.05it/s]
 28%|##7       | 1419/5148.0 [00:01<00:05, 692.93it/s]
 29%|##9       | 1503/5148.0 [00:01<00:04, 731.21it/s]
 31%|###       | 1589/5148.0 [00:02<00:04, 765.05it/s]
 32%|###2      | 1671/5148.0 [00:02<00:04, 779.87it/s]
 34%|###4      | 1757/5148.0 [00:02<00:04, 801.64it/s]
 36%|###5      | 1838/5148.0 [00:02<00:04, 751.72it/s]
 37%|###7      | 1923/5148.0 [00:02<00:04, 778.72it/s]
 39%|###9      | 2010/5148.0 [00:02<00:03, 802.31it/s]
 41%|####      | 2097/5148.0 [00:02<00:03, 819.46it/s]
 42%|####2     | 2180/5148.0 [00:02<00:03, 819.16it/s]
 44%|####3     | 2264/5148.0 [00:02<00:03, 824.85it/s]
 46%|####5     | 2347/5148.0 [00:02<00:03, 823.75it/s]
 47%|####7     | 2430/5148.0 [00:03<00:03, 825.47it/s]
 49%|####8     | 2513/5148.0 [00:03<00:03, 821.90it/s]
 50%|#####     | 2596/5148.0 [00:03<00:03, 706.51it/s]
 52%|#####1    | 2670/5148.0 [00:03<00:03, 639.60it/s]
 54%|#####3    | 2757/5148.0 [00:03<00:03, 696.39it/s]
 55%|#####5    | 2840/5148.0 [00:03<00:03, 731.79it/s]
 57%|#####6    | 2922/5148.0 [00:03<00:02, 751.63it/s]
 58%|#####8    | 3005/5148.0 [00:03<00:02, 771.12it/s]
 60%|#####9    | 3084/5148.0 [00:03<00:02, 774.30it/s]
 62%|######1   | 3169/5148.0 [00:04<00:02, 794.99it/s]
 63%|######3   | 3254/5148.0 [00:04<00:02, 808.65it/s]
 65%|######4   | 3339/5148.0 [00:04<00:02, 820.50it/s]
 67%|######6   | 3425/5148.0 [00:04<00:02, 831.52it/s]
 68%|######8   | 3509/5148.0 [00:04<00:01, 832.64it/s]
 70%|######9   | 3593/5148.0 [00:04<00:01, 833.14it/s]
 71%|#######1  | 3677/5148.0 [00:04<00:01, 833.01it/s]
 73%|#######3  | 3761/5148.0 [00:04<00:01, 754.61it/s]
 75%|#######4  | 3838/5148.0 [00:04<00:01, 704.50it/s]
 76%|#######5  | 3910/5148.0 [00:05<00:01, 686.72it/s]
 77%|#######7  | 3980/5148.0 [00:05<00:01, 671.27it/s]
 79%|#######8  | 4064/5148.0 [00:05<00:01, 715.48it/s]
 80%|########  | 4142/5148.0 [00:05<00:01, 733.54it/s]
 82%|########2 | 4227/5148.0 [00:05<00:01, 763.40it/s]
 84%|########3 | 4312/5148.0 [00:05<00:01, 785.73it/s]
 85%|########5 | 4394/5148.0 [00:05<00:00, 794.26it/s]
 87%|########7 | 4482/5148.0 [00:05<00:00, 818.89it/s]
 89%|########8 | 4567/5148.0 [00:05<00:00, 824.37it/s]
 90%|######### | 4652/5148.0 [00:05<00:00, 830.22it/s]
 92%|#########1| 4736/5148.0 [00:06<00:00, 811.69it/s]
 94%|#########3| 4818/5148.0 [00:06<00:00, 795.24it/s]
 95%|#########5| 4902/5148.0 [00:06<00:00, 805.89it/s]
 97%|#########6| 4986/5148.0 [00:06<00:00, 813.48it/s]
 99%|#########8| 5072/5148.0 [00:06<00:00, 824.82it/s]
100%|##########| 5148/5148.0 [00:06<00:00, 784.36it/s]

This 2-steps procedure will create a FreeWaterTensorFit object which contains all the diffusion tensor statistics free for free water contamination. Below we extract the fractional anisotropy (FA) and the mean diffusivity (MD) of the free water diffusion tensor.

FA = fwdtifit.fa
MD = fwdtifit.md

For comparison we also compute the same standard measures processed by the standard DTI model

dtimodel = dti.TensorModel(gtab)

dtifit = dtimodel.fit(data_small, mask=mask_small)

dti_FA = dtifit.fa
dti_MD = dtifit.md

Below the FA values for both free water elimination DTI model and standard DTI model are plotted in panels A and B, while the respective MD values are plotted in panels D and E. For a better visualization of the effect of the free water correction, the differences between these two metrics are shown in panels C and E. In addition to the standard diffusion statistics, the estimated volume fraction of the free water contamination is shown on panel G.

fig1, ax = plt.subplots(2, 4, figsize=(12, 6),
                        subplot_kw={'xticks': [], 'yticks': []})

fig1.subplots_adjust(hspace=0.3, wspace=0.05)
ax.flat[0].imshow(FA[:, :, 0].T, origin='lower',
                  cmap='gray', vmin=0, vmax=1)
ax.flat[0].set_title('A) fwDTI FA')
ax.flat[1].imshow(dti_FA[:, :, 0].T, origin='lower',
                  cmap='gray', vmin=0, vmax=1)
ax.flat[1].set_title('B) standard DTI FA')

FAdiff = abs(FA[:, :, 0] - dti_FA[:, :, 0])
ax.flat[2].imshow(FAdiff.T, cmap='gray', origin='lower', vmin=0, vmax=1)
ax.flat[2].set_title('C) FA difference')

ax.flat[3].axis('off')

ax.flat[4].imshow(MD[:, :, 0].T, origin='lower',
                  cmap='gray', vmin=0, vmax=2.5e-3)
ax.flat[4].set_title('D) fwDTI MD')
ax.flat[5].imshow(dti_MD[:, :, 0].T, origin='lower',
                  cmap='gray', vmin=0, vmax=2.5e-3)
ax.flat[5].set_title('E) standard DTI MD')

MDdiff = abs(MD[:, :, 0] - dti_MD[:, :, 0])
ax.flat[6].imshow(MDdiff.T, origin='lower', cmap='gray', vmin=0, vmax=2.5e-3)
ax.flat[6].set_title('F) MD difference')

F = fwdtifit.f

ax.flat[7].imshow(F[:, :, 0].T, origin='lower',
                  cmap='gray', vmin=0, vmax=1)
ax.flat[7].set_title('G) free water volume')

plt.show()
fig1.savefig('In_vivo_free_water_DTI_and_standard_DTI_measures.png')

From the figure, one can observe that the free water elimination model produces in general higher values of FA and lower values of MD than the standard DTI model. These differences in FA and MD estimation are expected due to the suppression of the free water isotropic diffusion components. Unexpected high amplitudes of FA are however observed in the periventricular gray matter. This is a known artefact of regions associated to voxels with high water volume fraction (i.e. voxels containing basically CSF). We are able to remove this problematic voxels by excluding all FA values associated with measured volume fractions above a reasonable threshold of 0.7:

FA[F > 0.7] = 0
dti_FA[F > 0.7] = 0

Above we reproduce the plots of the in vivo FA from the two DTI fits and where we can see that the inflated FA values were practically removed:

fig1, ax = plt.subplots(1, 3, figsize=(9, 3),
                        subplot_kw={'xticks': [], 'yticks': []})

fig1.subplots_adjust(hspace=0.3, wspace=0.05)
ax.flat[0].imshow(FA[:, :, 0].T, origin='lower',
                  cmap='gray', vmin=0, vmax=1)
ax.flat[0].set_title('A) fwDTI FA')
ax.flat[1].imshow(dti_FA[:, :, 0].T, origin='lower',
                  cmap='gray', vmin=0, vmax=1)
ax.flat[1].set_title('B) standard DTI FA')

FAdiff = abs(FA[:, :, 0] - dti_FA[:, :, 0])
ax.flat[2].imshow(FAdiff.T, cmap='gray', origin='lower', vmin=0, vmax=1)
ax.flat[2].set_title('C) FA difference')

plt.show()
fig1.savefig('In_vivo_free_water_DTI_and_standard_DTI_corrected.png')