Toward a High-Resolution Reconstruction of 3D Nerve Fiber Architectures and Crossings in the Brain Using Light Scattering Measurements and Finite-Difference Time-Domain Simulations

Unraveling the structure and function of the brain requires a detailed knowledge about the neuronal connections, i.e., the spatial architecture of the nerve fibers. One of the most powerful histological methods to reconstruct the three-dimensional nerve fiber pathways is 3D-polarized light imaging (3D-PLI). The technique measures the birefringence of histological brain sections and derives the spatial fiber orientations of whole human brain sections with micrometer resolution. However, the technique yields only a single fiber orientation for each measured tissue voxel even if it is composed of fibers with different orientations, so that in-plane crossing fibers are misinterpreted as out-of-plane fibers. When generating a detailed model of the three-dimensional nerve fiber architecture in the brain, a correct detection and interpretation of nerve fiber crossings is crucial. Here, we show how light scattering in transmission microscopy measurements can be leveraged to identify nerve fiber crossings in 3D-PLI data and demonstrate that measurements of the scattering pattern can resolve the substructure of brain tissue like the crossing angles of the nerve fibers. For this purpose, we develop a simulation framework that permits the study of transmission microscopy measurements — in particular, light scattering — on large-scale complex fiber structures like brain tissue, using finite-difference time-domain (FDTD) simulations and high-performance computing. The simulations are used not only to model and explain experimental observations, but also to develop new analysis methods and measurement techniques. We demonstrate in various experimental studies on brain sections from different species (rodent, monkey, and human) and in FDTD simulations that the polarization-independent transmitted light intensity (transmittance) decreases significantly (by more than 50%) with an increasing out-of-plane angle of the nerve fibers and that it is mostly independent of the in-plane crossing angle. Hence, the transmittance can be used to distinguish regions with low fiber density and in-plane crossing fibers from regions with out-of-plane fibers, solving a major problem in 3D-PLI and allowing for a much better reconstruction of the complex nerve fiber architecture in the brain. Finally, we show that light scattering (oblique illumination) in the visible spectrum reveals the underlying structure of brain tissue like the crossing angle of the nerve fibers with micrometer resolution, enabling an even more detailed reconstruction of nerve fiber crossings in the brain and opening up new fields of research.


I. INTRODUCTION
The human brain consists of a huge network of nerve fibers: Around 100 billion nerve cells are connected to 10,000 other nerve cells on average [1][2][3][4][5].Understanding the structure and function of the brain remains a key challenge for neuroscience.To figure out how brain function emerges from its structural organization, it is necessary to study the neuronal connections, i. e., the three-dimensional nerve fiber architecture of the brain.Developing a detailed network model of the brain, the so-called connectome [6,7], reveals connected brain regions and helps to identify important nerve fiber connections, which is a prerequisite for brain surgery.It also serves as a reference for fiber tractography algorithms, improving the interpretation of clinical data obtained from diffusion magnetic resonance imaging (MRI) [8][9][10].Finally, the connectivity of the nerve fibers exposes pathological changes in the brain's tissue structure, allowing to study neuro-degenerative diseases like Alzheimer's or Parkinson's disease and to develop new treatments and tools for improved diagnostics.To visualize and derive brain tissue properties and organization principles, light-microscopy techniques are widely used [11][12][13].
In this paper, we study how the scattering of light can be used to improve the interpretation of light microscopy images from fibrous tissue samples like brain tissue.For this purpose, we have developed a simulation framework that models light scattering in transmission microscopy measurements and allows an improved interpretation of the measured data by delivering 3D information about the underlying fiber structure.
Technological progress and new advances in tissue preparation and labeling have enabled the development of techniques that reveal the 3D nerve fiber architecture in both living and post-mortem brains [14], such as Optical Coherence Tomography [15][16][17], Micro-Optical Sec-tioning Tomography [18], Light-Sheet Microscopy [19][20][21][22][23] or Two-Photon Fluorescence Microscopy (TPFM) [23][24][25][26][27].While most of these techniques are limited to small sample sizes, the neuroimaging technique 3D Polarized Light Imaging (3D-PLI) [28,29] allows to study the nerve fiber architecture of whole post-mortem brains with microscopic resolution: unstained histological brain sections are illuminated by polarized light and the birefringence caused by the highly anisotropic structure of the nerve fibers is measured, thus revealing their spatial orientations [30].
Such neuroimaging techniques require special instruments that are not available in each laboratory.Many laboratories are equipped with simple transmission microscopes which only extract 2D information of the investigated tissue structures, making it difficult to reconstruct complex nerve fiber architectures across several brain sections.Here, we show that conventional bright-field transmission microscopy measurements can be used to obtain information about the 3D fiber structure of a sample, without need to change the experimental setup or to repeat measurements: the polarization-independent intensity of light that is transmitted through the sample (transmittance) depends on the orientation of the fibers with respect to the light beam, i. e., the transmittance of a brain section provides information about the out-of-plane inclination angle of the enclosed nerve fibers.These findings can also be transferred to other biological and non-biological samples with comparable fibrous structures, e. g., muscle fibers, collagen, or artificial fibers.
The correct reconstruction of nerve fiber crossings is a major challenge for many neuroimaging techniques and a prerequisite for the correct interpretation of clinical MRI data, allowing for better diagnostics and treatments of neuro-degenerative diseases.In standard 3D-PLI measurements, brain regions with in-plane crossing fibers cannot be distinguished from regions with out-ofplane fibers or regions with low fiber densities because they all yield small birefringence signals.So far, only the strength and phase of the birefringence signal are used to derive the spatial fiber orientations.The transmittance, which is the average value of the signal, has not been used for this purpose.Here, we demonstrate that the transmittance cannot only be used to gain information about the out-of-plane inclination of the fibers, but also to classify these regions and to identify crossing fibers that cannot unambiguously be determined with standard polarization microscopy techniques.
The transmittance is a measure of how much the light is attenuated when it passes through the brain tissue, i. e., it depends on tissue absorption as well as scattering of light.As the absorption coefficient of brain matter is small (less than 0.1 mm −1 [31,32]), the measured transmittance is expected to be mainly influenced by scattering.To study such complex light-tissue interactions at the microscopic level, we employed finitedifference time-domain (FDTD) simulations to compute the propagation of the light wave through the brain tissue sample [33][34][35][36][37]. FDTD simulations are a proven tool for studying for example light scattering in lithography applications [38][39][40] or nanostructures [41][42][43].They have also been applied to investigate microscopy measurements of non-biological and biological tissue samples [38,44,45], but not yet to brain tissue.One reason might be that simulations of tissue samples with dimensions of several micrometers are computationally very intense because the mesh size in the simulation needs to be much smaller than the wavelength.To still enable the investigation of larger samples like brain tissue, we used high-performance computing and a simplified simulation model for the optics of the imaging system and the inner structure of the nerve fibers.The developed simulation framework can easily be adapted to microscopy techniques with different optics (wavelength, polarization of light, numerical aperture, etc.) and tissue samples with comparable fibrous structures.
The paper is divided in an experimental and a simulation part: In Sec.II, we evaluate experimental data and develop techniques to obtain structural 3D information from transmittance images of samples with unknown substructure: we study measurement results obtained from brain sections of different species (rodent, monkey, human) and show that the transmittance decreases significantly (by more than 50%) with increasing out-ofplane inclination angle of the nerve fibers, using both 3D-PLI and TPFM measurements to access the fiber inclination.In Sec.III, we introduce the FDTD simulation framework and present the simulation results for artificial nerve fiber configurations with different inclination and crossing angles.The simulations show that the decrease in transmittance is mainly caused by isotropic light scattering and by the finite numerical aperture of the imaging system.The in-plane crossing angle of the fibers has no impact on the transmittance and can be determined from the respective scattering pattern.In Sec.IV, we combine the simulation results with experimental data and show that the transmittance can be used to distinguish between regions with in-plane crossing fibers, regions with out-of-plane fibers, and regions with low fiber density by combining the transmittance and the strength of the measured birefringence signal.

II. EXPERIMENTAL STUDIES
In this section, various experimental studies are presented that evaluate how the transmittance of brain sections depends on the out-of-plane inclination angle of the enclosed nerve fibers.For our studies, we mostly used 3D-PLI measurements as they provide both the transmittance and the three-dimensional nerve fiber orientations independently from each other [28,29].As both scattering and absorption contribute to the attenuation (transmittance) of light, differences in the transmittance might not only be caused by different fiber inclinations, but also by a different tissue composition or density.To investigate the inclination dependence independently from tissue composition or preparation, we performed our studies on different species, subjects, and brain sections.
The studies were conducted on healthy brains from mice, rats, vervet monkeys, and humans.All brains were obtained directly after death in accordance with legal and ethical requirements.The brains were deeply frozen, cut into sections of 60 µm thickness, embedded in a solution of 20% glycerin, and cover-slipped.A detailed description of the brain preparation can be found in Appendix A 1. The brain sections were measured with 3D-PLI, using a polarimeter with a numerical aperture of 0.15 and an object-space resolution of about 1.33 µm / px [28,29].The polarimeter consists of an LED light source, a rotating linear polarizer, a specimen stage containing the brain section, a circular polarization analyzer, and a camera which records the transmitted light intensity for different rotation angles {0 • , 10 • , . . ., 170 • } of the polarizer.More information about the 3D-PLI measurement can be found in Appendix A 2. The signals provide information about the spatial orientations of the nerve fibers.The amplitude of the signal (retardation | sin δ|) is related to the birefringence of the brain section and served here as a measure of the out-of-plane fiber inclination angle α, using δ ∝ cos 2 α [30].The transmittance was computed from the same signal without additional measurements by averaging the measured light intensities over all rotation angles.(For an ideal system, conventional transmission microscopy with unpolarized light should yield the same result.)To only consider effects caused by the brain tissue, the resulting transmittance values were normalized for each image pixel by the average transmitted light intensity without sample, yielding the normalized transmittance I T,N .

