Abstract
The diffusion-weighted-dependent attenuation of the MRI signal E(b) is extremely sensitive to microstructural features. The aim of this study was to determine which mathematical model of the E(b) signal most accurately describes it in the brain. The models compared were the monoexponential model, the stretched exponential model, the truncated cumulant expansion (TCE) model, the biexponential model, and the triexponential model. Acquisition was performed with nine b-values up to 2500 s/mm2 in 12 healthy volunteers. The goodness-of-fit was studied with F-tests and with the Akaike information criterion. Tissue contrasts were differentiated with a multiple comparison corrected nonparametric analysis of variance. F-test showed that the TCE model was better than the biexponential model in gray and white matter. Corrected Akaike information criterion showed that the TCE model has the best accuracy and produced the most reliable contrasts in white matter among all models studied. In conclusion, the TCE model was found to be the best model to infer the microstructural properties of brain tissue.
Introduction
Diffusion MRI (dMRI) is a powerful in vivo method for exploring biological microstructures. 1 The contrast of diffusion-weighted (DW) images is dependent on the magnitude of the diffusion-weighting factor (b-value). In clinical practice, especially for fast and early detection of ischemic stroke, DW images are generally acquired with a unique b-value (b = 1000 s/mm2). The geometric mean of DW images over three spatial directions is then used to calculate apparent diffusion coefficient (ADC) maps 2 with the help of a monoexponential model (MEM) formalism. As it uses the same limited range of b-values but is acquired for 6–60 spatial directions, 3 diffusion tensor imaging is also an MEM that includes the effect of water diffusion anisotropy 3 and leads to the calculation of a rotationally invariant form of the ADC, the mean diffusivity.
On the other hand, the higher the b-value, the more sensitive the DW images are to smaller absolute random displacements of water molecules, and therefore to the microscopic structures of brain tissue. 1 Therefore, high b-value dMRI consists of recording series of DW images with successively increasing b-values up to relatively high values (b ~2500 to ~12,000 s/mm2). The number of spatial direction measures acquired is adapted depending on the time constraints. The DW-dependent attenuation of the MRI signal obtained in each image pixel is related by Fourier transformation to the water displacement probability distribution (q-space method). In brain tissue, the standard deviation of this distribution has a characteristic length scale of tenths of μm,4,5 making it possible to investigate some microstructural changes related to cerebral insults. 5 However, Fourier transformation of the DW-dependent attenuation of the MRI signal (q-space) requires experimentally challenging conditions, whereas signal modeling is more flexible. Indeed, the DW-dependent attenuation of the MRI signal (designed as E(b)) acquired with intermediate or high b-values needs to be modeled numerically with a mathematical function describing its non-exponential decay, which is relevant to the postulated properties of the microscopic structures present in brain tissues.4,6–8 High b-value dMRI modeling and q-space are very sensitive to microscopic water displacements, highlighting white matter (WM) structures7,8 or cerebral insults in animals9–12 and in human brains,5,13 and thus having a potential clinical interest.
Two types of models are needed for the study of diffusion imaging. The first (signal models) describes empirically how the E(b) signal decays, whereas the second (tissue models) describes how the signal relates to the underlying tissue microstructure. However, for the former, several signal models have been developed to perform high b-value dMRI modeling, ie, to fit the E(b) signal4,6–8: the biexponential model (BEM),7,14–18 the triexponential model (TEM),3,7,12 the stretched exponential model (SEM),19–23 and the truncated cumulant expansion (TCE) model.1,24–26 The latter is closely related to the diffusional kurtosis imaging model.1,24–26 However, to date, no studies have compared the contrast and goodness-of-fit of all of these four different models on the same scans and the best model to be used in research imaging protocols remains a matter of debate. The aim of this study was to determine which mathematical model is the best to analyze the DW-dependent attenuation of the MRI signal of the human brain in a range of clinically accessible b-values, in terms of accuracy, contrast, and insensitivity to artifacts. Accuracy was studied with the Akaike information criterion (AIC), which is a reference statistical index of goodness-of-fit that makes it possible to test simultaneously the accuracy of different models potentially differing by their degrees-of-freedom.
Methods
Acquisition design
Subjects and patients
For the purpose of this study, we included 12 healthy human control subjects without memory dysfunction. The study protocol (ClinicalTrials.gov Identifier: NCT01686516) complied with the principles of the Declaration of Helsinki and was approved by the local ethics committee (South West France and Overseas French Departments I–II, France). Five subjects were young (26 ± 4 years) and seven were old (66 ± 6 years). To illustrate the clinical interest and the feasibility of high b-value dMRI acquisitions in hospital emergencies, a patient was scanned 4 days after a serious stroke infarct with a protocol identical to those employed for healthy subjects. Written informed consent was obtained from every subject/patient before participation.
MRI acquisitions
Magnetic resonance images were acquired with a Philips Achieva® 3 T magnet (maximum magnetic field gradient intensity: 40 mT/m) equipped with an eight-channel coil. DW spin echo images (31 slices covering the entire brain) were acquired with a pulsed gradient spin echo dMRI pulse sequence (A/8 = 34.5/25 ms, constant diffusion time td = 27 ms) with an echo-planar (EP) imaging readout (TR/TE = 4212/70 ms). DW images were recorded with an increasing range of b -values (0, 100, 200, 500, 750, 1000, 1500, 2000, and 2500 s/mm2, N b = 9) obtained by increasing the intensity of the magnetic field gradients. For each b -value, the scanner provided the geometric mean of DW images acquired in the read (R), phase (P), and slice (S) encoding directions. Spatial resolution used was 2 mm × 2 mm (matrix size 128 × 128) and slice thickness was 4 mm. Reconstruction for parallel imaging with sensitivity encoding for fast MRI (SENSE) was used with a SENSE factor of 2. Total dMRI scan duration was seven minutes, a duration compatible with the use of our protocol in hospital emergencies.
Image processing
Image registration
Otsu's method of histogram analysis 27 was used to extract the brain image from skull and background noise. To limit the effect of the patient's head motion 28 and susceptibility-related geometrical distortions of EP images, 28 Matlab® methods applying, respectively, soft linear registration and soft local, nonlinear registration of EP images were used.29–31 DW images were also visually checked for other remaining artifacts. Pixel-by-pixel division of registered DW images S(b) by S0 image (image acquired with no diffusion weighting) was performed to obtain images of the normalized DW-dependent attenuation of the MRI signal E(b) = S(b)/S0.
Segmentation and regions-of-interest (ROIs)
Segmentation accuracies of the packages FreeSurfer, SPM, and FSL have been previously compared, and for WM, a higher volumetric segmentation accuracy was found with the FSL package. 32 In this study, higher significance was given to WM compared to gray matter (GM). WM is less prone to partial volume effects of cerebrospinal fluid (CSF), so it is more reliable for unbiased diffusion model accuracy comparisons. EP images can be segmented, but care must be taken when doing so. 33 Segmentation of EP images was done with the FSL segmentation algorithm FAST (FMRIB's automatic segmentation tool), which provided binary masks that served as ROIs corresponding to CSF, GM, and WM.
Models of the E(b) signal
The parameter S0 was not considered as a free parameter in any of the models. To fit pixel-by-pixel the S(b)/S0 images, the following models were used:
Monoexponential model
This model of the E(b) signal approximates water diffusion in the tissue as Gaussian. It also has the following form for diffusion approximated as isotropic14,15:
Stretched exponential model
The SEM has the following form19–23:
Truncated cumulant expansion
The second-order TCE of the signal E(b) is a mathematical model closely related to the q-space method,24,25 but it does not use Fourier transformation as it has the following analytical formulation1,24–26:
In brain, Dapp exhibits higher values than ADC does.24–26 Apparent kurtosis Kapp is commonly interpreted as: (i) representative of the presence of restricted diffusion4,24 or (ii) representative of the presence of tissue microscopic heterogeneity.1,25 Indeed, in the latter case, it is possible to calculate the standard deviation of the local diffusivities (σ) in the following manner1,25:
The parameter σ represents the dispersion of the correlation function of the disorder existing in the biological medium, ie, σ is a parameter estimating the structural disorder frequency of barriers encountered during the diffusion of water molecules in the tissue, like those created by the spatial distribution of biomembranes. 1
Biexponential model
The BEM has the following form15–18:
This models the sum of the contribution (with no exchange) of two water pools (Ffast and Fslow), each associated, respectively, with expected “fast” and “slow” ADCs (Dfast and Dslow).15–18 The BEM has three degrees-of-freedom (Ndof = 3).
Triexponential model
To estimate the signal attenuation properties both at low
14
and high b-values,
15
three exponential functions are combined3,7,12 to express the TEM as:
Numerical implementation
E(b) images were fitted pixel-by-pixel over the full b-value range (0–2500 s/mm2) with an in-house Matlab® program implementing the models presented (MEM, SEM, TCE, BEM, TEM). Regression convergence of nonlinear systems described by equations (1), (2), (3), (5), and (6) was performed by simple least square minimization of the sum of squared residuals (SSR) of the fit, using 800 iterations of the Nelder–Mead simplex algorithm (fminsearch Matlab® function).
34
This procedure allowed the rapid calculation of the SSR maps and parametric maps corresponding to the free parameters of the models
Statistical analysis
Analysis of model contrasts
Parametric images were visually checked for artifacts. The contrast comparison between anatomically neighboring structures (ie, CSF compared to GM, GM compared to WM) for values of spatially averaged quantitative parameters
Spearman's rank correlation implemented in Matlab® was used to obtain the coefficient of correlation R between the values of different diffusion parameters measured in the WM for the 12 subjects, together with the corresponding risk probability. The significant risk probability level chosen was P < 0.05.
Analysis of goodness-of-fit
F-tests. To calculate the F-statistic associated with model comparisons,26,36,37 the approach proposed in the work of Yoshiura et al was used,
36
where the F-value is defined as:
AIC. The AIC has a particularly simple form when the E(b) signal is fitted by simple least square
38
:
Indeed, to be able to compare efficiently the accuracy of the models, an accurate penalty for models with high degrees-of-freedom needs to be taken into account. 38 The lower the AICc value, the better the model describes the experimental data obtained in the population sample. 38
Results
Contrast analysis
Parametric maps and corresponding contrasts are shown in Figure 1 (for MEM, SEM, and TCE) and Figure 2 (for BEM and TEM), and quantitative values are given in Table 1. The tissue contrasts between adjacent tissues are shown in the bar charts in Figures 1 and 2.
Parametric images (one slice) of the MEM (ADC images top) of the SEM (DDC and α maps bottom) and of the second-order TCE (Dapp, Kapp, and σ maps) in a representative (aged) subject Bar charts show the values of the free parameters of the models (mean ± the standard deviation in the population sample) Statistically significant tissue contrast differences for anatomically-neighboring tissues are also shown (*P<0 05 **P<0 01 ***P<0 001) Parametric images (one slice) of the BEM and of the TEM in a representative (aged) subject. Bar charts show the values of the free parameters of the model (mean ± the standard deviation in the population sample). Statistically significant tissue contrast differences for anatomically-neighboring tissues are also shown (*P<0.05, **P<0.01, ***P<0.001). Quantitative values of parameters of models of the DW-dependent attenuation of the MRI signal E(b), displayed as intersubject mean ± the intersubject standard deviation.

Artifacts and reproducibility
High sensitivity to CSF pulsatile motion artifacts was often observed in the SEM, BEM, and TEM parametric images (see Figs. 1 and 2, artifacts a round the ventricles in α, Ffast, Fslow, Dfast, Dslow, F2 maps).
Tissue contrasts.
ADC (Fig. 1), DDC (Fig. 1), Dfast (Fig. 2), D2 and D3 (Fig. 2) maps showed a significant intensity decrease from CSF to WM (bar charts of Figs. 1 and 2 and Table 1).
Parametric maps modeled with the BEM and the TEM, especially water fractions Ffast, Fslow, F1, F2, and F3, appeared noisy and did not exhibit any evident anatomical contrast (Fig. 2).
Maps of the parameter α (Fig. 1) showed a significant but weak contrast between GM and WM (P < 0.05). Conversely, a particularly strong and anatomically well-defined contrast (P < 0.001) appeared in Kapp maps, clearly revealing the structure of the WM and the GM (Fig. 1). Areas with higher intensities on Kapp maps appeared to be WM substructures, maybe areas of high tissular complexity (Fig. 1). In C maps, a strong contrast differentiating GM/WM tissues and CSF (P < 0.001) was observed (Fig. 1).
Contrast correlations
The correlation between the diffusion parameters obtained in WM was then examined. There was a strong correlation between DDC and ADC (R = 0.9580, P < 0.05), but not between DDC and D (R = 0.4964, P > 0.05) or between DDC and α(R = -0.1049, P ≫ 0.05).
However, there was a positive correlation between Dapp and ADC (R = 0.6130, P< 0.05), between Dapp and Kapp (R = 0.6224, P < 0.05), and between Da and σ (R = 0.9161, P < 0.05), as shown in Figure 3. The parameters extracted from BEM and TEM, especially the ADCs of the models, showed no correlation either with ADC or with Da (P > 0.05, data not shown).
Diffusion model contrast in stroke
The E(b) signal attenuation was also modeled in a stroke patient. The classical high T2-weighted signal and the decreased intensity of the ADC maps
2
were observed in the lesion (Fig. 4A and 4B, respectively). There was a similar decrease in DDC and decrease in Dapp in the entire lesion. In the stroke lesion infarcted WM (Fig. 4D), there was a strong increase in Kapp in infarcted WM (from 0.9 to 1.81) and a decrease in σ (from 0.46 to 0.42 μm2/ms). In the GM of the lesion, σ maps showed a greater decrease (from 0.61 to 0.43 μm2/ms). For BEM in infarcted GM (Fig. 4E), there were both a strong increase in Fslow (up to 69%) and a strong decrease in Dslow (from 0.57 to 0.18 μm2/ms). F1* maps representing the magnitude of the IVIM microcapillary perfusion were very noisy as in the healthy population sample, but the delineation of the infarcted area in F1* maps showed a larger extension of the lesioned area than in the ADC map (Fig. 4F, arrowhead).
Positive correlation between Dapp values and Kapp values (top), Dapp values and ADC values (middle), Dapp values and σ values (bottom), obtained for the 12 subjects in the white matter of the sample population. Spear man's coefficient of rank correlation (R) is also indicated for each linear regression (P < 0.05). an MRI slice of multi-parametric imaging in a patient with stroke (four days after the infarct), showing both grey and white matter insults in (A) T2-weighted images (T2W) acquired with TE/TR = 80/8624 ms; (B) maps of the MEM: ADC; (C) maps of the SEM; (D) maps of the second-order TCE model; (E) maps of the BEM; and (F) two maps of the TEM: F1* and D1*.

Models’ goodness-of-fit
For the sake of clarity, we present the modeled E(b) signals obtained in the WM of a representative subject (Fig. 5, left) and the corresponding SSR maps (Fig. 5, right). The SSR maps correspond to the anatomical location of the slice of the parametric maps presented in Figures 1 and 2. To be visually indicative of the compared goodness-of-fit, both the residual of the E(b) signal and the SSR maps were multiplied by Ndof, the degree-of-freedom of each model (Fig. 5). Figure 5 is also an illustration that an objective mathematical method is needed to differentiate which functions given by equations (2)–(5) are the best to describe the E(b) signal experimentally.
Left: the experimental diffusion-weighted-dependent attenuation of the MRI signal (the E(b) signal, empty triangles ▿) and the absolute value of the residual of the fit times the degree-of-freedom Ndof (dotted lines, black squares ▪, values multiplied by a factor 10) obtained in WM in a representative subject. Separate subfigures show the modeling by the MEM, by the SEM, by the second-order TCE model, by the BEM, and by the TEM, all represented by black dots (-•-). Right: SSR times Ndof maps of the models for an anatomical slice, which corresponds to the slice presented for the parametric images of the models in Figure 1 and Figure 2
F-tests
When analyzed with F-tests in all the segmented brain regions, the BEM, TEM, and TCE all performed better than the MEM for modeling the E(b) signal (P < 0.05). The BEM performed better than the TEM in all the brain areas (P < 0.05). The TCE was found to be a better model than the BEM in GM and WM (P < 0.05). The TCE and SEM have the same number of degrees-of-freedom (Ndof = 2), so they cannot be compared with the F-test method given in the work of Yoshiura et al, Equation (7). 36
AICc
AICc measured in the white matter of the entire population sample and assessing the best accuracy of the TCE model compared to other models.
In the representative fit accuracy provided, the residual values for bmax =2500 s/m2 was high for MEM, SEM, BEM, TEM but was low for the TCE model (Fig. 5), which illustrates the good convergence of this latter model for the highest b-values used in this work. However, for very low b-values (b = 100 and 200 s/mm2), a slight deviation of the TCE fit was observed (Fig. 5). Although the accuracy of the TCE model was the best in the WM, relatively high SSR values were also observed for the TCE fit in the GM area contaminated by CSF partial volume effects (Fig. 5, right).
Discussion
The main finding of the present study is that the TCE is the most appropriate mathematical function to model the DW-dependent attenuation of the MRI signal E(b) in the WM tissue. Our results showing the excellent goodness-of-fit of the TCE model in WM compared to the BEM model agree well with those in the previous study, 26 which used the geometric mean of DW images of five subjects, with N b = 16 and b-factors up to 2500 s/mm2 and compared only the BEM and TCE models. 26 The accuracy of these two models was compared using a ratio of F-values, a noise-independent method. 26 However, the F-test approach does not make it possible to compare more than two models simultaneously and efficiently. The additional strength of our methodology is that we estimated the models’ goodness-of-fit in a larger sample size and used a powerful statistical index (AICc). This allowed us to compare more than two models irrespective of their degree-of-freedom. We were thus able to compare the TCE model simultaneously with other models such as the SEM, BEM, and TEM. Our results showed an apparent discrepancy between the accuracies obtained with the F-test (Equation (7)) and AICc (equations (8) and (9)). Indeed, the latter method lead to a better accuracy of the MEM (N dof = 1) compared to the TEM (N dof = 6). This is because the Equation (9) has a term in 2N 2 dof + 4N dof whereas Equation (7) scales only in N dof , and because a larger penalty is applied to models with a high degree-of-freedom when the AICc approach is used (Table 2).
In our work, the TCE was found to be superior not only to the BEM but also to the TEM and the SEM in terms of accuracy, contrast, and insensitivity to artifacts. The TCE provided a strong contrast between GM and WM in Kapp maps, whereas α, Ffast, or Fslow and F1, F2, or F3 maps gave a poorly defined anatomical tissue contrast. The poor contrast and low reproducibility of the quantitative parameters of the BEM and the TEM and the poor fit stability of the TEM 7 may also be because of the fact that only one water pool might be present in the tissue,1,26 contrary to the assumptions of these models. It has been noted that overfitting errors may accumulate along the b -value range, especially for TEM, which is overdetermined. 7 It could be noted that the BEM is almost always used with CSF-suppressed dMRI sequences,7,16,17 thus avoiding the presence of CSF water, which could be mistaken for an additional water pool. On the contrary, the low-degree-of-freedom TCE model, which has no a priori hypothesis concerning the number of water pools, has previously been shown to be robust to infer WM properties both for acquisitions with or without CSF suppression. 39
However, for the TCE fit, relatively high SSR values were observed in the GM area contaminated by CSF partial volume effects, as has already been reported. 26 The mixing between partial volume effects of CSF water and tissue water was suspected to favor a large and non-physiologic bicompartmental attenuation of the E(b) signal. 26
A correlation between the related parameters ADC and Dapp was found in WM, and this strong linear relationship is an experimental confirmation that the TCE model adequately extends the MEM. Self correlation was also found between the TCE parameters Dapp, Kapp, and σ. This means that the parameters of the TCE model are self-consistent. In other words, there is a relation between quantity D (related to the cumulant of order one in b) and between quantity Kapp (related to the cumulant of order two in b). DDC also correlated well with ADC, but there was no correlation between DDC and α, the two SEM parameters. Furthermore, no significant correlations were found between ADC and the diffusion coefficients of BEM and TEM.
dMRI makes it possible to obtain signals related to tissue microstructure, but it has some drawbacks. In fact, the spatial resolution of EP images is lower than those of T1-weighted images (anatomical MRI). Consequently, segmentation of EP images is far from providing an accurate differentiation of the GM and of the subcortical WM like those obtained with the segmentation of T1-weighted images. 40 Owing to these intrinsic limitations of EPI-related methods, our study focused only on WM, which dominates by its large volume both the volume of the GM and the volume of the subcortical WM. Apart from this point, our study has the following limitations:
First, using more directions of diffusion measures might potentially refine the comparisons of accuracies. Nevertheless, considering that the microstructures present in tissue are of unknown nature, different tensorial forms of equations (2), 23 (5), 41 and (6) 3 with different degrees-of-freedom exist. On the other hand, for a given model, the degrees-of-freedom must be unequivocally defined in Equations (8) and (9). This further complicates the comparisons of models but will be addressed in a future study. Consequently, at present, we are using acquisition conditions very close to those in the study of Kiselev and Il'yasov, 26 where unidirectional fitting of the geometric average of DW images was used. This could potentially limit the study of WM, which has a non-null degree of anisotropy, 26 but a simple modeling design has the advantage of defining unequivocally the degrees-of-freedom of the models compared, a condition required for efficient accuracy comparisons.
Secondly, our conclusions about the accuracy of the models might be in part limited because model's accuracies were measured in a defined experimental context determined by the maximum b-value used (bmax = 2500 s/mm2). However, this b-value range is considered to be suitable for medical imaging in routine protocols24–26 and the aim of our study was to determine the best model of the E(b) signal in such conditions. The use of b-values higher than 2500 s/mm2 has the potential to better resolve the microstructural features of tissue but also has the disadvantage to decrease the signal-to-noise ratio. Some studies suggested that SEM20,22 or BEM13,16,17 are efficient models only when the b-value range is sufficiently high, and supposedly less effective for an intermediate b-value range. However, the precise b-value range where these models were considered as accurate was not explicitly nor consensually defined, as it was found to varies between studies (for SEM19–23 and for BEM15–18) and was not systematically compared to the theoretical b-value limit of modelling accuracy given by the radius of convergence 26 of their mathematical formulations. On the other hand, using the TCE model for a well-defined b-value range (0–2500 s/mm2) is consistent with its theoretical limit of accuracy. 26 Then, it is arguable that an adequate and flexible model that can accurately fit the E(b) signal over a high b-value range can also fit it accurately over an intermediate b-value range. This is at least the case for the more general cumulant expansion, which is adequate over an intermediate b-value range and can be extended to higher orders to describe the signal for higher b-values. 42
In addition, using data acquired with high b-values raises the important problem of noise, which is known to have a detrimental effect on model comparisons. 43 Indeed, the higher the b-value, the weaker the signal-to-noise ratio and also the higher the risk that the positive noise level is fitted by overdetermined mathematical functions. 42 In parallel imaging, the noise in DW images is strictly positive definite 43 and is known to follow a Rician distribution for SENSE acquisitions. 43 In addition, non-homogeneity of the variance of noise occurs locally in dMRI data acquired with SENSE owing to the reconstruction process. Model accuracy comparisons in the high b-value range, especially when using high-degree-of-freedom models, may benefit from appropriate local denoising. 44 Combined model accuracy comparisons and denoising method comparisons will be used to model dMRI images acquired at high b-values and high directional resolution in subsequent work.
Conclusion
In summary, this work confirms that the low-degree-of-freedom TCE model is superior to the BEM for modeling and analyzing accurately the DW-dependent attenuation of the MRI signal E(b) in the brain tissue. Additionally, the TCE model was also found to be superior not only to the BEM but also to the SEM and the TEM, at least in WM. This reinforces the hypothesis that the phenomenological TCE model is sensitive to accurate microscopic properties of tissues, whereas the SEM, BEM, and TEM are not.1,4,18,26 The parameter σ deducted from the TCE acquisition, an indicator of tissue microscopic heterogeneity, might prove particularly interesting as a new biomarker of microstructural changes occurring in stroke, as it would signal a change in biomembrane disorder. 1
Author Contributions
Conceived and designed the experiments: RN and BH. Analyzed the data: RN. Wrote the first draft of the manuscript: RN. Contributed to the writing of the manuscript: IS. Agreed on manuscript results and conclusions: RN, IS, and BH. Jointly developed the structure and arguments for the paper: RN, IS, and BH. Made critical revisions and approved the final version: RN, IS, and BH. All the authors reviewed and approved the final manuscript.
Footnotes
Acknowledgements
We thank F. Aubry, P. Celsis, and J. Pariente for their role in improving this manuscript.
