2.1.

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 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 $100\mu \mathrm{m}$, and the thickness of the dermal part was close to $500\mu \mathrm{m}$. Therefore, we obtained tissue disks of 12-mm diameter and height of $600\mu \mathrm{m}$. 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 $30\mu \mathrm{m}$) 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 ($\sim 1\mathrm{cm}\times 0.5\mathrm{mm}$) 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 $\sim 1\mathrm{cm}$ 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 $20\times$ objective with a field-of-view (FoV) of around $600\mu \mathrm{m}$. 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.

2.2.

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 $\mathbf{M}(z)$, which is dependent on optical path length $z$, is associated with a unique differential matrix $\mathbf{m}$. 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 $\mathbf{m}$ can be decomposed into G-antisymmetric ${\mathbf{m}}_{\mathrm{m}}$ and G-symmetric ${\mathbf{m}}_{\mathrm{u}}$ [where $\mathbf{G}=\mathrm{diag}(1,-1,-1,-1)$ is the Minkowski metric and T denotes matrix transposition]:¹⁹

The elements of matrix ${\mathbf{m}}_{\mathbf{m}}$ ($p_1$ through $p_6$) represent the elementary polarization properties—linear ($x\text{-}y$, $-45\mathrm{deg}$ $+45\mathrm{deg}$) and circular dichroism, linear ($x\text{-}y$, $-45\mathrm{deg}$ $+45\mathrm{deg}$) 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^13,19). The elements of ${\mathbf{m}}_{\mathrm{u}}$ matrix ($d_0$ through $d_9$) describe the depolarization properties of the medium. Diagonal terms ($d_7$ through $d_9$) represent the anisotropic depolarization coefficients, and the elements ($d_1$ through $d_6$) show the uncertainties of polarization properties.

The statistical meaning of coefficients of MM $\mathbf{M}$ of a continuous depolarizing medium implies that the depolarization is a result of a spatial or temporal averaging process over $\mathbf{M}$ when the polarization properties of medium (contained in differential matrix $\mathbf{m}$) fluctuate and matrix $\mathbf{M}$ varies. In such a case, it was demonstrated¹⁹ that the matrix ${\mathbf{m}}_{\mathrm{m}}$ represents mean values $⟨\mathbf{m}⟩$ of the polarization properties. The matrix ${\mathbf{m}}_{\mathrm{u}}$ contains mean square values of fluctuations of polarization properties, i.e., their variances (or uncertainties) $⟨\mathrm{\Delta }{\mathbf{m}}^2⟩$ and linearly depends on slab thickness $z$ [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 $\mathbf{m}$ from Eq. (3) into Eq. (1) and integrating the latter equation along $z$, 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 $\mathbf{m}$ of a homogeneous medium can be obtained from a simulated or experimentally measured MM $\mathbf{M}$ of a sample by computing the matrix logarithm, which can be represented as a sum of two matrices ${\mathbf{L}}_{\mathbf{m}}$ and ${\mathbf{L}}_{\mathbf{u}}$ of opposite G-symmetry:

Calculating the logarithm of Eq. (4) at $z=1$ (i.e., taking the thickness of the slab as unit one), we observe that the antisymmetric component ${\mathbf{L}}_{\mathbf{m}}$ and the symmetric component ${\mathbf{L}}_{\mathbf{u}}$, 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 ${\mathbf{L}}_{\mathbf{m}}$ and ${\mathbf{L}}_{\mathbf{u}}$ 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.

3.1.

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 $R_T=\sqrt{p_4^2+p_5^2}$ and dimensionless diagonal coefficient ${\alpha }_{22}$ of matrix ${\mathbf{L}}_{\mathbf{u}}$, which demonstrates the linear depolarization property, are shown in Fig. 2.

Fig. 2

Maps of (a) total linear retardance $R_T$ (in rad) and (b) dimensionless diagonal coefficient ${\alpha }_{22}$ of matrix ${\mathbf{L}}_{\mathbf{u}}$, 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 $600\mu \mathrm{m}$.

Three different zones are clearly distinguishable on both maps of total linear retardance $R_T$ and linear depolarization coefficient ${\alpha }_{22}$: (1) bare glass without tissue, (2) dermal, and (3) epidermal zones of human skin model. Microscope glass does not possess any measurable linear retardance ($R_T=0$) as well as does not depolarize transmitted light (${\alpha }_{22}=0$). 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 ${\alpha }_{22}$ 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 $p_3=p_6=0$, 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.

3.2.

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 $\epsilon$ 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 $\epsilon$ 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 $\epsilon$ 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 $M_{11}$—sample transmittance, $R_T$—total linear retardance, and ${\alpha }_{44}$—dimensionless diagonal element of matrix ${\mathbf{L}}_{\mathbf{u}}$, which represents circular depolarization at each pixel of an image, were used an input information. Furthermore, each parameter value was standardized using z-score [$X=(R_T-⟨R_T⟩)/{\sigma }_{R_T}$ $Y=({\alpha }_{44}-⟨{\alpha }_{44}⟩)/{\sigma }_{{\alpha }_{44}}$, $Z=(M_{11}-⟨M_{11}⟩)/{\sigma }_{M_{11}}$, where angle brackets denote mean value, $\sigma$ stands for a standard deviation] for preventing one parameter be dominant in data clustering. First, we arranged $2\times 2$ 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 $2\times 2$ neighborhood). Thus, we reduced the number of pixels to decrease the computational burden for subsequent segmentation. Then, we run the segmentation with the parameter $\mathrm{MinPts}=300$ and radius $\epsilon =0.2$ 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 $10\text{-}\mu \mathrm{m}$ 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 $Z$-score space (see text), $X$ – retardance, $Y$ – circular depolarization, $Z$ – 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 ${\alpha }_{22}$ or ${\alpha }_{33}$ representing linear depolarization in the framework $x\text{-}y$ axes and $\pm 45\mathrm{deg}$ axes, respectively, instead of parameter ${\alpha }_{44}$ representing circular depolarization.