A. Transmittance of flat and steep nerve fibers
In order to study how the transmittance depends on the out-of-plane inclination angle of the nerve fibers, the substructure of the investigated brain section needs to be taken into account.In large anatomical structures with densely packed nerve fibers, however, the inclination angles of the nerve fibers cannot be exactly determined by 3D-PLI or TPFM measurements.Instead, we compared the transmittance of flat nerve fibers (with inclination angles α < 45 • ) to the transmittance of steep nerve fibers (α > 45 • ) by investigating sections from brains that were cut along mutually orthogonal anatomical planes: one brain was cut along the coronal plane (dividing the brain into back and front), the other brain was cut along the sagittal plane (dividing the brain into left and right).As the sagittal plane is oriented perpendicular to the coronal plane, the transmittance of the same brain region can be evaluated for flat nerve fibers in one section plane and for steep nerve fibers in the other section plane.Since different brain sections from different specimens might not be comparable due to differences in the tissue structure, we only compared the transmittance values within the same brain section.
Figure 1 shows the normalized transmittance images I T,N of coronal and sagittal sections from rat and vervet monkey brains.The coronal (sagittal) section planes are indicated by blue (red) lines in the respective other brain section for reference.Selected brain structures were identified according to rat [46][47][48] and vervet [49][50][51] brain atlases.The transmittance values were evaluated in brain regions that have a relatively homogeneous tissue composition and that include predominantly flat nerve fibers (areas surrounded in green) or steep nerve fibers (areas surrounded in yellow).
The approximate orientation of the nerve fibers is known from the anatomy of the rat and the vervet brain as described in the atlases, and was confirmed by the retardation images | sin δ| shown below the transmittance images in Fig. 1: regions with flat nerve fibers show larger retardation values than regions with steep nerve fibers.The mean transmittance values and the standard deviation for the evaluated green and yellow areas can be found in Tab.I.In regions with flat nerve fibers, the mean transmittance values are larger (I T,N ∈ [0.18, 0.33]) than in regions with steep nerve fibers (I T,N ∈ [0.08, 0.15]).A region with flat (steep) nerve fibers which shows large (small) transmittance values in one section plane (coronal or sagittal), shows the opposite behavior in the corresponding orthogonal section plane.The difference is especially large when comparing the transmittance values of the corpus callosum (a massive fiber tract connecting the two hemispheres of the brain) and the cingulum (a C-shaped fiber structure running mostly perpendicular to the corpus callosum).In the coronal brain sections, the fibers of the cingulum (cg) run mostly perpendicular to the section plane and have about 55-67% lower transmittance values than the fibers of the corpus callosum (cc) which lie mostly within the section plane.In the sagittal brain sections, the situation is exactly the opposite: the transmittance values of the corpus callosum are about 50-54% less than the transmittance values of the cingulum.Fibers in the rat and vervet monkey brains show a similar pattern.
As expected, images obtained from conventional bright-field transmission microscopy with unpolarized light show similar effects, not only in vervet but also in human brain sections (see Fig. plemental Material): regions with steep (out-of-plane) nerve fibers appear darker than regions with flat (inplane) nerve fibers.

B. 3D-reconstruction of transmittance images
So far, single brain sections from different specimens were compared to each other.To study the transmittance across several consecutive brain sections, the transmittance images of 234 coronal sections from the right hemisphere of a vervet monkey brain were registered onto each other using in-house developed software tools (see Appendix A 2). Figure 3 shows the 3D-reconstructed transmittance images along three orthogonal anatomical planes: coronal (a), sagittal (b), and horizontal (c), as well as a detail of the 3D-volume (d).The white arrows point to the sagittal stratum -a white matter structure with nerve fibers that are oriented mostly perpendicular to the image plane (along the z-direction), as can be seen in the sagittal and horizontal planes.The structure appears much darker in the transmittance images than the surrounding regions.Thus, the observation that steep nerve fibers show lower transmittance values than flat nerve fibers is consistent across several consecutive brain sections.

C. Transmittance contrast of nerve fiber bundles in mutually orthogonal planes
In the previous studies, we considered bulk tissue with densely packed nerve fibers where the fiber inclinations cannot be exactly determined.In regions with distinct fiber bundles, however, the inclination angles can be estimated by manually evaluating the course of the fiber bundles in different section planes.For this purpose, we selected a structure in the rat brain that contains several distinct nerve fiber bundles with different, well-defined inclination angles -the so-called caudate putamen.To estimate the inclination angles of the nerve fibers, we evaluated the course of the bundles in mutually orthogonal section planes (coronal, sagittal, horizontal), see Fig. 4(a).As the brain sections were obtained from different brains and might differ in tissue composition, the transmittance images cannot be directly compared to each other.To still enable a comparison between the transmittance of flat and steep nerve fibers, the transmittance values in regions with fibers were compared to the transmittance values in regions with surrounding tissue for each brain section, assuming that the transmittance of the surrounding tissue does not depend much on the choice of the section plane.
To separate the fiber bundles from the surrounding tissue, we used the image contrast of the retardation images (see yellow lines in Fig. 4(b)).Figure 4(c) shows the corresponding histograms of the transmittance evaluated in regions with nerve fibers (pink) and in regions with surrounding tissue (cyan).As the coronal brain section contains mostly steep nerve fibers which yield low retardation values, the image contrast in the retardation image is not large enough to separate the fibers from the surrounding tissue.Therefore, we computed the histogram for the whole caudate putamen (area surrounded by yellow line) and assumed that the peak with lower (larger) transmittance belongs to nerve fibers (surrounding tissue).
From the minimum and maximum peak transmittance values of the histograms (pink and cyan numbers in Fig. 4(c)), we computed the transmittance contrast C ≡ (I T,max − I T,min )/(I T,max + I T,min ) between fiber bundles and surrounding tissue.For flat fiber bundles in the sagittal and horizontal brain sections, this contrast is much lower (5 • ≤ α ≤ 60 • : C ≈ 14-20%) than for steep fiber bundles in the coronal brain section (45 • ≤ α ≤ 85 • : C ≈ 62%).Assuming that the transmittance of the surrounding tissue is mostly independent of the fiber orientation, this demonstrates again that the transmittance values for steep nerve fibers are significantly lower than for flat nerve fibers.I T,N ) for the selected regions with nerve fibers (pink) and with surrounding tissue (cyan) in the caudate putamen.For the coronal brain section, the retardation image cannot be used to separate the fibers from the surrounding tissue because the fibers are oriented almost perpendicular to the image plane which leads to a small retardation signal and a small image contrast.Therefore, the histogram was computed over the whole selected region and the peak with lower (larger) transmittance was assumed to belong to fibers (surrounding tissue).The contrast values were computed from the respective peak values (numbers in pink and cyan) via: C = (max -min)/(max + min).The contrast for steep nerve fibers (coronal brain section) is much larger than for flat nerve fibers (sagittal and horizontal brain section).

D. Transmittance vs. inclination of nerve fiber bundles in TPFM images
The previous studies were only qualitative and the observed differences in the transmittance might also be caused by a different tissue composition, e. g., a different density of nerve fibers.To quantitatively study how the transmittance of a brain section depends on the inclination angles of the enclosed nerve fibers, the exact underlying fiber structure of the investigated brain section needs to be known.Therefore, we measured the caudate putamen of a coronal mouse brain section both with 3D-PLI and with TPFM to identify the inclination angles of the enclosed fiber bundles (see inset in Fig. 5(d) and Fig. S16 in the Supplemental Material).
The TPFM measurements were performed with a custom-made two-photon fluorescence microscope which achieves a resolution of 0.244 × 0.244 × 1 µm 3 and allows to perform an in-depth-scan of the brain section (see Appendix A 3 for more details).To obtain the inclination angles of the fiber bundles, the cross-sections of the bundles were determined in the first and the last slice of the TPFM image stack (cf.Fig. S16(d)).For each fiber bundle, the inclination angle was computed from the mid points of the corresponding cross-sections and from the thickness of the brain section (cf.Fig. S16(c)).
The fiber inclination and transmittance values were evaluated for 40 fiber bundles in the caudate putamen (see colored shapes in Fig. S16(b),(d)).The scatter plot in Fig. 5(d) shows the averaged transmittance values plotted against the determined fiber inclination angles.Although the values are broadly distributed, the scatter plot shows a clear tendency towards a decrease in transmittance with increasing fiber inclination angle.The values in orange belong to regions with lower fiber densities which might lead to overestimated transmittance values.However, the decrease in transmittance is also observed in regions with maximum fiber density (values in blue): while the mean transmittance values for flat nerve fibers (α < 50 • ) reach larger values (0.1 < I T,N < 0.2), the mean transmittance values for steep nerve fibers (α > 60 • ) are small (I T,N < 0.05).
All our experimental studies show that the transmittance of brain tissue decreases significantly (by more than 50%) with increasing out-of-plane inclination angle of the enclosed nerve fibers.

III. SIMULATION STUDIES
Although the experimental studies clearly show that the transmittance depends on the inclination angle of the nerve fibers, they do not provide enough information to describe this effect in detail.Based on the experimental results alone, it is not possible to make any predictions or draw conclusions for the interpretation of measured data.To model and better understand the observed transmittance effects, we performed numerical simulations on artificial nerve fiber configurations.This has the advantage that the exact underlying fiber structure, and thus the inclination angles of the nerve fibers, are known -also in bulk tissue with densely packed fibers.
As mentioned in Sec.I, the absorption coefficient of brain tissue is small so that the transmittance is expected to be mainly influenced by scattering.To study such complex light-matter interactions in microscopic detail, finite-difference time-domain (FDTD) algorithms are well suited.They discretize time and space, model the propagation of the light wave by approximating the spatial and temporal derivatives in Maxwell's curl equations by second-order central differences, and numerically compute the electromagnetic field components in space and time [33][34][35][36][37].As the mesh size of the spatial discretization needs to be much smaller than the wavelength (at most 25 nm), simulations of tissue samples with dimensions of several micrometers are computationally very intense.We have developed a simulation framework that allows for the first time to use FDTD simulations to study larger samples of fibrous tissue.For this purpose, we used a simplified simulation model for the brain tissue samples and the optics of the imaging system.We simulated the 3D-PLI measurement for various fiber configurations and evaluated the resulting transmittance values.
The artificial fiber configurations consist of about 700 fibers with uniformly distributed diameters between 1.0 µm and 1.6 µm and different fiber orientations.All fibers were generated in a volume of 30 × 30 × 30 µm 3 without intersections.The generation of the fiber configurations is described in Appendix B in more detail.Each fiber was represented by a simplified nerve fiber model, consisting of an inner axon with a constant radius and a surrounding myelin sheath with two layers, defined by different refractive indices (see Appendix C for motivation).
For the simulations of the 3D-PLI measurement, we used a conditionally stable FDTD algorithm to compute the propagation of the light wave through the tissue sample (artificial fiber configuration), described in more detail in Appendix D. The resulting electric field components were processed with analytical methods taking all optical components of the polarimeter into account, including the objective lens (with numerical aperture NA = 0.15) and the camera detector (see Appendix F).
The simulation studies were all performed for normally incident light with 550 nm wavelength and for the simulation parameters listed in Appendix E. One simulation run (volume of 30 × 30 × 30 µm 3 , mesh size of 25 nm) took about 8000 core hours on the supercomputer JUQUEEN (using an MPI grid of 16 × 16 × 16), allowing to perform many simulation runs with different parameters.The accuracy of the simulation results is discussed in Appendix G.

