Brain Structure MR Imaging Methods: Morphometry and Tractography

1. Introduction

Brain morphology is in constant change from the very beginning of the neurodevelopment in human beings. The characterization of the brain morphology and its biological implications on a specific subject is a complex task which requires efficient computational approaches. Radiology has traditionally assessed the main brain changes in different alterations from a macroscopic point of view, thus, not considering subtle changes as a results of neuronal plasticity. In conjunction with functional information, the structural neuroimaging methods have established as the key in the diagnosis of several central nervous system disorders, including tumours, neurodegenerative disorders and psychiatric diseases.

2. Brain morphometry

2.1. Introduction

Morphometry techniques use statistical methods to detect and to quantify subtle structural abnormalities that appear when comparing different populations. Nowadays, there are several methodologies which have been designed to achieve these goals. The fast evolution in terms of spatial resolution and signal-to-noise ratio in Magnetic Resonance (MR) scanners as well as the improvements on new imaging techniques and data processing algorithms, help to developing studies that increase the knowledge over many fields of neuroimaging. This section describes the scope of these new methodologies and the main processes related with their implementation.

2.2. Morphometric methods

The first method developed in order to measure anatomical differences was based on the manual delineation of brain structures and their analysis by defining regions of interest (ROI). Although its main advantage is the anatomical accuracy of the measures, there are some assumptions that should be taken into account, including high variability, poor reproducibility, the need for previous hypotheses about the anatomical areas and regions to study and computational requirements needed to study a large number of subjects.

In order to supply these restrictions, semiautomatic methods have been developed. These methodologies perform a fully computerized treatment of different brain areas, providing a reproducible way to define exploratory analysis without a priori knowledge about the spatial distribution of the potentially affected areas.

2.2.1. Deformation-based morphometry

Models based on deformation fields use the spatial transformations needed to register an image to a template. In this registration process, a three-dimensional nonlinear deformation map is generated, which contains the adjusted parameters obtained by the fitting process between both, the image and the template. The deformation-based morphometry (DBM) (Gaser et al., 2001) is therefore a useful methodology to find differences at the macroscopic level.

To obtain the deformation field, the algorithm is initialized and a first mesh is generated. At each iteration, this mesh is fitted to achieve the required target varying from low to high detail by a coarse to fine minimization strategy. This registration is followed by an estimation of the nonlinear deformations which are composed by a linear combination of 3D discrete cosine (DC) transform basis functions. Displacement vectors are then smoothed with and 8 x 8 x 8 mm Full Width at Half Maximum (FWHM) filter.

The statistical analysis of these parameters helps to determine whether there are specific differences between subjects. The deformation field provides information about both volume and position differences, and can be studied by analyzing the displacement vectors for each point or by quantifying the local signal variation. Multivariate statistical models are needed in order to make inferences about the differences between groups.

2.2.2. Tensor-based morphometry

Tensor-based morphometry (TBM) (Kipps et al., 2005) is a morphometric method which uses tensor magnitudes to identify regional changes in anatomical areas The estimation of these differences is based on the small variations that are generated when normalizing each voxel of an image (i_a, i_b, i_c) to a template reference (j_a, j_b, j_c).

By using the deformation fields, the determinants of the Jacobian matrix (J) can be estimated. This matrix is equivalent to a second-order tensor that provides univariate (point to point) information about how the brain shape varies from the original image to the template. This feature improves the use of the DBM method, because it avoids the use of the entire deformation field (multivariate approach) in order to determine if there are specific (local) differences between images.

For each voxel, the Jacobian matrix contains information about translation, rotation and shear transformations:

$$J=\left[\begin{array}{l}\partial i_a/\partial j_a\partial i_a/\partial j_b\partial i_a/\partial j_c\\ \partial i_b/\partial j_a\partial i_b/\partial j_b\partial i_b/\partial j_c\\ \partial i_c/\partial j_a\partial i_c/\partial j_b\partial i_c/\partial j_c\end{array}\right]$$E1

In order to perform the statistical analysis to detect differences between subjects, the J matrix can be decomposed into a rotation matrix (R) and symmetric positive-definite matrices (U or V), to satisfy the following axioms:

$$J=RU=VR$$E2

$$U=\sqrt{\left(J^{T}J\right)}$$E3

$$V=\sqrt{\left(J{J}^T\right)}$$E4

$$R=J{U}^{-1}=V^{-1}J$$E5

In the normalization process which applies rigid transformations, it holds that U = V = I, where I is the identity matrix. If U and V matrices are different from the matrix I, there is a change in the shape which can be encoded by the tensor E. For a given deformation, there are infinite ways to express the associated tensors in terms of an n-parameter:

$$\{\begin{array}{l}n\ne 0\Rightarrow E^n=n^{-1}(U^n-1)\\ \\ n=0\Rightarrow E^n=\ln (U)\end{array}$$E6

When n = 0, the obtained tensor is the Hencky tensor, which is useful to express local brain volume increases or decreases relative to the template. From these calculations, a comparison between images of many subjects can be done by extracting different variables (area, length and volume) and analyzing those using statistical models.

2.2.3. Diffeomorphic morphometry

This methodology is based on registering an image with a template using a flow field that encodes the geometric transformation required to normalize an image to another. A large deformation framework is used in order to conserve topology, obtaining a diffeomorphic and invertible deformation (Ashburner, 2007).

If there are two images A and B (with the same dimensions) and a function f which takes points from A and put those on B, then f can be considered as a translator; i. e. for each point of A provides the corresponding B-point. In order to maintain the diffeomorphic propriety, this function must be bijective; i. e. the relationship between A and B points must be 1 to 1 (a specific point of A only can be associated to a specific point of B and vice versa).

$$\{\begin{array}{l}f:A\to B\\ f^{-1}:B\to A\end{array}$$E7

where the transformation is diffeomorphic if there is a smooth bijective function f that transforms A into B and vice versa. The last step is to estimate a statistical model to detect significant changes between groups.

2.2.4. Voxel-based morphometry

The voxel-based morphometry (VBM) technique (Ashburner and Friston, 2000) is based on the normalization of several individual brains with to a specific template. These normalized images are voxel-by-voxel analyzed to detect variations of local tissue. Unlike other morphometric techniques, VBM is based on applying a mass univariate statistical analysis for each voxel. Typically, the brain is previously segmented into gray matter (GM), white matter (WM) and cerebrospinal fluid (CSF) maps. These calculations need a prior preprocessing to normalize the data in a common stereotactic space. The purpose of these processes is the minimization of the anatomical variability between different subjects, discounting macroscopic factors and allowing a statistical analysis to obtain subtle differences that can be attributed only to the anatomical variability between groups.

2.2.4.1. Signal heterogeneity

This step aims to minimize the bias field contained in the MR images. The lack of signal homogeneity, which may result from factors such as static magnetic field inhomogeneities, sensitivity of transmit and receiving coils and dielectric effect, directly affects the voxel intensities. In order to quantitatively evaluate the data, differences in the brightness between voxels of a particular region or area can be a source of bias for the algorithms convergence criteria (figure 1).

Figure 1.

Low-frequency bias field estimated from brain MR images

For these reasons, it is necessary to correct this inhomogeneity and there are several approaches, including modeling of field heterogeneity by DC basis functions (Ashburner and Friston, 2000), the use of Legendre polynomial basis functions (Brechbühler et al., 1996) or the Gaussian deconvolution on the histogram of the image (Sled et al., 1998). There are also methods which model the field by a linear combination of low frequency functions based on cubic B-splines adjusted by a cost function based on the intensity and the gradient of the image (Manjón et al., 2007).

2.2.4.2. Non-linear normalization

The nonlinear normalization (or warping) normalizes an image to a template by applying transformations that do not preserve the proportions of the original image. The main aim is to perform a deformation of an original image with a template to facilitate a high precise comparison within brain regions between different subjects.

The algorithm tries to reduce the difference between original and template images, using a standard least-squares minimization (Mean Square Error, MSE):

$$MSE={{\displaystyle \sum _i\left(f(x_i,y_i,z_i)-w\cdot g(x_i\text{'},y_i\text{'}{z}_i\text{'})\right)}}^2$$E8

where f(x_i, y_i, z_i) represents the value of the voxel i in the coordinate (x,y,z) of the original image f, g(x_i’, y_i’, z_i’) is the value of the voxel i in the coordinate (x’,y’,z’) of the template g and w represents a weighting factor.

2.2.4.3. Segmentation

The segmentation process aims to classify the MR brain images into GM, WM, CSF and other cortical and subcortical areas. Although there are many algorithms for brain segmentation, there is an efficient strategy commonly used by neuroimaging applications that in practice gives good results but theoretically is slightly away from the pure concept of segmentation. This method does not obtain the real tissue-intensity extracted from the image but a probability map for each class. Each voxel in these maps has a normalized brightness value in the range [0... 1], that reflects the probability of belonging to a particular tissue.

In order to identify and classify the different tissues, the algorithm analyzes the range of the brightness values of each voxel in the original image. If n is the number of bits of the image, then the intensity values can be assigned in the range [0... 2ⁿ-1]. For example, a coded image with 8 bits, has a brightness value between 0 and 255, with 0 black (no light) and 255 white (figure 2). With this approach, it is possible to represent images using cumulative graphs (histograms) in which each point represents the number of voxels with a given brightness level.

Figure 2.

Gray scale with 256 potential values (0 black, 255 white)

The intensities can be modeled by its mean and variance, adjusting the image histogram by a Gaussian-mixture function,. The algorithm performs a separate treatment for each tissue, assigning a different group (class) to each voxel (figure 3). Initially, these voxels are assigned to an initial value defined by a priori knowledge. Then, the algorithm obtains the total number of voxels in each group and their mean and variance. With these data, new iterations are recalculated:

$$p_{i,k}=\frac{1}{\sqrt{\left(2\pi c_k\right)}}\exp \left(\frac{-{(f(x_i)-v_k)}^2}{2c_k}\right)$$E9

where p_i,k represents the probability that the voxel i is assigned to the k-tissue, c_k is the variance of the tissue k, f(x_i) represents the brightness of the i-voxel in the image f and v_k is the mean of the k-tissue. With the new probabilities, the algorithm continues until either the convergence criterion is achieved or the fixed number of iterations is exceeded.

Figure 3.

Segmentation process. Top: original image. Bottom (from left to right): gray matter, white matter and cerebrospinal fluid probability maps.

2.2.4.4. Smoothing

The main purpose of the smoothing process is to increase the signal-to-noise ratio by reducing the high-frequency random noise. Additionally, smoothing involves other advantages such as increasing the normality of the data and the minimization of inter-subject anatomical differences. The smoothing kernel fixes the brightness of each voxel taking into account the Gaussian average of their adjacent voxels (neighbors). The filtered image is then blurred, mainly in edge and contour areas because the high-frequency signals are removed while the low-frequency bands are preserved (figure 4). The main parameter which defines the shape of the filter is the standard deviation (σ) expressed as the total amplitude at FWHM:

$$FWHM=\sigma \sqrt{\left(8\cdot \log (2)\right)}$$E10

Figure 4.

Smoothing of an image. From left to right: original image, 2D Gaussian kernel and smoothed image.

2.2.5. Statistics and results

Usually, the statistical analysis that follows the application of morphometric techniques is based on the General Linear Model (GLM) (Friston et al., 1995). This model allows statistical inferences selecting specific effects of interest in the study groups and is based on an equation that defines the measured signals by a linear combination of explanatory variables plus an error whose distribution is (assumed) Gaussian:

$$Y=X\beta +\epsilon$$E11

where Y represents the measured data, X models the design matrix, β represents the estimated parameters and ε is the error.

This structure allows the definition of a measured variable Y as a linear combination of explanatory variables plus an error. It is assumed that this error is independent and follows a Gaussian distribution with zero mean. The design matrix X is a model structure which includes covariates of interest that could potentially influence the results (age, sex, clinical scales, overall tissue volume,...).

By using a voxel-by-voxel approach, multiple statistical comparisons are tested. So, it is necessary to apply additional corrections to minimize the presence of false positives (type I errors). This problem can be solved by applying specific corrections to ensure the reliability of the results. In this sense, the Bonferroni correction based on setting the significance criteria to α / number of observations or the False Discovery Rate (FDR) (Genovese et al., 2002) that controls the fraction of false positives, can be used. The obtained maps are then colored and overlaid over a high resolution T1 image showing the morphometric differences between groups (figure 5).

Figure 5.

From left to right: original statistic map, colored statistic map and overlay over a T1 axial MR image.

2.3. Structured report

The final step of the morphometric procedure is to include all the information in a structured, concise and brief report (Marti-Bonnmati L., 2011). This report lists all the variables and numerical data calculated in the different processes:

• Parametric data

• The report should include the parameters used in the morphometric method (type of technique, normalization, segmentation, smoothing, templates…) and statistical information (type of test, thresholds, p-values,…).

• Volumetric measures

• Overall volumes of GM, WM and CSF and absolute (ml) and relative (%) values must be included. Furthermore, volumetric measurements for subcortical areas (for example, basal ganglia) are desirable. These values are compared with normal values (obtained from healthy subjects) after discounting potentially relevant sources of bias (age, sex, laterality,…)

• Figures, coordinates and labels for each area of interest

• The final report should also incorporate the significant areas showing differences between groups and their associated values and coordinates. If any, these areas should be overlaid onto a T1 template and detailed in a table which shows statistic values, location of the affected regions (including Brodmann areas) and cluster volumes.

3. White matter tractography

3.1. Introduction

The diffusion tensor magnetic resonance imaging (DT-MRI, DTI) technique is widely used nowadays to explore the anatomy of white matter tracts in the human brain in vivo. The DTI is a non-invasive technique that permits the visualization of white matter fiber bundles by the reconstruction of their trajectories in a voxel-by-voxel basis through the measurement of water diffusion along different directions. Water molecules movement in white matter is restricted by the axon and the longitudinal arrangement of myelin covering the axon. In these situations, where a main direction predominates over the others, the water molecules movement has a high anisotropy.

The DTI technique permits the acquisition of MR diffusion images with different orientations of the magnetic field gradients, thus, obtaining a set of images with information of the water movement directionalities for each anatomical cut. The number of gradient orientations is a key parameter in the acquisition of DTI data and, although it is mathematically enough to have 6 directions in order to calculate a tensor, a higher number of directions provide a higher directional resolution.

The computational processing of the DTI data permits the calculation of the orientation and fractional anisotropy (FA) voxelwise. In fact, FA parametric maps can be generated to depict the main orientation of the white matter structure. Computational algorithms specially designed for fiber tracking can be applied to the orientation and anisotropy data in order to reconstruct the trajectory of white matter tracts.

The DTI has a unique view of the tissue architecture of neurons and changes associated with various pathophysiological alterations. There is an increase in the use of this technique for the analysis of white matter alterations produced by tumours and the corresponding surgery planning. Also, the study of congenital abnormalities of the corpus callosum and cerebellum, epilepsy, schizophrenia and early and late Alzheimer's disease is being widely assessed by this technique (Catani M., 2006).

The DTI can be combined with other MRI techniques, such as conventional T1 and T2 images, MR perfusion studies or the results of the concentration of metabolites composing fiber bundles obtained from MR spectroscopy.

3.2. Principles of diffusion tensor MRI

The phenomenon of molecular thermal motion results in random movement of molecules in the three directions of space. These displacements are considered, in general, as translational motions of molecules characterized by Brownian nature. This movement or molecular diffusion in the human body takes place mostly between water molecules.

In some tissues of the human body, water molecules can present a free movement without barriers, also known as free diffusion, or a movement which is limited by the structure of the neighbouring tissues, known as restricted diffusion. The figure 6, shows both concepts.

In general, diffusion measurements express the effective displacement in space of the water molecules in a certain time interval (Le Bihan D., 1988). Although temperature modulates the molecular motion of water molecules (approximately 2.4% per degree Celsius) (Tofts PS., 2000), that thermal influence is not significant in the study of diffusion, as there are other biophysical properties that have a significant effect on the mobility of tissue water.

Figure 6.

Schematics of diffusion of water molecules in a free environment (left) and in a restricted environment (right).

If pure water is used as a reference standard, the average displacement of the water molecules in a range of about 50ms does not exceed 20µm. Because this dimension is comparable to the cell dimensions, there is a high probability that the water molecules also interact with intracellular components, hydrophobic membranes and macromolecules that impede the movement of water. Therefore, the "apparent" diffusion is several times lower than in the case of pure water. In biological systems, diffusion comprises a complex mixture of single thermal diffusion with exchange between the intracellular and extracellular compartments through cell membranes and tortuosity of the interstitial space, which is conditioned by cell size, organization and density clustering.

To understand diffusion and its quantification, it is assumed that in the initial time we have a group of molecules concentrated at one point. If we wait a time t, without exerting any action on the molecules, they will have expanded in the three dimensions following the Einstein's equation of diffusion:

$$r^2=6·D·t$$E12

where t is the time interval and r is the average radius of the distribution. As can be deduced, the diffusion coefficient D is expressed as units of distance squared per unit time. For use in radiology or clinical applications, it is usually expressed in mm²/s.

In order to study the physical diffusion properties explained above, MRI is the only imaging modality that allows visualization and calculation of molecular diffusion in vivo directly from molecular translational movement of water.

MR signal is sensitive to microscopic movements water molecules. During the de-phase of the spins after the radio frequency (RF) pulse, phase incoherencies appear in the spins relaxation due to thermal agitation of the water molecules, which accelerates the loss of spins synchronism and reduces the relaxation time. Moreover, the repeated movement of water molecules cause the nuclear spins displacement to other regions in which magnetic field differs from the original value, thus causing a frequency modulation of relaxation.

In order to quantify the displacement movement of the spins independently, a gradient in one direction can be applied immediately after the pulse. In this situation, the water molecules which have been moved in the direction of the gradient will be under a magnetic field be farther more of the original and therefore the signal drop faster.

The free diffusion approximation of light in the previous sections cannot be assumed in biological tissues, because sometimes, the movement of water molecules is restricted or defined to a certain direction. In the latter case, in which a molecule is most likely to move in one direction than another, one speaks of an anisotropic diffusion. The most obvious example (as will be seen below) takes place in the cerebral white matter, where the water molecules tend to move along axonal tracts of the different fascicles brain.

This anisotropy of diffusion can be characterized mathematically, considering a diffusion tensor in the following matrix form:

$$\left(\begin{array}{ccc}{D}_{xx}& D_{xy}& D_{xz}\\ D_{yx}& D_{yy}& D_{yz}\\ D_{zx}& D_{zy}& D_{zz}\end{array}\right)$$E13

Since the matrix is symmetric, ie Dxy = Dyx; Dyz = Dzy and dxz = Dzx, simply calculate 6 of the 9 parameters. Therefore, we can deduce that to extract directional properties of diffusion, it will require at least 6 different gradient directions.

An example of 6 acquisitions can be appreciated in figure 7.

Figure 7.

Diffusion images acquired in different magnetic field gradient orientations for the calculation of the diffusion tensor.

The eigenvector of the diffusion matrix provide the information about main orientation of the water molecules movement in each voxel. An example of orientation maps at different detail scales can be appreciated in figure 8:

From the diffusion matrix, the fractional anisotropy (FA) parameter can be calculated from the expression:

$F_f=\sqrt{\frac{{\left({\lambda }_1-\overline{\lambda }\right)}^2+{\left({\lambda }_2-\overline{\lambda }\right)}^2+{\left({\lambda }_3-\overline{\lambda }\right)}^2}{{\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2}}$E14

Being λ₁, λ₂ and λ₃ the eigenvalues of the diffusion matrix, and the average value of the eigenvalues. An example of the combination of both orientation and FA information in a voxel-by-voxel basis may be appreciated in figure 9.

Figure 8.

Orientation maps calculated from the diffusion matrix. From left to right: full brain map showing the vector field with the main orientations for each region. Detail of the vector field in a selected region. Voxel-by-voxel representation of the main diffusion orientation.

Figure 9.

Combined fractional anisotropy (FA) and orientation map. The level of intensity expresses the FA value, while the color indicates the main orientation (LR: left-right, AP: anterior-posterior, FH: foot-head).

3.4. Fiber segmentation

Different segmentation strategies exist for white matter tractography reconstructions of the fibers. The main segmentation techniques can be divided in seed based segmentation, regions of interest segmentation and white matter atlas segmentation.

3.4.1. Seed segmentation

This technique considers an initial point in a 3D space with a given orientation and FA. Thus, the algorithm will initiate a path by the neighbouring voxels showing similar orientations. This trajectory will be calculated until a too sharp angle exists between the orientation of the current voxel and the following. The fiber trajectory calculated will be reconstructed unless if it accomplishes also the condition of the minimum length, which is another parameter imposed in the segmentation to avoid the reconstruction of small fibers from random noise.

3.4.2. Regions of interest segmentation

This is the technique with a higher use in clinical routine nowadays. White matter fibers are reconstructed from regions of interest (ROIs) which are placed according to the user anatomical knowledge. This technique allows for the calculation of the fibers that pass through the ROIs that have been introduced. Exclusive ROIs can also be placed in order to avoid the reconstruction of fibers bundles which are adjacent to the one of interest.

An example of this technique can be appreciated in figure 10, where the uncinate fasciculus is reconstructed. Two ROIs are placed in order to exclusively reconstruct fibers crossing both regions.

Figure 10.

Segmentation of the uncinate fasciculus by the placement of two ROIs, in the frontal and temporal lobes.

3.4.3. Atlas based segmentation

The tracts segmentation using white matter atlas has a higher complexity. In general terms, the main basis of the method consists in the calculation of the orientation and FA maps for large series of subjects. All these data is anatomically co-registered and a final expert anatomical labelling is performed (O'Donnell LJ., 2007).

3.5. Extracted parameters

The main white matter tracts can be segmented routinely by ROI segmentation for clinical applications. The authors suggest the segmentation of the following white matter fasciculum according to experience with pre-surgical evaluation and study of neurodegenerative disorders:

• Corpus callosum

• Cingulate fasciculus

• Uncinate fasciculus

• Corticospinal fasciculus

• Inferior longitudinal fasciculus

• Superior longitudinal fasciculus

In figure 11, examples of fiber reconstructions in different pathologic conditions can be appreciated.

Figure 11.

White matter fasciculus reconstruction in different clinical cases. In a), main white matter fibers segmentation for the pre-surgical evaluation of a glioblastoma multiforme. The right superior longitudinal fasciculus, in blue, can be appreciated to be attached to the tumour periphery. In b), the reconstruction of a sectioned cingulum is observed in a patient after an emergency intervention due to an acute hydrocephaly. In c), corpus callosum fibers shortening due to advanced multiple sclerosis lesions.

In each reconstructed fiber bundle we can extract a set of parameters related to the microstructure:

• Fractional anisotropy (FA): its value ranges from 0 (pure isotropic) to 1 (highly anisotropic) and shows the degree of existence of a preferential diffusion direction within the voxel.

• Mean diffusivity (D): it is measured in mm²/s and expresses the degree of restriction to water molecules movement in a voxel. High D values reflect low degree of restriction to movement, while low D values show a restricted diffusion of molecules due to a higher cell density and reduced interstitial space.

• Number of fibers (NF): it is the total number of fibers that have been reconstructed in a certain fasciculum.

• Average length (L): it is mostly expressed in centimetres and provides the average length of the fibers of the reconstructed fasciculum.

3.6. Structured report

An adequate tractography report should be brief and concise (Marti-Bonmati L., 2011), and include:

• Parametric data: the results of the parameters presented in the anterior section (FA, D, NF, L) for each reconstructed white matter fasciculum. These values should be compared to values obtained in a large series of age-matched healthy subjects.

• Figures: representative figures of the main white matter tracts superimposed on anatomical images.

4. Conclusions and future challenges

The brain morphometry and tractography techniques have established as the main image processing methodologies for the characterization of brain structure in all types of central nervous system disorders. Although many centres benefit from their application to different clinical cases, large population studies have been mainly limited due to lack of standardization in the acquisition, processing and reporting techniques. The future challenges for these techniques have to be focused in multi-centre initiatives that facilitate the protocols sharing, the standardization of analysis procedures and the way this information is presented in adequate structured reports.