In our simulations, we tested up to five distinct diffusion times to generate b-values. Figure 1 illustrates the effects of the number of diffusion times on the parameter estimation at various SNR levels. We draw attention to a number of features in the plots. First, as SNR is increased from 5 to 20 (rows 1–3), the variability in the estimated parameters (${\displaystyle D_{\beta }}$, ${\displaystyle \beta }$, ${\displaystyle K^{\ast }}$) decreases. Second, increasing the number of distinct diffusion times used for parameter estimation decreases estimation variability, with the most significant improvement when increasing from one to two diffusion times (rows 1–3). Third, sampling with two distinct diffusion times provides more robust parameter estimates than sampling twice as many b-values using one diffusion time (compare 2 and 2′ violin plots, rows 1–3). Fourth, the last row (row 4) highlights the improvement in the coefficient of variation (CV) for each parameter estimate with increasing SNR. This result again confirms that the most pronounced decline of CV occurs when increasing from one to two diffusion times, and parameter estimations can be performed more robustly using DW-MRI data with a relatively high SNR.

Figure 1

Download asset Open asset

In Figure 2, we provide simulation results evaluating the choice of the two distinct diffusion times (assuming ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1<{\mathrm{\Delta }}_2}}$) by measuring the goodness-of-fit of the model. The smaller of the two diffusion times is stated along the abscissa, and the difference, i.e., ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2-{\mathrm{\Delta }}_1}}$, is plotted along the ordinate. The quality of fitting was measured using the coefficient of determination (the larger the ${\displaystyle R^2}$ value, the better the goodness-of-fit of the model) for each combination of abscissa and ordinate values. The conclusion from this figure is that ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1}}$ should be small, whilst ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2}}$ should be as large as practically plausible. Note, DW-MRI echo time was not considered in this simulation, but as ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2}}$ increases, the echo time has to proportionally increase. Because of the inherent consequence of decreasing SNR with echo time, special attention should be paid to the level of SNR achievable with the use of a specific ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2}}$. Nonetheless, our findings suggest that when ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1=8\text{ ms}}}$, ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2}}$ can be set as small as ${\displaystyle 21\text{ms}}$ to achieve an ${\displaystyle R^2>0.90}$ with an SNR as low as 5. If ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1}}$ is increased past ${\displaystyle 15\text{ms}}$, then the separation between ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1}}$ and ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2}}$ has to increase as well, and such choices benefit from an increased SNR in the DW-MRI data. The Connectome 1.0 DW-MRI data was obtained using ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1=19\text{ ms}}}$ and ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2=49\text{ ms}}}$, leading to a separation of ${\displaystyle 30\text{ms}}$. For this data, it is expected that with SNR = 5 an ${\displaystyle R^2}$ around 0.9 is feasible, and by increasing SNR to 20, the ${\displaystyle R^2}$ can increase to a value above 0.99.

Figure 2

Download asset Open asset

In Figure 3, we present the scatter plots of simulated ${\displaystyle K}$ values versus fitted ${\displaystyle K}$ values using simulated data with different number of diffusion times at various SNR levels. Four cases are provided, including fitting simulated data generated with ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1=19\text{ ms}}}$ (row 1) or ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2=49\text{ ms}}}$ (row 2), fitting data with both diffusion times (row 3), and fitting data with three diffusion times (row 4). ${\displaystyle R^2}$ values for each case at each SNR level are provided. This result verifies that sub-diffusion based kurtosis estimation (blue dots) improves using multiple diffusion times. The improvement in ${\displaystyle R^2}$ is prominent when moving from fitting single diffusion time data to two diffusion times data, especially when the data is noisy (e.g., SNR = 5 and 10). The improvement gained by moving from two to three distinct diffusion times is marginal (less than 0.01 improvement in ${\displaystyle R^2}$ value at SNR = 5 and 10, and no improvements for SNR = 20 data). Moreover, our simulation findings highlight the deviation away from the true kurtosis ${\displaystyle K}$ by using the traditional DKI method (orange dots), especially with kurtosis values larger than 1. Overall, fitting sub-diffusion model (Equation 3) to data with two adequately separated diffusion times can lead to robust estimation of mean kurtosis, via (Equation 9).

Figure 3

Download asset Open asset

In Figure 4, we show the time-dependence effect of the DKI metrics (${\displaystyle D_{DKI}}$ and ${\displaystyle K_{DKI}}$) after fitting the standard DKI model to our simulated data. In this fitting, we consider b-values of 0, 1000, 1400, and 2500 s/mm², and vary the diffusion time, as in Jelescu et al., 2022. We depict the parameter estimates, ${\displaystyle D_{DKI}}$ and ${\displaystyle K_{DKI}}$, from simulated data with added noise (SNR = 20) in (A) and (B) for the diffusion time between 10 and 110 ms. In (C) and (D), using data with no added noise, we illustrate the long-term fitting results and trends in the parameter estimates. In both (A) and (C), as diffusion time increases, ${\displaystyle D_{DKI}}$ decreases as expected for an effective diffusion coefficient. The mean ${\displaystyle D_{DKI}}$ (averaged over 1000 simulations) agrees with the true diffusivity ${\displaystyle D^{\ast }}$. When it comes to the kurtosis ${\displaystyle K_{DKI}}$, in the noiseless data setting (D), we see ${\displaystyle K_{DKI}}$ is converging to the true kurtosis value ${\displaystyle K^{\ast }}$ at large diffusion time, whilst in the noisy data setting (B), this trend is not obvious within the experimental diffusion time window. More explanations of the observed time-dependence of diffusivity and kurtosis are provided in the Discussion.

Figure 4

Download asset Open asset

Next, we sought to identify the minimum number and combination of b-values to use for mean kurtosis estimation based on the sub-diffusion model (Equation 3). The simulation results were generated using the Connectome 1.0 DW-MRI data b-value setting, which has 16 b-values, 8 from ${\displaystyle {\displaystyle \mathrm{\Delta }=19\text{ms}}}$ and 8 from ${\displaystyle {\displaystyle \mathrm{\Delta }=49\text{ms}}}$.

