Three-dimensional fiber orientation mapping of the human brain at micrometer resolution

TL;DR

This paper introduces a new method to map 3D fiber orientation in the human brain at micrometer resolution, enabling detailed study of brain connectivity.

Contribution

A novel computational method using PS-OCT to quantify volumetric 3D fiber orientation in full angular space is introduced.

Findings

The method measures polarization contrasts from two illumination angles to determine 3D optic axis orientation and true birefringence.

3D fiber orientation maps of entire coronal cerebrum sections and brainstem were generated with 10 μm in-plane resolution.

The approach reveals unprecedented details of fiber configurations in the human brain.

Abstract

The accurate measurement of three-dimensional (3D) fiber orientation in the brain is crucial for reconstructing fiber pathways and studying their involvement in neurological diseases. Comprehensive reconstruction of axonal tracts and small fascicles requires high-resolution technology beyond the ability of current in vivo imaging (e.g. diffusion magnetic resonance imaging). Optical imaging methods such as polarization-sensitive optical coherence tomography (PS-OCT) and polarization microscopy can quantify fiber orientation at micrometer resolution but have been limited to two-dimensional in-plane orientation or thin slices, preventing the comprehensive study of connectivity in 3D. In this work we present a novel method to quantify volumetric 3D orientation in full angular space with PS-OCT. We measure the polarization contrasts of the brain sample from two illumination angles of 0 and 15 degrees and apply a computational method that yields the 3D optic axis orientation and true birefringence. We further present 3D fiber orientation maps of entire coronal cerebrum sections and brainstem with 10 μm in-plane resolution, revealing unprecedented details of fiber configurations. We envision that our method will open a promising avenue towards large-scale 3D fiber axis mapping in the human brain as well as other complex fibrous tissues at microscopic level.

Figures

