An improved fiber tracking algorithm based on fiber assignment using the continuous tracking algorithm and two-tensor model

Abstract

This study tested an improved fiber tracking algorithm, which was based on fiber assignment using a continuous tracking algorithm and a two-tensor model. Different models and tracking decisions were used by judging the type of estimation of each voxel. This method should solve the cross-track problem. This study included eight healthy subjects, two axonal injury patients and seven demyelinating disease patients. This new algorithm clearly exhibited a difference in nerve fiber direction between axonal injury and demyelinating disease patients and healthy control subjects. Compared with fiber assignment with a continuous tracking algorithm, our novel method can track more and longer nerve fibers, and also can solve the fiber crossing problem. Research Highlights This study tested an improved fiber tracking algorithm based on fiber assignment using a continuous tracking algorithm and a two-tensor model. In this algorithm, tracking decisions are made by judging the estimation type of each voxel. This method can solve the cross-track problem. Abbreviations ADC, apparent diffusion coefficient; MD, mean diffusivity; FA, fractional anisotropy; RA, relative anisotropy

Introduction

Diffusion tensor imaging is an important complementary technique used in conjunction with diffusion weighted imaging. Diffusion weighted imaging utilizes the diffusion characteristics of water in tissue, in accord with water diffusion anisotropy. Diffusion tensor imaging can be used to detect the microstructure of tissue, including muscle, bone and white matter in the brain [1],[2] . It will be an important noninvasive way to explore the microstructure of the human brain and of high potential value in the diagnosis and treatment of nervous system diseases including multiple sclerosis, epilepsy, Alzheimer′s disease, and schizophrenia [3],[4],[5] etc. Diffusion tensor imaging has also been applied to preoperative neurosurgical visualization and the surgical planning of brain tumor surgery [6],[7],[8] .

Many algorithms have been used to perform fiber tracking, including fiber assignment with a continuous tracking algorithm [9] , the Runge-Kutta method [10] , and the tensor deflection algorithm [11] . The fiber assignment by continuous tracking algorithm and the Runge-Kutta method are versions of the streamline tracking algorithm, which takes the voxel principal direction as the fiber tracking direction. Tensor deflection uses the whole diffusion tensor matrix to determine the tracking direction. Fiber assignment by continuous tracking is currently the most commonly used method in clinical settings. However, the fiber assignment by continuous tracking algorithm cannot resolve the fiber crossing issue, which results in much shorter or sparser fibers. Several novel methods have recently been proposed. For example, Qazi et al [12] used a two-tensor extended streamline tractography technique to track the corticospinal tract with crossing fibers. Landman et al [13] used a Crossing Fiber Angular Regulation of Intra-voxel Structure approach to resolve fiber crossing, and set a robust framework. Each method has advantages and disadvantages, and there is currently no gold standard. The present study attempts to deal with this problem by combining the two-tensor model with a common algorithm.

Results

Quantitative analysis of subjects

A total of eight healthy subjects, two patients with axonal injury and seven patients with demyelination were included in the final analysis.

Basic parameter maps

In areas of large fiber bundles, diffusion anisotropy is strong, and the principal diffusion directions align well with the long axis of the bundle. This means that most white matter fibers can be reconstructed successfully based on vector fields. Here, we compare the vector field map of a normal subject′s brain with the brain vector field map of an axonal injury patient. The results revealed that the vector directions of the healthy brain were much more symmetrical compared to the axonal injury patient′s brain ([Figure 1]). {Figure 1}

As such, vector field maps can also be used to show the abnormal fiber directions indirectly.

Apparent diffusion coefficient (ADC), fractional anisotropy (FA) and relative anisotropy (RA) are three basic parameters that can be used in diffusion tensor imaging data processing. Maps of these parameters are shown in [Figure 2], calculated using custom-built "diffusion tensor imagingprocessTool" software. Due to skull removal, our maps are not affected by background noise. FA and RA maps show the contour of white matter much more clearly than ADC maps. {Figure 2}

