|Home | About | Journals | Submit | Contact Us | Français|
Diffusion tensor imaging (DTI) is sensitive to the directionally- constrained flow of water, which diffuses preferentially along axons. Tractography programs may be used to infer matrices of connectivity (anatomical networks) between pairs of brain regions. Little is known about how these computed connectivity measures depend on the scans’ spatial and angular resolutions. To determine this, we scanned 8 young adults with DTI at 2.5 and 3 mm resolutions, and an additional subject at 4 resolutions between 2–4 mm. We computed 70×70 connectivity matrices, using whole-brain tractography to measure fiber density between all pairs of 70 cortical and subcortical regions. Spatial and angular resolution affected the computed connectivity for narrower tracts (internal capsule and cerebellum), but also for the corticospinal tract. Data resolution affected the apparent role of some key structures in cortical anatomic networks. Care is needed when comparing network data across studies, and interpreting apparent disagreements among findings.
DTI and its variants (HARDI, DSI) are greatly advancing our understanding of neurological and psychiatric illness. DTI is widely used in international efforts to map brain connectivity in populations of thousands of subjects, such as the $30 million Human Connectome Project (www.humanconnectomeproject.org). Much of the effort to map brain networks has focused on DTI and its extensions, which can map axonal pathways and tracts in the living brain. Diffusion MRI is sensitive to the local direction and rate of water diffusion, at each location in the brain. Axonal pathways may be reconstructed using tractography methods to infer the most likely paths of tracts, based on tensor or higher-order (q-space) diffusion models . Other methods exist to map functional connectivity, but the meaning of connectivity is different – in resting-state fMRI and MEG/EEG, temporal correlations are measured between signals at pairs of locations in the brain. These correlation matrices are then thresholded to identify network “hubs” important for functional synchronization of brain activity.
Diffusion MRI connectivity mapping has broad applications in neurology and psychiatry for understanding disrupted patterns of brain connections in Alzheimer’s disease, autism  and neurogenetic disorders of childhood, as well as sex differences and genetic effects on connectivity . Much of the technical innovation focuses on q-space imaging, which enriches the information available on local directional and radial diffusion. The quest to improve the local diffusion model has led to elaborate q-space sampling schemes with large numbers of directional samples (HARDI), and/or multiple diffusion weightings (b-values) which are sensitive to non-monoexponential radial diffusion [4–6].
Several empirical and theoretical studies have modeled how signal to noise in DTI depends on the spatial and angular resolution of the scans [7–8]. Others have optimized the q-space sampling to boost SNR in clinically feasible scan times . Much less attention has been devoted to understanding how brain networks, and patterns of recovered connections, depend on the spatial and angular resolution of the scans. Spatial resolution affects even the simplest measures derived from DTI, such as fractional anisotropy – the most widely used measure of brain integrity; in larger voxels, FA measures can be greatly reduced by partial volume effects. Here we scanned 9 young adults to monitor how different scanning protocols affect the recovery of maps of brain connectivity.
We hypothesized that both spatial and angular resolution would affect measures of cortical and subcortical connectivity. We were interested in the anatomical scope and extent of these effects. Empirical data on these questions will help to determine how well scan data may be pooled or compared in multi-site DTI studies; many of these are now being planned or are underway. Those designing DTI protocols may also be interested to know how comparable their brain connectivity maps are likely to be, relative to independently collected data from other imaging centers.
Our first multi-resolution DTI dataset (called the “Mayo dataset”) was collected at the Mayo Clinic. 8 healthy subjects (age: 32.0 years ± 3.9 SD; 4 males) were scanned on a GE 3 Tesla brain MRI scanner with an 8-channel head coil. DWI data was collected using two protocols, summarized in Table 1.
For our second data set, a healthy male subject (32 years old) was scanned at multiple spatial resolutions (2×2×2, 2.5×2.5×2.5, 3×3×3 and 4×4×4 mm3) with axial DTI using the following acquisition parameters: TR=8000 ms, TE=83 ms, 128×128 matrix, 64 slices, b-value=1000 s/mm2, one baseline (b0) scan and 12 gradient directions. Although the number of diffusion directions for this second dataset is less than state-of-the-art (12), the large number of values of the spatial resolution makes it of interest.
We first eddy-corrected all EPI images using FSL (http://fsl.fmrib.ox.ac.uk/fsl). Geometric distortions due to magnetic susceptibility were corrected using a simultaneously collected field map, using the FSL “Prelude and Fugue” function. Non-brain regions were removed from a T2-weighted image in the corrected DWI dataset using the bet function in FSL.
Depending on the number of the gradient directions in each dataset, we used two different models to analyze the preprocessed DWI. For the 12-direction data, the standard diffusion tensor was fitted using the Diffusion Toolkit (http://trackvis.org/dtk/). More gradient directions were collected for the Mayo datasets (Table 1), so the orientation distribution function (ODF) for water diffusion was reconstructed voxelwise using the recently proposed Constant Solid Angle formula (CSA-ODF; Aganj et al., 2010) with 4th-order spherical harmonics (SH):
Here û is the arbitrary unit vector. Ẽ (û) := S(û)/S0, where S(û) is the diffusion signal, and S0 is the non diffusion-weighted signal (from the b0 scans). is the Laplace-Beltrami operator and FRT is the Funk-Radon transform (Funk, 1916). ODFs were reconstructed using 642 point samples, from a recursively subdivided icosahedral approximation of the unit sphere.
50 white matter (WM) regions of interest (ROIs; Table 2) were generated based on the ICBM young adult DTI-81 atlas (http://www.loni.ucla.edu/ICBM/). ROIs were transformed into our datasets’ native coordinate space using 9-parameter linear registration. Using FSL-FLIRT, we registered the atlas T2-weighted image to the b0 image of each individual dataset to obtain a transformation matrix. This was applied to all atlas WM ROIs to map them onto the target datasets.
White matter connectivity was then computed from a fast-marching based method . Given a seed region, the generalized fast-marching (GFM) method evolves a front outward according to the alignment of underlying fiber distributions. The GFM algorithm accept ts any spherical measure of fiber orientation – such as conventional tensors or ODFs – these are normalized to unit mass and sampled in directions coinciding with the displacement vectors to each voxel in the 98-neighborhood of a voxel under consideration (Fig. 1). The 98-neighborhood consists of all neighbors in the dilated 3×3×3 volume, plus those in the 5×5×5 volume whose displacement vectors do not coincide with any of those in the smaller neighborhood; this expansion of the neighborhood beyond the closest 26 voxels reduces the discretization error. At each iteration, the GFM algorithm adds one voxel to the advancing front; we decide which voxel to include based on a speed function that accounts for anisotropy and collinearity of fiber orientations at the source and destination voxels (Eq. 2):
where x is a voxel on the boundary of the front, xn is its n-th neighbor just outside the front, g is the mean generalized fractional anisotropy (GFA) between x and xn, un is the normalized displacement vector to the n-th neighbor, ue is the entrance direction with which the front moved to include x, and f(x,un) indicates the fiber orientation probability in the n-th direction at voxel x. The two terms in this equation are weighted differentially by the generalized fractional anisotropy (GFA) – in highly anisotropic regions, fiber orientation information is more reliable and the first term takes precedence, as it quantifies fiber orientation collinearity. In low anisotropy areas, the second term is favored, which prefers that the front continue to move in its current direction. Given this velocity, the front arrival time at each neighbor is computed simply from Eq. 3:
The voxel along the edge of the advancing front with the earliest arrival time is accepted into the front, and the neighbors of the newly-incorporated voxel are added to the list of potential transitions to consider in the next iteration. Marching continues in this manner until the front has evolved through the entire brain.
Connectivity is quantified by tracking a number of measures during the GFM evolution. The front arrival time is stored for each voxel along with the path distance taken to incorporate the voxel, the mean anisotropy along that path, the mean transition velocity of the front towards each voxel, and the minimum velocity step along the path to each voxel. Maps of all these attributes imply various features of “connectivity”, and are combined by standardizing the maps (subtracting the mean and dividing by the standard deviation across each subject’s map), performing principal component analysis (PCA), and extracting the projection of each voxel’s attributes onto the principal component vector. This score provides a robust “connectivity index” between WM ROIs.
Cortical connectivity matrices were computed as in . 35 cortical labels per hemisphere (Table 3) were automatically extracted from the same subjects’ T1-weighted structural MRI scans using FreeSurfer (http://surfer.nmr.mgh.harvard.edu/). As the software performs a linear registration, the resulting T1-weighted images and cortical models were aligned to the original T1 input image space and down-sampled to the space of the DWIs, using nearest neighbor interpolation (to avoid mixing of labels). To ensure tracts would intersect labeled cortical boundaries, labels were dilated with an isotropic box kernel of 5 voxels.
Fiber tracking was initiated by specifying three parameters: the threshold values for starting and stopping tracking, and the critical angle threshold for stopping tracking if the algorithm encounters a sharp turn in the fiber direction. Based on the diffusion tensor model, the Diffusion Toolkit (http://trackvis.org/dtk/) uses these parameters to generate 3D fiber tracts using the Fiber Assignment by Continuous Tracking (FACT) algorithm. With a brute-force reconstruction approach, we used all voxels in the volume as seed voxels to generate the fibers. After that, a spline filter was applied to each generated fiber, with units expressed in terms of the minimum voxel size of the dataset. After whole brain tractography, a 70×70 connectivity matrix was created. Each matrix element estimates the proportion of the total number of fibers, in that dataset, connecting each of the labels to the others.
We chose the P1 dataset (Table 1) for this test as it has the highest angular resolution (48 gradients). We downsampled the original dataset from 48 to k=47, 46… 16, 15 DWIs. Sub-sampling was based on maximizing the total angular distribution energy of the remaining set of k gradients, to optimize the uniformity of the spherical sampling (see ). We calculated white matter connectivity based on 4th-order SH CSA-ODFs (Eq. 1) for each subsampled dataset as in Section 2.3. Then we assessed how angular resolution affected white matter connectivity. (This angular subsampling slightly overestimates the consistency achievable in scans of the same subject across independent scanning sessions). But, by subsampling the same dataset, we can isolate the effect of angular sampling and model it, with other sources of variance held constant. Fig. 2 shows the standard deviation (a measure of instability) of the connectivity matrix elements among the connectivity maps calculated from subset 15 to subset 48; this standard deviation was computed in each of the 8 subjects, and then averaged across all 8, to infer general patterns. As expected, some of the thinnest (narrowest) fiber tracts – the cerebellar ICP and SCP, and the internal capsules – were strongly affected by altering the angular resolution. Even some of the major pathways, including the apparent connections of the cortico-spinal tract with the ACR, ALIC and SFO were also quite severely influenced (red colors, Fig. 2).
In general, changing the spatial resolution will also change the SNR for any DTI-derived measure. We chose the P2 dataset in Table 1 as the target, and gradually reduced its spatial resolution by downsampling its isotropic voxels of side 2.5 mm to 2.6, 2.7, …, 10 mm. Although other choices are possible, we chose linear interpolation to downsample the original P2 images to create each new image. For each downsampled subset, we calculated white matter connectivity based on 4th-order SH CSA-ODFs. Fig. 3 shows the standard deviation of connectivity matrix elements across connectivity maps calculated at all voxel sizes in the range 2.5, 2.6 …10 mm, averaged across all eight subjects. The computed WM connectivity in all tracts and all regions is affected by partial volume effects. Greatest differences were found in the connections of the medial lemniscus, cerebellar peduncles, internal capsules, which are among the thinnest tracts.
As 4th-order SH CSA-ODFs require at least 15 gradient directions to compute, we instead used the diffusion tensor model to evaluate white matter connectivity for the single subject scanned many times with only 12 gradient directions, but at multiple spatial resolutions. Fig. 4 shows the standard deviation of connectivity among connectivity maps calculated at 4 different isotropic spatial resolutions (2, 2.5, 3 and 4). These maps show more differences than those in Fig. 3, as signal averaging was used to boost the SNR for the scans with smaller voxels. By contrast, Fig. 3 is based on downsampling the exact same set of scans in 8 subjects, rather than performing new scans.
Fig. 5 shows the standard deviation of elements in cortical connectivity matrices for 70 ROIs in the 12-direction dataset, at 4 different spatial resolutions. The computed pattern of cortical connectivity heavily depends on the spatial resolution, with less apparent connectivity in scans with large voxels. The cortical connection between parahippocampal and fusiform gyri, and between corresponding structures in the left and right hemispheres were most affected by spatial resolution.
Overall, scans with larger voxels were prone to partial volume artifacts due to under-sampling, that causes a “loss” of some connections, especially for narrower tracts in the cerebellum and internal capsule. Differences were prominent throughout the brain. Care is needed when (1) interpreting anatomical connectivity patterns as objective measures of biological connectivity, (2) pooling data across scanners, and (3) comparing studies.
This project is funded by NIH grants R01 EB008432 and R01 EB007813 to P.T.