Digital histology with Mueller microscopy: how to mitigate an impact of tissue cut thickness fluctuations.

Mueller microscopy studies of fixed unstained histological cuts of human skin models were combined with an analysis of experimental data within the framework of differential Mueller matrix (MM) formalism. A custom-built Mueller polarimetric microscope was used in transmission configuration for the optical measurements of skin tissue model adjacent cuts of various nominal thicknesses (5 to 30  μm). The maps of both depolarization and polarization parameters were calculated from the corresponding microscopic MM images by applying a logarithmic Mueller matrix decomposition (LMMD) pixelwise. The parameters derived from LMMD of measured tissue cuts and the intensity of transmitted light were used for an automated segmentation of microscopy images to delineate dermal and epidermal layers. The quadratic dependence of depolarization parameters and linear dependence of polarization parameters on thickness, as predicted by the theory, was confirmed in our measurements. These findings pave the way toward digital histology with polarized light by presenting the combination of optimal optical markers, which allows mitigating the impact of tissue cut thickness fluctuations and increases the contrast of polarimetric images for tissue diagnostics.

Introduction

Studies of biological tissue with polarized light may bring important information on the tissue’s polarimetric properties, namely, depolarization power, retardance, and dichroism. The evolution of these properties with a disease (e.g., inflammation, degeneration, cancer, etc.) suggests using them as optical markers for diagnostics in clinical settings. However, using these parameters for diagnostics requires (i) understanding the fundamental processes of interaction of polarized light with tissue and (ii) finding the optimal set of optical parameters, which will increase the accuracy of diagnostics.

Almost all biological tissues scatter incident light, so they are usually highly depolarizing. Consequently, one needs to use Stokes–Mueller formalism to describe the interaction of polarized light with tissue. Within its framework, both incident and reflected (or transmitted) polarized light beams are described by real Stokes vectors. The corresponding transfer function of a sample is called a Mueller matrix (MM; real matrix). The experimental systems that measure all elements of MM of a sample are called Mueller polarimeters. These instruments may operate in a spectroscopic or angular-resolved mode, and may also provide both microscopic (few hundreds of ) and macroscopic (few ) polarimetric images of the sample.,,,

MM contains all information on polarimetric and depolarizing properties of a sample. However, a straightforward physical interpretation of MM elements is possible for a quite limited set of samples like basic polarimetric elements (polarizers and wave plates) and partial depolarizers. The most widely used method of nonlinear data reduction for the interpretation of general MM was proposed by Lu and Chipman. Despite widespread use of Lu–Chipman polar decomposition, it has some drawbacks when applied for the interpretation of MMs of biological samples. Lu–Chipman decomposition assumes a sequential appearance of basic optical effects (dichroism, retardance, and depolarization) along the path of the light beam. Obviously, this assumption does not hold for biological tissues, where all effects may appear simultaneously. It was shown that in transmission configuration, the most appropriate decomposition of MM, namely, logarithmic Mueller matrix decomposition (LMMD), is described within the framework of the differential formalism of fluctuating anisotropic media.,

Our prior studies of isotropic and anisotropic scattering phantoms, demonstrated the validity of LMMD in transmission configuration. Therefore, in this paper, we extend this approach to biological tissue models and report on the results of our studies of full-thickness skin equivalents with transmission Mueller microscopy and LMMD, discussing their potential diagnostics value. The algorithm of mitigating the impact of tissue thickness variations on polarimetric parameters of histological cuts is suggested. It is based on the theoretical predictions of Mueller differential formalism and Beer–Lambert law.

Materials and Methods

Tissue Samples and Experimental Set-Up

Human skin equivalents were produced in vitro from human cells and reflected the anatomy of human skin. The real skin can roughly be divided into three parts: epidermis, dermis, and subcutis (mostly fatty tissue, not included in our skin model). The skin models were generated from primary human skin cells (keratinocytes and fibroblasts)., The former cells differentiate in vitro and form an epidermis with the same anatomical layers as in vivo: stratum basale, stratum spinosum, stratum granulosum, and stratum corneum. The dermal part of the skin model consists of a collagen type 1 hydrogel with human primary fibroblasts. The real dermis can be divided into an uppermost part (stratum papillare) and a lower part (stratum reticulare). Stratum papillare of the dermis was not recreated in the skin model since it serves only as the mechanical interlocking of the epidermis and dermis. However, the typical cell sizes and shapes in this skin model are the same as the ones in real in vivo human skin.

