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.

Popular Summary

To understand the structure and function of the brain, we need to study the highly complex, 3D neuronal connections. To enable a correct reconstruction of the nerve fiber pathways, it is crucial to know the detailed substructure of the tissue, especially in regions with crossing fibers. Neuroimaging techniques such as 3D polarized light imaging (3D-PLI) reveal the fiber pathways with micrometer resolution, but fiber crossings still pose a major problem. Here, we show how light scattering can be used to detect fiber crossings in 3D-PLI images and to reveal their substructure.

In various experimental and simulation studies, we find that the transmitted light intensity strongly depends on the angle between the fibers and the direction of light propagation. It can be used not only to reveal 3D information but also to identify crossing fibers. Furthermore, we demonstrate that optical scattering reveals the substructure of brain tissue such as the crossing angles of nerve fibers. To explain our experimental observations, we develop a simulation framework for polarization microscopy that allows us to study light scattering on fibrous tissue models using high-performance computing.

Our results provide an improved reconstruction of nerve fiber pathways in the brain at neuronal scales and allow for a better understanding of the brain’s microstructure. We anticipate that scattering microscopy, which reveals difficult-to-access information about brain tissue at micrometer resolution, will become an integral part of future neuroimaging techniques. Also, the simulation framework and results can easily be generalized to other microscopy techniques and fibrous tissue samples, enabling applications beyond neuroscience.

Figure 1

Fiber orientation map of a coronal vervet monkey brain section, obtained from a 3D-PLI measurement with $1.33\mu \mathrm{m}$ pixel size. The enlarged area highlights three brain regions that all yield low birefringence signals: (i) gray matter with low fiber density, (ii) crossing nerve fibers (corona radiata), and (iii) steep out-of-plane fibers (fornix).

Figure 2

Normalized transmittance images $I_{T,N}$ of a coronal (a) and a sagittal (b) vervet brain section, obtained from a 3D-PLI measurement with $1.33\mu \mathrm{m}$ pixel size (the fiber orientation map in Fig. 1 is obtained from the same coronal brain section). For reference, the coronal (sagittal) section plane is indicated by a vertical blue (red) line in the respective other brain section. The enlarged areas on the right show anatomical brain regions with predominantly in-plane (out-of-plane) nerve fibers in green (yellow). Within each section, in-plane nerve fibers have larger transmittance values than out-of-plane fibers. $\mathrm{cg}=\mathrm{cingulum}$, $\mathrm{cc}=\mathrm{corpus}$ callosum, and $\mathrm{f}=\mathrm{fornix}$.

Figure 3

Transmittance of inclined fiber bundles. (a) Mean normalized transmittance values $\overline{I_{T,N}}$ plotted against the nerve fiber inclination angles $\alpha$ determined, respectively, from 3D-PLI and TPFM measurements of nerve fiber bundles in a mouse brain section (see Fig. S2 in Supplemental Material [55]). The values in blue belong to regions with similar (maximum) fiber density, and 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. (b) Simulated transmittance curves (mean transmittance $\overline{I_{T,N}}$ vs inclination $\alpha$) for a bundle of densely grown fibers (i) and a bundle with broad fiber orientation distribution (ii) for $\mathrm{NA}=1$ (purple curves) and $\mathrm{NA}=0.15$ (green curves). The transmittance curves are normalized by the mean transmittance values of the horizontal bundles, respectively. The simulations are performed with the parameters specified in Appendix pp6, using normally incident light with 550 nm wavelength. Both experimental and simulated data show that the transmittance decreases with an increasing fiber inclination (for $\mathrm{NA}=0.15$).

Figure 4

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\mu \mathrm{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 \delta |$) 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. [56] [(a) and (c) are adapted from Figs. 1 and 5B0 in [56] copyright 2015, with permission from Elsevier].

Figure 5