In Figure 5, we computed ${\displaystyle R^2}$ values for all combinations of choices when the number of non-zero b-values is 2, 3, and 4. We then plotted the ${\displaystyle R^2}$ values sorted in ascending order for each considered SNR level. Notably, as the number of non-zero b-values increase, the achievable combinations increase as well (i.e., 120, 560, and 1820 for the three different non-zero b-value cases). The colours illustrate the proportion of b-values used in the fitting based on ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1}}$ and ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2}}$. For example, 1:0 means only ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1}}$ b-values were used, and 2:1 means two ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1}}$ and one ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2}}$ b-values were used to generate the result. The b-value combinations achieving highest ${\displaystyle R^2}$ values are displayed in the inset pictures and the b-values are provided in Table 1.

Figure 5

Download asset Open asset

Our simulation results in Table 1 suggest that increasing the number of non-zero b-values from two to four improves the quality of the parameter estimation, also achievable by fitting the sub-diffusion model (Equation 3) to higher SNR data. The gain is larger by using higher SNR data than by using more b-values. For example, going from two to four non-zero b-values with SNR = 5 data approximately doubles the ${\displaystyle R^2}$, whereas the ${\displaystyle R^2}$ almost triples when SNR = 5 data is substituted by SNR = 20 data. Additionally, the use of ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1}}$ or ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2}}$ alone is not preferred (also see Figure 5), and preference is towards first using ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1}}$ and then supplementing with b-values generated using ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2}}$. At all SNR levels when only two non-zero b-values are used, one b-value should be chosen based on the ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1}}$ set, and the other based on ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2}}$. Moving to three non-zero b-values requires the addition of another ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1}}$ b-value, and when four non-zero b-values are used then two from each diffusion time are required. If we consider an ${\displaystyle R^2=0.90}$ to be a reasonable goodness-of-fit for the sub-diffusion model, then at least three or four non-zero b-values are needed with an SNR = 20. If SNR = 10, then three non-zero b-values will not suffice. Interestingly, an ${\displaystyle R^2}$ of 0.85 can still be achieved when SNR = 20 and two optimally chosen non-zero b-values are used.

Table 1 summarises findings based on having different number of non-zero b-values with ${\displaystyle R^2}$ values deduced from the SNR = 5, 10, and 20 simulations. We have chosen to depict five b-value combinations producing the largest ${\displaystyle R^2}$ values for the two, three, and four non-zero b-value sampling cases. We found consistency in b-value combinations across SNR levels. Thus, we can conclude that a range of b-values can be used to achieve a large ${\displaystyle R^2}$ value, which is a positive finding, since parameter estimation does not stringently rely on b-value sampling. For example, using three non-zero b-values an ${\displaystyle R^2\ge 0.90}$ is achievable based on different b-value sampling. Importantly, two distinct diffusion times are required, and preference is towards including a smaller diffusion time b-value first. Hence, for three non-zero b-values we find two b-values based on ${\displaystyle {\displaystyle {\mathrm{\Delta }}_1}}$ and one based on ${\displaystyle {\displaystyle {\mathrm{\Delta }}_2}}$. This finding suggests one of the ${\displaystyle {\mathrm{\Delta }}_1}$ b-values can be chosen in the range 50 to 350 s/mm², and the other in the range 1500 to 4750 s/mm². Additionally, the ${\displaystyle {\mathrm{\Delta }}_2}$ b-value can also be chosen in a range, considering between 2300 to 4250 s/mm² based on the Connectome 1.0 b-value settings. The b-value sampling choices made should nonetheless be in view of the required ${\displaystyle R^2}$ value. Overall, sampling using two distinct diffusion times appears to provide quite a range of options for the DW-MRI data to be used to fit the sub-diffusion model parameters. The suggested optimal b-value sampling in the last row of Table 1, primarily chosen to minimise b-value size whilst maintaining a large ${\displaystyle R^2}$ value, may be of use for specific neuroimaging studies, which will be used to inform our discussion on feasibility for clinical practice.

The benchmark mean kurtosis estimation in the brain is established using the entire b-value range with all diffusion encoding directions available in the Connectome 1.0 DW-MRI data. For two subjects in different slices, Figure 6 provides the spatially resolved maps of mean kurtosis computed using the sub-diffusion method (i.e., ${\displaystyle K^{\ast }}$) with one or two diffusion times, and using the standard method (i.e., ${\displaystyle K_{DKI}}$) considering the two distinct diffusion times. First, we notice a degradation in the ${\displaystyle K_{DKI}}$ image with an increase in diffusion time. Second, the use of a single diffusion time with the sub-diffusion model leads to ${\displaystyle K^{\ast }}$ values which are larger than either the ${\displaystyle K_{DKI}}$ values or ${\displaystyle K^{\ast }}$ values generated using two diffusion times. Third, the quality of the mean kurtosis map appears to visually be best when two diffusion times are used to estimate ${\displaystyle K^{\ast }}$. Superior grey-white matter tissue contrast (TC) was found for the ${\displaystyle K^{\ast }}$ map (${\displaystyle TC=1.73}$), compared to the ${\displaystyle K_{DKI}}$ maps (${\displaystyle TC=0.80}$ for the ${\displaystyle {\displaystyle \mathrm{\Delta }=19\text{ ms}}}$ dataset and ${\displaystyle TC=1.01}$ for the ${\displaystyle {\displaystyle \mathrm{\Delta }=49\text{ ms}}}$ dataset).

Figure 6

Download asset Open asset

In Figure 7, an error map (measured by root-mean-square-error, RMSE) from Subject 5 slice 74 was presented for fitting the sub-diffusion model to the DW-MRI data with two diffusion times. Sample parameter fittings in both b-space (Equation 3) and q-space (Equation 4) were provided for four representative white and grey matter voxels.

Figure 7

Download asset Open asset