These similarities and other functional properties, such as transporter expression, barrier function, etc., led to the use of such skin models as alternatives to animal models or human donor tissue. This is one of the reasons these models achieved regulatory acceptance by validation and adoption in the Organization for Economic Cooperation and Development guidelines for regulatory toxicological tests, e.g., skin irritation/corrosion (OECD TG 439). This means that these models are employed in Europe and other OECD countries to categorize substances for their potential to cause skin irritation and corrosion. Since human skin equivalents can be produced with less variability compared to that of human skin, these skin tissue models were chosen for our studies.

The models were grown in so-called cell culture inserts (Corning™ Snapwell™) with a diameter of 12 mm. The epidermis thickness was about , and the thickness of the dermal part was close to . Therefore, we obtained tissue disks of 12-mm diameter and height of . The grown tissue models were rinsed with phosphate-buffered salt solution and fixed with Roti®-Histofix 4% for 4 h at room temperature. Then, fixed samples were embedded in paraffin in an embedding machine. First, a disk of paraffin-embedded skin tissue model was cut along the diameter [see Fig. 1(a)]. Then, a set of adjacent histological cuts of different thickness (5, 10, 16, 20, and ) of around 1 cm in length and 0.5 mm in height were prepared from both parts of the disk using a microtome [see Figs. 1(b) and 1(c)]. Thereafter, the samples were deparaffinized for 20 min in Roticlear® and placed on a microscope glass slide [Fig. 1(d)]. There was no coverslip used in our studies.

Fig. 1

(a)–(e) Schematics of a preparation of the adjacent histological cuts of varying thickness; (e) photo of five microscope glass slides with two identical unstained cuts () of skin tissue model (see text) mounted on each slide. The nominal thicknesses of histological cuts are stated on paper labels.

The custom-built Mueller polarimetric liquid-crystal-based microscope operating in a visible wavelength range was used for the measurements of MM of thin samples in transmission configuration. The illumination arm of the setup consists of a white-light LED source, a set of lenses, and two diaphragms for independent control of beam divergence and size followed by a polarization state generator (PSG). The PSG is composed of a linear polarizer and two ferroelectric-liquid-crystal retarders (Meadowlark FPR-200-1550). A collimated beam of diameter and fixed polarization uniformly illuminate a sample. The transmitted light passed through an imaging lens (Thorlabs AC254-030-A-ML) followed by a polarization state analyzer (PSA) and a CCD camera (AV Stingray F-080B) coupled to a telephoto lens, which is adjusted to infinity. The sample is placed on the principal object plane; thus, a real space image of a sample is formed on the CCD detector. The arrangement of optical elements of PSA is reciprocal to that of a PSG. The wavelength of 533 nm was selected for our measurements by placing an interferential filter (spectral bandwidth of 20 nm) after a PSA. The measurements of histological cuts of skin models were performed with a objective with a field-of-view (FoV) of around . More details on the experimental setup can be found in Ref. 10.

A CCD camera was placed directly in a transmitted beam path. When switching from sample-to-sample, the overall signal registered by CCD showed significant variations. The transmitted intensity was higher for thinner samples compared to thicker ones, in accordance with Beer–Lambert law. To avoid the saturation problem, we used the measurement protocol described below. Due to the technical characteristics of the CCD detector, the polarimetric measurements were performed within a given intensity range to ensure the linearity of CCD response. Therefore, the integration time of a CCD was adjusted for every sample to get a well-balanced signal level for all 16 images needed to measure the corresponding MM. This procedure helped us avoid both over- and under-exposure. It is worth noting that while all histological cuts were relatively thin and transmitted a significant fraction of the direct light, the scattering of light produced noticeable effects in depolarization properties due to the incoherent summation of direct and scattered light signal on CCD.

Logarithmic Decomposition of Mueller Matrix

