Comparison of 3D orientation distribution functions measured with confocal microscopy and diffusion MRI.

MRI Acquisition

All animal procedures were approved by the Vanderbilt University Animal Care and Use Committee. Diffusion MRI experiments were performed on an adult squirrel monkey brain that had been perfusion fixed with physiological saline followed by 4% paraformaldehyde. Prior to fixation, the brain had undergone micro-electrode array recording experiments for an unrelated study. Due to this, there was slight atrophy in one hemisphere near the motor and pre-motor areas, resulting in assymetric hemispheres in both histology and MRI. However, this does not impact the results of our study, which is focused on characterizing white matter structure by determining the underlying fiber orientation distribution. The brain was then immersed in 4% paraformaldehyde for 3 weeks. The brain was transferred into a phosphate-buffered saline medium for 24 hours and scanned on a Varian 9.4 T, 21 cm bore magnet using a multi-shot multi-slice spin echo EPI sequence (TR = 6.7s; TE = 42ms; δ = 8ms; Δ = 27ms; max gradient strength = 30G/cm; voxel size = 400um isotropic; partial Fourier = .75; NEX = 5). A 30-direction diffusion-sampling scheme based on an electrostatic repulsion algorithm ( Jones et al., 1999 ) was used to acquire 30 diffusion-weighted images at a b-value of 3200 s/mm 2 , and 2 additional images were collected with b=0. This set of data was used for calculating diffusion tensors using a weighted linear least squares fit. Next, a 90-direction scheme was used to acquire diffusion weighted-images at a b-value of 6400 s/mm 2 , and 6 additional images at b=0. From this data set, the MRI-FOD was estimated using constrained spherical deconvolution with the damped Richardson-Lucy algorithm ( Dell'acqua et al., 2010 ) and fit to 8 th order spherical harmonic (SH) coefficients. MRI data processing was done using the high angular resolution diffusion imaging (HARDI) toolbox for MATLAB, available at http://neuroimagen.es/webs/hardi_tools/ .

Histological Procedures

After imaging, the brain was sectioned on a cryomicrotome at a thickness of 80um in the coronal plane and mounted on glass slides. Using a Canon EOS20D (Lake Success, NY, USA) digital camera with a zoom lens of 70-300 mm, the tissue block was digitally photographed prior to cutting every other section, resulting in a 3D “block-face” volume with a through-plane resolution of 160um. The tissue sections were mounted on glass slides and stained following the procedures outlined in ( Budde and Frank, 2012 ). Briefly, tissue sections were rinsed in PBS and dehydrated through graded ethanol solutions. The fluorescent lipophilic dye, “DiI”, (1,1’-dioctadecyl-3,3,3’3’-tetramethylindocarbocyanine percholarate) in 100% ethanol (.25mg/mL) was rinsed over sections for 1 minute. The stained sections were then rehydrated through graded ethanol solutions, and coverslipped with Fluoromount-G mounting medium.

Confocal Acquisition

All histological data were collected using an LSM 710 inverted confocal microscope (Carl Zeiss, Inc. Thornwood, NY. USA). For all selected tissue slices, confocal acquisition consists of two protocols: [1] creating a 2D montage of the entire tissue and [2] constructing a 3D high-resolution image in a selected region of interest. The 2D montage ( Fig. 1A ) consists of approximately 600-900 individual tiles acquired using a 10x oil objective at a resolution of 0.80μm 2 , which are stitched together using Zeiss software, ZEN 2010. Acquisition for a single slice takes approximately 30 minutes. To correct for image inhomogeneity and tiling effects in the image, we found it useful to increase the zoom feature to 1.5x or higher at the expense of collecting more tiles. This 2D montage is used for image registration, and for localizing the 3D high-resolution region of interest. Prior to 3D z-stack acquisition, two steps are performed. First, tissue thickness in the z-dimension is determined by adjusting the focal plane depth to determine where fluorescence begins and ends. This thickness is used to correct orientation estimates for tissue shrinkage (see Histological FOD below). Second, it is necessary to increase the laser output as deeper layers are imaged due to the increases in light scatter and absorption at greater tissue depths (see Confocal Pre-Processing ). The laser power is adjusted for approximately 5 different depths ranging from the coverslip to the end of the tissue, at each step ensuring that the image intensity range will cover the full 8-bit depth from 0-255 units. The LSM 710 interpolates the laser output between depths. The 3D z-stack ( Fig. 1B ) is then collected using a 63x oil objective at a nominal resolution of 0.18μm×0.18μm×0.42μm. Typical acquisition time to acquire the entire section thickness with an in-plane field of view of 1.6mm × 1.6mm (equivalent to 16 MRI voxels) is approximately 8 hours. The through-plane resolution is the “optimal” slice-thickness, calculated from the LSM710 software based on a 1.0 Airy unit pinhole diameter and an excitation wavelength of 543nm. Stitching, again, is performed using ZEN 2010 software to create a single 3D z-stack. Finally, all confocal data are converted from the LSM file format to TIFF images and imported into MATLAB for further processing.

Confocal Pre-processing

The aim here is to extract the histological-FODs from the 3D z-stacks in areas equivalent to the size of an MR voxel. To do this, we use structure tensor analysis to obtain an orientation estimate for every pixel in the 3D z-stack that is occupied by a fiber. Prior to structure tensor analysis, four sources of anisotropy inherent to confocal microscopy must be accounted for ( Fig. 1C ). Three corrections are performed directly on the confocal z-stack prior to structure tensor analysis, and the final correction performed post-analysis. The first is an attenuation of the image intensity as a function of tissue depth. This effect is caused by light scatter and absorption which decreases the intensity of excitation light penetrating to the deeper layers of the tissue, and consequently, the fluorescence of these layers. Because structure tensor analysis is based on image intensity gradients, this artifact could result in a bias in fiber orientation estimates ( Khan et al., 2015 ). The attenuation correction is performed in the Confocal Acquisition stage described above. Increasing the laser power for deeper layers generates a z-profile that has a relatively constant mean intensity in each x-y plane containing fibers. The second source of anisotropy arises from the confocal microscope’s point spread function (PSF). The PSF is the 3D diffraction pattern resulting from the systems response to an infinitely small point source of light. This diffraction pattern is known to be nearly three times wider through-plane than in-plane ( Pawley and Masters, 1996 ), leading to anisotropic blurring of the image; in-plane structures will be better resolved than those oriented through-plane. To deblur the confocal data, we use the iterative Lucy-Richardson algorithm ( Biggs and Andrews, 1997 ) and a computed theoretical model of the confocal microscope’s PSF ( Pawley and Masters, 1996 ). This model takes into account various confocal parameters including the numerical aperture, refractive index, wavelength of light, and the acquired image resolution. The Lucy-Richardson deconvolution algorithm is a maximum-likelihood approach to find the statistically most likely image, given the blurred image and assuming Poisson noise ( Biggs, 2010 ; Biggs and Andrews, 1997 ), which is an appropriate noise model of the photon-counting process of confocal imaging ( Pawley and Masters, 1996 ). The final pre-processing step is to correct for the anisotropic acquisition resolution. This ensures that fibers oriented laterally in the image will contain an equivalent number of pixels per length as fibers oriented axially. Interpolation to isotropic resolution is accomplished using cubic interpolation.

Structure Tensor Analysis

The structure tensor was introduced in the late 1980’s for point and edge detection ( Bigun and Granlund, 1987 ; Harris and Stephens, 1988 ), and has since become popular in image processing and computer vision, with applications including texture analysis and materials science ( Axelsson, 2008 ; Krause et al., 2010 ). This analysis technique is applied to our entire 3D confocal image, f(x,y,z). The structure tensor ( Köthe, 2003 ) is based on the gradient of f: [1] ∇ f σ = ( f x , f y , f z ) T which is calculated with Gaussian derivative filters: [2] f x = g x , σ ∗ f , f y = g y , σ ∗ f , f z = g z , σ ∗ f where denotes the convolution operation and g x,σ , g y,σ , and g z,σ are the spatial derivatives in the x, y, and z-direction, respectively, of a 3D Gaussian with standard deviation σ: [3] g σ ( x , y , z ) = 1 ( 2 π σ 2 ) 3 e − ( x 2 + y 2 + z 2 ) 2 σ 2 For illustration purposes, we show this step as performed on a simulated cylindrical fiber ( Fig. 2A ), representative of the neuronal structures seen in the 3D confocal image (similar illustrations appear in ( Khan et al., 2015 ) and ( Arseneau, 2006 )). Ideally, the image gradients are orthogonal to the fibers at all points ( Fig. 2B ). Next, an object known as the gradient square tensor, is calculated for each point in the image by taking the dyadic product of the gradient vector with itself: [4] G S T ( x , y , z ) σ = ∇ f σ ∇ f σ T = ( f x 2 f x f y f x f z f x f y f y 2 f y f z f x f z f y f z f z 2 ) Each tensor element is averaged over a local neighborhood to create the pixel-wise structure tensor. For spatial averaging, we choose a 3D Gaussian filter with standard deviation ρ: [5] S T ρ ( ∇ f σ ) = g ρ ∗ ( ∇ f σ ∇ f σ T ) This results in a 3-by-3 symmetric, semi-positive definite, rank-two tensor. Much like the diffusion tensor, this matrix will have three positive eigenvalues, and can be visualized as an ellipsoid ( Fig. 2C ). In DTI, one is typically interested in the largest eigenvalue and eigenvector, which points in the direction of greatest diffusion, and is usually assumed to be parallel to the primary structure orientation in the MR voxel. However, in structure tensor analysis, the image intensity gradients are strongest perpendicular to the fibers, which means the largest two eigenvectors will also be perpendicular to the fiber bundles. Hence, we make the assumption that the direction of minimal intensity variation is parallel to the fiber orientation at each pixel, a direction given by the eigenvector corresponding to the smallest eigenvalue. The certainty in estimated fiber orientation can be described by the Westin-measure ( Westin et al., 2002 ) defining how planar the structure tensor is: [6] C p = λ 2 − λ 3 λ 1 where λ 1 , λ 2 , and λ 3 are the primary, secondary, and tertiary eigenvalues of the structure tensor. This value varies from 0 to 1 and will be large in areas, like that depicted in Fig. 2C , where the first two eigenvalues are much larger than the third. This measurement is used to threshold the confocal image, so voxels with low certainties are not included in the final orientation distribution. For the results presented in this paper, the spatial derivatives were calculated using a Gaussian with standard deviation σ = 1μm, and spatial averaging performed using a Gaussian with standard deviation ρ = 2.5μm (these values were chosen based on comparisons to distributions of manually traced fibers – see section 3.1).

Histological-FOD

After a fiber orientation has been extracted for all pixels in the image ( Fig. 1C ) one final correction for anisotropy must be performed. It is known that tissue samples may shrink due to processing, sectioning, and staining ( Eltoum et al., 2001 ; Woods AE, 1994 ). These effects are mainly a result of fixation and dehydration in alcohol solutions during the staining procedure ( Wehrl et al., 2015 ; Williams et al., 1997 ). We use the thickness measurement before acquisition of each 3D z-stack to perform a geometric correction to the orientation estimate for every pixel in the image by assuming linear shrinkage in the through-plane (z) direction. Once all estimated vectors have been appropriately re-oriented to account for tissue shrinkage, the results are thresholded using both image intensity and the certainty value. This yields an orientation estimate for every pixel in our z-stack that is occupied by a fiber. A histogram representing the histological-FOD is then created as a function of polar and azimuthal angle, where the orientation estimates are placed into bins that cover constant solid angles over a sphere. This FOD is fit to high order (20) SH coefficients, and throughout this paper is displayed as a three dimensional glyph ( Fig. 1D ) in the same way that the MRI-FOD’s are typically displayed.

Image Registration

In order to make a quantitative comparison of the histological-FOD and the MRI-FOD, the data must be aligned and oriented appropriately. A multi-step registration procedure ( Choe et al., 2011 ) was used to align histology to MRI data. The first step is registration of the 2D confocal montage to the corresponding block face image using mutual information based 2D linear registration followed by 2D nonlinear registration using the adaptive bases algorithm (ABA) ( Rohde et al., 2003 ). Next, all block face photographs were assembled into a 3D block volume, which is registered to the MRI b=0 image using a 3D affine transformation followed by 3D nonlinear registration with ABA. Given the location of the 3D z-stack in the 2D confocal montage, we can use the combined deformation fields to determine the MRI signal from the same tissue volume. The MRI signal of interest is analyzed in MRI native space. As described above, we derive the tensor using a WLLS fit, and estimate the MRI-FOD using constrained spherical deconvolution. The final step is to transform the diffusion tensor and MRI-FOD to histological space to facilitate comparisons with the histological FOD. For the tensor, we apply the preservation of principal directions (PPD) strategy ( Alexander et al., 2001 ) twice: once to transform the tensor from MRI-space to block-space, and again to transform to histological space. For the MRI-FOD, we choose the approach developed in Hong et al. ( Hong et al., 2009 ). This method takes into account rotation, scaling, and shearing effects of the spatial transformations, and can be applied to any orientation distribution on a sphere. After these corrections, both the histological-FOD and the MRI-FOD are in histological space, and quantitative analysis can be performed.

Histological Procedures

After imaging, the brain was sectioned on a cryomicrotome at a thickness of 80um in the coronal plane and mounted on glass slides. Using a Canon EOS20D (Lake Success, NY, USA) digital camera with a zoom lens of 70-300 mm, the tissue block was digitally photographed prior to cutting every other section, resulting in a 3D “block-face” volume with a through-plane resolution of 160um. The tissue sections were mounted on glass slides and stained following the procedures outlined in ( Budde and Frank, 2012 ). Briefly, tissue sections were rinsed in PBS and dehydrated through graded ethanol solutions. The fluorescent lipophilic dye, “DiI”, (1,1’-dioctadecyl-3,3,3’3’-tetramethylindocarbocyanine percholarate) in 100% ethanol (.25mg/mL) was rinsed over sections for 1 minute. The stained sections were then rehydrated through graded ethanol solutions, and coverslipped with Fluoromount-G mounting medium.