A. Simulated transmittance of inclined fibers
To better understand the inclination dependence of the transmittance that we observed in our experimental studies, we generated an artificial bundle of densely grown fibers (see Fig. 5(a) and Appendix B 1) for different inclination angles α = {0 • , 10 • , . . ., 90 • }, and computed the transmittance from a simulated 3D-PLI measurement.
Figure 5(b) shows the resulting scattering patterns (i.e., the intensity per wave vector angle θ k ) of light transmitted through the sample for inclination angles α = 0 • and 70 • .The white circles represent steps of ∆θ k = 10 • , from 0 • (center) to 90 • (outer circle).The transmittance images and scattering patterns for all inclination angles can be found in Fig. S17 in the Supplemental Material.
For flat fibers (α < 45 • ), the light is mostly scattered under angles perpendicular to the principal axis of the fiber bundle (i.e., along the y-axis).For intermediate inclination angles, the light is scattered more and more in the direction of the fibers (i.e., in the positive x-direction).For an inclination angle of 70 • , the light is broadly scattered in almost all directions (see Fig. 5(b)).
Due to the numerical aperture of the employed imaging system (NA ≈ 0.15), light scattered under angles larger than arcsin(NA) ≈ 8.6 • does not contribute to the measured transmittance images.To study the effect of the finite numerical aperture on the measured transmittance values, we simulated the imaging system without aperture (NA = 1) considering light scattered under all angles, and with aperture (NA = 0.15) considering only Figure 5(c) shows the resulting transmittance curves (mean values of the simulated transmittance images plotted against the inclination angles of the fiber bundle) for NA = 1 (orange curves) and NA = 0.15 (blue curves).The solid curves were obtained from the bundle of densely grown fibers (i) which has similar fiber orientations (the mode angle difference between the local fiber orientation vectors and the predominant orientation of the fiber bundle is less than 10 • ).The dashed curves were obtained for a bundle with broad fiber orientation distribution (ii) which contains many different fiber orientations (the mode angle difference is about 25 • , see Appendix B 2).To enable a better comparison between horizontal fiber bundles (α = 0 • ) and vertical fiber bundles (α = 90 • ), all curves were divided by the mean transmittance value of the horizontal bundle, respectively.Figure S18(a),(c) in the Supplemental Material shows the (normalized) transmittance curves in separate figures.
For NA = 1, the transmittance for steep fibers is similar to or even slightly larger than the transmittance for flat fibers.For NA = 0.15, the transmittance decreases significantly between α = 30 • and α = 70 • .This implies that the observed decrease in transmittance is caused by the finite numerical aperture of the imaging system: for steep fibers, the light is scattered almost uniformly in all possible directions (cf.Fig. 5(b) for α = 70 • ) so that the detected transmitted light intensity becomes minimal.In simulation studies with polarized light, we could show that the decrease in transmittance is independent of the direction of polarization (see Fig. 14(a) in Appendix G) which suggests that the decrease is caused by isotropic (not by anisotropic) scattering of light.
While the transmittance for the bundle of densely grown fibers becomes minimal at α = 70 • for NA = 0.15 (the transmittance is 90% less than the transmittance for the horizontal bundle) and the transmittance for vertical fibers (α = 90 • ) is only about 25% less than for horizontal fibers (α = 0 • ), the transmittance for the bundle with broad fiber orientation distribution decreases monotonically with increasing fiber inclination angle and becomes minimal for vertical fibers (the transmittance for vertical fibers is more than 80% less than for horizontal fibers).Due to the broad fiber orientation distribution, the vertical bundle contains many fibers with inclinations between 60 • and 70 • , which explains why the minimum transmittance is shifted to larger inclination angles.
Especially for the bundle with broad fiber orientation distribution, the simulated transmittance curves (Fig. 5(c)) show a similar behavior as the measured transmittance values in the scatter plot (Fig. 5(d)).
As mentioned in the beginning, the brain sections used for the 3D-PLI measurements are embedded in glycerin solution.With increasing time after this tissue embedding, we observed that the brain sections become more and more transparent, i. e., the transmittance increases (see Fig. 6).To enable optimal transmittance contrasts, the brain sections were therefore measured directly after the embedding and the simulations were performed assuming that the refractive indices of the nerve fibers correspond to given literature values (see Appendix C).With increasing time after the tissue embedding, the glycerin solution presumably soaks into the myelin sheaths which surround the axons.As the glycerin solution has a lower refractive index than the myelin lipids, the effective refractive index of the myelin sheaths is therefore expected to decrease.Figure S18(b) in the Supplemental Material shows the resulting transmittance curves for the bundle of densely grown fibers with a reduced myelin refractive index: the mean transmittance for steep fibers (α = 70-80 • ) is only about 30% less than the mean transmittance for horizontal fibers, and the absolute transmittance values become larger.This corresponds to the experimental observation that the transmittance increases with increasing time after the tissue embedding.

B. Simulated transmittance of crossing fibers
Our previous simulation studies have shown that the transmittance strongly depends on the out-of-plane fiber inclination angle.Hence, the transmittance could be used to distinguish out-of-plane fibers from in-plane crossing fibers, which both yield small birefringence signals and can to date not be distinguished in standard 3D-PLI measurements.To study this in more detail, we simulated the transmittance of horizontal (in-plane) crossing fibers for different crossing angles and compared the results to the transmittance of steep (out-of-plane) fibers.
The horizontal crossing fibers were generated as sep-  The scattering patterns of separate and interwoven crossing fiber bundles look similar for all crossing angles.The underlying fiber configuration, i. e., the crossing angle of the fiber bundles, is clearly visible in all scattering patterns.The mean transmittance values for the interwoven fiber bundles are up to 13% larger than those for the separate fiber bundles and in both cases mostly independent of the crossing angle.
Figure 7(b) shows the scattering pattern for three mutually orthogonal, interwoven fiber bundles.The fiber configuration is similar to the horizontal 90 • -crossing, interwoven fiber bundles, but one third of the fibers is oriented in the z-direction (see Appendix B 2 and Fig. 10(c)).This configuration leads to lower transmittance values than the horizontal crossing fibers.
Figure 7(c) shows the scattering pattern for a vertical fiber bundle (bundle with broad fiber orientation distribution and α = 90 • , cf.Fig. 5(c)(ii)).The scattering pattern looks clearly different from the scattering pattern of the horizontal crossing fibers and the mean transmittance is much lower.bundles, and the bundle with broad fiber orientation distribution for α = 90 • .
The transmittance curves of horizontal crossing fibers are similar for separate and interwoven fiber bundles.
While the mean transmittance of the separate crossing fibers corresponds more or less to the mean transmittance of the horizontal fiber bundle for χ = 0 • , the transmittance values of the interwoven crossing fibers slightly increase with increasing crossing angle (by max.11%).
For all simulated fiber bundles that contain vertical or steep fibers, the mean transmittance values (densely dotted lines) are more than 26% less than for the horizontal crossing fibers (solid lines).For interwoven crossing fibers, the transmittance value is reduced by more than one half when the horizontal crossing fibers are combined with a vertical fiber bundle (orthogonal bundles).For the vertical bundle with broad fiber orientation distribution and the steep bundle of densely grown fibers (with α = 70 • ), the difference between the transmittance values is especially large: the transmittance is about 80-90% less than for the horizontal crossing fibers.
Our simulations of crossing fiber bundles have shown that the transmittance for horizontal fibers is mostly independent of the crossing angle between the bundles and much larger than the transmittance for vertical fibers.This suggests that the transmittance values can be used to distinguish between horizontal crossing and vertical fibers in 3D-PLI measurements, and to detect vertical fibers within fiber crossings.

