2.1.

Tissue Preparation

The Animal Research Ethics Committee of the Montréal Heart Institute approved all surgical procedures in accordance with the Canadian Council on Animal Care recommendations. Five C57Bl/6 mouse brains were used for this study; the brains are part of another ongoing investigation. The tissue preparation procedure followed the methodology presented by Ragan et al.⁵ Briefly, mice were anaesthetized under 2% isoflurane and perfused with 20 ml phosphate buffered saline and then by a mixture of 4% paraformaldehyde with 1% gadolinium. Each mouse brain was skull extracted and embedded in 4% agarose cylindrical blocks for serial imaging.

2.2.

dMRI Acquisition

One of the brains was imaged ex vivo with a standard three-dimensional (3-D) spin echo diffusion MRI sequence,²¹ using an Agilent 7 Tesla scanner equipped with $600\mathrm{mT}/\mathrm{m}$ gradients and a custom-built 1-loop cylindrical coil. Sequence parameters were: $\mathrm{TE}=0.021\mathrm{s}$, $\mathrm{TR}=0.4\mathrm{s}$, 30 gradient-encoding directions with $b=2079\mathrm{s}/{\mathrm{mm}}^2$ and five acquisitions with $b=0$, $\delta =5\mathrm{ms}$, $\mathrm{\Delta }=12\mathrm{ms}$, $\text{gradient amplitude}=320\mathrm{mT}/\mathrm{m}$, $\mathrm{FOV}=19.2\times 12.8\times 12.8\mathrm{mm}$, and an acquisition matrix of $128\times 96\times 96$ giving a resolution of $150\times 133\times 133\mu \mathrm{m}$, for a total acquisition time of 36 h. The diffusion MRI data preprocessing consisted of three steps: (1) correction of the field homogeneity artefacts,²² (2) reduction of the Rician noise bias,²³ and (3) volume resampling to 0.133 mm isotropic resolution.²⁴ Motion artefacts corrections were unnecessary due to the ex vivo nature of the sample. Then, in-house implementations^25,26 of diffusion tensor imaging and high angular resolution diffusion imaging reconstructions were performed using the Dipy library.²⁷ The FA was computed from the local diffusion tensors with a nonnegative least square method. The constrained spherical deconvolution of Dipy^25,28 was used to reconstruct the fiber orientation distribution functions (fODF). The principal directions of diffusion in each voxel and the AFD_max^17,18 were extracted from the fODF. AFD_max is the maximal value of the fODF on the sphere, and it can be interpreted as the maximum of the apparent fiber density (AFD). Finally, the number of fiber orientations (NuFO) within a voxel was computed with the method presented by Dell’Acqua et al.²⁶ and using a data-driven threshold set to 1.5 times the AFD_max values in the ventricles. The purpose of this threshold in the ventricles was to remove noisy peaks of the fODF, which are simply high-frequency peaks with low fODF amplitude. The ventricles are used because NuFO should be 0 there.

2.3.

SS-OCT and Image Acquisition

An in-house automated high-throughput histology setup combined with a swept-source optical coherence tomography (SS-OCT) microscope was used to image the agarose embedded rodent brains (Fig. 1). The SS-OCT system was reported in previous publications.^30,31 The swept-source laser was operated at a central wavelength of ${\lambda }_0=1310\mathrm{nm}$ with a tuning bandwidth of $\mathrm{\Delta }\lambda =100\mathrm{nm}$ (Axsun, 1310 Swept Source Engine). A cropped Blackman apodization function was used to reduce the side lobes caused by the rectangular shape of the swept-source spectrum. The $3\times$ microscope objective (Thorlabs, LSM04 Scan Lens) was enclosed in a 3-D printed watertight immersion chamber ending with a glass coverslip. This immersion chamber had two purposes: protect the scanning lens from the water and biological tissue debris created by the slicing process and impose a constant air/water column in the sample arm. The measured SS-OCT system resolutions in water were $r_{x,y}=8\mu \mathrm{m}$ laterally and $r_z=10\mu \mathrm{m}$ axially. Its sensitivity was 98.5 dB, and its sensitivity roll-off along depth was $-0.01\mathrm{dB}/\mu \mathrm{m}$.

Fig. 1

(a) Main components of the serial massive histology setup. (b) Parameters of the motor displacement model. (c) Volume rendering of a reconstructed brain (visualization made with the Blender software²⁹).

