On the need for bundle-specific microstructure kernels in diffusion MRI.

Probing microstructure with diffusion magnetic resonance imaging (dMRI) on a scale orders of magnitude below the imaging resolution relies on biophysical modelling of the signal response in the tissue. The vast majority of these biophysical models of diffusion in white matter assume that the measured dMRI signal is the sum of the signals emanating from each of the constituent compartments, each of which exhibits a distinct behaviour in the b-value and/or orientation domain. Many of these models further assume that the dMRI behaviour of the oriented compartments (e.g. the intra-axonal space) is identical between distinct fibre populations, at least at the level of a single voxel. This implicitly assumes that any potential biological differences between fibre populations are negligible, at least as far as is measurable using dMRI. Here, we validate this assumption by means of a voxel-wise, model-free signal decomposition that, under the assumption above and in the absence of noise, is shown to be rank-1. We evaluate the effect size of signal components beyond this rank-1 representation and use permutation testing to assess their significance. We conclude that in the healthy adult brain, the dMRI signal is adequately represented by a rank-1 model, implying that biologically more realistic, but mathematically more complex fascicle-specific microstructure models do not capture statistically significant or anatomically meaningful structure, even in extended high-b diffusion MRI scans.

Introduction

Probing tissue microstructure with in vivo magnetic resonance imaging (MRI) on a scale orders of magnitude below the imaging resolution relies on biophysical modelling of the signal response in intra- and extra-cellular space. In diffusion-weighted MRI (dMRI), the signal is sensitized to the thermal motion of water molecules on length scales that commensurate with various cell dimensions (Le Bihan et al., 1986; Callaghan et al., 1988). Hence, dMRI provides a non-invasive, yet indirect, probe into the tissue micro-environment, despite the coarse spatial resolution of MRI (Alexander et al., 2019; Novikov et al., 2019). Analysis of dMRI data on the macroscopic (voxel) level has shown sensitivity to several microstructural changes, associated with development, ageing, and disorder of the human brain, and of white matter in particular (Neil et al., 2002; Horsfield and Jones, 2002; Sundgren et al., 2004; Sullivan and Pfefferbaum, 2006). However, a more specific interpretation of the signal in terms of brain microstructure relies on the development and validation of a suitable biophysical model (Jelescu and Budde, 2017).

The vast majority of these biophysical white matter models derive from a model of multiple non-exchanging compartments that contribute linearly to the measured dMRI signal. Each compartment is assumed to exhibit a distinct behaviour in the b-value and/or orientation domain, characterised by a kernel or response function that encapsulates the expected signal behaviour that would be measured for that compartment in isolation. Many models further impose the constraint that the various compartments in a single voxel share the same orientation structure; for example, the extra-cellular space may be modelled as a distinct compartment with its own signal characteristics, but oriented identically to the intra-axonal compartment. Under this assumption, this shared orientation structure can be captured in a single orientation distribution , and the microstructural characteristics can be modelled by a single voxel-level kernel , whose signal response is the mixture of the constituent compartment responses.

Nearly all of the most common biophysical white matter models hence subsume an overarching spherical convolution model (Tournier et al., 2004; Anderson, 2005; Jespersen et al., 2007), in which the dMRI signal

The voxel-wise microstructure kernel describes the signal response in straight, parallel white matter fascicles, as an axially-symmetric function of the b-value and of the cosine θ of the angle between the fibre orientation and the diffusion gradient direction . This voxel-level response can encapsulate the mixture of all compartments that share the same orientational structure, as well as any isotropic compartments. The (unnormalized) orientation distribution function captures the direction and dispersion of the different fascicles in the voxel, and is fundamental in tractography.

Different biophysical models adopt different constraints and functional forms for and , many of which have been heavily debated in the community (Jelescu and Budde, 2017; Novikov et al., 2018a), but can ultimately be written in the form of Eq. (1). The microstructure kernel is typically modelled as a mixture of restricted, hindered and fast water compartments, each parameterized by their axial and radial diffusivities and signal fractions (Behrens et al., 2003, 2007; Kroenke et al., 2004; Assaf and Basser, 2005; Jespersen et al., 2007, 2010; Alexander et al., 2010; Fieremans et al., 2011; Panagiotaki et al., 2012; Zhang et al., 2012; Ferizi et al., 2015; Jelescu et al., 2016; Kaden et al., 2016; Reisert et al., 2017; Novikov et al., 2018b). Other approaches have adopted empirical forms for the response kernel, in some cases invariant across the brain (Tournier et al., 2004; Anderson, 2005; Tournier et al., 2007; Dell’Acqua et al., 2007; Jian and Vemuri, 2007; Descoteaux et al., 2009; Jeurissen et al., 2014; Christiaens et al., 2017). The orientation structure is often represented as an orientation distribution function (ODF) in the spherical harmonics (SH) basis (Frank, 2002), as it is the most general representation and reduces the convolution in (1) to a multiplication of corresponding SH coefficients (Tournier et al., 2004, 2007; Jespersen et al., 2007, 2010; Reisert et al., 2017; Novikov et al., 2018b). Other models have selected functional descriptors of the orientation structure, for instance in the form of δ-functions or Watson-distributions representing one or several fibre populations and their dispersion (Behrens et al., 2003, 2007; Zhang et al., 2012; Sotiropoulos et al., 2012).

