PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of blackwellopenThis ArticleFor AuthorsLearn MoreSubmit
Nmr in Biomedicine
 
NMR Biomed. 2017 September; 30(9): e3734.
Published online 2017 June 23. doi:  10.1002/nbm.3734
PMCID: PMC5563694

Diffusion MRI microstructure models with in vivo human brain Connectome data: results from a multi‐group comparison

Abstract

A large number of mathematical models have been proposed to describe the measured signal in diffusion‐weighted (DW) magnetic resonance imaging (MRI). However, model comparison to date focuses only on specific subclasses, e.g. compartment models or signal models, and little or no information is available in the literature on how performance varies among the different types of models. To address this deficiency, we organized the ‘White Matter Modeling Challenge’ during the International Symposium on Biomedical Imaging (ISBI) 2015 conference. This competition aimed to compare a range of different kinds of models in their ability to explain a large range of measurable in vivo DW human brain data. Specifically, we assessed the ability of models to predict the DW signal accurately for new diffusion gradients and b values. We did not evaluate the accuracy of estimated model parameters, as a ground truth is hard to obtain. We used the Connectome scanner at the Massachusetts General Hospital, using gradient strengths of up to 300 mT/m and a broad set of diffusion times. We focused on assessing the DW signal prediction in two regions: the genu in the corpus callosum, where the fibres are relatively straight and parallel, and the fornix, where the configuration of fibres is more complex. The challenge participants had access to three‐quarters of the dataset and their models were ranked on their ability to predict the remaining unseen quarter of the data. The challenge provided a unique opportunity for a quantitative comparison of diverse methods from multiple groups worldwide. The comparison of the challenge entries reveals interesting trends that could potentially influence the next generation of diffusion‐based quantitative MRI techniques. The first is that signal models do not necessarily outperform tissue models; in fact, of those tested, tissue models rank highest on average. The second is that assuming a non‐Gaussian (rather than purely Gaussian) noise model provides little improvement in prediction of unseen data, although it is possible that this may still have a beneficial effect on estimated parameter values. The third is that preprocessing the training data, here by omitting signal outliers, and using signal‐predicting strategies, such as bootstrapping or cross‐validation, could benefit the model fitting. The analysis in this study provides a benchmark for other models and the data remain available to build up a more complete comparison in the future.

Keywords: brain microstructure, Connectome, diffusion MRI, fornix, genu, model selection

Abbreviations used

CT
computerized tomography
CV
cross‐validation
DBF
diffusion basis function
DF
diffusion function
DT
diffusion tensor
DTI
diffusion tensor imaging
DW
diffusion‐weighted
EN
elastic net
ISBI
International Symposium on Biomedical Imaging
LASADD
Linear Acceleration of Sparse and Adaptive Diffusion Dictionary
LS
least‐squares
MRI
magnetic resonance imaging
PDD
principal diffusion directions
RSI
restriction spectrum imaging
ROI
region of interest
SFM
sparse fascicle model
SNR
signal‐to‐noise ratio
SSE
sum of squared errors
TE
echo time
WM
white matter

1. INTRODUCTION

Diffusion‐weighted (DW) magnetic resonance imaging (MRI) can provide unique insights into the microstructure of living tissue and is increasingly used to study the microanatomy and development of normal functioning tissue as well as its pathology. Mathematical models for analysis and interpretation have been crucial for the development and translation of DW‐MRI. Even though diffusion tensor imaging (DTI),1 which is based on a simple Gaussian model of the DW‐MRI signal, has shown promise in clinical applications,2 e.g. Alzheimer's disease,3 multiple sclerosis4 or brain tumors,5 a much wider variety of DW‐MRI models has been proposed to extract more information from the DW signal.

Models generally fall between two extremes: ‘models of the tissue’ and ‘models of the signal’. Models of the tissue6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17 describe the underlying tissue microstructure in each voxel explicitly with a multi‐compartment approach.18, 19, 20 Models of the signal focus on describing the DW signal attenuation without describing the underlying tissue composition that gives rise to the signal explicitly.21, 22, 23, 24, 25, 26, 27, 28, 29 Other approaches fall between these two classes and include some features of the tissue, such as the distribution of fibre orientations, but often describe the signal from individual fibres without modelling the fibre composition explicitly.30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40

Despite this explosion of DW‐MRI models, a broad comparison on a common dataset and within a common evaluation framework is lacking, so little is understood about which models are more plausible representations or explanations of the signal. Panagiotaki et al.18 established a taxonomy of diffusion compartment models and compared 47 of them using data from the fixed corpus callosum of a rat acquired on a pre‐clinical system. Later, Ferizi et al.39 performed a similar experiment using data from a live human subject, while Ferizi et al.41, 42 explored a different class of models, which aim to capture fibre dispersion. Rokem et al.43 compared two classes of models using cross‐validation and test–retest accuracy. All these studies18, 43, 44 aim to evaluate variations with specific classes of models with all other variables of the parameter estimation pipeline (i.e. noise model, fitting routine, etc.) fixed. While this provides fundamental insight into which compartments are important in compartment models, questions remain about the broader landscape of models; in particular, which classes of models explain the signal best and how strongly performance depends on the choice of parameter‐estimation procedure.

Publicly organized challenges provide a unique opportunity to bring a research community together to gain a quantitative and unbiased comparison of a diverse set of methods applicable to a particular data‐processing task. Such publicly organized challenges have helped to establish a common ground for the evaluation of competing methods in a variety of imaging‐related tasks, e.g. in brain MR image registration45 and segmentation.46 In DW‐MRI, public challenges have focused on recovering synthetic intra‐voxel fibre configurations47 or evaluating tractography techniques48, 50 and have been very successful at driving research and translation forward. Another interesting comparison of reconstruction methods using DW‐MRI data was based on the signal acquired from a physical phantom.49 Here we report on such a community‐wide challenge to model the variation of DW‐MRI signals at the voxel level in the in vivo human brain.

Modelling the diffusion signal is a key step in realizing practical and reliable quantitative imaging techniques based on diffusion MRI. The challenge in the area is to extract the salient features from the diffusion signal and relate them to the principal features of the underlying tissue (e.g. in the case of brain white matter (WM) the fibre orientation, axonal packing and axonal size). Three distinct questions arise.

  1. Given the richest possible dataset that samples the space of achievable measurements as widely as possible, which mathematical model can capture best the intrinsic variation of the acquired signal?
  2. Which tissue features can be derived from the model?
  3. What subset of those features can we estimate from limited acquisition time on a standard clinical scanner and what dataset best supports such estimates?

The intuition gained from (i) is generalizable over a wide range of applications, while (ii) and (iii) are highly dependent on the MRI study design and the available hardware. Therefore, our challenge focuses on question (i), as an understanding of (i) is necessary to inform (ii) and (iii). To that end, we acquire the richest possible dataset using the most powerful hardware available and the most motivated subject available (UF). Specifically, we use the Connectome scanner,51 which is unique among human scanners in having 300 mT/m gradients, rather than 40 mT/m as is typical of state‐of‐the‐art human scanners. Preclinical work by Dyrby et al.13 highlights the benefits of such strong gradients and the first results from the Connectome scanner42, 52, 53, 54 are now starting to verify those findings.

This kind of model comparison, based on prediction error, is a common and crucial part of the development of any statistical model‐based estimation applications. Burnham and Anderson55 explain how and why such comparisons should be performed to reject models that are theoretically plausible but not supported by the data. To that end, we used a uniquely rich dataset acquired on the Connectome system42 composed of around 5000 points in q space with, for each shell, a unique combination of gradient strength, diffusion time, pulse width and echo time. This offers the opportunity for the comparison of the many different types of models within a common framework, over a very wide range of measurement space. Using this rich dataset, we organized the White Matter Modeling challenge, held during the 2015 International Symposium on Biomedical Imaging (ISBI) in New York. The goal of the challenge was to evaluate and compare the models in two different tissue configurations that are common in the brain: (1) a WM region of interest where fibres are relatively straight and parallel, specifically the genu of the corpus callosum; and (2) a region in which the fibre configuration is more complex, specifically the fornix. Challenge participants had access to three‐quarters of each whole dataset; the participating models were evaluated on how well they predicted the remaining ‘unseen’ part of the data. As announced before the challenge, the final ranking was based exclusively on the performance on the genu data. In this article, however, we include results from both the genu and the fornix.

The article is organized as follows. We first describe in section 2 the experimental protocol, data post‐processing and preparation of the training and testing data for the challenge. We then present the methods for ranking the models and tabulate the various models involved in the competition succinctly. We report the challenge results in section 3 and discuss these results in section 4; a more detailed description of the models follows in the Appendix.

2. MATERIAL AND METHODS

2.1. The complete experiment protocol

