Spin movement observed with diffusion MRI is typically analyzed by means of the Ensemble Average Propagator (EAP), denoted by $$$P$$$. In the ideal case (under the narrow pulse assumption) $$$P(\textbf{r})$$$ represents the voxel-average probability that a spin undergoes a translation $$$\textbf{r}$$$ in the diffusion time $$$\Delta$$$, and it is related to the diffusion MRI signal $$$S$$$ of a pulsed-gradient spin-echo sequence through a Fourier transform¹: $$S(\textbf{q})=\int_\Omega\!e^{i\textbf{q}\cdot\textbf{r}}P(\textbf{r})\,dr.$$ Here $$$\textbf{q} = \gamma\delta\textbf{g}$$$ with $$$\gamma$$$, $$$\delta$$$, and $$$\textbf{g}$$$ respectively the gyromagnetic constant, the gradient pulse width, and the gradient.

The propagator provides access to quantitative information regarding the underlying diffusion process, such as the degree of Gaussianity. In this abstract we propose a method to estimate the orientation-dependent maximum range of spin movement, a quantity that can be of interest for example when searching for dominant axon orientations in the brain. We illustrate our approach by means of an analytical example and on brain diffusion MRI data.

It follows from the Shapley-Folkman-Starr Theorem² that (finite speed) spatially homogeneous stochastic processes, and thus the associated EAP, have a convex support function. If we write $$$\Omega$$$ for this support function (which is simply the set of points where $$$P(\textbf{r})$$$ is non-zero), then we find

$$H(\textbf{q}):=\log\left(S(-i\textbf{q})\right)=\|\textbf{q}\|\log\left[\left(\int_\Omega\!\exp\left(\frac{\textbf{q}}{\|\textbf{q}\|}\cdot\textbf{r}+\frac{1}{\|\textbf{q}\|}\log P(\textbf{r})\right)^{\|\textbf{q}\|}\right)^{\frac{1}{\|\textbf{q}\|}}\right]\leq\sup_{\textbf{r}\in\Omega}\textbf{q}\cdot\textbf{r}+\log P(\textbf{r}).$$

The extension of $$$S$$$ to complex-valued arguments is well-defined if $$$P$$$ has a moment-generating function (a mild assumption that is implicitly present in e.g. generalized Diffusion Tensor Imaging (DTI)), and the function $$$H$$$ is called the cumulant-generating function. In what follows we consider the large $$$\|\textbf{q}\|$$$ limit where the term $$$\log P(\textbf{r})$$$ can be neglected and the inequality reduces to an approximate equality. We then have the relation

$$H(\textbf{q})=\sup_{\textbf{r}\in\partial\Omega}\textbf{q}\cdot\textbf{r}.$$

From this expression we can infer the boundary $$$\partial\Omega$$$, the maximum range of spin movement, as the unit level set of the algebraic dual $$$F(\textbf{r}):=H^*(\textbf{r})$$$.

We reconstruct the function $$$H$$$ by fitting a generalized DTI model³ up to order $$$4$$$ to the signal, and substituting $$$\textbf{q}\to-i\textbf{q}$$$ before radially linearizing the function in the neighborhood of the set defined by $$$\log S(-i\textbf{q})=1$$$ (or generally some arbitrary constant). A numerical approximation of $$$F(\textbf{r})$$$ was achieved by evaluating the equality⁴ $$$H(\textbf{q})=F\left(\frac{1}{2}\nabla^2H^2(\textbf{q})\cdot\textbf{q}\right)$$$ for $$$246$$$ approximately uniformly distributed directions.

Practical constraints on the gradient hardware mean this method may significantly underestimate the range of movement even at high $$$b$$$-values, as illustrated by the toy example shown in Fig. 1. Instead we find an estimate of an observational boundary, indicating a range outside of which particles seldom traverse. We test this on a slice of a Human Connectome Project data set⁵. Estimated (voxel-wise) maximum displacements are shown in Fig. 2. The difference between the largest and smallest maximum displacements is depicted in Fig. 3, together with a standard FA map to illustrate the resemblance. Fig. 4 shows the maximum displacements per orientation as glyphs, while Fig. 5 shows the maximum relative displacements (i.e., the voxel-wise minimum subtracted for each orientation).

We introduced a new method that allows one to extract position- and orientation-dependent characteristic lengths that describe expected largest spin displacements (Fig. 1). The method relies on existing signal models, and can be applied to e.g. generalized DTI³ or diffusion kurtosis imaging⁶. The key assumption in the method is approximate spatial homogeneity of the voxel, and the few voxels with (severe) partial volume effects will significantly benefit from constrained reconstruction⁷.

We show in Fig. 2 that our method produces values within physically realistic and expected ranges in HCP data. Maps obtained by looking at the largest differences highlight areas where there is strong diffusion in a single direction, and show a strong resemblance to (DTI-based) fractional anisotropy. We note that local orientations along which movement is largest (Figs. 4 and 5) align relatively well with expected fiber orientations, and future work will focus on comparisons with e.g. spherical deconvolution results⁸.