Validation of DTI Tractography-Based Measures of Primary Motor Area Connectivity in the Squirrel Monkey Brain

Tracer Injection

The experiment was carried out on a New World squirrel monkey. The brain of this species shares many features of functional organization and microstructural complexity with humans, although the cortical surface is less gyrated. The motor cortex, for example, is completely exposed on the brain surface, so it is easily accessible for electrical microstimulation and injection of tracers. Other advantages of New World monkeys are that the cortical connections of the primary motor cortex (M1) are well known from previous studies [33] and their smaller brains reduce the time of histological processing and anatomical data analysis.

Under general anesthesia using aseptic techniques, a bidirectional tracer, biotinylated dextran amine (BDA; Molecular Probes Inc., Eugene, OR) was injected (as a 10% solution in phosphate buffer) into left hemisphere M1. Pressure injections of BDA were carried out using a 2 ul Hamilton syringe. Eight injections (1 µl/site; indicated by the black squares in Fig. 1) were made to cover a large M1 region representing the forearm and identified by intracortical microstimulation (see the stimulation sites labeled ‘Fa’ in Fig. 1). After each injection, the needle was left in the brain for 5–10 minutes and then retracted stepwise to avoid leakage of the tracer along the needle track. Surgical, microstimulation and injection procedures have been described in detail elsewhere [33]. After surgery, the monkey was allowed to recover from the procedure, giving the tracer sufficient time to be transported along axons to all regions connected to M1.

All animal procedures were approved by the Vanderbilt University Animal Care and Use Committee, and followed guidelines of National Institutes of Health for the care and use of laboratory animals. The monkey was housed in adjoining individual primate cages allowing social interactions with other monkeys and received extensive environmental enrichment such as play objects that are changed for variability, novel food objects and foraging devices. Normal food, treats and fruit were provided to the animal by trained staff. Analgesics were given following the procedure and in response to any evidence of pain (behavioral abnormality, loss of appetite, tissue edema, or muscle spasms). The monkey was closely monitored in the week following the surgical procedures.

Ex vivo MRI Data Acquisition

One week after surgery, the monkey was given a lethal dose of barbiturate, and perfused through the heart. All blood was rinsed out with physiological saline (0.9% NaCl) followed by fixative (4% paraformaldehyde). The brain was removed from the skull and stored in buffered saline overnight. The next day, the brain was scanned on a 9.4 Tesla Varian scanner (Varian Inc., Palo Alto, CA). First, T2-weighted structural images were acquired by running a standard gradient echo multi-slice (GEMS) sequence with full brain coverage (TR = 963 ms, TE = 4 ms, flip angle = 20°, slice gap = 0 mm, voxel size = 300×300×300 µm³, data matrix = 128×128×192, SNR≈50). Then diffusion weighted imaging was performed using a pulsed gradient spin echo (PGSE, [2]) multishot spinwarp imaging sequence with the same FOV as the structural images (TR = 5.2s, TE = 26 ms, number of diffusion gradient directions = 31, b = 0, 1200s/mm², voxel size = 300×300× 300 µm³, data matrix = 128×128×192, number of acquisitions = 10, SNR≈25, scanning time≈50 hr). The b value used in this experiment was lower than is optimal for diffusion studies in fixed tissue [34], due to hardware limitations. A low b value decreases the available diffusion contrast-to-noise ratio (CNR) in the image data, which has the same effect as higher image noise. To compensate for this shortcoming, we extended the scan time to 50 hours, which yielded a CNR comparable to in vivo human studies (equivalent to an in vivo study with mean diffusivity = 0.7×10⁻³ mm²/s and SNR≈20).

Block/micrograph Data Acquisition

Following MRI scanning, and before sectioning, the entire brain was placed in 30% sucrose for cryoprotection. Two days later the brain was cut serially on a microtome in a coronal plane to produce 50 µm thick frozen sections. All sections were collected in phosphate buffer, but prior to cutting every third section (i.e., at 150 µm intervals), the frozen tissue block was photographed using a Canon digital camera (image resolution = 50 µm/pixel, image size = 3330×4000pixels, number of images = 286), mounted above the microtome to facilitate intermodality image registration (see section 2.4) [35].

Sections were divided into six series. Every sixth thin section was processed for BDA histochemistry [36], producing a series of 74 sections covering all regions connected to M1. These sections were photographed under 0.5× magnification using a Nikon DXM1200F digital camera mounted on a Nikon E-800 microscope (image resolution = 7 µm/pixel, image size = 3840×3072 pixels). More than one photograph was needed to cover the entire section. Those component photographs were merged into a panoramic micrograph using Adobe Photoshop CS3 (San Jose, CA) and the background of each panorama was extended to a net image size of 6660×6660 pixels. The resulting stack of images defined a “standard micrograph space” which is the target space of all higher resolution microscopic data. Another series of sections was stained for Nissl substance to identify the cytoarchitectonic borders of M1 and other cortical regions connected to M1. Four series of remaining sections were used for other purposes.