Different algorithms of decomposition of MM have been extensively studied, and several methods (e.g., Lu–Chipman, reverse, symmetrical and differential) were proposed for the MM data interpretation. Among them, a logarithmic decomposition method developed for transmission geometry is the one that considers all optical properties as continuously distributed within the volume of medium. It makes LMMD particularly suitable for the studies of biological tissue in a transmission configuration. We briefly summarize the key steps of LMMD below. Within the framework of differential matrix formalism of a fluctuating anisotropic medium, the transmission MM is described by the following equation:

The MM , which is dependent on optical path length , is associated with a unique differential matrix . This matrix is constant for both nondepolarizing and depolarizing media, which are homogeneous along the light propagation direction. For a depolarizing medium, the differential matrix can be decomposed into G-antisymmetric and G-symmetric [where is the Minkowski metric and T denotes matrix transposition]:

The elements of matrix ( through ) represent the elementary polarization properties—linear (, ) and circular dichroism, linear (, ) and circular retardance (recall that the dichroic and retardance elementary properties are proportional, respectively, to the imaginary and real parts of the linear and circular anisotropies, i.e., of the differences of the linear and circular complex refractive indices of the medium,). The elements of matrix ( through ) describe the depolarization properties of the medium. Diagonal terms ( through ) represent the anisotropic depolarization coefficients, and the elements ( through ) show the uncertainties of polarization properties.

The statistical meaning of coefficients of MM of a continuous depolarizing medium implies that the depolarization is a result of a spatial or temporal averaging process over when the polarization properties of medium (contained in differential matrix ) fluctuate and matrix varies. In such a case, it was demonstrated that the matrix represents mean values of the polarization properties. The matrix contains mean square values of fluctuations of polarization properties, i.e., their variances (or uncertainties) and linearly depends on slab thickness [see Eq. (3), brackets refer to the spatial averaging in the transverse plane]. If the medium is assumed to be homogeneous in the longitudinal direction of light propagation, then

Substituting the statistical representation of differential matrix from Eq. (3) into Eq. (1) and integrating the latter equation along , we obtain

It follows from Eq. (4) that mean values of polarimetric properties scale up linearly with thickness while the depolarization properties evolve quadratically with thickness. The differential matrix of a homogeneous medium can be obtained from a simulated or experimentally measured MM of a sample by computing the matrix logarithm, which can be represented as a sum of two matrices and of opposite G-symmetry:

Calculating the logarithm of Eq. (4) at (i.e., taking the thickness of the slab as unit one), we observe that the antisymmetric component and the symmetric component , respectively, equal the mean values and (half) the variances of the polarization properties, accumulated over the slab thickness:

It is worth to recall that each element of and matrices has a straightforward physical interpretation in terms of polarimetric properties of a sample., The experimental validation of logarithmic decomposition of MM on biological tissue samples, namely, a linear dependence of polarization properties with thickness and quadratic dependence of depolarization properties will be described in the next sections.

Results and Discussion

Logarithmic Decomposition of Mueller Matrix Images

The set of five skin model histological cuts of varying thickness was measured with a Mueller microscope. All experimental MM images of histological cuts were processed pixelwise by applying LMMD. The corresponding maps of total linear retardance and dimensionless diagonal coefficient of matrix , which demonstrates the linear depolarization property, are shown in Fig. 2.

Fig. 2

Maps of (a) total linear retardance (in rad) and (b) dimensionless diagonal coefficient of matrix , calculated from the experimental MM data by applying pixelwise LMMD. Labels represent the mean thickness of histological cuts measured with a stylus profilometer. FoV is about .

