{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Mean signal diffusion kurtosis imaging (MSDKI)\n\nDiffusion Kurtosis Imaging (DKI) is one of the conventional ways to estimate\nthe degree of non-Gaussian diffusion (see `example_reconst_dki`)\n[Jensen2005]_. However, a limitation of DKI is that its measures are highly\nsensitive to noise and image artefacts. For instance, due to the low radial\ndiffusivities, standard kurtosis estimates in regions of well-aligned voxel may\nbe corrupted by implausible low or even negative values.\n\nA way to overcome this issue is to characterize kurtosis from average signals\nacross all directions acquired for each data b-value (also known as\npowder-averaged signals). Moreover, as previously pointed [NetoHe2015]_,\nstandard kurtosis measures (e.g. radial, axial and standard mean kurtosis)\ndo not only depend on microstructural properties but also on mesoscopic\nproperties such as fiber dispersion or the intersection angle of crossing\nfibers. In contrary, the kurtosis from powder-average signals has the advantage\nof not depending on the fiber distribution functions [NetoHe2018]_,\n[NetoHe2019]_.\n\nIn short, in this tutorial we show how to characterize non-Gaussian diffusion\nin a more precise way and decoupled from confounding effects of tissue\ndispersion and crossing.\n\nIn the first part of this example, we illustrate the properties of the measures\nobtained from the mean signal diffusion kurtosis imaging (MSDKI)[NetoHe2018]_\nusing synthetic data. Secondly, the mean signal diffusion kurtosis imaging will\nbe applied to in-vivo MRI data. Finally, we show how MSDKI provides the same\ninformation than common microstructural models such as the spherical mean\ntechnique [NetoHe2019]_, [Kaden2016b]_.\n\nLet's import all relevant modules:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport matplotlib.pyplot as plt\n\n# Reconstruction modules\nimport dipy.reconst.dki as dki\nimport dipy.reconst.msdki as msdki\n\n# For simulations\nfrom dipy.sims.voxel import multi_tensor\nfrom dipy.core.gradients import gradient_table\nfrom dipy.core.sphere import disperse_charges, HemiSphere\n\n# For in-vivo data\nfrom dipy.data import get_fnames\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.io.image import load_nifti\nfrom dipy.segment.mask import median_otsu" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Testing MSDKI in synthetic data\n\nWe simulate representative diffusion-weighted signals using MultiTensor\nsimulations (for more information on this type of simulations see\n`example_simulate_multi_tensor`). For this example, simulations are\nproduced based on the sum of four diffusion tensors to represent the intra-\nand extra-cellular spaces of two fiber populations. The parameters of these\ntensors are adjusted according to [NetoHe2015]_ (see also\n`example_simulate_dki`).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mevals = np.array([[0.00099, 0, 0],\n [0.00226, 0.00087, 0.00087],\n [0.00099, 0, 0],\n [0.00226, 0.00087, 0.00087]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the acquisition parameters of the synthetic data, we use 60 gradient\ndirections for two non-zero b-values (1000 and 2000 $s/mm^{2}$) and two\nzero bvalues (note that, such as the standard DKI, MSDKI requires at least\nthree different b-values).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Sample the spherical coordinates of 60 random diffusion-weighted directions.\nn_pts = 60\ntheta = np.pi * np.random.rand(n_pts)\nphi = 2 * np.pi * np.random.rand(n_pts)\n\n# Convert direction to cartesian coordinates.\nhsph_initial = HemiSphere(theta=theta, phi=phi)\n\n# Evenly distribute the 60 directions\nhsph_updated, potential = disperse_charges(hsph_initial, 5000)\ndirections = hsph_updated.vertices\n\n# Reconstruct acquisition parameters for 2 non-zero=b-values and 2 b0s\nbvals = np.hstack((np.zeros(2), 1000 * np.ones(n_pts), 2000 * np.ones(n_pts)))\nbvecs = np.vstack((np.zeros((2, 3)), directions, directions))\n\ngtab = gradient_table(bvals, bvecs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simulations are looped for different intra- and extra-cellular water\nvolume fractions and different intersection angles of the two-fiber\npopulations.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Array containing the intra-cellular volume fractions tested\nf = np.linspace(20, 80.0, num=7)\n\n# Array containing the intersection angle\nang = np.linspace(0, 90.0, num=91)\n\n# Matrix where synthetic signals will be stored\ndwi = np.empty((f.size, ang.size, bvals.size))\n\nfor f_i in range(f.size):\n # estimating volume fractions for individual tensors\n fractions = np.array([100 - f[f_i], f[f_i], 100 - f[f_i], f[f_i]]) * 0.5\n\n for a_i in range(ang.size):\n # defining the directions for individual tensors\n angles = [(ang[a_i], 0.0), (ang[a_i], 0.0), (0.0, 0.0), (0.0, 0.0)]\n\n # producing signals using Dipy's function multi_tensor\n signal, sticks = multi_tensor(gtab, mevals, S0=100, angles=angles,\n fractions=fractions, snr=None)\n dwi[f_i, a_i, :] = signal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that all synthetic signals were produced, we can go forward with\nMSDKI fitting. As other Dipy's reconstruction techniques, the MSDKI model has\nto be first defined for the specific GradientTable object of the synthetic\ndata. For MSDKI, this is done by instantiating the MeanDiffusionKurtosisModel\nobject in the following way:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "msdki_model = msdki.MeanDiffusionKurtosisModel(gtab)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "MSDKI can then be fitted to the synthetic data by calling the ``fit`` function\nof this object:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "msdki_fit = msdki_model.fit(dwi)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the above fit object we can extract the two main parameters of the MSDKI,\ni.e.: 1) the mean signal diffusion (MSD); and 2) the mean signal kurtosis (MSK)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "MSD = msdki_fit.msd\nMSK = msdki_fit.msk" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "kurtosis (MK) from the standard DKI.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dki_model = dki.DiffusionKurtosisModel(gtab)\ndki_fit = dki_model.fit(dwi)\n\nMD = dki_fit.md\nMK = dki_fit.mk(0, 3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we plot the results as a function of the ground truth intersection\nangle and for different volume fractions of intra-cellular water.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig1, axs = plt.subplots(nrows=2, ncols=2, figsize=(10, 10))\n\nfor f_i in range(f.size):\n axs[0, 0].plot(ang, MSD[f_i], linewidth=1.0,\n label='$F: %.2f$' % f[f_i])\n axs[0, 1].plot(ang, MSK[f_i], linewidth=1.0,\n label='$F: %.2f$' % f[f_i])\n axs[1, 0].plot(ang, MD[f_i], linewidth=1.0,\n label='$F: %.2f$' % f[f_i])\n axs[1, 1].plot(ang, MK[f_i], linewidth=1.0,\n label='$F: %.2f$' % f[f_i])\n\n# Adjust properties of the first panel of the figure\naxs[0, 0].set_xlabel('Intersection angle')\naxs[0, 0].set_ylabel('MSD')\naxs[0, 1].set_xlabel('Intersection angle')\naxs[0, 1].set_ylabel('MSK')\naxs[0, 1].legend(loc='center left', bbox_to_anchor=(1, 0.5))\naxs[1, 0].set_xlabel('Intersection angle')\naxs[1, 0].set_ylabel('MD')\naxs[1, 1].set_xlabel('Intersection angle')\naxs[1, 1].set_ylabel('MK')\naxs[1, 1].legend(loc='center left', bbox_to_anchor=(1, 0.5))\n\nplt.show()\nfig1.savefig('MSDKI_simulations.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: MSDKI_simulations.png\n :align: center\n\n MSDKI and DKI measures for data of two crossing synthetic fibers.\n Upper panels show the MSDKI measures: 1) mean signal diffusivity (left\n panel); and 2) mean signal kurtosis (right panel).\n For reference, lower panels show the measures obtained by standard DKI:\n 1) mean diffusivity (left panel); and 2) mean kurtosis (right panel).\n All estimates are plotted as a function of the intersecting angle of the\n two crossing fibers. Different curves correspond to different ground truth\n axonal volume fraction of intra-cellular space.\n\nThe results of the above figure, demonstrate that both MSD and MSK are\nsensitive to axonal volume fraction (i.e. a microstructure property) but are\nindependent to the intersection angle of the two crossing fibers (i.e.\nindependent to properties regarding fiber orientation). In contrast, DKI\nmeasures seem to be independent to both axonal volume fraction and\nintersection angle.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reconstructing diffusion data using MSDKI\n\nNow that the properties of MSDKI were illustrated, let's apply MSDKI to in-vivo\ndiffusion-weighted data. As the example for the standard DKI\n(see `example_reconst_dki`), we use fetch to download a multi-shell\ndataset which was kindly provided by Hansen and Jespersen (more details about\nthe data are provided in their paper [Hansen2016]_). The total size of the\ndownloaded data is 192 MBytes, however you only need to fetch it once.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fraw, fbval, fbvec, t1_fname = get_fnames('cfin_multib')\n\ndata, affine = load_nifti(fraw)\nbvals, bvecs = read_bvals_bvecs(fbval, fbvec)\ngtab = gradient_table(bvals, bvecs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before fitting the data, we perform some data pre-processing. For illustration,\nwe only mask the data to avoid unnecessary calculations on the background of\nthe image; however, you could also apply other pre-processing techniques.\nFor example, some state of the art denoising algorithms are available in DIPY_\n(e.g. the non-local means filter `example-denoise-nlmeans` or the\nlocal pca `example-denoise-localpca`).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "maskdata, mask = median_otsu(data, vol_idx=[0, 1], median_radius=4, numpass=2,\n autocrop=False, dilate=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have loaded and pre-processed the data we can go forward\nwith MSDKI fitting. As for the synthetic data, the MSDKI model has to be first\ndefined for the data's GradientTable object:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "msdki_model = msdki.MeanDiffusionKurtosisModel(gtab)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The data can then be fitted by calling the ``fit`` function of this object:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "msdki_fit = msdki_model.fit(data, mask=mask)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's then extract the two main MSDKI's parameters: 1) mean signal diffusion\n(MSD); and 2) mean signal kurtosis (MSK).\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "MSD = msdki_fit.msd\nMSK = msdki_fit.msk" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For comparison, we calculate also the mean diffusivity (MD) and mean kurtosis\n(MK) from the standard DKI.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dki_model = dki.DiffusionKurtosisModel(gtab)\ndki_fit = dki_model.fit(data, mask=mask)\n\nMD = dki_fit.md\nMK = dki_fit.mk(0, 3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's now visualize the data using matplotlib for a selected axial slice.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "axial_slice = 9\n\nfig2, ax = plt.subplots(2, 2, figsize=(6, 6),\n subplot_kw={'xticks': [], 'yticks': []})\n\nfig2.subplots_adjust(hspace=0.3, wspace=0.05)\n\nim0 = ax.flat[0].imshow(MSD[:, :, axial_slice].T * 1000, cmap='gray',\n vmin=0, vmax=2, origin='lower')\nax.flat[0].set_title('MSD (MSDKI)')\n\nim1 = ax.flat[1].imshow(MSK[:, :, axial_slice].T, cmap='gray',\n vmin=0, vmax=2, origin='lower')\nax.flat[1].set_title('MSK (MSDKI)')\n\nim2 = ax.flat[2].imshow(MD[:, :, axial_slice].T * 1000, cmap='gray',\n vmin=0, vmax=2, origin='lower')\nax.flat[2].set_title('MD (DKI)')\n\nim3 = ax.flat[3].imshow(MK[:, :, axial_slice].T, cmap='gray',\n vmin=0, vmax=2, origin='lower')\nax.flat[3].set_title('MK (DKI)')\n\nfig2.colorbar(im0, ax=ax.flat[0])\nfig2.colorbar(im1, ax=ax.flat[1])\nfig2.colorbar(im2, ax=ax.flat[2])\nfig2.colorbar(im3, ax=ax.flat[3])\n\nplt.show()\nfig2.savefig('MSDKI_invivo.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure::MSDKI_invivo.png\n :align: center\n\n MSDKI measures (upper panels) and DKI standard measures (lower panels).\n\nThis figure shows that the contrast of in-vivo MSD and MSK maps (upper panels)\nare similar to the contrast of MD and MSK maps (lower panels); however, in the\nupper part we insure that direct contributions of fiber dispersion were\nremoved. The upper panels also reveal that MSDKI measures are let sensitive\nto noise artefacts than standard DKI measures (as pointed by [NetoHe2018]_),\nparticularly one can observe that MSK maps always present positive values\nin brain white matter regions, while implausible negative kurtosis values are\npresent in the MK maps in the same regions.\n\n## Relationship between MSDKI and SMT2\nAs showed in [NetoHe2019]_, MSDKI captures the same information than the\nspherical mean technique (SMT) microstructural models [Kaden2016b]_. In this\nway, the SMT model parameters can be directly computed from MSDKI.\nFor instance, the axonal volume fraction (f), the intrisic diffusivity (di),\nand the microscopic anisotropy of the SMT 2-compartmental model [NetoHe2019]_\ncan be extracted using the following lines of code:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "F = msdki_fit.smt2f\nDI = msdki_fit.smt2di\nuFA2 = msdki_fit.smt2uFA" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig3, ax = plt.subplots(1, 3, figsize=(9, 2.5),\n subplot_kw={'xticks': [], 'yticks': []})\n\nfig3.subplots_adjust(hspace=0.4, wspace=0.1)\n\nim0 = ax.flat[0].imshow(F[:, :, axial_slice].T,\n cmap='gray', vmin=0, vmax=1, origin='lower')\nax.flat[0].set_title('SMT2 f (MSDKI)')\n\nim1 = ax.flat[1].imshow(DI[:, :, axial_slice].T * 1000, cmap='gray',\n vmin=0, vmax=2, origin='lower')\nax.flat[1].set_title('SMT2 di (MSDKI)')\n\nim2 = ax.flat[2].imshow(uFA2[:, :, axial_slice].T, cmap='gray',\n vmin=0, vmax=1, origin='lower')\nax.flat[2].set_title('SMT2 uFA (MSDKI)')\n\nfig3.colorbar(im0, ax=ax.flat[0])\nfig3.colorbar(im1, ax=ax.flat[1])\nfig3.colorbar(im2, ax=ax.flat[2])\n\n\nplt.show()\nfig3.savefig('MSDKI_SMT2_invivo.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure::MSDKI_SMT2_invivo.png\n :align: center\n\n SMT2 model quantities extracted from MSDKI. From left to right, the figure\n shows the axonal volume fraction (f), the intrinsic diffusivity (di), and\n the microscopic anisotropy of the SMT 2-compartmental model [NetoHe2019]_.\n\nThe similar contrast of SMT2 f-parameter maps in comparison to MSK (first panel\nof Figure 3 vs second panel of Figure 2) confirms than MSK and F captures the\nsame tissue information but on different scales (but rescaled to values between\n0 and 1). It is important to note that SMT model parameters estimates should\nbe used with care, because the SMT model was shown to be invalid NetoHe2019]_.\nFor instance, although SMT2 parameter f and uFA may be a useful normalization\nof the degree of non-Gaussian diffusion (note than both metrics have a range\nbetween 0 and 1), these cannot be interpreted as a real biophysical estimates\nof axonal water fraction and tissue microscopic anisotropy.\n\n\n### References\n.. [Jensen2005] Jensen JH, Helpern JA, Ramani A, Lu H, Kaczynski K (2005).\n Diffusional Kurtosis Imaging: The Quantification of\n Non_Gaussian Water Diffusion by Means of Magnetic Resonance\n Imaging. Magnetic Resonance in Medicine 53: 1432-1440\n.. [NetoHe2015] Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015).\n Exploring the 3D geometry of the diffusion kurtosis tensor -\n Impact on the development of robust tractography procedures and\n novel biomarkers, NeuroImage 111: 85-99\n.. [NetoHe2018] Henriques RN, 2018. Advanced Methods for Diffusion MRI Data\n Analysis and their Application to the Healthy Ageing Brain\n (Doctoral thesis). Downing College, University of Cambridge.\n https://doi.org/10.17863/CAM.29356\n.. [NetoHe2019] Neto Henriques R, Jespersen SN, Shemesh N (2019). Microscopic\n anisotropy misestimation in spherical\u2010mean single diffusion\n encoding MRI. Magnetic Resonance in Medicine (In press).\n doi: 10.1002/mrm.27606\n.. [Kaden2016b] Kaden E, Kelm ND, Carson RP, Does MD, Alexander DC (2016)\n Multi-compartment microscopic diffusion imaging. NeuroImage\n 139: 346-359.\n.. [Hansen2016] Hansen, B, Jespersen, SN (2016). Data for evaluation of fast\n kurtosis strategies, b-value optimization and exploration of\n diffusion MRI contrast. Scientific Data 3: 160072\n doi:10.1038/sdata.2016.72\n\n.. include:: ../links_names.inc\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 }