## Abstract

The neuroimaging technique three-dimensional polarized light imaging (3D-PLI) provides a high-resolution reconstruction of nerve fibres in human post-mortem brains. The orientations of the fibres are derived from birefringence measurements of histological brain sections assuming that the nerve fibres—consisting of an axon and a surrounding myelin sheath—are uniaxial birefringent and that the measured optic axis is oriented in the direction of the nerve fibres (macroscopic model). Although experimental studies support this assumption, the molecular structure of the myelin sheath suggests that the birefringence of a nerve fibre can be described more precisely by multiple optic axes oriented radially around the fibre axis (microscopic model). In this paper, we compare the use of the macroscopic and the microscopic model for simulating 3D-PLI by means of the Jones matrix formalism. The simulations show that the macroscopic model ensures a reliable estimation of the fibre orientations as long as the polarimeter does not resolve structures smaller than the diameter of single fibres. In the case of fibre bundles, polarimeters with even higher resolutions can be used without losing reliability. When taking the myelin density into account, the derived fibre orientations are considerably improved.

## 1. Introduction

Unravelling the architecture and connectivity of nerve fibres in the human brain is one of the greatest challenges in neuroscience. Over the past years, several methods have been developed to reconstruct the human connectome [1–3]. The neuroimaging technique three-dimensional polarized light imaging (3D-PLI) has been employed to reconstruct the three-dimensional architecture of nerve fibres in human post-mortem brains with a resolution of a few micrometres [4,5]. 3D-PLI enables the investigation of the pathways of long-range fibre bundles as well as single fibres and thus serves as a bridging technology between the macroscopic and the microscopic scale.

The spatial orientations of the nerve fibres are derived by transmitting polarized light through histological brain sections in a polarimeter and measuring their birefringence. To relate the measured signal to the fibre orientation, an *effective model* of birefringence is used which assumes that the fibre density is constant over the whole brain section [4] and that the measured optic axis indicates the predominant fibre orientation [5,6]. This assumption is based on various experimental studies on white matter which show that the average birefringence of parallel nerve fibres is negatively uniaxial and that the measured optic axis is oriented along the length of the fibres [7–10].

The majority of nerve fibres in the brain consist of an axon and a surrounding myelin sheath. The cytoplasm of the axon contains tubular polymers (microtubules) and neurofilaments running along the length of the axon [11,12]. The myelin sheath is formed by oligodendrocytes (glial cells) which are spirally wrapped around the axon. The cell membranes are bimolecular layers consisting of lipid molecules and membrane proteins. The membrane proteins are embedded in the bilayer or attached to the membrane surface [13–15], whereas the lipid molecules are oriented radially to the fibre axis [15–17]. The cell organelles of the axon and the protein framework of the myelin sheath lead to a weak positive birefringence with respect to the longitudinal fibre axis [7,8,10,17–20]. The anisotropic structure of the lipid molecules causes a positive birefringence with respect to the radial fibre axis [7,8,15,18,21].

The effective model of uniaxial negative birefringence that is currently used in 3D-PLI seems reasonable for sufficiently low optical resolutions. However, it might no longer be valid if the anisotropic molecular structure of the nerve fibres is resolved. In this paper, we investigated the limitations of the effective model in terms of the optical resolution of the polarimeter using numerical simulations. The simulations were performed with a modified version of SimPLI [6], a simulation method that models the birefringence of the fibres with the Jones matrix calculus and allows data to be generated from synthetic fibre constellations that is comparable to experimental data. In order to study and understand the most dominant effects that generate the birefringence signals in 3D-PLI, the anisotropic molecular structure of the nerve fibres was described by a simplified birefringence model with radial optic axes (*microscopic model*) and the effective model of uniaxial negative birefringence by a birefringence model with axial optic axes (*macroscopic model*). To investigate the limitations of the effective model, the transition between the microscopic and the macroscopic model was investigated depending on the optical resolution of the imaging system.

## 2. Three-dimensional polarized light imaging

The neuroimaging technique 3D-PLI determines the orientation of nerve fibres in post-mortem brains at the micrometre scale. The principles of 3D-PLI have been explained in detail by Larsen *et al.* [22] and Axer *et al.* [4,5]. This section describes the measurement and data analysis procedures that are relevant for this study.

### 2.1. Measurement

To determine the orientation of the nerve fibres, a post-mortem brain—obtained from a body donor in accordance with ethical requirements—is fixed in buffered formaldehyde for several months, frozen and cut with a cryotome into histological sections of 70 µm, which are measured with a polarimeter. For the 3D-PLI measurement, two state-of-the-art polarimeters with different optical resolutions and sensitivities are employed: the large-area polarimeter (LAP) has a pixel size of 64 µm and is mainly used for single-shot images of whole human brain sections. The polarizing microscope (PM) has a pixel size of 1.33 µm (i.e. down to small axonal diameters), which enables complex fibre constellations to be disentangled.

The LAP contains a pair of crossed linear polarizers and a quarter-wave retarder (with its fast axis adjusted at an angle of −45° with respect to the transmission axis of the first linear polarizer) (figure 1*a*). The employed light source emits incoherent, non-polarized, diffusive light with a peak wavelength of 525 nm. During the measurement, the polarizers and the quarter-wave retarder are rotated simultaneously around the stationary tissue sample. For each rotation angle *ρ* = 0°, 10°, … , 170°, the transmitted light intensity is recorded by a CCD camera so that a series of 18 images is acquired.