The large combinatorial space of different model choices for and , not to mention the many different fitting strategies, has led to a plethora of microstructure imaging techniques. The quest for the right model on the effective length scale in dMRI therefore critically depends on validating the underlying model assumptions in vivo (Novikov et al., 2018a). Most validation studies have focused on the functional form of the microstructure kernel (Dhital et al., 2018; Veraart et al., 2019; Tax et al., 2019), but until now the assumed voxel-wise uniqueness of this kernel and its corresponding ODF has not been evaluated.

Indeed, many of the most general models rely on the widely adopted assumption that, within an individual voxel, the microstructure-specific aspects can be expressed as a single kernel that can adequately represent all fascicles in that voxel. This implicit assumption of the convolution model in Eq. (1) imposes a clear simplification of the biology given that microstructural features might be fascicle-specific and that the coarse voxel size in dMRI inevitably introduces partial volume effects. Consider for example the fibre configuration sketched in Fig. 1, where two non-parallel fibre bundles have sufficiently different microstructural properties leading to them having different kernels and . In this case, the total signal can generally not be described as a single convolution. To describe such a system requires a more generalized model with multiple convolutions, in this instance , unless either both kernels or both ODFs are identical up to a normalization factor.

Fig. 1

Sketch of a voxel where the microstructure is not described by a single rotation-invariant kernel, as most models assume. The blue, vertical white matter fascicles are modelled as the spherical convolution of microstructure kernel and ODF . The orthogonal fascicles in orange are also modelled as a convolution (of and ), but with more restrictive diffusion resulting in a different microstructure kernel. This work validates if in vivo adult brain dMRI data can adequately be described as the convolution of a single microstructure kernel and a single ODF .

Note that, as per the equation above, multi-compartment models reduce to a single convolution if all directional compartments are coupled by a single ODF. Any isotropic compartments are trivially absorbed into a single kernel too. As such, the ball-and-stick model (Behrens et al., 2003, 2007), NODDI (Zhang et al., 2012), and other common 2- and 3-compartment models (Novikov et al., 2018b; Reisert et al., 2017), all abide by the single-convolution assumption. Multi-tissue constrained spherical deconvolution (MT-CSD) (Jeurissen et al., 2014) does too, as long as all but one response functions are isotropic. On the other hand, multi-tensor models (Tuch et al., 2002) and the CHARMED model (Assaf and Basser, 2005) where intra- and extra-axonal compartments are uncoupled in orientation are examples of models that do not assume the single-convolution form.

In this work, we will evaluate how well the in vivo dMRI signal in adult white matter is described by a single per-voxel response in even the most complex fibre configurations. Our approach is based on a voxel-wise decomposition into spherical harmonics and radial orthonormal components that, under the single-convolution model and in the absence of noise, offers a rank-1 descriptor of the dMRI signal (Christiaens et al., 2019). Testing the model rank in vivo can therefore validate the assumption of a single response per voxel, without imposing any prior constraints or functional form on the response or on the ODF. We will evaluate the relative effect size of higher-order signal components in rich multi-shell data of 4 healthy adult volunteers, and determine their significance with permutation testing. Source code and data are publicly available at https://osf.io/mjbnz.

Materials and methods

Data and preprocessing

Multi-shell dMRI data were acquired in 4 healthy adults (3 male, 1 female) on a 3T Philips Achieva system with 32-channel head coil. The multi-shell encoding scheme was designed to comprise 6 shells at , 1600, 2500, 3600, 4900, and (linear in ) with 70 isotropically-distributed directions each and 30 volumes, all interleaved for minimal duty cycle (Hutter et al., 2017). Each volume is acquired with isotropic resolution 2.5 ​mm, matrix size slices, , , multiband factor 2, and SENSE factor 2. We also acquired a matched single-band reference scan with AP-PA and LR-RL phase encoding.