One healthy volunteer was scanned over two non‐stop 4 h sessions. The imaged volume comprised twenty 4 mm thick whole‐brain sagittal slices covering the corpus callosum left–right. The image size was 110×110 and the in‐plane resolution 2×2 mm2. 45 unique and evenly distributed diffusion directions (taken from http://www.camino.org.uk) were acquired for each shell, with both positive and negative polarities; these directions were the same in each shell. We also included 10 interleaved b=0 measurements, leading to a total of 100 measurements per shell. Each shell had a unique combination of Δ={22,40,60,80,100,120} ms, δ={3,8} ms and |G|={60,100,200,300} mT/m (see Table 1). The measurements were randomized within each shell, whereas the gradient strengths were interleaved. We inspected the images visually and did not observe any obvious shifts from gradient heating. The minimum possible echo time (TE) for each gradient duration and diffusion time combination was chosen to enhance signal‐to‐noise ratio (SNR) and potential estimation of compartment‐specific relaxation constants. The SNR of b=0 images was 35 at TE = 49 ms and 6 at TE = 152 ms. The SNR was computed by assessing the signal mean and noise variance across the selected WM voxels on multiple b=0 images. In both cases these estimates matched reasonably well. More details about the acquisition protocol can be found in Ferizi et al.42

Table 1

The scanning protocol used, acquired in ~8 hours over two non‐stop sessions. The protocol has 48 shells, each with 45 unique gradient directions (‘blip‐up‐blip‐down’)

2.2. Post‐processing

All post‐processing was performed using Software Library (FSL).56 The DW images were corrected for eddy current distortions separately for each combination of δ and Δ using FSL's E d d y module (www.fmrib.ox.ac.uk/fsl/eddy) with its default settings. The images were then co‐registered using FSL's F n i r t package. As the 48 shells were acquired across a wide range of TEs, over two days, we chose to proceed in two steps. First, within each quarter of the dataset (different day, different δ) we registered all the b=0 images together. We then applied these transformations to their intermediary DW images, using a trilinear resampling interpolation. The second stage involved co‐registering the four different quarters. To help the co‐registration, especially between the two days images that required some through‐plane adjustment as well, we omitted areas of considerable eddy‐current distortions by reducing the number of slices from 20 to 5 (i.e. leaving two images either side of the mid‐sagittal plane) and reducing the in‐plane image size to 75×80.

2.3. Training and testing data

The data for this work originated from two regions of interest (ROIs), each containing 6 voxels (see Figure Figure1).1). The first ROI was selected in the middle of the genu in the corpus callosum, where the fibres are mostly straight and coherent. The second ROI's fibre configuration is more complex: it lies in the body of fornix, where two bundles of fibres bend and bifurcate.

nbm3734-fig-0001

We only consider two ROIs, each containing six voxels from the genu in the corpus callosum, where the fibres are approximately straight and parallel, and from the fornix, where the configuration of fibres is more complex

The dataset was split into two parts: the training dataset and the testing dataset. The training dataset was fully available for the challenge participants. The testing dataset was retained by one of the organizers (UF). The DW signal of the training dataset (36 shells, with acquisition parameters shown in black in Table 1) was provided together with the gradient scheme on the challenge website (http://cmic.cs.ucl.ac.uk/wmmchallenge/). This data was used by the participants to estimate their DW‐MRI model parameters. The signal attenuation in the testing dataset (12 shells, with acquisition parameters shown in red in Table 1) was kept unseen. It contained one shell, chosen at random, from each TE‐specific set of four shells (i.e of the same combination of δ and Δ). The challenge participants were then asked to predict the signal for the corresponding gradient scheme. They were free to use as much or as little of the training data provided as they wished to predict the signal of the test dataset for the six voxels in each ROI.

Figure Figure22 shows the DW signal attenuation for each shell in the genu dataset, with stars in the legend indicating which shells were left out for testing. In this plot, a small number of data appear as ‘outliers’ (two such data are shown with arrows in the bottom‐left subplot of Figure Figure2).2). Specifically, we counted about 10 of them among more than 4812 measurements, most of them being in the b=300 s/mm2 shell. Since these outliers appear to be specific to the b=300 s/mm2 shell and are not in other shells with similar b value, we attribute them to a momentary twitching of the subject rather than more systematic effects, such as perfusion.

nbm3734-fig-0002

Diffusion‐weighted signal from the genu ROI, averaged over the six voxels. Across each column and row, the signal pertains to one of the gradient strengths or pulse times δ used; in each subplot, the six shells shown in different colours ...

Similarly, Figure Figure33 shows the signal for the fornix region, with the signal over the six voxels averaged out.

nbm3734-fig-0003

Diffusion‐weighted signal from the fornix ROI, averaged over the six voxels. The legend's b value is in s/mm2 units. Testing shells are marked with a star. On the x‐axis is the cosine of the angle between the applied diffusion gradient ...

2.4. Model ranking

Models were evaluated and ranked based on their ability to predict the unseen DW signal accurately. Specifically, the metric used was the sum of square differences between the hidden signal and the predicted signal, corrected for Rician noise:57

SSE=1Ni=1N(S˜iSi2+σ2)2σ2
(1)

where N is the number of measurements, S˜i is the ith measured signal, S i its prediction from the model and σ the noise standard deviation.

2.5. Competing models

Here we give a short summary of the competing models. Additionally, Table 2 provides a summary of their key characteristics. More details are included in the Appendix.

  • Ramirez‐Manzanares: a dictionary‐based technique that accounts for multiple fibre bundles and models the distribution of tissue properties (axon radius, parallel diffusivity) and the orientation dispersion of fibres.
  • Nilsson: a multi‐compartment model that models isotropic, hindered and restricted diffusion and accounts for varying (T 1, T 2) relaxation times for each compartment.58
  • Scherrer a multi‐compartment model in which each compartment is modelled by a statistical distribution of 3‐D tensors.16
  • Ferizi1 and Ferizi2: two three‐compartment models that account for varying T 2 relaxation times for each compartment. As regards the intracellular compartment, Ferizi1 models the orientation dispersion by using dispersed sticks as one compartment; Ferizi2 uses a single radius cylinder instead.42
  • Poot: a three‐compartment model comprising an isotropic diffusion compartment, a tensor compartment and a model‐free compartment in which an Apparent Diffusion Coefficient (ADC) is estimated for each direction independently. T 2 relaxation times are also estimated for each compartment.59
  • Rokem: a combination of the sparse fascicle model43 with restriction spectrum imaging60 that describes the signal arising from a multi‐compartment model in a densely sampled spherical grid, using L1 regularization to enforce sparsity.
  • Eufracio: an extension of the Diffusion Basis Function (DBF) model that accounts for multiple b‐value shells.
  • Loya‐Olivas1 and Loya‐Olivas2: two models based on the Linear Acceleration of Sparse and Adaptive Diffusion Dictionary (LASADD) technique. Loya‐Olivas1 uses the DBF signal model, while Loya‐Olivas2 uses a three‐compartment tissue model. The optimization uses linearized signal models to speed up computation and sparseness constraints to regularize.
  • Alipoor: a model of fourth‐order tensors, corrected for T 2‐relaxation across different shells. A robust LS fitting was applied to mitigate influence of outliers.
  • Sakaie: a two‐compartment model of restricted and hindered diffusion with angular variation. A simple exclusion scheme based on the b=0 signal intensity was applied to remove outliers.
  • Fick: a spatio‐temporal signal model to represent 3‐D diffusion signal simultaneously over varying diffusion time. Laplacian regularization was applied during the fitting.61
  • Rivera: a regularized linear regression model of diffusion encoding variables. This is intentionally built as a simplistic model to be used as a baseline for model comparison.

Table 2

Summary of the various diffusion models evaluated. Tissue models are models that include an explicit description of the underlying tissue microstructure with a multi‐compartment approach. In contrast, signal models focus on describing the DW signal ...

While the challenge organizers also had competing models (Ferizi1, Ferizi2 and Scherrer), only Ferizi had access to the hidden data. The hidden data were never used to tune the results of his models.

3. RESULTS

Figure Figure44 shows the averaged prediction error in each ROI (top subplot is for the genu, bottom subplot is for the fornix) and the corresponding overall ranking of the participating models in the challenge. The first six models in the genu ranking performed similarly, each higher ranked model marginally improving on the prediction error. The prediction error clearly increased at a higher rate for the subsequent models. In the fornix dataset, the prediction error was higher than in the genu. For both datasets, the first six models were the same, albeit permuted. Most of the models performed similarly in terms of ranking in both genu and fornix cases, i.e. Nilsson (second in genu/first in fornix), Scherrer (third/second) and Ferizi_2 (fourth/fourth). Others performed significantly better in one of the cases, with Ramirez‐Manzanares (first/sixth) being the most notable.

nbm3734-fig-0004

Overall ranking of models by sum‐of‐squared‐errors (SSE) metric over all voxels in genu (top) and fornix (bottom) ROIs. The colors represent different ranges of b‐value shells

Figure Figure44 also details the prediction error for different ranges of b values in the unseen dataset. Models inevitably vary in their prediction capabilities; some models perform better within a given b‐value range but are penalized more in another. Across the models, as the figure shows, the ranking between models was dominated by the signal prediction accuracy for b values between 750 and 1400s/mm2; specifically, the shell that has the largest weight on this error is the b=1100 s/mm2 one. The top‐ranking models, nevertheless, were better at predicting the signal for higher b‐value images as well. The prediction performance of lower b‐value images (<750s/mm2) in the genu was less consistent across ranks. For example, the models of Rokem and Sakaie outperformed most of the higher ranking models in this low b‐value range. The fornix is a more complex region than the genu, hence the performance across the shells is less consistent. In the fornix, the prediction errors were generally larger than in the genu across all b values for all models, except Rivera's, which showed the opposite effect. The prediction errors of the b=0 images were also larger than in the genu, especially for the highly ranked models of Poot and Ferizi. The prediction errors in other b‐value shells followed the overall ranking of the models more closely.

Figure Figure55 shows the prediction error for each voxel independently. In the genu plot, the best performing models had high consistency of low prediction errors across all individual voxels. These were followed by the models with consistent larger prediction error in all voxels. Most of the lowest ranking models not only had largest prediction errors, they also showed large variations in prediction performance. For example, while the model of Loya‐Olivas2 was competitive in voxel 5, it ranked low due to large prediction errors in voxels 4 and 6. The results in the fornix show a lower consistency of prediction errors between the voxels than in the genu. Specifically, two voxels (3 and 4) showed substantially larger prediction errors and were likely responsible for much of the overall ranking.