Fiber tracking results

[Figure 3] shows fiber tracking results from a healthy subject and an axon injury patient in the midbrain area. [Figure 3]A0 is part of a fluid-attenuated inversion recovery (FLAIR) image from a healthy subject, while [Figure 3]A1-A3 shows data from the axonal injury patient with different rectangular regions of interest, and [Figure 3]B0-D3 are the fiber tracking results. The parameters used during the tracking process strongly influence the results [14],[15] .{Figure 3}

The parameters determine when fiber tracking should terminate (i.e. the stopping criteria). Therefore, we set three important parameters, the FA value threshold, the deflection angle threshold, the iteration step to the experimental values 0.2, 50° and 0.5, respectively, for both the traditional method and our own algorithm, before comparison. The fiber assignment by continuous tracking algorithm results are shown in yellow ([Figure 3]B0-B3) for comparison with our new algorithm results in red ([Figure 3]C0-C3). The results revealed that our method can track more fibers and longer fibers than the traditional method. From a clinical point of view, the axonal injury patient ([Figure 3]A1-D1) exhibited fewer fibers than the healthy control ([Figure 3]A0-D0). Substantially greater differences were revealed in the comparison between the normal side ([Figure 3]A2-D2) and abnormal side ([Figure 3]A3-D3) of the midbrain area in the same axonal injury patient.

The comparison shown in [Figure 4], between a healthy control and a demyelinating disease patient, was produced using a similar method to [Figure 3]. Our method ([Figure 4]C0-C1) was found to be more effective for fiber tracking compared to the traditional method ([Figure 4]B0-B1). Moreover, the demyelinating disease patient exhibited fewer fibers compared to the healthy control in the internal capsule area ([Figure 4]A1-D1), consistent with clinical data.{Figure 4}

Discussion

The currently proposed algorithm for detecting crossing tracts in the white matter of the human brain is based on the mixed tensor model, in accord with physiologically realistic features in a brain connectivity network. Information transfer between processing units is the basis of human brain function. White matter is the key medium determining information transfer, constituting a network of microscopic cellular wire-like axons, which carry electric signals from their own cell bodies to the synapses of other cells. Diffusion tensor tractography enables the visualization of this network, providing an effective tool for exploring abnormal functional connectivity in the brain. Many diseases are related to white matter abnormalities, such as diffuse axonal injury and demyelinating disease.

The present results revealed that the tracking fibers of patients with demyelinating disease or diffuse axonal injury are very different from those of a healthy control, consistent with the disease physiology.

Regarding fiber tracking methodology, our improved method (fiber assignment by continuous tracking algorithm based on mixed tensor models) is superior to the traditional algorithm in the following ways. First, because the method involves the two-tensor model, the direction of a voxel theoretically does not only depend upon its main direction, but also on the type of estimation. For line estimation, there is only one direction; while for the plane estimation, two directions are tracked. Although the processes underlying estimation and decision-making are complicated, this method can be used to solve cross-track problem, avoiding the missing or shortening of fibers. Second, in practice, the tracking results in our experimental data are consistent with theory. Third, when we change the tracking parameter threshold within acceptable values, such as tracking step (0 < step < 0.8), FA value (0.15 < FA value < 0.25) and deflection angle (45° < angle < 70°), our improved algorithm always produced better results than the traditional algorithm. However, this algorithm has its own limitation. The improved algorithm is time consuming, taking approximately twice as long as the traditional algorithm to finish the tracking, due to increased complexity. In future studies, we will perform algorithm optimization, such as reducing loop times and rewriting "if" sentences. Since diffusion tensor imaging is the only approach to noninvasively detect the architecture of white matter and the performance of the particular tracking method used has a strong effect on the tracking results, our improved algorithm is potentially valuable for improving the analysis of nervous system diseases.

Subjects and Methods

Design

Molecular imaging algorithm study.

Time and setting

This study was performed at the Department of Radiology, Xiamen Second Hospital, China from February to July 2011.

Subjects