The dMRI data were first processed with patch-level image denoising based on random matrix theory in complex (magnitude and phase) data, which accounts for the noise correlations in most common dMRI examinations (Cordero-Grande et al., 2019). Since this step operates in the complex domain, effects of Rician bias are minimized in subsequent magnitude images. Mean SNR of the raw data (before denoising) was measured to be at in white matter. The data was subsequently filtered for Gibbs ringing suppression (Kellner et al., 2016). The reverse phase encoding reference scan was used to estimate the B0 field map (Andersson et al., 2003), which was then used as part of subject motion and image distortion correction (Andersson and Sotiropoulos, 2016). Brain masks were segmented with FSL (Smith, 2002).

Voxel-wise decomposition

The spherical convolution (1) is conventionally projected onto the orthonormal basis of real, symmetric spherical harmonics of order and index . Here, we use a linear least-squares fit in the SH basis of order in all experiments. In this basis, and without loss of generality, convolution with a single response reduces to a multiplicationwhere are SH coefficients representing the dMRI signal in shell b and are SH coefficients representing the ODF. The axially-symmetric response kernel is represented by coefficients in the basis of zonal spherical harmonics , or equivalently by scaled coefficients that encompass the convolution normalization factor to simplify the notation.

With multi-shell data, we can arrange the voxel-wise signal coefficients of corresponding SH order in matriceswhere k shells are laid out in rows and their SH coefficients are laid out in columns. The SH convolution (2) can then be written as the outer product of vectors (indexed in b) and (indexed in m). As such, with a single per-voxel response and in the absence of noise, we haveand by construction, the matrices have matrix rank (Christiaens et al., 2019). If, on the other hand, the signal is not adequately described by a single rotation-invariant response, we can introduce a generalized model with multiple kernels and their corresponding ODFs:and where are generally of rank . Therefore, the matrix rank of can reveal the presence of multiple microstructure kernels for different orientations in the voxel.

In the special case where all kernels are identical up to a scale factor (), equation (5) reduces to the single-convolution model (with ) due to distributivity of the convolution and matrix product operations. Similarly, when all ODFs are identical up to a scale factor () the generalized model also reduces to the form of (4) (in this case with ), which can for instance arise when all fascicles are parallel. Furthermore, we also note that in all voxels because matrices have only one column. This property reflects the observation that isotropic compartments can always be absorbed in a single per-voxel kernel: the SH band offers no capability to disentangle multiple independent kernels and ODFs.

In this work, we use the singular value decomposition (SVD) of matrices as a means to determine the number of significant microstructure components. Specifically,where and are the left and right singular vectors and are the singular values. This decomposition is a voxel-wise specialization of the spherical harmonics and radial decomposition (SHARD) representation introduced in Christiaens et al. (2019), and will generally be of full rank due to inevitable noise and biological variability in the data. By comparing this SVD to (4), it is clear that the leading left and right singular vectors and retrieve — up to a normalization factor — the predominant microstructure response and ODF respectively. The subsequent components span the orthogonal space, i.e., the axes of variability not captured in the leading component. The associated singular values offer a measure of the effect size of this unexplained variance, thus making it possible to determine which components contain meaningful structure.

Relative signal explained

The squared singular values measure the signal power of each component i in the SH band of order . In order to retrieve a single measure of each component’s total effect size, we compute the signal power across all harmonic bands

Its square root measures the root-mean-squared (RMS) signal in each component; it has the same units as the signal. The ratiothen provides a measure of the signal power explained in the leading component. Under the rank-1 model, this ratio is expected to be near 100%. Mapping this measure across the brain thus enables to evaluate the proportion of the signal that is explained in the single-convolution model. Conversely, measures the local residual signal not explained with a single spherical convolution.

Permutation testing

In order to evaluate if there is evidence for multiple (fascicle-specific) kernels in a voxel, we wish to assess the effect size of additional components in relation to a distribution of unstructured residuals. We formulate the null hypothesis as

Components contain no meaningful structure.

We wish to evaluate the null hypothesis with as few prior assumptions as possible. To this end, we deploy a permutation testing strategy in which the null distribution is sampled numerically using residual bootstrapping. This process destroys any structure in the residuals of a rank-1 model by inducing random permutations and applying the appropriate leverage. Specifically, we first calculate the rank-1 representation in each voxel and its residuals to the dMRI signal . Subsequently, we sample bootstrap instances, each of which is generated with the following procedure:

Apply a random permutation to the residuals and correct for leverage with a factor h (see below):

Construct bootstrap instance ;

Fit bootstrap instance into the SH basis, construct matrices and their singular values ;

Store the component effect sizes of the bootstrap instance.

