Scalable mapping of myelin and neuron density in the human brain with micrometer resolution

Optical coherence tomography (OCT) is an emerging 3D imaging technique that allows quantification of intrinsic optical properties such as scattering coefficient and back-scattering coefficient, and has proved useful in distinguishing delicate microstructures in the human brain. The origins of scattering in brain tissues are contributed by the myelin content, neuron size and density primarily; however, no quantitative relationships between them have been reported, which hampers the use of OCT in fundamental studies of architectonic areas in the human brain and the pathological evaluations of diseases. Here, we built a generalized linear model based on Mie scattering theory that quantitatively links tissue scattering to myelin content and neuron density in the human brain. We report a strong linear relationship between scattering coefficient and the myelin content that is retained across different regions of the brain. Neuronal cell body turns out to be a secondary contribution to the overall scattering. The optical property of OCT provides a label-free solution for quantifying volumetric myelin content and neuron cells in the human brain.

The cells, dendrites, and axons in the human brain are structured into cytoarchitectonic and myeloarchitectonic areas, based on cell type, size, density, and the density of myelin sheath surrounding the axons. Those structural components are the substrate for cognitive competencies and the specific locations of neuropathological processes [1][2][3] . Despite significant advances in imaging technology in the past decades, our understanding of human brain structures at 1-100 μm scale, in which neurons are organized into functional cohorts, is still limited. Quantitative features such as cell and myelin density have only been reported in a small number of subjects and over a small region of the brain [4][5][6][7][8][9] . Currently, Gallyas Silver stain and Nissl stain are two of the standard histology methods to study myelin content and Neurons in the human brain 10,11 . Despite their ubiquity, complex procedures have to be taken to apply these methods. One needs to cut the brain into tens of micron thickness slices and mount the slice on a glass slide, which induces inevitable tissue damage and distortions. The slices need to be stained and excessive pigment needs to be washed, which is labor intensive and subject to error and variability. Strict control of digitization such as illumination power and camera exposure time are pivotal for downstream quantitative analysis. After imaging individual slices, tremendous efforts are paid to reconstruct the volume, making it challenging for large scale study.
Optical coherence tomography (OCT) has been widely used in imaging the brain, which allows volumetric reconstruction of multiple cubic centimeters of tissue with mesoscopic to microscopic resolution. OCT utilizes the backscattered light to acquire tissue microstructural information 12 , which has proven to be useful in revealing cancerous tissue boundaries 13,14 , 3D vascular structures [15][16][17] , fiber tracts, and individual neurons and laminar structures across the cerebral cortex in rat and human brain 18,19 . OCT also allows quantification of tissue scattering. Fitting the OCT depth profile to a nonlinear model allows the calculation of tissue optical properties such as the attenuation coefficient and the back-scattering coefficient 20,21 . Myelin and neuron scattering have been described as the origins of tissue scattering. Wang et al. 20 found that the scattering coefficient is higher in white matter and subcortical nuclei regions with highly myelinated fibers, compared to less myelinated grey matter. Srinivasan et al. 19 found that myelinated fiber tracts are highly scattering while cell bodies have a lower scattering coefficient. Despite these investigations, quantitative correlations between tissue optical properties and these structural components have yet to be investigated.
Here, we report our work on quantifying the relationship between tissue optical properties and myelin content and neuron density in the human brain using automated serial sectioning OCT 22 . We established a computational model of the scattering coefficient with myelin content and cell density as the origins of scattering. By using the ground truth of Gallyas Silver stain and Nissl stain, we showed that the scattering coefficient has a strong linear relationship with the myelin content across different regions of the human brain. We also found that in grey matter, the cell body scattering serves as a secondary contribution to the overall tissue scattering and that the scattering coefficient has a moderate correlation with cell density. Our study provides a novel method for measuring myelin content and neuron density of the human brain tissues in a scalable sample size. The study also has important implications in evaluating brain diseases. As demyelination and neuron loss are two of the pathological hallmarks in neurodegenerative diseases such as Alzheimer's disease and Chronic Traumatic Encephalopathy (CTE) [23][24][25][26][27][28] , characterization of the optical property in diseased and normal brains will advance our understanding of pathological evolutions and their impact on complex functions.