3.3.

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 $30\mu \mathrm{m}$). 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 $500\mu \mathrm{m}$ (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) $20\times$ microscopic image of histological cut of tissue, nominal thickness: $16\mu \mathrm{m}$ and (b) measured 2-D profile of tissue cuts, nominal thickness: $16\mu \mathrm{m}$.

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.

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 ${\alpha }_{22}$, ${\alpha }_{33}$, and ${\alpha }_{44}$ 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 $R_T=\sqrt{p_4^2+p_5^2}$ (${\mathrm{p}}_4$ and ${\mathrm{p}}_5$ are linear retardance along $x\text{-}y$ axes and linear retardance along $\pm 45\mathrm{deg}$ axes, respectively) and total linear dichroism $D_T=\sqrt{p_1^2+p_2^2}$ 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 $Y$ axis for $R_T$ 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 ${\alpha }_{22}$, ${\alpha }_{33}$, and ${\alpha }_{44}$ vary quadratically with thickness [see Figs. 5(c)–5(e)]. The absolute values of ${\alpha }_{44}$ are larger compared to absolute values of ${\alpha }_{22}$ and ${\alpha }_{33}$ 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.

3.4.

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, $M_{11}$ element of MM represents a transmittance $I/I_0$ of tissue sample (without the glass). It follows from the Beer–Lambert law:

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 $(k,l)$ 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 ${\mu }_T^{k,l}$, $A^{k,l}$, and $B^{k,l}$ 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 ($5\text{-}\mu \mathrm{m}$ nominal thickness) calculated using Eqs. (9)–(12): (a) $R_T/\ln M_{11}$, (b) ${\alpha }_{44}/{\ln }^2(M_{11})$, (c) $R_T^2/{\alpha }_{44}$, and (d) ${\alpha }_{44}/R_T^2$. FoV is about $600\mu \mathrm{m}$. All these ratios are invariant under tissue thickness fluctuations.

We assume that the distributions of ratios $R_T/\ln (M_{11})$, ${\alpha }_{44}/{\ln }^2(M_{11})$, $R_T^2/{\alpha }_{44}$ and ${\alpha }_{44}/R_T^2$ values, which are thickness invariant, should become more peaked compared to the distributions of both $R_T$ and ${\alpha }_{44}$ 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 $H(X)=-\sum _{x_i\in C}p(x_i)(\ln p(x_i))$, where $X$ is discrete random variable, $X\in$ set $C$, $p(X)$ 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.

For a dermal layer of skin tissue model cuts, the smallest values of entropy correspond to the distribution of the parameter ${\alpha }_{44}/R_T^2$ compared to other parameters calculated from Eqs. (9)–(12). The entropy values of ${\alpha }_{44}$ distributions also drop for thick ($>10\mu \mathrm{m}$) dermis layers, but the values of ${\alpha }_{44}$ 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.