Three different zones are clearly distinguishable on both maps of total linear retardance and linear depolarization coefficient : (1) bare glass without tissue, (2) dermal, and (3) epidermal zones of human skin model. Microscope glass does not possess any measurable linear retardance () as well as does not depolarize transmitted light (). The dermal part of a skin model histological cut demonstrates strong retardance. This effect is related to tissue anisotropy due to the presence of the aligned collagen fibers. On all images (see Fig. 2), there is a thin layer of epidermis underneath of dermis. The former layer does not show any retardance for all values of thickness of histological cuts. This is an expected result because the epidermal layer does not contain any aligned collagen fibers. The values of were always lower within the zone of epidermal layer compared to the zone of dermal layer, which means that the epidermis is more depolarizing compared to the dermis. We attribute this fact to stronger scattering of light by cells in the epidermis compared to scattering of light on aligned collagen fibers and scarce fibroblasts in a dermal layer. As expected, there was no circular retardance and circular dichroism observed for all tissue cuts [coefficients , Eq. (2)]. It can be explained by the fact that optical activity in biological tissues is related to the presence of chiral molecules (e.g., glucose), which were absent in the studied skin model tissue. We have focused on the analysis of dermal layer of skin model cuts, because this zone of tissue possesses both polarization and depolarization properties contrary to epidermis layer, which depolarizes light only. Thus, to verify the predictions of differential formalism on thickness dependence of polarization and depolarization properties of fluctuating anisotropic media, we need to estimate the mean values of polarization and depolarization properties of dermal layer.

Image Segmentation

To delineate epidermal and dermal zones of skin tissue model on microscopic images, we have applied the MATLAB subroutine density-based spatial clustering of applications with noise (DBSCAN), which makes use of two control parameters for data clustering, namely, radius and number of neighbors “MinPts.” A dataset of points is defined in a parametric space. For an arbitrary initial data point, the algorithm finds neighbors within a circle with radius and assigns them to the same cluster. For any neighbor point P having a predefined number MinPts (or more) of data points within a circle with radius and center P, the cluster is expanded by adding those points as well. If the number of points in the neighborhood of data point P is less than the threshold value MinPts, the point is considered to be a noise. Contrary to K-means clustering algorithm, DBSCAN method does not require one to fix the number of clusters. It makes DBSCAN one of the most commonly used and cited clustering algorithms well adapted for the image segmentation problem.

The results of clustering depend on the input information we use. In this study, the values of —sample transmittance, —total linear retardance, and —dimensionless diagonal element of matrix , which represents circular depolarization at each pixel of an image, were used an input information. Furthermore, each parameter value was standardized using z-score [ , , where angle brackets denote mean value, stands for a standard deviation] for preventing one parameter be dominant in data clustering. First, we arranged blocks of neighboring pixels in a “super pixel” and defined its value by a bilinear interpolation (an output pixel value is a weighted average of pixels in the nearest neighborhood). Thus, we reduced the number of pixels to decrease the computational burden for subsequent segmentation. Then, we run the segmentation with the parameter and radius and obtained three well separated zones: (1) bare glass, (2) dermal layer, and (3) epidermal layer and random outliers (or noise part). The results of the segmentation for a histological cut of nominal thickness are presented in Fig. 3. Black markers represent the noise part, and blue, red, and green markers correspond to bare glass, dermal, and epidermal layers of the tissue, respectively. One can notice the presence of thin zone rendered in green above the dermal layer in Fig. 3(b). Despite being classified as epidermis, it makes a part of dermis. Most probably, the edge part of dermal layer has a different thickness because of cutting artifacts, which, in turn, alters all optical parameters used for the image segmentation.

Fig. 3

(a) Clustering results in XYZ -score space (see text), – retardance, – circular depolarization, – total transmitted intensity. Black markers: noise, blue: bare glass, red: dermis, green: epidermis and (b) corresponding image segmentation: noise, bare glass, dermis, and epidermis zones are rendered in black, blue, red, and green, respectively.

After selecting the group of pixels corresponding to the dermal layer of skin model tissue cuts, the mean values of polarization and depolarization properties were calculated over the pixels of dermal layer for all tissue cuts. It is worth to mention that clustering results are almost the same if we chose as input information the values of parameter or representing linear depolarization in the framework axes and axes, respectively, instead of parameter representing circular depolarization.

Thickness Dependence of Polarization and Depolarization Properties

The thickness of the sample should be known for a correct assessment of the dependence of polarimetric properties on thickness. We have used a Stylus Profilometer (Bruker DektakXT) to measure the thickness of tissue cuts and check if it matches the nominal values of thickness (5 to ). The number of depth scans for a generation of a 2-D image was set to 10, and the width of the scanning area was fixed at (close to the FoV of Mueller microscope). The resulting 2-D depth profile provides information on homogeneity and uniformity of sample thickness (Fig. 4).