Quantitative findings for kurtosis are provided in Table 2. The analysis was performed for sub-cortical grey matter (scGM), cortical grey mater (cGM) and white matter (WM) brain regions. For specifics we refer the reader to the appropriate methods section. The table entries highlight the differences in mean kurtosis when computed using the two different approaches. The trend for the traditional single diffusion time approach is that an increase in ${\displaystyle {\displaystyle \mathrm{\Delta }}}$ results in a slight decrease in the mean ${\displaystyle K_{DKI}}$, and an increase in the coefficient of variation (CV) for any region. For example, the mean ${\displaystyle K_{DKI}}$ in the thalamus reduces from 0.65 to 0.58, whilst the CV increases from 30% to 39%. As much as 30% increase in CV is common for scGM and cGM regions, and around 10% for WM regions. The CV based on the ${\displaystyle K^{\ast }}$ value for each region is less than the CV for ${\displaystyle K_{DKI}}$ with either ${\displaystyle {\displaystyle \mathrm{\Delta }=19\text{ ms}}}$ or ${\displaystyle \mathrm{\Delta }=49\text{ ms}}$.

Figure 8 presents the distributions of the fitted parameters (${\displaystyle D}$ and ${\displaystyle K}$) in specific brain regions, based on the sub-diffusion model (Panel A) and the standard DKI model with ${\displaystyle {\displaystyle \mathrm{\Delta }=19\text{ ms}}}$ (Panel B). The distributions are coloured by the probability density. Yellow indicates high probability density, light blue indicates low probability density. In each subplot, the diffusivity is plotted along the abscissa axis and the kurtosis is along the ordinate axis. Results for the standard DKI model with ${\displaystyle {\displaystyle \mathrm{\Delta }=49\text{ ms}}}$ are qualitatively similar, so are not shown here. From Panel (B), we see an unknown nonlinear relationship between the DKI pair, ${\displaystyle D_{DKI}}$ and ${\displaystyle K_{DKI}}$, in all regions considered. By comparison, the sub-diffusion based ${\displaystyle K^{\ast }}$ and ${\displaystyle D^{\ast }}$ (Panel A) are less correlated with each other, indicating ${\displaystyle D^{\ast }}$ and ${\displaystyle K^{\ast }}$ carry distinct information, which will be very valuable for characterising tissue microstructure. Tables 3 and Table 4 present the ${\displaystyle D_{\beta }}$ and β values used to compute ${\displaystyle D^{\ast }}$ and ${\displaystyle K^{\ast }}$ in Table 2 and Figure 8.

Figure 8

Download asset Open asset

Results for reductions in diffusion encoding directions to achieve different levels of SNR with the purpose of shortening acquisition times will be benchmarked against the ${\displaystyle K^{\ast }}$ maps in Figure 6 and the ${\displaystyle K^{\ast }}$ values reported in Table 2.

Figure 9 presents the qualitative findings for two subjects generated using all, optimal and sub-optimal b-value samplings with SNR = 6 (3 non-collinear directions, 6 measurements), 10 (8 non-collinear directions, 16 measurements), and 20 (32 non-collinear directions, 64 measurements) DW-MRI data. The quality of the mean kurtosis map improves with increasing SNR, and also by optimising b-value sampling. Optimal sampling at SNR = 10 is qualitatively comparable to the SNR = 20 optimal sampling result, and to the benchmark sub-diffusion results in Figure 6.

Figure 9

Download asset Open asset

Quantitative verification of the qualitative observations are provided in Table 5. Significant differences in brain region-specific mean kurtosis values occur at the SNR = 6 level, which are not apparent when SNR = 10 or 20 data with optimal b-value sampling were used. The average errors are relative errors compared to the benchmark kurtosis values reported in Table 2. When using optimal b-values, average errors range from 22% to 43% at SNR = 6, from 13% to 43% at SNR = 10, and from 8% to 20% at SNR = 20, across brain regions. When using sub-optimal b-values, average errors range from 47% to 57% at SNR = 6, from 24% to 102% at SNR = 10, and from 27% to 72% at SNR = 20. Also, the brain region-specific CV for mean kurtosis was not found to change markedly when SNR = 10 or 20 data were used to compute ${\displaystyle K^{\ast }}$. The result of reducing the SNR to 6 leads to an approximate doubling of the CV for each brain region. These findings confirm that with optimal b-value sampling, the data quality can be reduced to around the SNR = 10 level, without a significant impact on the region-specific mean kurtosis estimates derived using the sub-diffusion model.

Figure 10 summarises the intraclass correlation coefficient (ICC) distribution results (μ for mean; σ for standard deviation) for the specific brain regions analysed. The two sets of ICC values were computed based on all DW-MRI (i.e., SNR = 20; subscript A) and the SNR = 10 optimal b-value sampling scheme (subscript O). As the value of μ approaches 1, the inter-subject variation in mean kurtosis is expected to greatly outweigh the intra-subject scan–rescan error. The value of μ should always be above 0.5, otherwise parameter estimation cannot be performed robustly and accurately, and values above 0.75 are generally accepted as good. The ${\displaystyle {\mu }_A}$ values for all regions were in the range 0.76 (thalamus) to 0.87 (caudate), and reduced to the range 0.57 (thalamus) to 0.80 (CC) when optimal sampling with SNR = 10 was used to estimate the ${\displaystyle K^{\ast }}$ value. Irrespective of which of the two DW-MRI data were used for ${\displaystyle K^{\ast }}$ estimation, the value of μ was greater than or equal to 0.70 in 20 out of 24 cases. The ${\displaystyle {\mu }_O}$ was less than 0.70 for only the thalamus, putamen, and pallidum. The loss in ICC by going to SNR = 10 data with optimal b-value sampling went hand-in-hand with an increase in ${\displaystyle \sigma }$, which is not unexpected, since the uncertainty associated with using less data should be measurable. Overall, ${\displaystyle {\mu }_A}$, ${\displaystyle {\mu }_O}$, and ${\displaystyle {\sigma }_A}$, ${\displaystyle {\sigma }_O}$, were fairly consistent across the brain regions, suggesting the DW-MRI data with SNR = 10 is sufficient for mean kurtosis estimation based on the sub-diffusion framework.

Figure 10

Download asset Open asset

