Apparent propagator anisotropy from single‐shell diffusion MRI acquisitions

Purpose The apparent propagator anisotropy (APA) is a new diffusion MRI metric that, while drawing on the benefits of the ensemble averaged propagator anisotropy (PA) compared to the fractional anisotropy (FA), can be estimated from single‐shell data. Theory and Methods Computation of the full PA requires acquisition of large datasets with many diffusion directions and different b‐values, and results in extremely long processing times. This has hindered adoption of the PA by the community, despite evidence that it provides meaningful information beyond the FA. Calculation of the complete propagator can be avoided under the hypothesis that a similar sensitivity/specificity may be achieved from apparent measurements at a given shell. Assuming that diffusion anisotropy (DiA) is nondependent on the b‐value, a closed‐form expression using information from one single shell (ie, b‐value) is reported. Results Publicly available databases with healthy and diseased subjects are used to compare the APA against other anisotropy measures. The structural information provided by the APA correlates with that provided by the PA for healthy subjects, while it also reveals statistically relevant differences in white matter regions for two pathologies, with a higher reliability than the FA. Additionally, APA has a computational complexity similar to the FA, with processing‐times several orders of magnitude below the PA. Conclusions The APA can extract more relevant white matter information than the FA, without any additional demands on data acquisition. This makes APA an attractive option for adoption into existing diffusion MRI analysis pipelines.


Isotropic Ensamble Average Propagator
In the main document, section 3.1, we explicitly calculate the inner product that defines the PA by using the simplification in eq. (9), yielding an anisotropy metric related to the PA for a specific shell, namely the APA. After [1], this metric requires the definition of an isotropic signal equivalent to the mono-exponential model. The rationale behind [1] is that the EAP can be averaged over the directional coordinates to obtain the closest isotropic signal to the original one. By using we have: where r = 1 and R = Rr. The inner integral in R 3 is usually computed in spherical coordinates, so that dq = q 2 du dq and a straightforward manipulation yields: where J 1/2 stands for Bessel's function of the first kind with index 1/2. The meaning of eq. (3) is that the isotropic EAP defined as the directional average of the original EAP becomes the Bessel transform of the isotropic diffusion signal defined as the directional average of the original diffusion signal. In other words, the isotropic equivalent to the EAP corresponds to a diffusion signal that is indeed computed in the same manner: By assuming a mono-exponential model of E(q) itself, we can compute: • The (squared) 2 norm of E(q): (5) • The (squared) 2 norm of E I (q): • The inner product between E(q) and E I (q): • The cosine between both two signals, which implicitly defines the PA: Using this complete approach, we achieve a a closed form of the APA that involves the computation of a quadratic form on the measured values at each voxel. Pursuing an analogous formulation to that in AMURA [2], in the main document we propose an alternative formulation leading to a linear computation with negligible deviations from the model. To that end, instead of using eq. (4) for the isotropic equivalent, we use the following expression instead:: with

Practical implementations
Note that the computation of E I 2 = E, E I in eq. (4) requires evaluating a double surface integral in the orientation variables u and v. In a practical implementation, such integrals are computed based on spherical harmonics expansions. In precise terms, let {u n } N n=1 be the set of the N acquired gradients within the measured shell; let {Y j (u)} M j=1 be the set of the M first (low order) spherical harmonics (typically M < N ). The coefficients {c j } M j=1 of a given orientation function, S(u), in this basis will be fitted as a Laplacian-regularized least squares problem: where the M × 1 vector C stacks the coefficients c j ; the N × 1 vector S stacks the measurements of the orientation function, S(u n ); the N × M matrix B stacks the values of the basis functions, B n,j = Y j (u n ); λ is the Laplace-Beltrami regularization parameter so that the M × M matrix L contains the eigenvalues of spherical harmonics for the Laplacian (we fix it to constant value of 0.006 in all cases). Since the 0−th order spherical harmonic encodes the DC component of the signal, the integral of the orientation function over the unit sphere reduces to a scaled version of its c 0 coefficient. For example, eq. (5) can be estimated as: where s stands for the first row of B T B + λL 2 −1 B T and the N × 1 vectorD stacks the N values of D(u n ) −3/2 . In order to compute eqs. (6) and (7), we arrange a N × N matrixD whose entries areD n 1 ,n 2 = (D(u n 1 ) + D(u n 2 )) −3/2 . Then: where the left product with s stands for the outermost integral in eq. (6) in the variable v for each constant u (each column u n 2 ), meanwhile the right product with s T stands for the innermost integral in eq. (6) in the variable u for each constant v (each row u n 1 ). This way, the computational complexity of computing simple surface integrals remains linear with the number of sampled gradients, O(N ), while the complexity of double surface integrals becomes O(N 2 ).

Comparison of both APA implementations
[ Figure 1 about here.] In order to compare the results given by the simplified APA implementation proposed in the main document to the complete implementation described in this supplementary material, an experiment is carried out. For the sake of comparison, we will denote full-APA (F-APA) to the complete implementation derived from eq. (8), and simply APA to the fast approach used along the paper, derived from eq. (9). For a visual comparison we use three slices from the HCP volume MGH1007 using a single shell for b=3000 s/mm 2 . We consider the absolute error as the quality measure: Error A mask is used in order to limit the measurement of the error to the white matter area. Both measures are bounded in the interval [0,1] and so will be the error. A visual comparison for slices (42, 52, 65) is shown in Fig. S1. The average error for the white matter is calculated for the whole volume at b=3000 s/mm 2 and b=5000 s/mm 2 and it can be found in Table S1. Note that the average error is smaller than 1% for both shells, which implies that both measures are practically the same and both provide very similar results. In addition, from the visual comparison in Fig. S1, we can also conclude that the error is uniformly distributed in the white matter, which, from a practical view point, it can be seen as a very small bias in the measure. Some of the experiments with real data from the main document were also redone for F-APA with no noticeable differences from those with APA.
[ Table 1 about here.] Although both methods provide similar results, in this work we have opted for APA, for a matter of simplicity. As previously stated, the computational complexity of computing simple surface integrals is O(N ), while the complexity of double surface integrals becomes O(N 2 ). Since APA is only based on simple integrals, similar results can be obtained in reduced computation time. In order to test it on real data, we measure the execution times of computing the metrics over the MGH1007 volume previously described. The experiment is run on a quad-core Intel(R) Core(TM) i7-4770K 3.50GHz processor under Ubuntu Linux 16.04 So using one single shell in MATLAB R2013b without multi-threading. The results are reported in Table S2. As a consequence of the simplicity on the implementation of APA, it shows execution times 200 faster than F-APA.
Thus, since the implementation error is so small and the gain in execution time is so high, we have opted to use the fast implementation thorough the whole paper. [