UNC-Utah NA-MIC framework for DTI fiber tract analysis

Diffusion tensor imaging has become an important modality in the field of neuroimaging to capture changes in micro-organization and to assess white matter integrity or development. While there exists a number of tractography toolsets, these usually lack tools for preprocessing or to analyze diffusion properties along the fiber tracts. Currently, the field is in critical need of a coherent end-to-end toolset for performing an along-fiber tract analysis, accessible to non-technical neuroimaging researchers. The UNC-Utah NA-MIC DTI framework represents a coherent, open source, end-to-end toolset for atlas fiber tract based DTI analysis encompassing DICOM data conversion, quality control, atlas building, fiber tractography, fiber parameterization, and statistical analysis of diffusion properties. Most steps utilize graphical user interfaces (GUI) to simplify interaction and provide an extensive DTI analysis framework for non-technical researchers/investigators. We illustrate the use of our framework on a small sample, cross sectional neuroimaging study of eight healthy 1-year-old children from the Infant Brain Imaging Study (IBIS) Network. In this limited test study, we illustrate the power of our method by quantifying the diffusion properties at 1 year of age on the genu and splenium fiber tracts.


INTRODUCTION
Since its invention in the 1980s diffusion weighted imaging has become increasingly popular for the analysis of brain pathologies and development. It is now standard practice for investigators to collect diffusion weighted images (DWI) concurrently with other modalities of interest without any expertise on how to preprocess or analyze DWI. To address this growing wealth of DWI data, we have created a coherent framework of tools for the pre-processing and analysis of DWI and diffusion tensor images (DTI) that is accessible for the non-technical user. It is the aim of this paper to describe our DTI processing workflow and to provide example data processed with this framework for increased workflow clarity.
Diffusion weighted images capture the movement of water molecules. When locally unconstrained, such water molecules move by Brownian motion in an isotropic spherical fashion, for example as is the case within the brain lateral ventricles. However when locally confined, the diffusion of water molecules is constrained by local boundaries, such as cell membranes. Within the human brain white matter (WM), largely composed of neuronal axons, the diffusion of water is more ellipsoid or anisotropic in shape due to the confines of the internal microstructure and the myelin surrounding neuron axons.
Characterizing the properties of diffusion can inform us about the integrity of WM micro-structure and provides insight on the mechanisms and progression of disease. To date DTI has been used in research to further the understanding of pathology in many neurologic diseases including multiple sclerosis (Zhu et al., 2013), Alzheimer's (Patil et al., 2013), and Parkinson's (Du et al., 2012), as well as normal neurodevelopment (Geng et al., 2012), and aging (Lebel et al., 2012). The most common properties measured in DWI/DTI are axial (AD, λ || ), radial (RD, λ ⊥ ), and mean diffusivity (MD), as well as fractional anisotropy (FA). Axial diffusivity has been shown to decrease in the case of cell death, while radial diffusivity will increase in the case of myelin damage (Song et al., 2003). Alterations in FA reflect changes of AD, RD, and MD and describe the overall pattern of diffusion.
Several analysis frameworks exist for the study of DWI/DTI data, typically falling into one of three categories: region of interest (ROI) based analysis, voxel based analysis (VBA), and quantitative tractography. ROI driven analysis can be performed by registering individual subject DWI images either to their corresponding anatomical images or by registering the DWI of all individuals to a prior existing atlas. In the first type of ROI analysis, regional segmentations from structural (T1w, T2w) MRIs of the same subject are propagated via linear or deformable co-registration Frontiers in Neuroinformatics www.frontiersin.org of the structural images to the diffusion image data and mean statistics are collected (Shi et al., 2013). However, the structural WM parcellations are commonly lobar based and thus cannot separate the measurements from different fiber tracts located in the same cortical lobe/region. In the second type of ROI analysis coregistration with a prior WM fiber atlas is performed (Faria et al., 2010). This type of analysis again requires deformable registration and results in a single measure per fiber region, potentially combining regions that show significantly different fiber situations. The positive aspects of both of these types of region based analysis is that processing is relatively simple and it is robust against imperfect registration. On the other hand, the analysis of large ROIs means non-specific findings as ROIs can contain multiple fiber tracts and multiple fiber situations (single, fanning, crossing fiber situations). This form of analysis thus leads to limited localization of findings. Another DTI processing framework is VBA. This style of analysis requires registration of each subject into a study specific reference space and then all voxels are tested individually for significance (Liu et al., 2009). While this technique allows for whole brain analysis and is suitable for hypothesis generation, this analysis assumes absolutely perfect registration on a voxel by voxel basis and is thus highly sensitive to the deformable co-registration. Another form of VBA is FSL's Tract Based Spatial Statistics (TBSS). Here a WM skeleton is created for the brain, maximal FA values are projected to the WM skeleton, and voxel wise analysis of the skeleton is performed . Due to the projection to the WM skeleton, TBSS is less sensitive to registration accuracy compared to standard VBA. While FSL's TBSS tool provides a coherent framework that works well with imperfect registration, it does not have an explicit tract representation and instead provides a skeletal voxel representation that cannot be uniquely linked to individual fibers throughout the brain. Further, the use of maximal FA values renders TBSS sensitive to DWI artifacts.
Our framework falls into the category of quantitative tractography, where anatomically informed curvilinear regions are used to analyze diffusion at specific locations all along fiber tracts. Diffusion property profiles are computed as weighted averages taking into account the neighborhood at each step along the fiber bundles . This form of analysis results in highly localized statistics that can be visualized back on the individual fiber bundles.
In this paper, we propose the use of our 3D Slicer based UNC-Utah NA-MIC DTI framework for a coherent atlas fiber tract based DTI analysis. Most steps of the framework utilize graphical user interfaces (GUI) to simplify interaction and provide accessibility for non-technical researchers. At the heart of the framework is a set of tools anchored around the multi-purpose image analysis platform 3D-Slicer. Several workflow steps are handled via external modules called from Slicer in order to provide an integrated approach. Our workflow starts with data conversion from DICOM format, followed by an automatic as well as interactive DWI and DTI quality control (QC) steps. Though underappreciated in the field, appropriate and thorough QC is a must for every DTI study. Our framework then centers around a DTI atlas that is either provided as a prior template or computed directly as an unbiased average atlas from the study data via deformable atlas building . Fiber tracts are defined via interactive label map tractography on that atlas followed by consecutive interactive fiber cleaning steps. Individual subject DTI fiber tract profiles of FA, MD, RD, and AD are extracted automatically using the atlas mapping information. These fiber tract profiles are then analyzed using our statistics toolbox (FADTTS; Zhu et al., 2011). The local statistical results are then mapped back on to the fiber bundles for visualization within 3D Slicer (Pieper et al., 2004).
This framework has been developed over several years and we continuously improve the current toolset. Within the past year we have added entropy-based detection of vibration artifacts as part of the automatic QC and an interactive atlas building tool called DTIAtlasBuilder. All tools that are part of the UNC-Utah NA-MIC DTI framework are open source and available on NITRC, providing the neuroimaging field with a transparent and coherent toolset for DTI fiber tract based analysis.
We illustrate the use of our framework on a small sample, cross sectional neuroimaging study of eight healthy 1-year-old children from the Infant Brain Imaging Study (IBIS) Network (http://www.ibisnetwork.org/UNC). In this paper, we demonstrate our method by quantifying the diffusion properties within the genu and splenium fiber tracts. The corresponding raw data, intermediate files, and final results for data presented in this manuscript have been made freely available on NITRC.

MATERIALS AND METHODS
Our framework for the fiber tract based analysis of DTI images is composed of four essential segments: (1) QC, (2) Atlas Creation, (3) Interactive Tractography, and (4) Statistical Analysis. The tool workflow can be seen in Figure 1 and an overall conceptual schematic of the process is displayed in Figure 2. All of the tools referenced in the description of our workflow can be utilized as stand-alone command line applications to facilitate scripting and grid computing, or interactively as part of 3D Slicer as external modules. While the first segment of our framework is termed QC, it is critical to perform quality assessment of the processing at every step in the framework to ensure the data that is analyzed is correct.

DWI quality control -DTIPrep
Analysis starts by converting the raw DWI from DICOM to NRRD format using a tool called DicomToNrrdConverter in 3D Slicer. This tool is capable of computing DWI via the b-matrix information stored in Siemens DICOM headers, as well as recognizing many other manufacturer specific tags. A tool called DTIPrep (Liu et al., 2010) is then used to check for and correct several common artifacts found in DWI. In summary, DTIPrep performs image information checking, slice-wise artifact detection, creates a baseline average, corrects for subject motion and eddy current induced deformations, and finally detects dominant direction artifacts . To carry out the automatic QC a protocol must be created that contains all parameters of the DTIPrep QC. Protocols are fully modifiable for individual studies and include a significant number of parameters for every step. DTIPrep gives the option of creating a default protocol from a given scan. On completion of the automatic QC it is further Frontiers in Neuroinformatics www.frontiersin.org
Step 2 is atlas creation.
Step 4 is parameter profile creation and statistical analysis.
advised to perform additional visual QC. Here artifacts that may be missed by the automatic QC step can be detected and gradients still containing artifacts or motion can be excluded from the final DWI ( Figure 3A). After visual QC, the DTI tensor and property images are then automatically estimated from the DWI data via the dtiestim and dtiprocess modules.

DTI quality control -3D slicer
In addition to checking the DWI's for artifacts, the resulting DTI must also undergo visual QC using 3D Slicer. This QC step is to ensure the correct directional encoding, sufficient signal-to-noise (SNR) in the FA data, and to test the overall quality of the data by performing landmark-based/"fiducial" tractography. Visual DTI QC begins by loading the FA image and scrolling through the slices for each plane (axial, sagittal, coronal) to assess the apparent level of SNR for the volume. Next, the DTI volume is analyzed to determine if the header information encoding directionality for the volume is correct. The DTI image is colored by tract direction with tracts moving right to left colored red, tracts moving anterior to posterior colored green, and tracts moving inferior to superior colored blue. The corpus callosum (CC) should appear red, the cingulum green, and the corticospinal tract should be blue. Also special attention is directed to detect dominant directional artifacts if low FA regions (such as the gray matter) display a single color throughout most of the DTI image, as seen in Figure 3B. Visualization of the primary direction via line glyphs is the next step to ensure the correctness of the diffusion measurement frame. Glyphs should follow the direction of the main fiber tracts. If the glyphs are incorrect (due to incorrect DICOM information or scanner software malfunction), then the corrected measurement frame must be determined and edited in the volume header.
Once the header information is confirmed to be correct for the volume, fiducial based tractography will provide information whether tractography is feasible, and anatomically accurate. Fiducial tractography is commonly performed for the major fiber tracts (CC genu and splenium, cingulum, and corticospinal tract) as well as any tracts of particular interest for the neuroimaging study at hand. In the resulting fiber tracts, it is important to observe (a) sufficient quantity of fibers delineating a fiber tract, (b) whether fibers completely track between the expected terminal regions, and whether the fiber location is anatomically accurate throughout the length of the fiber tract. If artifacts remain in the volume, then fibers may break prematurely, appear warped, or incorrectly localized, or if the data has a very low SNR the fiber bundles of interest may be impossible to track. If tracking is not feasible or anatomically incorrect, then exclusion of this dataset is suggested.

Brain mask creation, editing, and application -bet2/AutoSeg
The first step in atlas building is to extract out the brain regions for each DTI image. To do this, a brain mask must be computed, edited, and applied to the DWIs in order to create a skull-stripped DTI. Many tools are available for brain mask creation, such as the "Brain Mask from DWI"module in 3D Slicer, bet2 within DTIPrep, or a tissue segmentation based approach in AutoSeg (Van Ginneken et al., 2007). Brain mask creation by bet2 can be computed Frontiers in Neuroinformatics www.frontiersin.org FIGURE 2 | Schematic of DTI processing workflow. The input raw DWI data first undergoes vigorous quality control before the final conversion to DTI for Atlas building. Interactive tractography is performed on the Atlas.
From these fiber bundles, diffusion parameter profiles are extracted and analyzed along with other study variables to determine significant interactions. using the baseline image or the IDWI (isotropic DWI), while both are used jointly to compute the brain mask with AutoSeg. Mask editing is performed also within 3D Slicer, but any editing tools with Slicer compatible format (NRRD, Nifti, GIPL,Analyze) can be employed such as itk-SNAP (Yushkevich et al., 2006). The edited brain mask can then be applied to the DWIs, creating a skull stripped DTI (referred to as "Original DTI" in this paper) ready for atlas building.

Atlas creation -DTIAtlasBuilder
The main atlas building step is performed with a tool called DTI-AtlasBuilder, detailed in Figure 4. This program allows the user to create an atlas image as an unbiased average of all co-registered DTI images in the study. Atlases are created iteratively starting from affine atlases (step 1) and moving to deformable diffeomorphic atlases (step 2, 3, and 4). The registration is done in three steps: (1) affine registration with the Slicer General Registration module (called BRAINSFit; Johnson et al., 2007), (2) unbiased diffeomorphic atlas computation (Joshi et al., 2004) with the GreedyAtlas module within AtlasWerks (Joshi et al., 2004), and (3) a refinement step via symmetric diffeomorphic registration with the Slicer DTI registration module DTI-Reg using the Advanced Normalization Tools (ANTS; Avants et al., 2008;Klein et al., 2009). The final, refinement step significantly improves registration accuracy in our experience, most likely attributable to the use of local normalized cross-correlation as the image similarity metric. A final step will combine multiple transforms to compute the overall transformation (GlobalDisplacementFields) from the original DTI images to the Final DTI Atlas. This final transform is required for the analysis of fiber tract properties in the individual subject space.
In the first step, DTIAtlas Builder will generate scalar FA images from the skull stripped DTIs. The affine atlas building starts by registering subject images to a pre-exiting FA template in MNInormative space or to the first subject in case no prior template is specified. In order to avoid biasing the subsequent analysis, the image intensities are calibrated with the template's intensities using standard histogram quantile normalization available in ITK. The affine transforms are applied to the original DTIs and the transformed images are averaged to create the Affine DTI Atlas (see Figure 4). Each image's affine transform is then updated by recomputing the affine registration in a second loop to the Affine DTI Atlas (rather than to the initial template as in the first loop). This transform, Affine Transform 1B, is then applied to the original DTI to create affinely co-registered DTI's for each subject.
In the second step of atlas building, the unbiased diffeomorphic atlas is computed from the affinely registered FAs via diffeomorphic, fluid flow based atlas building. This results in deformation fields from the affine space to the unbiased atlas space. These deformation fields are then applied to the affinely registered DTIs and the transformed images are averaged to create a Greedy Diffeomorphic DTI Atlas (Figure 4).
The last stage of atlas building utilizes DTI-Reg to refine the diffeomorphic registration computed in the previous step.

Frontiers in Neuroinformatics www.frontiersin.org
To start, FA images from the original subject DTIs are histogram normalized to the Greedy Diffeomorphic FA Atlas, and then diffeomorphically registered to it. The resulting deformation fields are then applied to the original subject DTI and the resulting diffeomorphically registered DTIs are averaged to create the Final DTI Atlas. The deformation fields are then updated with another loop of registrations with DTI-Reg in order to coregister the images with the Final DTI Atlas (rather than with the Greedy Diffeomorphic FA Atlas as in the first loop). The overall Global Displacement Fields are finally computed by composing the affine transform from step 1, as well as the affine and deformation transform from this last step. These displacement fields will be used to map the fiber tracts created in atlas space to the individual subject spaces for analysis. The final atlas co-registered DTI images are computed by applying these final deformations fields. It is noteworthy that this deformation of the DTI images is unlike the processing in TBSS/FSL, where FA images are deformed. Furthermore, all DTI image interpolations are performed using standard Log Euclidean based resampling (Arsigny et al., 2006). With the atlas building complete, it is important to perform QC on the Affine, Greedy Diffeomorphic, and Final DTI Atlases, simply accessible with a single button within DTI-AtlasBuilder. In this QC step we can visually determine the point to point correspondence between the subject images and the atlas to determine if computed registration transforms are appropriate. The progression of atlas sharpness can be seen in Figure 4C.
Software path configuration for the tool can be performed automatically or manually, and the settings can be saved for future use. This tool also has a "no GUI"/command line option, given that the data file, parameter file, and software configuration file are available (all can be determined and saved via the GUI). There is also direct support for the GRID software LSF (Load Sharing Facility, Platform Computing, Inc).

Label map single tensor tractography -3D slicer
While atlas building is fully automatic, the definition of fiber tracts of interest needs to be performed interactively. Fiber tractography only needs to be performed once on the Final DTI Atlas built in the prior step. In our framework, fiber tractography is performed via the label map tractography module or the full brain tractography module within 3D Slicer. If using label map tractography, each fiber tract has a separate anatomically informed label created in the editor module and is tracked individually (Figure 5). Alternative fiber tracking tools can be employed as long as the resulting fibers are stored in standardized vtk-based file format. Resulting DTI fiber bundles require post-processing via ROIs to remove unwanted, erroneous fibers.

Coarse fiber bundle editing with ROIs -3D slicer
The results of label map tractography may contain fibers in addition to the fiber bundle of interest. Coarse cleaning is performed with interactive 3D regions of interest (ROIs) in 3D Slicer's tractography display module. Here fibers to be kept or removed are selected using ROIs iteratively. ROIs should also be utilized for separating tracts that run within the left and right hemispheres in order to allow separate analysis of the right and left portions (for example the left and right hemispheric portions of the fornix).

Fine fiber bundle editing -FiberViewerLight
Fiber tracts often need additional, finer editing via clustering and cleaning within FiberViewerLight. Fiber length based cleaning and several fiber clustering algorithms (center of gravity, Hausdorff distance, mean distance, normalized cut) to further delineate which fibers belong in the bundle of interest and which ones do not (Figure 6). The decision for inclusion or exclusion is based on expected fiber tract anatomy. Once fiber bundles are deemed appropriate, diffusion property profiles for each fiber tract and subject image are extracted as described in the next section.

Diffusion property profiles -DTIAtlasFiberAnalyzer
As a final step prior to statistical analysis, diffusion properties along the fiber tracts, called fiber property profiles, are extracted. For that purpose, a fiber tract parameterization needs to be established, which is achieved via arclength parameterization (Corouge et al., 2005) starting from each fiber's intersection with an "origin" plane. While DTIAtlasFiberAnalyzer automatically determines an appropriate origin plane, investigators can choose to interactively determine that plane in order to associate a specific anatomical fiber location to this origin within FiberViewerLight (see Figure 6B).
The parameter profile generation and processing is performed with a tool called DTIAtlasFiberAnalyzer, shown in Figure 7. The profile extraction operates on the original DTI data in order to avoid the need for tensor deformation by deforming the parameterized fibers into original space via the previously computed Global Deformation Field. At the deformed fiber locations the tool can then calculate the different fiber property profiles along each fiber in the bundle (FA, RD, AD, MD, GA, Fro, lambda 2, and lambda 3) in subject space. Figure 7B shows example FA profiles of the genu calculated and visualized with DTIAtlasFiberAnalyzer for this study. Visualization of these Frontiers in Neuroinformatics www.frontiersin.org fiber property profiles for each subject versus the fiber bundle in atlas space allows an alternative assessment of the registration quality to atlas space ( Figure 6B). Inspecting the atlas fiber bundle parameterization is an additional QC step to confirm correct parameterization ( Figure 7C). An error in parameterization usually requires a change in origin plane to correct the problem.
The fiber property profiles are the input for the analysis tool called Functional Analysis of Diffusion Tensor Tract Statistics (FADTTS; Zhu et al., 2011) that performs statistical analysis on the fiber bundle profiles. This tool has the ability to compute group difference and correlational analyses (Figure 8) used to analyze the data for this paper are available with the sample data at www.nitrc.org.

Statistical analysis of diffusion properties along the fiber tract -FADTTS
Using the diffusion profiles obtained from DTIAtlasFiberAnalyzer, it is possible to test how these diffusion properties relate to other variables of interest within the sample data. FADTTS performs several statistics tests including a multi-variate coefficient model, weighted least squares estimation, and functional principal component analysis. These tests produce both global and local test statistics with confidence bands and the local values are corrected for multiple comparisons with false discovery rate (FDR). Several statistical plots are commonly created, a sample of which can be seen in Figure 8. From the beta plots, the direction, magnitude, and significance of the relationship of interest can be assessed by the beta sign, distance from zero, and the significance marker (a filled circle) on the curve, respectively (Figures 8C,D).

Mapping covariate p-values on fiber bundles for visualization -MergeStatWithFiber
Upon computation of test statistics with FADTTS, the statistics can be mapped to their corresponding location on the fiber bundle using a tool called MergeStatWithFiber. The input for MergeS-tatWithFiber is the comma-separate-value (csv) file resulting from FADTTS with the addition of arclength as the first column, and a header row containing column labels. The result of this merge is then visualized in 3D Slicer ( Figure 8E). The color map shown is determined by the significance threshold set in MergeStatWith-Fiber, as such that all points on the fiber that have a p-value above the significance threshold are colored blue, and all points with values below the significance threshold are colored from cyan to red, with points indicated in red signifying the areas of greatest significance. If desired, this color map can be altered within Slicer to another color table. MergeStatWithFiber can be used as a command line tool or through the final tab in the DTIAtalsFiberAnalyzer GUI (Figure 7).

RESULTS
For our example data we created an unbiased diffeomorphic atlas with eight 12-month-old IBIS participant scans, and performed interactive tractography of the genu and splenium on this atlas. Diffusion parameter profiles were created for these two fiber bundles and these profiles were further analyzed with FADTTS. While each of the participants were categorized as 12-month-old for the IBIS study, the interaction between a subject's actual age in weeks (ranging from 11.8 to 14.5 weeks) and the DTI properties was analyzed to illustrate the use of our tools. On the genu, three areas showed significant interaction between Age and FA, but no areas of significant interaction were found on the splenium (data not shown). Figure 8C shows the graphical output of this Frontiers in Neuroinformatics www.frontiersin.org result from FADTTS in the blue line for Age with the filled in circles indicating significance. Figure 8E portrays the same information of significance mapped back on to the original genu fiber bundle.

DISCUSSION AND CONCLUSION
With the use of diffusion tensor imaging (DTI) on the rise, the field of neuroimaging is in need of a coherent paradigm for localized fiber tract based DTI analysis. Here, we illustrate the UNC-Utah NA-MIC framework for DTI fiber tract analysis assessing the genu and splenium in 1-year-old from the UNC arm of the IBIS study.
To demonstrate the use of our framework we analyzed the relationship of subject age in weeks with diffusion tensor properties. For the genu there were three locations of significant interaction as visible in Figure 8. The splenium demonstrated no areas of significant interaction. To facilitate utilization of our workflow, all of the data, intermediate files, scripts, and tools are freely available on NITRC.
Over the past few years we have created and tested our tools to provide maximum usability and a detailed analysis of diffusion at each point along fiber bundles. As these are tools for research, the tools in our framework are constantly evolving as we improve their function and incorporate new options. For example, we are currently enhancing DTIPrep to allow for simulation based error estimation, enhancing DTIAtlasBuilder to improve its unbiased atlas building step, as well as providing a GUI (non-Matlab) based interface to FADTTS directly within DTIAtlasFiberAnalyzer.
While all tools employed here are individually available online (see Software and example dataset below), all needed modules of this framework are directly available for download via the 3D Slicer extension manager (version 4.3 and above; Figure 9). We hope our workflow will fill the void in the field for the many neuroimaging investigators that need to analyze their DTI data with a localized fiber based approach.