Imaging was achieved by sequentially cutting thin tissue slices (around $200\mu \mathrm{m}$) with a vibrating blade⁵ and by moving the sample under the microscope objective with a motorized stage (Zaber, T-LSR150B). At each motor position, the sampling beam was raster scanned over the objective’s field of view using galvanometric mirrors. An OCT A-line was acquired for each sampling point, thus resulting in a mosaic of volumetric OCT tiles for each tissue slice. After a slice acquisition, the sample was moved axially using a motorized jack (Thorlabs, L490MZ/M), and this process was repeated until the whole tissue was sliced and imaged. Thousands of OCT volumetric tiles were acquired, each covering a size of $\mathrm{\Delta }x\times \mathrm{\Delta }y\times \mathrm{\Delta }z=2.5\times 2.5\times 0.8{\mathrm{mm}}^3$ with an anisotropic sampling resolution of $4.9\times 4.9\times 6.5\mu {\mathrm{m}}^3$. The volumes were assembled using the postprocessing reconstruction method presented in Sec. 2.4.

The SS-OCT full depth range was 5 mm, but only 0.8-mm A-lines were kept for the data reconstruction. This range spanned depths between 195 to $975\mu \mathrm{m}$ from the zero-delay position of the OCT. This range usually comprises water and tissue. Indeed, due to variations in cutting depths, the water–tissue interface was typically located between 395 and $570\mu \mathrm{m}$ from the zero-delay position. The upper and lower limits of the 0.8-mm depth range were chosen (1) to ensure that there was an overlap between consecutive tissue slices, (2) to remove the auto-correlation artefacts and internal reflections due to the coverslip, (3) to remove the areas affected by complex ambiguity artefacts, and (4) to remove deeper sample areas where the signal was dominated by multiple scattered photons or by noise.

2.4.

Volume Reconstruction

The data reconstruction consisted of three steps: (1) tile registration, (2) $XY$ slices stitching, and (3) $Z$ slices stitching.

First, a displacement model of the motorized sample stage was used to estimate each tile position within the mosaic. The displacement parameters used were the tile overlaps $(O_x,O_y)$, the angle $\varphi$ between the $X$ and $Y$ motor axes, and the angle $\theta$ between the laser scanning reference frame and the motorized stage reference frame [Fig. 1(b)]. Given a tile at the grid position $(i,j)$ within the mosaic, its associated motor position $\overrightarrow{D}(i,j)$ was given by

Finally, the displacement model described by the extracted parameters $(\widehat{\theta },\widehat{\varphi },{\widehat{O}}_x,{\widehat{O}}_y)$ was used to compute the $XY$ tile positions for each tissue slice. The advantages of using this model instead of directly computing the tile positions via image registration are that it is more robust to noise and it gives accurate position at multiple resolutions even when there is missing information. Moreover, it takes advantage of the motors’ high displacement resolution ($<1\mu \mathrm{m}$).

The OCT volumes were then stitched within each tissue slice using weighted average linear blending over neighboring tiles overlap areas

Fig. 2

Laplace blending example. (a) A tile $I(x,y)$ (red) is added to the mosaic $M(x,y)$ (blue). The Laplace equation is solved over the overlap area (violet). (b) Laplace blending weights solved for this geometry. The contour lines traced inside the overlap area indicate the extent of the modified blending weights for different values of the blending width parameter $w$.

The parameter $w$ was used to control the blending seams visibility. When $w=0$, there is no transition between the mosaic and the new tile. The extent of the blending area increases as $w$ becomes larger. When $w=1$, the blending is done over the whole overlap area. In other words, $w$ controls the steepness of the blending. A value of $w=0.3$ was used for the reconstruction presented in this study. Image blending acts as a spatial compounding that reduces the speckle noise in the overlap areas only. This results in a grid-shaped artifact in the reconstructed slices as no denoising occurs inside the tiles. Using the modified blending weights ${\alpha }^{\prime }(x,y)$, the spatial extent of the grid-shaped artifact was reduced.