Simulated transmittance of crossing fibers. (a),(b) Separate and interwoven fiber bundles with crossing angle $\chi$: The upper figures show the generated bundles before cropping. The lower figures show the bundles after being cropped to a volume of $30\times 30\times 30{\mu \mathrm{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\times 30\times 30{\mu \mathrm{m}}^3$. The white dotted lines indicate the main directions of the two horizontal fiber bundles in the $xy$ plane, and the third fiber bundle is oriented in the $z$ direction. All fiber configurations are generated from 700 fibers with diameters between 1.0 and $1.6\mu \mathrm{m}$. (d) Mean transmittance values for different crossing angles $\chi$ shown for in-plane crossing (solid curves) and out-of-plane fiber configurations (densely dotted curves). For better comparison, the values are divided by the mean transmittance value of the corresponding horizontal fiber bundle for $\chi =0°$. The mean normalized transmittance values $\overline{I_{T,N}}$ are computed from a simulated 3D-PLI measurement (with numerical aperture $\mathrm{NA}=0.15$). The simulations are performed with the parameters specified in Appendix pp6, using normally incident light with 550 nm wavelength. 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. 3] for $\alpha =70°$ and 90° and for a vertical fiber bundle with broad fiber orientation distribution [see Fig. 3]. The mean transmittance is mostly independent of the crossing angle and larger than the mean transmittance of out-of-plane fibers.

Figure 6

Combined analysis of transmittance and retardation images allowing the classification of brain regions with low birefringence signals in 3D-PLI measurements. The figure shows the normalized transmittance image ($I_{T,N}$) and the retardation image ($|\sin \delta |$) of a coronal section through the right hemisphere (occipital lobe) of a vervet monkey brain [cf. Fig. 13 in Appendix pp2-s2] obtained from a 3D-PLI measurement with $1.33\mu \mathrm{m}$ pixel size. The transmittance in the region with maximum retardation (orange ellipse) is used as a reference value: Regions with small retardation values and notably lower transmittance values (regions surrounded by a yellow line) are expected to contain steep (out-of-plane) fibers. Regions with small retardation values and similar transmittance values (cyan) are expected to contain flat (in-plane) crossing fibers. Regions with small retardation values and larger transmittance values (purple) belong to regions with low fiber density, i.e., regions with a large amount of unmyelinated axons or surrounding tissue.

Figure 7

Simulated scattering patterns for different artificial nerve fiber constellations: (a) densely grown fiber bundle [cf. Fig. 3] with different inclination angles $\alpha$ and (b) in-plane crossing fibers [separate and interwoven bundles; cf. Figs. 5 and 5] with different crossing angles $\chi$. The scattering patterns show the underlying substructure, e.g., the crossing angle of the fibers (indicated by the black lines around the patterns).

Figure 8

Scattering measurement of two optic tracts, extracted from a 30-$\mu \mathrm{m}$-thin section of a human optic chiasm, placed on top of each other with a crossing angle of approximately 80°. The measurement is performed as described in Appendix pp1-s5 with a pixel size in object space of about $6.5\mu \mathrm{m}$ and a mask with a rectangular hole which illuminates the sample under a fixed polar angle of $49.1°$. During the measurement, the mask is rotated by angles of $\{0°,15°,\dots ,345°\}$ around the center of the sample (starting on top and rotating clockwise; see the compass rose on top). The upper image shows the transmitted light intensity recorded by the camera averaged over all rotation angles of the mask. For evaluation, a region of $10\times 10\mathrm{pixels}$ is selected in each of the two fiber bundles [(1) and (2)] and in the crossing region (3). The selected regions are indicated by a yellow square in the upper image, and the in-plane orientations of the nerve fibers (derived from visible tissue structures) are indicated by green and purple lines. The graphs show the average transmitted light intensity $I$ in the three evaluated regions plotted against the rotation angle of the mask (black curves). The red curves show the transmitted light intensity of the middle pixel ($1\times 1\mathrm{pixels}$) in the selected regions. The dashed vertical lines indicate the orientations of the fibers as marked in the upper image. The dash-dotted lines in (3) indicate the position of the measured peaks.

Figure 9

Scattering measurement of a coronal vervet brain section in comparison to simulated scattering patterns. (a) The scattering measurement is performed as described in Appendix pp1-s5 with a pixel size in object space of $6.3\mu \mathrm{m}$ and a mask with a circular hole with a fixed polar angle of illumination of $47.6°$ and rotation angles of $\{0°,22.5°,\dots ,337.5°\}$. Three regions with parallel in-plane fibers (corpus callosum), parallel out-of-plane fibers (fornix), and in-plane crossing fibers (corona radiata) are selected for evaluation. The images in the top row show the fiber orientation maps obtained from a 3D-PLI measurement with $1.33\mu \mathrm{m}$ pixel size (the fiber orientation map of the whole brain section is shown in Fig. 1). The selected regions ($10\times 10\mathrm{pixels}$) are marked by a small square, and the in-plane orientations of the nerve fibers (derived from visible structures and/or anatomical knowledge) are marked by purple and green lines. The second row shows the average transmitted light intensity $I$ in the evaluated regions plotted against the rotation angle of the mask (black curves). The red curves show the transmitted light intensity of the middle pixel ($1\times 1\mathrm{pixels}$) in the selected regions. The dashed vertical lines indicate the orientations of the fibers as marked in the fiber orientation maps. (b) The bottom row shows artificial fiber bundles and the simulated scattering patterns (cf. Fig. 7): parallel in-plane fibers, parallel out-of-plane fibers with 70° inclination, and interwoven crossing fibers with a 90° crossing angle. The green and purple lines around the scattering patterns indicate the main orientations of the artificial fibers. The white dashed circle indicates the angle under which the sample is illuminated in the measurement ($47.6°$). The graphs above show the intensity of the corresponding scattering pattern, evaluated along the dashed circle in the clockwise direction starting at the top (see the white arrow in the scattering pattern). To account for the finite size of the hole of illumination, a Gaussian blur with a diameter of 8° in angular space is applied to the scattering patterns before evaluation.

Figure 10

Schematic setup of the scattering measurement. (a) Side view of the setup, consisting of an LED panel with a mask, a specimen stage with a sample, and a camera. The numbers without brackets belong to the scattering measurement of the vervet brain section (Sec. 4c), and the numbers in brackets to the measurement of the human optic chiasm (Sec. 4b). (b) Top view of the mask used for the scattering measurement in Sec. 4c (circular hole with 3.5 cm diameter and 11.5 cm distance from the center, yielding an illumination angle of $47.6°$ with respect to the center of the sample). During the measurement, the mask is rotated in steps of $22.5°$ around the center in the clockwise direction. [For the scattering measurement in Sec. 4b, a mask with a rectangular hole ($2.4\times 4{\mathrm{cm}}^2$) and steps of 15° are used.].

Figure 11

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, and selected anatomical structures are labeled (see the legend). The upper two rows of each panel show the normalized transmittance images $I_{T,N}$, and the third row shows the retardation images $|\sin \delta |$, both obtained from 3D-PLI measurements with $1.33\mu \mathrm{m}$ pixel size. The images in the first row show the whole brain section, and the images in the second and third rows 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.

Figure 12

Bright-field transmission microscopy 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 ($\mathrm{GM}=\text{gray}$ matter). (b) Bright-field transmission microscopy image of the same region with $0.91\mu \mathrm{m}$ pixel size. (c) Normalized transmittance image of the same region obtained from a 3D-PLI measurement with $1.33\mu \mathrm{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 nonsteep fibers; the transmitted light intensity images obtained from bright-field transmission microscopy (b) and 3D-PLI (c) look similar.

Figure 13

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\mu \mathrm{m}$ pixel size. The brain is cut along the coronal plane ($xy$ plane), and the resulting brain sections are 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$) planes. 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.

Figure 14

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\mu \mathrm{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 \delta |$) of the enlarged views. The image contrast is 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 is computed over the whole selected region, and the peak with lower (larger) transmittance is assumed to belong to fibers (surrounding tissue). The contrast values are computed from the respective peak values (numbers in pink and cyan) via $\mathcal{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).

Figure 15

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 or cytoplasmic and an extracellular space of about 3 nm. (d) Each cell layer (two lipid bilayers with separating cytoplasm) is considered as one “myelin layer” with an effective refractive index $n_m=1.47$ (blue), the extracellular space is considered to be filled with glycerin solution (“glycerin layer”) with a refractive index $n_g=1.37$ (yellow). The myelin and glycerin layers are assumed to contribute $3/4$ and $1/4$ to the overall myelin sheath thickness $t_{\text{sheath}}$, respectively. (e) Nerve fibers are modeled with double myelin layers with thickness $t_m=(3/7)t_{\text{sheath}}$ and a single separating glycerin layer with thickness $t_g=(1/7)t_{\text{sheath}}$. The myelin sheath thickness contributes approximately one-third to the overall fiber radius ($t_{\text{sheath}}=0.35r$). The inner axon is modeled with a radius $r_{\mathrm{ax}}=0.65r$ and a refractive index $n_{\mathrm{ax}}=1.35$.

Figure 16

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 and retarder (dark gray), sample (light gray), and objective lens, detector, and camera (blue). The table and figure on the right-hand side show how the optical elements are modeled by FDTD simulations: The incoherent and diffusive light source (LED) with peak wavelength $\widehat{\lambda }$ and full width at half maximum (FWHM) is modeled by performing several simulation runs with plane waves that have different wavelengths ($\lambda$) and angles of incidence ($\phi$, $\theta$). The modeled light source emits coherent light that is circularly polarized. The tissue sample is represented by an artificial fiber architecture, the rotating analyzer by a rotated Jones matrix [with rotation matrix $R(\rho )$]. The numerical aperture (NA) of the imaging system is modeled by considering only wave vector angles ${\theta }_k<\arcsin (\mathrm{NA})$. The spherical microlenses of the camera detector are modeled by performing a moving average over the area of the microlens with radius $r=1.33\mu \mathrm{m}/2$.

Figure 17

Error estimation for different numbers of myelin layers. (a) Dimensions of the simulation volume ($xy$ or $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 $\mathit{\Delta }$. All fibers are modeled with a diameter of $1\mu \mathrm{m}$, consisting of an inner axon (green) with a diameter of $0.65\mu \mathrm{m}$ and a surrounding myelin sheath with a thickness of $0.175\mu \mathrm{m}$. The myelin sheath is composed of alternating layers of myelin (blue) and glycerin (yellow), and the myelin layers are 3 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 or 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. 15. (c) Normalized transmittance images, corresponding mean values $\overline{I_{T,N}}$, and transmittance profiles on the right (middle image pixels evaluated along the $y$ axis; see the white dashed lines) obtained from 3D-PLI simulations with different numbers $L$ of myelin layers and Yee mesh sizes defined in (b). The simulations are performed for normally incident light with 550 nm wavelength and simulation parameters specified in Table 2. The profiles with nonitalic 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 $\mathit{\Delta }$ (relative to the glycerin layer thickness $t_g$) and the transmittance image with a realistic myelin sheath. The values for ARDM (blue) and RMAD (orange) are computed using Eqs. (h1) and (h2); the values surrounded in red belong to fibers with double myelin layers ($L=2$) and 25 nm mesh size, which are used for the simulation studies in Secs. 3 and 4.

Figure 18

Simulated transmitted light intensity for different simulation parameters. Bundle of densely grown fibers [cf. Fig. 3] simulated for different inclination angles $\alpha$: (a) Mean transmitted light intensity for light polarized along the $x$ axis ($\overline{I_x}$) or along the $y$ axis ($\overline{I_y}$) plotted against $\alpha$. The simulations are performed for a numerical aperture $\mathrm{NA}=0.15$ using a normally incident plane wave with 550 nm wavelength and simulation parameters specified in Table 2. (b) Transmittance curves (mean transmittance value $\overline{I_{T,N}}$ plotted against $\alpha$) obtained from 3D-PLI simulations for normally incident light with different wavelengths $\lambda =\{545,550,555\}\mathrm{nm}$ and a Yee mesh size of $\mathit{\Delta }=25\mathrm{nm}$. The simulations for $\lambda =550\mathrm{nm}$ are also performed for diffusive light (with angles of incidence $\{\theta =0°\}$, $\{\theta =3°$; $\phi =4\times 90°\}$). The curves are normalized by the mean transmittance value of the horizontal bundle, respectively. The solid curves are computed for the numerical aperture of the imaging system ($\mathrm{NA}=0.15$), and the dashed curves are computed for $\mathrm{NA}=1$. The black crosses belong to simulations with $\mathrm{NA}=0.15$, $\lambda =550\mathrm{nm}$, and $\mathit{\Delta }=12.5\mathrm{nm}$.