Previous attempts have been made in optimising the DW-MRI acquisition protocol for mean kurtosis estimation based on the traditional, single diffusion time kurtosis model (Hansen et al., 2013; Hu et al., 2022; Poot et al., 2010; Hansen et al., 2016). Considerations have been made towards reducing the number of b-shells, directions per shell, and sub-sampling of DW-MRI for each direction in each shell. Our findings suggest that robust estimation of kurtosis cannot be achieved using the classical model for mean kurtosis, as highlighted previously (Yang et al., 2022; Ingo et al., 2014; Ingo et al., 2015; Barrick et al., 2020). A primary limitation of the traditional method is the use of the cumulant expansion resulting in sampling below a b-value of around ${\displaystyle 2500\text{s}/{\text{mm}}^2}$ (Jensen et al., 2005; Jensen and Helpern, 2010), and using the sub-diffusion model this limitation is removed (Yang et al., 2022). Our simulation findings and experiments using the Connectome 1.0 data confirm that mean kurtosis can be mapped robustly and rapidly using the sub-diffusion model applied to an optimised DW-MRI protocol. Optimisation of the data acquisition with the use of the sub-diffusion model has not been considered previously.

Mean kurtosis values can be generated based on having limited number of diffusion encoding directions (refer to Figure 9 and Table 5). Given that each direction for each b-shell takes a fixed amount of time, then a four b-shell acquisition with six directions per shell will take 25 times longer than a single diffusion encoding data acquisition (assuming a single b-value=0 data is collected). The total acquisition time for the diffusion MRI protocol for the Connectome 1.0 data was 55 min, including 50 b = 0 s/mm² scans plus seven b-values with 32 and nine with 64 diffusion encoding directions (Tian et al., 2022). This gives a total of 850 scans per dataset. As such, a single 3D image volume took ${\displaystyle 3.88\text{s}}$ to acquire. Conservatively allowing ${\displaystyle 4\text{s}}$ per scan, and considering SNR = 20 data (i.e., 64 directions) over four b-values and a single b-value = 0 scan, DW-MRI data for mean kurtosis estimation can be completed in ${\displaystyle 17\text{min}8\text{s}}$ (${\displaystyle R^2=0.96}$). At SNR = 10 (i.e., 32 directions), DW-MRI data with the same number of b-values can be acquired in ${\displaystyle 4\text{min}20\text{s}}$ (${\displaystyle R^2=0.91}$). If an ${\displaystyle R^2=0.85}$ (SNR = 10) is deemed adequate, then one less b-shell is required, saving an additional ${\displaystyle 64\text{s}}$. We should point out that even though 2-shell optimised protocols can achieve ${\displaystyle R^2=0.85}$ with SNR = 20, this is not equivalent in time to using 3-shells with SNR = 10 (also ${\displaystyle R^2=0.85}$). This is because 4×${\displaystyle 4\times }$ additional data are required to double the SNR (equivalent to acquiring an additional 4-shells). However, only one extra b-shell is required to convert 2-shell data to 3-shells with SNR = 10. Our expected DW-MRI data acquisition times are highly feasible clinically, where generally neuroimaging scans take around ${\displaystyle 15\text{min}}$ involving numerous different MRI contrasts and often a DTI acquisition.

An early study on estimating mean kurtosis demonstrated the mapping of a related metric in less than ${\displaystyle 1\text{min}}$ over the brain (Hansen et al., 2013). Clinical adoption of the protocol lacked, possibly since b-values are a function ${\displaystyle \delta }$, ${\displaystyle G}$, and ${\displaystyle \mathrm{\Delta }}$. Hence, different b-values can be obtained using different DW-MRI protocol settings, leading to differences in the mapping of mean kurtosis based on the data (we showed the ${\displaystyle {\displaystyle \mathrm{\Delta }}}$ effect in Figure 6 and Table 2). Our findings suggest this impediment is removed by sampling and fitting data with b-values across two distinct diffusion times. Nonetheless, we should consider what might be an acceptable DW-MRI data acquisition time.

A recent study on nasopharyngeal carcinoma investigated reducing the number of b-shell signals based on fixing diffusion encoding directions to the three Cartesian orientations (He et al., 2022). The 3-shell acquisitions took ${\displaystyle 3\text{min}2\text{s}}$ to acquire, whilst the 5-shell data required 5 min 13 s. They investigated as well the impact of using partial Fourier sampling, i.e., reducing the amount of data needed for image reconstruction by reducing k-space line acquisitions for each diffusion encoded image. Their benchmark used 5-shells (200, 400, 800, 1500, 2000 s/mm²), and found partial Fourier sampling with omission of the 1500 s/mm² b-shell produced acceptable results. With this acquisition the scan could be completed in 3 min 46 s, more than 2 times faster than the benchmark 8 min 31 s. Our proposed 3-shells acquisition with an ${\displaystyle R^2=0.85}$ (see SNR = 10 results in Table 1) executable under ${\displaystyle 4\text{min,}}$ is therefore inline with current expectations. Note, at the ${\displaystyle R^2=0.85}$ level the ICC for the different brain regions were in the range 0.60–0.69, and these were not formally reported in Figure 10. This level of reproducibility is still acceptable for routine use. We should additionally point out that we used the Subject 1 segmentation labels, after having registered each DW-MRI data to the Subject 1 first scan. This approach results in slight mismatch of the region-specific segmentations when carried across subjects, inherently resulting in an underestimation of ICC values.

Less than ${\displaystyle 4\text{min,}}$ DW-MRI data acquisitions can potentially replace existing data acquisitions used to obtained DTI metrics, since even the estimation of the apparent diffusion coefficient improves by using DW-MRI data relevant to DKI (Veraart et al., 2011b; Wu and Cheung, 2010). Additionally, it is increasingly clear that in certain applications the DKI analysis offers a more comprehensive approach for tissue microstructure analysis (Li et al., 2022b; Liu et al., 2022; Huang et al., 2022; Guo et al., 2022b; Goryawala et al., 2023; Guo et al., 2022a; Wang et al., 2022; Li et al., 2022a; Maralakunte et al., 2023; Hu et al., 2022; Zhou et al., 2023; Li et al., 2022c). As such, multiple b-shell, multiple diffusion encoding direction DW-MRI acquisitions should be used for the calculation of both DTI and DKI metrics.