The last data reconstruction step was to stitch the tissue slices together to get a complete 3-D volume. The depth translations between consecutive tissue slices were computed as follows. First, the whole OCT reflectivity tissue slices were smoothed using a 2-D $XY$ Gaussian kernel of standard deviation $k=25\mu \mathrm{m}$. Then, the image structure boundaries within each mosaic were enhanced by computing the normalized 2-D $XY$ image gradient magnitude. Using the normalized gradient maps, the cross correlation between each depth position of the $n^{\prime }\mathrm{th}+1$ slice and the last depth position of the $n$’th slice were computed. The depth associated with the maximum image gradient cross correlation was selected as the optimal translation between the $n$’th and $n^{\prime }\mathrm{th}+1$ slices. This procedure was repeated for all adjacent slices.

The use of the image gradient for this registration process was motivated by the OCT depth-dependent signal attenuation. Indeed, the tissue contrast at the bottom of a slice is attenuated by all the structures encountered by the sampling beam. These scattering structures are removed by the slicing process, which increases the signal measured at corresponding positions in the next slice. When using the OCT reflectivity to compute the depth translation between slices, the algorithm tended to align the top tissue layers together as they exhibited similar contrasts solely due to the fact that they are not affected by signal attenuation. Instead, the normalized image gradient modulus delineates the brain and fiber tracts boundaries and is not affected by signal attenuation.

Another aspect that needs to be considered when stitching the tissue slices together is the roughness and denivelation of the water/tissue interface: the slicing process does not always result in a clear cut and the tissue deformations can introduce holes in the reconstructed data. This is taken into consideration by detecting the water/tissue interface depth for each A-line and creating a data mask to keep only the pixel under the interface in the assembled volume. Finally, the tissue slices were stitched together using the 3-D blending weights given by solving the Laplace equation ${\nabla }^2\alpha (x,y,z)=0$ and the water/tissue masks. The 3-D Laplace equation was solved with a 7-stencil finite difference approach with a Jacobi iteration scheme. The boundary conditions were defined in the same way as for the 2-D case explained above, except that the domain of the moving volume being added was constrained by the tissue mask.

2.5.

Optical Attenuation Estimation

Vignetting effects were compensated for in each tile prior to using the stitching algorithm. The tissue tiles containing agarose were first identified by applying the Li threshold method³⁴ to separate them into background and foreground volumes. The background mosaic tiles were averaged together, and their $Z$-axis average intensity projection (AIP) was computed. This agarose AIP was used to compensate for the 2-D vignetting effect introduced by both optic misalignment and beam scanning.³⁵

As described in van Leeuwen et al.,³⁶ the sampling beam of a single-mode fiber reduces to a Gaussian beam, and its confocal point-spread function (PSF) adds a depth-dependent intensity contrast to each A-line. To compensate for this effect, multiple A-lines in agarose were averaged [see Fig. 3(a)]. Then, a tissue model affected by the axial component of the Gaussian beam PSF $h(z)$ was fitted on the average agarose profile to find the beam parameters: the focal plane depth $z_f$ and the effective Rayleigh length $z_R$

Fig. 3

(a) Example of an axial PSF extraction from an average A-line within agarose. (b) Schematic of the A-line signal extrapolation method used to reduce attenuation coefficients biases. (c) AIPs over $250\mu \mathrm{m}$ of a horizontal slice from a reconstructed mouse brain. (Left) Original OCT reflectivity contrast, (middle) reflectivity contrast compensated for the depth-dependent signal attenuation, and (right) the computed attenuation coefficient map. The white arrow indicates the direction of the MOA. The tissue was cut along the coronal plane with the cerebellum facing the microscope.

A single PSF profile was synthesized using the extracted parameters and used to normalize each A-line in the dataset. This method assumes that the depth-dependent signal attenuation in agarose is negligible compared with the signal roll-off introduced by the microscopes axial PSF.

The signal attenuation within each A-line was then computed. The common approach to tissue attenuation estimation is to fit an exponential function on an average A-line and to report the extracted signal decay with depth obtained from this regression. The fitted equation is usually a single-scattered photon model and assumes that the attenuation coefficient is constant over the A-line scanning range. In addition to the loss of axial resolution, another drawback of this approach is the artefact that appears when the A-line traverses inhomogeneous tissues (e.g., as shown for deeply embedded fiber bundles using PS-OCT).¹² To address these limitations, depth-resolved local tissue attenuation coefficients were estimated using the method introduced by Vermeer et al.,³⁷ which assumes that detected photons have experienced a single scattering event. A second assumption of this model is that the backscattering photons measured are a fixed fraction of the attenuated photons. Following this model, the attenuation at a given depth $i$ is given by