The histogram of determines the null distribution of the corresponding component i. The number of instances where the actual effect size in the data exceeds that of the randomized bootstrap, i.e., where , determines the voxel-wise p-valuesof components . then offers the primary metric of potential structure beyond the single-convolution model.

Correction for leverage is needed to rescale the variance of the residuals to the effective degrees of freedom (Cook and Weisberg, 1982). The mean squared residuals of a least-squares fit of a model with κ degrees of freedom to ν data points scales as , where is the noise variance in the data. Hence, the root-mean-squared error of a least squares fit underestimates the actual noise level, and an inverse scaling is applied to compensate for this in permutation testing. For a well conditioned system, a leverage factor achieves the appropriate scaling (Davidson, 2006). In our case, the total no. degrees of freedom in the model ( is the no. shells) and the no. data points , so . We also verified in simulations that this definition approximates the effective noise scaling in a rank-1 SVD.

The per-voxel p-values are subsequently corrected for multiple comparisons using the Benjamini-Hochberg false discovery rate (FDR) procedure (Benjamini and Hochberg, 1995). FDR procedures control the proportion of false positives relative to the number of positive results, and are generally less stringent than family-wise error rate procedures such as Bonferroni correction. The specific FDR procedure used here (Benjamini-Hochberg, which assumes independent or positively correlated tests) is itself less stringent than other FDR approaches (Benjamini and Yekutieli, 2001). We adopt significance level and .

Note that the test statistic in components is affected by SH terms only, since the SH term is always rank-1. Therefore, only anisotropic signal content (the contrast in the angular domain) contributes to the statistical evaluation.

Results

Fig. 2a displays the RMS signal power in the leading component. The total RMS power in the other components is shown in Fig. 2b. We observe that the leading component captures the main signal in the parenchyma. The other components, while not free of anatomical structure, primarily capture the overall noise modulation and image artefacts including a posterior fat shift artefact and flow or pulsation effects near the brain stem. Fig. 2c shows the ratio R of signal power explained in the leading component across the brain. We find that the rank-1 model explains of the signal in the parenchyma, and observe no consistent reduction of this ratio in regions of crossing fibres. Similar figures for the other subjects are provided in supplementary materials and support this observation. A histogram of the % signal power explained in the rank-1 representation across all 4 subjects, shown in Fig. 3, confirms that the effect of additional components in a generalized model is small.

Fig. 2

Signal power (ratio) in the leading component in Subject 4. (a) The root-mean-squared (RMS) signal power in the principal voxel-wise SHARD component. (b) The total RMS signal power in components 2 to r, i.e., . (c) The ratio R of the signal power explained in the first component. Supplementary Figures for the other subjects are provided online.

Fig. 3

Histogram of the signal power explained in the first component across the full brain and across all 4 subjects (blue), and its survival function (; orange curve). The plot shows that the data is well represented by the leading component alone.

To determine if these additional components — however small — contain meaningful structure, we study the results of the permutation testing in Fig. 4. Fig. 4a shows the voxel-wise p-value of component 2 as estimated with the non-parametric residual bootstrapping procedure. We observe that in scattered voxels with little discernible anatomical structure. After correcting for multiple comparisons in Fig. 4b, voxels with a significant 2nd component are increasingly sparse, with consistent anatomical structure only occurring at the interface between the ventricles and the corpus callosum where we do not expect multi-fascicle contributions. Visual inspection of the preprocessed dMRI data showed evidence of ventricular CSF pulsation, and previous studies have also pointed to effects of physiological noise in these regions (Walker et al., 2011). In components , no voxels survive multiple comparisons. Corresponding figures for the other subjects are provided as supplementary material. Generally, permutation testing did not find evidence for higher-order components, indicating that the single-convolution model is an appropriate description for these data.

Fig. 4

Permutation testing in Subject 2 to determine if component 2 contains any meaningful structure. (a) p-value of a voxel-wise residual bootstrapping test. (b) Significant voxels (; shown in red) overlaid onto the mean image. Significance is determined after correcting for multiple comparisons using the Benjamini-Hochberg False Discovery Rate procedure with . Supplementary Figures for the other subjects are provided online.

Discussion

This paper has introduced a model-independent means of evaluating the single convolution assumption implicit to most biophysical models for white matter diffusion MRI, i.e., the assumption that the signal is sufficiently described with a single ODF and a single rotation-invariant response function per voxel. Our approach is independent of the functional form of the microstructure response and independent of the complexity of the fibre orientation distribution, relying instead on a matrix decomposition shown to be rank-1 under the single-convolution model. Results in adult brain data show that the effect size of any signal components beyond the single-convolution model is small ( of the signal) and that these additional components are not statistically significant when corrected for multiple comparisons. Therefore, even though different fascicles in a voxel may have different microstructural characteristics, we find no evidence of their observability in the diffusion MRI signal that would warrant the added complexity of multi-fascicle models in healthy adult white matter across the broad range of b-values sampled in these data.