Results
Optical property maps resembling quantitative histology. Figure 1 shows the histology and optical property maps for somatosensory cortex. Gallyas Silver stain (Fig. 1a) exhibits contrast among the supragranular layers which consist of pyramidal neurons, numerous stellate neurons and sparse axons (indicated by blue region), infragranular layers with large pyramidal neurons and axon bundles that connect to the subcortical structures (green region), and the white matter (red region), which mainly consists of highly myelinated axon bundles and glial cells. The Gallyas Optical Density (OD) map ( Fig. 1b) demonstrates that the supragranular layers have the lowest OD value, followed by the infragranular layers which have intermediate amount of myelinated axon bundles. The white matter exhibits the highest OD value due to the highly myelinated and densely packed axonal bundles. In addition, smaller features can be seen in the Gallyas OD map as well, such as the thin band of denser myelin content at the upper right region of layer IV (red arrow), possibly due to the high-density fibers in the outer band of Baillarger 29 . The Nissl stain and Cellular Occupation per Area (COPA) map (Fig. 1c,d) show contrast for cell bodies. The external granular layer (layer II) and the external pyramidal layer (layer III) exhibit the highest COPA value (red region), especially in the lower left part of the sample (black arrow), probably due to the higher neuron density. The internal granular layer (layer IV), internal pyramidal layer (layer V), and the fusiform layer (layer VI) present alternating contrasts (small green arrows). The white matter generally exhibits a low value of COPA, which is mainly attributed to the glia cells.
The µ s map (Fig. 1e) strongly resembles the Gallyas OD map. The white matter shows highest µ s because of the highly scattering myelin sheath. As the sparse axon branches into the cortex, µ s decreases accordingly. The supragranular layer shows the lowest µ s , due to the lack of myelin content. In addition, the thin band feature at the upper right region (red arrow) found in Gallyas OD map can also be seen in the µ s map. Apart from that, the infragranular layers (IV, V, VI) show additional laminar structures (small green arrows) similar to the COPA map but not in Gallyas OD map. Overall, the µ s map seems to be strongly correlated with myelin content and slightly modulated by the neuron scattering. The µ ′ b map (Fig. 1f) offers another feature dimension to aid in discriminating tissue types. It is noticeable that µ ′ b varies within the white matter, possibly highlighting fibers oriented within the image plane (green arrow). This is possibly because the fibers oriented within the imaging plane direct more back-scattered photons to the detectors than the fibers oriented through the imaging plane. Consequently, the µ ′ b map offers potential information about fiber orientation. The ratio of µ ′ b /µ s map (Fig. 1g) www.nature.com/scientificreports/ provides another useful feature to distinguish structures in the brain. Similar to the µ ′ b map, the ratio of µ ′ b /µ s map also highlights the region with fibers oriented within the imaging plane (green arrow). In addition, the ratio of µ ′ b /µ s map also highlights region of higher value in the superficial layers indicated by the blue arrow, the cause of such contrast requires further investigation. The AIP image (Fig. 1h) is a nonlinear function of the µ s and µ b maps, which provides an overall view of the tissue structure for reference.
The relationship of scattering coefficient and myelin content. The similarity between the µ s and Gallyas OD maps indicates that myelin is a crucial factor contributing to the brain tissue scattering. To quantitatively inspect the relationship, we plot µ s versus Gallyas OD in selective ROIs, which covered all the laminar layers as well as the white matter for the five brain regions. The scatter plots indicate a strong linear relationship between scattering coefficient and Gallyas OD, which is consistent with the Mie theory 30 . Therefore, we fit the data with a linear model and presented the results in Fig. 2a Remarkably, the five samples share similar linear relationship as indicated by the slope parameter ( k 1 ) of the fitting results, although individual samples have distinct patterns of µ s and Gallyas OD distributions. In the cerebellum (Fig. 2a), data points from the white matter, granular layer and molecular layer form three discrete clusters. However, they all follow a shared linear function ( k 1 = 5.79) in which higher Gallyas OD is associated with higher µ s . Similar patterns are observed in the somatosensory cortex (Fig. 2c), where the six cortical layers and the white matter group into three clusters and share a co-linear relationship ( k 1 = 7.48). Supragranular layers (layer I, II and III) form a cluster with the lowest Gallyas OD and µ s , infragranular layers (layer IV, V and VI) form a cluster with intermediate values, and the white matter cluster exhibits the highest values. The other two cortical regions of SupFrontal and BA21 (Fig. 2d-e) only present two discrete clusters. The supragranular layers and infragranular layers in both regions fall into a single cluster with low Gallyas OD and µ s . The white matter tracts show high values close to those of somatosensory cortex. Interestingly, the SupFrontal displays a within-cluster trend in both grey matter and white matter, suggesting a myelin gradient across cortical layers. Regardless, the slope parameters ( k 1 = 8.57 for SupFrontal and k 1 = 6.24 for BA21) demonstrate a similar relationship as those revealed in the cerebellum and the somatosensory cortex. In the hippocampus (Fig. 2b), due to complex anatomical structures, data points from different layers form a continuous distribution. For example, the fimbria, white matter and fornix show gradually increasing Gallyas OD and µ s values, while having large overlaps with the CA4 and dentate gyrus. Despite the different distribution pattern from other tissues, the fitting result is comparable with a slope parameter k 1 = 8.26.
An extraordinary linear relationship between Gallyas OD and µ s is revealed in all the samples (P < 10 -10 ) with high correlation coefficients (PCC > 0.85 for all brain regions, see Fig. 3b). The slope between Gallyas OD and µ s falls in a narrow range of 5.8 to 8.6. Slope variations are possibly due to different staining backgrounds, as well  Fig. 2f we combined all the data from the 5 samples and fit a single linear function, which reveals an average slope of 7.2 (correlation coefficient = 0.936, P < 10 -60 ). These results provide evidence that the linear relationship between optical scattering and myelin content holds in different regions of the brain. Therefore, scattering coefficient may potentially be used as an objective measurement of the myelin content in the human brain.
The relationship of scattering coefficient and joint myelin content and neuron density. The results of the linear model in session 2.2 suggest that myelin is a strong predictor of scattering coefficient in both grey and white matter of the human brain. In the cortex, neuronal cell body is another factor that should be taken into consideration when interpreting the scattering coefficient of the brain tissue. In this section, we first re-evaluate the relationship of µ s and Gallyas OD by including the factor of COPA. Next, we quantify the correlation between µ s and COPA to examine the cell body contribution in scattering. Figure 3a compares the slopes for Gallyas OD ( k 1 ) in univariate regression where only Gallyas OD was considered, and multivariate regression where COPA was included as well. We found that k 1 values were highly comparable between these two models with differences less than 10% in all brain regions Besides, the PCC of Gallyas OD in the multivariate regression (Fig. 2d) was found higher than 0.9 for four out of five samples (P < 10 -10 for all), with a slightly lower PCC of 0.6 in the hippocampus. Compared to the univariate regression, the partial correlation coefficient in the multivariate regression has only minor difference in most samples. The consistency of k 1 in the two analyses reinforces that the linear relationship between µ s and Gallyas OD largely preserves even when cell body scattering is taken into account. The high PCC in multivariate model consolidate the finding in session 2.2 that µ s is strongly correlated with myelin content across brain regions. In addition, when evaluating the R 2 of Pearson's correlation and NRMSE (Fig. 3c,d) for the regression, we found that adding COPA into the model only results in small improvements, which indicates that myelin content is a dominant factor to scattering coefficient. Figure 4 elaborates the correlation between µ s and COPA in the multivariate regression model. The slope for COPA ( k 2 ) exhibits large variations across the 5 different samples (Fig. 4a). In somatosensory cortex, Supfrontal and BA21, a positive k 2 was found, while in cerebellum and hippocampus, k 2 was negative. Two-sided t-test reveals a significant positive correlation only in the somatosensory region (P < 0.01, PCC = 0.45). The correlations   . 4b), but further statistical test fails to find a significance (P = 0.11 for SupFrontal and P = 0.35 for BA21). These results suggested that the neuronal scattering is a small contribution to the overall scattering coefficient and the effect varies across the brain. The negative correlation in the cerebellum and hippocampus was counterintuitive. However, it should be noted that the size of neurons in densely packed layers such as the granular cells in the cerebellum and the dentate gyrus in the hippocampus is much smaller than that of the other layers, which leads to a reduced scattering coefficient. The intercept in the multivariate regression exhibits large variation as well (Fig. 4c). As the intercept in the regression encompasses the unmodeled components to the tissue scattering, such as the extracellular matrix, the small intercepts in hippocampus, somatosensory, and BA21 indicate a negligible contribution from these remaining components. However, in cerebellum and SupFrontal cortex, we found a significant intercept (P < 0.001), indicating substantial scattering components remained. Overall, the fitting results for COPA in the model are coherent with findings in session 2.2 and 2.3 that µ s is dominant by the myelin factor. Neuronal cell body, however, only plays a secondary contribution to tissue scattering in the human brain.
Scattering coefficient differentiating neocortex from allocortex. As shown in Fig. 2, µ s and Gallyas OD exhibit distinctive cross-layer patterns among the 5 samples. As both metrics are strong predictors of myelin content, we examined the mean µ s and Gallyas OD in the white matter (red ROIs and red dots on scatter plots of Fig. 2) for all the samples. We found that µ s distribution exhibits a similar pattern as in Gallyas OD (Fig. 5a,b). Hippocampus and cerebellum have significantly lower Gallyas OD and µ s compared to the other three samples in the cortex, indicating a lower myelin content in the two regions of the brain. Anatomically, cerebellum and hippocampus belong to the allocortex while somatosensory, SupFrontal and BA21 belong to the neocortex. The two types of cortices possess different developmental trajectories of myelination (Miller et al. 36 ). The mean COPA in the grey matter (blue dots on scatter plots of Fig. 2) of the allocortex is significantly higher than that of the neocortex (Fig. 5c). Our results suggest that in addition to histology, the optical property obtained by OCT serves as a viable tool to differentiate the neocortex from the allocortex, with a distinction resulting from underlying myelin content.
The relationship of back scattering and joint myelin content and neuron density. We analyzed the relationships of µ ′ b and µ ′ b /µ s with respect to the COPA and Gallyas OD, respectively. For cerebellum, somatosensory cortex, BA21 and SupFrontal, there is a strong linear relationship between µ ′ b and Gallyas OD (Supplemental Figs. 5, 8), which is likely a result of the strong dependency between µ s and µ ′ b . However, there is no clear relationship between the ratio map of µ ′ b /µ s and Gallyas OD or COPA (Supplemental Figs. 6, 7).