These estimated attenuation coefficient maps were used to compute an attenuation correction bias field $B(x,y,z)$ profile for each A-line using Beer–Lambert’s law. The original OCT tissue slices were normalized by these attenuation correction fields during the slice stitching process, which allowed for reducing signal intensity discontinuities between consecutive tissue slices [Fig. 3(c)]

2.6.

OCT Brain Template Creation

An average OCT brain template was created to be able to compare the contrasts and segmentation results of each individual brain (see Fig. 4). As will be discussed in Sec. 2.8, this OCT template was also used by the multimodal registration of the OCT brains with the dMRI data. The template creation methodology is based on existing algorithms in the literature^39,40 and uses the ANTs.²⁰ The normalization process was done with a subset of the assembled brains ($n=4$). Furthermore, to get a symmetric template, the normalization also included the brain’s mirror version obtained by flipping them along the midsagittal plane. Each brain was first aligned to the Allen mouse brain reference template⁴¹ using rigid and affine transforms and a mutual information metric. The average volume obtained by this coarse alignment was used as an initial reference template for the normalization procedure. The normalization was done iteratively, using the average template of iteration $k$ as a reference for the $k+1$ iteration. Each iteration consisted of a sequence of rigid, affine, and symmetric normalization diffeomorphic (SyN)⁴² registrations done at four scales each and using the mutual information metric. A total of $M=20$ normalization iterations were done. Finally, the OCT mouse brain templates were registered to the Allen mouse brain template using a sequence of rigid, affine, and SyN registrations. For computational reasons, the template creation was done at a subsampled resolution of $50\mu \mathrm{m}\mathrm{per}\text{voxel}$. These OCT mouse brain templates will be publicly available on our laboratory website.⁴³

Fig. 4

(a) Horizontal slices of the OCT reflectivity brains used for the template creation. (b) Mirrored version of the brain shown in (a) that were also used by the template creation algorithm. (c) Resulting OCT reflectivity brain template after $\mathrm{M}=20$ iterations. All brains were downsampled to an isotropic resolution of $50\mu \mathrm{m}\mathrm{per}\text{voxel}$.

2.7.

OCT Brain Image Segmentation

Given the observed visual contrast between white and gray matter, the Atropos algorithm,⁴⁴ an $n$-tissue segmentation method developed for MRI data and part of ANTs, was used for segmenting the brain tissue from the combination of the OCT reflectivity and attenuation contrasts. The reflectivity and attenuation contrasts were combined to form a 4-D matrix, the first three dimensions being the $X$-, $Y$-, and $Z$-axes and the fourth dimension being a two-length vector with elements being the reflectivity and attenuation contrasts. A Markov random field with a radius of $1\times 1\times 1$ voxel and weight $w=0.3$ was used to regularize the segmented volumes. The mathematical expression of the multivariate $n$-tissue segmentation method is described elsewhere.⁴⁴ All brains were registered to the Allen mouse brain common coordinate framework (CCF)⁴¹ prior to the tissue classification procedure, which allowed visual comparisons of the segmented tissue classes with the Allen mouse brain atlas. This segmentation procedure was applied to all individual brains and to the OCT brain template.

The segmentation algorithm was initialized using $K$-means clustering combined with a Gaussian mixture model (GMM).⁴⁵ The optimal number of classes $n$ was determined by computing the Bayesian information criterion (BIC) for varying values of $n$. The GMM associated with the minimum BIC value was chosen.⁴⁶ Using this criterion, the optimal number of classes was found to be $n=5$. Thus, the tissue was segmented into five classes corresponding to (1) low and (2) high attenuation gray matter, white matter fibers (3) parallel and (4) orthogonal to the direction of the MOA, and (5) other white matter regions exhibiting high attenuation and reflectivity contrasts. These tissue class names were chosen based on the visual comparisons of the segmented labels with the Allen mouse brain atlas structures and based on the reflectivity and attenuation statistics per tissue label (see the results section and Fig. 7). To restrict the segmentation process to tissue voxels only, the brain and ventricles were masked prior to the Atropos algorithm application. The brain tissue was segmented with a global intensity threshold determined by the Otsu method combined with a morphological watershed algorithm. The ventricles were excluded from this segmentation process by delineating them using the ITK-SNAP⁴⁷ tool prior to the Atropos algorithm application.