In the given experimental setup, this conclusion stands for healthy in vivo adult brain diffusion MRI data acquired on a clinical system. Whether the single-convolution assumption is also sufficient to describe ex vivo and post mortem data, pediatric brain data, or pathology, remains an open question. Similarly, while our acquisition protocol was designed to collect rich multi-shell data well beyond clinically feasible scan time (47 ​min), higher SNR, even longer acquisitions or stronger MRI gradient systems might reveal evidence for additional components. In these cases, the presented method can be used to evaluate if the single-convolution assumption is appropriate for the subject group and data at hand. Given that the data in this study supersedes what is considered clinically practical, the evidence suggests that we can build on the single-convolution model to extract clinically applicable biomarkers of adult white matter microstructure.

Of note is that the voxel-wise decomposition is also directly applicable to generalized diffusion encoding waveforms, and to diffusion time and echo time dependent contrasts. Indeed, any orientation-dependent contrast can be projected in the SH basis and incorporated as additional rows in matrices (3), subject to the same rank-1 decomposition in the single-convolution model. Including generalized diffusion encoding or time dependence might hence add complementary contrast to detect and recover multiple kernels and ODFs, and may be the subject of future work. Here, we prioritised linear diffusion MRI across an extended b-value range as the most clinically relevant protocol to investigate at this point.

With the null result in mind, we now shine some light on how different the kernels need to be to break the single-convolution assumption at a given SNR. In Fig. 5, we show contour plots to demonstrate the SNR required to detect differences in standard model parameters between both fascicles (see also Appendix A: Simulations). With an asymptotic SNR after denoising, the adult brain data in this work only allows the detection of strong differences between kernels, either in the intra-axonal fraction and diffusivity or in combinations of several parameters. This finding highlights that the single-convolution model is sufficient in most practical applications, even if the biological tissue properties are fascicle-specific.

Fig. 5

Simulations of the observability of kernel differences between two equally-weighted fascicles in a crossing. The curves plot the significance levels and of the bootstrapping test as a function of signal-to-noise ratio (SNR) and of differences in “standard model” parameters. For instance, when changing the intra-axonal fraction f, two kernels are simulated as . The presence of multiple kernels is detected when their differences are large enough and when SNR if sufficiently high.

Finally, while the presented model-free decomposition has a clear advantage encompassing a wide array of microstructure models, specific functional forms and constraints in model-based approaches may boost their statistical power. Multi-tensor models (Tuch et al., 2002; Scherrer and Warfield, 2010) are a case in point, where our generalized model is constrained to Gaussian responses and to ODFs in the form of SH δ-functions pointing along each fascicle direction. Whilst the multi-tensor model often reveals different diffusivities for the different fascicles in a voxel, these are driven in part by the fibre dispersion in each fascicle and in part by the non-mono-exponential signal response of the microstructure. Other examples include MT-CSD approaches with multiple anisotropic ODFs (Pietsch et al., 2019) that constrain the per-tissue responses to be invariant across the brain, and tractography-informed approaches (Reisert et al., 2014; Girard et al., 2017) that regularize microstructure parameters along fibre tracts. By using a local decomposition and deliberately avoiding functional forms, this work generalizes across a broad class of microstructure models. We can, however, not exclude that non-negativity or other sparsity-enforcing constraints on the ODF may reveal evidence for fascicle-specific microstructure kernels.

If, based on the evidence in this work, we accept the single-convolution assumption as valid, the proposed local SHARD decomposition has rank ​= ​1. This has interesting consequences for various applications. For instance, we could achieve voxel-level dMRI denoising, and potentially also outlier removal, by suppressing all but the leading component. Another interesting application can be in microstructure imaging based on rotation-invariant signal features (Reisert et al., 2017; Novikov et al., 2018b). Not only is the single-convolution assumption a prerequisite for using rotation-invariants in this way, the coefficients of the rank-1 decomposition can themselves form a set of rotation-invariant features useful in fitting a biophysical signal model.

Conclusion

While microstructural features are likely fascicle-specific, the diffusion MRI data in the white matter is well described by the spherical convolution of a single fibre orientation distribution and a single response function per voxel, regardless of their form and complexity. Our results suggest that we can build upon the current models of diffusion in adult white matter, despite the implicit single convolution assumption, and that a more complex model would not significantly improve the goodness-of-fit.