nbm3734-fig-0005

Sum‐of‐squared‐errors (SSE) per voxel for each model in genu and fornix. The size of rectangles represent the SSE value per voxel

Finally, we report in Figures Figures6,6, ,7,7, ,88 and and99 an illustration of the quality of fit of each model to four representative shells, including the b=1100s/mm2 shell mentioned above; Figures Figures66 and and77 concern the genu data and Figures Figures88 and and99 are for fornix data.

nbm3734-fig-0006

Genu signal for the group consisting of the best seven from 14 models. We show only four (of twelve) representative shells; these are shown by blue stars, while red circles denote the model‐predicted data. The best models are listed first. The ...

nbm3734-fig-0007

Genu signal for the second group of 14 models. Raw testing data are shown by blue stars, while red circles denote the model‐predicted data. The x‐axis is the cosine of the angle between G and n

nbm3734-fig-0008

Fornix signal for the group consisting of the best 7 from 14 models. We show only four (of twelve) representative shells; these are shown by blue stars, while red circles denote model‐predicted data. The best models are listed first. The x‐axis ...

nbm3734-fig-0009

Fornix signal for the second group of 14 models. Raw testing data are shown by blue stars, while red circles denote the model‐predicted data. The x‐axis is the cosine of the angle between G and n

4. DISCUSSION

The challenge set out to compare the ability of various kinds of models to predict the diffusion MR signal from WM over a very wide range of measurement parameters – exploring the boundaries of possible future quantitative diffusion MR techniques. The 14 challenge entries were a good representation of the many available models that are proposed in the literature. The acquired data aimed to cover the broadest spectrum of experimental parameters possible. The participating models use a variety of fitting routines and modelling assumptions, providing additional insight into the effects of algorithmic and modelling choices during parameter estimation. Although the set of methods tested is not sufficient to make a full comparison of each independent feature (diffusion model, noise model, fitting routine, etc.) and the number of combinations prohibits an exhaustive comparison, the results of the challenge do reveal some important trends.

In contrast with earlier model comparisons,18, 43, 44 the results provide new insight into which broad classes of model explain the signal best and what features of the estimation procedure are important. This information is very timely, as recent model‐based diffusion MRI techniques, such as NODDI,15 SMT,17, 40 DIAMOND,16 DKI62 and LEMONADE,63 are starting to become widely adopted in clinical studies and trials. Despite their success, intense debate continues in the field about applicability of different models and fitting routines.64, 65 The insights from this challenge provide key pointers to the important features of the next‐generation of front‐line imaging techniques of this type. Moreover, the data and evaluation routines remain available to form the basis of an expanding ranking of models and fitting routines and a benchmark for future model development.

4.1. Main conclusions

The first insight is on the type of model used. Signal models do not necessarily outrank tissue models; indeed, using our dataset, models of the signal (Alipoor, Sakaie, Fick, Rivera) ranked on average lower than models of the tissues, despite their theoretical ability to offer more flexibility in describing the raw signal. This is quite surprising, as the current perception within the field is that, generally, we can capture the signal variation much better through a functional description of the signal (signal models) rather than via a biophysical model of the tissue (tissue models). The former generally consist of bases of arbitrary complexity, whereas the latter are generally very parsimonious models that rely on extremely crude descriptions of tissue (e.g. white matter as parallel impermeable cylinders). The results suggest that the flexibility of signal models can rapidly lead to overfitting. However, the tissue models can explain the signal relatively well even with just a few parameters (compare the quality‐of‐fit plots of the Rivera model in Figure Figure77 with the signal prediction of the top models in Figure Figure6:6: the higher the b value, the worse the prediction of the linear signal model). Certain underlying assumptions may cause the signal models to perform less well than expected. For example, they are often designed to work with data with a single diffusion time and do not generalize naturally to incorporate the additional dimension (although see Fick et al.61 for some steps towards generalization). Many of the tissue models, on the other hand, naturally account for finite δ, varying diffusion times and gradient strength (e.g. the Ramirez‐Manzanares, Nilsson and Ferizi models in our collection). We cannot draw any conclusion about the benefits of an adjustable number of parameters in a model, because of the limited number of models in our study that do this and because the models differ in a range of other aspects.

The second insight concerns the choice of noise modelling. Despite the fact that SNR at b=0 and T E=152 ms falls to about 6, use of the Rician noise model does not appear to be a significant benefit in predicting unseen signal; here, however, we do not investigate the effect on estimated model parameters, which may still benefit from the more accurate noise model. In this challenge, most participants used non‐linear least‐squares or maximum‐likelihood optimization. Additional regularization of the objective function (Eufracio & Rivera/Lasso, Rokem/Elastic Net, Fick/Laplacian) appeared to have little benefit over non‐regularized optimization.

The third observation is about removing signal outliers. Five of the eleven models preprocessed the training data by clearing out outliers, including the top two models. We tried this procedure with two good models that did not use such a procedure, Ferizi1 and Ferizi2, and observed that it did not affect the ranking, though it did improve the prediction error marginally. This is understandable, considering the relatively little weight these apparent outliers have on the total number of measurements (10 points from a 4812‐strong dataset). Additionally, specific strategies for predicting the signal, e.g. bootstrapping or cross‐validation, as used by the top two models of Ramirez‐Manzanares and Nilsson, may also help the model ranking.

4.2. Limitations and future directions

Although this challenge provides several new insights into the choice of model and fitting procedure for diffusion‐based quantitative imaging tools, it has a number of limitations that future challenges might be designed to address. One limitation of the study is that we use a very rich acquisition protocol that is not representative of common or clinical acquisition protocols. In particular, we cover a very wide range of b values and the data acquisition (protocol) we use consists of many TEs, unlike many other multi‐shell diffusion datasets that use a fixed TE. As stated in the Introduction, our intention is to sample the measurement space as widely as possible to support the most informative models possible. Varying the TE makes it possible to probe compartment‐specific T 2 (the decay of which Ferizi et al.42 finds to be monoexponential at the voxel level), an investigation that would be impossible with a single TE. However, the good performance of DIAMOND also shows that a model with fixed δ and Δ can still capture the signal variation in multi‐TE datasets and that, while the majority of the full data was ignored in each of the reconstructions, its prediction error compared favourably with other techniques.

We use the unique human Connectome scanner51 to acquire a dataset with gradients of up to 300mT/m, which is not readily available in most current MR machines. However, previous preclinical work by Dyrby et al.13 suggests that high diffusion gradients enrich the signal, which helps model fitting and comparison. Future challenges might be designed that focus on explaining the signal and estimating parameters from data more typical of clinical acquisitions.

Assessing the prediction performance on unseen data as in this challenge is different from assessing the fitting error: it implicitly penalizes models that overfit the data. However, since most of the missing shells lie in between other shells (in terms of b values, TEs, etc.), the quality of signal extrapolation was not assessed. We get a glimpse of this from Figure Figure4,4, where the SSE is unevenly distributed between the b values. Here, the shell that bore the largest error is the b=1100 s/mm2 one; see also Figures Figures66 and and7.7. Of all ‘unseen’ shells, this shell combines the lowest Δ and highest |G|, placing it on the edge of the range of the measurement space sampled. Such a b‐value shell combines high signal magnitude with high sensitivity, i.e. the gradient of signal against b‐value is highest in this range, which makes it hard to predict. (We stress that this observation is in the context of the wider multi‐shell acquisition, and is not to be seen in isolation for its potential impact on single‐shell acquisition methods.) On the other hand, the variability of prediction errors in the b<750 s/mm2 range could arise from the varying sensitivity of different models to the free water component, which is challenging to estimate as it can easily be confounded with hindered water, or physiological effects, which are mostly observable in this low b‐value range. Future work can take this further, by selecting unseen shells outside the min–max range of experimental parameters. This is likely to penalize more complex models that overfit the data even more strongly.

We did not take into account the computational demand of each model, and this might limit the generalization of the results. Models that use bootstrapping generally have a higher computational burden and may not be feasible for large datasets, e.g. whole brain coverage.

The dataset used in this challenge is specific to one subject who underwent a long‐duration acquisition, which adds to the question of generalizability. The subsequent preprocessing of the data is also a factor to bear in mind: the registration of two 4h datasets, across such a broad range of echo times, poses its own challenges for certain non‐homogenous regions in the brain, such as the fornix (compared with, for example, the relatively large genu). Thus the results may be somewhat subject‐specific and may be affected by residual alignment errors.

Another limitation is that we only look at isolated voxels inside the corpus callosum and the fornix. Questions still remain about which models are viable even in the most coherent areas of the brain with the simplest geometry, so we believe our focused challenge on well‐defined areas is an informative first step necessary before extending the idea to the whole of the white matter, which would make for an extremely complex challenge. We note, however, recent work by Ghosh et al.66 that illustrates such an approach with Human Connectome Project (HCP) data.

We focused here on comparing models based on their ability to predict unseen data. Although models that reflect true underlying tissue structure should explain the data well, we cannot infer in general that models that predict unseen data better are mechanistically closer to the tissue than those that do not. As we discuss in the Introduction, the main power of evaluating models in terms of prediction error is to reject models that cannot explain the data. Thus, while the identification of parsimonious models that explain the data certainly has great benefit, further validation is necessary through comparison of the parameters that they estimate with independent measurements, e.g. obtained through microscopy (our challenge makes no attempt to assess the integrity of parameter estimates themselves, but future challenges might use such performance criteria). Models can be evaluated to some extent by sanity checking the realism of their fitted parameter values, as in for example Jelescu et al.64 or Burcaw et al.67 However, obtaining accurate ground‐truth values for quantitative evaluation remains a hard and yet unsolved problem for diffusion MRI in general. In particular, histology can only roughly approximate the in vivo ground truth and introduces its own set of challenges in sample preparation, acquisition and biophysical interpretation.12, 13, 65, 68, 69, 70, 71 This challenge highlights the need for improved model comparison and validation methods.