IV. COMBINATION OF EXPERIMENTAL AND SIMULATION STUDIES
In this section, we combine the results from the experimental studies and the simulation studies to develop a classification for brain regions with low birefringence signals that cannot be distinguished by 3D-PLI measurement.The simulations in Sec.III B have shown that the transmittance does not depend on the crossing angle between in-plane nerve fibers.First, we verify this prediction by investigating the optic chiasm of a hooded seal [52] -a region that contains fibers with crossing angles of around 90 • in the image plane (see Fig. 8).
While the retardation values in the region with crossing fibers (region B in orange) are broadly distributed (the birefringence signals of crossing fibers cancel out), the transmittance values in this region show a similar distribution as in a region with mostly parallel fibers (region A in blue), see histograms in Fig. 8(d).
The peak transmittance value of region B is slightly lower than in region A because the number of fibers in the crossing region (two crossing bundles) is larger than in the region with parallel fibers (one bundle).Thus, the transmittance depends on the tissue density, but not on the crossing angles between the nerve fibers -as predicted by the simulations in Sec.III B.
To demonstrate that the transmittance can be used to classify regions with small birefringence signals (i. e., small retardation values obtained from 3D-PLI measurements) into regions with in-plane crossing fibers, regions with steep fibers, and regions with low fiber density, the predictions obtained from the simulation studies in Sec.III were applied to experimental data (see Fig. 9).
As the transmittance depends on absorption, the region with maximum absorption was determined as a reference: The retardance δ of brain tissue increases with decreasing fiber inclination angle α and with increasing thickness d of birefringent tissue components (δ ∝ d ∆n cos 2 α, where ∆n is the birefringence of the tissue [28]).Assuming that a brain section contains all possible nerve fiber configurations, the region with maximum retardation signal | sin δ| max (orange ellipse in Fig. 9A) is therefore expected to contain mostly horizontal parallel fibers (α ≈ 0 • ) with a high fiber density (max.d ∆n) and thus to cause a maximum of absorption.Regions with even lower transmittance values are accordingly expected to contain steep (out-of-plane) fibers which increase the scattering and thus the attenuation of light.
By comparing the normalized transmittance values (I T,N ) of regions with small retardation values to the transmittance of the region with maximum retardation (I ref ≡ I T,N (| sin δ| max )), the regions can be classified into three categories (see Fig.For regions with slightly lower or larger transmittance values, an unambiguous classification is not possible.Provided that the region with maximum retardation has the largest tissue absorption, lower transmittance values can only be caused by steep fibers.Similar transmittance values, however, could also be caused by a small number of steep fibers, and larger transmittance values could be caused by a small number of in-plane crossing fibers (or a smaller number of steep fibers).A classification by means of retardation and transmittance values can therefore only serve as an indication of the underlying fiber configuration and should always be considered in addition to individual tissue characteristics.As the transmittance depends on the tissue preparation, the combined analysis of transmittance and retardation should only be performed section-wise.Brain atlases and 3Dreconstructed images (cf.Fig. 3) validate the classification of regions in Fig. 9.

V. DISCUSSION AND CONCLUSION
Conventional bright-field transmission microscopy measurements of fibrous tissue samples provide usually only 2D information about the underlying fiber architecture.This considerably limits the application of these methods in the analysis of three-dimensional fiber structures, e. g., when studying the complex architecture of nerve fibers in the brain.In this paper, we show both in experimental and simulation studies that the scattering of light enables a more enhanced interpretation of transmission microscopy images, providing additional structural information about the three-dimensional fiber architecture in brain tissue samples.First, we exploited various techniques to derive threedimensional structural information from brain tissue samples with unknown substructures.Our experimental studies on brain sections from different species (rodent, monkey, and human, see Sec.II) have shown that the polarization-independent transmitted light intensity (transmittance) significantly decreases with increasing out-of-plane inclination angle of the enclosed nerve fibers (by more than 50%).Using finite-difference time-domain (FDTD) simulations, we could successfully model this effect and show that the decrease in transmittance is caused by polarization-independent (isotropic) light scattering and by the finite numerical aperture of the imaging system (see Sec. III A).Polarizationdependent light scattering which leads to diattenuation (polarization-dependent attenuation of light) cannot explain the observed transmittance effect because the diattenuation of brain tissue was shown to be small [53,54].
Furthermore, we could use the FDTD simulations to explain the increasing transparency of brain tissue samples with increasing time after tissue embedding (see Fig. 6): when the embedding solution soaks into the surrounding myelin sheaths of the nerve fibers, this leads to an equalization of the effective refractive indices and thus to reduced scattering, which increases the transparency of the tissue.
Note that both measured and simulated transmittance values depend on many parameters such as the homogeneity of brain tissue or the density of nerve fibers, which are not easily accessible.As long as the exact underlying tissue structure is not known, our simulation results can therefore only serve as a qualitative prediction for the interpretation of measured data.
As the simulated samples are only characterized by their geometry and refractive indices, biological and non-biological samples with comparable fibrous structures (e. g., muscle fibers, collagen, artificial fibers) are expected to show similar transmittance effects.To increase the transmittance contrast between flat and steep fiber structures, the embedding solution should have a different refractive index than the fibers (cf.Fig. S18(a)-(b) in the Supplemental Material).
The observed transmittance effects are mostly independent of the polarization so that standard transmission microscopy techniques can be used to obtain 3D information about the underlying fiber configurations, without need to change the experimental setup or to repeat measurements.This greatly enhances the application of conventional bright-field transmission microscopy which is available in many laboratories and so far only used to gain 2D structural information.
When generating a detailed model of the nerve fiber architecture in the brain, the reconstruction of brain regions with crossing nerve fibers poses a major challenge.With 3D-Polarized Light Imaging (3D-PLI), it is generally not possible to distinguish brain regions with in-plane crossing fibers from regions with out-of-plane fibers or from regions with low fiber densities, because they all yield low birefringence signals.Using FDTD simulations, we could show that the transmittance of inplane fiber configurations does not depend on the crossing angle between the fibers (see Sec. III B).Applying the predictions of our simulation studies to experimen-tal data, we could demonstrate that a combined analysis of transmittance and retardation (strength of the birefringence signal) enables to distinguish between these regions (see Sec. IV).Our simulations also revealed that the transmittance can be used to detect out-of-plane fibers in regions with in-plane crossing fibers, which is not possible with current techniques.The combined analysis of transmittance and retardation images can also be applied to past 3D-PLI measurements in order to validate and -if necessary -correct the reconstructed fiber orientations.
Apart from the fiber inclination, the transmitted light intensity reveals much more information about the underlying tissue structure when studying the exact pattern of the scattered light.Our simulation studies in Sec.III B have shown that the scattering pattern can be used, for example, to identify the crossing angle of nerve fiber bundles, which is not easily accessible with current measurements.How the scattering pattern is related to the exact underlying fiber structure and tissue homogeneity will be addressed in future studies.Major features of the scattering pattern (like the fiber crossing angle) can be determined by simply placing an aperture between light source and sample and measuring the transmitted light intensity for different positions of the aperture.
The FDTD simulations proved to be a valuable and reliable tool in many aspects: they allow to better understand the interaction of polarized light with brain tissue, to find explanations for the observed transmittance effects, to make general predictions, and to improve the measurement procedure and analysis.In contrast to previous top-down simulation approaches of 3D-PLI that model the optical properties of the nerve fibers by series of Jones matrices (simPLI [30,52]), the FDTD simulations solve Maxwell's equations and allow to model more complex effects like the scattering of light, but they require much more computing time.
Most nerve fibers in the brain are surrounded by a so-called myelin sheath, which consists of multiple layers with 3-5 nm thickness (see Appendix C).If the exact layered structure of the myelin sheath is modeled, the mesh size in the simulations can be at most 3 nm.In this case, the simulation of a single nerve fiber with 1 µm diameter consumes almost 290 000 core hours (see Appendix G 1).To enable the simulation of larger tissue samples with various simulation parameters, we developed a simplified nerve fiber model with double myelin layers and a simplified model for the imaging system (the incoherent and diffusive light source was modeled by monochromatic light with normal incidence).We could show that these simplified models still reproduce the observed transmittance effects and that our results are not sensitive to small changes in the simulation parameters (see Appendix G and [55]), so our model is a good compromise between accuracy and computing time.The developed simulation framework can easily be adapted to microscopy techniques with different optics (numerical aperture, wavelength, polarization, etc.) and to other species and tissue types.Further brain tissue components like glial cells can easily be added to the simulation model [56].
In summary, we have developed and successfully applied a versatile simulation framework for transmission microscopy measurements of fibrous tissue samples that allows to study light scattering in larger samples like brain tissue, using finite-difference time-domain simulations.We have demonstrated both in experimental and simulation studies on various brain tissue samples that the polarization-independent transmitted light intensity (transmittance) provides information about the 3D orientation of the enclosed fibers, allowing to use simple bright-field transmission microscopy to study three-dimensional fiber structures.Finally, we could show that the transmittance can be used to classify brain regions with low birefringence signals without changing the experimental setup or repeating measurements.This enables a more enhanced interpretation of threedimensional nerve fiber architectures in the brain.old), as well as on brain sections from vervet monkeys (African green monkey: Chlorocebus aethiops sabaeus, male, between one and two years old), rats (Wistar, male, three months old), mice (C57BL/6, male, six months old), and a hooded seal [52].All animal procedures were approved by the institutional animal welfare committee at Forschungszentrum J ülich GmbH, Germany, and are in accordance with European Union (National Institutes of Health) guidelines for the use and care of laboratory animals.The human brain was acquired in accordance with the local ethic committee of the University of Rostock, Germany.A written informed consent of the subject is available.
The brains were removed from the skull within 24 hours after death, immersed in a buffered solution of 4% formaldehyde for several weeks, immersed for several days in solutions of 10% and 20% glycerin combined with 2% Dimethyl sulfoxide for cryoprotection, dipped in cooled isopentane for several minutes, and deeply frozen.The frozen brains were cut with a cryostat microtome (Leica Microsystems, Germany) at a temperature of -30 • C into sections of 60 µm.The brain sections were mounted on cooled glass slides, embedded in 20% glycerin solution, covered by a cover glass, sealed with lacquer, and weighted for several hours to prevent the development of air bubbles.The sections were measured one day after embedding to obtain optimal transmittance images.

3D-PLI measurement
The 3D-PLI measurements were performed with a high-resolution Polarizing Microscope (PM) manufactured by Taorad GmbH, Germany.The microscope has been used in previous 3D-PLI studies to measure the three-dimensional nerve fiber orientations at high resolution [28,29,59,60].The light source consists of a single white LED (IntraLED 2020+ operated at 24 W) with integrated K öhler illumination and a bandpass filter, generating a wavelength spectrum λ = (550 ± 5) nm.Further components are a rotatable linear polarizer, a specimen stage, a circular analyzer (quarter-wave retarder combined with linear polarizer), and a CCD camera (monochrome RETIGA-4000R camera by QImaging with Kodak KAI-04022-ABA image sensor) which records an image for each rotation angle ρ = {0 • , 10 • , . . ., 170 • } of the polarizer, yielding a series of 18 images.The microscope is equipped with a motorized specimen stage (Märzhäuser, Germany) which performs a translational scan of the brain section in tiles of 2.7 × 2.7 mm 2 .To allow for stitching, the tiles were measured with an overlap of 30% on all sides.The objective lens (Nikon TL Plan Fluor EPI P 5x) has a 5× magnification and a numerical aperture of 0.15.The resolution in object space is about 1.33 µm / px.
The transmittance and retardation images in Secs.II and IV were computed as described in Axer et al. [28,29] by performing a discrete harmonic Fourier anal-ysis on the measured light intensities I(ρ) per image pixel: I(ρ) = a 0 + a 2 cos(2ρ) + b 2 sin(2ρ).The transmittance I T corresponds to the average over all 18 images and was computed from the Fourier coefficient of order zero (I T = 2 a 0 ), the retardation | sin δ| corresponds to the amplitude of the intensity signal and was computed from the Fourier coefficients of order zero and two | sin δ| = (a 2 2 + b 2 2 ) 1/2 /a 0 , where δ is the phase shift induced by the birefringent brain tissue.The transmittance images were normalized by the transmittance image measured without sample, yielding normalized transmittance images (I T,N ).
Images of several consecutive brain sections (see Fig. 3) were registered onto each other using in-house developed software tools based on the software packages ITK, elastix, and ANTs [61][62][63][64][65] which perform linear and non-linear transformations.As undistorted reference volume, aligned blockface images were used: a picture of the brain block surface (blockface image) was taken every time before sectioning, a pattern of ARTag markers [66] was used to determine the position of the brain block in two-dimensional space [67].

TPFM measurement
The TPFM measurements were performed with a custom-made two-photon fluorescence microscope [23,68] at the European Laboratory for Non-Linear Spectroscopy (LENS), University of Florence, Italy.The microscope is equipped with a mode-locked titanium-sapphire laser with a wavelength of 800 nm which is coupled into a scanning system based on a pair of galvanometric mirrors.The laser is focused onto the sample by a waterimmersion 25 × objective lens (LD LCI Plan-Apochromat 25x/0.8Imm Corr DIC M27).The lateral displacement of the sample was realized by a motorized xy-stage (enabling tile-wise scanning of the sample).The axial displacement (along the z-axis) was realized by a closedloop piezoelectric stage.The fluorescence signals were collected by two photomultiplier tubes, enabling to detect red and green fluorescence.The setup achieves a resolution of 0.244 × 0.244 × 1 µm 3 .The sample was measured in tiles of 250 × 250 µm 2 , with an overlap of 10% to allow for stitching.