Micrograph-space to DTI-space Registration

The standard micrograph space was registered to the DTI image data using a multi-step procedure. First, every BDA-labeled section was down-sampled (to 256×256 pixels) and registered to the down-sampled photograph (256×256 pixels) of the corresponding tissue block using a 2D affine transformation followed by a 2D non-rigid transformation, semi-automatically calculated via the Thin-Plate Spline algorithm [37]. Next, all the down-sampled block photographs were assembled into a block volume and all the non-diffusion weighted MRI images were similarly stacked into a DTI volume. Then the block volume was registered to the DTI volume using a 3D affine transformation followed by a 3D non-rigid transformation automatically calculated via the Adaptive Bases Algorithm [38]. The multi-step registration described here was very similar to the registration procedure validated in an earlier study [39], which showed that the accuracy of the overall registration was approximately one MRI voxel (0.3 mm). The deformation fields produced by all the above registration steps were saved in order to transfer other data from the micrograph space (i.e., higher resolution microscopic data, described in section 2.5) into DTI space.

Fiber and Density Distribution Map (DDM) Reconstruction

Our ultimate goal in this part of the study was to reconstruct histological and DTI-tractography-based DDMs in DTI space and render them on the same white-gray matter (WGM) interface. Prior to calculating DDMs, we reconstructed the 3D WGM interface in DTI space: the T2w image volume was rigidly registered to the DTI volume; white matter was segmented using a variational level set approach [40] on the aligned T2w images; the 3D outline of white matter was extracted to generate triangle meshes representing the WGM interface. In addition, the BDA-injection region and different cortical projection regions were segmented manually in the standard micrograph space based on architecture revealed by Nissl- and BDA-labeled sections [33] and then transferred to DTI space.

To resolve each individual BDA-labeled fiber, 4× micrographs (0.87 µm/pixel) covering the WGM boundary were acquired using a Nikon D1R1 digital camera mounted on a Nikon AZ100 M microscope. The BDA-labeled fibers were segmented by a series of morphological processes: top-hat filtering was performed to correct uneven illumination, global thresholding was implemented to extract fibers, and a set of statistical properties (e.g., area of object and ratio of perimeter to area) of each extracted object was calculated and used to remove the non-fiber objects. For each micrograph, points along the WGM boundary were selected manually and fit to a curve with 6 pixels (∼5 µm) width to represent the histological WGM boundary. Segmented fibers that overlapped the fitted boundary were extracted automatically. The centroids of the overlapping fibers were rigidly transferred to the standard micrograph space, slice by slice.

Next, the standard size micrograph (6660×6660 pixels) was divided into 256×256 square units and the number of centroids in each unit area, n_i,j (i,j = 1,2,…,256), was determined. These n_i,j values served as the intensity of the 2D DDM. The DDMs were transferred from micrograph space to block space then to DTI space using the deformation field associated with the registration procedure described earlier (see section 2.4). Immediately following each spatial transformation step, the Jacobian matrix of the corresponding deformation field was calculated and then used to compensate the density change caused by the raw geometric transformation. Finally, to facilitate visualizing the density information in 3D, we mapped these gray-matter-distributed densities onto the 3D WGM interface in DTI space. Thus, ideally, the value of a certain voxel on the 3D WGM interface in DTI space represented the number of BDA-labeled fibers whose centroids were located on the WGM boundary in high resolution micrograph space.

To identify the locations of BDA-labeled somas in the micrographic volume, we exhaustively plotted the centers of these somas and the outline of the brain tissue for each BDA-labeled section under a 6.3× objective of a Zeiss microscope equipped with the vector graphing software Igor Pro 2.0 (Wavemetrics, Inc., Lake Oswego, OR). Using Canvas 11 (ACD Systems, Victoria, British Columbia, Canada), each vector graph was converted to a bitmap image with resolution high enough to resolve any two adjacent soma markers. The bitmap image was rigidly aligned with the 0.5× micrograph of the corresponding BDA-labeled section by matching the outline of the brain tissue and then was cropped into a standard size image. Next, each labeled soma was segmented by simple thresholding and the centroid of the soma was extracted automatically; the procedures (including gridding, counting, and calculating, then transferring and compensating the DDM) were similar to those for processing BDA-labeled fiber data. Finally, we mapped these gray-matter-distributed densities onto the 3D WGM interface in DTI space: the density value of every gray matter voxel was added to the closest voxel located on the 3D WGM interface. Thus, ideally, the value for each voxel at the interface was the sum of the number of BDA-labeled somas located vertically (i.e., normal to the surface) beyond this voxel.

