Emmanuel Caruyer

Alessandro Daducci

Maxime Descoteaux

Jean-Christophe Houde

Jean-Philippe Thiran

Ragini Verma

Phantomas: a flexible software library to simulate diffusion MR phantoms

Target audience-We introduce Phantomαs, an open-source software library to create realistic phantoms in structural and diffusion MRI. This is intended as a tool for researchers working in method development or translational research in the field of diffusion-weighted MRI. The software allows evaluation of the whole tractography processing pipeline: acquisition design, preprocessing, local reconstruction and tractography algorithm. Purpose-Diffusion MRI has been an active area of research for the last two decades. However, while the data processing pipeline from DWIs to fiber tracking involves several steps, it is still difficult to evaluate the influence of a specific noise reduction method, a specific representation of the angular diffusion information, or a specific tracking algorithm, on fiber tracking results. By starting from a ground truth configuration of fiber bundle geometries, Phantomαs is designed to enable quantitative assessment of each step involved in the pipeline, using metrics based on successful reconstruction of the true connectivity profile of the phantom, such as defined in Tractometer 1 . Complementarily, for the evaluation of local reconstruction accuracy, local metrics 2 can be computed to measure the accuracy of angular models. Existing softwares 8,9 present similar objectives, yet Phantomas offers an original perspective in the computation of volume fraction (VF), signal attenuation and T1/T2 images. Methods-Geometric primitives. The phantom mimics a portion of the brain, with white matter fiber bundles connecting different cortical regions. Each fiber bundle is represented by a tube of constant radius wrapped around its centerline. The centerlines are defined as 3rd-order piecewise polynomials, controlled by user-defined control points. The centerlines are further constrained to be normal to the surface, which represents the gray-white matter interface. An example of a fiber configuration is represented in Fig. 1a. We also introduce some water-filled regions, represented as spheres, and a shell of water-filled material surrounding the phantom, intended to mimic the cerebrospinal fluid around the brain tissue. Tissue volume fraction. In order to simulate MR signal on a grid of given resolution, the first step is to identify the VF of each tissue type: white matter (WM), gray matter (GM) and cerebrospinal fluid (CSF). This is done in a hierarchical fashion: first, for each geometric primitive, we identify a map of voxels intersected by the corresponding object. Then for each voxel either partially filled by an object, or intersected by many objects, we evaluate the VF on a grid subdivision of the given voxel. This two-step procedure allows for an efficient computation of the VF for the whole dataset. Signal relaxation. The influence of longitudinal and transverse relaxation on the signal decay is also considered. Relaxation times reported in vivo 6 vary across individuals and across brain regions. We simulate maps of T1 and T2 relaxation times of the same resolution as the target image, using Gaussian distributed samples with the following mean and standard deviations (in ms): GM T1=1331±13, T2=110±2; WM T1=832±10, T2=79.6±0.6; CSF T1=3500±100, T2=250±10. A correlation ratio of 0.2 is introduced between direct neighboring voxels, which are not statistically independent. The image correlation matrix is factorized using Cholesky decomposition to simulate random image; since the dimension of the problem is large, but the covariance matrix is sparse, efficient algorithms 3 are used to solve the factorization problem. Fiber orientation distribution and signal attenuation. As for the simulation of tissue VF, the estimation of the fiber orientation distribution (FOD) in a given voxel starts by the computation of fiber orientation samples, corresponding to the closest orientation of each bundle on a regular subdivision of the voxel. Then, we use kernel density estimation with a symmetric Von-Mises Fischer kernel to estimate a continuous FOD. The signal attenuation in WM is computed based on a composite hindered and restricted model of diffusion 4 , where the extra-axonal compartment is modeled by a diffusion tensor and the intra-axonal compartment is modeled by a cylinder 7 . Since this model has cylindrical symmetry, the attenuation is computed as the convolution of the FOD by a template diffusion profile. This is done to save computational complexity, since the computation of the cylindrical restricted compartment is time and resource consuming. Results-We simulated a configuration of 27 fiber bundles, on a 70x70x70 grid, with 3 spherical regions of free water, represented in Fig. 1.a. The range of fiber configurations includes fiber bending, crossings, and kissing fibers. As a result, a range of intravoxel white matter configurations were created, with a total of 3695 voxels in pure WM, out of which 1192 contain more than one (and up to 5) fiber bundles. Besides, 15482 voxels contain a fraction of WM, and a fraction of GM and/or CSF. Crossing angles range from 0 to 90°, as reported in Fig. 1.c. The simulation of 64 volumes of DWIs (b = 3000s/mm²) were completed in about 15 minutes on a standard laptop. An example of fiber reconstruction using MRtrix 5 probabilistic tracking on constrained spherical deconvolution is reported in Fig. 1.d.