Discussion and conclusions
Previous studies have demonstrated the ability of OCT to differentiate cortical laminar structures and to identify fiber tracts and subcortical nuclei 20,31 . Here, by using serial sectioning OCT, we have quantitatively investigated the contribution of structural components to the optical property in human brain samples and established a model of the scattering coefficient with regard to myelin content and neuron density. We have found that the scattering coefficient is strongly correlated with the myelin content (P < 10 -10 , PCC > 0.85) and that linear relationship is retained across different regions of the brain. The domination of myelin content in tissue scattering is reasonable considering the high index of refraction of myelin (n = 1.47) with respect to the surrounding aqueous environment (n = 1.35). The results from our study suggest that the optical property can be used a robust predictor for myelin content of the human brain. This strong correlation between scattering coefficient and myelin content has important implications in neurodegenerative diseases. It has been shown that the breakdown of myelin sheath is an indication of pathological abnormality and can result from several neurodegenerative diseases such as multiple sclerosis 32-34 , Alzheimer's disease 25,35 and chronic traumatic encephalopathy 23,24 . The quantitative measurement of myelin content could potentially be useful in characterizing the degree of demyelination in pathological brain samples. In addition, we have shown that scattering coefficient enabled the differentiation of various brain regions, such as neocortex and allocortex that have distinct myelination trajectories in development 36 . Previous studies have revealed that the degree of myelination and the order of maturation in the brain is associated with the vulnerability to psychiatric disorders [37][38][39] . Systematic characterization of myelination in the brain may provide a new avenue to map out the regional vulnerability to a range of brain disorders. Neuronal cell body turns out to be a secondary contribution to the overall scattering, and the correlation varies across different brain regions. In somatosensory cortex we found a significantly positive correlation (P < 0.01), indicating a strong laminar structure with differed neuron density and size, while in other brain regions we observed negative or moderately positive correlations. The lack of major contribution made by cell bodies to scattering in brain tissues have been reported in previous studies. Kalashnikov et al. 40 found light scattering from the neuron body contributes less than 10% of the observed backscattering signal when using cultured Hela neurons. Besides, Magnain et al. 41 found out that myelin density and fiber orientation could disrupt the identification of neuron cell bodies by using an optical coherence microscopy (OCM). As a result, the weak correlation revealed by the current linear model is not unexpected. Indeed, most of the OCT studies on brain cancer samples have reported difficulties in differentiating cancer of various stages from normal grey matter merely by scattering coefficient 42,43 . Our studies might provide an explanation for those challenges because cell body contribution is only a minor factor for light scattering in the brain. Further improvement of our scattering model may increase the sensitivity to neuron scattering. In our model, we assumed the same k1 and intercept between the grey matter and white matter ( Fig. 2 and Supplementary Fig. 9). However, as the refractive index of the extracellular space in grey and white matter could be different, allowing parameter tuning might result in a better fitting. When formulating the relationship between scattering coefficient and COPA, we assumed the phase function Qs to be constant for all neuron bodies (Eq. 4, "Method"), which bears a drawback if there's a large variation in the neuron size. Considering the dependency between phase function and scatterer size, a nonlinear model might improve the performance for correlating scattering coefficient with neuronal cell bodies.
Compared to histological methods, serial sectioning OCT offers a new window for quantifying myeloarchitecture and neuroarchitecture in human brain. Histological stains have been standard methods for studying myelin content and cell density. However, the outcome of histology heavily depends on the concentration of the contrast agent, the pH and temperature, and bleaching procedure for undesired pigment, which may vary across studies. Histology also requires a series of manual processes that are prone to human errors. Consequently, variations from slice to slice are inevitable 44 . In addition, sample must be cut into thin slices before being stained and imaged, which introduces tissue damage and distortions that are challenging to correct during 3D reconstruction. Serial sectioning OCT, on the contrary, uses intrinsic optical properties of tissue that does not depend on external contrast agents. The scattering coefficient is insensitive to system setup, incident power, and acquisition parameters. In volumetric imaging, the images are acquired on block-face prior to slicing, avoiding the vast majority of distortions incurred by tissue cutting and mounting. As a result, serial sectioning OCT generates images with consistent qualities across slices and samples. The technique, being automated and robust, is favorable to expand to larger sample size. Thus, serial sectioning OCT provides an attractive solution for quantifying volumetric myelin content and neuron cells in the human brain.
A few future directions in this field can be pointed out. First, another optical property named birefringence may directly relate to myelin content in the brain 45 . With polarization-sensitive OCT, we can measure the birefringence of myelinated fibers in addition to the scattering coefficient 46 . Inclusion of both parameters will enable a model for predicting and synthesizing myelin content in the human brain. Second, the local index of refraction may serve as another optical property to quantify neuronal characteristics 19 . Besides, high-resolution OCM has proven to be able to visualize neurons in brain tissue. Magnain et al. 41 used OCM to identify neurons that were validated by co-registered Nissl stain images in human entorhinal cortex. In addition to OCT, other imaging techniques such as two photon microscopy (2PM) could be useful to quantify neuron density as well. With elongated depth of profile, 2PM is able to cover a volume of tissue with high volume rate and generates an AIP of autofluorescence signals [47][48][49][50] . Hence, one of the future directions is to use multimodal techniques to obtain accurate measurement of cell and myelin content in scalable human brain samples. Lastly, the scattering coefficient measurement and correlation with myelin content could potentially provide valuable information in in vivo applications such as imaging guided neurosurgery. In a few related studies, Ben Arous et al. 51 have applied deep optical coherence microscopy in post-surgery mouse brain imaging for fiber tracking. Almasian et al. 52 have used the attenuation coefficient extracted from OCT to differentiate normal and glioma tissue during human brain surgery. OCT imaging of fresh brain tissues has revealed the intricate architecture of the underlying  18,53 . Future investigation of quantitative OCT with fast fitting algorithm will allow us to measure scattering coefficients of the brain in real time that has great potential in neurosurgical and neuromodulation guidance. Several empirical considerations ought to be clarified in this study. First, in the effort of correlating to histology, we thoroughly examined the quality of histological images and used the slices that have consistent staining intensity as ground truth. Yet there might be minor variations of staining that were not normalized among different samples, which could be one of the reasons for the variations of the slope and intercept values observed in the fitting results. Second, we formulated COPA based on the assumption that neuron bodies do not overlap on Nissl images, which may result in an underestimation in regions with high neuron populations. As an alternative, we also calculated the OD of Nissl stain and correlated it with the optical properties. The results were not significantly different from using COPA. Third, to obtain COPA we used an empirical thresholding method to segment the neuron body while excluding the smaller glia cells, the accuracy of which may depend on brain regions and staining quality. In the future, a deep learning based classifier may improve the segmentation 9 . Lastly, we acknowledge that there are more sources of scatterers in addition to myelin content and neuron that contribute to the tissue scattering. In this study, we used a simplified two-regressor linear model to include only myelin content and neuron and left the other factors as residuals of the model. Although our results showed strong evidence that myelin content is a major contributor, the intercepts that indicated the secondary source contribution varied with brain regions. The regions with larger secondary contributions such as the cerebellum or SupFrontal cortex will benefit from a model that accounted for contributions from other scattering components. A more complete model incorporating multiple variables into the scattering model, including contributions from protein, lipid bubbles, fiber structures, and extracellular space is a promising future direction to further improve accuracy. This type of multivariate Mie scattering model can be built by integrating refractive index, size and concentration of the various components into the existing model. For example, the water content of the fibers is associated with fiber density. Scattering changes when light travels from myelin sheath to water or vice versa where the index of refraction changes. Alternation of water content due to pathology changes the local density of axonal pack and may reflect on scattering coefficient. In addition, the geometry of fibers including the ratio of myelin sheath wrapping and thickness with respect to the total axonal diameter may affect the scattering coefficient. Further investigation of a Mie scattering model incorporating the layout of cylindrical scatterers can help to better understand this mechanism 54 . The orientation of fibers may affect the backscattering coefficient, which could be simulated by constructing the scattering model with off-axis fiber bundles. The release or accumulation of proteins, lipids, and other molecules may contribute to the overall scattering coefficient. Precise quantification of molecular concentration, index of refraction, and size is needed to build a molecule-specific Mie scattering model, which may be obtained by histological staining and super-resolution imaging 44,55,56 .
In conclusion, we have demonstrated the use of optical scattering obtained by serial sectioning OCT to measure myelin content and neuron density of human brain tissues. The scattering coefficient has a strong linear relationship with the myelin content across different brain regions, which promotes a robust label-free measurement for myelin with substantially reduced labor efforts comparing to traditional histology. The scattering coefficient was also moderately modulated by the neuronal cell bodies, the precise measurement of which requires complementary quantification or imaging techniques. Our approach has great potentials for enabling large-scale investigation of myeloarchitecture in various brain regions as well as studies of neurodegenerative processes in pathological brain sample.