Finally, the tissue labels obtained from the OCT brain template segmentation were used by the Atropos algorithm as a priori information to guide the segmentation in other assembled brains. An advantage of this approach is a reduction in the effect of acquisition variability and artefacts on the segmentation results by imposing an a priori bias to the classification problem. This was done with the fifth mouse brain that was also acquired with dMRI. The volume was first coregistered to the symmetric OCT template. Then, the Atropos algorithm was used as described above, but by replacing the initialization step by the OCT brain template’s tissue label map and using a prior weight of $w=0.15$. The average reflectivity and attenuation coefficients per tissue classes were then computed and reported in Fig. 7.

2.8.

OCT and dMRI Coregistration and Comparison

In qualitative estimates of the attenuation coefficient maps, white matter with different orientations was apparent in images suggesting a compensation for fiber orientation (e.g., see the anterior commissure in Fig. 9). To validate this, an assembled brain was coregistered to a FA map obtained from a dMRI acquisition of the same animal (Fig. 5). Registration templates were used to align the two imaging modalities together as this method was found to be more robust than directly aligning the individual brains.⁴⁰ To summarize the registration algorithm, the OCT reflectivity and attenuation maps were aligned to the OCT brain templates. Then, the precomputed transforms from the OCT brain template to the Allen mouse brain were applied. Similar registration steps were done with the dMRI metrics. The FA volume was first aligned to a publicly available dMRI mouse brain template.⁴⁸ This FA template was then coregistered to the Allen mouse brain. Thus, both the OCT brain volumes and the dMRI data were mapped to the Allen mouse CCF, which allowed voxelwise comparisons of the contrasts. All registrations were done with ANTs with a sequence of rigid, affine, and SyN transforms at four scales and using the mutual information similarity metric. The coregistration was done at a resampled isotropic resolution of $50\mu \mathrm{m}$ to match the resolution of the OCT brain template.

Fig. 5

Example of the multimodal registration algorithm. The OCT brain template generated in this work and a publicly available dMRI mouse brain template⁴⁸ were used as intermediary registration templates. Each arrow indicates the registration of a given volume onto another (for example the OCT brain was mapped to the OCT brain template and then mapped to the Allen mouse brain).

The angles between the direction of the MOA of the OCT brain and the fiber bundle orientations given by the water diffusion data were computed. The OCT signal dependence on fiber orientation could then be investigated. Other MRI metrics were also compared with the reflectivity and attenuation contrasts: the FA, the AFD_max, and the number of fibers orientations within a voxel (NuFO). Plots of the average OCT reflectivity and attenuation signal as a function of these dMRI metrics were generated as follows. First, each dMRI metric was separated into 32 equidistant intervals. Then, 3-D masks indicating the spatial location in the brain of the voxels included in each dMRI metric interval were created. These masks were applied to the OCT reflectivity and attenuation maps. The plots in Fig. 8 report the average reflectivity and attenuation values for each dMRI metric value interval. The error intervals are expressed as the standard error of the mean ($\mathrm{SEM}=\sigma /\sqrt{N}$), with $N$ being the number of voxels in each interval and $\sigma$ being the standard deviation of the reflectivity/attenuation values for these voxels.

3.1.

Volumes Acquisition and Reconstruction

Typically, an acquisition of a mouse brain using a $3\times$ microscope objective generated a mosaic of shape $7\times 5\times 70$ at a sampling resolution of $4.8\times 4.8\times 6.5\mu {\mathrm{m}}^3$. This represents around 320 GB of disk space. Figure 6 shows a reconstructed coronal brain tissue slice at the original microscope’s sampling resolution. Due to the large amount of data, the OCT volumes were resampled prior to reconstruction at an isotropic resolution of $25\mu \mathrm{m}\mathrm{per}\text{voxel}$, which resulted in a final brain volume size of around 300 MB. Figure 3(c) shows a horizontal section within an assembled brain at a resolution of $25\mu \mathrm{m}/\text{voxel}$ with three contrasts: (1) OCT reflectivity, (2) compensated OCT reflectivity, and (3) estimated signal attenuation. The compensated OCT map was obtained by dividing the OCT reflectivity signal by the cumulated attenuation for a given depth (Eq. 12). This normalization removed the signal attenuation caused by the sampling depth and by highly scattering structures traversed by the sampling beam. The remaining contrast in the compensated OCT map was due to local differences of refractive index in the brain microstructure. Indeed, the attenuation compensation increased the backscattering signal under highly scattering structures (e.g., under the corpus callosum or the arbor vitae). It also uniformized the backscattering contrast in homogeneous areas such as in the neocortex. For the white matter fibers, the backscattering contrast was higher due to the layered structure of the myelin sheaths. Also, Leahy et al.¹³ proposed that the backscattering coefficient changes based on the fiber tract orientation in relation to the direction of the incident beam. Our data corroborate this model, as illustrated by the anterior commissure olfactory and temporal limbs in Fig. 9.