DTIStudio [41] software was used to perform tensor fitting and fiber tracking (using the FACT algorithm, [6]) over the whole brain (tensor fitting method: standard linear-fitting; tracking parameters: start fractional anisotropy, FA = 0.1, stop FA = 0.2 and stop angle = 70°). Streamlines penetrating both of the following two ROIs were selected: (i) the BDA-injection region transferred from the micrograph space and (ii) the WGM interface just underneath (i). The above deterministic tracking scheme is referred to as the ‘DS’ scheme below. The DDM for this tracking method was produced by counting the number of the streamlines passing through or terminating within the voxels at the WGM interface. In addition, the streamline terminals distributed within the cortical regions connected to M1 were also mapped onto the interface and then counted for each interface voxel to produce a DTI ‘streamline terminal’ DDM.

Probabilistic fiber tracking was performed with the FMRIB’s Diffusion Toolbox (FDT) v2.0 from the commonly-used FMRIB’s Software Library (FSL, http://www.fmrib.ox.ac.uk/fsl). First, the bedpostx command was run to calculate the distributions of diffusion parameters at each brain voxel using the Markov Chain Monte Carlo sampling method [10] with automatic relevance determination (ARD) weight assigned to 1 (default) and 0.5, respectively. The probtrackx command was used to generate probabilistic streamlines (using the following parameters: sample number = 10000, curvature threshold = 0.2, modified Euler streamlining = on, step length = 0.1 mm and distance correction = on). The BDA-injection region was used as the seed mask, the WGM interface underneath the injection region as the waypoint mask and the gray matter mask (excluding the injection region) as the termination mask. These probabilistic tracking schemes are referred to as ‘FSL1’ (ARD weight = 1) and ‘FSL2’ (ARD weight = 0.5), respectively. The voxel value in the resulting dataset represents the number of streamline samples passing through this voxel times the expected length of these samples. The probabilistic-tractography-derived DDM was generated by keeping the values of those voxels located at the WGM interface and setting the remaining voxels to zero. Again, the DDMs were rendered on the 3D WGM interface in DTI space.

Quantitative Analysis

The number of DTI-tractography-derived streamlines has been proposed as a measure of cortico-cortical connectivity [18], based on the hypothesis that the number of streamlines is proportional to the number of axons connecting cortical regions. To test this hypothesis for inter-regional connectivity, the numbers of BDA-labeled fibers (N_B) passing into cortical ROIs and the numbers of DTI streamlines (N_D) projecting to the same ROIs were quantified. These data were fit to a linear model:

(1)

where a and b are, respectively, the estimated intercept and slope of the regression line, and e is the error term. If a = 0 lay within the 95% of confidence interval of the intercept, then the data were fit to a second tier linear model with intercept equal to zero:

(2)

The Pearson’s correlation coefficient of the fit, r, was calculated to indicate the degree of correlation between tractography and histology fiber data.

To verify the hypothesis in the case of intra-regional connectivity, we calculated the Pearson’s correlation coefficient, r_p, between tractography and histology fiber DDMs within each cortical projection region on a voxel-by-voxel basis:

(3)

where B_i and D_i are, respectively, the numbers of BDA-labeled fibers and DTI streamlines distributed in the i ^th voxel located at the WGM interface directly underneath this projection region, n is the total number of these voxels, and the bar represents the mean.

The original seed region covered the gray matter in the injection region, but it was found that DTI tractography was more effective when the seed region was intentionally extended one or two voxels into the white matter under the WGM interface (see below). To examine the effect of extending the depth, d_w, of the seed region into the white matter on quantitative metrics, we calculated all correlations using seed masks with d_w = 0, 1 and 2 voxels ( = 0, 0.3 and 0.6 mm).

In addition, since obtaining the DDM of BDA-labeled somas was less time-consuming, we studied the feasibility of using soma number to represent histological connectivity. We therefore calculated the relationship between numbers of somas and fibers across all projection regions. We also calculated the association between numbers of somas and numbers of DS streamline-terminals across all the projection regions.

Fiber Orientation Comparison

To elucidate the reasons for discrepancy between tractographic and histological connectivity, we compared the tensor orientation with BDA-labeled fiber orientation in the high resolution micrograph space. Based on a multi-step registration procedure [39], the locations of diffusion tensors were transferred from the DTI space into the high resolution histological space and the tensor orientations were corrected based on the preservation of principal direction algorithm [42].

M1-connected Cortices

BDA injections in the forearm representation area of the primary motor cortex reveal somatotopically distributed connections with the ipsi−/contralateral supplementary motor areas (iSMA/cSMA); ipsi−/contralateral anterior cingulate cortex (iAC/cAC); ipsi−/contralateral premotor cortex (iPM/cPM); ipsilateral M1 excluding the injection region (iM1ex); contralateral M1 (cM1); ipsi−/contralateral anterior parietal cortex (iPA/cPA); ipsi−/contralateral posterior parietal cortex (iPP/cPP); anterior part of upper bank of ipsilateral lateral sulcus –a combination of parietal ventral and rostral (iPVR) areas; posterior part of upper bank of ipsi−/contralateral lateral sulcus cortex–secondary somatosensory cortex (iS2/cS2). Labeling of ipsi−/contralateral prefrontal cortex (iPF/cPF) and contralateral PVR (cPVR) were not observed in this case. The above regions served as our regions of interest (ROIs). The borders of these ROIs were identified based on their cytoarchitectural features seen in Nissl preparations [43].

Detection and Counting of Cross-interface BDA-labeled Fibers

Figure 2 illustrates the pipeline for semi-automatic detection of cross-interface BDA-labeled fibers in the high resolution micrograph space. The WGM boundary (yellow curve shown in Fig. 2C) was estimated by fitting a series of manually-placed markers (red dots in Fig. 2B) on 4× micrographs covering a specific ROI. Interface-crossing fibers (fibers with color contours shown in Fig. 2D) were detected automatically by performing an “AND” operation on the fiber and boundary binary masks. Locations of the interface-crossing fibers were defined by the centroids of these fibers, and the number of cross-interface fibers was counted (red numbers in Fig. 2D). It was verified that shifting the WGM matter boundary towards white matter or gray matter by 50 pixels (∼44 µm) did not significantly change the number of interface-crossing fibers. The deviation of semi-automated counting from manual counting as the “gold standard” was also studied. Four to five sections for each ROI were randomly drawn for verification. The deviations for all the ROIs were in the range of 3%–5%.

Tensor Fitting Results

The DWI data were fit to single-tensors using DTIStudio. The primary diffusion direction and FA were calculated in each voxel. Figure 3 shows a color-coded map of primary diffusion direction and the FA scalar map of the same coronal slice.

Correlation of Inter-ROI Connectivity Strength

To test the hypothesis that the number of streamlines, N_D, connecting two cortical regions is proportional to the number of anatomical fibers, N_B, connecting those regions, we identified the number pair for connections between M1 and each ROI and plotted them in Fig. 4 as well as displayed them as connectivity backbones, shown in Fig. 5. The points lying on the horizontal axis (N_D = 0) and vertical axis (N_B = 0) in Fig. 4 indicate false negative and false positive connections, respectively. Approximately 75–98% of streamlines reached true positive ROIs for our predefined regions. By performing linear regression based on model (1) and calculating the 95% confidence interval for the intercept, we found that all nine confidence intervals included the origin. Therefore, we fit to the constrained linear model (2) to generate the regression lines in Fig. 4. The Pearson correlation coefficient r and corresponding p value for each case are shown in Table 1. Because the data span many orders of magnitude, the Pearson correlation coefficient (and corresponding p) is dominated by the largest values. For this reason, we also calculated the Spearman rank correlation coefficient, r_s, and corresponding p value, which depends only on the relative rank of connection strengths measured by the two methods. The values of r_s in Table 1 are quite low and do not reach statistical significance. To focus on the strongest connections, we chose the ten regions with the highest number of BDA-labeled fibers (i.e., number of BDA fibers>100) and calculated the Spearman rank correlation coefficient, (r_s)₁₀, and corresponding p value for this subset of the data.

Comparison across the three tractography schemes with d_w = 0 or 0.3 mm (distance of zero or one voxel) indicates that the linear correlations for the two probabilistic schemes are significantly lower than those of the deterministic scheme. When d_w = 0.6 mm, however, the linear correlations for these three schemes are quite similar and deterministic tractography gave no false negative results. In addition, under any fixed d_w, the FSL1 and FSL2 schemes have similar r value, but FSL2 has a smaller number of false negative ROIs than FSL1 does.

The influence of d_w on the three tractography schemes, shown in the rows of Fig. 4A–C and first nine rows of Table 1 are different. As d_w increases, the number of false negative ROIs decreases except for the FSL2 scheme. For the FSL1 and FSL2 schemes, as d_w increases, the r value increases significantly. For the DS scheme, r value increases less strongly than for the FSL1 and FSL2 schemes.

The proportional relation between the number of BDA-labeled somas and the number of streamline-terminals, shown in Fig. 4D and Table 1 row 10 is similar to the relationship between BDA-labeled fibers and streamlines in Fig. 4A and Table 1 row 1, except the number of false negative ROIs were more for the soma-terminals case. The proportional relationship between BDA-labeled soma and fibers, shown in Fig. 4E and bottom 2 rows in Table 1, is statistically significant.

Fig. 5 shows the backbones of the BDA-labeled, DS and FSL2 derived inter-regional connectivity. There are no edges connecting IR to iPF, cPF and cPVR shown in Fig. 5A but there are edges between IR to those regions in Fig. 5BC, which indicates iPF, cPF and cPVR are false positive regions detected by DS and FSL2 with d_w = 0.6 mm.

Correlation of Intra-ROI Connectivity Spatial Distribution

To visualize the spatial distribution of fibers (or somas) as well as streamlines (or streamline-terminals) from the ROI scale down to the voxel scale, the density distribution of BDA somas, BDA interface-crossing fibers, DTI streamline-terminals and DTI interface-crossing streamlines were rendered on the 3D WGM interface in DTI space (Fig. 6). Comparison between Fig. 6B and Fig. 6C–E indicated that distribution patterns of BDA and DTI fibers (DS and FSL schemes) are quite different. Table 2 lists the quantitative evaluation of the agreement: the Pearson’s correlation coefficient, r_p, of anatomical and DTI derived distributions in each individual ROI. Overall, the significant linear correlations for all the ROIs were lower than 0.5. The FSL2 scheme provided stronger linear correlations than FSL1 did, but comparable correlations to the DS scheme for most ROIs. The agreement in iPM is higher than the other projection regions due to higher sensitivity of DTI along this fiber pathway.

We also investigated whether deeper seed regions could help to achieve better distribution similarity by increasing the probability that streamlines could bypass crossing fibers directly beneath the WGM interface. Results for d_w = 0, 0.3 and 0.6 mm for each tractography scheme are shown in rows of Table 2. For the DS scheme, no significant trend in correlations was found as d_w increased (p>0.99 in one-way ANOVA). For the FSL1 and FSL2 schemes, correlations for ipsilateral ROIs did not have statistically significant differences (p>0.75 and p>0.72) although they seemed to have subtly increasing trends as d_w increased. No test was performed for contralateral ROIs, since the correlations were weak in those cases.

Sources of Discrepancies between Histological and DTI Connectivity

The challenges facing diffusion MRI tractography are well known, in principle. Limitations due to image artifacts and noise [26], [44], [45], partial volume averaging [46], and angular resolution [47] have been studied both experimentally and theoretically. Partial volume averaging was an obvious source of error in this study, as illustrated in Figures 7 and 8. Although BDA-labeled fibers connecting the injection region to the contralateral hemisphere run left/right in the plane of Figure 7, the principal diffusion direction (corresponding to most fibers in these voxels) is oriented anterior/posterior just under the injection region (Fig. 7B) and superior/inferior in deeper white matter (Fig. 7C). The disagreement between the BDA fibers bound for the contralateral hemisphere and the dominant fiber orientation prevents most DTI streamlines from reaching the midline and the other hemisphere. In fact, the strong anterior/posterior orientation of tensors immediately under the injection region is partially responsible for the strong (false positive) connectivity to the ipsilateral prefrontal cortex (Fig. 6) measured with both the DS and FSL2 schemes.

Probabilistic tracking methods based on the multi-tensor model are predicted to be better able to resolve crossing fibers, but the histology vs. FSL fiber correlations (Table 2) are not significantly higher than the histology vs. DS fiber correlations in this study. Figure 8AB shows that FSL2 is also limited by partial volume averaging in this data set. The FSL2 streamlines do not reach the WGM interface under the contralateral PM cortex. Figure 8C displays the two fiber orientations estimated by the FSL2 scheme in each white matter voxel immediately underneath the contralateral PM cortex. The vast majority of axons in these voxels travel parallel to the WGM interface and strongly bias the FSL-estimated fiber orientations also to be nearly parallel to the WGM interface. All six FSL cases (two ARD values, three seed region depths) were subject to similar limitations (although it is possible that the estimated orientations could be more accurate at higher SNR and/or b-value [48]). In summary, partial volume averaging of diffusion tensor data is likely a major source of error in our cortical connectivity measurements for both deterministic and probabilistic tracking algorithms.