To achieve ${\displaystyle R^2>0.92}$ for estimating kurtosis ${\displaystyle K^{\ast }}$, it is necessary to have four b-values, e.g., two relatively small b-values (350 s/mm² using ${\displaystyle {\displaystyle \mathrm{\Delta }=19\text{ ms}}}$, and 950 s/mm² using ${\displaystyle \mathrm{\Delta }=49\text{ ms}}$, both with ${\displaystyle G=68\text{mT}/\text{m}}$) and two larger b-values of around 1500 s/mm² and 4250 s/mm² (using ${\displaystyle {\displaystyle \mathrm{\Delta }=19\text{ ms}}}$ and ${\displaystyle \mathrm{\Delta }=49\text{ ms}}$, respectively, both with ${\displaystyle G=142\text{mT}/\text{m}}$) (see bottom row in Table 1). If two or three non-zero b-values are considered sufficient (with ${\displaystyle R^2=0.85}$ or ${\displaystyle 0.90}$), then the larger ${\displaystyle \mathrm{\Delta }}$ needs to be used to set the largest b-value to be 2300 s/mm², and the other(s) should be set using the smaller ${\displaystyle {\displaystyle \mathrm{\Delta }}}$. For the two non-zero b-values case, the b-value from the smaller ${\displaystyle \mathrm{\Delta }}$ should be around 800 s/mm². For the three non-zero b-values case, the b-values from the smaller ${\displaystyle \mathrm{\Delta }}$ would then need to be 350 s/mm² and 1500 s/mm². Interestingly, b-value = 800 s/mm² lies around the middle of the 350 s/mm² to 1500 s/mm² range. The additional gain to ${\displaystyle R^2=0.92}$ can be achieved by splitting the b-value with the larger ${\displaystyle \mathrm{\Delta }}$ into two, again with ${\displaystyle 2300\text{s}/{\text{mm}}^2}$ near the middle of the two new b-values set. In addition, the separation between ${\displaystyle {\mathrm{\Delta }}_1}$ and ${\displaystyle {\mathrm{\Delta }}_2}$ needs to be as large as plausible, as can be deduced from the simulation result in Figure 2, but attention should be paid to signal-to-noise ratio decreases with increased echo times (He et al., 2022).

A recent study on optimising quasi-diffusion imaging (QDI) considered b-values up to 5000 s/mm² (Spilling et al., 2022). Whilst QDI is a non-Gaussian approach, it is different to the sub-diffusion model, but still uses the Mittag-Leffler function and involves the same number of model parameters. The DW-MRI data used in their study was acquired with a single diffusion time. Nevertheless, particular points are worth noting. Their primary finding was parameter dependence on the maximum b-value used to create the DW-MRI data. They also showed the accuracy, precision, and reliability of parameter estimation are improved with increased number of b-shells. They suggested a maximum b-value of 3960 s/mm² for the 4-shell parameter estimation. Our results do not suggest a dependence of the parameter estimates on maximum b-value (note, if a maximum b-value dependence is present, benchmark versus optimal region specific results in Tables 2 and 5 should show some systematic difference; and we used the real part of the DW-MRI data and not the magnitude as commonly used). These findings potentially confirm that ${\displaystyle \mathrm{\Delta }}$ separation is an important component of obtaining a quality parameter fit from which mean kurtosis is deduced (Figure 2).

Commonly available human clinical MRI scanners are capable of ${\displaystyle 40\text{mT/m}}$ gradient amplitudes. Recently, the increased need to deduce tissue microstructure metrics from DW-MRI measurements has led to hardware developments resulting in 80 mT/m gradient strength MRI scanners. These were initially sought by research centres. The Connectome 1.0 scanners achieve 300 mT/m gradient amplitudes (Tian et al., 2022), which in turn allow large b-values within reasonable echo times. The Connectome 2.0 scanner is planned to achieve 500 mT/m gradient amplitudes (Huang et al., 2021), providing data for further exploration of the ${\displaystyle {\displaystyle (q,\mathrm{\Delta })}}$ space by providing a mechanism for increasing our knowledge of the relationships between the micro-, meso-, and macro-scales. Whilst there are three Connectome 1.0 scanners available, only one Connectome 2.0 scanner is planned for production. Hence, robust and fast methods utilising existing 80 mT/m gradient systems are needed, and the Connectome scanners provide opportunities for testing and validating methods.

The b-value is an intricate combination of ${\displaystyle \delta }$, ${\displaystyle G}$, and ${\displaystyle {\displaystyle \mathrm{\Delta }}}$. An increase in any of these three parameters increases the b-value (note, ${\displaystyle \delta }$ and ${\displaystyle G}$ increase b-value quadratically, and ${\displaystyle \mathrm{\Delta }}$ linearly). An increase in ${\displaystyle \delta }$ can likely require an increase in ${\displaystyle \mathrm{\Delta }}$, the consequence of which is a reduction in signal-to-noise ratio as the echo time has to be adjusted proportionally. Partial Fourier sampling methods aim to counteract the need to increase the echo time by sub-sampling the DW-MRI data used to generate an image for each diffusion encoding direction (Zong et al., 2021; Heidemann et al., 2010). The ideal scenario, therefore, is to increase ${\displaystyle G}$, as for the Connectome 1.0 and 2.0 scanners.

Based on our suggested b-value settings in Table 1, a maximum b-value of around 4250 s/mm² is required for robust mean kurtosis estimation (assume ${\displaystyle R^2>0.90}$ is adequate and achieved via four non-zero b-values and two distinct Δs, and we should highlight that b-value = 1500 s/mm² (${\displaystyle {\displaystyle \mathrm{\Delta }=19\text{ ms}}}$) and b-value = 4250 s/mm² (${\displaystyle {\displaystyle \mathrm{\Delta }=49\text{ ms}}}$) were achieved using ${\displaystyle G=142\text{mT/m}}$). Considering 80 m/Tm gradient sets are the new standard for MRI scanners, an adjustment to ${\displaystyle \delta }$ to compensate for ${\displaystyle G}$ is needed (recall, ${\displaystyle q=\gamma \delta G}$). Hence, for a constant ${\displaystyle q}$, ${\displaystyle G}$ can be reduced at the consequence of proportionally increasing ${\displaystyle \delta }$. This then allows MRI systems with lower gradient amplitudes to generate the relevant b-values. For example, changing the ${\displaystyle \delta }$ from 8 ms to 16 ms would result in halving of ${\displaystyle G}$ (using a maximum ${\displaystyle G}$ of 142 m/TM from the Connectome 1.0 can achieve the optimal b-values; hence halving would require around 71 m/TM gradient pulse amplitudes). Increasing ${\displaystyle \mathrm{\Delta }}$ above 49 ms to create larger b-value data is unlikely to be a viable solution due to longer echo times leading to a loss in SNR. SNR increases afforded by moving from 3 T to 7 T MRI are most likely counteracted by an almost proportional decrease in ${\displaystyle T_2}$ times (Pohmann et al., 2016), in addition to 7 T MRI bringing new challenges in terms of increased transmit and static field inhomogeneities leading to signal inconsistencies across an image (Kraff and Quick, 2017).