a](#F1)). To implement the tilted illumination, we either mounted the sample on a rotation stage or incorporated a rotation stage into the PS-OCT scan head. The physical effect of tilting the beam or the sample is the same ^[20](#R20)^. For the normal illumination, one volumetric acquisition covered a field of view (FOV) of 3.5 mm × 3.5 mm. The FOV was 3.5 mm × 2 mm for the 15 deg tilted illumination. Automatic tile-scan was used to cover the whole sample, with a 20% overlap between tiles in each direction of the raster scan. For volumetric imaging of large-scale samples, a 150 μm thick slice

a](#F2)), respectively. The long axis of the fibers was positioned along the x-axis of the laboratory frame, creating a 0 deg in-plane orientation. The obtained PS-OCT measurements were obtained from the normal and a 15 deg tilted illuminations in each experimental condition, and 3D axis orientation was estimated using the proposed method above. As expected, our optimization yielded the same estimated in-plane orientations of 0 deg ([Fig. 2b](#F2)), whereas distinct through-plane orientations with a difference of approximately 30 deg between conditions 1 and 2, and 1 and 3 ([Fig. 2c](#F2)). We

a](#F3)). A linear fitting of the estimated through-plane orientation with respect to the preset tilting angle revealed a slope of 0.97, very close to the theoretical value of 1. As an example, we showed two through-plane orientation maps reconstructed as the stage was positioned at 0 deg and 50 deg, in [Fig. 3b](#F3). The mean through-plane orientation of the fibers was measured to be 18.6 deg and 73.6 deg, respectively. In the second experiment, we imaged another region of the sample by tilting the stage at 10 deg, 30 deg and 50 deg angles and estimated the corresponding through-plane orient

a](#F4)) when the sample was oriented at −90, −60, −30 and 0 deg in the XY plane. Meanwhile, the through-plane orientation estimate was consistent among the different in-plane orientation setups ([Fig. 4b](#F4)). To evaluate our method quantitatively, we analyzed the fiber orientation distributions (FODs) of a 3 mm × 3 mm ROI as shown in [Fig. 4a](#F4) and [Fig. 4b](#F4). The mean estimated in-plane orientation was −82.8, −57.2, −24.8 and 8.5 deg ([Fig. 4c](#F4)), accurately following the in-plane rotation. Note that there was a small angular offset (average of 5.9 deg) between the long axis o

a](#F6)), where multiple colors in the in-plane orientation map ([Fig. 6b](#F6)) indicated the confluence of fibers bundles from different cortical regions into the superficial white matter. This low apparent birefringence is not due to low myelin content. Instead, it is caused by a group of fibers going through plane, proved by the high through-plane angle at the same locations ([Fig. 6d](#F6)). These dark spots are diminished in the true birefringence map ([Fig. 6c](#F6)). It is also intriguing to find that the splenium of the corpus callosum exhibits much higher true birefringence compared

f](#F7)). To demonstrate the advantages of measuring 3D axis orientation by PS-OCT, we performed volumetric imaging of a slab of a human brainstem at the level of the medulla oblongata ([Fig. 7a](#F7)-[7d](#F7)), following dMRI scanning of the whole brainstem and cerebellum ([Fig. 7e](#F7)). The tissue went through PS-OCT imaging with both normal and 15 deg tilted illumination and sectioning as previously reported for large volumetric reconstruction ^[17](#R17)^. The total reconstructed volume was approximately 30×20×3 mm^3^. The brainstem was positioned in an axial plane, i.e., left-right (re

a](#F8)) present the sophisticated fiber network in the human brainstem. The geometry of the fiber tracts was very well preserved without the errors caused by inter-slice alignment and the slice damage typically seen in traditional histological images including PLI. To leverage the volumetric 3D axis orientation for connectivity studies, we tracked the trajectory of the fibers using off-the-shelf dMRI tractography toolkit in selective fiber bundles ([Fig. 8b](#F8)), including the solitary tract (sol), the amiculum of the olive (ami, ROI 1), the pyramidal tract (py, ROI 2), and the spinal trige

Full text

Introduction

The fiber tracts constituting the white matter in the human brain form an intricate network of connections that support neuronal communication and brain function. Deciphering the structural connectivity of the brain is crucial to understand normal behavior and what goes wrong in neurological disorders. Diffusion weighted magnetic resonance imaging (dMRI) detects the preferential direction of water diffusion, which is parallel to the fiber bundles; it allows fiber pathways to be imaged in vivo and has been widely used in clinical settings such as neurosurgical planning and targeting for deep brain stimulation ^1^. Despite tremendous advancements over the past decade ^2,3^, the mm-resolution of dMRI limits its ability to resolve complex fiber configurations, particularly in the regions where multiple fiber bundles from different directions intersect ^4^. Microscopic methods are necessary to unveil the ground truth fiber architectures and to guide the validation of new dMRI techniques ^4,5^. Optical methods have used various contrast mechanisms to image myelinated axons at micrometer resolution, such as fluorescent labeling ^6^, third harmonic generation ^7^ and anti-stokes Raman scattering ^8^. However, these techniques do not provide direct measures of fiber orientation that can be used to reconstruct white-matter tracts and study structural connectivity in the brain. Instead, fiber orientation is computationally derived from the acquired images, using techniques such as structure tensor analysis, which has been applied to histological staining ^9,10^, fluorescence microscopy ^11^, and label-free optical coherence tomography (OCT) images ^12^.

Imaging techniques using light polarization take advantage of the birefringent nature of the myelin sheath surrounding axons to generate quantitative and label-free contrasts for myelinated fibers ^13^. Birefringence originates in the differential refraction between two orthogonal axes as light travels through the anisotropic materials. The fiber geometry causes optical anisotropy with negative birefringence and an optic axis orientation aligned with the physical fiber direction. In recent years, methods that rely on birefringence, such as polarized light imaging (PLI) ^14^ and automatic serial sectioning polarization-sensitive OCT (PS-OCT) ^15,16^, have been applied to map microscopic fiber orientation in large-scale brain samples. PLI uses multiple oblique views to derive the 3D fiber orientations in human brain sections, after cutting of the sample into 50–100 μm thick slices ^17^. Following fiber tracts across sections requires inter-slice registration with microscopic accuracy, which is an extremely challenging task ^17,18^. In serial sectioning PS-OCT, the blockface of the sample is imaged prior to sectioning, thus avoiding section-to-section distortions. The conventional optic axis orientation measurement of PS-OCT from one illumination angle quantifies an in-plane fiber orientation ( $[eqn]$ , Fig 1a), which is a projection of the 3D fiber axis onto a plane orthogonal to the beam. Information on the through-plane fiber orientation ( $[eqn]$ , Fig 1a) is not captured. However, fiber tracts in the brain run in all different directions in 3D; therefore, it is impossible to track them without the full 3D axis orientation. Another drawback of traditional PS-OCT is that it measures an apparent birefringence that has a dependence on the through-plane axis. Without removing the dependency, the measured birefringence does not truly reflect the structural anisotropy, which is an important indicator of myelin integrity in neurological disorders ^19^. One approach to obtain the through-plane orientation is to use multiple illumination angles. The in-plane optic axis orientations obtained by different illumination angles allow the calculation of the through-plane orientation via a geometric derivation. Our previous study suggested a minimum of three beams were required to accurately recover the 3D orientation of muscle fibers ^20^. Similarly, various illumination configurations have been proposed in PS-OCT and validated in simple birefringent samples of tendons and muscles ^20–23^. Those configurations have been using widely spaced imaging directions, and either excluding or assuming optic axes co-planar to the illumination beam directions, therefore limiting the range of the measurable through-plane axis. Furthermore, 3D fiber orientation measurements using PS-OCT in the human brain have not been reported.

In this work, we present the estimation of 3D axis orientation in human brain tissue that allows a complete mapping of the fiber axis within the 360-degree spherical space with a novel procedure and a computational framework. The estimation is enabled by only two illumination beams separated by a small angle of 15 deg. This approach formulates a birefringence vector composed of both optic axis orientation and true birefringence, which is estimated by minimizing the error from the two measured birefringence values and in-plane optic axis orientations. We conduct a systematic validation on the performance of our approach and show the 3D fiber orientation maps in various large human brain samples, including coronal sections of a whole hemisphere and a medulla sample.

Results

Reconstruction of 3D axis orientation via two illumination angles

Traditional PS-OCT only measures the apparent birefringence and in-plane optic axis orientation, whereas the through-plane component of the 3D axis is lost. Recent PS-OCT work has reported that the true 3D optic axis can be recovered by the measurements of apparent birefringence ^21^ or optic axis orientations ^20^ at variable illumination angles. Both the measured apparent birefringence and the optic axis have a dependence on the illumination angle and the measurements from different incident angles are linked by the through-plane orientation of the fibers. Geometrically, the 3D optic axis lies along the intersection of multiple planes, each of which is formed by an incident beam and its corresponding optic axis orientation. However, it is necessary to also account for the effective birefringence values, in order to manage the case when the optic axis is co-planar to the illumination beams, and hence the two planes coincide. Our forward model uses both optic axis and birefringence to capture these physical relationships.

To estimate the 3D optic axis via two illuminations, we first define the true birefringence vector $[eqn]$ , where $[eqn]$ is the true birefringence, $[eqn]$ is the in-plane orientation, and $[eqn]$ is the through-plane orientation (see the insert in Fig. 1a). We imaged the sample with normal illumination and illumination tilted by $[eqn]$ deg about the y-axis with our automatic serial sectional PS-OCT system (Fig. 1a and see Methods for details). Under the normal illumination, we obtained the apparent birefringence $[eqn]$ and orientation measurement $[eqn]$ (Fig. 1b); whereas under the tilted illumination, we measured the apparent birefringence $[eqn]$ and orientation measurement $[eqn]$ (Fig. 1c). By defining estimated birefringence vectors $[eqn]$ the two illuminations, we formulated the optimization problem in Eq. 1 as finding the birefringence vector $[eqn]$ that minimizes the sum of the mean square errors between the estimated and the measured apparent birefringence vectors under the two illuminations (see Methods for details):

$$[eqn]$$

After estimating the true birefringence vector $[eqn]$ , we then obtained the true birefringence $[eqn]$ and through plane orientation $[eqn]$ (Fig. 1d). The data acquisition and image reconstruction procedures are demonstrated in Fig. 1, using a homogeneous fibrous sample in the human brain.

Recovering full-range through-plane axis orientation with high accuracy

We used a section of a human corpus callosum in which the fibers run in a mostly uniform direction to validate the accuracy of the 3D fiber orientation estimation. To create different through-plane angles, we positioned the sample on a horizontal surface (0 deg) and two triangular wedges made of agarose gel (−30 and 30 deg, Fig. 2a), respectively. The long axis of the fibers was positioned along the x-axis of the laboratory frame, creating a 0 deg in-plane orientation. The obtained PS-OCT measurements were obtained from the normal and a 15 deg tilted illuminations in each experimental condition, and 3D axis orientation was estimated using the proposed method above. As expected, our optimization yielded the same estimated in-plane orientations of 0 deg (Fig. 2b), whereas distinct through-plane orientations with a difference of approximately 30 deg between conditions 1 and 2, and 1 and 3 (Fig. 2c). We further examined the through-plane orientation in a smaller region of interest (ROI) and plotted the fiber orientation distributions (FOD). We found the mean of the through-plane FODs to be 9.2 deg on 0 deg surface (no wedge), −26.9 deg on −30 deg wedge and 37.7 deg on 30 deg wedge (Fig. 2d). The difference between the three estimates agreed with the known difference between the inclination angles of the wedges. The estimates showed a positive offset consistently with respect to the preset wedge angles, indicating the presence of the fibers oriented at a small through-plane angle at the 0 deg wedge setting (about 6.7 deg on average). It is important to underline that our method was able to distinguish the positive and negative sign of the through-plane orientations (green vs. magenta wedges in Fig. 2), which is critical in connectivity analysis as 30 deg and −30 deg fibers run into completely different directions.

We then used another block of the corpus callosum sample to characterize the range of our through-plane orientation measurement. As the white matter fiber tracts in the brain could run in all directions in spherical space, it is important for our method to capture the full +/−90 deg range of through-plane angles. We performed two experiments on different regions of the sample. The sample was positioned on a tilting stage to artificially vary the through-plane angle. In the first experiment, the stage was tilted counterclockwise from 0 deg to 50 deg in 10 deg steps. Our estimated through-plane orientation successfully followed the preset tilting angles (Fig. 3a). A linear fitting of the estimated through-plane orientation with respect to the preset tilting angle revealed a slope of 0.97, very close to the theoretical value of 1. As an example, we showed two through-plane orientation maps reconstructed as the stage was positioned at 0 deg and 50 deg, in Fig. 3b. The mean through-plane orientation of the fibers was measured to be 18.6 deg and 73.6 deg, respectively. In the second experiment, we imaged another region of the sample by tilting the stage at 10 deg, 30 deg and 50 deg angles and estimated the corresponding through-plane orientation. The slope of a linear fitting between the estimation and the pre-set tilting angle was 1.21, and the estimated through-plane orientation at 0 deg tilting was 26.8 deg. The slight deviation of the fitted slope from 1 was possibly attributed to the increased noise at the large tilting angle, as the estimated through-plane orientation was close to 90 deg. Nevertheless, we demonstrated that our method can measure a full range of through-plane orientation angles. Taken together with Fig. 2, which shows a successful estimation of the sign of the through-plane orientation, these results suggest that our method is capable of recovering the full spherical range of the 3D fiber orientation in the human brain.

Robust estimation under variety of in-plane orientations

A previous study showed that the accuracy of the through-plane axis was highly dependent on the in-plane orientation; as a result, a minimum of three illumination angles was required to recover the through-plane axis reliably ^20^. The method that we propose here uses an optimization algorithm that eliminates the need for the third illumination angle, while retaining a viable through-plane orientation estimation. To validate the robustness of the estimation under a variety of in-plane orientation scenarios, we mounted a corpus callosum sample on a horizontal surface while rotating it counterclockwise from 0 to 90 deg in the XY plane by 30 deg steps. We expected that through-plane orientation should remain constant around 0 deg while the in-plane orientation should increase by 90 deg through rotation. We obtained the measurements via the two illumination angles, normal and tilted by 15 deg (Fig. 4). Our method estimated the in-plane orientations accurately (Fig. 4a) when the sample was oriented at −90, −60, −30 and 0 deg in the XY plane. Meanwhile, the through-plane orientation estimate was consistent among the different in-plane orientation setups (Fig. 4b). To evaluate our method quantitatively, we analyzed the fiber orientation distributions (FODs) of a 3 mm × 3 mm ROI as shown in Fig. 4a and Fig. 4b. The mean estimated in-plane orientation was −82.8, −57.2, −24.8 and 8.5 deg (Fig. 4c), accurately following the in-plane rotation. Note that there was a small angular offset (average of 5.9 deg) between the long axis of the fibers and the x-axis of the laboratory frame. The through-plane orientation estimation remained consistently close to 0 deg (2.8, −1.2, −2.4 and −1.6 deg) across the different in-plane orientation settings (Fig. 4d). This suggests that our method with only two illumination incidences can measure the 3D axis of all fibers robustly regardless of their orientation in the XY plane.

We demonstrated the advantage of a strategy using an xy-axis swap at 45 deg for through-plane axis estimation (Fig. 4, see details in Methods), by comparing to the results without xy-swap (using Eq. 1, 2, 3 for all in-plane orientations). In the latter case, both the in-plane (Fig. S2a) and through-plane (Fig. S2b) orientation showed high errors at −90 deg rotation (see Fig. S1 for details). We also compared our optimization result against an earlier method using a geometric derivation (Fig. S2c, see details in reference ^16^). Using only one tilted beam about the y-axis, the geometric derivation method resulted in errors in through-plane estimation as the in-plane angle reached −60 deg. The error continued to grow as the in-plane angle increased to −90 deg.

Estimating true birefringence of myelin content

Birefringence is a measure of myelin concentration in the brain and has important indications in neurodegenerative disorders involving myelin degradation. However, it is known that traditional PS-OCT measures an apparent birefringence that depends on the optic axis orientation of the tissue: the larger the through-plane orientation $[eqn]$ , the smaller the apparent birefringence $[eqn]$ ^14^. This is in contrast to true birefringence $[eqn]$ , which depends only on structural anisotropy and not axis orientation. Estimation of true birefringence requires removing the effect of the through-plane angle. Using our computational framework, we also investigated the estimation of true birefringence in the corpus callosum sample. We computed the mean of the apparent birefringence and true birefringence within the ROIs from the three inclination setups as shown in Fig. 2a and 2b. As expected, the apparent birefringence was reduced in the ±30 deg inclination setups (4.9 × 10^−4^ and 4.1 × 10^−4^) compared to the 0 deg setup (5.8 × 10^−4^), showing the dependency on the through-plane angle (Fig. 5). By jointly estimating the through-plane angle and true birefringence, we found that the latter remained consistent across the three setups, regardless of the through-plane angle (mean value: 6.2 × 10^−4^), thus providing an unbiased measure of the underlying myelin content in this corpus callosum sample.

Mapping 3D fiber orientation of a full coronal section in the human brain

To investigate the feasibility of 3D orientation measurements in a large-scale sample, we imaged a coronal section of a human brain hemisphere, with a surface area of ~6 cm × 9 cm. We first acquired dMRI data of the sample followed by PS-OCT. For PS-OCT imaging, the tissue surface was flat-faced, and imaged with both normal and 15 deg tilted illuminations to recover the 3D orientations. We generated two maps of apparent birefringence and in-plane orientation obtained with traditional single-illumination PS-OCT (Fig. 6a-6b), and three additional maps acquired with our method, including true birefringence, through-plane orientation and 3D axis orientation (Fig. 6c-6e). As a comparison, we also show dMRI tractography in the same plane as the PS-OCT images (Fig. 6f).

Overall, the 3D orientation map of PS-OCT (Fig. 6e) is consistent with dMRI (Fig. 6f) and the known neuroanatomy of white matter architecture in this region (Fig. 6) ^24^. However, the much higher resolution of PS-OCT can be easily appreciated in the images. In the 3D orientation map, medial-lateral/right-left (red) and inferior-superior (blue) directions are in-plane and anterior-posterior (green) is through-plane. The splenium of the corpus callosum is situated on the superiomedial border of the lateral ventricle and contains highly birefringent fibers oriented in the anterior-posterior (green) and medial-lateral (red) directions. Lateral to the lateral ventricle, three white matter tracts could be clearly delineated. The tapetum (Tp) borders the lateral edge of the lateral ventricle and contains fibers coursing inferior-superior (blue). The sagittal stratum (SSt) sits lateral to the Tp and contained fibers running primarily anterior-posterior (green). Lateral to the SSt, the inferior longitudinal fasciculus (ILF) consists of mostly inferior-superior oriented fibers (blue). Throughout the superficial white matter, short association fibers (e.g., U-fibers), are seen wrapping around the sulci and extending down the gyri into gray matter. Gyral fibers show predominantly in-plane orientations (red, blue), while small bundles of through plane fibers (green) are consistently observed at the base of sulci. The high-resolution 3D fiber orientation map of PS-OCT depicts a diverse, heterogeneous organization of fiber bundles in deeper white matter. The smooth transition of the colored fiber orientation along tracts indicates realistic continuity that one might expect in the densely packed fascicular highways of the deep white matter.

Compared to traditional, apparent birefringence, true birefringence revealed the unbiased myelin distribution in the brain by removing the effect of the through-plane angle. In the apparent birefringence maps, dark spots were consistently present in fibers right beneath sulci (white arrows in Fig. 6a), where multiple colors in the in-plane orientation map (Fig. 6b) indicated the confluence of fibers bundles from different cortical regions into the superficial white matter. This low apparent birefringence is not due to low myelin content. Instead, it is caused by a group of fibers going through plane, proved by the high through-plane angle at the same locations (Fig. 6d). These dark spots are diminished in the true birefringence map (Fig. 6c). It is also intriguing to find that the splenium of the corpus callosum exhibits much higher true birefringence compared to the rest of the white matter in this section (red arrows in Fig. 6c), as effect that is confounded in the apparent birefringence map by the high through-plane angles.

Volumetric 3D orientation and tractography of a human brainstem

The micrometer resolution PS-OCT reveals complex configurations of small fiber tracts at a voxel size 5,000–1,000,000 times smaller than the dMRI. Particularly in the brainstem, compact fiber tracts that connect small nuclei run in all directions, making it extraordinarily challenging to resolve small fibers (Fig. 7f). To demonstrate the advantages of measuring 3D axis orientation by PS-OCT, we performed volumetric imaging of a slab of a human brainstem at the level of the medulla oblongata (Fig. 7a-7d), following dMRI scanning of the whole brainstem and cerebellum (Fig. 7e). The tissue went through PS-OCT imaging with both normal and 15 deg tilted illumination and sectioning as previously reported for large volumetric reconstruction ^17^. The total reconstructed volume was approximately 30×20×3 mm^3^. The brainstem was positioned in an axial plane, i.e., left-right (red), and anterior-posterior (green) being the x and y axis of our laboratory frame, respectively. It is well known that there are large numbers of fibers running through the axial plane in the brainstem, as well as substantial in-plane fibers connecting the small nuclei. The true birefringence map indicated a high myelin content in the human medulla with the presence of various nuclei. Using the maps generated by PS-OCT, we annotated some of the major structures on the high-resolution images based on the Paxinos Atlas (Fig. 7f) ^25^ and inspected the orientation of the fiber tracts (Fig. 7b-7d). For example, the direction of internal arcuate fibers seen in the reticular formation (RF) running from the cuneate (Cu) and gracile nuclei (not shown in Fig. 7) to the contralateral medial lemniscus (ml), is consistent with the mostly in plane fibers in the RF and ml. The amiculum of the olive (ami) is a capsule of myelinated fibers surrounding the inferior olivary (IO) running both in-plane and through-plane. The spinal trigeminal tract (sp5), as well as the inferior cerebellar peduncle (icp), run inferior to superior, connected with other structures via numerous input and output fibers throughout their length in the axial plane. Finally, two other main fiber tracts, the pyramidal (pyr) and solitary tract (sol) run inferior to superior and thus appear mostly blue in the color rendering of the 3D fiber orientation.

We demonstrated the advantage of volumetric imaging with our method by reconstructing 19 consecutive block-face sections with a 10 μm in-plane resolution and 150 μm slice thickness. The 2.85 mm volumetric rendering of true birefringence (Supplementary Video 1), in-plane orientation (Supplementary Video 2), through-plane orientation Supplementary Video 3), and 3D axis orientation (Fig. 8a) present the sophisticated fiber network in the human brainstem. The geometry of the fiber tracts was very well preserved without the errors caused by inter-slice alignment and the slice damage typically seen in traditional histological images including PLI. To leverage the volumetric 3D axis orientation for connectivity studies, we tracked the trajectory of the fibers using off-the-shelf dMRI tractography toolkit in selective fiber bundles (Fig. 8b), including the solitary tract (sol), the amiculum of the olive (ami, ROI 1), the pyramidal tract (py, ROI 2), and the spinal trigeminal tract (sp5, ROI 3). We then visualized the zoomed tracts in the latter three ROIs. The ami forms a circumferential layer surrounding the inferior olivary nucleus (ION) while intervening with the ION from different directions. It is interesting to see that complex fiber configurations occur throughout white matter regions. For example, the fibers in py and sp5 form a twisting and inter-weaving pattern in small bundles while they travel through the superior-inferior axis. Moreover, the bundles in py and sp5 reveal differed inter-weaving angles that are smaller in py than sp5. Those twisting patterns might facilitate a large number of axonal groups reaching their widely spread terminals through these compact white matter pathways and forming a certain topology in their trajectories. The microscopic details revealed by our technique can provide unprecedented information to understand the complicated fiber connections in the brainstem that play a central role in functions such as modulating wakefulness in human consciousness ^26^.