Bright-field transmission microscopy
The bright-field transmission microscopy images were obtained from ZEISS Axio Imager Vario.The microscope is equipped with a white microLED which emits unpolarized light with wavelengths between 400 nm and 750 nm.The objective lens (Plan Apochromat 5x) has a 5× magnification and a numerical aperture of 0.16.The resolution in object space is about 0.91 µm / px.

Appendix B: Generation of artificial fiber configurations 1. Densely grown fiber bundle
The bundle of densely grown fibers (see Fig. 5(a)) was generated by in-house developed software: N = 700 circles with uniformly distributed diameters (d ∈ [1.0, 1.6] µm) were randomly uniformly placed in the xy-plane (in an area of 45 × 30 µm 2 ).The circles were initialized with a random speed (max.0.1 µm displacement per step) and collided with each other (assuming elastic collision with particle mass r 2 ) until a solution was reached without collision in the xy-plane.To obtain well-distributed fibers, the previous step was repeated 250 times before the positions of the circles were stored.To obtain a 3D fiber volume, the circle positions were stored while incrementing the z-position by 1 µm per step.To generate fiber bundles with different inclination angles, the resulting bundle of densely grown fibers was rotated around the y-axis with respect to the center position and cropped to a volume of 30 × 30 × 30 µm 3 .To prevent fibers from touching each other after discretization, all fiber diameters were reduced by 5%.In the resulting fiber bundle, about 60% of the volume is filled with fibers.

Inhomogeneous fiber bundles
Inhomogeneous fiber bundles, like the bundle with broad fiber orientation distribution (Fig. 5(c)(ii)) or crossing fibers (Fig. 7(a),(b)), were generated by in-house developed software [69] which allows collision control in 3D.Starting from well-distributed straight fibers with N = 700 and d ∈ [1.0, 1.6] µm (obtained after 250 steps as described in the previous section), the fibers were divided iteratively into segments of 2-5 µm and assigned a random displacement in the x-, y-, and z-direction.The resulting fiber segments were split or merged until the length of each segment was again between 2-5 µm, ensuring that the maximum angle between adjacent segments was less than 20 • .When a collision between two segments was detected, the segments were exposed to a small repelling force and the previous step was repeated until no more collisions were detected.To prevent fibers from touching each other, all fiber diameters were reduced by 5%.The resulting fiber bundle was cropped to a volume of 30 × 30 × 30 µm 3 .The fiber bundles were generated from different configurations of straight fibers and different random displacements: • Bundle with broad fiber orientation distribution (Fig. 5(c)(ii)): The fiber bundle was generated from a bundle of straight horizontal fibers in the xdirection and a maximum random displacement of 10 µm.In the resulting fiber bundle, about 33% of the volume is filled with fibers.To generate fiber bundles with different inclination angles, the re-sulting bundle was rotated around the y-axis with respect to the center position.
• Separate crossing fiber bundles (Fig. 7(a)): The bundle of straight horizontal fibers in the x-direction was divided in an upper and a lower bundle of thickness z/2, respectively.The upper bundle was rotated around the z-axis about the center position by an angle +χ/2, the lower bundle was rotated by an angle −χ/2, resulting in two separate bundles with crossing angle χ (cf.Fig. 10a).The resulting fibers were used as input for the algorithm with a maximum displacement of 1 µm.Depending on the crossing angle of the resulting fiber bundle, between 40-50% of the volume is filled with fibers.
• Interwoven crossing fiber bundles (Fig. 7(a)): Each fiber layer in the z-direction of the straight horizontal fiber bundle (oriented in the x-direction) was rotated alternately by ±χ/2 (cf.Fig. 10(b)).The resulting fibers were used as input for the algorithm with a maximum displacement of 1 µm.Depending on the crossing angle of the resulting fiber bundles, between 40-50% of the volume is filled with fibers.
• Mutually orthogonal, interwoven fiber bundles (Fig. 7(b)): The straight horizontal fiber bundle (oriented in the x-direction) was divided in three types of alternating layers: one layer was rotated + 45 • around the z-axis, one − 45 • around the zaxis, and one +90 • around the y-axis, yielding two horizontal fiber bundles in the xy-plane and one vertical fiber bundle oriented along the z-axis.The resulting fibers were used as input for the algorithm with a maximum displacement of 1 µm.In the resulting fiber bundle (cf.Fig. 10(c)), ca.32% of the volume is filled with fibers.For the simulation studies in Sec.III, a simplified model was used to represent the myelin sheath (see Fig. 11(d)): Each cell layer (two lipid bilayers with separating cytoplasm) was considered as one myelin layer with an effective refractive index n m = 1.47 (blue), the extracellular space was considered to be filled with glycerin solution (glycerin layer) with a refractive index n g = 1.37 (yellow).Assuming that the extracellular space increases when being embedded in glycerin, the myelin and glycerin layers were assumed to contribute 3/4 and 1/4 to the overall myelin sheath thickness t sheath , respectively.The refractive index of the cytoplasmic layer was neglected in this model.
The myelin sheath thickness contributes approximately one third to the overall fiber radius r [76].Hence, the myelin sheath thickness was chosen to be t sheath = 0.35 r and the radius of the inner axon r ax = 0.65 r.The refractive index of the axon (green) was chosen to correspond to the refractive index of cytoplasm (n ax = 1.35).The myelin sheath was modeled as double myelin layers with thickness t m = (3/7) t sheath each and a single glycerin layer with thickness t g = (1/7) t sheath separating the myelin layers.Interruptions of the myelin sheath (nodes of Ranvier) and the small space between axon and myelin sheath (periaxonal space [71]) were neglected in this model.

Appendix D: FDTD algorithm
The propagation of the polarized light wave through the brain tissue sample (nerve fiber configuration) was simulated by a massively parallel 3D Maxwell solver based on a conditionally stable finite-difference timedomain (FDTD) algorithm [33].The algorithm computes the electromagnetic field components numerically by discretizing space and time and approximating Maxwell's curl equations by finite differences: The discretization is realized with a cubic Yee grid [77] (each electric field component is surrounded by four magnetic field components and vice versa) and a leapfrog time-stepping scheme.The spatial and temporal derivatives in Maxwell's curl equations are approximated by second-order central differences.For more details, see Menzel et al. [34].
The simulations were performed with the software TDME3D TM [35,78] -a massively parallel threedimensional FDTD Maxwell Solver, Copyright EMBD (European Marketing and Business Development BVBA).The software solves Maxwell's equations for arbitraryshaped objects that are illuminated by arbitrary incident plane waves and that consist of linear, isotropic, lossy materials with known permeability, permittivity, and conductivity.For the FDTD simulations, a combined algorithmic approach was used: In free space, Yee's algorithm was applied.To compute the interaction of the light with brain tissue, an unconditionally stable Lie-  Trotter-Suzuki product formula approach was used.This results in a computationally efficient but conditionally stable algorithm.For more information, see De Raedt [79].The simulations were performed on the supercomputer JUQUEEN [58] at Forschungszentrum J ülich GmbH, Germany.

Appendix E: Simulation parameters
Table II lists the parameters that were used for the simulation studies in Sec.III.
All fiber configurations were generated in a volume of 30 × 30 × 30 µm 3 .As described in Appendix C, each fiber was modeled by an inner axon and a surrounding myelin sheath with two layers and different refractive indices (see Fig. 11).The surrounding medium was assumed to be homogeneous with a refractive index n surr = n g = 1.37, which corresponds to the refractive index of gray brain matter as well as to the refractive index of the surrounding glycerin solution.To account for the fact that the brain sections are embedded in glycerin solution (see Appendix A 1), 0.5 µm thick layers of glycerin solution (with refractive index n g = 1.37) were added at the bottom and on top of the sample, yielding a medium with dimensions 30 × 30 × 31 µm 3 .The dimensions of the simulation box were chosen to be 30 × 30 × 35 µm 3 to leave some space for light source and detection planes.The simulation volume was surrounded by uniaxial perfectly matched layer (UPML) absorbing boundaries of 1 µm thickness, thick enough to prevent light from being reflected back into the simulation volume.The different components of the sample were simulated as dielectrics with real refractive indices (as described in Appendix C).Absorption was neglected because the absorption coefficients of brain tissue are small [31,32].
The simulation studies were performed for a duration of 200 periods and a Courant factor of 0.8.The Yee mesh size was chosen to be ∆ = 25 nm.This mesh size is just large enough to account for the double myelin layers of the nerve fiber model (see Appendix C: the glycerin layer for fibers with 1 µm diameter is 25 nm thick).
The light source was modeled as plane monochromatic wave.The simulation studies were performed for normally incident and coherent light with left-handed circular polarization and a wavelength of 550 nm (corresponding to the peak wavelength of the employed light source, see Appendix A 2).Using an MPI grid of 16 × 16 × 16 on JUQUEEN, each simulation run (i.e., the calculation of one configuration, one wavelength, and x × y × z = 30 × 30 × 30 µm 3 fiber radius: r ∼ 0.5 µm axon: r ax = 0.65 r, n ax = 1.35 myelin sheath: double myelin layers: t m = 3 7 t sheath , n m = 1.47 single glycerin layer: t g = 1 7 t sheath , n g = 1.37 one angle of incidence) consumed between 7000-8000 core hours, required a minimum memory between 260-360 GB, and lasted between 1:45-2:00 hours.

Appendix F: Computation of the transmitted light intensities
Figure 12 shows how the 3D-PLI measurement was modeled by means of FDTD simulations.For the simulations, a mathematically equivalent polarimetric setup of the employed microscope was considered in which the sample is illuminated by (left-handed) circularly polarized light and analyzed by a rotating linear polarizer (analyzer).
The computation of the transmitted light intensities consists of several steps: 1.) Maxwell Solver: After passing the polarizing filters in front of the sample (see Fig. 12 on the left), the light wave is left-handed circularly polarized.The propagation of the light wave through the sample was com-puted by TDME3D as described in Appendix D. The resulting light wave is represented by a superposition of monochromatic plane waves with different wave vectors k and real amplitudes E 0,k : where r and t are the spatial and temporal coordinates, ω is the angular frequency, φ is the phase, and A k and B k are defined as: A k = E 0,k cos φ and B k = E 0,k sin φ.
Note that every index k denotes a different wave vector k and is not related to the wave number k = 2π/λ (the wavelength of the transmitted light waves is the same as for the ingoing light wave).

