This MATLAB toolbox contains all necessary functions for parameter estimation of the Standard Model (SM) of diffusion in white matter1. Check our recent paper for details on this implementation and on the Standard Model in general. Below we provide instructions on how to run the toolbox. See the 'example.m' script that performs the parameter estimation in an example dataset.
For a python implementation of SMI check our TMI package here.
Over the last 15-20 years, multiple approaches aimed to model the physics of water diffusion in white matter have relied on similar assumptions. This led to the unifying framework dubbed Standard Model (SM) of diffusion in WM as formulated in (Novikov et al., 2018)2. In a nutshell, this model disentangles signal contributions from different structures, i.e. compartments, present in a white matter voxel.
Briefly, axons (and possibly glial processes) are represented by impermeable zero-radius cylinders (the so-called “sticks”) arranged in locally coherent fiber fascicles. The diffusion in the extra-axonal space of each fascicle is assumed to be Gaussian and described by an axially symmetric diffusion tensor. The third, optional tissue compartment is the cerebro-spinal fluid. Such multicomponent fascicles (also called kernel) are distributed in a voxel according to an arbitrary fiber orientation distribution function (ODF). All fascicles in a voxel are assumed to have the same compartment fractions and diffusivities, and differ from each other only by orientation.
The SM encompasses a number of WM models made of anisotropic Gaussian compartments with axons represented by sticks (Kroenke et al., 2004; Jespersen et al., 2007, 2010; Fieremans et al., 2011; Zhang et al., 2012; Sotiropoulos et al., 2012; Jensen et al., 2016; Jelescu et al., 2016a; Kaden et al., 2016; Reisert et al., 2017; Novikov et al., 2018; Veraart et al., 2018, to mention a few). From the SM point of view, earlier models impose constraints either on compartment parameters or the functional form of the fiber ODF; such constraints improve robustness but may introduce biases into the estimation of remaining parameters.
For more details please look at our recent publication: Reproducibility of the Standard Model of diffusion in white matter on clinical MRI systems, Volume 257, 15 August 2022, 119290, (2022), NeuroImage, or Section 3 in this review by Novikov et al. (2019).
- Sufficient coarse-graining for Gaussian diffusion at the fiber bundle level (thus, no time-dependence)
- All fiber bundles inside a voxel have the same microstructural properties
- Negligible water exchange between compartments in a bundle
- No assumptions on the fibre bundle's diffusivities
- No assumptions on the functional form of the ODF
Note that this model does not apply to gray matter, where exchange and the presence of soma cannot be ignored.
- This implementation of the SM supports as input a 4D array [nx x ny x nz x N] of diffusion-weighted data, meaning a 3D spatial arrangement of voxels with N diffusion measurements arranged along 4th dimention.
- Protocol files (FSL format) are needed: a b-value (.bval) file, a [1 x N] vector, and a directions file (.bvec), a [3 x N] array.
- A noise map, this must be a 3D array of size [nx x ny x nz]. This can be previously estimated using MPPCA, see this Matlab implementation.
On the b-value units: We recommend microstructural units [milliseconds / (squared micrometers)]. Note that typical scanner units are b=1000 [seconds / (squared millimeters)] which equal 1 [milliseconds / (squared micrometers)], thus to convert from scanner units to microstructural units one should divide by 1000.
Input data can also have:
- Multiple B-tensor shapes (specified by β). This input must be a [1 x N] vector. Only axially symmetric b-tensors are supported. β is a unitless scalar between -0.5 and 1 that indicates the B-tensor shape (see figure below).
- Multiple echo times (TE, assumed to be in [milliseconds]). This input must be a [1 x N] vector.
In this general scenario, each measurement is thus fully specified by: a b-value (b), a unit direction (u) (axis of symmetry of B), a b-tensor shape (β), and TE. See the figure and equation below to understand how these parameters make a b-tensor B:
- If no β is supplied (this input can be absent or an empty array) the code assumes β=1 (linear tensor encoding, LTE, measurements).
- If no TE is supplied (this input can be absent or an empty array), the code assumes that the TE is the same across all measurements. In this case compartmental T2 values will not be outputted and water fractions will be T2-weighted.
Standard Model parameters (diffusion) + compartmental T2 values (only if multiple TE data was provided). See figure below for some examples:
- Fractions (f, fw) and anisotropy (p2,p4) are adimensional.
- Diffusivities are returned in microstructural units [squared micrometers / milliseconds] (independently of input units).
- Compartmental T2 values are returned in [milliseconds].
- Note that compartmental T2 maps will only be outputed if variable TE data was used.
- Rotational invariants for each diffusion shell will also be an output. To normalize them and remove T2 contrast, simply divide by b0 for fixed TE data or by proton density for variable TE data.
We recommend using the default options but the code provides users with some flexibility. Recommended inputs:
-
Diffusion data (4D array) + protocol information + mask (binary 3D array) + noise map (3D array).
- If a mask is not provided the fit will be performed in all voxels present in the first 3 dimensions of the 4D array.
- If a noise map is not provided (not recommended), it will be estimated using the repetitions of the non-diffusion-weighted images.
- We recommend preprocessing your raw data with DESIGNER before doing SM estimation. DESIGNER outputs a robust noise estimation, otherwise MPPCA can be used to robustly estimate the noise level.
-
The current SMI implementation is written in Matlab, future work may translate it to other languages.
Example usage1
We provide three example datasets here that were preprocessed with DESIGNER. These contain
- Multiple b-values
- Multiple b-values and tensor shapes
- Multiple b-values, tensor shapes, and echo times.
See the example file 'example.m' with the SMI fit. Note that data is provided in nifti format, we used this nifti toolbox but feel free to use any and modify the lines where the data is loaded.
% =====================================================================
% Minimal usage:
% Add SMI.m to the path, e.g.:
addpath('/Documents/SantiagoCoelho/Git/SMI')
% Load dwi, protocol, and mask
% Specify protocol information
options.b = bval;
options.dirs = dirs;
% Specify mask and noise map
options.mask = logical(mask);
options.sigma = sigma;
% Run SM fitting (dwi is a 4D array)
[out] = SMI.fit(dwi,options);
% =====================================================================
% More advanced usage:
% Add SMI.m to the path, e.g.:
addpath('/Documents/SantiagoCoelho/Git/SMI')
% Load dwi, protocol, and mask
% Specify protocol information
options.b = bval;
options.beta = beta;
options.dirs = dirs;
options.TE = TE;
% Specify mask and noise map
options.mask = logical(mask);
options.sigma = sigma;
% Specify options for the fit
options.compartments = {'IAS','EAS','FW'}; % The order does not matter
options.NoiseBias = 'Rician'; % use 'None' for zero-mean noise
options.MLTraining.bounds = [0.05, 1, 1, 0.1, 0, 50, 50, 0.05; 0.95, 3, 3, 1.2, 0.5, 150, 120, 0.99];
% The order is: [f, Da, Depar, Deperp, fw, T2a, T2e, p2] (If data has
% fixed TE then the T2a and T2e priors are simply ignored)
% Run SM fitting (dwi is a 4D array)
[out] = SMI.fit(dwi,options);
% =====================================================================
First, the rotational invariants of the diffusion signal are estimated. This is done using a least squares estimator. If data has a non-negligible rician bias, we suggest to correct it during the fitting (see advanced options below).
Then, the SM parameters are estimated from the rotational invariants of the signal (up to L=4 as default option). Unlike conventional parameter estimation approaches which rely on an analytical forward model, e.g. maximum likelihood, here we use data-driven machine learning (ML) regression. This is done by applying a sufficiently flexible nonlinear regression to training data generated with the forward model of interest, considering a wide distribution of model parameters, the noise level, and the protocol that was used. Then, such regression is applied to the data of interest.
For typical SNR values found in clinical dMRI experiments, the optimal regression, i.e. minimizing MSE of the training data, can be achieved already by a cubic polynomial, as it captures all relevant degrees of freedom in the data represented by the set of rotational invariants.
The code provides some additional flexibility:
- Rician bias correction: to de-bias the DWI before the spherical harmonics fit.
- Variable number of compartments: 'IAS', 'EAS', 'FW'. This will be extended but at the moment the only two options are {'IAS', 'EAS'} (default) or {'IAS', 'EAS', 'FW'}.
- User-defined parameter distributions for the training data (for the machine learning estimator that performs RotInvs -> kernel).
- Output spherical harmonic decomposition of the ODF for fiber tracking (normalized for using it with MRtrix3).
The Standard Model is very complex and this is why noise propagates into the model parameters nonlinearly. This results in the kernel diffusivities and additional compartments being very challenging to estimate. If you have multiple TE and b-tensor shapes your chances of getting accurate and precise parameters are much better but if you only have two-shell data then you will likely only get reliable axonal fraction and p2.
- For conventional 2-shell datasets (b=1000,2000 s/mm^2) set 'options.compartments' to {'IAS', 'EAS'} (no FW compartment).
- For ex-vivo data modify the training distribution bounds to account for the decreased diffusivity values.
- This paper in Imaging Neuroscience discusses how to get the most from 2-shell datasets using multiple approaches for microstructure estimation (SMI, WMTI, NODDI, SMT).
- This paper soon to be published in Human Brain Mapping explores how variable TE helps determine the free-water fraction.
Do not hesitate to reach out to Santiago.Coelho@nyulangone.org (or @santicoelho in Twitter) for feedback, suggestions, questions, or comments1.
A US patent contains some of the related developments.
% Authors: Santiago Coelho (santiago.coelho@nyulangone.org), Jelle Veraart, Els Fieremans, Dmitry Novikov
% Copyright (c) 2022 New York University
%
% Permission is hereby granted, free of charge, to any non-commercial entity ('Recipient') obtaining a
% copy of this software and associated documentation files (the 'Software'), to the Software solely for
% non-commercial research, including the rights to use, copy and modify the Software, subject to the
% following conditions:
%
% 1. The above copyright notice and this permission notice shall be included by Recipient in all copies
% or substantial portions of the Software.
%
% 2. THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
% NOT LIMITED TO THE WARRANTIESOF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
% IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BELIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
% WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF ORIN CONNECTION WITH THE
% SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% 3. In no event shall NYU be liable for direct, indirect, special, incidental or consequential damages
% in connection with the Software. Recipient will defend, indemnify and hold NYU harmless from any
% claims or liability resulting from the use of the Software by recipient.
%
% 4. Neither anything contained herein nor the delivery of the Software to recipient shall be deemed to
% grant the Recipient any right or licenses under any patents or patent application owned by NYU.
%
% 5. The Software may only be used for non-commercial research and may not be used for clinical care.
%
% 6. Any publication by Recipient of research involving the Software shall cite the references listed
% below.
%
% REFERENCES:
% - Coelho, S., Baete, S., Lemberskiy, G., Ades-Aron, B., Barrol, G., Veraart, J., Novikov, D.S., Fieremans, E., 2022. Reproducibility of the Standard Model of diffusion in white matter on clinical MRI systems. Neuroimage. doi: 10.1016/j.neuroimage.2022.119290. Epub 2022 May 8. PMID: 35545197; PMCID: PMC9248353.
% - Novikov, D.S., Veraart, J., Jelescu, I.O., Fieremans, E., 2018. Rotationally-invariant mapping of scalar and orientational metrics of neuronal microstructure with diffusion MRI. NeuroImage 174, 518 – 538.
% - Reisert, M., Kellner, E., Dhital, B., Hennig, J., Kiselev, V.G., 2017. Disentangling micro from mesostructure by diffusion MRI: A Bayesian approach. NeuroImage 147, 964 – 975.
Footnotes
-
Please cite these works if you use the SMI toolbox in your publication:
- Coelho, S., Baete, S., Lemberskiy, G., Ades-Aron, B., Barrol, G., Veraart, J., Novikov, D.S., Fieremans, E., 2022. Reproducibility of the Standard Model of diffusion in white matter on clinical MRI systems. Neuroimage. doi: 10.1016/j.neuroimage.2022.119290. Epub 2022 May 8. PMID: 35545197; PMCID: PMC9248353.
- Dmitry S. Novikov, Jelle Veraart, Ileana O. Jelescu, Els Fieremans, Rotationally-invariant mapping of scalar and orientational metrics of neuronal microstructure with diffusion MRI, NeuroImage, Volume 174, 2018, Pages 518-538.
- Marco Reisert, Elias Kellner, Bibek Dhital, Jürgen Hennig, Valerij G. Kiselev, Disentangling micro from mesostructure by diffusion MRI: A Bayesian approach, NeuroImage, Volume 147, 2017, Pages 964-975.
-
The SM name originates from a tongue-in-cheek association with the Standard Model in particle physics, as both encompass a fair bit of previous modeling effort from various groups. Here “standard” refers to common assumptions among modeling approaches and does not imply “exact”. ↩