Fig. 6

(a) AIP of an OCT reflectivity tissue slice reconstructed at the original anisotropic sampling resolution of $4.8\times 4.8\times 6.5\mu {\mathrm{m}}^3$. Neither the axial PSF nor the $XY$ vignetting effect were compensated for in this example to better represent the raw data aspect. Each volumetric tile has a field of view of $2.5\times 2.5\times 0.8{\mathrm{mm}}^3$ and an overlap with their neighbors of around 0.5 mm. The arrows indicate the mammillothalamic tracts (red), the columns of the fornix (blue), and the corpus callosum and the striatum (green). (b) Coronal, sagittal, and horizontal slices within the delineated yellow region in (a), showing fibers in the striatum and the corpus callosum. The tissue was sliced along the coronal plane, with the cerebellum facing the microscope objective.

The stitching algorithms were implemented in Python and were integrated as a suite of tools in a private galaxy server⁴⁹ running in a Docker container⁵⁰ on an $8\text{cores}/32\mathrm{GB}$ RAM computer (64 bits). This implementation strategy was chosen to allow multiple users to launch data reconstructions simultaneously. The galaxy server can be installed on a cluster engine to provide more processing power and thus reduce computing time if needed. The tissue preparation and acquisition process took $\sim 10\mathrm{h}$ (including perfusion, surgery, agarose embedding, and slicing). The reconstruction with the galaxy server took another 3 h per brain. Tests were done with the stitching algorithm on a computing cluster (called Briarée and located at Calcul Québec) using 4 nodes of 12 cores each and 24 GB of RAM per node. The reconstruction time using the computing cluster was reduced to about 1 h for brains at the full resolution. Thus, the acquisition process using this high-throughput histology setup took between $1/2$ and 1 day overall to image and reconstruct a whole mouse brain in 3-D at an isotropic resolution of $25\mu \mathrm{m}\mathrm{per}\text{voxel}$.

3.2.

Attenuation Coefficients Estimation

The 3-D attenuation map resulted in a higher contrast for all white matter fibers in the brain (see Fig. 9), even those parallel to the optical axis. For example, in the anterior commissure the olfactory and temporal limbs exhibit low and high OCT reflectivity contrasts, respectively, but a high attenuation value compared with gray matter. The average OCT reflectivity $⟨R⟩$ and attenuation coefficients $⟨\mu ⟩$ measured for each tissue class are reported in Fig. 7. To summarize, the largest attenuation coefficient was found in the fibers orthogonal to the optical axis ($24\pm 5{\mathrm{cm}}^{-1}$), followed by the myelin-rich tissue ($20\pm 3{\mathrm{cm}}^{-1}$), the fibers parallel to the optical axis ($20\pm 5{\mathrm{cm}}^{-1}$), the high attenuation gray matter ($15\pm 2{\mathrm{cm}}^{-1}$), and the low attenuation gray matter ($11\pm 1{\mathrm{cm}}^{-1}$). For the reflectivity contrast, the gray matter classes and the fibers parallel to the direction of the sampling beam had similar reflectivity (around 0.23). The myelin rich tissue ($0.35\pm 0.05$) and the fibers orthogonal to the direction of the incident beam ($0.55\pm 0.1$) exhibited larger reflectivity values. These results were computed using the OCT template segmentation labels. Similar results were obtained with the individual brain segmentation labels. Interestingly, the histograms in Fig. 7 show that it was the combination of reflectivity and attenuation that allowed the separation of each tissue class. Thus, if the segmentation was only based on the reflectivity contrast, it would have not been possible to separate the gray matter and fibers parallel to the MOA as they all have similar reflectivity values. Similarly, a segmentation based solely on the attenuation coefficient contrast would not be able to separate the three white matter classes as they all have similar attenuation coefficients.