2.) Yee shift:
Before further processing, the electromagnetic field components were shifted in the x,y,z-direction to the middle of the corresponding Yee cell, respectively: where ∆x = ∆y = ∆z is the side length of the cubic Yee cell.
For each shift ∆ j in the direction j = {x, y, z}, the vector components A k,i and B k,i were recomputed as follows: (F7) After performing the shifts specified in Eq. (F3) to (F5), the resulting field vector is given by: (F8)

3.) Scattering pattern:
To study how much light is scattered under a certain angle (wave vector k), the scattering pattern was computed, i. e., the intensity per wave vector normalized by the ingoing light intensity (I 0 ) per image pixel (px):

Measurement Simulation
• sample size: (mm²-cm²) x 60 μmm • structures > 10 nm • refractive index: complex, (an)isotropic • sample size: < (100² μmm²) x 100 μmm • mesh size: 3-25 nm • refractive index: complex, isotropic • LED light source • wavelength: • angle of incidence: φ = 0°-360°, θ < 3°• plane wave • wavelength: fixed λ = {545, …, 555 nm} • coherent light • angle of incidence: fixed φ = {0°, 60°,…}, θ = {0°,…, 3°} • partial polarization (imperfect polarizers) • completely polarized light 12. Modeling of the 3D-PLI measurement.The figure and table on the left-hand side show the optical components of the polarimeter (the order of the polarizing filters is different than in the measurement, but the setup is mathematically equivalent): light source (green), polarizer/retarder (dark gray), sample (light gray), objective lens/detector/camera (blue).The table and figure on the right-hand side show how the optical elements were modeled by FDTD simulations: The incoherent and diffusive light source (LED) with peak wavelength λ and full-width at half-maximum (FWHM) was modeled by performing several simulation runs with plane waves that have different wavelengths (λ) and angles of incidence (ϕ, θ).The modeled light source emits coherent light that is circularly polarized.The tissue sample was represented by an artificial fiber architecture, the rotating analyzer by a rotated Jones matrix (with rotation matrix R(ρ)).The numerical aperture (NA) of the imaging system was modeled by considering only wave vector angles θ k < arcsin(NA).The spherical microlenses of the camera detector were modeled by performing a moving average over the area of the microlens with radius r = 1.33 µm / 2.

5.) Objective lens:
The objective lens was assumed to be ideal and both specimen and detector were assumed to lie within the corresponding focal planes of the lens.Thus, the propagation of the electromagnetic wave between sample and detector was assumed to be free and Ẽk (r, t, ρ) was evaluated at the z-position of the detection plane behind the sample (defined as z = 0): To account for the numerical aperture (NA) of the objective lens, only k-vectors were processed that fulfill: The employed imaging system has a numerical aperture of about 0.15, so only k-vectors with angles θ k ≤ 8.6 • were used for processing.

6.) Detector microlenses:
The camera sensor contains an array of spherical microlenses which bundle the light onto subjacent photodiodes for each image pixel.Assuming perfect microlenses and photodiodes that are completely covered by one microlens, respectively, the microlenses were modeled by applying a moving average over the area of the microlens.Instead of taking the magnification and the physical size of the microlenses into account, the microlenses were modeled with a diameter of 2 r 0 = 1.33 µm corresponding to the pixel size of the microscope in object space: To obtain the full image information (independent of the detector pixel position), no rasterizing was applied.

7.) Intensity:
In principle, the intensity detected by the camera sensor depends on the angle of incidence of the incident light: I cos θ k .As the numerical aperture is sufficiently small (NA = sin θ k ≈ 0.15 ⇔ cos θ k > 0.9886), the angle dependence was neglected, which enables to represent the intensity I(r, ρ) as Fourier series in ρ, as described below.
With this assumption, the light intensity recorded by the camera is given by the absolute squared value of the electric field vector.To compute the intensity at a certain point r in the image plane, the electric field vectors were summed over k and averaged over time: where FT −1 denotes the inverse discrete Fourier transform: The discrete Fourier transform (FT) is defined analogously.
To save computing time, the convolution in Eq. (F16) was replaced by a multiplication, making use of the convolution theorem: where the function J 1 (x) is the Bessel function of the first kind of order one, with k xy ≡ k 2 x + k 2 y and r 0 = 0.665 µm.
To simplify notation, the following abbreviations are defined: The intensity is then given by: The x-and y-components of the electric field vector Ẽk (r, t, ρ) behind the rotating analyzer were computed from E k (r, t) = A k cos(k • r − ωt) − B k sin(k • r − ωt) according to Eq. (F11).As the equation is linear in the xand y-components of E k (r, t), the x-and y-components of {A k , B k , E k } are transformed to { Ãk (ρ), Bk (ρ), Ẽk (ρ)} according to the same equation.As the Fourier transform is independent from ρ, Eq. (F11) also holds for the xand y-components of E (r) and Ẽ(r, ρ), yielding Fourier coefficients of order zero and two: where trigonometric identities have been used: cos 2 x = where e m (r) and f m (r) are functions of the inverse discrete Fourier transforms: Thus, the transmitted light intensity I(r, ρ) can be written in terms of a Fourier series: where the Fourier coefficients a m (r) and b m (r) are computed from the six inverse discrete Fourier transforms defined above: E x (r), E y (r), X x (r), X y (r), Y x (r), Y y (r).
For non-normally incident light (k x 0 or k y 0), the transmitted light intensity contains Fourier coefficients of order four (cf.Eq. (F27)).
Using Eq. (F30), the light intensity was computed for arbitrary rotation angles ρ and normalized by the ingoing light intensity per image pixel: In the experiment, the measured light intensities are normalized by the light intensities measured without specimen to compensate for filter inhomogeneities.This image calibration could be modeled by performing an additional simulation run without sample.To save computing time, the simulated light intensities were simply normalized by I 0 (without considering the imaging system) and only relative values were used for the comparison between measured and simulated light intensities.The Fourier coefficient of order zero a 0,N (r) obtained from the simulated normalized transmitted light intensity I N (r, ρ) was used to compute the simulated transmittance images I T,N (r): Figure S19 in the Supplemental Material summarizes the most important steps of computing the transmitted light intensities for 3D-PLI simulations.The computation was carried out in Python (version 2.7.6) using the NumPy package (version 1.12.1)[82,83].To obtain the intensity at a certain pixel position (x,y), the inverse discrete Fourier transform was computed in two dimensions by means of the Fast Fourier Transform [84].To enable an efficient use of the FFT, the number of grid points in x and y (N x and N y ) were set to be a multiple of two: (F35)

Appendix G: Error estimation of simulation results
When modeling the optical components of the imaging system, the limitations of the simulation software need to be taken into account: the simulated light wave is completely polarized and coherent, the materials are characterized by isotropic refractive indices, and size and resolution of the simulated geometries are limited due to finite computing time.
Using completely polarized light for the simulations implies that the optical elements are assumed to be ideal (unpolarized light source, ideal polarizing filters, no polarization-sensitivity of the camera).For the employed polarizing microscope, these assumptions are reasonable because the optical components are of high quality.Moreover, the transmittance can be considered to be mostly independent from the polarization properties of the imaging system.
The simulation studies in Sec.III were performed for a reduced sample size (30 × 30 × 30 µm 3 ) and 200 periods.Simulations with larger sample sizes in x/y and more periods yielded similar results [85].
To further reduce computing time, the simulation studies were performed for a simplified nerve fiber model (axon surrounded by double myelin layers, cf.Fig. 11(e)), a Yee mesh size of 25 nm, and normally incident light with 550 nm wavelength.To estimate the accuracy of the simulation results, the transmittance images were simulated for different numbers of myelin layers L, different Yee mesh sizes ∆, different wavelengths λ, and different angles of incidence θ.To study the influence of one simulation parameter at once, only one simulation parameter was varied while all other simulation parameters were chosen as in Tab.II (with normally incident light and 550 nm wavelength).
To estimate the accuracy of the resulting transmittance images, the absolute relative difference between the mean values (ARDM) and the relative mean absolute difference (RMAD) between the images were computed: In this notation, the "image" refers to the transmittance image for which the absolute relative difference is computed (obtained, e. g., from simulations with different Yee mesh sizes).The "reference image" is the transmittance image used for comparison (obtained, e. g., from the simulation with minimum mesh size).The symbol represents the average over all image pixels.As the simulation studies mostly investigate the mean transmittance values, the ARDM is a direct measure for the accuracy of the simulation results, while the RMAD is a measure for the reliability of the ARDM as an error estimate.