Discussion

Directly imaging connective pathways in the human brain is crucial for neurosurgical planning and intraoperative guidance. In this work, we presented a new approach for estimating 3D fiber orientation and true birefringence of human brain samples using PS-OCT, enabling investigation of connectivity and myelin content at the microscale. Traditional PS-OCT only measures a projection of 3D fiber orientation and an apparent birefringence that depends on through-plane angle of fiber tracts. Our method measures fiber orientations in full spherical orientation space and birefringence free of the orientation bias to reveal the true structural anisotropy of myelinated fibers. By formulating a birefringence vector composed of true birefringence value and 3D optic axis orientation vector, we verified that our method could recover accurate estimates of a large range of through-plane angles at different in-plane orientations. We showed 3D fiber orientation maps in full sections of a human cerebral hemisphere and a brainstem. The 10 μm resolution maps revealed intricate fiber organization far beyond the capabilities of non-invasive imaging techniques such as dMRI. To our knowledge, this is the first thorough characterization of 3D fiber axis measurement in the human brain with microscopic spatial resolution, using PS-OCT based methods.

Our major contribution is the accurate and robust reconstruction of 3D axis orientation in full spherical space, using only two illumination beams. Previously, quantifying 3D orientation by PS-OCT via multiple beams has been proposed based on either apparent birefringence deviation or geometrical transformation of in-plane orientations ^20,21^. A common caveat of those methods is that they did not show the estimation for the full angular space, a critical factor for studying structural connectivity in the human brain as the fiber tracts run in all different directions. Early studies suffered from the limitation that only the axis lying within the 2D plane spanned by the illumination beams could be recovered ^21^. Later methods removed this limitation, but required a third illumination to complete the estimation in the 3D space, substantially increasing the acquisition time and hardware complexity ^20,22^. The birefringence vector defined in this work makes use of both apparent birefringence and in-plane orientation metrics in PS-OCT. This formulation obviates the requirement for the third illumination angle and allows the retrieval of the vector from only two sets of measurements. This implementation significantly reduces the engineering complexity and the acquisition time for large-scale samples. Another work utilized a dual-angle PS-OCT with a large tilt angle, based on geometric reasoning but failed in our hands for axes co-planar to the two illumination beams ^23^. Here we significantly advanced the 3D axis orientation reconstruction with a small tilt angle of 15 deg as the second illumination, which reduces the intensity loss of light caused by refraction and yields a high signal-to-noise-ratio (SNR). The small tilting angle also simplifies image registration between the illuminations. Importantly, we showed that our method was able to recover the axis orientation in full spherical space (Fig. 2–4) with high angular precision (Fig. 3). Therefore, the proposed PS-OCT method enables quantitative characterization of comprehensive fiber configurations in the human brain.

