{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Estimating diffusion time dependent q-space indices using qt-dMRI\nEffective representation of the four-dimensional diffusion MRI signal --\nvarying over three-dimensional q-space and diffusion time -- is a sought-after\nand still unsolved challenge in diffusion MRI (dMRI). We propose a functional\nbasis approach that is specifically designed to represent the dMRI signal in\nthis qtau-space [Fick2017]_. Following recent terminology, we refer to our\nqtau-functional basis as $q au$-dMRI. We use GraphNet regularization --\nimposing both signal smoothness and sparsity -- to drastically reduce the\nnumber of diffusion-weighted images (DWIs) that is needed to represent the dMRI\nsignal in the qtau-space. As the main contribution, $q au$-dMRI provides\nthe framework to -- without making biophysical assumptions -- represent the\n$q au$-space signal and estimate time-dependent q-space indices\n($q au$-indices), providing a new means for studying diffusion in\nnervous tissue. $q au$-dMRI is the first of its kind in being\nspecifically designed to provide open interpretation of the\n$q au$-diffusion signal.\n\n$q au$-dMRI can be seen as a time-dependent extension of the MAP-MRI\nfunctional basis [Ozarslan2013]_, and all the previously proposed q-space\ncan be estimated for any diffusion time. These include rotationally\ninvariant quantities such as the Mean Squared Displacement (MSD), Q-space\nInverse Variance (QIV) and Return-To-Origin Probability (RTOP). Also\ndirectional indices such as the Return To the Axis Probability (RTAP) and\nReturn To the Plane Probability (RTPP) are available, as well as the\nOrientation Distribution Function (ODF).\n\nIn this example we illustrate how to use the $q au$-dMRI to estimate\ntime-dependent q-space indices from a $q au$-acquisition of a mouse.\n\nFirst import the necessary modules:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.data.fetcher import (fetch_qtdMRI_test_retest_2subjects,\n read_qtdMRI_test_retest_2subjects)\nfrom dipy.reconst import qtdmri, dti\nimport matplotlib.pyplot as plt\nimport numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Download and read the data for this tutorial.\n\n$q\\tau$-dMRI requires data with multiple gradient directions, gradient\nstrength and diffusion times. We will use the test-retest acquisitions of two\nmice that were used in the test-retest study by [Fick2017]_. The data itself\nis freely available and citeable at [Wassermann2017]_.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fetch_qtdMRI_test_retest_2subjects()\ndata, cc_masks, gtabs = read_qtdMRI_test_retest_2subjects()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "data contains 4 qt-dMRI datasets of size [80, 160, 5, 515]. The first two are\nthe test-retest datasets of the first mouse and the second two are those of the\nsecond mouse. cc_masks contains 4 corresponding binary masks for the corpus\ncallosum voxels in the middle slice that were used in the test-retest study.\nFinally, gtab contains the qt-dMRI gradient tables for the DWIs in the dataset.\n\nThe data consists of 515 DWIs, divided over 35 shells, with 7 \"gradient\nstrength shells\" up to 491 mT/m, 5 equally spaced \"pulse separation shells\"\n(big_delta) between [10.8-20] ms and a pulse duration (small_delta) of 5ms.\n\nTo visualize qt-dMRI acquisition schemes in an intuitive way, the qtdmri module\nprovides a visualization function to illustrate the relationship between\ngradient strength (G), pulse separation (big_delta) and b-value:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure()\nqtdmri.visualise_gradient_table_G_Delta_rainbow(gtabs[0])\nplt.savefig('qt-dMRI_acquisition_scheme.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: qt-dMRI_acquisition_scheme.png\n :align: center\n\nIn the figure the dots represent measured DWIs in any direction, for a given\ngradient strength and pulse separation. The background isolines represent the\ncorresponding b-values for different combinations of G and big_delta.\n\nNext, we visualize the middle slices of the test-retest data sets with their\ncorresponding masks. To better illustrate the white matter architecture in the\ndata, we calculate DTI's fractional anisotropy (FA) over the whole slice and\nproject the corpus callosum mask on the FA image.:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "subplot_titles = [\"Subject1 Test\", \"Subject1 Retest\",\n \"Subject2 Test\", \"Subject2 Retest\"]\nfig = plt.figure()\nplt.subplots(nrows=2, ncols=2)\nfor i, (data_, mask_, gtab_) in enumerate(zip(data, cc_masks, gtabs)):\n # take the middle slice\n data_middle_slice = data_[:, :, 2]\n mask_middle_slice = mask_[:, :, 2]\n\n # estimate fractional anisotropy (FA) for this slice\n tenmod = dti.TensorModel(gtab_)\n tenfit = tenmod.fit(data_middle_slice, data_middle_slice[..., 0] > 0)\n fa = tenfit.fa\n\n # set mask color to green with 0.5 opacity as overlay\n mask_template = np.zeros(np.r_[mask_middle_slice.shape, 4])\n mask_template[mask_middle_slice == 1] = np.r_[0., 1., 0., .5]\n\n # produce the FA images with corpus callosum masks.\n plt.subplot(2, 2, 1 + i)\n plt.title(subplot_titles[i], fontsize=15)\n plt.imshow(fa, cmap='Greys_r', origin='lower', interpolation='nearest')\n plt.imshow(mask_template, origin='lower', interpolation='nearest')\n plt.axis('off')\nplt.tight_layout()\nplt.savefig('qt-dMRI_datasets_fa_with_ccmasks.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: qt-dMRI_datasets_fa_with_ccmasks.png\n :align: center\n\nNext, we use qt-dMRI to estimate of time-dependent q-space indices\n(q$\\tau$-indices) for the masked voxels in the corpus callosum of each dataset.\nIn particular, we estimate the Return-to-Original, Return-to-Axis and\nReturn-to-Plane Probability (RTOP, RTAP and RTPP), as well as the Mean Squared\nDisplacement (MSD).\n\nIn this example we don't extrapolate the data beyond the maximum diffusion\ntime, so we estimate $q\\tau$ indices between the minimum and maximum\ndiffusion times of the data at 5 equally spaced points. However, it should the\nnoted that qt-dMRI's combined smoothness and sparsity regularization allows\nfor smooth interpolation at any $q\\tau$ position. In other words, once\nthe basis is fitted to the data, its coefficients describe the the entire\n$q\\tau$-space, and any $q\\tau$-position can be freely recovered.\nThis including points beyond the dataset's maximum $q\\tau$ value\n(although this should be done with caution).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "tau_min = gtabs[0].tau.min()\ntau_max = gtabs[0].tau.max()\ntaus = np.linspace(tau_min, tau_max, 5)\n\nqtdmri_fits = []\nmsds = []\nrtops = []\nrtaps = []\nrtpps = []\nfor i, (data_, mask_, gtab_) in enumerate(zip(data, cc_masks, gtabs)):\n # select the corpus callosum voxel for every dataset\n cc_voxels = data_[mask_ == 1]\n # initialize the qt-dMRI model.\n # recommended basis orders are radial_order=6 and time_order=2.\n # The combined Laplacian and l1-regularization using Generalized\n # Cross-Validation (GCV) and Cross-Validation (CV) settings is most robust,\n # but can be used separately and with weightings preset to any positive\n # value to optimize for speed.\n qtdmri_mod = qtdmri.QtdmriModel(\n gtab_, radial_order=6, time_order=2,\n laplacian_regularization=True, laplacian_weighting='GCV',\n l1_regularization=True, l1_weighting='CV'\n )\n # fit the model.\n # Here we take every 5th voxel for speed, but of course all voxels can be\n # fit for a more robust result later on.\n qtdmri_fit = qtdmri_mod.fit(cc_voxels[::5])\n qtdmri_fits.append(qtdmri_fit)\n # We estimate MSD, RTOP, RTAP and RTPP for the chosen diffusion times.\n msds.append(np.array(list(map(qtdmri_fit.msd, taus))))\n rtops.append(np.array(list(map(qtdmri_fit.rtop, taus))))\n rtaps.append(np.array(list(map(qtdmri_fit.rtap, taus))))\n rtpps.append(np.array(list(map(qtdmri_fit.rtpp, taus))))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The estimated $q\\tau$-indices, for the chosen diffusion times, are now\nstored in msds, rtops, rtaps and rtpps. The trends of these\n$q\\tau$-indices over time say something about the restriction of\ndiffusing particles over time, which is currently a hot topic in the dMRI\ncommunity. We evaluate the test-retest reproducibility for the two subjects by\nplotting the $q\\tau$-indices for each subject together. This example\nwill produce similar results as Fig. 10 in [Fick2017]_.\n\nWe first define a small function to plot the mean and standard deviation of the\n$q\\tau$-index trends in a subject.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def plot_mean_with_std(ax, time, ind1, plotcolor, ls='-', std_mult=1,\n label=''):\n means = np.mean(ind1, axis=1)\n stds = np.std(ind1, axis=1)\n ax.plot(time, means, c=plotcolor, lw=3, label=label, ls=ls)\n ax.fill_between(time,\n means + std_mult * stds,\n means - std_mult * stds,\n alpha=0.15, color=plotcolor)\n ax.plot(time, means + std_mult * stds, alpha=0.25, color=plotcolor)\n ax.plot(time, means - std_mult * stds, alpha=0.25, color=plotcolor)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start by showing the test-retest MSD of both subjects over time. We plot the\n$q\\tau$-indices together with $q\\tau$-index trends of free\ndiffusion with different diffusivities as background.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# we first generate the data to produce the background index isolines.\nDelta_ = np.linspace(0.005, 0.02, 100)\nMSD_ = np.linspace(4e-5, 10e-5, 100)\nDelta_grid, MSD_grid = np.meshgrid(Delta_, MSD_)\nD_grid = MSD_grid / (6 * Delta_grid)\nD_levels = np.r_[1, 5, 7, 10, 14, 23, 30] * 1e-4\n\nfig = plt.figure(figsize=(10, 3))\n# start with the plot of subject 1.\nax = plt.subplot(1, 2, 1)\n# first plot the background\nplt.contourf(Delta_ * 1e3, 1e5 * MSD_, D_grid, levels=D_levels, cmap='Greys',\n alpha=.5)\n\n# plot the test-retest mean MSD and standard deviation of subject 1.\nplot_mean_with_std(ax, taus * 1e3, 1e5 * msds[0], 'r', 'dashdot',\n label='MSD Test')\nplot_mean_with_std(ax, taus * 1e3, 1e5 * msds[1], 'g', 'dashdot',\n label='MSD Retest')\nax.legend(fontsize=13)\n# plot some text markers to clarify the background diffusivity lines.\nax.text(.0091 * 1e3, 6.33, 'D=14e-4', fontsize=12, rotation=35)\nax.text(.0091 * 1e3, 4.55, 'D=10e-4', fontsize=12, rotation=25)\nax.set_ylim(4, 9.5)\nax.set_xlim(.009 * 1e3, 0.0185 * 1e3)\nax.set_title(r'Test-Retest MSD($\\tau$) Subject 1', fontsize=15)\nax.set_xlabel('Diffusion Time (ms)', fontsize=17)\nax.set_ylabel('MSD ($10^{-5}mm^2$)', fontsize=17)\n\n# then do the same thing for subject 2.\nax = plt.subplot(1, 2, 2)\nplt.contourf(Delta_ * 1e3, 1e5 * MSD_, D_grid, levels=D_levels, cmap='Greys',\n alpha=.5)\ncb = plt.colorbar()\ncb.set_label('Free Diffusivity ($mm^2/s$)', fontsize=18)\n\nplot_mean_with_std(ax, taus * 1e3, 1e5 * msds[2], 'r', 'dashdot')\nplot_mean_with_std(ax, taus * 1e3, 1e5 * msds[3], 'g', 'dashdot')\nax.set_ylim(4, 9.5)\nax.set_xlim(.009 * 1e3, 0.0185 * 1e3)\nax.set_xlabel('Diffusion Time (ms)', fontsize=17)\nax.set_title(r'Test-Retest MSD($\\tau$) Subject 2', fontsize=15)\nplt.savefig('qt_indices_msd.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: qt_indices_msd.png\n :align: center\n\nYou can see that the MSD in both subjects increases over time, but also slowly\nlevels off as time progresses. This makes sense as diffusing particles are\nbecoming more restricted by surrounding tissue as time goes on. You can also\nsee that for Subject 1 the index trends nearly perfectly overlap, but for\nsubject 2 they are slightly off, which is also what we found in the paper.\n\nNext, we follow the same procedure to estimate the test-retest RTAP, RTOP and\nRTPP over diffusion time for both subject. For ease of comparison, we will\nestimate all three in the same unit [1/mm] by taking the square root of RTAP\nand the cubed root of RTOP.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Again, first we define the data for the background illustration.\nDelta_ = np.linspace(0.005, 0.02, 100)\nRTXP_ = np.linspace(1, 200, 100)\nDelta_grid, RTXP_grid = np.meshgrid(Delta_, RTXP_)\nD_grid = 1 / (4 * RTXP_grid ** 2 * np.pi * Delta_grid)\nD_levels = np.r_[1, 2, 3, 4, 6, 9, 15, 30] * 1e-4\nD_colors = np.tile(np.linspace(.8, 0, 7), (3, 1)).T\n\n# We start with estimating the RTOP illustration.\nfig = plt.figure(figsize=(10, 3))\nax = plt.subplot(1, 2, 1)\nplt.contourf(Delta_ * 1e3, RTXP_, D_grid, colors=D_colors, levels=D_levels,\n alpha=.5)\n\nplot_mean_with_std(ax, taus * 1e3, rtops[0] ** (1 / 3.), 'r', '--',\n label='RTOP$^{1/3}$ Test')\nplot_mean_with_std(ax, taus * 1e3, rtops[1] ** (1 / 3.), 'g', '--',\n label='RTOP$^{1/3}$ Retest')\nax.legend(fontsize=13)\nax.text(.0091 * 1e3, 162, 'D=3e-4', fontsize=12, rotation=-22)\nax.text(.0091 * 1e3, 140, 'D=4e-4', fontsize=12, rotation=-20)\nax.text(.0091 * 1e3, 113, 'D=6e-4', fontsize=12, rotation=-16)\nax.set_ylim(54, 170)\nax.set_xlim(.009 * 1e3, 0.0185 * 1e3)\nax.set_title(r'Test-Retest RTOP($\\tau$) Subject 1', fontsize=15)\nax.set_xlabel('Diffusion Time (ms)', fontsize=17)\nax.set_ylabel('RTOP$^{1/3}$ (1/mm)', fontsize=17)\n\nax = plt.subplot(1, 2, 2)\nplt.contourf(Delta_ * 1e3, RTXP_, D_grid, colors=D_colors, levels=D_levels,\n alpha=.5)\ncb = plt.colorbar()\ncb.set_label('Free Diffusivity ($mm^2/s$)', fontsize=18)\n\nplot_mean_with_std(ax, taus * 1e3, rtops[2] ** (1 / 3.), 'r', '--')\nplot_mean_with_std(ax, taus * 1e3, rtops[3] ** (1 / 3.), 'g', '--')\nax.set_ylim(54, 170)\nax.set_xlim(.009 * 1e3, 0.0185 * 1e3)\nax.set_xlabel('Diffusion Time (ms)', fontsize=17)\nax.set_title(r'Test-Retest RTOP($\\tau$) Subject 2', fontsize=15)\nplt.savefig('qt_indices_rtop.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: qt_indices_rtop.png\n :align: center\n\nSimilarly as MSD, the RTOP is related to the restriction that particles are\nexperiencing and is also rotationally invariant. RTOP is defined as the\nprobability that particles are found at the same position at the time of both\ngradient pulses. As time increases, the odds become smaller that a particle\nwill arrive at the same position it left, which is illustrated by all RTOP\ntrends in the figure. Notice that the estimated RTOP trends decrease less fast\nthan free diffusion, meaning that particles experience restriction over time.\nAlso notice that the RTOP trends in both subjects nearly perfectly overlap.\n\nNext, we estimate two directional $q\\tau$-indices, RTAP and RTPP,\ndescribing particle restriction perpendicular and parallel to the orientation\nof the principal diffusivity in that voxel. If the voxel describes coherent\nwhite matter (which it does in our corpus callosum example), then they describe\nproperties related to restriction perpendicular and parallel to the axon\nbundles.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# First, we estimate the RTAP trends.\nfig = plt.figure(figsize=(10, 3))\nax = plt.subplot(1, 2, 1)\nplt.contourf(Delta_ * 1e3, RTXP_, D_grid, colors=D_colors, levels=D_levels,\n alpha=.5)\n\nplot_mean_with_std(ax, taus * 1e3, np.sqrt(rtaps[0]), 'r', '-',\n label='RTAP$^{1/2}$ Test')\nplot_mean_with_std(ax, taus * 1e3, np.sqrt(rtaps[1]), 'g', '-',\n label='RTAP$^{1/2}$ Retest')\nax.legend(fontsize=13)\nax.text(.0091 * 1e3, 162, 'D=3e-4', fontsize=12, rotation=-22)\nax.text(.0091 * 1e3, 140, 'D=4e-4', fontsize=12, rotation=-20)\nax.text(.0091 * 1e3, 113, 'D=6e-4', fontsize=12, rotation=-16)\nax.set_ylim(54, 170)\nax.set_xlim(.009 * 1e3, 0.0185 * 1e3)\nax.set_title(r'Test-Retest RTAP($\\tau$) Subject 1', fontsize=15)\nax.set_xlabel('Diffusion Time (ms)', fontsize=17)\nax.set_ylabel('RTAP$^{1/2}$ (1/mm)', fontsize=17)\n\nax = plt.subplot(1, 2, 2)\nplt.contourf(Delta_ * 1e3, RTXP_, D_grid, colors=D_colors, levels=D_levels,\n alpha=.5)\ncb = plt.colorbar()\ncb.set_label('Free Diffusivity ($mm^2/s$)', fontsize=18)\n\nplot_mean_with_std(ax, taus * 1e3, np.sqrt(rtaps[2]), 'r', '-')\nplot_mean_with_std(ax, taus * 1e3, np.sqrt(rtaps[3]), 'g', '-')\nax.set_ylim(54, 170)\nax.set_xlim(.009 * 1e3, 0.0185 * 1e3)\nax.set_xlabel('Diffusion Time (ms)', fontsize=17)\nax.set_title(r'Test-Retest RTAP($\\tau$) Subject 2', fontsize=15)\nplt.savefig('qt_indices_rtap.png')\n\n\n# Finally the last one for RTPP.\nfig = plt.figure(figsize=(10, 3))\nax = plt.subplot(1, 2, 1)\nplt.contourf(Delta_ * 1e3, RTXP_, D_grid, colors=D_colors, levels=D_levels,\n alpha=.5)\n\nplot_mean_with_std(ax, taus * 1e3, rtpps[0], 'r', ':', label='RTPP Test')\nplot_mean_with_std(ax, taus * 1e3, rtpps[1], 'g', ':', label='RTPP Retest')\nax.legend(fontsize=13)\nax.text(.0091 * 1e3, 113, 'D=6e-4', fontsize=12, rotation=-16)\nax.text(.0091 * 1e3, 91, 'D=9e-4', fontsize=12, rotation=-13)\nax.text(.0091 * 1e3, 69, 'D=15e-4', fontsize=12, rotation=-10)\nax.set_ylim(54, 170)\nax.set_xlim(.009 * 1e3, 0.0185 * 1e3)\nax.set_title(r'Test-Retest RTPP($\\tau$) Subject 1', fontsize=15)\nax.set_xlabel('Diffusion Time (ms)', fontsize=17)\nax.set_ylabel('RTPP (1/mm)', fontsize=17)\n\nax = plt.subplot(1, 2, 2)\nplt.contourf(Delta_ * 1e3, RTXP_, D_grid, colors=D_colors, levels=D_levels,\n alpha=.5)\ncb = plt.colorbar()\ncb.set_label('Free Diffusivity ($mm^2/s$)', fontsize=18)\n\nplot_mean_with_std(ax, taus * 1e3, rtpps[2], 'r', ':')\nplot_mean_with_std(ax, taus * 1e3, rtpps[3], 'g', ':')\nax.set_ylim(54, 170)\nax.set_xlim(.009 * 1e3, 0.0185 * 1e3)\nax.set_xlabel('Diffusion Time (ms)', fontsize=17)\nax.set_title(r'Test-Retest RTPP($\\tau$) Subject 2', fontsize=15)\nplt.savefig('qt_indices_rtpp.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: qt_indices_rtap.png\n :align: center\n.. figure:: qt_indices_rtpp.png\n :align: center\n\nAs those of RTOP, the trends in RTAP and RTPP also decrease over time. It can\nbe seen that RTAP$^{1/2}$ is always bigger than RTPP, which makes sense as\nparticles in coherent white matter experience more restriction perpendicular to\nthe white matter orientation than parallel to it. Again, in both subjects the\ntest-retest RTAP and RTPP is nearly perfectly consistent.\n\nAside from the estimation of $q\\tau$-space indices, $q\\tau$-dMRI\nalso allows for the estimation of time-dependent ODFs. Once the Qtdmri model\nis fitted it can be simply called by qtdmri_fit.odf(sphere,\ns=sharpening_factor). This is identical to how the mapmri module functions,\nand allows to study the time-dependence of ODF directionality.\n\nThis concludes the example on qt-dMRI. As we showed, approaches such as qt-dMRI\ncan help in studying the (finite-$\\tau$) temporal properties of diffusion\nin biological tissues. Differences in $q\\tau$-index trends could be\nindicative of underlying structural differences that affect the time-dependence\nof the diffusion process.\n\n.. [Fick2017]_ Fick, Rutger HJ, et al. \"Non-Parametric GraphNet-Regularized\n Representation of dMRI in Space and Time\", Medical Image Analysis,\n 2017.\n.. [Wassermann2017]_ Wassermann, Demian, et al. \"Test-Retest qt-dMRI datasets\n for 'Non-Parametric GraphNet-Regularized Representation of dMRI in\n Space and Time' [Data set]\". Zenodo.\n https://doi.org/10.5281/zenodo.996889, 2017.\n\n" ] } ], "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.9.13" } }, "nbformat": 4, "nbformat_minor": 0 }