Fig. 4

(a)  microscopic image of histological cut of tissue, nominal thickness: and (b) measured 2-D profile of tissue cuts, nominal thickness: .

The mean values of thickness averaged over 10 profilometer scans for the histological cuts of skin tissue models are presented in Table 1. Mean values differ significantly from the nominal ones, and this difference becomes larger for thicker samples. The slices of tissue-containing paraffin blocks were cut by the microtome with micrometer-controlled precision. The analysis of data from Table 1 suggests that the deparaffinization of tissue slides induces significant variations in the thickness of tissue cuts.

Table 1

Average thickness of tissue cuts measured with a stylus profilometer.

Nominal thickness (μm)	Mean measured thickness (μm)	Standard deviation (μm)
5	3.0	0.8
10	6.2	1.1
16	7.6	2.3
20	10.1	4.2
30	10.5	2.9

To verify the dependence of polarization and depolarization properties of tissue cuts on thickness, we plotted the averaged values of polarization and depolarization parameters of dermal layer versus thickness of histological cuts measured with stylus profilometer (see Fig. 5).

Fig. 5

Thickness dependence plots of (a) total linear retardance (in rad), (b) total linear dichroism, (c)–(e) dimensionless depolarization coefficients , , and averaged over a dermal layer of the histological cuts. Experimental data are shown by black markers, solid red lines represent linear (a, b) and quadratic (c)–(e) fit curves.

As predicted by the theory [Eqs. (2) and (4)], a total linear retardance ( and are linear retardance along axes and linear retardance along axes, respectively) and total linear dichroism depend linearly on thickness [see Figs. 5(a) and 5(b)]. The presence of linear dichroism can be explained by the scattering on nonspherical scatterers like elongated collagen fibers. While intercept of linear regression curve with axis for is equal to zero, this is not the case for linear dichroism linear regression curve. We attribute this effect to the scattering of transmitted light on a rough surface of tissue. Surface scattering of an anisotropic medium does not affect the retardance values but also contributes to the increase of values of linear dichroism. The values of anisotropic depolarization coefficients , , and vary quadratically with thickness [see Figs. 5(c)–5(e)]. The absolute values of are larger compared to absolute values of and for the same tissue thickness. It means that the linear polarization of incident light is preserved better compared to the circular one. This is an indication of Rayleigh scattering regime.

Mitigating the Impact of Tissue Cut Thickness Fluctuations

The pathological changes of tissue (cancer, fibrosis, inflammation, etc.) will affect both polarization and depolarization properties of tissue. An ultimate goal of digital pathology consists of delineating the abnormal zones of a microscope image of histological cut using the maps of optimal optical markers providing the highest contrast. As we have shown above, both polarization and depolarization parameters of anisotropic scattering media vary with tissue thickness because of changing optical path length. Therefore, controlling the thickness of histological cuts is one of the crucial issues for accurate diagnostics. However, in practice, it is impossible to measure the real depth profile of histological cuts with profilometer as we did in these studies, because for a standard histology analysis, the tissue cuts are mounted on a microscope glass slide and protected by a coverslip (i.e., the tissue is “sandwiched” between two glasses).

We explore several approaches to eliminate the impact of local variations of tissue thickness on its measured polarization and depolarization properties. During the calibration of Mueller microscope, a bare glass was used as the reference sample. Since MM of a bare glass was included in the calibration data, element of MM represents a transmittance of tissue sample (without the glass). It follows from the Beer–Lambert law: where and are the intensities of input and output light beam, respectively, is a sum of absorption coefficient and scattering coefficient of the medium, is the thickness of a sample. The total fluence for intensity of the input light beam was controlled by the exposure time. During the measurements of tissue cuts of different thickness, the exposure time was varied to prevent the saturation of the detected signal, as was described in Sec. 2.1. That is why, values for different tissue cuts were rescaled to match an exposure time (250 ms) used in the calibration process. Finally, applying Eq. (7) pixelwise to image, one can produce a microscopic image of the optical density of studied tissue cut.