Different numbers of myelin layers
To estimate the accuracy of the simplified nerve fiber model, a straight single fiber with reduced simulation volume (see Fig. 13(a)) was simulated for different numbers L of myelin layers with thickness t m (and L − 1 separating glycerin layers with thickness t g = t m /3) as well as for a realistic model of the myelin sheath consisting of 43 thin layers (see Fig. 13(b)).The Yee mesh size was chosen to be small enough to resolve all geometric features: For most samples, the mesh size was chosen to be one third of the glycerin layer thickness (∆ = t g /3).Fibers with two myelin layers (L = 2) were also simulated for larger mesh sizes (∆ = t g /2 = 12.5 nm and ∆ = t g = 25 nm).The realistic myelin sheath was simulated for ∆ = t g = 3 nm, consuming 288 358 core hours on JUQUEEN (using an MPI grid of 64 × 64 × 16).
Figure 13(c) shows the corresponding transmittance images, mean values, and line profiles obtained from 3D-PLI simulations with normally incident light and λ = 550 nm for the straight single fibers shown in Fig. 13(b).The mean values and line profiles for L ≥ 1 look similar.For better comparison, Fig. 13(d) shows the absolute relative differences (ARDM and RMAD) between the transmittance images with L = {0, 1, 2, 3, 4, 5} and the transmittance image with realistic myelin sheath.The relative differences decrease with increasing number of myelin layers L and with decreasing mesh size ∆.A fiber with two or more myelin layers and a mesh size ∆ = t g /3 yields similar transmittance values as the fiber with realistic myelin sheath.With increasing mesh size, the relative differences increase.For a fiber with double myelin layers and a mesh size ∆ = 12.5 nm (25 nm), the differences are: ARDM ≈ 1.2% (2.3%) and RMAD ≈ 1.6% (2.8%).For a mesh size of 25 nm, the differences are still smaller than for a fiber without or with a single myelin layer.Thus, a fiber with double myelin layers and a mesh size of 25 nm is a good compromise between accuracy and computing time and was used for all simulation studies in Sec.III.In interesting cases, the simulations were repeated for a reduced mesh size (∆ = 12.5 nm), see black crosses in Fig. 14(b).The simulations were performed for normally incident light with 550 nm wavelength and simulation parameters specified in Tab.II.The profiles with non-italic labels belong to the displayed transmittance images.(d) Relative differences between the transmittance images with different numbers L of myelin layers and different mesh sizes ∆ (relative to the glycerin layer thickness t g ) and the transmittance image with realistic myelin sheath.The values for ARDM (blue) and RMAD (orange) were computed using Eq.(G1) and Eq.(G2), the values surrounded in red belong to fibers with double myelin layers (L = 2) and 25 nm mesh size, which were used for the simulation studies in Sec.III.
were performed for a Yee mesh size of 25 nm.For some inclination angles, the simulations for NA = 0.15 were repeated for a smaller mesh size (∆ = 12.5 nm), see black crosses.The transmittance curves for different wavelengths and for diffusive light (obtained from simulation runs with different angles of incidence) look all very similar.The maximum difference between the normalized transmittance values is less than 0.03.In addition, the simulations with smaller mesh size (black crosses) yield similar results as the simulations with larger mesh size (curves), the maximum difference between the normalized transmittance values is only about 0.005.Thus, the transmittance curves for the bundle of densely grown fibers are not sensitive to small changes in wavelength, angle of incidence, or mesh size, which demonstrates the validity of the simulation results.).The corpus callosum (cc) contains mainly out-of-plane nerve fibers, while the cingulum (cg) contains mainly in-plane nerve fibers (cf.Fig. 1).The mean transmitted light intensity of the selected region in the corpus callosum is about 53% (a) and 32% (b) less than the mean transmitted light intensity of the selected region in the cingulum, respectively.This demonstrates that the same effects (less transmitted light intensity in regions with out-of-plane nerve fibers than in regions with in-plane nerve fibers) can also be observed in conventional transmission microscopy which uses unpolarized light.The electromagnetic field components behind the sample (represented by a set of vectors {A k , B k , k}, green box) were computed by TDME3D and shifted to the middle of the corresponding Yee cell (with side length ∆).To study how much light is scattered under a certain angle (wave vector k), the scattering pattern I k was computed, i. e., the intensity per wave vector normalized by the ingoing light intensity I 0 per image pixel (orange box).The spherical microlenses of the camera detector were modeled by applying a moving average over the area of the microlens with radius 0.665 µm (J 1 is the Bessel function of the first kind of order one).The numerical aperture of the imaging system was modeled by considering only waves with directions of propagation k that fulfill θ k < arcsin(NA) ≈ 8.63 • .Applying a 2D inverse discrete Fast Fourier Transform (FFT), the transmitted light intensities were computed and normalized by the ingoing light intensity per pixel (red box).

FIG. 1 .
FIG. 1. Transmittance and retardation images of coronal and sagittal brain sections for a rat (a) and a vervet monkey (b).The coronal (sagittal) section planes are indicated by blue (red) lines in the respective other brain section for reference, selected anatomical structures are labeled (see legend).The upper two rows of each panel show the normalized transmittance images I T,N , the third row shows the retardation images | sin δ|, both obtained from 3D-PLI measurements with 1.33 µm pixel size.The images in the first row show the whole brain section, the images in the second and third row show an enlarged view.The selected regions in yellow (green) belong to steep (flat) nerve fibers, which appear dark (bright) in the transmittance and retardation images.

FIG. 2 .
FIG. 2. Bright-field transmission microscopy (see Appendix A 4) vs. 3D-PLI measurement of a human brain section (right occipital lobe).(a) Photograph of the brain block surface before sectioning (blockface image).The enlarged region shows the sagittal stratum, a white matter structure that runs mostly perpendicular to the section plane (GM = gray matter).(b) Bright-field transmission microscopy image of the same region with 0.91 µm pixel size.(c) Normalized transmittance image of the same region obtained from a 3D-PLI measurement with 1.33 µm pixel size.(d) Corresponding retardation image.Regions with steep out-of-plane fibers (sagittal stratum) show lower transmittance (and retardation) values than the neighboring regions with non-steep fibers; the transmitted light intensity images obtained from bright-field transmission microscopy (b) and 3D-PLI (c) look similar.

FIG. 3 .
FIG. 3. 3D-reconstructed normalized transmittance images (I T,N ) of the right hemisphere of a vervet monkey brain (234 consecutive sections from the occipital lobe) obtained from 3D-PLI measurements with 1.33 µm pixel size.The brain was cut along the coronal plane (xy-plane), the resulting brain sections were registered onto each other in the z-direction.(a)-(c) Cross-sections of the 3D-volume shown along the coronal (xy), sagittal (xz), and horizontal (yz) plane.The colored lines indicate the position of the displayed xy-, xz-, and yz-planes.(d) Detail of the 3D-volume.The white arrows point to the sagittal stratum -a white matter structure that runs mostly perpendicular to the image plane (along the z-direction) and which appears much darker in the transmittance images than the surrounding tissue.

FIG. 4 .
FIG. 4. Transmittance contrast of nerve fiber bundles in mutually orthogonal anatomical planes.(a) Normalized transmittance images (I T,N ) of a coronal, sagittal, and horizontal rat brain section obtained from 3D-PLI measurements with 1.33 µm pixel size.The colored lines indicate the approximate position of the section planes.The enlarged views show the region of the caudate putamen and the maximum and minimum angles under which the nerve fiber bundles are oriented with respect to the section planes.(b) Corresponding retardation images (| sin δ|) of the enlarged views.The image contrast was used to select regions with fibers and with surrounding tissue in the caudate putamen (yellow lines).(c) Histograms of the transmittance values (I T,N) for the selected regions with nerve fibers (pink) and with surrounding tissue (cyan) in the caudate putamen.For the coronal brain section, the retardation image cannot be used to separate the fibers from the surrounding tissue because the fibers are oriented almost perpendicular to the image plane which leads to a small retardation signal and a small image contrast.Therefore, the histogram was computed over the whole selected region and the peak with lower (larger) transmittance was assumed to belong to fibers (surrounding tissue).The contrast values were computed from the respective peak values (numbers in pink and cyan) via: C = (max -min)/(max + min).The contrast for steep nerve fibers (coronal brain section) is much larger than for flat nerve fibers (sagittal and horizontal brain section).

FIG. 5 .
FIG. 5. Transmittance of inclined fiber bundles.(a) 3D view and cross-sections through mid-planes for an artificial bundle of densely grown fibers shown exemplary for an inclination angle α = 45 • .(b) Light-scattering patterns obtained from 3D-PLI simulations for the bundle of densely grown fibers with α = 0 • and 70 • .To compute the mean normalized transmittance values I T,N for the numerical aperture of the imaging system (NA = 0.15), only wave vector angles θ k ≤ 8.6 • were considered (see red circles).The white circles represent steps of ∆θ k = 10 • .(c) Simulated transmittance curves (mean transmittance I T,N vs. inclination α) for the bundle of densely grown fibers (i) and for a bundle with broad fiber orientation distribution (ii) for NA = 1 (orange curves) and NA = 0.15 (blue curves).The transmittance curves were normalized by the mean transmittance values of the horizontal bundles, respectively.The simulations were performed with the parameters specified in Appendix E, using normally incident light with 550 nm wavelength.(d) Mean normalized transmittance values I T,N plotted against the nerve fiber inclination angles α determined respectively from 3D-PLI and TPFM measurements of nerve fiber bundles in a mouse brain section (see Fig. S16 in the Supplemental Material).The values in blue belong to regions with similar (maximum) fiber density, the values in orange belong to regions with variable fiber density in which the transmittance might be overestimated.The error bars indicate the standard error of the mean for the evaluated transmittance values.Both experimental and simulated data show that the transmittance decreases with increasing fiber inclination.

FIG. 6 .
FIG.6.Normalized transmittance images I T,N of the sagittal rat brain section in Fig.1(a) obtained from a 3D-PLI measurement one day after tissue embedding and 16 months later.With increasing time after the tissue embedding, the brain section becomes more transparent.

Figure 7 ( 3 /FIG. 7 .
FIG. 7. Simulated transmittance of crossing fibers.(a)-(c) 3D views and light-scattering patterns for (a) horizontal crossing fibers (separate bundles in the upper row, interwoven bundles in the lower row) with crossing angle χ, (b) three mutually orthogonal, interwoven fiber bundles, and (c) a vertical fiber bundle with broad fiber orientation distribution (cf.Fig.5(c)(ii)).The mean normalized transmittance values I T,N were computed from a simulated 3D-PLI measurement (with numerical aperture NA = 0.15).The simulations were performed with the parameters specified in Appendix E, using normally incident light with 550 nm wavelength.(d) Mean transmittance values for different crossing angles χ shown for in-plane crossing (solid curves) and out-ofplane fiber configurations (densely dotted curves).For better comparison, the values were divided by the mean transmittance value of the corresponding horizontal fiber bundle for χ = 0 • .Apart from the fiber configurations shown in this figure, the mean transmittance values are also displayed for the bundle of densely grown fibers (see Fig.5(c)(i)) for α = 70 • and 90 • .While the scattering patterns show the crossing angles of in-plane crossing fibers (see (a)), the mean transmittance is mostly independent of the crossing angle and larger than the mean transmittance of out-of-plane fibers (see (d)).
9): 1.I T,N I ref : regions with notably lower transmittance values are expected to contain steep (out-ofplane) fibers (see yellow arrows and regions surrounded by a yellow line), 2. I T,N ∼ I ref : regions with similar transmittance values are expected to contain flat (in-plane) crossing fibers (see cyan arrows), 3. I T,N I ref : regions with notably larger transmittance values are expected to have a lower fiber density (see magenta arrows).

FIG. 8 .
FIG. 8. Crossing nerve fibers in the optic chiasm of a hooded seal: (a) brain tissue before sectioning, (b) unnormalized transmittance and retardation images of the middle brain section obtained from 3D-PLI measurements with 1.33 µm pixel size, (c) schematic drawing of the optic chiasm consisting of optic tracts (o.t.) and optic nerves (o.n.),(d) normalized histograms of the transmittance image (I T ) and retardation image (| sin δ|) for a region with mostly parallel fibers (blue) and a region with nearly 90 • -crossing fibers (orange).Unlike the retardation, the transmittance does not depend on the crossing angles of the nerve fibers, only on the tissue density.More information about the sample can be found in Dohmen et al. [52] (Figs.(a) and (c) were adapted from Figs. 1B and 5B0 in [52] Copyright (2015), with permission from Elsevier).