The imaging principle works as follows: the quarter-wave retarder transforms the linearly polarized light from the first polarizer into circularly polarized light. The birefringent brain tissue induces an additional phase shift so that the outgoing light is elliptically polarized. The fraction of light that then passes the second linear polarizer depends on the local orientation of the optic axis of the birefringent tissue, which is assumed to coincide with the local fibre orientation.

The polarimetric set-up of the PM is slightly different to the set-up of the LAP (the order of the optical elements is reversed and only the first linear polarizer is rotatable). However, the imaging principle and the signal analysis are similar [4] so that the following considerations are only described for the LAP.

### 2.2. Signal analysis

The measured light intensity of an individual pixel describes a sinusoidal curve across the acquired image series, which depends on the orientation of the fibres within this pixel (figure 1*b*). A physical description of the measured light intensity profile can be derived with the Jones matrix calculus [23,24], assuming that the light is coherent and completely polarized and that the optical elements are linear. For simplicity, the derivation is shown for a single pixel at a certain rotation angle *ρ*.

In the Jones matrix calculus, all optical elements in the polarimeter are represented by Jones matrices (cf. figure 1*a*). The Jones matrices of the crossed linear polarizers are given by [25]
2.1

The Jones matrix of a wave retarder that is rotated by an angle *ψ* in counterclockwise direction and induces along the fast axis a phase shift *δ* between the two orthogonal components of the light wave is given by [25]
2.2

In the experimental set-up, the fast axis of the quarter-wave retarder is rotated by −45° with respect to the axis of the first linear polarizer. Thus, the quarter-wave retarder can be described by the Jones matrix of a rotated wave retarder as given in equation (2.2) with a rotation angle of *ψ* = −45° and a phase shift of *δ* = 90°:
2.3

Under the assumption that the birefringence of the brain tissue can locally be described as negatively uniaxial with the optic axis indicating the predominant fibre direction (effective model), the brain tissue can locally be represented by a wave retarder that introduces a phase shift *δ* along the fast axis (fibre axis). During the measurement, the two polarizers and the quarter-wave retarder are rotated simultaneously around the specimen stage in counterclockwise direction by a rotation angle *ρ*. For simplicity, the equivalent case is considered in which the brain tissue is rotated by an angle (−*ρ*) in counterclockwise direction while the other optical elements are fixed. Thus, the brain tissue can be described by the Jones matrix of a rotated wave retarder as given in equation (2.2) with phase shift *δ* and rotation angle *ψ* = *φ* − *ρ*, where *φ* denotes the in-plane orientation of the optic axis:
2.4

When light with an electric field vector *E*_{0} passes through the 3D-PLI set-up, the resulting output beam with electric field vector *E*_{T} can be described by multiplication of the associated Jones matrices. As the Jones matrix calculus cannot be used to describe the non-polarized light emitted by the employed light source, the Jones vector is used to describe the horizontally polarized light after the first linear polarizer (cf. figure 1*a*):
2.5

Using , the transmitted light intensity is calculated, yielding a sinusoidal intensity profile (figure 1*b*):
2.6where corresponds to the transmitted light intensity averaged over all rotation angles (here referred to as *transmittance*) and to the peak-to-peak amplitude of the normalized sinusoidal intensity profile (here referred to as *retardation*). The phase shift *δ* is given by (appendix A)
2.7where *λ* is the wavelength of the incident light, *t* the thickness of the brain section, Δ*n* the local birefringence of the sample and *α* the local out-of-plane inclination angle of the fibre. Thus, the intensity profile in equation (2.6) is a direct measure of the spatial fibre orientation defined by the direction angle *φ* and the inclination angle *α* (figure 1*c*).

In order to compute transmittance, direction and retardation, the intensity profile is fitted by means of a discrete harmonic Fourier analysis [4,26]. The inclination angle *α* is calculated from the measured retardation by rearranging equation (2.7). The direction and inclination angles are combined to a unit vector indicating the local fibre orientation in three dimensions. Putting all unit vectors of several adjacent brain sections together, a three-dimensional volume of vectors is created and the fibre tracts are reconstructed with streamline algorithms.

## 3. Simulation of three-dimensional polarized light imaging using the Jones matrix formalism

### 3.1. Simulation model

3D-PLI derives the nerve fibre orientations based on the fact that the average birefringence of parallel fibres is negatively uniaxial [7–10] and assuming that the orientation of the measured optic axis corresponds to the local fibre orientation. To investigate the limitations of this effective birefringence model, a straight single fibre and a hexagonal bundle of straight parallel fibres were simulated and the birefringence of the fibres was modelled according to a microscopic and a macroscopic model for different optical resolutions of the simulated imaging system.

#### 3.1.1. Microscopic model

