Automated computation of nerve fibre inclinations from 3D polarised light imaging measurements of brain tissue

The method 3D polarised light imaging (3D-PLI) measures the birefringence of histological brain sections to determine the spatial course of nerve fibres (myelinated axons). While the in-plane fibre directions can be determined with high accuracy, the computation of the out-of-plane fibre inclinations is more challenging because they are derived from the amplitude of the birefringence signals, which depends e.g. on the amount of nerve fibres. One possibility to improve the accuracy is to consider the average transmitted light intensity (transmittance weighting). The current procedure requires effortful manual adjustment of parameters and anatomical knowledge. Here, we introduce an automated, optimised computation of the fibre inclinations, allowing for a much faster, reproducible determination of fibre orientations in 3D-PLI. Depending on the degree of myelination, the algorithm uses different models (transmittance-weighted, unweighted, or a linear combination), allowing to account for regionally specific behaviour. As the algorithm is parallelised and GPU optimised, it can be applied to large data sets. Moreover, it only uses images from standard 3D-PLI measurements without tilting, and can therefore be applied to existing data sets from previous measurements. The functionality is demonstrated on unstained coronal and sagittal histological sections of vervet monkey and rat brains.

To better understand the function of the brain and to treat neurological diseases, a detailed reconstruction of the intricate and densely grown nerve fibre network is needed. Diffusion magnetic resonance imaging allows to study the course of nerve fibre pathways in vivo and for whole brain volumes 1,2 , but the resolution is not sufficient to reconstruct individual nerve fibres 3,4 . Microscopy techniques like optical coherence tomography 5,6 , light-sheet microscopy 7,8 , or two-photon fluorescence microscopy 9,10 allow an in-depth scan of the sample and provide microscopic resolution in all three dimensions. However, they can only be applied to smaller brain volumes and need advanced image processing to follow the course of nerve fibres, which is especially challenging in regions with densely packed fibres. The microscopy technique 3D polarised light imaging (3D-PLI) 11,12 analyses the spatial course of nerve fibre pathways by transmitting polarised light through unstained histological brain sections and measuring their birefringence (optical anisotropy). In contrast to the other techniques, 3D-PLI requires sectioning of the brain and advanced registration algorithms to recover the original brain volume. But it allows to retrieve the three-dimensional orientations of the nerve fibres also in densely grown regions and for large brain sections, making 3D-PLI a powerful approach for analysing the brain's nerve fibre network in microscopic detail.
The three-dimensional nerve fibre orientations are derived from the measured birefringence signals, which are mainly caused by myelin-a stack of cell membranes with highly ordered molecular structure that surrounds most axons in the white matter 13 . (In the following, the term nerve fibre will solely be used to refer to myelinated axons.) The orientation of the nerve fibres within the section plane can be determined with high degree of accuracy, provided the fibre bundles have a well-defined orientation 14 . The computation of the out-of-plane fibre orientation (inclination) requires a more sophisticated analysis: The fibre inclination is directly related to the Methods Experimental methods. We demonstrate the functionality of our algorithm on coronal and sagittal sections of primate and rodent brains. To validate the automatically computed nerve fibre inclinations, we compared our results from 3D-PLI to in-depth-tissue scans from two-photon fluorescence microscopy.
Preparation of brain sections. The brain sections were obtained from three healthy rats (Wistar, male, 3 months old) and two healthy vervet monkeys (male, between 1 and 2 years old). All animal procedures were approved by the institutional animal welfare committee at Forschungszentrum Jlülich GmbH, Germany, and were in accordance with the European Union and National Institutes of Health guidelines for the use and care of laboratory animals and in compliance with the ARRIVE guidelines. Euthanasia of rats was carried out under controlled isoflurane inhalation followed by decapitation. The monkey brains were obtained in accordance with the Wake Forest Institutional Animal Care and Use Committee (IACUC #A11-219). Euthanasia procedures conformed to the AVMA Guidelines for the Euthanasia of Animals, using ketamine/pentobarbital anesthesia followed by perfusion with phosphate buffered saline and fixation with 4% paraformaldehyde.
The brains were removed from the skull within 24 h after death, fixed with 4% buffered formaldehyde for several weeks, cryo-protected with 20% glycerine and 2% dimethyl sulfoxide, deeply frozen, and coronally or sagittally cut with a cryostat microtome (Polycut CM 3500, Leica, Microsystems, Germany) into sections of 60 µm. A coronal vervet monkey brain section (no. 512), a sagittal vervet monkey brain section (no. 374), a coronal rat brain section (no. 185), and a sagittal rat brain section (no. 176) were selected for evaluation. The sections were mounted on glass slides, embedded in 20% glycerine solution, cover-slipped, sealed with lacquer, and measured with 3D-PLI afterwards.
3D polarized light imaging (3D-PLI). The 3D-PLI measurements were performed with the polarising microscope (LMP-1, Taorad GmbH, Germany), as described in Menzel et al. 14,19 : Incoherent light with a wavelength of about 550 nm was transmitted through a linear polariser, the brain section, and a circular analyser (quarter-wave retarder and linear polariser), see Fig. 1a. The first polariser was rotated in steps of 10° and for each rotation angle ( ρ = {0 • , 10 • , . . . , 170 • } ) an image was recorded with a CCD camera (Qimaging Retiga 4000R), yielding a pixel size in object space of about 1.33 µm. For each image pixel, the intensity values of the resulting image series form a sinusoidal signal with respect to the rotation angle (see Fig. 1b), which can be described as 11 www.nature.com/scientificreports/ The average of the signal (transmittance T ) is a measure of the tissue attenuation, caused by absorption and scattering. The phase of the signal (direction ϕ ) indicates the in-plane direction angle of the nerve fibres within the brain section. The amplitude of the normalised signal (retardation R = | sin δ| ) indicates the strength of birefringence and is related to the out-of-plane angle of the nerve fibres (inclination α ) via δ ∼ cos α 2 (see below). Figure 1c illustrates the definition of ϕ and α , Fig. 1d shows an example of the resulting transmittance and retardation images.
Two-photon fluorescence microscopy (TPFM). To validate the automatically computed nerve fibre inclinations, we compared the results from 3D-PLI to two-photon fluorescence microscopy (TPFM) measurements of the same brain section. The TPFM measurements were performed on a coronal rat brain section ( 6 × 8 tiles in the caudate putamen, cf. Fig. 6b) as described in Menzel et al. 19 , with a custom-made two-photon fluorescence microscope 9,22,23 consisting of a model-locked titanium-sapphire laser with a wavelength of 800 nm and a waterimmersion 25× objective lens (LD LCI Plan-Apochromat 25x/0.8 Imm Corr DIC M27), achieving a resolution of 0.244 × 0.244 × 1 µm 3 . The fluorescence signals were collected by two photomultiplier tubes, detecting red and green fluorescence, respectively. TPFM allows a full in-depth scan of the 60 µm thin brain section, yielding one image every z = 1 µm. The axial displacement was realised with a closed-loop piezoelectric stage. A motorised xy-stage enabled tile-wise scanning of the sample. The sample was measured in tiles of 250 × 250 µm 2 , with an overlap of 10%, and the resulting images were stitched together using the software TeraStitcher 24 . Due to the slightly different autofluorescence signal of different tissue components, the nerve fibre bundles can be manually traced across different slices of a TPFM image stack. To determine the inclination angle of a fibre bundle, the cross-sections of the fibre bundle were determined in the first and last slices of the image stack that still show the fibre bundle. The inclination angle α was computed geometrically from the midpoint positions of the crosssections, taking the number of images (z-depth) between the first and last slice into account. 230 fibre bundles in the caudate putamen were selected for evaluation. Only well-separated fibre bundles with well-defined fibre inclination were chosen. To enable a comparison between the geometrically computed TPFM inclinations and the automatically computed 3D-PLI inclinations, the transmittance image of the 3D-PLI measurement was registered onto the maximum intensity projection of the TPFM image stack using affine transformations, and the individual fibre bundles were masked in the registered transmittance image. To avoid registration artefacts, the inverse transformation was applied to the resulting fibre mask, allowing to locate the fibre bundles in the original (non-registered) inclination image and compare the 3D-PLI inclination to the corresponding TPFM inclination. To ensure that the tissue had not deformed between measurements, the maximum intensity projection of the TPFM image stack was compared to the transmittance image, ensuring that visible borders of nerve fibre bundles overlap. (1) (1 + sin δ sin (2(ρ − ϕ))). where is the wavelength of the incident light, d m the thickness of the birefringent tissue (myelin), and n the birefringence of the material. As these parameters are usually not precisely known 11 , the fibre inclinations can only be computed based on certain assumptions.
Unweighted model. Assuming that n d m does not change much within the measured brain section, the fibre inclination angle can be computed for each image pixel via: where R ref is the retardation expected to occur in regions with densely packed, parallel in-plane nerve fibres Transmittance-weighted model. Usually, the amount of nerve fibres (and thus the amount of birefringent material d m ) changes notably within a brain section. When a region with a large amount of nerve fibres is used as reference, the unweighted model leads to an over-estimation of the fibre inclination angles in other regions with a smaller amount of nerve fibres. To take the different amounts of nerve fibres and hence the different retardation values into account, the transmittance can be used as a reference. Due to the high scattering coefficient of myelin, regions with many nerve fibres appear darker than regions with less nerve fibres. When light with intensity I 0 is transmitted through a material with thickness d and attenuation coefficient µ , the resulting light intensity T is attenuated according to the Lambert-Beer law: In brain tissue, light passes through birefringent tissue (myelin) with attenuation coefficient µ m , and non-birefringent tissue (e.g. cell bodies) with attenuation coefficient µ c . For now, we assume that both attenuation coefficients do not depend on the exact tissue composition or the fibre inclination. The maximum thickness of birefringent tissue in the brain section is defined as d Combined model. As mentioned in the Introduction, the transmittance is only a good measure of the amount of nerve fibres (myelinated axons) if the attenuation of myelin is dominant compared to the attenuation of other tissue components. In regions with a high degree of myelination (HM-regions), the attenuation can be considered to be mostly caused by the scattering of myelin so that the fibre inclinations can be computed with the transmittance-weighted model. However, in regions with a low degree of myelination (LM-regions), the transmittance is dominated by other tissue components such as cell bodies so that the unweighted model should be used instead to compute the fibre inclinations.
To account for different tissue compositions, we defined P HM ∈ [0, 1] as the probability that a region is highly myelinated and computed the nerve fibre inclinations using a linear combination of both models, the transmittance-weighted model for HM-regions (Eq. 7) and the unweighted model for LM-regions (Eq. 3): Determination of the degree of myelination. To determine the degree of myelination ( P HM ) for all image pixels of a measured brain section, we first used a binary classification to separate LM-regions ( P HM = 0 ) from HM-regions ( P HM = 1 ). Based on this initial separation, we then defined transition zones with 0 < P HM < 1.
Binary classification. We followed the classification introduced in Benning et al. 21 to separate regions with low and high myelination (LM-and HM-regions). Note that the authors used the terms "grey matter" and "white matter", while we use the terms LM-and HM-regions to avoid confusion with the anatomical definitions. The separation of regions was performed by analysing the transmittance (T) and retardation (R) images from 3D-PLI.
While all regions with high birefringence ( R > R thres ) can be considered as HM-regions, not all regions with low birefringence ( R ≤ R thres ) can be considered as LM-regions because regions with crossing fibres, out-of-plane fibres, or a smaller number of nerve fibres also result in smaller birefringence values. As in-plane crossing and out-of-plane fibres have similar or even lower transmittance values than regions with parallel in-plane fibres 19 , the transmittance in a region with many parallel in-plane fibres can be used as a reference ( • HM-regions (high myelination): All three threshold parameters ( R thres , T thres , T back ) were computed from points of maximum curvature in the retardation and transmittance histograms (see blue/black vertical lines in Fig. 2b). To prevent the background from having any influence on the computed fibre inclinations, the retardation and transmittance of all background pixels were set to their minimum and maximum value, respectively, before analysing the corresponding histograms. Before computing the transmittance histogram, the values were normalised to [0,1] and a circular median filter with 5 px radius was applied to the transmittance image. We selected a radius of 5 px (and not 10 px as in Benning et al.) because a larger radius leads to more clouding artefacts in the resulting inclination map, while a smaller radius leads to increased noise and makes non-relevant details like cells more visible.
As in Benning et al. 21 , the threshold parameters R thres and T back were computed as the points of maximum curvature behind and before the biggest peak in the retardation and transmittance histogram, respectively. In contrast to Benning et al. 21 who used 128 bins, we decided to use histograms with 256 bins in order to increase accuracy. However, histograms with 256 bins might show several small peaks or a plateau (cf. Fig. 2b on the very left), which makes the determination of the threshold parameters difficult. To ensure that no outliers were chosen, we first considered a histogram with 64 bins and computed the threshold parameters (points of maximum curvature) in a range corresponding to the full-width at half-maximum (FWHM) of the peaks, using larger ranges as Benning et al. 21 to account for thinner peaks: 20 × FWHM behind the retardation peak and 10 × FWHM before the transmittance peak. Then, we considered histograms with 128 bins and recomputed the threshold parameters in a range of bins close to the previous results [ 2 · (result − 1), 2 · (result + 1) ]. The last step was repeated for histograms with 256 bins. In this way, the computed threshold parameters stay close to the previous results while being more accurate.
The threshold parameter T thres was computed as the point of maximum curvature between T ref and T back . The value for T ref was obtained by computing the average transmittance of the region with the highest retardation values, which is expected to contain the highest amount of parallel in-plane fibres (cf. Eq. 2). To avoid outliers, the union-find connected components algorithm by Oliveira et al. 26 was used to identify a region of connected pixels (0.009-0.011% of the image size) that has the highest retardation values.
Transition zones. When the inclination angles are computed separately for HM-and LM-regions, artificially sharp borders might occur, e.g. in regions where nerve fibres from the white matter spread out into the cortex. To account for regions that do not clearly belong to one or the other tissue type, we defined transition zones between HM-and LM-regions. www.nature.com/scientificreports/ To determine the width of these transition zones, the threshold parameters R thres and T thres were considered, which have a major impact on the separation of HM-and LM-regions (cf. striped blue areas in Fig. 2b): In a transition zone, slight changes of the threshold parameters cause the region to be identified as a different tissue type. To find out whether a separation is susceptible to small changes, the threshold parameters were recomputed several times from slightly different retardation and transmittance histograms that were generated by bootstrapping (random sampling with replacement): Pixels were randomly drawn from the retardation and transmittance images; the same pixel could be selected multiple times. To identify the best compromise between computing time and accuracy, the algorithm was run for different sample sizes (100%, 25%, 4%, 1% of all image pixels) for 500 iterations. From the resulting values for R thres and T thres , the positive/negative means ( R thres+/− and T thres+/− ) were computed, respectively. The positive (negative) mean is the mean of all recomputed threshold values that are larger (smaller) than the original threshold value. The whole process was repeated 10 times and the average values for each iteration were forwarded to further analysis. Figure 2c shows the positive and negative means for the different sample sizes and up to 500 iterations for a coronal rat brain section. The different sample sizes yield similar curves. The smaller the sample size, the more the values differ from the ones obtained with 100% sample size. After 200 iterations (vertical dashed line), the values for 100% and 25% sample size (cyan and magenta curves) do not change much and the maximum difference between 25% and 100% sample size is still within one histogram bin ( 1/256 = 0.0039 ). Therefore, a sample size of 25% and 200 iterations were selected for computing R thres+/− and T thres+/− , which define the width of the transition zone between LM-and HM-regions.
HM-probability map. To realize a smooth transition between regions with low myelination ( P HM = 0 ) and high myelination ( P HM = 1 ), the width of the transition zone was described by a sigmoid function (cf. top image in Fig. 3c).
The result is an HM-probability map which indicates the probability P HM that a region is highly myelinated (HM-region). In this representation, regions with P HM > 0.95 were considered as HM-regions, regions with While there is a sharp separation between fibre bundles and surrounding tissue in the putamen (yellow arrow), there is a smooth transition for fibres spreading from the white matter into the cortex (magenta arrow). The slightly lower P HM -values in the corona radiata (cr), which is a region rich of white matter tracts, are due to the steep and crossing nerve fibres which impair the classification of regions (see "Discussion").  As out-of-plane nerve fibres have lower transmittance values than in-plane nerve fibres due to an increased amount of scattering (see Menzel et al. 19 ), all (median-filtered) transmittance values smaller than T M (transmittance of the region with the highest retardation values, expected to contain densely packed nerve fibres) were clipped to this value, so that T ≥ T M . Figure 4b shows the median-filtered transmittance image (enlarged view of the region marked by the yellow rectangle) without this correction (top) and with correction (bottom). Without correction, the transmittance values of the cingulum (cg), which contains mostly very steep out-of-plane nerve fibres, are notably darker than those of the corpus callosum (cc), which contains mostly in-plane nerve fibres. The same applies to the fornix (f) and the optic tract (ot).
Software. The HM-probability maps and the nerve fibre inclinations were computed with the specifically developed open-source software PLImig. The algorithm is parallelised and uses GPU acceleration, allowing for a fast and reproducible computation of nerve fibre inclinations in large brain sections. For more information, the reader is referred to our GitHub page (https:// github. com/ 3d-pli/ PLImig). Depending on the available graphics memory, the images were processed in a different number of chunks to avoid overflow. The coronal and sagittal rat brain sections (417 MB and 843 MB, using 32-bit float values) were computed using an NVIDIA RTX 3070 and 8 GB of RAM. The coronal and sagittal vervet brain sections (6 GB and 6.5 GB) were processed on the supercomputer JURECA 27 . Using 512 GB RAM, one GPU (NVIDIA A100 40GB) and 128 CPU cores on one node, the computing time for the coronal vervet brain section (34,669 × 46,341 pixels) was 15 min for the HM-probability map and 5 min for the rest.  www.nature.com/scientificreports/ to whole brain section, (ii) transmittance-weighted model ( α HM ) applied to whole brain section, (iii) unweighted and transmittance-weighted models applied to LM-and HM-regions individually, and (iv) unweighted and transmittance-weighted models applied to LM-and HM-regions using the HM-probability map combining both formulas in the transition zone.

Comparison of different models for computing the fibre inclination angle in
When applying the unweighted model to the whole image, using the highest retardation of the brain section as reference ( R ref ,LM = R ref ), the fibre inclinations in low myelinated regions (e.g. cortex) are highly overestimated (Fig. 5i). On the other hand, when applying the transmittance-weighted model to the whole image ( R ref ,HM = R ref ), the fibre inclinations in these regions are under-estimated, generating a lot of saturated values at α = 0 • in black (Fig. 5ii). When applying the two different models to HM-and LM-regions individually, using their highest retardation values as reference ( R ref ,HM , R ref ,LM ), the estimated fibre inclinations become more realistic throughout the entire brain section. The contrast between out-of-plane fibres (cingulum = cg) and in-plane fibres, both in regions with densely packed nerve fibres (corpus callosum = cc) and in regions with bulked fibre architecture (cortex), becomes clearly visible (Fig. 5iii). However, the sharp separation of HM-and LM-regions leads to an underestimation of the inclination (dark pixels) in regions where fibres extend to cortical areas (see yellow arrow in (III)), even worsening where sharp bending is involved (cyan arrows in (III)). The fourth model using the HM-probability map to combine the formulas for HM-and LM-regions in the transition zone reduces the underestimation considerably, leading to a more continuous course of fibres from deep white matter regions into the cortex (see yellow and cyan arrows in (IV)), corresponding to the anatomical expectation.

Comparison of fibre inclinations obtained from 3D-PLI and TPFM.
To quantitatively evaluate the automatically computed nerve fibre inclinations, the caudate putamen of a coronal rat brain section was measured both with 3D-PLI and with two-photon fluorescence microscopy (TPFM) as described in the Methods.
As TPFM allows an in-depth scan of the sample, it allows to directly assess the three-dimensional course of nerve fibre bundles in the caudate putamen (cf. Fig. 6c), so that it can serve as a comparison for the nerve fibre orientations determined by 3D-PLI. 230 fibre bundles with well-defined shape and orientation were selected for evaluation (Fig. 6a,b). The average 3D-PLI inclination of each bundle was plotted against the TPFM inclination, which was geometrically reconstructed from the corresponding bundle (see Methods). Figure 6d,e show the results for the in-plane fibre direction angle ϕ and the out-of-plane fibre inclination angle α for the four different scenarios displayed in Fig. 5.
While the direction angles computed from 3D-PLI correspond mostly to those obtained from TPFM (see Fig. 6d), the computed 3D-PLI inclination angles are generally larger than those obtained from TPFM (see Fig. 6e). When applying the unweighted model to the whole brain section (i), the computed 3D-PLI inclinations are highly over-estimated, especially for lower fibre inclinations ( < 50 • ). The computed 3D-PLI inclinations range between 50 • and 70 • , mostly independent of the actual inclination of the nerve fibres. When applying the transmittance-weighted model to the whole brain section (ii), the computed 3D-PLI inclinations become much more similar to those obtained from TPFM. However, even when taking their standard deviations ( �α ≈ 30 • ) into account, most values do not reach the TPFM inclinations. When applying the unweighted and transmittanceweighted models individually to LM-and HM-regions regions (iii), a global shift of the 3D-PLI inclinations www.nature.com/scientificreports/ of about −5 • can be observed. When using a linear transition between the models (iv), the computed 3D-PLI inclinations are still over-estimated, but the difference to the fibre inclinations derived from the TPFM measurements becomes slightly smaller at higher inclinations.
Fibre inclinations in coronal and sagittal brain sections. Figure 7 shows the HM-probability maps and the corresponding inclination maps (computed from Eq. 8) for coronal and sagittal sections of vervet monkey and rat brains. The HM-probability maps reliably separate regions with low myelination (black) from the rest. Individual nerve fibres are visible in the putamen (Pu) of the coronal vervet brain section, and in the caudate putamen (CPu) of the coronal and sagittal rat brain sections. The HM-probability values in the white matter are slightly reduced (light green) for in-plane crossing fibres (cr) and steep out-of-plane fibres (coronal: cg, f; sagittal: cc), yielding  19 , Fig. S2c). The threedimensional fibre orientation (dashed yellow line, defined by the in-plane direction angle ϕ and the out-ofplane inclination angle α ) was geometrically computed from the centres of gravity of the bundle cross-sections mentioned above (yellow shapes). (d) Direction angles of the fibre bundles derived from the aligned 3D-PLI and TPFM measurements. The grey bars denote the positive and negative standard deviation when averaging the pixel values over the selected regions in the 3D-PLI direction image. (e) Inclination angles of the fibre bundles derived from the aligned 3D-PLI and TPFM measurements. The 3D-PLI inclinations were computed for four different scenarios (cf. Fig. 5): (i) unweighted model ( α LM ) applied to whole brain section, (ii) transmittanceweighted model ( α HM ) applied to whole brain section, (iii) unweighted and transmittance-weighted models applied to LM-and HM-regions individually, (iv) unweighted and transmittance-weighted models applied to LM-and HM-regions using a linear combination of both formulas in the transition zone. www.nature.com/scientificreports/ slightly lower estimations of the inclination angles (cf. Fig. 6e). The inclination maps show a good distinction of regions with in-plane and out-of-plane nerve fibres. Although the coronal and sagittal sections belong to different brains, the different section planes (coronal/sagittal) clearly show an inverse pattern: In the coronal section plane, the computed inclination angles of the corpus callosum (cc) are much lower than those in the cingulum (cg) or the fornix (f). In the sagittal section plane, it is exactly the other way around.

Discussion
When analysing measurement data from 3D-polarised light imaging (3D-PLI) of histological brain sections, the determination of the out-of-plane nerve fibre inclinations has been a major challenge, requiring anatomical knowledge and effortful manual adjustment of parameters. Here, we introduced an automatic computation of the fibre inclinations, by evaluating the retardation and transmittance images (normalised amplitude and average of the measured 3D-PLI signal). As these images are generated during standard 3D-PLI measurements, past measurements can be easily evaluated to improve the interpretation of fibre orientations in 3D-PLI data.
In particular, no special equipment like a tilting stage is required to compute the nerve fibre inclinations. In contrast to the manual adjustment, no expert knowledge is required, the results for each measured brain section are reproducible, and large data sets can be analysed. In this way, we were able to provide a reliable estimate of the nerve fibre inclinations, which can still be adjusted by an expert in a second step.
In addition to the fibre inclinations, we provided an HM-probability map, indicating the probability that a region is highly myelinated. Benning et al. 21 proposed a binary separation of regions with low and high myelination. We improved the existing algorithm (taking different features/plateaus of the histograms into account and excluding background tissue when determining the threshold parameters). Furthermore, we introduced the important concept of transition zones. The classification of brain regions into regions with high myelination (HM-regions), low myelination (LM-regions), and transition zones is not only relevant when applying different inclination formulas. Another possible application might be Independent Component Analysis (ICA), which has been used for noise and artefact removal in 3D-PLI images 28 , especially in regions with low myelination (low birefringence signals) 21 . The HM-probability map allows to determine different degrees of myelination, helping to improve the ICA in those regions and enabling enhanced tissue segmentation.
The HM-probability maps and fibre inclinations were computed with a specifically developed open-source software, which is parallelised and uses GPU acceleration to efficiently process large images. While the computation of the fibre inclinations takes less than 5 min for the coronal vervet brain section (34 669 × 46 341 px), the generation of the HM-probability map takes about 15 min. The sampling consumes a major part of the computing time, but it is much more efficient than alternative approaches, e.g. defining transition zones by the local variance of pixel values (requiring to study many surrounding pixels).
The choice of the reduced sample size and number of iterations (cf. Fig. 2c) has a large impact on the computing time. When using 100% instead of 25% sample size, four times more computing time (i e. more than 1 h) is required. To further reduce computing time, the number of iterations can be reduced to 100 without noticeable deviation of the threshold parameters from their asymptotic optimum (cf. Fig. 2c). Also, a more advanced bootstrapping approach (requiring smaller sample sizes) could be employed. The program could be extended to run on multiple GPUs, using the chunking algorithm that is already present. This would allow to process even larger brain sections and ideally achieve a linear speed-up.
While the software was shown to yield reliable results and correctly reproduce the overall inclination of major nerve fibre structures, several aspects should to be taken into account. First of all, the transmittance-weighted model used to compute the fibre inclinations in the HM-regions depends on the transmittance values. When the transmittance image has not enough contrast (e.g. in brain sections measured several days/weeks after tissue embedding) or when the transmittance image shows irregularities (caused e.g. by an inhomogeneous distribution of the embedding glycerine solution or tile/stitching artefacts), this leads to artefacts in the resulting inclination map, which would not be visible when solely using the birefringence signal (retardation) for computing the inclination, as in the LM-regions. Hence, the software should only be used when the transmittance image contrast is free of artefacts. Another point to take into account is that the algorithm uses the transmittance and retardation histograms to compute the parameters in the inclination formula. A reliable computation of parameters only works when the histograms include a broad variety of tissue components as well as various nerve fibre structures (in-plane and out-of-plane fibres).
While the computed inclinations can serve as an estimate of the real out-of-plane angle in homogeneous nerve fibre tissue, they are highly over-estimated in regions with in-plane crossing nerve fibres due to the reduced birefringence signal 15 . As shown in Menzel et al. 19 , the transmittance image can serve as an indicator whether large inclination angles belong to out-of-plane nerve fibres or are caused by in-plane crossing fibres or a low fibre density. Regions with high inclination values can be considered as reliable if the transmittance in these regions is notably lower than in regions with in-plane nerve fibres. When comparing the transmittance image (Fig. 4a) with the inclination image (Fig. 7a) for the coronal vervet brain section, it becomes apparent that cingulum and fornix fulfill this criterion (i e. their transmittance is notably lower than the transmittance of the corpus callosum which contains mostly in-plane nerve fibres). The corona radiata is dotted with regions that show substantially darker transmittance values than the corpus callosum (caused by the superior longitudinal fascicle -a steep fibre tract crossing this region), but it also contains regions with similar transmittance values, indicating that the high inclination values in these areas cannot be considered as reliable and are most probably caused by in-plane crossing nerve fibres. As the threshold parameters in the transmittance histogram are very susceptible to small changes, we decided to keep an automated detection of crossing nerve fibres for future algorithms.
So far, the nerve fibre inclinations were computed by applying a transmittance-weighted model to the whole brain section and manually adjusting the parameters. The models presented here (with and without transmittance www.nature.com/scientificreports/ weighting in HM-and LM-regions, Eqs. 3 and 7) lead to a major enhancement of the computed fibre inclinations in regions with low myelination (cf. Fig. 5(ii) and (iii)). The computation of an HM-probability map allowed for the first time to consider transition zones, reducing step artefacts at the boundaries between HM-and LM-regions (cf. Fig. 5(III) and (IV)): The course of nerve fibres from highly myelinated regions into the cortex becomes much more continuous (yellow arrows), even border regions with sharp bending that were clearly visible as dark stripes become more softened (cyan arrows). However, the dark fissures of the inclination at the borders are still present in some regions and still need to be improved by another tuning of R thres+/− and T thres+/− . As shown in Fig. 7, the consistency of nerve fibre inclinations in whole coronal and sagittal brain sections of different species (monkeys and rodents) validates the reliability of the algorithms applied. The findings presented in Fig. 6 show that the different models have different impacts on the computed inclinations of nerve fibre bundles (in the caudate putamen of a rat brain). When using the unweighted model (i), the computed inclinations are mostly independent of the real underlying fibre inclination. The transmittanceweighted model (ii) introduces the proportionality of the computed inclinations (they increase with increasing fibre inclination) in particular for values below 50 • . The introduction of the HM-probability map (iv), necessary to generate valid fibre inclinations in regions with low myelination and to obtain a smooth transition at the boundaries, has only a small impact on the computed fibre inclinations of the nerve fibre bundles: (iv) almost shows a similar behaviour as (ii). Hence, the linear combination of both models (Eq. 8) to compute the nerve fibre inclinations in regions with different degrees of myelination can be used without changing the inclination values in highly myelinated regions. Although the transmittance-weighted inclinations (ii-iv) are much more realistic than those computed with the unweighted model (i), the comparison of 3D-PLI and two-photon fluorescence microscopy (TPFM) shows a strong over-estimation of the computed fibre inclinations ( �α ≈ 30 • ). However, the inclinations obtained from 3D-PLI might also be larger than those obtained from TPFM because the TPFM measurements have been performed several weeks afterwards so that the thickness of the brain section might have decreased due to dehydration. For this reason, a direct comparison of 3D-PLI and TPFM values is critical. However, as the affine image registration of 3D-PLI to TPFM measurements did not reveal large tissue deformations, maintaining also a high affinity of the corresponding direction angles, major deviations of the fibre orientations are unlikely.
In conclusion, we have introduced a versatile software that allows for the first time to automatically compute the out-of-plane angles of nerve fibres obtained from standard 3D-PLI measurements of brain tissue (without tilting). So far, these nerve fibre inclinations were computed by applying a transmittance-weighted model to the whole brain section, requiring anatomical knowledge and effortful manual adjustment of parameters. By distinguishing areas with different degrees of myelination and defining transition zones, we were able to apply a regionally specific computation of fibre inclinations that accounts for subtle changes in the tissue composition. The developed algorithm is parallelised and uses GPU acceleration, allowing for a fast and reproducible computation of nerve fibre inclinations in large brain sections, thus greatly improving the analysis and interpretation of 3D-PLI data also for past measurements.

Data availability
All data supporting the findings of this study are included in the provided figures.

Code availability
The software PLImig that was used to compute the HM-probability maps and the nerve fibre inclinations in Figs. 3, 5 and 7 is publicly available on GitHub (https:// github. com/ 3d-pli/ PLImig). Data analysis was performed with the open-source software Fiji (https:// fiji. sc/ Fiji). License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.