Robust, fast and accurate mapping of diffusional mean kurtosis
eLife assessment
This paper proposes a valuable new method for the assessment of the mean kurtosis for diffusional kurtosis imaging by utilizing a recently introduced subdiffusion model. The evidence supporting the claims that this technique is robust and accurate in brain imaging is solid; however, there is a need to include a summary of the clear limitations.
https://doi.org/10.7554/eLife.90465.3.sa0Valuable: Findings that have theoretical or practical implications for a subfield
 Landmark
 Fundamental
 Important
 Valuable
 Useful
Solid: Methods, data and analyses broadly support the claims with only minor weaknesses
 Exceptional
 Compelling
 Convincing
 Solid
 Incomplete
 Inadequate
During the peerreview process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments
Abstract
Diffusional kurtosis imaging (DKI) is a methodology for measuring the extent of nonGaussian diffusion in biological tissue, which has shown great promise in clinical diagnosis, treatment planning, and monitoring of many neurological diseases and disorders. However, robust, fast, and accurate estimation of kurtosis from clinically feasible data acquisitions remains a challenge. In this study, we first outline a new accurate approach of estimating mean kurtosis via the subdiffusion mathematical framework. Crucially, this extension of the conventional DKI overcomes the limitation on the maximum bvalue of the latter. Kurtosis and diffusivity can now be simply computed as functions of the subdiffusion model parameters. Second, we propose a new fast and robust fitting procedure to estimate the subdiffusion model parameters using two diffusion times without increasing acquisition time as for the conventional DKI. Third, our subdiffusionbased kurtosis mapping method is evaluated using both simulations and the Connectome 1.0 human brain data. Exquisite tissue contrast is achieved even when the diffusion encoded data is collected in only minutes. In summary, our findings suggest robust, fast, and accurate estimation of mean kurtosis can be realised within a clinically feasible diffusionweighted magnetic resonance imaging data acquisition time.
Introduction
Diffusionweighted magnetic resonance imaging (DWMRI) over a period of more than 30 years has become synonymous with tissue microstructure imaging. Measures of how water diffuses in heterogeneous tissues allow indirect interpretation of changes in tissue microstructure (Le Bihan and JohansenBerg, 2012). DWMRI has predominantly been applied in the brain, where properties of white matter connections between brain regions are often studied (Lebel et al., 2019), in addition to mapping tissue microstructural properties (Tournier, 2019). Applications outside of the brain have clinical importance as well, and interest is growing rapidly.
Generally, DWMRI involves the setting of diffusion weightings and direction over which diffusion is measured. Whilst diffusion tensor imaging (DTI) can be performed using a single diffusion weighting, a socalled bshell, and at least six diffusion encoding directions (Le Bihan et al., 2001), other models tend to require multiple bshells each having multiple diffusion encoding directions. Diffusional kurtosis imaging (DKI) is a primary example of a multiple bshell, multiple diffusion encoding direction method (Jensen et al., 2005). DKI is considered as an extension of DTI (Jensen and Helpern, 2010; Hansen et al., 2013; Veraart et al., 2011b), where the diffusion process is assumed to deviate away from standard Brownian motion, and the extent of such deviation is measured via the kurtosis metric. Essentially, the increased sampling achieved via DKI data acquisitions allows more complex models to be applied to data (Van Essen et al., 2013; Shafto et al., 2014), in turn resulting in metrics of increased utility for clinical decision making.
Recent clinical benefits of using kurtosis metrics over other DWMRIderived measures have been demonstrated for grading hepatocellular carcinoma (Li et al., 2022b), prognosing chronic kidney disease (Liu et al., 2022), differentiating parotid gland tumours (Huang et al., 2022), measuring response to radiotherapy treatment in bone tumour (Guo et al., 2022b) and glioblastoma (Goryawala et al., 2023), identifying tissue abnormalities in temporal lobe epilepsy patients with sleep disorders (Guo et al., 2022a) and brain microstructural changes in mild traumatic brain injury (Wang et al., 2022), monitoring of renal function and interstitial fibrosis (Li et al., 2022a), detecting the invasiveness of bladder cancer into muscle (Li et al., 2022d), aiding management of patients with depression (Maralakunte et al., 2023), delineating acute infarcts with prognostic value (Hu et al., 2022), predicting breast cancer metastasis (Zhou et al., 2023), diagnosing Parkinson’s disease (Li et al., 2022c), amongst others reported prior and not listed here.
The routine use of DKI in the clinic has nonetheless lagged due the inability to robustly estimate the kurtosis metric (Veraart et al., 2011a; Tabesh et al., 2011; Kuder et al., 2012; Henriques et al., 2021). A known requirement for estimating kurtosis in DKI is to restrict the maximum bvalue to $2000{\displaystyle 3000{\displaystyle \text{s}/{\text{mm}}^{2}}}$ for brain studies (Jensen et al., 2005; Jensen and Helpern, 2010; Poot et al., 2010), with the optimal maximum bvalue found to be dependent on tissue type (Poot et al., 2010). This suggests that the traditional kurtosis model is less accurate at representing the diffusion signal at large bvalues. Moreover, multiple bshell, multiple direction highquality DWMRI data can take many minutes to acquire, which poses challenges for clinical imaging protocols involving a multitude of MRI contrasts already taking tens of minutes to execute. Reduction of DKI data acquisition times through parallel imaging, optimisation of bshells and directions have been investigated (Zong et al., 2021; Heidemann et al., 2010; Zelinski et al., 2008), and DWMRI data necessary for DKI analysis has been shown to supersede the data required for DTI (Veraart et al., 2011b). Therefore, an optimised DKI protocol can potentially replace clinical DTI data acquisitions without adversely affecting the estimation of DTI metrics.
For DKI to become a routine clinical tool, DWMRI data acquisition needs to be fast and provides a robust estimation of kurtosis. The ideal protocol should have a minimum number of bshells and diffusion encoding directions in each bshell. The powder averaging over diffusion directions improves the signaltonoise ratio (SNR) of the DWMRI data used for parameter estimation. Whilst this approach loses out on directionality of the kurtosis, it nonetheless provides a robust method of estimating mean kurtosis (Henriques et al., 2021), a metric of significant clinical value.
Instead of attempting to improve an existing modelbased approach for kurtosis estimation, as has been considered by many others, we considered the problem from a different perspective. In view of the recent generalisation of the various models applicable to DWMRI data (Yang et al., 2022), the subdiffusion framework provides new, unexplored opportunities, for fast and robust kurtosis mapping. We report on our investigation into the utility of the subdiffusion model for practically useful mapping of mean kurtosis.
Results
Simulation studies were conducted to establish the requirements on the number of diffusion times and the separation between them for accurate estimation of mean kurtosis based on the subdiffusion model augmented with random Gaussian noise, following (Equation 10). Testing and validation were performed using the human Connectome 1.0 brain dataset (Tian et al., 2022). The 2 × 2 × 2 mm^{3} resolution data were obtained using two diffusion times $({\displaystyle \mathrm{\Delta}=19,49\phantom{\rule{thinmathspace}{0ex}}\text{ms})}$ with a pulse duration of $\displaystyle \delta =8\phantom{\rule{thinmathspace}{0ex}}\text{ms}$ and G = 31, 68, 105, 142, 179, 216, 253, and 290 mT/m, respectively, generating bvalues = 50, 350, 800, 1500, 2400, 3450, 4750, and 6000 s/mm^{2}, and bvalues = 200, 950, 2300, 4250, 6750, 9850, 13500, and 17800 s/mm^{2}, according to $b\text{value}=(\gamma \delta G{)}^{2}(\mathrm{\Delta}\delta /3)$. Up to 64 diffusion encoding directions per bshell were set. The traditional method for mean kurtosis estimation was implemented (producing $K}_{DKI$), which is limited to the use of DWMRI generated using a single diffusion time (Jensen et al., 2005; Jensen and Helpern, 2010; Veraart et al., 2011a; Poot et al., 2010), alongside our implementation based on the subdiffusion model (Equation 3), wherein mean kurtosis $K}^{\ast$ is computed as a function of the subdiffusion model parameter β (refer to (Equation 9)) using either a single or multiple diffusion times.
Multiple diffusion times for robust and accurate mean kurtosis estimation
In our simulations, we tested up to five distinct diffusion times to generate bvalues. 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 ($D}_{\beta$, $\beta$, $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 bvalues 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 DWMRI data with a relatively high SNR.
In Figure 2, we provide simulation results evaluating the choice of the two distinct diffusion times (assuming $\displaystyle {\mathrm{\Delta}}_{1}<{\mathrm{\Delta}}_{2}$) by measuring the goodnessoffit of the model. The smaller of the two diffusion times is stated along the abscissa, and the difference, i.e., $\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 $R}^{2$ value, the better the goodnessoffit of the model) for each combination of abscissa and ordinate values. The conclusion from this figure is that $\displaystyle {\mathrm{\Delta}}_{1}$ should be small, whilst $\displaystyle {\mathrm{\Delta}}_{2}$ should be as large as practically plausible. Note, DWMRI echo time was not considered in this simulation, but as $\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 {\mathrm{\Delta}}_{2}$. Nonetheless, our findings suggest that when $\displaystyle {\mathrm{\Delta}}_{1}=8\text{ms}$, $\displaystyle {\mathrm{\Delta}}_{2}$ can be set as small as $21\text{ms}$ to achieve an ${R}^{2}>0.90$ with an SNR as low as 5. If $\displaystyle {\mathrm{\Delta}}_{1}$ is increased past $15\text{ms}$, then the separation between $\displaystyle {\mathrm{\Delta}}_{1}$ and $\displaystyle {\mathrm{\Delta}}_{2}$ has to increase as well, and such choices benefit from an increased SNR in the DWMRI data. The Connectome 1.0 DWMRI data was obtained using $\displaystyle {\mathrm{\Delta}}_{1}=19\text{ms}$ and $\displaystyle {\mathrm{\Delta}}_{2}=49\text{ms}$, leading to a separation of $30\text{ms}$. For this data, it is expected that with SNR = 5 an $R}^{2$ around 0.9 is feasible, and by increasing SNR to 20, the $R}^{2$ can increase to a value above 0.99.
In Figure 3, we present the scatter plots of simulated $K$ values versus fitted $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 {\mathrm{\Delta}}_{1}=19\text{ms}$ (row 1) or $\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). $R}^{2$ values for each case at each SNR level are provided. This result verifies that subdiffusion based kurtosis estimation (blue dots) improves using multiple diffusion times. The improvement in $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 $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 $K$ by using the traditional DKI method (orange dots), especially with kurtosis values larger than 1. Overall, fitting subdiffusion model (Equation 3) to data with two adequately separated diffusion times can lead to robust estimation of mean kurtosis, via (Equation 9).
Timedependence in DKI metrics
In Figure 4, we show the timedependence effect of the DKI metrics ($D}_{DKI$ and $K}_{DKI$) after fitting the standard DKI model to our simulated data. In this fitting, we consider bvalues of 0, 1000, 1400, and 2500 s/mm^{2}, and vary the diffusion time, as in Jelescu et al., 2022. We depict the parameter estimates, $D}_{DKI$ and $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 longterm fitting results and trends in the parameter estimates. In both (A) and (C), as diffusion time increases, $D}_{DKI$ decreases as expected for an effective diffusion coefficient. The mean $D}_{DKI$ (averaged over 1000 simulations) agrees with the true diffusivity $D}^{\ast$. When it comes to the kurtosis $K}_{DKI$, in the noiseless data setting (D), we see $K}_{DKI$ is converging to the true kurtosis value $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 timedependence of diffusivity and kurtosis are provided in the Discussion.
Towards fast DKI data acquisitions
Next, we sought to identify the minimum number and combination of bvalues to use for mean kurtosis estimation based on the subdiffusion model (Equation 3). The simulation results were generated using the Connectome 1.0 DWMRI data bvalue setting, which has 16 bvalues, 8 from $\displaystyle \mathrm{\Delta}=19\phantom{\rule{thinmathspace}{0ex}}\text{ms}$ and 8 from $\displaystyle \mathrm{\Delta}=49\phantom{\rule{thinmathspace}{0ex}}\text{ms}$.
In Figure 5, we computed $R}^{2$ values for all combinations of choices when the number of nonzero bvalues is 2, 3, and 4. We then plotted the $R}^{2$ values sorted in ascending order for each considered SNR level. Notably, as the number of nonzero bvalues increase, the achievable combinations increase as well (i.e., 120, 560, and 1820 for the three different nonzero bvalue cases). The colours illustrate the proportion of bvalues used in the fitting based on $\displaystyle {\mathrm{\Delta}}_{1}$ and $\displaystyle {\mathrm{\Delta}}_{2}$. For example, 1:0 means only $\displaystyle {\mathrm{\Delta}}_{1}$ bvalues were used, and 2:1 means two $\displaystyle {\mathrm{\Delta}}_{1}$ and one $\displaystyle {\mathrm{\Delta}}_{2}$ bvalues were used to generate the result. The bvalue combinations achieving highest $R}^{2$ values are displayed in the inset pictures and the bvalues are provided in Table 1.
Our simulation results in Table 1 suggest that increasing the number of nonzero bvalues from two to four improves the quality of the parameter estimation, also achievable by fitting the subdiffusion model (Equation 3) to higher SNR data. The gain is larger by using higher SNR data than by using more bvalues. For example, going from two to four nonzero bvalues with SNR = 5 data approximately doubles the $R}^{2$, whereas the $R}^{2$ almost triples when SNR = 5 data is substituted by SNR = 20 data. Additionally, the use of $\displaystyle {\mathrm{\Delta}}_{1}$ or $\displaystyle {\mathrm{\Delta}}_{2}$ alone is not preferred (also see Figure 5), and preference is towards first using $\displaystyle {\mathrm{\Delta}}_{1}$ and then supplementing with bvalues generated using $\displaystyle {\mathrm{\Delta}}_{2}$. At all SNR levels when only two nonzero bvalues are used, one bvalue should be chosen based on the $\displaystyle {\mathrm{\Delta}}_{1}$ set, and the other based on $\displaystyle {\mathrm{\Delta}}_{2}$. Moving to three nonzero bvalues requires the addition of another $\displaystyle {\mathrm{\Delta}}_{1}$ bvalue, and when four nonzero bvalues are used then two from each diffusion time are required. If we consider an ${R}^{2}=0.90$ to be a reasonable goodnessoffit for the subdiffusion model, then at least three or four nonzero bvalues are needed with an SNR = 20. If SNR = 10, then three nonzero bvalues will not suffice. Interestingly, an $R}^{2$ of 0.85 can still be achieved when SNR = 20 and two optimally chosen nonzero bvalues are used.
Table 1 summarises findings based on having different number of nonzero bvalues with $R}^{2$ values deduced from the SNR = 5, 10, and 20 simulations. We have chosen to depict five bvalue combinations producing the largest $R}^{2$ values for the two, three, and four nonzero bvalue sampling cases. We found consistency in bvalue combinations across SNR levels. Thus, we can conclude that a range of bvalues can be used to achieve a large $R}^{2$ value, which is a positive finding, since parameter estimation does not stringently rely on bvalue sampling. For example, using three nonzero bvalues an ${R}^{2}\ge 0.90$ is achievable based on different bvalue sampling. Importantly, two distinct diffusion times are required, and preference is towards including a smaller diffusion time bvalue first. Hence, for three nonzero bvalues we find two bvalues based on $\displaystyle {\mathrm{\Delta}}_{1}$ and one based on $\displaystyle {\mathrm{\Delta}}_{2}$. This finding suggests one of the $\mathrm{\Delta}}_{1$ bvalues can be chosen in the range 50 to 350 s/mm^{2}, and the other in the range 1500 to 4750 s/mm^{2}. Additionally, the $\mathrm{\Delta}}_{2$ bvalue can also be chosen in a range, considering between 2300 to 4250 s/mm^{2} based on the Connectome 1.0 bvalue settings. The bvalue sampling choices made should nonetheless be in view of the required $R}^{2$ value. Overall, sampling using two distinct diffusion times appears to provide quite a range of options for the DWMRI data to be used to fit the subdiffusion model parameters. The suggested optimal bvalue sampling in the last row of Table 1, primarily chosen to minimise bvalue size whilst maintaining a large $R}^{2$ value, may be of use for specific neuroimaging studies, which will be used to inform our discussion on feasibility for clinical practice.
Benchmark mean kurtosis in the brain
The benchmark mean kurtosis estimation in the brain is established using the entire bvalue range with all diffusion encoding directions available in the Connectome 1.0 DWMRI data. For two subjects in different slices, Figure 6 provides the spatially resolved maps of mean kurtosis computed using the subdiffusion method (i.e., $K}^{\ast$) with one or two diffusion times, and using the standard method (i.e., $K}_{DKI$) considering the two distinct diffusion times. First, we notice a degradation in the $K}_{DKI$ image with an increase in diffusion time. Second, the use of a single diffusion time with the subdiffusion model leads to $K}^{\ast$ values which are larger than either the $K}_{DKI$ values or $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 $K}^{\ast$. Superior greywhite matter tissue contrast (TC) was found for the $K}^{\ast$ map ($TC=1.73$), compared to the $K}_{DKI$ maps ($TC=0.80$ for the $\displaystyle \mathrm{\Delta}=19\text{ms}$ dataset and $TC=1.01$ for the $\displaystyle \mathrm{\Delta}=49\text{ms}$ dataset).
In Figure 7, an error map (measured by rootmeansquareerror, RMSE) from Subject 5 slice 74 was presented for fitting the subdiffusion model to the DWMRI data with two diffusion times. Sample parameter fittings in both bspace (Equation 3) and qspace (Equation 4) were provided for four representative white and grey matter voxels.
Quantitative findings for kurtosis are provided in Table 2. The analysis was performed for subcortical 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 \mathrm{\Delta}$ results in a slight decrease in the mean $K}_{DKI$, and an increase in the coefficient of variation (CV) for any region. For example, the mean $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 $K}^{\ast$ value for each region is less than the CV for $K}_{DKI$ with either $\displaystyle \mathrm{\Delta}=19\text{ms}$ or $\mathrm{\Delta}=49\text{ms}$.
Figure 8 presents the distributions of the fitted parameters ($D$ and $K$) in specific brain regions, based on the subdiffusion model (Panel A) and the standard DKI model with $\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 \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, $D}_{DKI$ and $K}_{DKI$, in all regions considered. By comparison, the subdiffusion based $K}^{\ast$ and $D}^{\ast$ (Panel A) are less correlated with each other, indicating $D}^{\ast$ and $K}^{\ast$ carry distinct information, which will be very valuable for characterising tissue microstructure. Tables 3 and Table 4 present the $D}_{\beta$ and β values used to compute $D}^{\ast$ and $K}^{\ast$ in Table 2 and Figure 8.
Reduction in DKI data acquisition
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 $K}^{\ast$ maps in Figure 6 and the $K}^{\ast$ values reported in Table 2.
Figure 9 presents the qualitative findings for two subjects generated using all, optimal and suboptimal bvalue samplings with SNR = 6 (3 noncollinear directions, 6 measurements), 10 (8 noncollinear directions, 16 measurements), and 20 (32 noncollinear directions, 64 measurements) DWMRI data. The quality of the mean kurtosis map improves with increasing SNR, and also by optimising bvalue sampling. Optimal sampling at SNR = 10 is qualitatively comparable to the SNR = 20 optimal sampling result, and to the benchmark subdiffusion results in Figure 6.
Quantitative verification of the qualitative observations are provided in Table 5. Significant differences in brain regionspecific mean kurtosis values occur at the SNR = 6 level, which are not apparent when SNR = 10 or 20 data with optimal bvalue sampling were used. The average errors are relative errors compared to the benchmark kurtosis values reported in Table 2. When using optimal bvalues, 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 suboptimal bvalues, 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 regionspecific CV for mean kurtosis was not found to change markedly when SNR = 10 or 20 data were used to compute $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 bvalue sampling, the data quality can be reduced to around the SNR = 10 level, without a significant impact on the regionspecific mean kurtosis estimates derived using the subdiffusion model.
Scan–rescan reproducibility of mean kurtosis
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 DWMRI (i.e., SNR = 20; subscript A) and the SNR = 10 optimal bvalue sampling scheme (subscript O). As the value of μ approaches 1, the intersubject variation in mean kurtosis is expected to greatly outweigh the intrasubject 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 $\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 $K}^{\ast$ value. Irrespective of which of the two DWMRI data were used for $K}^{\ast$ estimation, the value of μ was greater than or equal to 0.70 in 20 out of 24 cases. The $\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 bvalue sampling went handinhand with an increase in $\sigma$, which is not unexpected, since the uncertainty associated with using less data should be measurable. Overall, $\mu}_{A$, $\mu}_{O$, and $\sigma}_{A$, $\sigma}_{O$, were fairly consistent across the brain regions, suggesting the DWMRI data with SNR = 10 is sufficient for mean kurtosis estimation based on the subdiffusion framework.
Discussion
DWMRI allows the measurement of mean kurtosis, a metric for the deviation away from standard Brownian motion of water in tissue, which has been used to infer variations in tissue microstructure. Research on mean kurtosis has shown benefits in specific applications over other diffusion related measures derived from DWMRI data (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). Whilst many efforts have been made to optimise mean kurtosis imaging for clinical use, the limitations have been associated with lack of robustness and the time needed to acquire the DWMRI data for mean kurtosis estimation. The choice of the biophysical model and how diffusion encoding is applied are critical to how well kurtosis in the brain is mapped. Here, we evaluated the mapping of mean kurtosis based on the subdiffusion model, which allows different diffusion times to be incorporated into the data acquisition. Using simulations and the Connectome 1.0 public DWMRI dataset, involving a range of diffusion encodings, we showed that mean kurtosis can be mapped robustly and rapidly provided at least two different diffusion times are used and care is taken towards how bvalues are chosen given differences in the SNR level of different DWMRI acquisitions.
Reduction in scan time
Previous attempts have been made in optimising the DWMRI 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 bshells, directions per shell, and subsampling of DWMRI 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 bvalue of around $2500\text{s}/{\text{mm}}^{2}$ (Jensen et al., 2005; Jensen and Helpern, 2010), and using the subdiffusion 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 subdiffusion model applied to an optimised DWMRI protocol. Optimisation of the data acquisition with the use of the subdiffusion 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 bshell takes a fixed amount of time, then a four bshell acquisition with six directions per shell will take 25 times longer than a single diffusion encoding data acquisition (assuming a single bvalue=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^{2} scans plus seven bvalues 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 $3.88\text{s}$ to acquire. Conservatively allowing $4\text{s}$ per scan, and considering SNR = 20 data (i.e., 64 directions) over four bvalues and a single bvalue = 0 scan, DWMRI data for mean kurtosis estimation can be completed in $17\text{min}\phantom{\rule{thinmathspace}{0ex}}8\text{s}$ (${R}^{2}=0.96$). At SNR = 10 (i.e., 32 directions), DWMRI data with the same number of bvalues can be acquired in $4\text{min}\phantom{\rule{thinmathspace}{0ex}}20\text{s}$ (${R}^{2}=0.91$). If an ${R}^{2}=0.85$ (SNR = 10) is deemed adequate, then one less bshell is required, saving an additional $64\text{s}$. We should point out that even though 2shell optimised protocols can achieve ${R}^{2}=0.85$ with SNR = 20, this is not equivalent in time to using 3shells with SNR = 10 (also ${R}^{2}=0.85$). This is because 4×$4\times$ additional data are required to double the SNR (equivalent to acquiring an additional 4shells). However, only one extra bshell is required to convert 2shell data to 3shells with SNR = 10. Our expected DWMRI data acquisition times are highly feasible clinically, where generally neuroimaging scans take around $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 $1\text{min}$ over the brain (Hansen et al., 2013). Clinical adoption of the protocol lacked, possibly since bvalues are a function $\delta$, $G$, and $\mathrm{\Delta}$. Hence, different bvalues can be obtained using different DWMRI protocol settings, leading to differences in the mapping of mean kurtosis based on the data (we showed the $\displaystyle \mathrm{\Delta}$ effect in Figure 6 and Table 2). Our findings suggest this impediment is removed by sampling and fitting data with bvalues across two distinct diffusion times. Nonetheless, we should consider what might be an acceptable DWMRI data acquisition time.
A recent study on nasopharyngeal carcinoma investigated reducing the number of bshell signals based on fixing diffusion encoding directions to the three Cartesian orientations (He et al., 2022). The 3shell acquisitions took $3\text{min}2\text{s}$ to acquire, whilst the 5shell 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 kspace line acquisitions for each diffusion encoded image. Their benchmark used 5shells (200, 400, 800, 1500, 2000 s/mm^{2}), and found partial Fourier sampling with omission of the 1500 s/mm^{2} bshell 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 3shells acquisition with an ${R}^{2}=0.85$ (see SNR = 10 results in Table 1) executable under $4\text{min,}$ is therefore inline with current expectations. Note, at the ${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 DWMRI data to the Subject 1 first scan. This approach results in slight mismatch of the regionspecific segmentations when carried across subjects, inherently resulting in an underestimation of ICC values.
Less than $4\text{min,}$ DWMRI 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 DWMRI 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 bshell, multiple diffusion encoding direction DWMRI acquisitions should be used for the calculation of both DTI and DKI metrics.
DWMRI data acquisition considerations
To achieve ${R}^{2}>0.92$ for estimating kurtosis $K}^{\ast$, it is necessary to have four bvalues, e.g., two relatively small bvalues (350 s/mm^{2} using $\displaystyle \mathrm{\Delta}=19\text{ms}$, and 950 s/mm^{2} using $\mathrm{\Delta}=49\text{ms}$, both with $G=68\text{mT}/\text{m}$) and two larger bvalues of around 1500 s/mm^{2} and 4250 s/mm^{2} (using $\displaystyle \mathrm{\Delta}=19\text{ms}$ and $\mathrm{\Delta}=49\text{ms}$, respectively, both with $G=142\text{mT}/\text{m}$) (see bottom row in Table 1). If two or three nonzero bvalues are considered sufficient (with ${R}^{2}=0.85$ or $0.90$), then the larger $\mathrm{\Delta}$ needs to be used to set the largest bvalue to be 2300 s/mm^{2}, and the other(s) should be set using the smaller $\displaystyle \mathrm{\Delta}$. For the two nonzero bvalues case, the bvalue from the smaller $\mathrm{\Delta}$ should be around 800 s/mm^{2}. For the three nonzero bvalues case, the bvalues from the smaller $\mathrm{\Delta}$ would then need to be 350 s/mm^{2} and 1500 s/mm^{2}. Interestingly, bvalue = 800 s/mm^{2} lies around the middle of the 350 s/mm^{2} to 1500 s/mm^{2} range. The additional gain to ${R}^{2}=0.92$ can be achieved by splitting the bvalue with the larger $\mathrm{\Delta}$ into two, again with $2300\text{s}/{\text{mm}}^{2}$ near the middle of the two new bvalues set. In addition, the separation between $\mathrm{\Delta}}_{1$ and $\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 signaltonoise ratio decreases with increased echo times (He et al., 2022).
A recent study on optimising quasidiffusion imaging (QDI) considered bvalues up to 5000 s/mm^{2} (Spilling et al., 2022). Whilst QDI is a nonGaussian approach, it is different to the subdiffusion model, but still uses the MittagLeffler function and involves the same number of model parameters. The DWMRI 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 bvalue used to create the DWMRI data. They also showed the accuracy, precision, and reliability of parameter estimation are improved with increased number of bshells. They suggested a maximum bvalue of 3960 s/mm^{2} for the 4shell parameter estimation. Our results do not suggest a dependence of the parameter estimates on maximum bvalue (note, if a maximum bvalue 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 DWMRI data and not the magnitude as commonly used). These findings potentially confirm that $\mathrm{\Delta}$ separation is an important component of obtaining a quality parameter fit from which mean kurtosis is deduced (Figure 2).
Diffusion gradient pulse amplitudes
Commonly available human clinical MRI scanners are capable of $40\text{mT/m}$ gradient amplitudes. Recently, the increased need to deduce tissue microstructure metrics from DWMRI 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 bvalues 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 (q,\mathrm{\Delta})$ space by providing a mechanism for increasing our knowledge of the relationships between the micro, meso, and macroscales. 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 bvalue is an intricate combination of $\delta$, $G$, and $\displaystyle \mathrm{\Delta}$. An increase in any of these three parameters increases the bvalue (note, $\delta$ and $G$ increase bvalue quadratically, and $\mathrm{\Delta}$ linearly). An increase in $\delta$ can likely require an increase in $\mathrm{\Delta}$, the consequence of which is a reduction in signaltonoise 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 subsampling the DWMRI 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 $G$, as for the Connectome 1.0 and 2.0 scanners.
Based on our suggested bvalue settings in Table 1, a maximum bvalue of around 4250 s/mm^{2} is required for robust mean kurtosis estimation (assume ${R}^{2}>0.90$ is adequate and achieved via four nonzero bvalues and two distinct Δs, and we should highlight that bvalue = 1500 s/mm^{2} ($\displaystyle \mathrm{\Delta}=19\text{ms}$) and bvalue = 4250 s/mm^{2} ($\displaystyle \mathrm{\Delta}=49\text{ms}$) were achieved using $G=142\text{mT/m}$). Considering 80 m/Tm gradient sets are the new standard for MRI scanners, an adjustment to $\delta$ to compensate for $G$ is needed (recall, $q=\gamma \delta G$). Hence, for a constant $q$, $G$ can be reduced at the consequence of proportionally increasing $\delta$. This then allows MRI systems with lower gradient amplitudes to generate the relevant bvalues. For example, changing the $\delta$ from 8 ms to 16 ms would result in halving of $G$ (using a maximum $G$ of 142 m/TM from the Connectome 1.0 can achieve the optimal bvalues; hence halving would require around 71 m/TM gradient pulse amplitudes). Increasing $\mathrm{\Delta}$ above 49 ms to create larger bvalue 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 $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).
Relationship between subdiffusion based mean kurtosis $K}^{\ast$ and histology
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 subcortical 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 subdiffusion based mean kurtosis $K}^{\ast$ is consistent with published histology of normal human brain.
Timedependence of diffusivity and kurtosis
The timedependence of diffusivity and kurtosis has attracted much interest in the field of tissue microstructure imaging. Although our motivation here is not to map timedependent diffusion, we can nonetheless point out that the assumption of a subdiffusion model provides an explanation of the observed timedependence of diffusivity and kurtosis. From (Equation 4), the diffusivity that arises from the subdiffusion model is of the form $D}_{SUB}={D}_{\beta}{\overline{\mathrm{\Delta}}}^{\beta 1$, where $\overline{\mathrm{\Delta}}$ is the effective diffusion time, $D}_{SUB$ has the standard units for a diffusion coefficient $\text{mm}}^{2}\text{/s$, and $D}_{\beta$ is the anomalous diffusion coefficient (with units $\text{mm}}^{2}{\text{/s}}^{\beta$) associated with subdiffusion. In the subdiffusion framework (Equation 1), $D}_{\beta$ and $\beta$ are assumed to be constant, and hence $D}_{SUB$ exhibits a fractional powerlaw dependence on diffusion time. Then, following (Equation 8), $D}^{\ast$ is obtained simply by scaling $D}_{SUB$ by a constant $1/\mathrm{\Gamma}(1+\beta )$, and hence also follows a fractional powerlaw dependence on diffusion time. This timedependence effect of diffusivity was illustrated in our simulation results, Figure 4(A) and (C).
When it comes to kurtosis, the literature on the timedependence 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 subdiffusion 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, $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 timedependency 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 timedependence of kurtosis are quite mixed.
Furthermore, we summarise the benefits of using the subdiffusion based mean kurtosis measurement $K}^{\ast$. First, as shown in (Equation 9), subdiffusion based mean kurtosis $K}^{\ast$ is not time dependent, and hence has the potential to become a tissuespecific imaging biomarker. Second, the fitting of the subdiffusion model is straightforward, fast and robust, from which the kurtosis $K}^{\ast$ is simply computed as a function of the subdiffusion model parameter β. Third, the kurtosis $K}^{\ast$ is not subject to any restriction on the maximum bvalue, as in standard DKI. Hence, its value truly reflects the information contained in the full range of bvalues.
Extension to directional kurtosis
The direct link between the subdiffusion 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 DWMRI dataderived 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 DWMRI 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 DWMRI sampling over a large number of diffusion encoding directions within each bshell (Jensen and Helpern, 2010; Poot et al., 2010). As such, extension to directional kurtosis requires a larger DWMRI dataset acquired using an increased number of diffusion encoding directions. The number of bshells and directions therein necessary for robust and accurate mapping of directional kurtosis based on the subdiffusion 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 bshell in the DWMRI 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 rank4 kurtosis tensor is estimated from the DWMRI 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 rank2 diffusion tensor with six unique tensor entries is needed to be estimated, whilst in DKI in addition to the rank2 diffusion tensor, the kurtosis tensor is rank4 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 DWMRI data (including bvalue = 0) with different diffusion encoding properties have to be acquired (Jensen and Helpern, 2010). The traditional approach has been to set five distinct bvalues with 30 diffusion encoding directions within each bshell (Poot et al., 2010). Hence, to obtain the entries of the rank4 kurtosis tensor, much more DWMRI 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 DWMRI data acquisition times.
A kurtosis tensor framework based on the subdiffusion 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 rank2 β tensor with only six unknowns, requiring at least six distinct diffusion encoding directions. This type of approach can reduce the amount of DWMRI data to be acquired, and potentially serve as a viable way forward for the combined estimation of mean, axial, and radial kurtosis.
Kurtosis estimation outside of the brain
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 monoexponential model can be recovered from the subdiffusion equation by setting $\beta =1$. For this case, the product between the bvalue and fitted diffusivity has been reported to be insightful for bvalue 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 subdiffusion equation, we can investigate the size of $b{D}_{SUB}$ by analysing the four nonzero bvalue optimal sampling regime ($\mathrm{\Delta}}_{1$: 350 s/mm^{2} and 1500 s/mm^{2}; $\mathrm{\Delta}}_{2$: 950 s/mm^{2 }and 4250 s/mm^{2} from Table 1). Considering scGM, cGM, and WM brain regions alone, the rounded and dimensionless $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 $\mathrm{\Delta}}_{1$, and the other two are derived using $\mathrm{\Delta}}_{2$. Interestingly, the loglinear 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 $R}^{2$ across the entire brain, since β and $D}_{SUB$ are brain region specific and the most optimal sampling strategy should be β and $D}_{SUB$ specific. Whilst region specific sampling may provide further gains in the $R}^{2$ value, and improve ICC values for specific brain regions, such data would take a long time to acquire and require extensive postprocessing and indepth analyses.
Methods
Theory
Subdiffusion modelling framework
In biological tissue, the motion of water molecules is hindered by various microstructures, and hence the diffusion can be considerably slower than unhindered, unrestricted, free diffusion of water. The continuous time random walk formalism provides a convenient mathematical framework to model this subdiffusive behaviour using fractional calculus (Metzler and Klafter, 2000). The resulting probability density function $P(x,t)$ of water molecules at location $x$ (in units of mm) at time $t$ (in units of s) satisfies the time fractional diffusion equation:
where $\frac{{\mathrm{\partial}}^{\beta}}{\mathrm{\partial}{t}^{\beta}}$ is the time fractional derivative of order $\beta$ ($0<\beta \le 1$) in the Caputo sense, $D}_{\beta$ is the generalised anomalous diffusion coefficient with unit of $\text{mm}}^{2}/{\text{s}}^{\beta$, and the parameter β characterises the distribution of waiting times between two consecutive steps in the continuous time random walk interpretation. When $\beta =1$, the waiting times have finite mean; when $0<\beta <1$, the waiting times have infinite mean, leading to subdiffusion behaviour. The solution to the time fractional diffusion Equation 1 in Fourier space is:
where $E}_{\beta}(z)=\sum _{n=0}^{\mathrm{\infty}}\frac{{z}^{n}}{\mathrm{\Gamma}(1+\beta n)$ is the singleparameter MittagLeffler function, $\mathrm{\Gamma}$ is the standard Gamma function and by definition ${E}_{1}(z)=\mathrm{exp}(z)$. In the context of diffusion DWMRI, $k$ in Equation 2 represents the qspace parameter $q=\gamma G\delta$, $t$ represents the effective diffusion time $\overline{\mathrm{\Delta}}=\mathrm{\Delta}\delta /3$ and $p(k,t)$ represents the signal intensity $S(q,\overline{\mathrm{\Delta}})$, leading to the diffusion signal equation (Magin et al., 2020):
Defining $\displaystyle b={q}^{2}\overline{\mathrm{\Delta}}$, the DWMRI signal then can be expressed in terms of bvalues:
where
has the standard unit for a diffusion coefficient, s/mm^{2}.
Diffusional kurtosis imaging
The traditional DKI approach was proposed by Jensen and colleagues (Jensen et al., 2005; Jensen and Helpern, 2010) to measure the extent of nonGaussian diffusion in biological tissues using DWMRI data:
where $S$ is the signal for a given diffusion weighting $b$ (i.e., bvalue), $S}_{0$ is the signal when $b=0$, $D}_{DKI$ and $K}_{DKI$ are the apparent diffusivity and kurtosis. A major limitation of (Equation 6) is that it was developed based on the Taylor expansion of the logarithm of the signal at $b=0$, as such bvalues should be chosen in the neighbourhood of bvalue = 0 (Yang et al., 2022; Kiselev, 2010). Hence, to estimate diffusivity and kurtosis, Jensen and Helpern, 2010 suggested the use of three different bvalues (such as 0, 1000, 2000 s/mm^{2}) and the maximum bvalue should be in the range 2000 s/mm^{2} to 3000 s/mm^{2} for brain studies. Subsequently, the optimal maximum bvalue was found to be dependent on the tissue types and specific pathologies, which makes the experimental design optimal for a whole brain challenging (Chuhutin et al., 2017). The procedure for fitting kurtosis and diffusivity tensors is also not trivial, and a variety of fitting procedures are currently in use. We refer readers to the descriptive and comparative studies for detail on the implementation and comparison of methods (Veraart et al., 2011b; Chuhutin et al., 2017).
Mean kurtosis from the subdiffusion model
Yang et al., 2022 established that the traditional DKI model corresponds to the first two terms in the expansion of the subdiffusion model:
where diffusivity, $D}^{\ast$, and kurtosis, $K}^{\ast$, are computed directly via subdiffusion parameters $D}_{SUB$ and $\beta$:
where $D}_{SUB$ is defined in (Equation 5). Note the mean kurtosis expression in (Equation 9) was also derived by Ingo et al., 2015 using a different method. Their derivation follows the definition of kurtosis, $K=\u27e8{x}^{4}\u27e9/\u27e8{x}^{2}{\u27e9}^{2}3$, i.e., by computing the fourth moment $\u27e8{x}^{4}\u27e9$ and the second moment $\u27e8{x}^{2}\u27e9$ based on the subdiffusion Equation 1.
Connectome 1.0 human brain DWMRI data
The DWMRI 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 $2\times 2\times 2{\text{mm}}^{3}$ resolution data were obtained using two diffusion times $({\displaystyle \mathrm{\Delta}=19,49\phantom{\rule{thinmathspace}{0ex}}\text{ms})}$ with a pulse duration of $\delta =8\phantom{\rule{thinmathspace}{0ex}}\text{ms}$ and $G=31,68,105,142,179,216,253,290\phantom{\rule{thinmathspace}{0ex}}\text{mT/m}$, respectively, generating bvalues = 50, 350, 800, 1500, 2400, 3450, 4750, 6000 s/mm^{2} for $\mathrm{\Delta}=19\phantom{\rule{thinmathspace}{0ex}}\text{ms}$, and bvalues = 200, 950, 2300, 4250, 6750, 9850, 13500, 17800 s/mm^{2} for $\displaystyle \mathrm{\Delta}=49\phantom{\rule{thinmathspace}{0ex}}\text{ms}$, according to bvalue $=(\gamma \delta G{)}^{2}(\mathrm{\Delta}\delta /3)$. 32 diffusion encoding directions were uniformly distributed on a sphere for $b<2400\phantom{\rule{thinmathspace}{0ex}}\text{s}/{\text{mm}}^{2}$ and 64 uniform directions for $b\ge 2400\phantom{\rule{thinmathspace}{0ex}}\text{s}/{\text{mm}}^{2}$.
The FreeSurfer’s segmentation labels as part of the dataset were used for brainregionspecific analyses. Tian et al., 2022 provided the white matter averaged group SNR (23.10 ± 2.46), computed from 50 interspersed bvalue $=0\text{s}/{\text{mm}}^{2}$ images for each subject. Both magnitude and the real part of the DWMRI were provided. Based on an indepth analysis, the use of the real part of the DWMRI data was recommended, wherein physiological noise, by nature, follows a Gaussian distribution (Gudbjartsson and Patz, 1995).
Simulated DWMRI data at specific bvalues
DWMRI data were simulated to establish (i) the correspondence between actual versus fitted mean kurtosis using the traditional DKI and subdiffusion models based on various choices for $\mathrm{\Delta}$, and (ii) to investigate the impact of SNR levels and subsampling of bvalues on the mean kurtosis estimate. The DWMRI signal was simulated using the subdiffusion model (Equation 3) with random Gaussian noise added to every normalised DWMRI signal instance:
where $N(0,{\sigma}^{2})$ is white noise with mean of zero and standard deviation of $\sigma$ according to the normal distribution.
Two aspects influence $\sigma$ in the case of realvalued DWMRI data. These include the SNR achieved with a single diffusion encoding direction (i.e., Connectome 1.0 DWMRI data SNR was derived using only bvalue $=0\text{s}/{\text{mm}}^{2}$ data), and the number of diffusion encoding directions in each bshell across which the powder average is computed:
where $N}_{\mathrm{D}\mathrm{I}\mathrm{R}$ is the number of diffusion encoding directions for each bshell and assuming it is consistent across bshells. The $\sigma$ for the Connectome 1.0 data is approximately $1/(23.10\times 8)=0.0054$ based on 64 diffusion encoding directions. Achieving of SNR = 5, 10, and 20 for the simulation study can therefore be accomplished by changing only the $\sigma$ and keeping ${N}_{\mathrm{D}\mathrm{I}\mathrm{R}}=64$. As such, ${\sigma}_{\mathrm{S}\mathrm{N}\mathrm{R}=5}=0.0250$, ${\sigma}_{\mathrm{S}\mathrm{N}\mathrm{R}=10}=0.0125$, and ${\sigma}_{\mathrm{S}\mathrm{N}\mathrm{R}=20}=0.0063$.
Three simulation experiments were carried out at various SNR levels. The first simulation experiment was to examine the effect of the number of diffusion times on the accuracy of the parameter fitting for idealised grey matter and white matter cases. In (Equation 10) the choices of $\displaystyle {D}_{\beta}=3\times {10}^{4}\phantom{\rule{thinmathspace}{0ex}}{\text{mm}}^{2}/{\text{s}}^{\beta}$, $\beta =0.75$, and $\displaystyle {D}_{\beta}=5\times {10}^{4}\phantom{\rule{thinmathspace}{0ex}}{\text{mm}}^{2}/{\text{s}}^{\beta}$, $\beta =0.85$, were made for white matter and grey matter, respectively. These two distinct βs led to $K}^{\ast$ of 0.8125 and 0.4733 using (Equation 9). Diffusion times were chosen from the range ${\mathrm{\Delta}}_{i}\in [\delta ,\delta +50]$, where diffusion pulse length was set to $\delta =8\text{ms}$ to match the Connectome 1.0 data. A minimum required separation between any two $\mathrm{\Delta}$s was enforced, i.e., $30/(n1)$, where 30 corresponds with the $30\text{ms}$ difference between the $\mathrm{\Delta}$s for the Connectome 1.0 DWMRI data, and $n$ is the number of distinct diffusion times simulated. We considered as many as five distinct diffusion times. Simulations were conducted by randomly selecting sets of $\mathrm{\Delta}$s for 1000 instances, and then generating individual DWMRI simulated signals using (Equation 10), before fitting for $D}_{\beta$ and β, from which $K}^{\ast$ was computed using (Equation 9). For the case of two diffusion times, suggestions on the separation between them was given based on the goodnessoffit of the model.
The second simulation experiment was to investigate the effect of the number of diffusion times on the accuracy of parameter fitting. To generate simulated data reflecting observations in human data, $D}_{\beta$ and $\beta$ were restricted to ${D}_{\beta}\in [{10}^{4},{10}^{3}]$ in the unit of $\text{mm}}^{2}/{\text{s}}^{\beta$ and $\beta \in [0.5,1]$, corresponding to ${K}^{\ast}\in [0,1.7124]$. For the case of two diffusion times, suggestions on the separation between them were given based on the goodnessoffit of the model.
The third simulation experiment is to study bvalue subsampling under various SNR levels (SNR = 5, 10, and 20). We set the bvalues in the simulated data the same as those used to acquire the Connectome 1.0 dataset. At each SNR level, we selected combinations of two, three, and four bvalues, irrespective of the difference between them and the diffusion time set to generate the bvalue. Essentially, we explored the entire possible sets of bvalues for the three regimes, resulting in 120, 560, and 1820 combinations, respectively. A goodnessoffit measure for model fitting was used to make comparisons between the different bvalue combinations.
SNR reduction by downsampling diffusion encoding directions
We performed SNR reduction of the Connectome 1.0 DWMRI data by downsampling of diffusion directions in each bshell. The method of multiple subsets from multiple sets (PDMM) subsampling algorithm provided in DMRITool (Cheng et al., 2018) was applied to the bvectors provided with the Connectome 1.0 DWMRI data. Note, the bshells contained 32 directions if $b<2400\text{s}/{\text{mm}}^{2}$, and 64 directions if $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 noncollinear diffusion encoding directions in each bshell, 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 bvalue, and SNR = 6 had six).
Parameter estimation
Standard DKI model
The maximum bvalue used to acquire the DWMRI for standard DKI model fitting is limited to the range 2000 s/mm^{2} to 3000 s/mm^{2} due to the quadratic form of (Equation 6). We opted to use DWMRI data generated with bvalues $=50,350,800,1500,2400\text{s}/{\text{mm}}^{2}$ using $\mathrm{\Delta}=19\text{ms}$, and bvalues = 200, 950, 2300 s/mm^{2} using $\mathrm{\Delta}=49\text{ms}$. Note, the apparent diffusion coefficient, $D}_{DKI$ in (Equation 6) is time dependent, as can be deduced from (Equation 5) and (Equation 8). Thereby, standard DKI fitting can only be applied to DWMRI data generated using a single diffusion time. The model in (Equation 6) was fitted in a voxelwise manner to the powder averaged (i.e., geometric mean over diffusion encoding directions, often referred to as traceweighted) DWMRI data using the lsqcurvefit function in MATLAB (Mathworks, Version 2022a) using the trustregion reflective algorithm. Optimisation function specific parameters were set to TolFun = 10^{–4} and TolX = 10^{–6}. Parameters were bounded to the ranges of ${D}_{DKI}>0$ and $0<{K}_{DKI}\le 3$.
Subdiffusion model
For the single diffusion time case, the subdiffusion model in (Equation 3) was fitted to the powder averaged DWMRI data in a voxelwise manner using the same MATLAB functions as in the previous section. For each subject, spatially resolved maps of $D}_{\beta$ and $\beta$ were generated. The fitting strategy for multiple diffusion time data is to solve:
where $S}_{ij$ is the signal at the ith effective diffusion time $\overline{\mathrm{\Delta}}}_{i$ and the jth qspace parameter $q}_{j$, $SUB$ is the subdiffusion model (Equation 3) at $\overline{\mathrm{\Delta}}}_{i$ and $q}_{j$ for a given set of (${D}_{\beta},\beta$), $n$ is the number of diffusion times, and $J}_{i$ is the number of $q$values corresponding to $\overline{\mathrm{\Delta}}}_{i$ in data acquisition. This objective function allows incorporation of an arbitrary number of diffusion times, each having arbitrary number of $q$values. Parameters were bounded to the ranges of $0<\beta \le 1$ and ${D}_{\beta}>0$ for all the parameter fittings. Model parameters were found to be insensitive to the choice of initial values. Parameters $D}^{\ast$ and $K}^{\ast$ were computed analytically using the estimated $D}_{\beta$ and $\beta$ according to (Equation 8) and (Equation 9).
Goodnessoffit and regionbased statistical analysis
For all of the simulation results, the coefficient of determination, referred to as $R}^{2$, was used to assess how well the mean kurtosis values were able to be fitted. This value was calculated using a standard $R}^{2$ formula
where $\overline{K}}_{\text{Simulated}$ is the mean simulated $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 $TC=\left{\mu}_{WM}{\mu}_{GM}\right/\sqrt{{\sigma}_{WM}^{2}+{\sigma}_{GM}^{2}}$, where $\mu}_{WM$ and $\mu}_{GM$ are the mean parameter values in white and grey matters; and $\sigma}_{WM$ and $\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 regionsofinterest. We considered three categories of brain regions, namely subcortical 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 ttest was performed to test for significant differences in mean kurtosis population means under the optimal and suboptimal bvalue sampling regimes presented in Table 5 compared to the benchmark parameter values presented in Table 2.
Scan–rescan analysis using intraclass correlation coefficient (ICC)
For each of the six subjects, both the first (scan) and second (rescan) scan images were registered to the first scan images of Subject 1 using inbuilt MATLAB (Mathworks, Version 2022a) functions (imregtform and imwarp). We used 3D affine registration to account for distortions and warps common in DWMRI data. Cubic spline interpolation was applied to resample both the scan and rescan DWMRI data for each subject onto Subject 1’s first scan data grid. The FreeSurfer labels for Subject 1’s first scan were used for brain region analysis. For each voxel, the ICC measure was applied to assess scan–rescan reproducibility of mean kurtosis, as described by Duval et al., 2018 and Fan et al., 2021,
with ${s}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}}^{2}={\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}_{i}\left(1/2\cdot (({m}_{i,scan}{\overline{m}}_{i}{)}^{2}+({m}_{i,rescan}{\overline{m}}_{i}{)}^{2})\right)$ and ${s}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}^{2}={\mathrm{v}\mathrm{a}\mathrm{r}}_{i}({\overline{m}}_{i})$. Here, $m}_{i,scan$ is the kurtosis value measured in the voxel for subject i, and $\overline{m}}_{i$ is the average value between scan and rescan. An ICC histogram and the mean and standard deviation descriptive statistics were generated for all brain regions analysed.
Conclusion
The utility of DKI for inferring information on tissue microstructure was described decades ago. Continued investigations in the DWMRI 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 subdiffusion 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.
Data availability
The Connectome 1.0 human brain DWMRI data used in this study is part of the MGH Connectome Diffusion Microstructure Dataset (CDMD) (Tian et al., 2022), which is publicly available on the figshare repository: https://doi.org/10.6084/m9.figshare.c.5315474. MATLAB codes generated for simulation study, parameter fitting, and optimisation of bvalue sampling are openly available at https://github.com/mfarquhar/SubdiffusionDKI (copy archived at Farquhar, 2023).

figshareComprehensive diffusion MRI dataset for in vivo human brain microstructure mapping using 300 mT/m gradients.https://doi.org/10.6084/m9.figshare.c.5315474
References

Diffusiontime dependence of diffusional kurtosis in the mouse brainMagnetic Resonance in Medicine 84:1564–1578.https://doi.org/10.1002/mrm.28189

Single and multipleshell uniform sampling schemes for diffusion mri using spherical codesIEEE Transactions on Medical Imaging 37:185–199.https://doi.org/10.1109/TMI.2017.2756072

Scanrescan of axcaliber, macromolecular tissue volume, and gratio in the spinal cordMagnetic Resonance in Medicine 79:2759–2765.https://doi.org/10.1002/mrm.26945

SoftwareSubdiffusionDKI, version swh:1:rev:d1b3d199a21f9861bd103bb2bf413e91ca38c8fcSoftware Heritage.

Mapping early tumor response to radiotherapy using diffusion kurtosis imaging*The Neuroradiology Journal 36:198–205.https://doi.org/10.1177/19714009221122204

The Rician distribution of noisy MRI dataMagnetic Resonance in Medicine 34:910–914.https://doi.org/10.1002/mrm.1910340618

Experimentally and computationally fast method for estimation of a mean kurtosisMagnetic Resonance in Medicine 69:1754–1760.https://doi.org/10.1002/mrm.24743

Diffusion imaging in humans at 7T using readoutsegmented EPI and GRAPPAMagnetic Resonance in Medicine 64:9–14.https://doi.org/10.1002/mrm.22480

Toward more robust and reproducible diffusion kurtosis imagingMagnetic Resonance in Medicine 86:1600–1613.https://doi.org/10.1002/mrm.28730

Fast diffusion kurtosis imaging in acute ischemic stroke shows mean kurtosisdiffusivity mismatchJournal of Neuroimaging 32:941–946.https://doi.org/10.1111/jon.13000

Exponential analysis in physical phenomenaReview of Scientific Instruments 70:1233–1257.https://doi.org/10.1063/1.1149581

Diffusional kurtosis imaging: the quantification of nongaussian water diffusion by means of magnetic resonance imagingMagnetic Resonance in Medicine 53:1432–1440.https://doi.org/10.1002/mrm.20508

MRI quantification of nonGaussian water diffusion by kurtosis analysisNMR in Biomedicine 23:698–710.https://doi.org/10.1002/nbm.1518

BookThe cumulant expansion: an overarching mathematical framework for understanding diffusion NMRDiffusion MRI Oxford University Press.https://doi.org/10.1093/med/9780195369779.003.0010

7T: Physics, safety, and potential clinical applicationsJournal of Magnetic Resonance Imaging 46:1573–1589.https://doi.org/10.1002/jmri.25723

Advanced fit of the diffusion kurtosis tensor by directional weighting and regularizationMagnetic Resonance in Medicine 67:1401–1411.https://doi.org/10.1002/mrm.23133

Diffusion tensor imaging: concepts and applicationsJournal of Magnetic Resonance Imaging 13:534–546.https://doi.org/10.1002/jmri.1076

The diagnostic value of diffusion kurtosis imaging in Parkinson’s disease: a systematic review and metaanalysisAnnals of Translational Medicine 10:474.https://doi.org/10.21037/atm221461

Diffusion kurtosis imaging as an imaging biomarker for predicting prognosis in chronic kidney disease patientsNephrology Dialysis Transplantation 37:1451–1460.https://doi.org/10.1093/ndt/gfab229

Fractional calculus models of magnetic resonance phenomena: relaxation and diffusionCritical Reviews in Biomedical Engineering 48:285–326.https://doi.org/10.1615/CritRevBiomedEng.2020033925

Signaltonoise ratio and MR tissue parameters in human brain imaging at 3, 7, and 9.4 tesla using current receive coil arraysMagnetic Resonance in Medicine 75:801–809.https://doi.org/10.1002/mrm.25677

Optimal experimental design for diffusion kurtosis imagingIEEE Transactions on Medical Imaging 29:819–829.https://doi.org/10.1109/TMI.2009.2037915

Optimization of quasidiffusion magnetic resonance imaging for quantitative accuracy and timeefficient acquisitionMagnetic Resonance in Medicine 88:2532–2547.https://doi.org/10.1002/mrm.29420

Estimation of tensors and tensorderived measures in diffusional kurtosis imagingMagnetic Resonance in Medicine 65:823–836.https://doi.org/10.1002/mrm.22655

Diffusion MRI in the brain – theory and conceptsProgress in Nuclear Magnetic Resonance Spectroscopy 112–113:1–16.https://doi.org/10.1016/j.pnmrs.2019.03.001

Constrained maximum likelihood estimation of the diffusion kurtosis tensor using a Rician noise modelMagnetic Resonance in Medicine 66:678–686.https://doi.org/10.1002/mrm.22835

More accurate estimation of diffusion tensor parameters using diffusion Kurtosis imagingMagnetic Resonance in Medicine 65:138–145.https://doi.org/10.1002/mrm.22603

MR diffusion kurtosis imaging for neural tissue characterizationNMR in Biomedicine 23:836–848.https://doi.org/10.1002/nbm.1506

Theoretical models of the diffusion weighted MR signalNMR in Biomedicine 23:661–681.https://doi.org/10.1002/nbm.1520

Specific absorption rate studies of the parallel transmission of innervolume excitations at 7TJournal of Magnetic Resonance Imaging 28:1005–1018.https://doi.org/10.1002/jmri.21548
Article and author information
Author details
Funding
Australian Research Council (Discovery Early Career Research Award DE150101842)
 Qianqian Yang
Australian Research Council (Discovery Project Award DP190101889)
 Qianqian Yang
 Viktor Vegh
The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.
Acknowledgements
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 MRIbased 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.
Ethics
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/s41597021010926.
Version history
 Preprint posted:
 Sent for peer review:
 Reviewed Preprint version 1:
 Reviewed Preprint version 2:
 Version of Record published:
Cite all versions
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.
Copyright
© 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.
Metrics

 319
 views

 24
 downloads

 1
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
 Physics of Living Systems
Explaining biodiversity is a fundamental issue in ecology. A longstanding puzzle lies in the paradox of the plankton: many species of plankton feeding on a limited variety of resources coexist, apparently flouting the competitive exclusion principle (CEP), which holds that the number of predator (consumer) species cannot exceed that of the resources at a steady state. Here, we present a mechanistic model and demonstrate that intraspecific interference among the consumers enables a plethora of consumer species to coexist at constant population densities with only one or a handful of resource species. This facilitated biodiversity is resistant to stochasticity, either with the stochastic simulation algorithm or individualbased modeling. Our model naturally explains the classical experiments that invalidate the CEP, quantitatively illustrates the universal Sshaped pattern of the rankabundance curves across a wide range of ecological communities, and can be broadly used to resolve the mystery of biodiversity in many natural ecosystems.

 Chromosomes and Gene Expression
 Computational and Systems Biology
Genes are often regulated by multiple enhancers. It is poorly understood how the individual enhancer activities are combined to control promoter activity. Anecdotal evidence has shown that enhancers can combine subadditively, additively, synergistically, or redundantly. However, it is not clear which of these modes are more frequent in mammalian genomes. Here, we systematically tested how pairs of enhancers activate promoters using a threeway combinatorial reporter assay in mouse embryonic stem cells. By assaying about 69,000 enhancerenhancerpromoter combinations we found that enhancer pairs generally combine nearadditively. This behaviour was conserved across seven developmental promoters tested. Surprisingly, these promoters scale the enhancer signals in a nonlinear manner that depends on promoter strength. A housekeeping promoter showed an overall different response to enhancer pairs, and a smaller dynamic range. Thus, our data indicate that enhancers mostly act additively, but promoters transform their collective effect nonlinearly.