The microscopic model of birefringence considers the anisotropic molecular structure of a single nerve fibre. To investigate and understand the predominant effects generating the birefringence signals in 3D-PLI, a simplified model of birefringence was chosen for the simulations. As stated in §1, the average birefringence of parallel nerve fibres is negative with respect to the longitudinal fibre axis. Therefore, the positive birefringence of the axon and the myelin proteins is weak as compared to the birefringence effects of the myelin lipids [8,9,18,20,21]. Since the exact contribution of the different birefringence effects to the overall birefringence is unknown, the birefringence effects of the nerve fibres were modelled by considering only the anisotropic radial structure of the myelin sheaths: the fibres were simulated as hollow tubes (representing the myelin sheaths) with positive birefringence and radial optic axes (cf. figure 2*b*(ii)). The axons were considered to be non-birefringent.

#### 3.1.2. Macroscopic model

To compare the simulation results of the microscopic model with the effective model of uniaxial negative birefringence, a macroscopic model of birefringence was defined. According to the assumptions made in the effective model, a single nerve fibre was simulated as negatively birefringent with axial optic axes oriented along the length of the fibre (cf. figure 2*b*(i)). Dohmen *et al.* [6] used this simulation model to investigate the effect of crossing fibre constellations. As this study concentrates on straight parallel fibres, the macroscopic model only serves as a reference for the effective model to verify the simulations of the microscopic model. To ensure a better comparison with the microscopic model, the fibres were simulated as hollow tubes (and not as solid cylinders as in [6]).

### 3.2. Simulation method

The basic idea of the simulation method is to model the birefringent myelin sheaths as series of linear optical retarder elements which are represented by Jones matrices. By defining the direction of the optic axes (radial/axial), both the microscopic and the macroscopic model can be simulated. The simulation approach is based on the simulation tool SimPLI developed by Dohmen *et al.* [6]. For this study, the simulation tool was extended by the microscopic model and modified such that various fibre configurations with individual orientations, radii and myelin sheath thickness can be realized.

The simulation tool is based on several assumptions and simplifications: first of all, the use of the Jones matrix calculus requires linear optical elements and perfect polarizers (i.e. the outgoing light is assumed to be completely polarized). Another assumption is that the incident light can be described by parallel rays of light with straight optical pathways, i.e. the light is assumed to be non-diffusive and refraction, diffraction and scattering are neglected. For this study, a parallel and straight beam of light seems a reasonable approximation for the LAP because the imaging system has a small numerical aperture (the acceptance angle of the objective lens is less than 1°) so that the camera only captures light rays that are almost parallel to each other.

The simulation consists of several steps:

(1)

*Generation of synthetic nerve fibres in a three-dimensional volume*. The nerve fibres are modelled as hollow tubes representing the myelin sheaths (figure 2*a*). In order to approximate the geometry of the fibres, the simulation volume is discretized into small cubic volume elements (called*voxels*), as indicated schematically by the grid in figure 2*c*.(2)

*Generation of a three-dimensional vector field*. For sufficiently small voxel sizes, the birefringence of the myelin sheaths can approximately be described by assigning each myelin voxel*j*a unit vector that indicates the direction of the optic axis (*φ*,_{j}*α*) within the myelin sheath. In the macroscopic model, the vectors are oriented parallel to the fibre axis. In the microscopic model, the vectors are oriented radially to the fibre axis (figure 2_{j}*b*).(3)

*Generation of a synthetic 3D-PLI image series*. In order to model the birefringence effect of the myelin sheaths, each myelin voxel is represented by the Jones matrix of a rotated wave retarder. The retarder axis is aligned with the optic axis within the myelin voxel (figure 2*c*). The synthetic 3D-PLI image series is calculated analogously to the derivation of the sinusoidal intensity profile as given in equation (2.5), with*M*_{tissue}being replaced by the product of*N*matrices representing the myelin voxels along the optical path: 3.1The matrix is the Jones matrix of a rotated wave retarder as given in equation (2.2) and represents the*j*th myelin voxel. The rotation angle depends on the in-plane direction angle*φ*of the optic axis and the phase shift_{j}*δ*on the out-of-plane inclination angle_{j}*α*. The Jones matrices of the linear polarizers and the quarter-wave retarder are given by equations (2.1) and (2.3). For each rotation angle of the polarimeter (_{j}*ρ*= 0°, 10°, … , 170°), all Jones matrices along the optical path are multiplied (figure 2*c*), yielding a series of 18 synthetic 3D-PLI images with a sinusoidal intensity profile for each image pixel.

### 3.3. Simulation parameters

The choice of the simulation parameters was inspired by real experimental conditions. According to typical dimensions of large nerve fibres in human white matter [16,27–29], the diameter and the myelin sheath thickness of the simulated nerve fibres were chosen to be 15 µm and 2.5 µm, respectively (figure 3*a*). The fibres were generated in a simulation volume with dimensions *x* × *y* × *z* = 64 × 64 × 70 µm^{3}, corresponding to the pixel size of the LAP and the thickness of the brain section. The simulation volume was discretized into cubic voxels with a side length of Δ*x*_{sim}. In a preliminary study (see later, §4.1), the optimal voxel size was determined to be Δ*x*_{sim} = 0.1 µm, which was used for all following simulations. Note that the dimensions are given in micrometres to meet the experimental conditions. As only relative length scales matter for the qualitative simulation results, the units could be chosen arbitrarily.