Signal-to-noise-ratio (SNR) plays an important role in 3D orientation estimation.Especially in ex vivo human tissue imaging, the long fixation time significantly increases scattering, and therefore reduces the SNR and the degree of polarization over depth. In our procedure, we applied refractive index matching to our samples to increase the SNR in the PS-OCT measurements. Previously, we showed that index matching significantly improved the SNR and thus enhanced the accuracy of in-plane orientation and apparent birefringence measurements ^27^. It is noted that the index matching procedure may not be necessary in in vivo tissue imaging, where scattering is much smaller. Previous studies using fresh or lightly fixed birefringent tissue have demonstrated success of 3D axis measurement without refractive index matching.

Imaging connectional anatomy at the microscopic level is imperative for understanding the human brain. Although 3D fiber orientation has been previously computed from images of histological stains or measured with polarized light imaging, these approaches involve registration of a stack of 2D brain sections to recover volumetric structure, a challenging task due to the nonlinear distortions induced by sectioning and tissue mounting ^28^. This renders continuous tracing of 3D fibers across sections at microscopic scale impractical ^4^. In this work, we can achieve this by incorporating a serial sectioning PS-OCT technique ^15^ into 3D axis imaging, whose unique advantage is that the blockface of the tissue is imaged before it is cut, therefore eliminating tissue distortion and preserving volumetric inter-slice alignment with microscopic accuracy ^28^. We showed examples of 3D axis orientation maps in volumetric reconstructions of brainstem (Fig. 6–8), revealing the extraordinarily complex organization of fiber bundles at microscopic scale, including twisting and inter-weaving with different angles. Importantly, the reconstructed volume preserves the integrity of fiber pathways in 3D without aforementioned distortions. It is noted that, in the current study, we did not perform tractography for the entire brainstem sample, as standard dMRI tractography tools are not suitable for the full complexity of fiber architectures at the microscale. Future development involves building dedicated tools for microscopic tractography. Further development of serial PS-OCT and computational tools opens the potential of large-scale neuroanatomical connectivity studies at unprecedented details.