We assume that all skin model tissue cuts are homogeneous along with the incident light beam path (few microns scale), but tissue properties may vary over the imaged plane (FoV few hundreds of microns). Because of the linear dependence of retardance and quadratic dependence of depolarization on thickness, the following relations hold for each pixel of a microscopic image of histological cuts: where is the thickness of tissue cut, and and are linear and quadratic coefficients for a pixel . Consequently, the following quantities should not depend on the local thickness of tissue cut:

Equations (9)–(12) are invariant under tissue thickness fluctuations. Using these equations, we have calculated the thickness invariant microscopic maps for all histological cuts of different nominal thicknesses (see Fig. 6). While the values of , , and may still vary across the microscopic image, these variations are now related to the variations in tissue properties, not in tissue thickness.

Fig. 6

Maps of skin tissue model histological cut ( nominal thickness) calculated using Eqs. (9)–(12): (a) , (b) , (c) , and (d) . FoV is about . All these ratios are invariant under tissue thickness fluctuations.

We assume that the distributions of ratios , , and values, which are thickness invariant, should become more peaked compared to the distributions of both and values, which depend on both thickness fluctuations and fluctuations of tissue properties for a dermal layer of skin model tissue cuts. To check this assumption, we performed the statistical analysis of these distributions. We used the value of entropy , where is discrete random variable, set , is probability distribution function, as an “inversed” metric of distribution peakedness. Indeed, more peaked distributions are less undetermined; hence, their entropy should be lower compared to broader distributions. The calculated values of entropy are presented in Table 2.

Table 2

Entropy (in nats) of the distributions of parameters defined by Eqs. (9)–(12) for a dermal layer.

Real thickness (μm)	H(R_T)	H(α_44)	H(R_Tln(M_11))	H(α_44ln^2(M_11))	H(RT2α_44)	H(α_44RT2)
3.0	3.77	3.81	3.31	2.89	3.25	2.94
6.2	3.71	3.87	3.74	3.66	3.59	3.15
7.6	3.47	3.39	3.71	3.55	3.34	2.91
10.1	3.40	2.68	3.43	3.37	3.72	2.82
10.5	3.42	2.51	3.52	3.62	3.56	2.57

For a dermal layer of skin tissue model cuts, the smallest values of entropy correspond to the distribution of the parameter compared to other parameters calculated from Eqs. (9)–(12). The entropy values of distributions also drop for thick () dermis layers, but the values of depend both on tissue cut thickness and tissue polarimetric properties. We demonstrated that using the dependence of polarization and depolarization parameters on thickness, predicted by the differential MM formalism, one can produce the microscopic images with less fluctuations and higher contrast.

Conclusions

The experimental studies of histological skin tissue model cuts with the polarimetric Mueller microscope have confirmed the validity of phenomenological differential formalism of fluctuating anisotropic media for biological tissues. As per theoretical LMMD prediction, we have demonstrated that total linear retardance and total linear dichroism of dermal layer depend linearly on thickness while the depolarization parameters demonstrate quadratic dependence with thickness. The set of optical parameters including the circular depolarization and total linear birefringence (both derived from the logarithmic decomposition of MM of skin tissue model cuts) and the intensity of transmitted light (element ) was effectively used for the automated segmentation of microscopy images and delineation of the zones of bare glass, dermis, and epidermis.

An important issue, overlooked by many researchers working in the field of polarized light histology, appears to be the control and characterization of the real thickness of studied tissue cuts. The pathological changes of tissue (cancer, fibrosis, inflammation, etc.) will affect measured polarization and depolarization properties of a sample. However, changing thickness of tissue cut and, consequently, the optical path length will also affect these properties. Thus, for separating the contribution of both factors and reliable diagnostics of tissue with polarized light, the impact of the varying optical path length on polarization and depolarization optical markers of the specific disease has to be taken into account. We have proposed several approaches on using the linear and quadratic dependence on retardance and depolarization on thickness, respectively, combined with Beer–Lambert law for mitigating the impact of tissue thickness fluctuations and increasing the contrast of polarimetric images relevant for diagnostic purposes. Mueller microscopy studies of different types of tissue with pathologies will be the subject of our future work.