Since measuring the birefringence of the micrometre-thick brain sections is impossible with the employed set-ups and literature values are not given for the currently used preparation technique, an upper limit for the birefringence of the myelin sheaths Δ*n* was estimated: under the assumption that a brain section that is completely filled with a homogeneous birefringent material with in-plane optic axis (*α* = 0°) induces a maximum possible retardation the upper limit of the birefringence was calculated by rearranging equation (2.7): Note that the choice of Δ*n* only changes the overall magnitude of the retardation and does not affect the simulation results qualitatively. In the macroscopic (microscopic) model, the myelin voxels were simulated with axial (radial) optic axes and negative (positive) birefringence with respect to the optic axes.

The wavelength of the incident light was chosen to correspond to the peak wavelength of the LAP (*λ* = 525 nm). To study only the birefringence effect of the nerve fibres, the fibres were simulated without any absorption.

### 3.4. Simulation of the optical resolution

To investigate the effect of different optical resolutions on the measured 3D-PLI signal, the synthetic 3D-PLI image series were downsampled using the open-source image processing programme Fiji [30]: to account for the limited optical resolution of the polarimeter, the image series were first convoluted with a two-dimensional Gaussian filter with a standard deviation *σ*. Then, the effect of the spatial discretization of the CCD chip was modelled by resampling the resulting images with a sampling factor *f*_{s} (average when downsizing). To determine realistic parameters for *σ* and *f*_{s}, the imaging properties of the LAP were considered as a point of reference (appendix B).

Based on these considerations, the synthetic 3D-PLI image series were downsampled with different parameter sets (table 1), yielding images with different pixel sizes Δ*x*. The pixel size of the downsampled images was chosen such that a multiple of the pixel size corresponds to the side length of the simulation volume (Δ*x* = 64 µm/*n*, with *n* = 4, 8, 16, 32). The standard deviation was calculated as a linear function of the pixel size (*σ* = 0.714 Δ*x*, appendix B) and the sampling factor was calculated by dividing the pixel size of the high-resolution image series by the pixel size of the downsampled image . In the following, the optical resolution of the imaging system will be given in terms of the pixel size, which defines the set of downsampling parameters (*σ* and *f*_{s}) in table 1. Note that the simulation results will not change qualitatively as long as the ratio between the fibre dimensions and the downsampling parameters remains the same.

### 3.5. Calculation of the retardation curve

The determination of the inclination angle *α* is challenging for 3D-PLI because the peak-to-peak amplitude of the measured intensity profile is highly sensitive to noise and—among others—influenced by the density of myelinated nerve fibres (see below).

In the standard 3D-PLI analysis, the inclination angle is calculated from the measured intensity profile assuming that the brain tissue can locally be described by the effective model of uniaxial negative birefringence. In order to investigate whether the effective model can be used to extract the correct fibre inclinations, the retardation computed from equation (2.7) was compared to the retardation values derived from simulations using the macroscopic and the microscopic model (see §3.1). For that purpose, the retardation images were calculated for different fibre inclinations and different optical resolutions, respectively. For a better comparison between the retardation values of the single fibre and the fibre bundle, only the pixel in the centre of each (downsampled) retardation image was considered for evaluation. If pixels at other locations had been chosen, the retardation values of the single fibre would have been influenced by boundary effects that do not exist for the fibre bundle or real brain tissue which are completely filled with fibres. The retardation values from the centre of each downsampled retardation image were plotted against the corresponding inclination angle, yielding a *retardation curve* for each downsampling step. The retardation curves were compared to the normalized retardation curve of the effective model (cf. equation (2.7)), in the following referred to as the *theoretical curve*: .

To be able to compare different retardation curves, the retardation was normalized for each curve with the maximum retardation value, respectively: 3.2

As only birefringent material (mainly myelin) is responsible for the phase shift in equation (2.7), *t* describes not the thickness of the whole brain section but rather the *local myelin thickness t*_{m}, i.e. the combined thickness of myelin sheaths along the optical path. Due to the inhomogeneity of brain tissue, the local myelin density of a brain section is less than 100% i.e. the maximum possible retardation is . If the inclination is calculated under the assumption that the maximum possible retardation equals 1, the inclination angle will be overestimated. In order to obtain a more precise estimation of the inclination angle, a so-called *myelin density correction* was applied to the downsampled retardation images.

In the case of the macroscopic model, in which the optic axes within one fibre have the same orientations, *δ* scales linearly with *t*_{m}. In the case of the microscopic model, in which the optic axes within one fibre have different inclination angles, the upper limit of *δ* scales linearly with *t*_{m} as long as the optic axes of neighbouring myelin voxels have similar orientations (appendix C). Thus, the dependence on the myelin density can be eliminated to the greatest possible extent by multiplying the phase shift *δ* with a correction factor (*t*/*t*_{m}):
3.3In order to apply the myelin density correction to the downsampled retardation images, *t*_{m} was replaced by the combined thickness of myelin voxels along the optical path (after applying the Gaussian filter and resampling). The resulting retardation images were normalized according to equation (3.2), yielding .

## 4. Simulation results

### 4.1. Comparison of analytical and numerical solution

To estimate the accuracy of the simulation results for the microscopic model, a single fibre with radial optic axes and perpendicularly incident light (figure 3*b*) was generated for different voxel sizes (Δ*x*_{sim}) and the numerically computed phase difference between extraordinary and ordinary wave (Δ*Φ*_{num}) was compared to the analytical solution (Δ*Φ*_{ana}).