FIG. 9 .
FIG. 9. Combined analysis of transmittance and retardation images allowing to distinguish between brain regions with in-plane crossing and out-of-plane nerve fibers.The figure shows the normalized transmittance image (I T,N ) and the retardation image (| sin δ|) of a coronal section through the right hemisphere (occipital lobe) of a vervet brain (cf.Fig. 3(a)).The transmittance in the region with maximum retardation (orange ellipse) is used as a threshold value (I ref ≡ I T,N (| sin δ| max )).Regions with small retardation values and notably lower transmittance values (I T,N I ref , yellow arrows and regions surrounded by a yellow line) are expected to contain steep (out-of-plane) fibers.Regions with small retardation values and similar transmittance values (I T,N ∼ I ref , cyan arrows) are expected to contain flat (in-plane) crossing fibers.Regions with small retardation values and larger transmittance values (I T,N I ref , magenta arrows) belong to regions with low fiber density, i. e., regions with a large amount of unmyelinated axons or surrounding tissue.

FIG. 10 .
FIG. 10.Generation of crossing fibers.(a)-(b) Separate and interwoven fiber bundles with crossing angle χ: The upper figures show the generated bundles before cropping.The lower figures show the bundles after being cropped to a volume of 30 × 30 × 30 µm 3 .The white dotted line indicates the border between the upper and the lower bundle of the separate crossing fibers.(c) Three mutually orthogonal, interwoven fiber bundles cropped to a volume of 30 × 30 × 30 µm 3 .The white dotted lines indicate the main directions of the two horizontal fiber bundles in the xy-plane, the third fiber bundle is oriented in the z-direction.All fiber configurations were generated from 700 fibers with diameters between 1.0 µm and 1.6 µm.

FIG. 11 .
FIG. 11.Modeling of nerve fibers.(a) Schematic drawing of a nerve fiber (myelinated axon).(b) Cross-section through the nerve fiber showing the inner axon and the surrounding myelin sheath (formed by a type of glial cell which spirally wraps around the axon).(c) Schematic representation of the myelin structure consisting of several lipid bilayers (5 nm thick cell membranes)with an intracellular/cytoplasmic and an extracellular space of about 3 nm.(d) Each cell layer (two lipid bilayers with separating cytoplasm) was considered as one "myelin layer" with an effective refractive index n m = 1.47 (blue), the extracellular space was considered to be filled with glycerin solution ("glycerin layer") with a refractive index n g = 1.37 (yellow).The myelin and glycerin layers were assumed to contribute 3/4 and 1/4 to the overall myelin sheath thickness t sheath , respectively.(e) Nerve fibers were modeled with double myelin layers with thickness t m = (3/7) t sheath and a single separating glycerin layer with thickness t g = (1/7) t sheath .The myelin sheath thickness contributes approximately one third to the overall fiber radius (t sheath = 0.35 r).The inner axon was modeled with a radius r ax = 0.65 r and a refractive index n ax = 1.35.

•
FIG.12.Modeling of the 3D-PLI measurement.The figure and table on the left-hand side show the optical components of the polarimeter (the order of the polarizing filters is different than in the measurement, but the setup is mathematically equivalent): light source (green), polarizer/retarder (dark gray), sample (light gray), objective lens/detector/camera (blue).The table and figure on the right-hand side show how the optical elements were modeled by FDTD simulations: The incoherent and diffusive light source (LED) with peak wavelength λ and full-width at half-maximum (FWHM) was modeled by performing several simulation runs with plane waves that have different wavelengths (λ) and angles of incidence (ϕ, θ).The modeled light source emits coherent light that is circularly polarized.The tissue sample was represented by an artificial fiber architecture, the rotating analyzer by a rotated Jones matrix (with rotation matrix R(ρ)).The numerical aperture (NA) of the imaging system was modeled by considering only wave vector angles θ k < arcsin(NA).The spherical microlenses of the camera detector were modeled by performing a moving average over the area of the microlens with radius r = 1.33 µm / 2.

FIG. 13 .
FIG.13.Error estimation for different numbers of myelin layers.(a) Dimensions of the simulation volume (xy/yz-plane) used to simulate a straight single fiber with different numbers of myelin layers.(b) Cross-section through fibers with different numbers L of myelin layers and different Yee mesh sizes ∆.All fibers were modeled with a diameter of 1 µm, consisting of an inner axon (green) with a diameter of 0.65 µm and a surrounding myelin sheath with a thickness of 0.175 µm.The myelin sheath is composed of alternating layers of myelin (blue) and glycerin (yellow), the myelin layers are three times thicker than the glycerin layers.The realistic model of the myelin sheath contains 22 layers of 5 nm thick cell membranes (blue), interrupted by 3 nm thick alternating layers of cytoplasm (green) and surrounding glycerin solution (yellow), yielding a myelin sheath composed of 43 thin layers.The refractive indices are 1.35 for the axon/cytoplasm (green), 1.37 for the glycerin solution (yellow), and 1.47 for the myelin layers (blue).A motivation of the myelin sheath model and the corresponding refractive indices is shown in Fig.11.(c) Normalized transmittance images, corresponding mean values I T,N , and transmittance profiles on the right (middle image pixels evaluated along the y-axis, see white dashed lines) obtained from 3D-PLI simulations with different numbers L of myelin layers and Yee mesh sizes defined in (b).The simulations were performed for normally incident light with 550 nm wavelength and simulation parameters specified in Tab.II.The profiles with non-italic labels belong to the displayed transmittance images.(d) Relative differences between the transmittance images with different numbers L of myelin layers and different mesh sizes ∆ (relative to the glycerin layer thickness t g ) and the transmittance image with realistic myelin sheath.The values for ARDM (blue) and RMAD (orange) were computed using Eq.(G1) and Eq.(G2), the values surrounded in red belong to fibers with double myelin layers (L = 2) and 25 nm mesh size, which were used for the simulation studies in Sec.III.

FIG. 14 .
FIG.14.Simulated transmitted light intensity for different simulation parameters.Bundle of densely grown fibers (cf.Fig.5(a)) simulated for different inclination angles α: (a) Mean transmitted light intensity for light polarized along the x-axis (I x ) or along the y-axis (I y ) plotted against α.The simulations were performed for a numerical aperture NA = 0.15 using a normally incident plane wave with 550 nm wavelength and simulation parameters specified in Tab.II.(b) Transmittance curves (mean transmittance value I T,N plotted against α) obtained from 3D-PLI simulations for normally incident light with different wavelengths λ = {545, 550, 555} nm and a Yee mesh size of ∆ = 25 nm.The simulations for λ = 550 nm were also performed for diffusive light (with angles of incidence {θ = 0 • }, {θ = 3 • ; ϕ = 4 × 90 • }).The curves were normalized by the mean transmittance value of the horizontal bundle, respectively.The solid curves were computed for the numerical aperture of the imaging system (NA = 0.15), the dashed curves were computed for NA = 1.The black crosses belong to simulations with NA = 0.15, λ = 550 nm, and ∆ = 12.5 nm.

FIG. 15 .
FIG. 15.Supplementary Figure | 3D-PLI measurement vs. bright-field transmission microscopy of a sagittal vervet brain section (neighboring section to the one shown in Fig. 1(b), enlarged region): (a) Normalized transmittance image obtained from a 3D-PLI measurement (see Appendix A 2).(b) Bright-field transmission microscopy image obtained from ZEISS Axio Imager Vario with unpolarized light (see Appendix A 4).The corpus callosum (cc) contains mainly out-of-plane nerve fibers, while the cingulum (cg) contains mainly in-plane nerve fibers (cf.Fig.1).The mean transmitted light intensity of the selected region in the corpus callosum is about 53% (a) and 32% (b) less than the mean transmitted light intensity of the selected region in the cingulum, respectively.This demonstrates that the same effects (less transmitted light intensity in regions with out-of-plane nerve fibers than in regions with in-plane nerve fibers) can also be observed in conventional transmission microscopy which uses unpolarized light.

FIG. 16 .FIG. 17 .
FIG.16.Supplementary Figure | Transmittance vs. inclination of nerve fiber bundles.A coronal, 60 µm thin section from the left hemisphere of a mouse brain was measured with 3D-PLI (voxel size: 1.33 × 1.33 × 60 µm 3 ) and with TPFM (voxel size: 0.24 × 0.24 × 1 µm 3 ) for a region of 3 × 3 tiles in the caudate putamen: (a) normalized transmittance image of the coronal mouse brain section, (b) normalized transmittance image of the selected region in the caudate putamen, (c) single tile of a contrastenhanced TPFM image stack demonstrating the determination of the fiber inclination angle α, (d) middle slice of the TPFM image stack for the selected region in the caudate putamen, (e) mean normalized transmittance values I T,N plotted against the determined fiber inclination angles α.The inclination angles of 40 fiber bundles were determined from the center point positions of the cross-sections in the first and the last slice of the TPFM image stack (see colored shapes in Fig.(d)), the transmittance values were evaluated in the middle of the corresponding fiber bundles (see colored shapes in Fig.(b)).The blue (orange) data points in the scatter plot belong to regions that are (are not) completely filled with nerve fibers, cf.blue (orange) arrows in Figs.(b),(d).The error bars indicate the standard error of the mean for the evaluated transmittance values for each region.Although the values are broadly distributed, the scatter plot shows a clear tendency towards a decrease in transmittance with increasing fiber inclination, also in regions with maximum fiber density.

FIG. 18 ., e 2 ,FIG. 19 .
FIG. 18. Supplementary Figure | Simulated transmittance curves.Mean normalized transmittance values I T,N plotted against the inclination angles α for different fiber bundles and different myelin refractive indices n m : (a) bundle of densely grown fibers with n m = 1.47 (corresponds to given literature values of lipids), (b) bundle of densely grown fibers with n m = 1.39 (the reduced myelin refractive index models tissue with long embedding time), (c) bundle with broad fiber orientation distribution and n m = 1.47.The orange curves were computed from 3D-PLI simulations for the numerical aperture of the imaging system (NA = 0.15) and the blue curves for NA = 1.To enable a better comparison between horizontal and vertical fibers, the lower figures show the transmittance curves normalized by the mean transmittance of the horizontal fiber bundle, respectively.The 3D-PLI simulations were performed for a normally incident plane wave with 550 nm wavelength and simulation parameters specified in Tab.II.A reduced myelin refractive index (long embedding time) leads to larger transmittance values and a smaller decrease for steep fibers for NA = 0.15.While the transmittance curves for the bundle of densely grown fibers show a minimum for steep fibers (α = 70-80 • ), the mean transmittance values for the bundle with broad fiber orientation distribution decrease monotonically with increasing fiber inclination.

TABLE I .
Mean transmittance values (I

TABLE II .
Parameters for the simulation studies in Sec.III: expenses of one simulation run (computation of one fiber configuration, one wavelength, and one angle of incidence on JUQUEEN), dimensions of the simulation volume, and fiber properties (radius r, thickness t, refractive index n).