Method
Sample. Two human brains (mean age 53.5 ± 12.0 year old, 1 male and 1 female) were obtained from the Massachusetts General Hospital Autopsy Suite. The brains were neurologically normal without a previous diagnosis of neurological deficits. The tissues were fixed by immersion in 10% formalin for at least two months 57 . The post-mortem interval did not exceed 24 h. The samples in the study consisted of five anatomical regions from the two human brains, including the cerebellum, hippocampus, somatosensory cortex, superior frontal cortex, and middle temporal area 21. Each sample was embedded in 4% melted oxidized agarose and cross-linked with a borohydride-borate solution for block-face imaging 58 . Serial sectioning OCT system. The serial sectioning OCT system was previously described 22 . Briefly, the system consists of a spectral-domain OCT system to measure the optical properties of the sample as well as a customized vibratome for tissue slicing. The light source was a broadband super-luminescent diode with a center wavelength of 1300 nm and full width half maximum bandwidth of 170 nm, yielding an axial resolution of 3.5 μm in tissue. The spectrometer consisted of a 1024-pixel InGaAs line scan camera operating at an A-line rate of 47 kHz. The total imaging depth was estimated to be 1.5 mm. The sample arm used a 10 × water immersion objective (Zeiss, N-Achroplan), yielding a lateral resolution of 3.5 μm. The volumetric imaging covered a field of view (FOV) of 1.5 × 1.5 × 1.5 mm 3 . The voxel size, which was defined by the stepping size of galvo mirror scanning laterally and obtained by the total imaging depth divided by the number of pixels axially, was 2.9 μm isotropic. The sensitivity of the system was 105 dB. Consecutive image tiles were obtained to cover the entire area of the sample and a 50% overlap was used between adjacent tiles. A customized vibratome was mounted adjacent to the OCT to cut off a 50 μm thick slice of the tissue upon completion of the full area scan. The slices were retrieved for histological staining. It is noted that the comparison of OCT and histology images in this study was not conducted on the same slices but slices nearby, as OCT slices have been preserved for other uses.
Quantification of optical properties. Optical property maps of the five samples are estimated from the data, by using a nonlinear model to fit the OCT depth profile as described in 20 www.nature.com/scientificreports/ Law to model light propagation in tissue, which is a negative exponential function parameterized by scattering coefficient µ s and back-scattering coefficient µ b (Eq. 1). The OCT depth profile is also confined by a confocal function h(z) imposed by the objective lens and a system sensitivity function H(z) 20 imposed by the spectrometer as below: Note that µ ′ b is the relative back-scattering coefficient, which includes OCT system parameters such as power and spectrometer configuration and efficiency. Light absorption is ignored due to its small contribution in the near-infrared spectral range. We used a simplified Gaussian PSF function, which depends on the focus depth Z f and the effective Rayleigh range Z Rs of the objective: We fit the measured OCT profile using a nonlinear optimization tool with the least-square error in Matlab 60 to obtain the optimal solutions for parameters of µ s and µ b . In general, the fitting process tries to find the optimal value for four parameters:µ s , µ ′ b , Z f and Z Rs for each A-line. Without any fitting constraints there is a strong inter-dependency between these parameters 20 . To reduce this inter-dependency and to make the fitting more stable and robust, we spatially parameterized the Z f and Z Rs following the procedure described in Yang et al. 2020. The Z f and Z Rs values were pre-calibrated using an Intralipid phantom with a comparable scattering coefficient as the tissue sample. As a result, we have reduced the number of unknowns in our fitting model to two: µ s and µ ′ b . To reduce the noise and errors in estimating these optical properties, the OCT A-lines were averaged over a 30 × 30 μm 2 area before fitting. A ratio map of µ ′ b /µ s was obtained afterward. In reconstruction, we used ImageJ to stitch the AIP of each OCT image tile to generate the XY coordinates for overlapping tiles. The tiles were then blended using a customized Matlab code to remove artifacts caused by intensity non-uniformity across the field of view. The down-sampled (30 × 30 μm 2 , same as the volumetric averaging in fitting) stitching coordinates were used to stitch the µ s and µ ′ b results from previous fitting steps.
Quantitative histological analysis. Selected slices from the serial block-face scanning were processed with Nissl stain 61,62 for neuron body identification and Gallyas stain 10,36 for characterizing myelin content. The stained slices were digitized by an 80i Nikon Microscope (Micro-video Instruments, Avon, Massachusetts) with a 4 × objective. The pixel size was 1.9 μm. The Gallyas stain images were normalized in each channel and converted to mean Optical Density (OD) 44 to directly represent the myelin content: where I is the RGB vector of the histology image with I i representing the intensity of the red, green, or blue channel, respectively. Each RGB channel was normalized separately before converting to OD in log scale. Then the average of three channels was used as the mean OD that represents the myelin content. Based on Mie scattering theory, assuming the same myelin architecture generating a constant phase function in the brain, scattering coefficient is expected to be proportional to myelin density which is represented by Gallyas OD.
For the cellular scattering, according to Mie theory 63 , the scattering of sphere particles is related to the phase function Q s , the number density of the sphere N s , and the cross-section area of the sphere A s : Assuming that the phase function Q s is constant for all neuron bodies, the scattering coefficient becomes proportional to the product of cellular density N s and cellular cross-sectional area A s , which we define as the cellular occupation per area (COPA). We constructed the COPA map in Nissl stain images by calculating the total area occupied by cell bodies and divided it by the total area of the tissue in a small neighborhood. The computation was conducted in two steps: (1) we segmented the cell bodies in the Nissl stain image using a threshold-based method, which converts the pixels values inside cell body to be 1 and the others in extracellular space to be 0; (2) by calculating the ratio of number of pixels with value equal to 1 to the total number of pixels within a small neighborhood (a 200 × 200 µm box), we obtained the percentage of that local area occupied by cell bodies. The box moved across the entire image to form the COPA map.
Generalized linear model. Taking both myelin and neuronal cell bodies into account, we examined the quantitative relationship between the optical properties and the Gallyas OD and COPA to reveal the sources of tissue scattering. OCT block-face and the histology images were acquired on slices from the same brain region that were a few millimeters apart. We manually drew 60 to 90 Region of Interests (ROIs) on corresponding slices from both modalities for linear regression analysis. The area of ROIs ranges from 200 to 500 pixels. In larger areas with more uniform intensity, such as the grey matter in hippocampus, the area of ROIs were expanded to a few thousand pixels to get more precise measurement. The total number of ROIs depends on the number of distinct cortical layers or subdivisions. On average we drew 20 evenly distributed ROIs for each layer in the µ s , Gallyas OD and COPA maps, respectively. www.nature.com/scientificreports/ We built a multivariate, generalized linear model (GLM) for µ s to include Gallyas OD and COPA as predictors. In the grey matter, we included both Gallyas OD and COPA as contributing factors, whereas in white matter, we only considered Gallyas OD since contribution from neuronal cell bodies was neglectable. The GLM relationship can be mathematically described as, where the matrix X contains the measured Gallyas OD and COPA values and β contains the linear coefficients that will be estimated: where n is total number of ROIs, k 1 is the slope corresponding to the myelin contribution to the tissue scattering, k 2 is the slope corresponding to the neuronal contribution to the tissue scattering, and b is the contribution from other components, such as scattering from the extracellular matrix. The coefficients were calculated by a pseudoinverse solution as, We used R 2 of Pearson's correlation and Normalized Root Mean Square Error (NRMSE) for examining the goodness of the fitting. Partial correlation coefficient (PCC), two sided t-tests and multiple comparisons were conducted to reveal the significance of the contributions from Gallyas OD and COPA, respectively.
It should be noted that we set COPA in the white matter of all tissue types to be 0 before regression, as there are only glia cells and a few interstitial neurons 64 in the white matter, which won't contribute to scattering significantly. Leaving it unchanged, however, will bias the regression. In the multivariable regression, we also assumed that k 1 and b to be the same in grey and white matter, as we assume that the contribution from myelin and other extracellular components behave similarly, and thus use the same coefficients k 1 and b to describe them jointly.
Lastly, we also investigated the origins of contrast in the µ ′ b map and ratio of µ ′ b /µ s . Four linear models were used, namely, fitting µ ′ b as a linear function of Gallyas OD and COPA, respectively, and fitting the ratio of µ ′ b /µ s as a linear function of Gallyas OD and COPA, respectively. Due to the lack of strong correlation of µ ′ b or ratio of µ ′ b /µ s with respect to Gallyas OD and COPA (see session 2.4 and supplementary figures), no multivariate investigation was further performed.

Data availability
The datasets and analysis code generated during and/or analysed during the current study are available from the corresponding author on reasonable request.