Assuming that reflection and refraction effects can be neglected so that the associated extraordinary and ordinary waves follow the same pathway, Bear & Schmidt [18] derived an analytical expression for the phase difference:
4.1where *Γ* is the optical path length difference between the extraordinary and ordinary waves, *r*_{1} the radius of the whole nerve fibre (outer cylinder), *r*_{2} the radius of the non-birefringent axon (inner cylinder) and *r*_{0} the distance at which the light is incident perpendicular to the fibre axis (figure 3*b*).

In order to compute Δ*Φ*_{num}, the propagation of the ordinary and extraordinary waves was simulated separately: in the case of the ordinary wave, the light is polarized parallel to the longitudinal axis of the fibre. In the case of the extraordinary wave, the light is polarized perpendicular to the longitudinal axis of the fibre (figure 3*b*). The phase for both the ordinary wave (*Φ*_{o}) and the extraordinary wave (*Φ*_{e}) was calculated from the corresponding electric field vector *E*_{T} in equation (3.1):
4.2

The numerically computed phase difference Δ*Φ*_{num} = *Φ*_{e} − *Φ*_{o} was evaluated at various distances 0 < *r*_{0} < 5 µm away from the centre of the fibre and compared to the analytical solution given in equation (4.1), with *r*_{1} = 7.5 µm, *r*_{2} = 5 µm, *λ* = 525 nm and Δ*n* = 0.001875 (cf. §3.3). In order to study the impact of the spatial discretization on the accuracy of the numerical solution, the simulation was performed for various voxel sizes 1.50 µm > Δ*x*_{sim} > 0.06 µm. As a measure of consistency between the numerical and the analytical solution, the relative phase difference was calculated: (Δ*Φ*_{ana} − Δ*Φ*_{num})/Δ*Φ*_{ana}.

Figure 4*a*,*b* shows the relative phase difference plotted against *r*_{0} for various voxel sizes Δ*x*_{sim}. As can be seen, the numerical solution fluctuates around the analytical solution for voxel sizes of 0.5 µm and less. With smaller voxel sizes, the numerical solution approaches the analytical solution (indicated by the dashed black line). This behaviour is especially evident when considering the mean of the absolute relative phase difference for each voxel size (figure 4*c*): for a voxel size of Δ*x*_{sim} = 1.5 µm (corresponding to one-tenth of the fibre diameter), the mean absolute relative phase difference is about 12%. For Δ*x*_{sim} = 0.5 µm, it is about 6% and for Δ*x*_{sim} = 0.06 µm, it is only 0.8%. This demonstrates that the simulation tool produces correct results.

As a good compromise between computation time and accuracy, all following fibre simulations were performed with a voxel size of Δ*x*_{sim} = 0.1 µm (corresponding to 1/150 of the fibre diameter). For this voxel size, the relative phase difference is no more than 4% (figure 4*b*) and the mean relative phase difference is about 1.3% (figure 4*c*).

### 4.2. Simulation of a single fibre

In a preliminary study, the limitations of the effective model of uniaxial negative birefringence were first studied for a straight single fibre. The fibre was simulated according to both the macroscopic and the microscopic models with different inclination angles (*α* = 0°, 10°, … , 90°) and different optical resolutions. The dimensions of the fibre and the other simulation parameters were chosen as described in §3.3. The retardation curves were calculated from the downsampled retardation images (without/with myelin density correction) and normalized as described in §3.5. An example of downsampled and corrected retardation images can be found in appendix D.

Figure 5 shows the dimensions of the simulated single fibre and the corresponding retardation curves (continuous lines) for both simulation models and different optical resolutions (according to table 1). The theoretical retardation curve of the effective model is indicated by a dashed black line. In the case of the macroscopic model, the uncorrected retardation curves (figure 5*a*) are already very similar to the theoretical curve for all investigated optical resolutions. After the myelin density correction (figure 5*c*), all retardation curves match the theoretical curve exactly, independently of the optical resolution. In the case of the microscopic model (figure 5*b*,*d*), the retardation curves for a pixel size much smaller than the fibre diameter (Δ*x* *<* 2 µm) are inverted as compared to the theoretical curve for *α* < 90°, i.e. the microscopic and the macroscopic model yield totally different results. For intermediate pixel sizes (2 µm ≤ Δ*x* ≤ 8 µm), the retardation curves are non-monotonic, i.e. the assignment of the inclination angle is ambiguous. Finally, for pixel sizes larger than the fibre diameter (Δ*x* = 16 µm), the uncorrected retardation curve (figure 5*b*) is similar to the theoretical curve. After the myelin density correction (figure 5*d*), the retardation curve matches the theoretical curve almost exactly.

### 4.3. Simulation of a fibre bundle

In brain tissue, nerve fibres are usually organized in hexagonal close-packed fibre bundles [13]. In order to investigate the effect of fibre bundles on the 3D-PLI signal, a hexagonal bundle of straight parallel fibres with an inter-fibre spacing of 1 *μ*m was simulated (figure 6*e*). In order to obtain comparable results, the same dimensions and simulation parameters were chosen as for the single fibre.

Figure 6 shows the normalized retardation curves for both simulation models and different optical resolutions (according to table 1). The (downsampled) retardation images that were used to compute the corrected retardation curves of the microscopic model are shown in appendix D.