Figure 1 .

 1 Figure 1. a. An example of fiber bundle configuration, and spherical regions filled with free water. b. Corresponding slice of T1-weighted image. c. Ground truth fiber orientation distribution (zoom of the red square in 1.b). d. Fiber reconstructed with probabilistic streamline tracking using MRtrix 5 . Conclusions-We believe that Phantomαs will have a positive impact on the community for the validation and the development of novel methods in diffusion MRI. An early version of this tool was used successfully to create the testing dataset of the 2nd HARDI Reconstruction Challenge, organized at ISBI, San Francisco, CA, USA (2013). The software is released as a free and open-source software library written in Python and ANSI C, and is designed towards flexible use and integration with other software. References-1. M.-A. Côté et al., MedIA, 17(3), pp. 844-857, 2013. 2. A. Daducci et al., TMI, in press. 3. T. A. Davis and W. W. Hager, ACM Trans. Math. Software, 35, 2009. 4. Y. Assaf et al., Neuroimage 27(1), PP. 48-58, 2005. 5. J. Tournier et al., Int J Imag Syst Tech 22(1), pp. 53-66, 2012. 6. J. P. Wansapura et al., JMRI 9(4), pp. 531-538, 1999. 7. O. Söderman and B. Jönsson, JMR, A 117(1), pp. 94-97, 1995. 8. P. F. Neher et al., ISMRM, 20013. 9. T. G. Close et al., NeuroImage, 47(4), pp. 1288-1300, 2009.


Emmanuel Caruyer

Ragini Verma

Rotation-invariant measures for population study in HARDI

Target audience-We present a novel family of rotational-invariant measures for High Angular Resolution Diffusion Imaging 1 (HARDI) data. This research is motivated by need for the developing new biomarkers and thereby facilitating population-based studies in HARDI.

Purpose-From a spherical harmonics (SH) representation of the angular information in diffusion MRI, we derive rotational-invariant measures computed as homogeneous polynomials of the SH coefficients. When the SH is truncated to order , we find a total of 12 such invariants, which are expected to be sensitive to subtle changes in white matter structure, not captured by tensor fractional anisotropy (FA) or generalized fractional anisotropy (GFA). The method is generalizable to any spherical function, we report here the derivation and results based on the apparent diffusion coefficient 2 (ADC) profile.

Methods-Notations and problem statement. From the coefficients, of the truncated SH expansion of the ADC profile, we propose to search for all homogeneous polynomials, (where represents the vector of coefficients) that are invariant per rotation. We make use of the Wigner matrix, , associated to the rotation matrix . So for a given polynomial rank, t, we want to find all vectors of polynomial coefficients , such that , , . The dimensions for a given truncation order, L, and a given polynomial order, t, are N = (L + 1)(L + 2)/2, and . Homogeneous polynomial and linear transform. In order to solve the problem above, we try to derive the coefficients of the polynomial from the vector of coefficients . Albeit not reported here due to lack of space, we can find a linear transform, , such that , ,

. This means that applying a linear transform to the argument of the polynomial is equivalent to applying an equivalent transform to the coefficients of the polynomial. Then the problem resolves into finding such that , . Sufficient condition. Using arguments based on Euler angles decomposition and density of in , we can show that it is sufficient to solve the problem for the two particular rotations, and , of 1.0 rad about x and z axes, respectively. The problem reduces to a large linear system of equations; we use efficient solver for sparse systems 3 to find a solution. The set of solutions is further pruned to keep only a set of linearly and functionally independent polynomials.

Results-

The linear system gives 12 rotational invariant measures, after removing functionally dependent solutions. This is consistent with the dimension of the space of real symmetric SH: for L = 4, N = 15, and the 3 degrees of freedom of a rotation. These invariants are reported on Fig. 1.b. We fit an ADC profile in SH, L = 4, to images of a healthy subject scanned on a 3T Siemens 3T Verio TM scanner, 64 evenly distributed gradient directions, b = 3000 s/mm² and 2 non-weighted diffusion-weighted images. The maps reported on Fig. 1.b show a range of different contrasts, significantly different from FA and GFA, which suggests that they bring novel information on the "shape" of the ADC profile. Conclusions-We proposed a novel method to compute rotationally invariant measures in HARDI. The method is general and applies to the ADC, the ODF, or any other angular function represented by its SH coefficients. We believe that these measures will bring additional information on the angular complexity of diffusion functions, complementary to the indices already in use, based on tensor or on microstructure features. The results are complementary to recently published related methods 4,5 . The large number of measures will be used to create population classifiers, paving the way for pathology-specific biomarkers.