5. CONCLUSION

Challenges such as this have great value in bringing the community together and provide an unbiased comparison of wide‐ranging solutions to key data‐processing problems. They raise new insights and ideas, motivating more directed future studies. The data are publicly available for others to use, with more details of the dataset given on the Challenge website at http://cmic.cs.ucl.ac.uk/wmmchallenge/. On this website, an up‐to‐date ranking of the models will be available, where additional models can be added after the publication of the article and where the community will be able to evaluate further the impact of noise correction, compartment‐specific T 2 estimation, inter‐class model assumptions, e.g. tissue versus signal models, or indeed intra‐class model assumptions, e.g. whether cylinders or sticks are optimal models for the given dataset.42 This will provide an important benchmark for future models and parameter estimation routines.

ACKNOWLEDGEMENTS

Research reported in this manuscript was supported by the following organizations. EPSRC supported this work through grants EP/G007748, EP/L022680/1, EP/I027084/01, EP/M020533/1 and EP/N018702/1. U. Ferizi is also supported by the National Institute of Arthritis and Musculoskeletal and Skin Diseases (NIAMS) of the National Institute of Health (NIH) under award numbers R21AR066897 and R01AR067789. B. Scherrer was supported in part by NIHR01 NS079788, R01 EB019483, R01 EB018988 and BCH TRP Pilot and BCH CTREC K‐to‐R Merit Award. M. Nilsson is supported by the Swedish Strategic Research (SSF) Grant AM13‐ 0090. R. Deriche and R. Fick were supported in part by the European Research Council (ERC) under the Horizon 2020 research and innovation program (ERC Advanced Grant agreement No 694665 : CoBCoM).

The content is solely the responsibility of the authors and does not necessarily represent the official views of the funding bodies (EPSRC or NIH).

Lastly, we thank the ISBI 2015 challenge organizers, in particular Stephen R. Aylward (Kitware Inc., USA) and Badri Roysam (University of Houston, USA).

COMPETING MODELS

A.1.  Tissue models

A.1.1. Ramirez‐Manzanares (CIMAT, Mexico): Empirical Diffusion‐and‐Direction Distributions (ED3)

This work builds on the statistical modelling of the apparent diffusion coefficient72 and tackles the modelling of axon fibre dispersion in single15, 73 and multiple fibre bundle cases. The method estimates empirically (rather than imposes) the distribution of tissue properties (axon radius, parallel diffusion, etc.), as well as the orientational distribution of the bundles. The general framework is as follows:

  • estimation of mean principal diffusion directions (PDD) per axon bundle;
  • selection of a dense set of orientationally focused basis directions that capture the discrete non‐parametric fibre dispersion;
  • design of a dictionary of intra/extracellular synthetic DW signals, which are precomputed along the basis directions (see the DBF method in Ramirez‐Manzanares et al.74);
  • computation of the size compartments per diffusion atom of the dictionary (model fitting).

The PDDs are estimated from the Diffusion Tensor (DT) (single bundle case) and DBF74 (complex structure cases). The 120 orientations closest to the PDDs are selected from a set of 1000 evenly distributed orientations. The intra‐axonal signals S i are precomputed from the model in Van Gelderen et al.75 for restricted diffusion within a cylinder with radius R=1,2,,10μm and parallel diffusion d=1,1.1,,2.1μm2/ms. The extra‐axonal signals are generated as follows: S e from zeppelins with combinations of parallel and radial diffusion, d=1,1.1,,2.5μm2/ms and d=2,3,,8μm2/ms; isotropic diffusion compartment signals Siiso=expqτdiso for diso=2,2.1,2.2,,4μm2/ms; and the dot signal, which takes into account static proton density. The values of the dictionary atoms above were tuned by cross‐validation.76 The size compartments β[gt-or-equal, slanted]0 computed in the weighted non–negative LS formulation,

W[SS0TE(i=1NiβiiSiij=1NeβjeSjek=1NisoβkisoSkiso+βdot)]22
(A1)

indicate the atoms that explain the signal; the W weights are proportional to SNR. Overfitting is reduced by a bootstrap 77 procedure.

The cross‐validation experiments indicate that the reconstructions given by the robust fitting of this rich multi‐compartment diffusion dictionary allow us to predict accurately non‐acquired MR signals for different machine protocols. This is of most interest in the development of methods able to detect the complex microstructure heterogeneity associated with the different compartments within the voxels. The atoms with coefficients β>0 depict the empirical distributions and their orientations indicate non‐parametrical bundle–dispersion configurations (such as fanning or radially symmetrical). The recovered distributions reveal, for instance, an axon radius of between 1 and 4μm. One should take into account, however, that, since the heterogeneous intra/extra‐axonal T 2 relaxation feature is not modelled explicitly, the method may compensate for T 2 variations by using, for instance, large isotropic d iso coefficients to fit the signal accurately. For this reason, a direct interpretation of the fitted parameters may be misleading. The use of more specific models is a part of ongoing work.

A.1.2. Nilsson (Lund, Sweden): multi‐compartment model outlier rejection and separate fitting of b 0data

This multiple compartment model was developed specifically for the ISBI WM challenge and built up by relaxation‐weighted and time‐dependent diffusion tensors according to

S=S0iwieB:DieTE/T2i(1e(TRTE/2)/T1i)
(A2)

where B=bn2 and b=(γ δ g)2 t d. The diffusion time t d was corrected for rise times (ξ) according to t d=Δ−δ/3+ξ 3/30δ 2ξ 2/6δ. Each component was also described by a weight (w i) and relaxation times (T1i and T2i). The model featured three types of component, with either isotropic, hindered or restricted diffusion. Diffusion in the isotropic component was modelled by a single diffusion coefficient. The hindered and restricted components were modelled by cylinder‐symmetric tensors described by axial and radial diffusivities together with the polar and azimuth angles. In the restricted component, the apparent diffusion coefficient of the radial component depended on δ and Δ, as well as on the cylinder radius, according to Van Gelderen et al.78

Three modifications were performed to this very general model. First, to accommodate for potential bias in the b 0 images (which was the case for fornix data, where deviations of up to 20σ was observed), the prediction for b 0 data was obtained from the median of all signals acquired with identical TE instead of from Equation A2. Second, opposite direction acquisitions were rescaled by a free model parameter, in order to allow for potential gradient instabilities inducing differences between the directions and their opposite directions. Third, models were generated dynamically during fitting by randomly selecting up to four hindered components and up to three restricted components. One isotropic component was always included.

The model was first fitted to half of the diffusion‐weighted data (randomly selected), after which outliers were rejected (>2.5σ). Thereafter a second fit was performed. Both fit steps assumed Gaussian noise and utilized the ‘lsqcurvefit’ function in Matlab. The procedure was repeated 100 times for different randomly generated models.

To prepare for submission of the results, only the models that best predicted the hidden half of the data were selected, after which the median of the selected predictions was used for the final prediction.

A.1.3. Scherrer (Harvard, USA): distribution of anisotropic microstructural environments in diffusion compartment imaging (DIAMOND)

DIAMOND models the set of tissue compartments in each voxel by a finite sum of unimodal continuous distributions of diffusion tensors. This corresponds to a hybrid tissue model that combines biophysical and statistical modelling. As described by Scherrer et al.,16 the DW signal S k for a gradient vector g k and b value b k is modelled by

Sk=S0j=0Nfj1+bkgkTDj0gkκjκj

where S 0 is the non‐attenuated signal, N is the number of compartments, f j the relative fraction of occupancy of the jth compartment and κ j and Dj0 are respectively the concentration and expectation of the jth continuous tensor distribution. DIAMOND enables assessment of compartment‐specific diffusion characteristics such as the compartment FA (cFA), compartment RD (cRD) and compartment MD (cMD). It also provides a novel measure of microstructural heterogeneity for each compartment.

The estimation of a continuous distribution of diffusion tensors requires DW data acquired with the same timing parameters δ and Δ.16 To compare DIAMOND with other models in this dataset, we fitted one DIAMOND model separately for each {δ, Δ} group (i.e. for each TE group), leading to 12 DIAMOND models. One shell was missing in each TE group; we predicted its signal using the corresponding DIAMOND model. The model estimation was achieved as follows. We first computed the mean and standard deviation of S 0( μS0 and σS0) within each TE group and discarded DW signals with intensity larger than μS0+3σS0(simple artefact correction). We then estimated DIAMOND parameters as described in Scherrer et al.,16 considering Gaussian noise and cylindrical anisotropic compartments. For the genu, we considered a model with one freely diffusing and one anisotropic compartment; for the fornix, we considered a model with one freely diffusing compartment and two anisotropic compartments.

A.1.4. Ferizi_1 and Ferizi_2 (UCL, England)

This submission uses two three‐compartment models, as described in previous studies.39, 41 These models consist of (1) either a Bingham distribution of sticks or a cylinder for the intracellular compartment; (2) a diffusion tensor for the extracellular compartment; (3) an isotropic cerebrospinal fluid (CSF) compartment. The T 2 relaxation element is fitted beforehand to the (variable echo time) b=0 measurements. The signal model is as follows:

S=S0˜[fiexp(TET2i)Si+feexp(TET2e)Se+fcexp(TET2c)Sc]
(A3)

where f i, f e and f c are the weights of the intracellular, extracellular and third normalized compartment signals S i, S e and S c, respectively; the values of compartmental T 2 are indexed similarly; S0˜ is the proton density signal (which is TE‐independent and obtained from fitting to the b=0 signal). These models emerged from previous studies.41, 44 Here, however, a single white matter T 2 and separate compartmental diffusivities are additionally fitted.

There is a two‐stage model fitting procedure. The first step estimates the T 2 decay rate of tissue separately in each voxel, by fitting a bi‐exponential model to the b=0 intensity as a function of TE, in which one component is from tissue and the other from CSF. A preliminary analysis of voxels fully inside WM regions shows no significant departure from mono‐exponential decay; equal T 2 values are then assumed within the intra and extracellular compartments. When fitting the bi‐exponential model, the value of T 2 in CSF is fixed to 1000 ms (a more precise value of CSF is unlikely to be estimated with this protocol). Thus, for each voxel, the volume fraction of CSF, S0˜ and T 2 of the tissue are estimated. These three estimates are then fixed for all the subsequent model fits. Then, each model is fitted using the Levenberg–Marquardt algorithm with an offset Gaussian noise model. The model parameters obtained were similar to earlier estimates obtained using the full dataset,42 differing by between 5–10% from the original.

A.1.5. Poot (Erasmus, the Netherlands)

This submission uses a three‐compartment model, with for each compartment a different complexity of the diffusion model and an individual T 2 value. This model was developed specifically for the ISBI WM challenge and is the result of iteratively visualizing different projections of the residuals and trying to infer the maximum complexity that the rich data supports.

The first compartment models isotropic diffusion and, through the initialization procedure, it captures the fast diffusion components. The second compartment is modelled by a second‐order (diffusion) tensor and models intermediate diffusion strengths. The third compartment is model‐free, as the ADC is estimated for each direction independently. Each compartment additionally has an individual T 2 value and signal intensity at b=0, T E=0 (which could easily be translated into volume fractions). Hence, the complete model of a voxel in image j is given by

Sj(θ)=i=13AieTEjR2,iebjADCj,i=i=13eMi,jθ
(A4)

where S j is the predicted signal intensity of image j, A i is the non‐diffusion weighted signal intensity of compartment i at zero T E, T E is the echo time, R 2 is the reciprocal of the T 2 relaxation time of compartment i, b=(Δ−δ/3)δ 2|G|2 γ 2, with γ=42.5781 MHz/T, A D C j,1=c, ADCj,2=gjTDgj, ADCj,3=dhjT, where d is a vector with the ADC value of each orientation group and h j is a vector that selects the orientation group to which image j belongs (90 groups in total). Note that h j has at most one non‐zero element and that element has a value of one. As displayed in the rightmost part of Equation A4, the model can be written as a multiplication of matrices M i, containing all rows M i,j, with θ=[lnA1, R 2,1, c, lnA2, R 2,2, D 11, D 12, D 13, D 22, D 23, D 33, lnA3, R 2,3, d]T, which combines all 103 parameters into a single parameter vector. All parameters are estimated simultaneously from the 3311 measurements provided per voxel by a maximum‐likelihood estimator that assumes a Rician distribution of the measurements and simultaneously optimizes the noise level.59 Finally, the signal intensities of the ‘unseen’ data are predicted by substituting the estimate into Equation A4.

A.1.6. Rokem (Standford, USA): a restriction‐spectrum sparse fascicle model (RS‐SFM)

The sparse fascicle model (SFM)43 is a member of the large family of models that account for the diffusion MRI signal in the white matter as a combination of signals due to compartments corresponding to different axonal fibre populations (fascicles) and other parts of the tissue. Model fitting proceeds in two steps. First, an isotropic component is fitted. We model the effects of both the measurement echo time (TE) and the measurement b value on the signal. These are fitted as a log(TE)‐dependent decay with a low‐order polynomial function and a b‐value‐dependent multi‐exponential decay (also including an offset to account for the Rician noise floor). The residuals from the isotropic component are then deconvolved with the perturbations in the signal due to a set of fascicle kernels, each modelled as a radially symmetric (λ 2=λ 3) diffusion tensor. The putative kernels are distributed in a dense sampling grid on the sphere. Furthermore, restriction spectrum imaging (RSI)60 is used to extend the model, by adding a range of fascicle kernels at each sampling point, with different axial and radial diffusivities, capturing diffusion at different scales. To restrict the number of anisotropic components (fascicles) in each voxel and to prevent overfitting, the RS‐SFM model employs the Elastic N et algorithm (EN),79 which applies a tunable combination of L1 and L2 regularization on the weights of the fascicle kernels. We used elements of the SFM implemented in the dipy software library80 and the EN implemented in scikit‐learn.81 In addition, to account for differences in SNR, we implemented a weighted least‐squares strategy, whereby each signal's contribution to the fit was weighted by its TE, as well as the gradient strength used. EN has two tuning parameters, determining (1) the ratio of L1‐to‐L2 regularization and (2) the weight of the regularization relative to the least‐squares fit to the signal. To find the proper values of these parameters, we employed k‐fold cross‐validation,43 leaving out one shell of measurement in each iteration for cross‐validation. We determined that the tuning parameters with the lowest Least Squares Error (LSE)18 provide an almost‐even balance of L1 and L2 penalty with weak overall regularization. Because of the combination of a dense sampling grid (362 points distributed on the sphere) and multiple restriction kernels (45 per sampling point), the maximal number of parameters for the model is approximately 16 300, more than the number of data points. However, because regularization is employed, the effective number of parameters is much smaller, resulting in an active set of approximately 20 regressors.82 We have made the code to reproduce our results fully available at https://arokem.github.io/ISBI2015.

A.1.7. Eufracio (CIMAT, Mexico): diffusion basis functions for multi–shell scheme

This model is based on the Diffusion Basis Functions (DBF) model,74 a discrete version of the Gaussian Mixture Model for the sphere: ŝi=j=1mαjϕij+ϵ, with ŝi=si/s0, ϕij=exp(bqiTTjqi) and Tj=(χ1vjvjT+χ2I). The DBF model is reformulated by substituting ϕ ij and T j: ŝi=j=1mαjexpbiχ2giTgiexpbiχ1(vjTgi)2+ϵ. The first exponential can be defined as a scale factor that depends on the b values, βi=exp(biχ2qiTqi). In this way, the β i factors are associated with different b values, so the new model includes information for multi‐shell schemes. The coefficients α and the shell scale factor β are computed by solving the optimization problem:

minα,βcf(α,βc;λα,λβ)=BΦ˜αS22+λαα1+λββc0βc22s.t.1Tα=1,α0
(A5)

where B= diag(β c),

βc=1#CiCexp(biχ^2(qiTqi))