Following Table 2 (last column), white matter regions showed high kurtosis (0.87 ± 0.22), consistent with a structured heterogeneous environment comprising parallel neuronal fibres, as shown in Maiter et al., 2021. Cortical grey matter showed low kurtosis (0.40 ± 0.16). Subcortical grey matter regions showed intermediate kurtosis (0.60 ± 0.21). In particular, caudate and putamen showed similar kurtosis to grey matter, whilst thalamus and pallidum showed similar kurtosis properties to white matter. Histological staining results of the sub-cortical nuclei (Maiter et al., 2021) showed the subcortical grey matter was permeated by small white matter bundles, which could account for the similar kurtosis in thalamus and pallidum to white matter. These results confirmed that our proposed sub-diffusion based mean kurtosis ${\displaystyle K^{\ast }}$ is consistent with published histology of normal human brain.

The time-dependence of diffusivity and kurtosis has attracted much interest in the field of tissue microstructure imaging. Although our motivation here is not to map time-dependent diffusion, we can nonetheless point out that the assumption of a sub-diffusion model provides an explanation of the observed time-dependence of diffusivity and kurtosis. From (Equation 4), the diffusivity that arises from the sub-diffusion model is of the form ${\displaystyle D_{SUB}=D_{\beta }{\overline{\mathrm{\Delta }}}^{\beta -1}}$, where ${\displaystyle \overline{\mathrm{\Delta }}}$ is the effective diffusion time, ${\displaystyle D_{SUB}}$ has the standard units for a diffusion coefficient ${\displaystyle {\text{mm}}^2\text{/s}}$, and ${\displaystyle D_{\beta }}$ is the anomalous diffusion coefficient (with units ${\displaystyle {\text{mm}}^2{\text{/s}}^{\beta }}$) associated with sub-diffusion. In the sub-diffusion framework (Equation 1), ${\displaystyle D_{\beta }}$ and ${\displaystyle \beta }$ are assumed to be constant, and hence ${\displaystyle D_{SUB}}$ exhibits a fractional power-law dependence on diffusion time. Then, following (Equation 8), ${\displaystyle D^{\ast }}$ is obtained simply by scaling ${\displaystyle D_{SUB}}$ by a constant ${\displaystyle 1/\mathrm{\Gamma }(1+\beta )}$, and hence also follows a fractional power-law dependence on diffusion time. This time-dependence effect of diffusivity was illustrated in our simulation results, Figure 4(A) and (C).

When it comes to kurtosis, the literature on the time-dependence is mixed. Some work showed kurtosis to be increasing with diffusion time in both white and grey matter (Aggarwal et al., 2020), and in grey matter (Ianus et al., 2021), whilst others showed kurtosis to be decreasing with diffusion time in grey matter (Lee et al., 2020; Olesen et al., 2022; Jelescu et al., 2022). In this study, we provide an explanation of these mixed results. We construct a simulation of diffusion MRI signal data based on the sub-diffusion model (Equation 3) augmented with random Gaussian noise. Then we fit the conventional DKI model to the synthetic data. As shown in Figure 4(D), when there is no noise, ${\displaystyle K_{DKI}}$ increases with diffusion time in white matter, whilst decreasing with diffusion time in grey matter. When there is added noise, as shown in Figure 4(B), the time-dependency of kurtosis within the timescale of a usual MR experiment is not clear. This goes some way to explaining why the results in the literature on the time-dependence of kurtosis are quite mixed.

Furthermore, we summarise the benefits of using the sub-diffusion based mean kurtosis measurement ${\displaystyle K^{\ast }}$. First, as shown in (Equation 9), sub-diffusion based mean kurtosis ${\displaystyle K^{\ast }}$ is not time dependent, and hence has the potential to become a tissue-specific imaging biomarker. Second, the fitting of the sub-diffusion model is straightforward, fast and robust, from which the kurtosis ${\displaystyle K^{\ast }}$ is simply computed as a function of the sub-diffusion model parameter β. Third, the kurtosis ${\displaystyle K^{\ast }}$ is not subject to any restriction on the maximum b-value, as in standard DKI. Hence, its value truly reflects the information contained in the full range of b-values.

The direct link between the sub-diffusion model parameter β and mean kurtosis is well established (Yang et al., 2022; Ingo et al., 2014; Ingo et al., 2015). An important aspect to consider is whether mean β used to compute the mean kurtosis is alone sufficient for clinical decision making. Whilst benefits of using kurtosis metrics over other DW-MRI data-derived metrics in certain applications are clear, the adequacy of mean kurtosis over axial and radial kurtosis is less apparent. Most studies perform the mapping of mean kurtosis, probably because the DW-MRI data can be acquired in practically feasible times. Nonetheless, we can point to a few recent examples where the measurement of directional kurtosis has clear advantages. A study on mapping tumour response to radiotherapy treatment found axial kurtosis to provide the best sensitivity to treatment response (Goryawala et al., 2023). In a different study, a correlation was found between glomerular filtration rate and axial kurtosis is assessing renal function and interstitial fibrosis (Li et al., 2022a). Uniplor depression subjects have been shown to have brain region specific increases in mean and radial kurtosis, whilst for bipolar depression subjects axial kurtosis decreased in specific brain regions and decreases in radial kurtosis were found in other regions (Maralakunte et al., 2023). This selection of studies highlights future opportunities for extending the methods to additionally map axial and radial kurtosis.