Figure 1 .

 1 Figure 1. a. FA and GFA, for reference. b. The 12 ADC-based polynomial invariants, for increasing truncation order, L, and increasing polynomial rank, t.


Emmanuel Caruyer

Luke Bloy

Birkan Tunç

Jérémy Lecoeur

Varsha Shankar

Ragini Verma

A comparative study of 16 tractography algorithms for the corticospinal tract: reproducibility and subject-specificity

Target audience-We present a comparative study of fiber tractography algorithms on a test-retest dataset. This work and the conclusions it leads to will guide the choice of a tracking algorithm for applications ranging from surgical planning to group-based connectivity analysis.

Purpose-We evaluate 16 different tracking algorithms, using their either freely available implementations or in house software. We compare their ability to produce reproducible tracking, as well as to capture subject-specific tract geometry.

Methods-Data acquisition. We scanned 9 healthy subjects (3F/6M, 31.25 ± 4.2y) at three time points separated by a 2 weeks interval. Diffusion-weighted EPI were acquired on a Siemens 3T Verio TM scanner, b = 1000 s/mm², on 64 evenly distributed directions plus 2 b = 0 images. Pre-processing. The diffusion-weighted images were regularized and corrected for Rician noise 1 . The diffusion tensor was estimated using linear least squares fitting, and fractional anisotropy (FA) maps were computed. Deformation field from and to the JHU Eve template 2 were computed using non-linear registration 3 of the FA images. Fiber tracking. The fiber tracking of left and right corticospinal tracks (CST) was performed in subject-space. When applicable, the tractography was initiated in the cerebral peduncle (CP), using 5 seeds per voxel in that region of interest (ROI). The step size was set to 0.2mm, and the maximum turning angle to 35°. Tracking was terminated whenever underlying FA fell below 0.15. The tract was a posteriori filtered to keep only streamlines going through the posterior limb of the internal capsule (PLIC), and decimated to keep only 300 streamlines per track. The regions of interest were taken from Eve labels, warped to subject space using the above-computed deformation field. Tracking algorithms are streamline and probabilistic streamline on constrained spherical deconvolution (MRtrix 5 ), FACT, second-order Runge-Kutta, streamline and tensorline on tensor (Diffusion toolkit [START_REF] Wang | ISMRM[END_REF] ), Gibbs globals tracking (MITK [START_REF] Reisert | [END_REF] ), 4th order Runge-Kutta and tensorline on tensor (Camino 8 ), tensorline, streamline, tensor deflection, bootstrap streamline, bootstrap tensor deflection on tensor and streamline on CSA ODF (in-house software), and unscented Kalman filter on multi-tensor model [START_REF] Malcom | MICCAI[END_REF] . Statistical analysis. Tract similarity was computed using a Gaussian process framework 4 . To this end, normalized tract probability images were computed in subject space, and warped back to template space. As a summary of the similarity matrix, we computed the ratio of the average inter-subject distance over the intra-subject distances.

Results-All the results reported here were obtained using a specific instance of the tracking parameters for each method. When applicable, the parameters were chosen to be as consistent as possible across the methods being compared. Reproducibility and specificity. All the methods tested reported an inter-subject average distance at least 3 times higher than its intra-subject counterpart (Fig. 1.a). Among the methods giving a lower score, we can identify two different trends. On one hand, some methods such as MITK Gibbs tracking, are more conservative, and therefore provide a very high intra-subject similarity, at the cost of loosing inter-subject discriminative power. On the other hand, some methods such as Camino Euler tracking are more exploratory, and on average captures differences between subject, at the cost of a lower reproducibility. Interestingly, the probabilistic streamline tracking in MRTrix is comparable with other methods in terms of reproducibility. Finally, the similarity matrix which presents the most well-defined clusters "visually" was obtained with DTK tensorline method. Conclusions-We compared a total of 16 tracking methods in terms of their ability to capture subject-specific white matter fiber path geometries. For the specific tracking parameters chosen for this study, we can identify 9 methods giving globally superior results (Fig 1 .a). A more comprehensive study, including possible variation of the parameters, and extension to a high angular resolution diffusion imaging dataset, will be undertaken in the future.

Figure 1 .

 1 Figure 1. a. Ratio of the average inter-subject distance over the average intra-subject distance. b. similarity matrix for selected methods; repetitions of the same subject are grouped. c. Result of left CST tracking for selected methods.