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.
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.1–5 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).6 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 mode7,8 and may also provide both microscopic (few hundreds of ) and macroscopic (few ) polarimetric images of the sample.2–4,5,9,10
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.11 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.12 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.13,14
Our prior studies of isotropic and anisotropic scattering phantoms10,15 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).16,17 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 43918). 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.
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]:19
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 medium13,19). 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 demonstrated19 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, then19
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.13,19 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.
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.
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.20
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.
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).
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.
Average thickness of tissue cuts measured with a stylus profilometer.
|Nominal thickness (μm)||Mean measured thickness (μm)||Standard deviation (μm)|
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).
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)].21 The presence of linear dichroism can be explained by the scattering on nonspherical scatterers like elongated collagen fibers.22 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.23 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: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:
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.
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,24 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.
Entropy (in nats) of the distributions of parameters defined by Eqs. (9)–(12) for a dermal layer.
|Real thickness (μm)||H(RT)||H(α44)||H(RTln(M11))||H(α44ln2(M11))||H(RT2α44)||H(α44RT2)|
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.
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.
The authors have no relevant financial interests in the manuscript and no other potential conflicts of interest to disclose.
H.R.L. and T.S.H.Y. gratefully acknowledge the funding from the doctoral school “Interfaces” of École Polytechnique, France. P.L. acknowledges the funding of his internship at LPICM, École Polytechnique, France, by the Chinese Scholarship Council.
Hee Ryung Lee is a PhD student at Applied Optics and Polarimetry Group of the Laboratory of Physics of Interfaces and Thin Films, École Polytechnique, France. He received his MS degree in physics from Kyung Hee University (Seoul, Korea) and École Polytechnique (Palaiseau, France) in 2017. His current research interests include biophotonics, Mueller polarimetry, and polarized Monte Carlo modeling. He is a member of the SPIE student chapter at École Polytechnique, France.
Pengcheng Li is a PhD candidate at Laboratory of Optical Imaging and Sensing, Department of Physics, Tsinghua University, China. His research topics are Monte-Carlo simulation of polarized photon scattering in turbid media, Mueller matrix transformation theory (orientation invariant parameters of Mueller matrix), and its application in distinguishing different sources of anisotropies.
Thomas Sang Hyuk Yoo received his PhD degree in optics in 2018 from the École Polytechnique, France. He received his MS degree in engineering from Kyung Hee University (Seoul, Korea) and École Polytechnique (Palaiseau, France) in 2015. His current research interests include optical metrology and image processing.
Enric Garcia-Caurel received his PhD in physics from the University of Barcelona, Spain, in 2001. Since 2003, he has been a researcher at the Laboratory of Physics of Interfaces and Thin Films at the École Polytechnique (France), and has been an associated researcher at Synchrotron SOLEIL (France) since 2015. He is the author of four patents and more than 90 papers in peer-reviewed journals and international conferences. His current research interests include ellipsometry, polarimetry, materials science, and biophotonics. He is a SPIE senior member and advisor of the OSA–SPIE student chapter at the École Polytechnique.
Tatiana Novikova is a head of the Applied Optics and Polarimetry Group at the Laboratory of Physics of Interfaces and Thin Films, École Polytechnique, France. She received her PhD degree in applied mathematics from the Institute of Mathematical Modeling (Moscow, Russia) in 1999 and HDR (habilitation) in physics from the University Paris XI (Orsay, France) in 2015. Her research activities are focused on Mueller polarimetry and its applications for biomedical diagnostics and metrology in microelectronics, computational optics of polarized light, interaction of electromagnetic waves with matter, and physics of weakly ionized plasmas. She is an OSA senior member and SPIE member.