In the case of the macroscopic model, the uncorrected retardation curves (figure 6*a*) are very similar to the theoretical retardation curve of the effective model (dashed black line) for all investigated optical resolutions. As compared to the retardation curves of the single fibre (figure 5*a*), the retardation curves of the fibre bundle are closer to the theoretical curve. After the myelin density correction (figure 6*c*), the curves are almost identical. In the case of the microscopic model, the uncorrected retardation curves (figure 6*b*) are also closer to the theoretical curve as compared to the uncorrected retardation curves of the single fibre (figure 5*b*). The myelin density correction (figure 6*d*) makes only a small difference, especially for low optical resolutions. For the simulated fibre bundle, the transition between the microscopic and the macroscopic model already occurs for pixel sizes larger than the fibre radius (Δ*x* ≥ 8 µm).

## 5. Discussion

In 3D-PLI, the fibre orientations are derived under the assumption that the brain tissue can (locally) be described as a homogeneous and uniaxial birefringent material with the optic axis indicating the predominant fibre direction. Furthermore, the density of myelinated fibres is assumed to be the same for the whole brain section. In this paper, the limitations of this effective birefringence model have been studied for the first time. For that purpose, a single fibre and a hexagonal fibre bundle (with diameters *d*) were simulated based on the Jones matrix calculus, employing a microscopic and a macroscopic model of birefringence and different optical resolutions (defined by the pixel size Δ*x* as given in table 1).

The transition between the two models is apparent when analysing the retardation curves: for high optical resolutions the radial optic axes of the microscopic model are resolved. In this case, the optic axes are oriented perpendicular to the longitudinal fibre axis so that the retardation curves are inverted as compared to the macroscopic model and fibres with high inclination angles are interpreted as flat fibres. The zero retardation value for *α* = 90° is an artefact arising from the fact that the retardation is evaluated at the centre of the retardation image which—in the case of vertical fibres—contains no myelin (cf. figure 8, upper right corner). For intermediate optical resolutions (Δ*x* *<* *d*), there is a transition zone between the microscopic and the macroscopic model so that an unambiguous assignment between retardation and inclination is not possible. For sufficiently low optical resolutions (single fibre: Δ*x* > *d*; fibre bundle: Δ*x* *>* *d*/2), the microscopic and the macroscopic model yield similar results (figures 5*d* and 6*d*) so that the effective model of uniaxial negative birefringence can be used to compute the fibre inclinations.

Thus, for the simulated fibre bundle (consisting of five fibre layers with *d* = 15 µm), the effective model can be used to interpret LAP measurements (Δ*x*_{LAP} = 64 µm > *d*/2), but not to interpret PM measurements (Δ*x*_{PM} = 1.33 µm < *d*/2). However, the diameters of the simulated fibres represent an upper estimate of typical fibre diameters in the human brain. The diameters of myelinated nerve fibres range from 0.3 to 15 µm [16,27–29] and the majority of the fibres (e.g. 80% in the corpus callosum [27]) have diameters of 1 µm or below so that the condition Δ*x*_{PM} > *d*/2 is still fulfilled. In addition, fibre diameters much smaller than 15 µm implicate that the measured brain section (with thickness 70 µm) contains much more fibre layers than the simulated fibre bundle. A comparison between the simulated single fibre and the fibre bundle suggests that the more fibre layers located along the optical path, the smaller the minimum pixel size for which the effective model is still valid. To verify this hypothesis, the limitations of the effective model should also be studied in terms of the number of fibre layers along the optical path. However, a larger number of fibre layers also increases the probability that fibres with different spatial orientations are measured within the same volume, which poses a major challenge for 3D-PLI [6]. In future studies, the limitations of the effective model should therefore also be investigated for non-parallel fibre structures.

The simulations have shown that—in regions with parallel fibre structures—the effective model of uniaxial negative birefringence is valid for the employed optical set-ups. For imaging systems with very high optical resolutions, the effective model needs to be reconsidered. Even if the optical resolution is too high to extract the correct fibre inclinations, 3D-PLI remains a valuable neuroimaging technique as the image contrasts of transmittance and retardation still provide detailed structural information on the two-dimensional nerve fibre architecture in large histological brain sections.

The effective model that is currently used for the data analysis in 3D-PLI does not only assume parallel fibre structures, but also a uniform myelin density. The simulations have shown that the retardation signal is considerably influenced by the myelin density, which impairs the reconstructed fibre orientations. It could be demonstrated that the estimation of the fibre inclination is considerably improved by the myelin density correction which incorporates the local myelin thickness of the examined tissue into the calculation of the inclination angle. While the correction has a large effect on the retardation curves of the single fibre, the effect is smaller for the fibre bundle which is much more homogeneous than the single fibre. Thus, the myelin density correction is especially useful for regions with an inhomogeneous density of myelinated nerve fibres (e.g. for transition zones between white and grey matter). In the case of the microscopic model, the correction does not work as well as for the macroscopic model because the retardation also depends on the direction of the radially oriented optic axes in the myelin sheath, but it is still a considerable improvement. In order to incorporate the myelin density correction into the 3D-PLI signal analysis, the local myelin thickness *t*_{m} of the sample needs to be determined. The intensity values of the transmittance image seem to be a good measure of the local myelin thickness in brain tissue [5].