Fig. 7

(a) Scatterplot of the OCT reflectivity and attenuation contrasts by tissue classes obtained from the OCT brain template segmentation results. For each label, $N=750$ sample voxels were selected randomly. (b) Average FA computed for each segmented tissue label. (c) OCT reflectivity and (d) attenuation coefficient histograms per tissue classes. The average values per label are expressed as $\text{mean}\pm \text{standard deviation}$. The label colors are the same as shown in Fig. 9.

The largest attenuation coefficient variations were within both fiber bundle classes. One cause of this variability is that the fibers in each group are not precisely orthogonal or parallel to the direction of the MOA, thus introducing intermediary attenuation values. In addition, the Markov random field regularization and the resolution downsampling can introduce partial volume effects that mix the contrasts from different tissue components. Other sources of attenuation coefficients variations include coregistration errors, artefacts introduced by the attenuation model, tissue slicing defects, and global brain orientation variations.

3.3.

Correlation Between dMRI and OCT

The correlations between the diffusion MRI metrics and the optical contrasts are shown in Fig. 8. These comparisons were made on a single assembled brain that was coregistered to its FA map. This map was measured for the same animal after its sacrifice and before skull extraction and histological acquisition. A threshold value of $\mathrm{FA}=0.4$ was used to separate the voxels into tissue and white matter tracts. This FA threshold was chosen to extract both the single-direction fiber tracts ($\mathrm{FA}>0.6$) and the crossing fiber areas, which are associated with FA values between 0.4 and 0.6.²⁶

Fig. 8

Correlation between the dMRI metrics and the OCT contrasts. The average reflectivity (a) and attenuation coefficient (d) are plotted as a function of FA for different fiber bundle orientations. Also, these OCT contrasts are reported as a function of the angle between the fibers and the direction of the MOA (b), the AFD_max (e), and the number of fibers orientations within a voxel (NuFO, f). The error intervals represent the standard error of the mean ($\sigma /\sqrt{N}$). (c) Main tissue orientations given by the dMRI acquisition. The white arrow represents the direction of the MOA given by coregistering the assembled brain with the FA volume. The average values in (e) and (f) were computed within the fiber bundles only by imposing a threshold of $\mathrm{FA}>.4$.

The optical contrast variation with fiber bundle orientation as measured by dMRI is shown in Fig. 8(b). Consistent with observations reported previously in the literature,¹³ the OCT signal was higher for fibers orthogonal to the direction of the MOA. For example, in Fig. 6(a), the mammillothalamic tracts (red arrows) and columns of the fornix (blue arrows) are parallel to the MOA and appear darker than the surrounding tissue, as opposed to the corpus callosum and striatum (green arrows), which contain fibers parallel to the MOA. This characteristic was also observed for the attenuation coefficients. This can be explained in part by the attenuation model assumption that states that the measured backscattered signal is a fixed fraction of the attenuated light. This assumption results in a linear relationship between reflectivity and attenuation, which can be observed in the Fig. 7 scatterplot.

The multimodal comparison of OCT with diffusion MRI also revealed that reflectivity and attenuation in fibers orthogonal to the direction of the incident beam increase with FA [Fig. 8(a)]. For fibers that are parallel to the incident beam, a negative correlation is observed between these measures. This result can be explained with partial volume effects and the orientation dependence of the OCT signal in fiber bundles. Each fiber bundle orientation population within a voxel will exhibit a different orientation-dependent OCT contrast [Fig. 8(b)]. For crossing fiber areas, low reflectivity (parallel) and high reflectivity (orthogonal) fibers are mixed and contribute to the measured OCT signal in different proportions. As FA becomes larger, the fibers within a bundle will be better aligned. If the bundle main orientation is orthogonal to the MOA, the reflectivity signal will increase as the proportion of low reflectivity parallel fibers in the bundles diminishes. This partial volume effect is consistent with the observed negative correlation of OCT contrast with the number of fibers orientations within a voxel [NuFO, Fig. 8(f)]. Some causes of FA variations in white matter are the presence of crossing fibers, fiber dispersion, fiber bundle density, or microstructural architecture of the cellular membranes.⁵¹ These fiber bundles characteristics should be considered when interpreting OCT contrasts in white matter.