Notably, estimates for axial and radial kurtosis require directionality of kurtosis to be resolved, resulting in DW-MRI sampling over a large number of diffusion encoding directions within each b-shell (Jensen and Helpern, 2010; Poot et al., 2010). As such, extension to directional kurtosis requires a larger DW-MRI dataset acquired using an increased number of diffusion encoding directions. The number of b-shells and directions therein necessary for robust and accurate mapping of directional kurtosis based on the sub-diffusion model is an open question.

There are three primary ways of determining mean kurtosis. These include the powder averaging over diffusion encoding directions in each shell, and then fitting the model, as in our approach. A different approach is to ensure each b-shell in the DW-MRI data contains the same diffusion encoding directions, and then kurtosis can be estimated for each diffusion encoding direction, after which the average over directions is used to state mean kurtosis. Lastly, the rank-4 kurtosis tensor is estimated from the DW-MRI data, from which mean kurtosis is computed directly. The latter two approaches are potential candidates for extending to axial and radial kurtosis mapping. Note, in DTI a rank-2 diffusion tensor with six unique tensor entries is needed to be estimated, whilst in DKI in addition to the rank-2 diffusion tensor, the kurtosis tensor is rank-4 with 15 unknowns, resulting in 21 unknowns altogether (Hansen et al., 2016). As such, DKI analysis for directional kurtosis requires much greater number of diffusion encoding directions to be sampled than DTI. This automatically means that at least 22 DW-MRI data (including b-value = 0) with different diffusion encoding properties have to be acquired (Jensen and Helpern, 2010). The traditional approach has been to set five distinct b-values with 30 diffusion encoding directions within each b-shell (Poot et al., 2010). Hence, to obtain the entries of the rank-4 kurtosis tensor, much more DW-MRI data is needed in comparison to what is proposed for mean kurtosis estimation in this study. Estimation of the tensor entries from this much data is prone to noise, motion and image artifacts in general (Tabesh et al., 2011), posing challenges on top of long DW-MRI data acquisition times.

A kurtosis tensor framework based on the sub-diffusion model where separate diffusion encoding directions are used to fit a direction specific β is potentially an interesting line of investigation for the future, since it can be used to establish a rank-2 β tensor with only six unknowns, requiring at least six distinct diffusion encoding directions. This type of approach can reduce the amount of DW-MRI data to be acquired, and potentially serve as a viable way forward for the combined estimation of mean, axial, and radial kurtosis.

Although our study has been focusing on mean kurtosis imaging in the human brain, it is clear that DKI has wide application outside of the brain (Li et al., 2022b; Liu et al., 2022; Huang et al., 2022; Guo et al., 2022b; Li et al., 2022a; Zhou et al., 2023). Without having conducted experiments elsewhere, we cannot provide specific guidelines for mean kurtosis imaging in the breast, kidney, liver, and other human body regions. We can, however, point the reader in a specific direction.

The classical mono-exponential model can be recovered from the sub-diffusion equation by setting ${\displaystyle \beta =1}$. For this case, the product between the b-value and fitted diffusivity has been reported to be insightful for b-value sampling (Yablonskiy and Sukstanskii, 2010), in accordance with a theoretical perspective (Istratov and Vyvenko, 1999). It was suggested the product should approximately span the (0, 1) range. Considering our case based on the sub-diffusion equation, we can investigate the size of ${\displaystyle b{D}_{SUB}}$ by analysing the four non-zero b-value optimal sampling regime (${\displaystyle {\mathrm{\Delta }}_1}$: 350 s/mm² and 1500 s/mm²; ${\displaystyle {\mathrm{\Delta }}_2}$: 950 s/mm² and 4250 s/mm² from Table 1). Considering scGM, cGM, and WM brain regions alone, the rounded and dimensionless ${\displaystyle b{D}_{SUB}}$ values are (0.09, 0.38, 0.20, 0.88), (0.13, 0.55, 0.30, 1.34), and (0.05, 0.23, 0.11, 0.48), respectively, and note that in each case the first two effective sampling values are based on ${\displaystyle {\mathrm{\Delta }}_1}$, and the other two are derived using ${\displaystyle {\mathrm{\Delta }}_2}$. Interestingly, the log-linear sampling proposed in Istratov and Vyvenko, 1999 is closely mimicked by the effective sampling regime (scGM: −2.42, −1.62, −0.97, −0.13; cGM: −2.06, −1.20, −0.60, 0.29; WM: −2.94, −2.23, −1.48, −0.73; by sorting and taking the natural logarithm). This analysis also informs on why it may be difficult to obtain a generally large ${\displaystyle R^2}$ across the entire brain, since β and ${\displaystyle D_{SUB}}$ are brain region specific and the most optimal sampling strategy should be β and ${\displaystyle D_{SUB}}$ specific. Whilst region specific sampling may provide further gains in the ${\displaystyle R^2}$ value, and improve ICC values for specific brain regions, such data would take a long time to acquire and require extensive post-processing and in-depth analyses.

The DW-MRI dataset provided by Tian et al., 2022 was used in this study. The publicly available data were collected using the Connectome 1.0 scanner for 26 healthy participants. The first seven subjects had a scan–rescan available. We evaluated qualitatively the seven datasets, and chose the six which had consistent diffusion encoding directions. Subject 2 had 60 instead of 64 diffusion encoding directions, and hence, was omitted from this study. The ${\displaystyle 2\times 2\times 2{\text{mm}}^3}$ resolution data were obtained using two diffusion times ${\displaystyle ({\displaystyle \mathrm{\Delta }=19,49\text{ ms})}}$ with a pulse duration of ${\displaystyle \delta =8\text{ms}}$ and ${\displaystyle G=31,68,105,142,179,216,253,290\text{mT/m}}$, respectively, generating b-values = 50, 350, 800, 1500, 2400, 3450, 4750, 6000 s/mm² for ${\displaystyle \mathrm{\Delta }=19\text{ ms}}$, and b-values = 200, 950, 2300, 4250, 6750, 9850, 13500, 17800 s/mm² for ${\displaystyle {\displaystyle \mathrm{\Delta }=49\text{ ms}}}$, according to b-value ${\displaystyle =(\gamma \delta G{)}^2(\mathrm{\Delta }-\delta /3)}$. 32 diffusion encoding directions were uniformly distributed on a sphere for ${\displaystyle b<2400\text{s}/{\text{mm}}^2}$ and 64 uniform directions for ${\displaystyle b\ge 2400\text{s}/{\text{mm}}^2}$.