The purpose of this study was to explore and understand the most dominant effects that generate the birefringence signals in 3D-PLI. To fully understand the physical processes behind 3D-PLI and to improve the interpretation of the reconstructed fibre orientations, a direct comparison between simulation and experiment is required. The long-term aim should be to develop a simulation tool of 3D-PLI that considers all relevant effects needed for reproducing the experimental results. To this end, the simulation model should be extended step by step and the relevant effects should be identified.

Although the simulations show that the simplified microscopic model can already be used to explain the effective negative birefringence of parallel nerve fibres, future studies should include the positive birefringence of the axon and investigate how this modification changes the transition between the microscopic and the macroscopic model.

So far, only straight and parallel fibres have been investigated. To provide more realistic fibre models, the fibres should be simulated with varying fibre diameters, myelin sheath thickness and spatial orientations. As fibres with different spatial orientations pose a major challenge for 3D-PLI [6], future studies should focus on investigating inhomogeneous, non-parallel fibre structures. To enable a direct comparison with the experiment, the simulated fibre configurations should be based on experimentally determined fibre structures.

In addition to a more realistic fibre model, the propagation of light should also be simulated more realistically. In this study, the incident light was described by a parallel beam of light. However, in the experiment, the employed light source emits diffusive light, i.e. the sample is illuminated by light with slightly different angles of incidence. As the measured birefringence signals depend on the angle between the light wave and the nerve fibres, a non-zero angle of incidence changes the retardation curves. For the LAP, which has a small numerical aperture, the effect can presumably be neglected. However, for systems with higher optical resolutions and higher numerical apertures, the effect of a Gaussian distribution of incident angles should be investigated further.

Moreover, the simulations were based on the Jones matrix calculus which is only applicable to completely polarized and coherent light. As the light source emits incoherent light and the polarizers are not perfect, the Jones matrices should be replaced by Müller matrices [31], which enable partially polarized and incoherent light to be studied.

Finally, the assumption of a linear optical pathway is a great simplification. The refractive index of the myelin sheath is higher than the refractive indices of the inner axon and the surrounding tissue [9,32,33], which will cause refraction/reflection at the interfaces and scattering of light. In future studies, the effects of refraction and scattering on the measured birefringence signal should be investigated in more detail. As the used simulation tool (SimPLI) is based on a matrix calculus, other simulation approaches will be required to investigate such nonlinear pathways.

## 6. Conclusion

In this study, we laid a theoretical foundation for 3D-PLI. The effective model of uniaxial negative birefringence, which is currently used to compute the nerve fibre orientations from experimental data, has been validated for the first time. Using simulations based on the Jones matrix calculus, we have shown that the effective model can be used for the employed optical set-ups, i.e. as long as the polarimeter does not resolve structures smaller than the diameter of single nerve fibres. The developed Jones matrix formalism for simulating 3D-PLI has proved to be a powerful tool to gain a deeper theoretical understanding of the physical processes behind 3D-PLI and to better interpret the experimental data. The simulations enable not only validation of the computational model of the fibre reconstruction, but also optimization of the experimental set-up and the measurement method.

## Authors' contributions

M.M. substantially contributed to the conception and design of the study as well as to the acquisition, analysis and interpretation of the simulated data. She carried out the simulations as well as the analytical calculations and drafted the manuscript. K.M. participated in the design of the study, contributed to theoretical considerations and to the interpretation of the simulated data, and helped draft the manuscript. H.D.R. contributed to the interpretation of the simulated data, to the theoretical considerations and to the revision of the manuscript. J.R. conducted experimental measurements, helped transfer the measurement results to the simulation and revised the manuscript. K.A. contributed to the anatomical content of the study and to the revision of the manuscript. M.A. coordinated the study, participated in the conception and design, contributed to the analysis and interpretation of the simulated data, and helped draft the manuscript. All authors read the final manuscript and gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by the Helmholtz Association portfolio theme ‘Supercomputing and Modeling for the Human Brain’ and by the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 604102 (Human Brain Project).

## Acknowledgements

The authors thank Melanie Dohmen for the introduction to the simulation software SimPLI.

## Appendix A. Derivation of the phase shift

When polarized light passes through the birefringent brain section, it is split into an ordinary and an extraordinary wave which both experience different refractive indices. The refractive index *n*_{e} that the extraordinary wave experiences when passing through the birefringent tissue under an angle *θ* with respect to the optic axis, is given by [34]:
A 1
A 2where *n*_{o} is the ordinary refractive index and the principal extraordinary refractive index of the brain tissue.

The birefringence of biological tissue ( [35]) is small as compared to the values of the refractive indices *n*_{o} and *n*_{E} (*n* = 1.3–1.5 [33]). Therefore, a Taylor expansion can be applied to the function
A 3in Δ*n* = 0:
A 4

The same expansion can be done for in . With these Taylor expansions, equation (A 2) can be written as A 5

Choosing a coordinate system in which the light propagates in the *z*-direction and the brain tissue lies in the *xy*-plane, the optic axis (oriented in the direction of the nerve fibres) makes an angle *θ* with the *z*-axis, i.e. the out-of-plane inclination angle of the fibre is *α* = 90°−*θ*. With this definition follows: .