and C is the set of indices grouped by different b values (#C is the number of elements in it). The regularization term weighted by λ α demands sparseness and the term weighted by λ β prevents overfitting. The problem in Equation A5 is solved in three steps. First, the active atoms are predicted (α i>0) with α˜=argminαf(α,βc;λα,λβ). Second, the active atoms are corrected with α=arg min{αi}:α˜i>0f(α,βc;0,λβ). Finally, the factors β c are updated with βc=argminβcf(α,βc;λα,λβ). To solve each step, the active sets algorithm for quadratic programming is used.

To train the model for the WMM'15 data, Equation A5 is solved for each voxel with the training data to find the optimal weights α j and scale factors β c that best reproduce the training data. For this challenge, the β c factors are grouped by the 36 training shells and the method parameters are set by hand: λ α=0.5, λ β=0.02, χ 1=9.5×10−4and χ 2=5×10−5. To predict the unseen signal at each voxel, the reformulated model is used with the optimal weights α j and the 12 scale factors for the unseen β c are calculated by interpolation with the 36 optimal β c of the training data.

A.1.8. Loya‐Olivas_1 and Loya‐Olivas_2 (CIMAT, Mexico): Linear Acceleration of Sparse and Adaptive Diffusion Dictionary (LASADD)

LASADD is a multi–tensor based technique to adapt dynamically the diffusion functions (DFs) dictionary to a DW–MRI signal.83, 84 The method changes the size and orientation of relevant diffusion tensors (DTs). The optimization algorithm uses a special DT expression and assumptions to reduce the computational cost.

The one‐compartment version (LASADD–1C) is based on a DBF multi–tensor model:74 si=j=1nαjϕi,j, where si=si/s0i, ϕi,j=expbigiTTjgi, α j>0 and j=1nαj=1. LASADD expresses the DT as

Tj=χ1jvjvjT+χ2jI,
(A6)

where χ{1,2}j are scalars associated with the eigenvalues, v j is the principal diffusion direction (PDD) and I is the identity matrix. The algorithm iterates three steps, like Aranda et al.:85, 86 Predict, Correct, and Generate, until convergence. Prediction selects the relevant DFs using LASSO to regulate the number to choose. Correction adjusts the volume fraction, size and orientation of the DTs. Taking advantage of the DT expression and Taylor first‐order series approximation of the exponential, the optimizations are reduced to bounded least‐squares problems, which are solved by a Projected Gauss–Seidel scheme. Generation controls the overestimation of fibres by adding to the basis the DTs resulting from combining two and three DFs for the new iteration.

An extra refinement to the computed results, named LASADD–3C, splits each detected DF into three compartments:87 intracellular (IC), extracellular (EC) and CSF. The multi‐tensor model is si=j=1nαjICψi,j+j=1nαjECθi,j+αCSFωi with j=1nαjIC+αjEC+αCSF=1. The ψ i,j models the directional IC compartment diffusion for each fibre bundle using TjIC=χ0jvjvjT. The EC compartment with hindered diffusion uses the representation (A6) for θ i,j. The isotropic diffusion ω i uses T CSF=χ 3 I. This stage keeps the PDDs fixed and only adjusts the α's and χ's of the three compartments.

The parameters of the models were estimated using the training dataset: the b values using the equation by Stejskal and Tanner88 and the S 0 values as the median of the gradient–free signals with equal echo time per voxel. The initial basis comprises 33 PDDs distributed in the unitary sphere. The bounds χ {0,1}[set membership][1,39]×10−4 and χ{2,3}[1,9]×104mm2/s and the LASSO regularization parameter (equals 1.7) was tuned by hand such as to provide the minimum error. The best multi‐tensorial model for both algorithms was used for each voxel to predict the corresponding unseen data.

A.2.  Signal models

A.2.1. Alipoor (Chalmers, Sweden)

The diffsuion MRI signal is modelled as a fourth‐order symmetric tensor as proposed by Özarslan and Mareci.89 Let g i=[x i y i z i] and ai=[zi44yizi36yi2zi24yi3ziyi44xizi312xiyizi212xiyi2zi4xiyi3 6xi2zi212xi2yizi 6xi2yi2 4xi3zi4xi3yixi4]T be a gradient encoding direction and the corresponding design vector, respectively. The diffusion signal is then described by

S(gi)=S0expTET2exp(btTai)
(A7)

where S(g i) is the measured signal when the diffusion sensitizing gradient is applied in the direction g i, S 0 is the observed signal in the absence of such a gradient, b is the diffusion‐weighting factor and tR15 contains the distinct entries of a fourth‐order symmetric tensor. Note that d(g i,t)=d(g i) is used for simplification. Given measurements in N>15 different directions, the least‐squares (LS) estimate of the diffusion tensor is t^=(GTG)1GTy, where G is an N×15 matrix defined by G=[a 1 a 2(...)a N]T and yi=b1ln(S(gi)/S0). We use the weighted LS tensor estimation method in Alipoor et al.90 to mitigate the influence of outliers.

To estimate the diffusion signal for a given acquisition protocol with T E=T E x, b=b x and δ=δ x, the two non‐diffusion‐weighted measurements with the closest T E to T E x (among measurements with δ=δ x) are used to estimate T 2 and S 0 for each voxel. Then, data from the closest shell to b x(among shells with δ=δ x) are used to estimate the tensor describing the underlying structure.

A.2.2. Sakaie–Tatsuoka–Ghosh (Cleveland, USA): an empirical approach

As the extent of q space in the dataset is unusually comprehensive, we chose a simple, generic approach to gain intuition. Visual inspection suggested use of a restricted and hindered component, each with angular variation:

Si=ATEi(fRi+(1f)exp(biDi))
(A8)

where S i is the predicted signal for a signal acquired with T E i, b i. ATEi is the median signal at a given TE with no diffusion weighting. Fitted parameters are f, the volume fraction of R i, the restricted component, and D i, the diffusivity. R i and D i are modelled as spherical harmonics with real, antipodal symmetry91 with maximum degree 4. The model has 31 fit parameters for each voxel. Data were fitted using using a non‐linear least‐squares algorithm (lsqcurvefit, MATLAB). Prior to the fit, data points with non‐zero bvalue that had a signal higher than the median of the b=0 signal plus 1.4826 times the median absolute deviation were excluded. Shells with a normalized median signal smaller than that of shells with lower bvalues were also excluded. Normalization was performed by dividing by the median of the b=0 signal with the same TE.

A.2.3. Fick (INRIA, France): a spatio‐temporal functional basis to represent the diffusion MRI signal

We use our recently proposed spatio‐temporal (3D+t) functional basis61 to represent the diffusion MRI signal simultaneously over three‐dimensional wave vector q and diffusion time τ. Based on Callaghan's theoretical model of spatio‐temporal diffusion in pores,92 our basis represents the 3D+t diffusion signal attenuation E(q,τ) as a product of a spatial and temporal functional basis as

E(q,τ)=N=0Nmax{jlm}o=0Omaxc{jlmo}Sjlm(q,us)To(τ,ut)
(A9)

where T o is our temporal basis with basis order o and S jlm is the spatial isotropic MAP‐MRI basis29 with radial and angular basis orders j, l and m. Here, Nmax and Omax are the maximum spatial and temporal order of the bases, which can be chosen independently. We formulate the bases themselves as

Sjlm(q,us)=4πil(2π2us2q2)l/2e2π2us2q2×Lj1l+1/2(4π2us2q2)Ylm(u)To(τ,ut)=exp(utτ/2)Lo(utτ)
(A10)

with u s and u t the spatial and temporal scaling functions, Ylm the spherical harmonics and L o a Laguerre polynomial. We calculate the spatial scaling u s by fitting an isotropic tensor to the TE‐normalized signal attenuation E(q,·) for all q. Similarly, we compute u t by fitting an exponential eutτ/2 to E(·,τ) for all τ. We fit our basis using Laplacian‐regularized least squares in the following steps. We first denote Ξ i(q,τ,u s,u t)=S jlm(i)(q,u s)T o(i)(τ,u t) with i[set membership]{1…N coef}, with N coef the number of fitted coefficients. We then construct a design matrix QRNdata×Ncoef with Qik=SNi(A,qk)Toi(τk,ut). The signal is then fitted as c=argminc||yQc||2+λ U(c) with y the measured signal, c the fitted coefficients and λ the weight for our analytic Laplacian regularization U(c). We used generalized cross‐validation93 to find the optimal regularization weighting λ in every voxel. In our submitted results, we used a spatial order of 8 and a temporal order of 4, resulting in 475 fitted coefficients.

A.2.4. Rivera (CIMAT, Mexico): baseline method: robust regression

We regard this very simplistic model as a baseline for other model‐based methods. It assumes as little information as possible from the diffusion signal. The vector of independent variables is x i=[g i,|G|ii,δ i,T E i,b i], containing the gradient strength g, the echo time T E and the b value b. Given signal s i, we then estimate the parameters of the linear regression model:

s=Xθ+ϵ
(A11)

where θ[set membership]R 23 is the unknown vector of coefficients, ϵ is the residual error and

X=x,x|2,Δδ,ΔTE,Δb,δTE,δb,TEb,1

is the matrix design (x|2 is obtained from squaring each element of the matrix x). To account for outliers, we estimate θ with a weighted (robust) least‐squares approach using the Lasso regularization:

θt+1=arg minθWt(Xθy)22+λθ1
(A12)

where W 0 is the identity matrix and each subsequent W is computed via

Wt+1=diag(vit+1wit+1)
(A13)

with outlier weighting in ωit+1=κ2/(κ2+(yiXθit+1)2) through κ, an arbitrary parameter that controls the outlier sensitivity. The protocol weight

vit+1=meanjΩi{wjt+1}andΩi={j:TEj=TEi,|Gj|=|Gi|}
(A14)

computes a confidence factor for the complete protocol.

Equations A12 and A13 are iterated three times. The final estimated signal is computed using (A11), using the protocol of the unseen signal.

Notes

Ferizi U, Scherrer B, Schneider T, et al. Diffusion MRI microstructure models with in vivo human brain Connectome data: results from a multi‐group comparison. NMR in Biomedicine. 2017;30:e3734 https://doi.org/10.1002/nbm.3734

Notes

* Uran Ferizi, Benoit Scherrer and Torben Schneider joint first co‐authors.

REFERENCES

1. Basser PJ, Mattiello J, LeBihan D. Estimation of the effective self‐diffusion tensor from the NMR spin echo. J Magn Reson, Ser B. 1994;103(3):247–254. [PubMed]
2. Assaf Y, Pasternak O. Diffusion tensor imaging (DTI)‐based white matter mapping in brain research: A review. J Mol Neurosci: MN. 2008;34(1):51–61. [PubMed]
3. Rose SE, Chen F, Chalk JB, et al. Loss of connectivity in Alzheimer's disease: An evaluation of white matter tract integrity with colour coded MR diffusion tensor imaging. J Neurol, Neurosurg Psychiatry. 2000;69(4):528–530. [PubMed]
4. Werring D, Brassat D, Droogan A, et al. The pathogenesis of lesions and normal‐appearing white matter changes in multiple sclerosis. Brain. 2000;123(8):1667–1676. [PubMed]
5. Price S, Jena R, Burnet N, et al. Improved delineation of glioma margins and regions of infiltration with the use of diffusion tensor imaging: An image‐guided biopsy study. Amer J Neuroradiol. 2006;27(9):1969–1974. [PubMed]
6. Mitra PP, Sen PN, Schwartz LM, Le Doussal P. Diffusion propagator as a probe of the structure of porous media. Phys Rev Lett. 1992;68(24):3555. [PubMed]
7. Behrens T, Woolrich M, Jenkinson M, et al. Characterization and propagation of uncertainty in diffusion‐weighted MR imaging. Magn Reson Med. 2003;50(5):1077–1088. [PubMed]
8. Behrens T, Berg HJ, Jbabdi S, Rushworth M, Woolrich M. Probabilistic diffusion tractography with multiple fibre orientations: What can we gain? NeuroImage. 2007;34(1):144–155. [PubMed]
9. Assaf Y, Basser PJ. Composite hindered and restricted model of diffusion (CHARMED) MR imaging of the human brain. NeuroImage. 2005;27(1):48–58. [PubMed]
10. Jespersen SN, Kroenke CD, Østergaard L, Ackerman JJH, Yablonskiy DA,Modeling dendrite density from magnetic resonance diffusion measurements. Neuroimage. 2007;34(4):1473–1486. [PubMed]
11. Jespersen SN, Bjarkam CR, Nyengaard JR, et al. Neurite density from magnetic resonance diffusion measurements at ultrahigh field: Comparison with light microscopy and electron microscopy. Neuroimage. 2010;49(1):205–216. [PubMed]
12. Alexander DC, Hubbard PL, Hall MG, et al. Orientationally invariant indices of axon diameter and density from diffusion MRI. NeuroImage. 2010;52(4):1374–1389. [PubMed]
13. Dyrby TB, Hall MG, Ptito M, et al. Contrast and stability of the axon diameter index from microstructure imaging with diffusion MRI. Magn Reson Med. 2013;70(3):711–721. [PubMed]
14. Sotiropoulos SN, Behrens TE, Jbabdi S. Ball and rackets: Inferring fiber fanning from diffusion‐weighted MRI. NeuroImage. 2012;60(2):1412–1425. [PubMed]
15. Zhang H, Schneider T, Wheeler‐Kingshott CA, Alexander DC. NODDI: Practical in vivo neurite orientation dispersion and density imaging of the human brain. NeuroImage. 2012;61(4):1000–1016. [PubMed]
16. Scherrer B, Schwartzman A, Taquet M, Sahin M, Prabhu SP, Warfield SK. Characterizing brain tissue by assessment of the distribution of anisotropic microstructural environments in diffusion‐compartment imaging (DIAMOND). Magn Reson Med. 2016;76(3):963–977. https://doi.org/10.1002/mrm.25912 [PubMed]
17. Kaden E, Kelm ND, Carson RP, Does MD, Alexander DC. Multi‐compartment microscopic diffusion imaging. NeuroImage. 2016;139:346–359. [PubMed]
18. Panagiotaki E, Schneider T, Siow B, Hall MG, Lythgoe MF, Alexander DC. Compartment models of the diffusion MR signal in brain white matter: A taxonomy and comparison. NeuroImage. 2012;59(3):2241–2254. [PubMed]
19. Stanisz GJ, Wright GA, Henkelman RM, Szafer A. An analytical model of restricted diffusion in bovine optic nerve. Magn Reson Med. 1997;37(1):103–111. [PubMed]
20. Nilsson M, van Westen D, Ståhlberg F, Sundgren PC, Lätt J. The role of tissue microstructure and water exchange in biophysical modelling of diffusion in white matter. Magn Reson Mater Phys, Biol Med. 2013;26(4):345–370. [PMC free article] [PubMed]
21. Liu C, Bammer R, Moseley ME. Generalized diffusion tensor imaging (GDTI): A method for characterizing and imaging diffusion anisotropy caused by non‐Gaussian diffusion. Isr J Chem. 2003;43(1‐2):145–154.
22. Tuch DS. Q‐ball imaging. Magn Reson Med. 2004;52(6):1358–1372. [PubMed]
23. Jensen JH, Helpern JA, Ramani A, Lu H, Kaczynski K. Diffusional kurtosis imaging: The quantification of non‐Gaussian water diffusion by means of magnetic resonance imaging. Magn Reson Med. 2005;53(6):1432–1440. [PubMed]
24. Descoteaux M, Angelino E, Fitzgibbons S, Deriche R. Regularized, fast, and robust analytical Q‐ball imaging. Magn Reson Med. 2007;58(3):497–510. [PubMed]
25. Assemlal HE, Tschumperlé D, Brun L. Efficient and robust computation of pdf features from diffusion MR signal. Med Image Anal. 2009;13(5):715–729. [PubMed]
26. Yablonskiy DA, Sukstanskii AL. Theoretical models of the diffusion weighted MR signal. NMR Biomed. 2010;23(7):661–681. [PubMed]
27. Kiselev VG. The cumulant expansion: An overarching mathematical framework for understanding diffusion NMR. In: Derek KJ, ed. Diffusion MRI: Theory, Methods, and Applications OUP USA; 2011:152–168.
28. Özarslan E, Koay C, Basser P. Simple harmonic oscillator based estimation and reconstruction for one‐dimensional q‐space MR. In: Proc. Intl. Soc. Mag. Reson. Med, Vol. 16; 2008:35.
29. Özarslan E, Koay CG, Shepherd TM, et al. Mean apparent propagator (MAP) MRI: A novel diffusion imaging method for mapping tissue microstructure. NeuroImage. 2013;78:16–32. [PubMed]
30. Tournier JD, Calamante F, Gadian DG, Connelly A. Direct estimation of the fiber orientation density function from diffusion‐weighted MRI data using spherical deconvolution. Neuroimage. 2004;23(3):1176–1185. [PubMed]
31. Alexander DC. Maximum entropy spherical deconvolution for diffusion MRI. Biennial International Conference on Information Processing in Medical Imaging: 19th International Conference, IPMI 2005. Glenwood Springs, CO, USA. Berlin: Springer; 2005:76–87. [PubMed]
32. Anderson AW. Measurement of fiber orientation distributions using high angular resolution diffusion imaging. Magn Reson Med. 2005;54(5):1194–1206. [PubMed]
33. Jian B, Vemuri BC. Multi‐fiber reconstruction from diffusion MRI using mixture of Wisharts and sparse deconvolution Information Processing in Medical Imaging. Berlin: Springer; 2007:384–395. [PMC free article] [PubMed]
34. Jian B, Vemuri BC, Özarslan E, Carney PR, Mareci TH. A novel tensor distribution model for the diffusion‐weighted MR signal. NeuroImage. 2007;37(1):164–176. [PubMed]
35. Sakaie KE, Lowe MJ. An objective method for regularization of fiber orientation distributions derived from diffusion‐weighted MRI. NeuroImage. 2007;34(1):169–176. [PubMed]
36. Dell'Acqua F, Rizzo G, Scifo P, Clarke RA, Scotti G, Fazio F. A model‐based deconvolution approach to solve fiber crossing in diffusion‐weighted MR imaging. IEEE Trans Biomed Eng. 2007;54(3):462–472. [PubMed]
37. Jbabdi S, Sotiropoulos S, Savio A, Graña M, Behrens T. Model‐based analysis of multishell diffusion MR data for tractography: How to get over fitting problems. Magn Reson Med. 2012;68(6):1846–1855. [PubMed]
38. Rathi Y, Michailovich O, Laun F, Setsompop K, Grant PE, Westin CF. Multi‐shell diffusion signal recovery from sparse measurements. Med Image Anal. 2014;18(7):1143–1156. [PubMed]
39. Ferizi U, Schneider T, Panagiotaki E, et al. A ranking of diffusion MRI compartment models with in vivo human brain data. Magn Reson Med. 2014;72(6):1785–1792. [PubMed]
40. Kaden E, Kruggel F, Alexander DC. Quantitative mapping of the per‐axon diffusion coefficients in brain white matter. Magn Reson Med. 2016;75(4):1752–1763. https://doi.org/10.1002/mrm.25734 [PubMed]
41. Ferizi U, Schneider T, Tariq M, Wheeler‐Kingshott C, Zhang H, Alexander D. The importance of being dispersed: A ranking of diffusion MRI models for fibre dispersion using in vivo human brain data In: Mori K, editor; , Sakuma I, editor; , Sato Y, editor; , Barillot C, editor; , Navab N, editor. , eds. Medical Image Computing and Computer‐Assisted Intervention – MICCAI 2013, Lecture Notes in Computer Science, vol. 8149. Berlin: Springer; 2013:74–81. [PubMed]
42. Ferizi U, Schneider T, Zhang H, Wheeler‐Kingshott CA, Alexander DC. White matter compartment models for in vivo diffusion MRI at 300 mT/m. Neuroimage. 2015;118:468–483. https://doi.org/10.1016/j.neuroimage.2015.06.027 [PubMed]
43. Rokem A, Yeatman JD, Pestilli F, et al. Evaluating the accuracy of diffusion MRI models in white matter. PLoS One. 2015;10(4):e0123272. [PubMed]
44. Ferizi U. Compartment Models and Model Selection for in‐vivo Diffusion‐MRI of Human Brain White Matter [PhD thesis]. London: University College London; 2014.
45. Klein A, Andersson J, Ardekani BA, et al. Evaluation of 14 nonlinear deformation algorithms applied to human brain MRI registration. NeuroImage. 2009;46(3):786–802. [PubMed]
46. Mendrik A, Vincken K, Kuijf H, et al. MRBrainS challenge: Online evaluation framework for brain image segmentation in 3T MRI scans. Comput Intell Neurosci. 2015;2015:813696 https://doi.org/10.1155/2015/813696 [PubMed]
47. Daducci A, Canales‐Rodriguez EJ, Descoteaux M, et al. Quantitative comparison of reconstruction methods for intra‐voxel fiber recovery from diffusion MRI. IEEE Trans Med Imaging. 2014;33(2):384–399. [PubMed]
48. Fillard P, Descoteaux M, Goh A, et al. Quantitative evaluation of 10 tractography algorithms on a realistic diffusion MR phantom. Neuroimage. 2011;56(1):220–234. [PubMed]
49. Ning L, Laun F, Gur Y, et al. Sparse reconstruction challenge for diffusion MRI: Validation on a physical phantom to determine which acquisition scheme and analysis method to use? Med Image Anal. 2015;26(1):316–331. [PubMed]
50. Pujol S, Wells W, Pierpaoli C, et al. The DTI challenge: Toward standardized evaluation of diffusion tensor imaging tractography for neurosurgery. J Neuroimaging. 2015;25(6):875–882. [PubMed]
51. Setsompop K, Kimmlingen R, Eberlein E, et al. Pushing the limits of in vivo diffusion MRI for the Human Connectome project. NeuroImage. 2013;80:220–233. [PubMed]
52. McNab JA, Edlow BL, Witzel T, et al. The human connectome project and beyond: Initial applications of 300 mT/m gradients. NeuroImage. 2013;80:234–245. [PubMed]
53. Duval T, McNab J, Setsompop K, et al. In vivo estimation of axon diameter in the human spinal cord using 300 mT/m gradients. International Society of Magnetic Resonance in Medicine (ISMRM) 24th Annual Meeting & Exhibition Milan, Italy; 2014.
54. Huang SY, Nummenmaa A, Witzel T, et al. The impact of gradient strength on in vivo diffusion MRI estimates of axon diameter. NeuroImage. 2015;106:464–472. [PubMed]
55. Burnham KP, Anderson DR. Model Selection and Multimodel Inference: A Practical Information‐Theoretic Approach. Berlin: Springer; 2002.
56. Jenkinson M, Bannister P, Brady M, Smith S. Improved optimization for the robust and accurate linear registration and motion correction of brain images. NeuroImage. 2002;17(2):825–841. [PubMed]
57. Jones DK, Basser PJ. Squashing peanuts and smashing pumpkins': How noise distorts diffusion‐weighted MR data. Magn Reson Med. 2004;52(5):979–993. [PubMed]
58. Nilsson M, Lätt J, Ståhlberg F, Westen D, Hagslätt H. The importance of axonal undulation in diffusion MR measurements: A Monte Carlo simulation study. NMR Biomed. 2012;25(5):795–805. [PubMed]
59. Poot D, Klein S. Detecting statistically significant differences in quantitative MRI experiments, applied to diffusion tensor imaging. IEEE Trans Med Imaging. 2015;34(5):1164–1176. [PubMed]
60. White NS, Leergaard TB, D'Arceuil H, Bjaalie JG, Dale AM. Probing tissue microstructure with restriction spectrum imaging: Histological and theoretical validation. Hum Brain Mapp. 2013;34(2):327–346. [PubMed]
61. Fick R, Wassermann D, Pizzolato M, Deriche R. A unifying framework for spatial and temporal diffusion in diffusion MRI International Conference on Information Processing in Medical Imaging. Berlin: Springer; 2015:167–178. [PubMed]
62. Fieremans E, Jensen JH, Helpern JA. White matter characterization with diffusional kurtosis imaging. Neuroimage. 2011;58(1):177–188. [PubMed]
63. Novikov DS, Veraart J, Jelescu IO, Fieremans E. Mapping orientational and microstructural metrics of neuronal integrity with in vivo diffusion MRI. arXiv preprint arXiv:1609.09144. 2016 Sep 28.
64. Jelescu IO, Veraart J, Adisetiyo V, Milla SS, Novikov DS, Fieremans E. One diffusion acquisition and different white matter models: How does microstructure change in human early development based on WMTI and NODDI? NeuroImage. 2015;107:242–256. [PubMed]
65. Jelescu IO, Veraart J, Fieremans E, Novikov DS. Degeneracy in model parameter estimation for multi‐compartmental diffusion in neuronal tissue. NMR Biomed. 2016;29(1):33–47. [PubMed]
66. Ghosh A, Zhang H, Alexander DC. To be dispersed or not to be dispersed: A study using HCP data International Society of Magnetic Resonance in Medicine (ISMRM) 24th Annual Meeting & Exhibition. Singapore; 2016.
67. Burcaw LM, Fieremans E, Novikov DS. Mesoscopic structure of neuronal tracts from time‐dependent diffusion. NeuroImage. 2015;114:18–37. [PubMed]
68. Barazany D, Basser PJ, Assaf Y. In vivo measurement of axon diameter distribution in the corpus callosum of rat brain. Brain. 2009;132(5):1210–1220. [PubMed]
69. Fieremans E, Burcaw LM, Lee HH, Lemberskiy G, Veraart J, Novikov DS. In vivo observation and biophysical interpretation of time‐dependent diffusion in human white matter. NeuroImage. 2016;129:414–427. [PubMed]
70. Guglielmetti C, Veraart J, Roelant E, et al. Diffusion kurtosis imaging probes cortical alterations and white matter pathology following cuprizone induced demyelination and spontaneous remyelination. Neuroimage. 2016;125:363–377. [PubMed]
71. Kelm ND, West KL, Carson RP, Gochberg DF, Ess KC, Does MD. Evaluation of diffusion kurtosis imaging in ex vivo hypomyelinated mouse brains. NeuroImage. 2016;124:612–626. [PubMed]
72. Yablonskiy DA, Bretthorst LG, Ackerman JJ. Statistical model for diffusion attenuated MR signal. Magn Reson Med. 2003;50(4):664–669. [PubMed]
73. Axer H, Axer M, Krings T, Graf Diedrich v, Keyserlingk D. Quantitative estimation of 3‐d fiber course in gross histological sections of the human brain using polarized light. J Neurosci Methods. 2001;105(1):121–131. [PubMed]
74. Ramirez‐Manzanares A, Rivera M, Vemuri BC, Carney P. Diffusion basis functions decomposition for estimating white matter intravoxel fiber geometry. IEEE Trans Med Imaging. 2007;26(8):1091–1102. [PubMed]
75. Van Gelderen P, DesPres D, van Zijl P, Moonen CT. Evaluation of restricted diffusion in cylinders. phosphocreatine in rabbit leg muscle. J Magn Reson, Ser B. 1994;103(3):255–260. [PubMed]
76. Stone M. Cross-validatory choice and assessment of statistical predictions. J R Stat Soc Ser B (Methodological). 1974;1:111–147.
77. Efron B. Bootstrap methods: Another look at the jackknife. Ann Stat. 1979:1–26.
78. Vangelderen P, DesPres D, Vanzijl P, Moonen C. Evaluation of restricted diffusion in cylinders. phosphocreatine in rabbit leg muscle. J Magn Reson, Ser B. 1994;103(3):255–260. [PubMed]
79. Zou H, Hastie T. Regularization and variable selection via the elastic net. J R Stat Soc Ser B Stat Methodol. 2005;67(2):301–320.
80. Garyfallidis E, Brett M, Amirbekian B, et al. DIPY, a library for the analysis of diffusion MRI data. Front Neuroinform. 2014;8:8. [PubMed]
81. Pedregosa F, Varoquaux G, Gramfort A, et al. SCIKIT‐LEARN: Machine learning in Python. J Mach Learning Res. 2011;12:2825–2830.
82. Zou H, Hastie T, Tibshirani R, et al. On the ‘degrees of freedom’ of the lasso. Ann Stat. 2007;35(5):2173–2192.
83. Loya‐Olivas A, Aranda R, Rivera M. LASADD: linear acceleration method for adapting diffusion dictionaries. In: International Society of Magnetic Resonance in Medicine (ISMRM) 23rd Annual Meeting and Exhibition; 2015; Toronto, Ontario, Canada.
84. Loya AK. Algoritmos para solución de modelos no lineales en DW‐MRI [Master's thesis]. Guanajuato, Mexico: Centro de Investigación en Matemáticas; 2015.
85. Aranda R, Ramirez‐Manzanares A, Rivera M. Recovering detailed intra‐voxel white matter structure by using an adaptive diffusion dictionary. In: International Society of Magnetic Resonance in Medicine (ISMRM) 23rd Annual Meeting and Exhibition; Toronto, Ontario, Canada:2015.
86. Aranda R, Ramirez‐Manzanares A, Rivera M. Sparse and Adaptive Diffusion Dictionary (SADD) for recovering intra‐voxel white matter structure. Med Image Anal. 2015;26(1):243–255. [PubMed]
87. Sherbondy AJ, Rowe MC, Alexander DC. Microtrack: An algorithm for concurrent projectome and microstructure estimation Medical Image Computing and Computer‐Assisted Intervention – MICCAI 2010, Lecture Notes in Computer Science, vol. 6361. Berlin: Springer; 2010:183–190. [PubMed]
88. Stejskal EO, Tanner JE. Spin diffusion measurements: Spin echoes in the presence of a time‐dependent field gradient. J Chem Phys. 1965;42(1):288–292. https://doi.org/10.1063/1.1695690
89. Özarslan E, Mareci TH. Generalized diffusion tensor imaging and analytical relationships between diffusion tensor imaging and high angular resolution diffusion imaging. Magn Reson Med. 2003;50(5):955–965. [PubMed]
90. Alipoor M, Gu IY, Mehnert AJ, Lilja Y, Nilsson D. On high order tensor‐based diffusivity profile estimation Engineering in Medicine and Biology Society (EMBC), 2013 35th Annual International Conference of the IEEE. Osaka, Japan; 2013:93–96. [PubMed]
91. Alexander D, Barker G, Arridge S. Detection and modeling of non‐Gaussian apparent diffusion coefficient profiles in human brain data. Magn Reson Med. 2002;48(2):331–340. [PubMed]
92. Callaghan PT. Pulsed‐gradient spin‐echo NMR for planar, cylindrical, and spherical pores under conditions of wall relaxation. J Magn Reson, Ser A. 1995;113(1):53–59.
93. Craven P, Wahba G. Smoothing noisy data with spline functions. Numer Math. 1978;31(4):377–403.

Articles from Wiley-Blackwell Online Open are provided here courtesy of Wiley-Blackwell, John Wiley & Sons