Diffusion MRI has been widely used in clinical and preclinical studies. However, the millimeter resolution of dMRI is insufficient in regions where multiple microscopic fascicles converge, forming complicated structures like crossing, branching, or fanning ^29^; microscopic imaging needed to validate and to guide dMRI technologies in those challenging regions ^4,5^. PS-OCT has been used to evaluate the accuracy of dMRI orientation obtained with different q-space sampling schemes, spatial resolutions, and orientation reconstruction methods in human white matter samples ^30,31^. The ability to measure 3D axis orientation in large-scale samples will enable full validation of dMRI.

Birefringence has been used as a biomarker in myelin involved neurological disorders, such as multiple sclerosis ^32^, stroke ^33^, cerebral amyloid angiopathy ^34^, and Alzheimer’s disease ^35^. Previous studies have found myelin defects in those diseases and an associated birefringence reduction. The alteration of birefringence is an indication of changes in the integrity and alignment of myelin sheaths. One drawback of traditional PS-OCT is that it measures an apparent birefringence that depends on the fiber orientation and hence may yield biased estimates of myelin content (Fig. 5). In this work, we introduced true birefringence, which is free of the influence of through-plane fiber orientation and hence yields an unbiased measure of myelin content (Fig. 5). The results from the whole coronal section showed a significantly higher true birefringence in corpus callosum, indicating a high myelin content, which is supported in previous studies using different imaging modalities (e.g., see red arrows in Fig. 6). Our method to estimate the birefringence potentially advances our understanding of myelin distribution in the brain, and hence merits further investigation and validation. It is worthwhile mentioning that the birefringence estimation is also dependent on the objective used in the microscope that depicts different scales of myelin arrangement ^36^. With similar resolution and slice thickness, the birefringence value has been reported 2–4 × 10^−4^ in the corpus callosum of the rodent brains using either PS-OCT or polarization microscopy ^37,38^. The ability to quantify myelin content in tens to hundreds of cubic centimeters of human brain tissue promises to advance our understanding of the mechanism of various neurodegenerative diseases.