The FreeSurfer’s segmentation labels as part of the dataset were used for brain-region-specific analyses. Tian et al., 2022 provided the white matter averaged group SNR (23.10 ± 2.46), computed from 50 interspersed b-value ${\displaystyle =0\text{s}/{\text{mm}}^2}$ images for each subject. Both magnitude and the real part of the DW-MRI were provided. Based on an in-depth analysis, the use of the real part of the DW-MRI data was recommended, wherein physiological noise, by nature, follows a Gaussian distribution (Gudbjartsson and Patz, 1995).

We performed SNR reduction of the Connectome 1.0 DW-MRI data by downsampling of diffusion directions in each b-shell. The method of multiple subsets from multiple sets (P-D-MM) subsampling algorithm provided in DMRITool (Cheng et al., 2018) was applied to the b-vectors provided with the Connectome 1.0 DW-MRI data. Note, the b-shells contained 32 directions if ${\displaystyle b<2400\text{s}/{\text{mm}}^2}$, and 64 directions if ${\displaystyle b\ge 2400\text{s}/{\text{mm}}^2}$. We consider SNR = 20 to be the full dataset. The SNR = 10 data was constructed by downsampling to eight non-collinear diffusion encoding directions in each b-shell, and three were required for the SNR = 6 data. In the downsampled data, each diffusion encoding direction was coupled with the direction of opposite polarity (i.e., SNR = 10 had sixteen measurements for each b-value, and SNR = 6 had six).

For all of the simulation results, the coefficient of determination, referred to as ${\displaystyle R^2}$, was used to assess how well the mean kurtosis values were able to be fitted. This value was calculated using a standard ${\displaystyle R^2}$ formula

where ${\displaystyle {\overline{K}}_{\text{Simulated}}}$ is the mean simulated ${\displaystyle K}$ value.

For human data, we computed the region specific mean and standard deviation for each subject, and reported the weighted mean and pooled standard deviation along with the coefficient of variation (CV), defined as the ratio of the standard deviation to the mean. The weights were the number of voxels in the associated regions in each subject. Following Barrick et al., 2020, the tissue contrast (TC) is computed as ${\displaystyle TC=\left|{\mu }_{WM}-{\mu }_{GM}\right|/\sqrt{{\sigma }_{WM}^2+{\sigma }_{GM}^2}}$, where ${\displaystyle {\mu }_{WM}}$ and ${\displaystyle {\mu }_{GM}}$ are the mean parameter values in white and grey matters; and ${\displaystyle {\sigma }_{WM}}$ and ${\displaystyle {\sigma }_{GM}}$ are the standard deviations of parameter values. Higher TC values indicate greater tissue contrast.

Human data were analysed voxelwise, and also based on regions-of-interest. We considered three categories of brain regions, namely sub-cortical grey matter (scGM), cortical grey matter (cGM), and white matter (WM). The scGM region constituted the thalamus (FreeSurfer labels 10 and 49 for left and right hemisphere), caudate (11, 50), putamen (12, 51), and pallidum (13, 52). The cGM region was all regions (1000–2999) and separately analysed the fusiform (1007, 2007) and lingual (1013, 2013) brain regions, whilst the WM had white matter fibre regions from the cerebral (2, 41), cerebellum (7, 46), and corpus callosum (CC; 251–255) areas. The average number of voxels in each region were 3986 (scGM), 53326 (cGM), 52121 (WM), 1634 (thalamus), 831 (caudate), 1079 (putamen), 443 (pallidum), 2089 (fusiform), 1422 (lingual), 48770 (cerebral WM), 2904 (cerebellum WM), and 447 (CC). For each brain region, a t-test was performed to test for significant differences in mean kurtosis population means under the optimal and sub-optimal b-value sampling regimes presented in Table 5 compared to the benchmark parameter values presented in Table 2.

The utility of DKI for inferring information on tissue microstructure was described decades ago. Continued investigations in the DW-MRI field have led to studies clearly describing the importance of mean kurtosis mapping to clinical diagnosis, treatment planning and monitoring across a vast range of diseases and disorders. Our research on robust, fast, and accurate mapping of mean kurtosis using the sub-diffusion mathematical framework promises new opportunities for this field by providing a clinically useful, and routinely applicable mechanism for mapping mean kurtosis in the brain. Future studies may derive value from our suggestions and apply methods outside the brain for broader clinical utilisation.

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

Qianqian Yang and Viktor Vegh acknowledge the financial support from the Australian Research Council (ARC) Discovery Project scheme (DP190101889) for funding a project on mathematical model development and MRI-based investigations into tissue microstructure in the human brain. Qianqian Yang also acknowledges the support from the ARC Discovery Early Career Research Award (DE150101842) for funding a project on new mathematical models for capturing heterogeneity in human brain tissue. Authors also thank the members of the Anomalous Relaxation and Diffusion Study (ARDS) group for many interesting discussions involving diffusion MRI.

This study uses a publicly available MGH Connectome Diffusion Microstructure Dataset. The institutional review board of Mass General Brigham approved the study protocol and ensured that all relevant ethical regulations were complied with. Written informed consent was obtained from all participants. More details about the MGH Connectome Diffusion Microstructure Dataset are available at https://doi.org/10.1038/s41597-021-01092-6.

1. Preprint posted: June 28, 2023
2. Sent for peer review: June 29, 2023
3. Reviewed Preprint version 1: November 2, 2023
4. Reviewed Preprint version 2: August 14, 2024
5. Version of Record published: October 7, 2024

You can cite all versions using the DOI https://doi.org/10.7554/eLife.90465. This DOI represents all versions, and will always resolve to the latest one.

© 2023, Farquhar et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.