Another source of OCT contrast variation seems to be the fiber bundle density, as shown by the negative correlation of both reflectivity and attenuation with the AFD_max [Fig. 8(e)]. The observed decrease of the OCT signal with density is counterintuitive and will require further investigation. A hypothesis to explain this effect is that the scattering anisotropy of fiber bundles increases with fiber density, in which case forward scattering becomes prominent, thus reducing the perceived reflectivity of the fiber bundle. The decrease of attenuation with density could be due to the use of an attenuation model that assumes a constant fraction of backscattered photons to attenuated photons. Future work focusing on fiber bundles exhibiting various density values is planned, including Monte Carlo simulations that can help better understand the origin of the OCT contrast.

3.4.

Tissue Segmentation

The tissue classification was first performed on all individual brains and on the OCT brain template without using prior information to guide the segmentation. The segmentation results for the brain template are reported in Fig. 9. Based on visual inspection, the five segmented tissue classes for each brain were similar. Quantitatively, an average Dice similarity coefficient of $0.82\pm 0.03$ was obtained when comparing the segmentation results among all the individual brains. This Dice score indicates that the segmentation method is accurate when similar classification initialization was used for each brain. The segmentation differences among individual samples can be caused by coregistration errors, morphological differences among the individual brains, or differences in cutting artifacts, for example.

Fig. 9

Five-tissue segmentation results using an OCT/attenuation mouse brain template ($n=4$). (a)–(c) A coronal slice and (d)–(f) a horizontal slice. (a) and (d) The OCT reflectivity contrast, (b) and (e) the estimated attenuation coefficients, and (c) and (f) the segmentation results for the same slices. The white arrow indicates the direction of the MOA. The red and green arrows in (d) point to the temporal and olfactory limbs of the anterior commissure, respectively. The green and yellow labels represent white matter fiber bundles parallel and orthogonal to the MOA, the dark and light-blue labels contain mostly gray matter tissues, and the red label represents other tissues with high myelin content (e.g., medulla, pons, and midbrain).

The optical contrast dependencies on fiber orientations were used by the segmentation algorithm to distinguish between fiber tracts and myelin-rich brain tissues. Indeed, the average FA values per tissue label [Fig. 7(b)] show that the fiber bundles orthogonal (yellow) and parallel (green) to the MOA exhibit a larger FA signal than those in the myelin-rich tissue (red). For example, the medulla (Fig. 9, red) contains many densely packed crossing fibers that result in high reflectivity and attenuation contrasts and in a low water diffusion FA, as opposed to the spinocerebellar tract (Fig. 9, green) that was classified as a parallel fiber tract due to its remaining orientation-dependent contrast. The low FA values computed in each tissue class can be caused by misclassification errors or by crosstalk between each label. This arises from the design of the segmentation method, which only considers the combined reflectivity and attenuation contrasts distributions without prior morphological information. As a consequence, different tissue types exhibiting similar contrasts can be classified into the same class, or voxels inside a single morphological structure are classified into different labels. An example of this misclassification can be seen in Fig. 9 where part of the corpus callosum was labeled as orthogonal fibers (yellow) and as myelin-rich tissue (red).

The segmentation results reveal that another source of OCT contrast seems to be related to the density of myelinated fibers in brain tissue, as defined by the ratio of neuronal cell bodies to neurites (myelinated axons and dendrites) volume fractions. Using optical coherence microscopy, Srinivasan et al.⁵² showed that the cell bodies exhibit lower scattering coefficients than myelinated fibers, which increase the signal attenuation in myelin rich tissue. This contrast mechanism could be responsible for the distinct classification of the thalamus and hypothalamus (Fig. 9, light blue), the medulla (Fig. 9, red), and the neocortex (Fig. 9, dark blue). Indeed, the scatterplot in Fig. 7(a) reveals that each of these tissue areas is associated with increasing reflectivity and attenuation values. Similarly, using polarization-sensitive OCT, Liu et al¹⁶ showed that the attenuation coefficient was significantly different among the MLs, the GLs, and the white matter regions of mice’s cerebellar slices. In histology, the granular and molecular layers are delineated based on their neuronal cell body densities, with the molecular layer being sparse and the granular layer being densely packed.⁵³ This relationship between the OCT contrast and neurite density will be evaluated in future work.