Our method was able to map fiber networks in the human brainstem, one of the most complex in the human brain, in which small fiber tracts form an intricate 3D network ^39^. It is worth emphasizing that our method can be applied to quantitatively image the 3D axis orientation and true birefringence in other fibrous tissues as well, such as cartilage ^21^, capsular ligament ^40^, enamel ^41^, uterus ^42^ and heart ^23^. Characterizing the 3D fiber organization in these tissues will play an essential role in understanding the normal functions and the pathological impact in diseases.

One limitation of the current study is that we applied our method on the en-face orientation data that integrated the information along the depth, assuming a constant axis of fibers bundles within a small volume intersected by the two beams. In future work, we will develop methods to extract local birefringence and orientation ^21,43–45^ to further obtain depth-resolved 3D axis orientation. We also plan to develop an automatic serial PS-OCT system with two beam illuminations for large-scale ex vivo human brain imaging. The combination of quantitative 3D axis orientation and serial sectioning PS-OCT opens new avenues towards imaging connectional anatomy in the healthy and diseased brain at unprecedented scales.

Methods

Human brain samples

We used samples excised from two neurologically normal brains from the Massachusetts General Hospital (MGH) Autopsy Suite for this study (age at death: 70 years and 46 years; 1 male and 1 female) and one neurologically normal brainstem sample from a subject enrolled in the Université de Tours body donation program (age at death: 76 years, male). Before death, the participants provided written proof of their willingness to become body donors and consented to the use of their bodies for educational or research purposes. The post-mortem intervals did not exceed 24 hours. The brains were fixed by immersion in 10% formalin and were cut into smaller blocks. We cut two blocks in the corpus callosum region from one of the MGH samples and one coronal slab from the other MGH sample. We applied refractive index matching to all the samples by incubations in 2,2-thiodiethanol solutions as described in our previous work ^27^. The ex vivo imaging procedures are approved by the Institutional Review Board of the MGH.