A total of eight healthy subjects (average age, 54.50 years), two axonal injury patients and seven patients with a demyelinating disease (average age, 55.78 years) were included. Healthy subjects were volunteers from the Department of Radiology, Xiamen Second Hospital, China (without nervous system disease). Axonal injury and demyelination patients were recruited from the out-patient clinic or inpatients at the Xiamen Second Hospital. They were diagnosed by two attending physicians from the Department of Image and an experienced doctor in accord with clinical symptoms, magnetic resonance imaging and biochemistry tests.

Subjects gave written informed consent before participating. The examination took 15 minutes, and did not cause any harm to participants, in accord with the Helsinki Declaration.

Methods

Data acquisition

This data set was acquired using dual spin-echo DW-EPI with 13 gradient directions and 28 slices (whole-brain coverage) using an 8NVHEAD coil and a GE SIGNA HDe 1.5 T scanner (GE Healthcare, Bethesda, MD, USA). The acquisition parameters were as follows: b = 1 000 s/mm 2 , flip angle = 90°, echo time = 107 ms, repetition time = 10 000 ms, slice thickness = 4.0 mm, spacing = 1.0 mm, field of view = 26 cm × 26 cm, matrix = 128 × 128, number of excitations = 2.0. The time of diffusion tensor imaging data acquisition was approximately 5 minutes. T2-weighted images were acquired, and we used the FLAIR sequence shown in [Figure 5]. Parameters: repetition time = 8 602 ms, echo time = 123.8 ms, inversion time = 2 100 ms, slice thickness = 4.0 mm, spacing = 1.0 mm, field of view = 24 cm × 24 cm, matrix = 288 × 160, number of excitations = 1.0. The procedure lasted approximately 3 minutes. {Figure 5}

Data preprocess

Preprocessing raw data were necessary to remove effects caused by movement or noise after acquisition. There are several traditional preprocessing steps using statistical parametric mapping software, as follows: realign-slice timing-normalizing-smoothing. In performing data preprocessing, it is important to maintain the originality of the data. However, it is unclear to what degree our preprocessing method is able to achieve this. Statistical parametric mapping preprocessing might be an appropriate method, but is time-consuming. As such, we chose a simpler "fast preprocessing" method to keep maximum data originality. Our fast preprocessing method is described as follows. First, transfer the DICOM data into an analyzed data format (".img, .hdr" in MRIcro software; www.mricro.com). Second, smoothing is performed to reduce the data noise. Third, skull removal is used to remove background noise, to avoid analyzing data outside the brain when performing subsequent fiber tracking. [Figure 6] shows data before and after preprocessing. Compared with [Figure 6]A, [Figure 6]B exhibits less image noise.{Figure 6}

Diffusion tensor imaging principles

The motion of water molecules is either free or restricted depending on tissue structure. Water diffusion happens in three dimensions. If the tissue contains many fibers, the velocity of water diffusion is fastest along the direction of the fiber, and slowest in the direction perpendicular to it. This characteristic of diffusion is called anisotropy [16] . While fibers are fewer, water diffuses isotropically. Due to the above characteristics, water molecules are diffused in an anisotropic way in the white matter of the brain, which contains abundant fibers, under most circumstances. In 1994, Basser et al [17] proposed the concept of diffusion tensor, assuming the shape of water diffusion as an ellipsoid. In geometry, the ellipsoid is expressed as a three-dimensional matrix; the diffusion tensor matrix D. The signal S and D can be defined as follows,

Where S and S 0 are the signal intensity and the original signal intensity without diffusion weighting, respectively, while b is a factor describing the gradient timing and strength, and g is a unit vector representing the direction of a diffusion gradient.

[INLINE:2]

The three axes of the ellipsoid correspond directly to the three descending eigenvalues: l1 , l2 , l3 , and their corresponding eigenvectors are e1 , e2 and e3 . Some commonly used quantitative indicators, such as mean diffusivity (MD), FA and RA, can be derived from eigenvalues. These parameters, which have been used to evaluate white matter disease in multiple sclerosis, Alzheimer′s disease, epilepsy and brain tumors, are important in clinical application [18],[19] .