Thus, when the light passes through a brain section of thickness *t*, the extraordinary wave experiences a phase shift with respect to the ordinary wave which depends on the inclination angle of the optic axis:
A 6This is the formula of the phase shift as given in equation (2.7).

## Appendix B. Derivation of the downsampling parameters

In previous measurements, the optical resolution of the LAP was investigated by employing a *USAF* test chart which contains line pairs (lp) with different spacings [36]. From the measured line intensity profiles, the Michelson contrast was computed as
B 1where *I*_{max} corresponds to the mean intensity of the maxima and *I*_{min} to the mean intensity of the (local) minima in the line intensity profile (cf. figure 7*b*). The largest number of line pairs per millimetre that can just be resolved (according to the Rayleigh criterion) was determined to be 5.66 lp mm^{−1}, which corresponds to a width per line pair of *l*_{LAP} = 176.7 µm and a contrast of . A width per line pair of 157.5 µm yields a Michelson contrast of 11.2%.

According to these measurement results, a test image with three lines (pixel size: 0.1 µm) and a line width of was created, and the downsampling procedure (Gaussian filter and resampling) was applied to the test image (figure 7*a*). The sampling factor was calculated by dividing the pixel size of the test image by the pixel size of the LAP: *f*_{s,LAP} = 0.1 µm/64 µm. To reproduce the measured contrast of the line intensity profile (figure 7*b*), a Gaussian filter with a standard deviation of *σ*_{PM} = 45.7 µm was applied. To avoid boundary effects and ensure a symmetric line intensity profile, the dimensions of the image (1216 × 1216 µm) were chosen such that the downsampled image consists of an odd number of pixels (19 × 19 px).

Based on the determined parameters for the LAP (Δ*x*_{LAP}, *σ*_{LAP}, *f*_{s,LAP}), downsampling parameters for imaging systems with other optical resolutions (table 1) were derived: the ratio between the pixel size of the LAP image and the determined standard deviation of the two-dimensional Gaussian filter is Analogous measurements of the PM yield a similar ratio between pixel size and standard deviation [36]. Assuming that this ratio is the same for all simulated imaging systems, the standard deviation of the two-simensional Gaussian filter was calculated from the pixel size Δ*x* of the resulting downsampled image:
B 2

After applying the Gaussian filter, the synthetic image series (with pixel size Δ*x*_{sim}) was resampled with a sampling factor of
B 3yielding a downsampled image series with pixel size Δ*x*.

## Appendix C. Dependence of the phase shift on the local myelin thickness

Each myelin voxel of a simulated nerve fibre is represented by the Jones matrix of a rotated wave retarder as defined in equation (2.2). Depending on what kind of model is used (macroscopic or microscopic), the retarder axis is either oriented parallel or radially to the fibre axis (figure 2*b*).

In the macroscopic model, the optic axes of the myelin voxels are all oriented in the fibre direction (*φ*, *α*) so that the voxels can be described by the same Jones matrix *M*_{δ}(*β*) with phase shift *δ* and . When the light propagates through *N* voxels of myelin, the multiplication of the *N* corresponding Jones matrices yields (using equation (2.2) and *R*(*β*)*R*(−*β*) = ):
C 1Thus, the *N* myelin voxels with thickness Δ*t* (along the optical path) and phase shift *δ* can be replaced by one myelin voxel with side length and phase shift:
C 2In other words, the phase shift *δ* (and for small *δ* also the retardation ) scales linearly with the combined thickness of myelin voxels (*N*Δ*t*), i.e. with the local myelin thickness *t*_{m}.

In the microscopic model, the optic axes of the myelin voxels along the optical path all have different orientations (*φ _{j}*,

*α*

_{j}) (figure 2

*c*). If the optic axes of neighbouring myelin voxels have a similar direction ( and ), the multiplication of the

*N*Jones matrices of the voxels can be simplified. For , one can define and the multiplication of a pair of rotation matrices yields: C 3using a first-order approximation in

*η*

_{21}.

The multiplication of two Jones matrices yields (with ): C 4

The multiplication of four Jones matrices yields (ignoring terms in the order of ): C 5where the elements of the secondary diagonal are given by C 6 C 7

If the number of myelin voxels (i.e. the number of matrices *M _{j}*) is small, the elements of the secondary diagonal in the resulting matrix can be neglected for . If the number of myelin voxels is large, the arguments of the exponential functions in the secondary diagonal will take all possible values and cancel each other for . In both cases, the multiplication of

*N*Jones matrices yields: C 8

Thus, the *N* myelin voxels with thickness Δ*t* and phase shift *δ _{j}* can be replaced by one myelin voxel with thickness (

*N*Δ

*t*=

*t*

_{m}) and phase shift: C 9given that the optic axes of neighbouring myelin voxels have similar directions.

The analytical considerations have shown that the phase shifts of individual voxels add together in both simulation models: in the macroscopic model, the phase shift scales linearly with the local myelin thickness *t*_{m}. In the microscopic model, this is only true for the upper limit of the phase shift. The dependence on the local myelin thickness is taken into account in the myelin density correction (see §3.5).

## Appendix D. Retardation images of the fibre bundle

- Received August 17, 2015.
- Accepted September 11, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.