System, data acquisition and analysis

We used a custom-built automatic serial sectioning PS-OCT (as-PSOCT) system ^15,16^ to validate our 3D orientation method (Fig. 1a). The system consists of a spectral domain PS-OCT centered at 1300 nm, motorized xyz translational stages, and a vibratome for tissue sectioning. The PS-OCT is a free-space system. A quarter wave plate is placed in the sample arm to ensure circularly polarized light on the sample. Using a Jones matrix formalism, the retardance and the optic axis orientation are obtained by the amplitude and the phase of the complex signals, respectively, on the two polarization channels. Automatic imaging and sectioning of brain blocks is controlled by custom-built software for coordinating data acquisition, xyz stage translation, and vibratome sectioning. The samples were first flat faced with the vibratome and then imaged using a scan lens (OCT-LSM3, Thorlabs, Newton, NJ). The range of imaging depth was 2.6 mm with an axial resolution of 4.2 μm in tissue.

We imaged the sample with normal illumination and illumination tilted by 15 deg about the y-axis (Fig. 1a). To implement the tilted illumination, we either mounted the sample on a rotation stage or incorporated a rotation stage into the PS-OCT scan head. The physical effect of tilting the beam or the sample is the same ^20^. For the normal illumination, one volumetric acquisition covered a field of view (FOV) of 3.5 mm × 3.5 mm. The FOV was 3.5 mm × 2 mm for the 15 deg tilted illumination. Automatic tile-scan was used to cover the whole sample, with a 20% overlap between tiles in each direction of the raster scan. For volumetric imaging of large-scale samples, a 150 μm thick slice was removed from the tissue surface by the vibratome to expose the deeper region until the whole volume was imaged.

We obtained measurements of the apparent birefringence and the optic axis orientation from the two illumination angles. The apparent birefringence was calculated by performing a linear regression of the retardance profile along 150 μm in depth and extracting the slope of the linear fit. En-face optic axis orientation images were determined by the peak of a histogram formed by binning the measured orientation values into 5 deg intervals along 150 μm depth. The apparent birefringence and the en-face orientation data were then used to compute 3D axis orientation and true birefringence.

Image registration between the normal and titled illuminations

We applied a demon affine registration algorithm ^46,47^ to register the images between the normal and titled illumination incidences. We first registered the apparent birefringence maps from the two incidences and obtained a transformation field. Then we applied this transformation field to register the optic axis orientation maps.

Formulation of 3D orientations in normal and tilted spherical coordinate system

To estimate the 3D optic axis via two illuminations, we define a birefringence vector $[eqn]$ , where the magnitude $[eqn]$ is the true birefringence value and the unit vector $[eqn]$ represents the 3D orientation. The in-plane ( $[eqn]$ ) and through-plane ( $[eqn]$ ) orientations are obtained by vector calculus of $[eqn]$ in spherical coordinates. Let $[eqn]$ be the unit vector parallel to the illumination beam. We use two illumination beams, one along the z-axis, $[eqn]$ , and one rotated about the y-axis by $[eqn]$ , $[eqn]$ . We obtain two measurements of the apparent optic axis orientation $[eqn]$ and $[eqn]$ , and two measurements of the apparent birefringence $[eqn]$ and $[eqn]$ . This yields an apparent birefringence vector $[eqn]$ under the normal illumination $[eqn]$ (Fig. 1b, left panel) and $[eqn]$ under the tilted illumination $[eqn]$ (Fig. 1c, left panel). It is noted that the apparent birefringence values in each of the two frames are related to the true birefringence ^48^ as $[eqn]$ and $[eqn]$ , where $[eqn]$ and $[eqn]$ are the polar angles $[eqn]$ in the two illumination frames.