[INLINE:3]

Meanwhile, we can perform fiber tracking according to eigenvectors [20] . We express the fiber assignment by continuous tracking algorithm and our own algorithm below.

Fiber assignment by continuous tracking algorithm

In each voxel, the fiber assignment by continuous tracking algorithm takes the corresponding eigenvector of the biggest eigenvalue as the fiber tracking direction. This is complimented in a continuous coordinate space. According to this rule, given a start point′s coordinate is xt , its corresponding eigenvector is e(→)1t , Δ t is the step, then the next point′s coordinate x t+1 can be expressed as,

[INLINE:4]

All the points tracked into a line are then connected, forming a fiber according to the start point. Accordingly, we tracked all fibers in a given region of interest. The criteria during the tracking process should be noted here.

First, the FA value threshold setting was typically set to 0.2, which was regarded as the boundary value between white matter and gray matter, because the FA values of gray matter were usually lower than 0.2 [21],[22] .

Second, the deflection angle threshold setting was typically set to 45°-70°, because the probability of sharp turns between two adjacent neurons is close to zero in the white matter of human brain [15] .

This method, and the computation involved, was simple and easy to understand. This computation can effectively perform fiber visualization in the region with significant diffusion anisotropy. Nevertheless, when a region with fiber crossing is encountered, the tracking is partially finished or interrupted, resulting in shorter or sparser fibers.

Two-tensor model

The similarity of the diffusion ellipsoid was separately categorized into line model, plane model or sphere models. In general, most voxels contained one main fiber direction, which is the line model. However, more than one main fiber direction was exhibited when there were crossing fibers. Using the single-tensor model (Equation 1) to fit this type of diffusion data leads to errors. When two crossing fibers (plane model) were taken into account, the two-tensor model [23] was defined as follows,

[INLINE:5]

Where, s0 and l3 were calculated from the single-tensor fitting, f was the fraction of the first tract, qi was the angle between two principal directions in the plane and l1 was the principal diffusivity. Here, we could solve Equation 7 using Levenberg-Marquet method in MATLAB.

Based on the two-tensor model, this method significantly improved the ability to track smaller and longer fibers in the brain. Each voxel was described using a two-tensor model and performed fiber tracking, but this would require a lot of time and a high level of computer performance. In addition, most of the voxels contained one main fiber direction. Due to reasons above, we combined the single-tensor model, two-tensor model and fiber assignment by continuous tracking algorithm to perform fiber tracking.

Fiber assignment by continuous tracking algorithm based on mixed tensor models

A voxel was estimated into the line (λ1 > > λ2≈λ3 ), plane (λ1≈λ2 > > λ3 ), or sphere (λ1≈λ2≈λ3 ) model according to its eigenvalues. Our new algorithm steps:

First, calculate all voxels′ tensors using single-tensor model.

Second, judge the type of estimation of each voxel, and divide them into three groups:

(1) Line estimation: using single tensor fitting.

(2) Plane estimation: using two-tensor fitting.

(3) Sphere estimation: no fitting, tracking is ended here.

After the regulation described above, we obtained tensors for all voxels. We performed fiber tracking using our fiber assignment by continuous tracking algorithm based on mixed tensor models from a starting point. Departing from the traditional fiber assignment by continuous tracking algorithm, we performed estimation-type judgments for each voxel before tracking the next point. Thus, for line estimation, we tracked the next point according to this voxel′s principal direction; when performing plane estimation, we tracked the next two points according to this voxel′s two principal directions. This method can solve the cross-track problem and improve the accuracy of fiber tracking. [Figure 7] shows the algorithm flowchart.{Figure 7}

All algorithms are coded with MATLAB (Version R2009b; MathWorks, Natick, Massachusetts, USA) environment, and the software named "diffusion tensor imaging processing tool" based on the Graphical User Interface, developed to post-process data efficiently.