We define an estimated birefringence vector $[eqn]$ , which in turn defines its estimated in-plane orientations and polar angles, $[eqn]$ and $[eqn]$ , respectively, in the two frames. Therefore, the estimated apparent birefringence vector $[eqn]$ under the normal illumination $[eqn]$ can be written as,

$$[eqn]$$

where the expression in the parenthesis projects the birefringence vector into the plane orthogonal to $[eqn]$ , thereby scaling it by $[eqn]$ . The square root adds the required second $[eqn]$ . Similarly, the estimated apparent birefringence vector $[eqn]$ under the tilted illumination $[eqn]$ can be written as,

$$[eqn]$$

Strategy of xy-axis swapping for robust 3D axis estimation

The above computational model works well for the majority of the fiber orientations. However, a problem arises when applying Eqs. 1 – 3 to the fibers with in-plane orientations running along the y-axis (i.e., $[eqn]$ ). Since the optic axis measurement ranges between −90 deg and 90 deg, the distribution of the noise at $[eqn]$ forms a U-shape (bounded by −90 deg and 90 deg, see details in Fig. S1), which deviates from the close to normal distribution of the noise at intermediate angles in Eq. 1. Consequently, the computational estimation is sub-optimal. To overcome this problem, we interchanged the x-axis and y-axis for those fibers in post-processing, to effectively convert the optic axis orientation to close to 0 deg. As a result, the noise distribution is returned back to an approximately normal one. The tilted illumination in the swapped xy-axis is equivalent to tilting about x-axis by $[eqn]$ , which yields the apparent birefringence vector as $[eqn]$ . Correspondingly, the estimated apparent birefringence vector is written as,

$$[eqn]$$

where $[eqn]$ is the new direction vector of the tilted illumination beam about the x-axis. Therefore, we combine the two sets of equations and adapt a criterion that Eqs. 1, 2, and 3 are used to retrieve fibers with optic axis $[eqn]$ < 45 deg and Eqs. 1, 2, and 4 are used to estimate fibers with optic axis $[eqn]$ ≥ 45 deg.

Optimization of 3D orientation

The optimization algorithm takes the registered en-face orientation ( $[eqn]$ and $[eqn]$ ) and apparent birefringence ( $[eqn]$ and $[eqn]$ ) from normal and tilted incidences as inputs, converts them into the birefringence vector $[eqn]$ and $[eqn]$ , and outputs the estimated $[eqn]$ , which yields the true birefringence value $[eqn]$ and the 3D orientation vector $[eqn]$ (see details in previous section). The 3D orientation vector is further converted into in-plane and through-plane orientations for better visualization. We used the function fminsearch in Matlab (Mathworks, Natick, MA) to minimize the difference between the estimated and the measured birefringence vectors. The initialization values were set to be $[eqn]$ and the optimization stopped based on an error tolerance of 1e-20 in the fminsearch function.

Quantifying fiber orientation distribution

To quantitatively analyze the orientations obtained from the uniform corpus callosum sample, we constructed fiber orientation distributions (FODs) by generating histograms of the orientation map in every 0.5 mm × 0.5 mm region using a bin width of 5 deg. We found the primary peak and full width at half maximum (FWHM) of each FOD, then calculated the averaged peak direction and FWHM over the FODs in a 2.5 mm × 2.5 mm or 3 mm × 3 mm region of interest (ROI). Pixels with low birefringence values were excluded from the analysis.

Microscopic tractography

Tractography was generated for the volumetric brainstem data, using a DTI model in Diffusion Toolkit ^49^. The estimated 3D axis orientation was used to track the fiber using the streamline algorithm with an angular threshold of 60 degree and no threshold on the intensity. Due to the high resolution of serial sectioning PS-OCT, fiber tracking on the full 3D data would generate an extremely large data size beyond the ability of any visualization tools. As a result, we manually segmented four white matter structures in both brainstem hemispheres and performed tractography on each segmented structure of interest. The tracts were visualized in TrackVis (Version 0.6.1; trackvis.org).

Diffusion weighted MRI and tractography

The coronal slab sample was packed in fomblin and scanned in a small-bore 4.7 T Bruker BioSpin MRI scanner (maximum gradient strength 660 mT/m) using a diffusion-weighted 3D echo-planar imaging sequence with the following parameters: 0.25 mm isotropic resolution, TR = 500 ms, TE = 47.8 ms, 20 segments, 5 averages, $[eqn]$ , $[eqn]$ , b-value = 4,000 s/mm^2^, one b=0 volume and 12 volumes with non-colinear gradient directions. After denoising ^50^ and correction for eddy-current distortions ^51^, deterministic tensor tractography ^29^ was performed, seeding in every white-matter voxel.

The brainstem sample was packed in fombin and scanned in a 3T Siemens Trio MRI scanner using a diffusion-weighted SSFP sequence with the following parameters: 0.3 mm isotropic resolution, TR = 31.19 ms, TE = 25.6 ms, flip angle = 35°, bandwidth = 149 Hz/Px, b-value = 20,000 s/mm^2^, one b=0 volume and 30 gradient directions. Diffusion tensor imaging (DTI) model was used to extract the fiber orientation.