Complementary analysis of Mueller-matrix images of optically anisotropic highly scattering biological tissues

Using optical techniques for tissue diagnostics (so-called ‘optical biopsy’) has been a subject of extensive research for many years. Various groups have been exploring different spectral and/or imaging modalities (e.g. diffuse reflectance spectroscopy, autofluorescence, Raman spectroscopy, optical coherence tomography (OCT), polarized light microscopy, etc.) for biomedical applications. In this paper, we report on using multi-wavelength imaging Mueller polarimetry combined with an appropriated image post-processing for the detection of tissue malignancy. We investigate a possibility of complementary analysis of Mueller matrix images obtained for turbid tissue-like scattering phantoms and excised human normal and cancerous colorectal tissue samples embedded in paraffin. Combined application of correlation, fractal and statistical analysis was employed to assess quantitatively the polarization-inhomogeneous scattered fields observed at the surface of tissue samples. The combined analysis of the polarimetric images of paraffin-embedded tissue blocks has proved to be an efficient tool for the unambiguous detection of tissue malignant transformation. A fractal structure was clearly observed at spatial distributions of depolarization of light scattered in healthy tissues in a visible range of spectrum, while corresponding distributions for cancerous tissues did not show such dependence. We demonstrate that paraffin does not destroy a fractal structure of spatial distribution of depolarization. Thus, the loss of fractality in spatial distributions of depolarization for cancerous tissue is related to the structural changes in the tissue sample induced by cancer itself and, therefore, may serve as a marker of the disease. The obtained results emphasize that a combined use of statistical, correlation and fractal analysis for the Mueller-matrix image post-processing is an effective approach for an assessment of variations of optical properties in turbid tissue-like scattering media and biological tissues, with a high potential to be transferred to clinical practice for screening cancerous tissue samples.


Background
Mueller polarimetry is the cutting edge optical technique widely used in various applications, associated with the studies of heterogeneous complex structures, including composite anisotropic materials, polymers, biological tissues, etc. [1][2][3][4][5].This technique analyzes the changes of polarization states of incident light reflected/transmitted/ scattered by a sample.In particular, a polarized light, which propagates through biological tissue, may change its state of polarization due to optical anisotropy of tissue and even become partially or fully depolarized because of scattering within a tissue.Typically, phase anisotropy (linear and circular birefringence), amplitude anisotropy (linear and circular dichroism) and scattering properties of tissues (depolarization power) are altered by structural and/or biochemical tissue malformation changes, e.g.due to inflammation, tumor growth, etc.Therefore, the physical parameters, which quantitatively define these changes, are considered as valuable markers for the diagnostic purpose [6][7][8][9][10][11][12].
The real 4 × 4 Mueller matrix contains complete information about polarization and depolarization properties of the scattering medium [13].In fact, the straightforward interpretation of matrix elements in term of optical properties suitable for quantitative characterization of studied medium can be performed for very limited class of samples (e.g.homogenous isotropic scattering media).Probing a multi-component highly heterogeneous structure of real biological tissue with polarized light produces very complex optical response due to the anisotropic scattering that lead to non-linear dependence of the elements of measured Mueller matrix on tissue optical properties [14].A phenomenological approach is used for a decomposition of measured Mueller matrix into the set of Mueller matrices describing the basic optical properties, such as dichroism, retardation and depolarization.Well-developed Mueller matrix algebra describes the numerous types of decompositions algorithms, including sum, product serial-parallel ones and others [15][16][17][18][19][20][21].The selection of a particular decomposition algorithm depends on the specific application.In present study of tissue phantoms and real tissues specimens we applied Lu-Chipman polar decomposition approach [18] and then complementary applied statistical, correlation and fractal analysis utilizing images of depolarization parameter Δ.
Most of modern polarimetric imaging systems [7][8][9][10][11] present the measured Mueller matrix data in form of two-dimensional (m × n) distributions of matrix elements, known also as Mueller matrix images (MMI) and denoted here as q(m × n).
For the quantitative evaluation of two-dimensional distributions q(m × n) the statistical moments of the first (Z 1 ), second (Z 2 ), third (Z 3 ), and fourth (Z 4 ) orders are used [22]: where N = m × n is the number of pixels in the lightsensitive zone of digital cameras.
The conventional autocorrelation approach [22] implemented for the assessment of spatial homogeneity by calculating the autocorrelation functions of each n th line of the image is used: The notation Δx stands for the 'pitch' of the coordinate (x) variation in the distribution q Ã ðxÞ ¼ qðxÞ−q , where q is the mean value, which is defined over the entire line.
The autocorrelation function is defined by summing partial dependencies: For the quantitative assessment of the relationships (3), the statistical moments of the second (K 2 ) and forth (K 4 ) orders, which characterize the FWHM and the sharpness of the autocorrelation function K(Δx) peak, respectively, were calculated.
The assessment of scale self-similarity of distributions q(m × n) is based on the fractal analysis, which includes the following steps [22]: (i) Calculation of power spectra J(q) using a discrete Fourier transform of the corresponding autocorrelation function; (ii) Determination of the log-log dependence of logJ(q) − log(ν), where ν = l −1 is the spatial frequency, and l is the size of the structural element in q(m × n); (iii)Dependencies logJ(q) − log(ν) are approximated by the least square fit with curves V(η) The straight segments of these curves serve to determine the slope angles η and calculate fractal dimensions.
The fractal classification of distributions q(m × n) takes place.If the dependencies V(η = const) are linear in the range of 2-3 decades of l sizes, the MMI are considered to be fractal.If there are several constant slope angles of V(η) curves, the distributions q(m × n) become multifractal.If there are no stable constant values of slope angle η within the whole range of sizes l, the distributions q(m × n) are believed to be random.
For the quantitative assessment of the logarithmic dependencies logJ(q) − log(ν) a second-order statistical moment was proposed in [22].
In current paper, we demonstrate new opportunities for differentiation of depolarizing samples by using combined application of statistical, correlational and fractal analysis of the image of depolarization parameter Δ.

Methods
Tissue phantoms and paraffin-embedded blocks of excised normal and cancerous colorectal tissues in vitro have been used in the experiments.The tissue phantoms contain rutile (TiO 2 ) particles (mean size 0.53 μm, standard deviation 0.01 μm, confirmed by electron microscopy measurements), embedded in 1 mm thick PVC-based host matrix (2 cm × 2 cm), with different concentrations: 1.5, 3.0, and 6.0 mg/ml.The concentration values were selected to assure the scattering coefficients μ s of the phantoms vary providing change from single scattering regime (optically thin layer) to multiple scattering regime (optically thick layer).From spectrophotometric measurements the intrinsic absorption coefficient was found to be negligible, therefore, all three tissue phantoms were considered as pure scattering ones.The details on phantom fabrication and characterization are described in [23].
The paraffin-embedded blocks of excised colorectal human tissue were prepared according to standard protocol utilized for the pathology analysis.First, the tissue specimens were fixed in formalin for several hours, then they were dehydrated by immersion in alcohol, cleared by organic solvent to remove alcohol and finally were infiltrated with paraffin wax.When molten paraffin solidified, it provided a support matrix for tissue.Then thin sectioning of two selected blocks was performed.A conventional histological analysis of stained tissue sections by pathologist did not find malignancy on the histological cut from the first block and confirmed the presence of malignancy in tissue from the second block.
The experimental measurements with phantoms were performed in the transmission mode with Mueller polarimetric microscope, whereas the measurements with the thick tissue blocks were performed in reflection mode by using multi-wavelength imaging Mueller-matrix polarimeter.The detailed description and specifications of the experimental system used in current study are presented in [19,24], and omitted here for brevity.We applied Lu-Chipman decomposition of experimental Mueller matrices.Neither significant diatteniation nor retardance were found for any class of studied samples.The depolarization (Δ) maps were used for further analysis.
The set of chosen samples provides the possibility of conducting a comparative analysis of the dependencies of statistical (Eq.( 1)), correlational (see Eqs. (2-3)) and fractal parameters (Eq.( 4)).These parameters characterize the spatial distribution of depolarization parameter Δ values in terms of different scattering multiplicities of the model samples with controlled parameters ("phantoms") and samples of real partially-depolarizing tissues with various pathological conditions.
The analysis of depolarizing properties of tissue phantoms 1-3 in the frameworks of the statistical, correlation and fractal approaches reveals that with the increase in number of scattering events along the light propagation path (i.e. increase of scattering coefficient from 2.5 mm − 1 to 10 mm − 1 ) the maximum of the distribution shifts to larger value of depolarization (see Figs. 1, 1.4-1.6).The dispersion of Δ values distribution increases monotonically with the increase of mean value of depolarization (see Figs. 1, 1.4-1.6).This fact correlates well with the measurements of the degree of depolarization, averaged over an emerging beam area [19].On the other hand, the different phantom samples are characterized by the "individual" half-width, asymmetry and sharpness of the peak of distribution N(Δ).
As one can see the depolarization Δ(m × n) maps of tissue phantoms 1 and 2 are spatially heterogeneous.This is indicated by the fast decaying dependencies of the autocorrelation functions K Δ (Δx) (see Figs. 1, 1.7 and 1.8).This fact may be related to coordinate heterogeneity of the sample structure that was described in details early [22].Consequently, different scattering multiplicity appears in different parts of the phantom layer.Optically, this is displayed in the coordinate fluctuations of the values of depolarization degree.With the increase of scattering multiplicity, these fluctuations decrease as it is obvious from the K Δ (Δx) peak's FWHM raise for the sample 3 (see Fig. 1, 1.6).This tendency is confirmed by the correlation analysis of the depolarization map of the partially depolarizing sample 3. Two-dimensional distribution of depolarization parameter Δ for this sample is characterized by a larger range of values.As a result, the autocorrelation function K Δ (Δx) shows smooth decreasing dependency (see Fig. 1, 1.9) [22].
Scale self-similar structure of the depolarization maps of phantom samples 1-3 essentially depends on the scattering multiplicity.Therefore, for the samples 1 and 2 the logarithmic dependencies logJ(Δ) − log l −1 do not have a constant slope angle η of the approximating curves V(η) in the entire range of the structural elements dimensions l of the depolarization maps (see Figs. 1, 1.10, 1.11).Distribution Δ(m × n) for the sample 3 is multi-fractal (two values of η) (see Fig. 1, 1.12).Thus, the increase of scattering multiplicity leads to almost equiprobable contribution of the different-scale structural elements of the phantom samples to the polarization structure of the optical object field formation.Along with this, the raise of the depolarization parameter value is accompanied by an increase of the values for low (~10 μm) and medium (~100 μm) sizes l of structural elements in the dependencies logJ(Δ) − log l −1 .Apparently, this regularity is conditioned by the growth of the high-frequency component in the distribution Δ(m × n) due to the increasing scattering multiplicity of optically-thick phantom sample (scattering coefficient of 10 mm −1 ).The results of the statistical, correlational and fractal analysis of Δ(m × n) for the wavelength 0.55 μm are shown in the Table 1.
Analysis of the obtained results revealed that the following parameters are the most sensitive to the increase of scattering multiplicity within the volume of phantom samples with the scattering coefficient varying from 2.5 mm − 1 to 10 mm − 1 : the dispersion of distribution Δ(m × n) (increases up to 24 times); the 4th order statistical moment, which characterizes the sharpness of the peak of distribution of depolarization parameter Δ values (decreases 20 times).We found following behavior of the parameters with the increase of the scattering coefficient from 2.5 mm −1 to 10 mm −1 : Thus, a high sensitivity of the set of statistical (Z i ), correlational (K j ), and fractal (D f ) parameters to the changes of two-dimensional distributions of depolarization degree Δ of structurally similar samples with different scattering coefficients was identified.These findings open the perspectives for the quantitative characterization of partially depolarizing optically anisotropic biological tissues, which are much more complex samples compared to optically isotropic tissue phantoms.
The peaks of spatial distributions of depolarization parameter Δ for the samples of paraffin-embedded tissue blocks are localized within the range 0.25 ≤ Δ ≤ 0.85 (see Fig. 2a-c, panels 3, 4).The comparison of distributions of depolarization Δ for the phantom samples (see Figs. 1, 1.4-1.6)and biological tissue reveals that latter distributions are multi-modal with several peaks.With the increase of wavelength the peaks of histograms of depolarization parameter Δ are shifted to larger values for both healthy and cancerous samples of paraffin-embedded tissue blocks (see Fig. 2a-c), panels 3, 4).Besides, the range of Δ values expands.This fact  indicates the 'longwave' increase of depolarization of the radiation rearranged by the biological samples.We attribute this fact to deeper penetration of longer wavelengths into the tissue.Increased optical path results in larger number of scattering events, which randomize the polarization of incident light.Furthermore, the depolarization parameter Δ of cancerous paraffin-embedded tissues is higher than that of the normal tissue sample.It is worth to mention that an opposite trend on depolarization was observed in the experiments with fresh thick tissue specimens (colon, uterine cervix) [8,24] when epithelial surface of tissue was imaged.In the above-mentioned measurement configuration [8,24] an imaging plane was orthogonal to the plane of tissue histological cuts seen and analyzed by pathologist.
As one can see, spatial depolarization distributions Δ(m × n) for the paraffin-embedded tissue blocks have complex coordinate-heterogeneous topographic structure (see Fig. 2a-c, panels 1, 2) compared to the similar maps obtained for the tissue phantoms (see Figs. 1, 1.1-1.3).Obviously, this difference is related to the morphological structure of biological tissue.Quantitatively, the topographic heterogeneity in the spatial distribution of depolarization Δ(m × n) is illustrated by the trend of autocorrelation functions K Δ (Δx) achieved for healthy and cancerous paraffin-embedded tissue blocks (see Fig. 2a-c, panels 5, 6), which have a larger half-width K 2 ↑ and less sharp peak K 4 ↓ compared to the phantom samples due to a larger range of values of the depolarization parameter Δ (see Tables 1 and 2).
The analysis of data presented in Table 2 revealed the most sensitive parameters for the differentiation of normal and cancerous tissue: 1 st group (non-colored boxes): the difference between the parameter values for cancerous and non-cancerous tissues does not exceed 25-45 percent; 2 nd group (highlighted in yellow): the difference between the corresponding values of parameters varies between 50 and 100 percent; 3 rd group (highlighted in green): the difference between the corresponding values of parameters is two to five-fold.
Spectral analysis shows that with the increase of probe beam wavelength, coordinate uniformity of the distributions Δ(m × n) grows for both healthy and cancerous samples due to the depolarization enhancement.From the quantitative point of view, this trend is illustrated by the increase of the FWHM and decrease of the peak sharpness of the dependencies K Δ (Δx) (see Fig. 2a-c, panels 5, 6).As it can be seen from Table 2, along with the wavelength increase (λ↑) the parameter K 2 increases and the parameter K 4 decreases.
Comparative analysis of the logarithmic dependencies logJ(Δ) − log l −1 calculated for the distributions of depolarization Δ(m × n) of biological samples, reveals significant differences for normal and cancerous tissues.It has been shown, that the distributions of depolarization Δ(m × n) for normal tissue are fractal in all spectral range (see Fig. 2a-c, panels (7)).For cancerous sample, the distributions of depolarization Δ(m × n) are random in the region of medium sizes (~100 μm − 300 μm) (see Fig. 2a-c, panels (8), rectangular boxes).The discovered pattern is in a good correlation with the data obtained by Mueller-matrix mapping of optically-thin histological cuts of different organ tissues (prostate, cervix and uterine wall) [22] where the 'oncological destruction' of fractality of two-dimensional structure of depolarization maps is associated with the formation of new fibrillar networks.

Conclusions
The combined application of statistical, correlation and fractal analysis for the quantitative assessment of polarization-inhomogeneous scattered fields observed at the surfaces of isotropic scattering tissue phantoms and Table 2 Statistical (Z i ), correlation (K j ) and fractal (D f ) parameters introduced for quantitative assessment of depolarization Δ(m × n) of light scattered in the normal and cancerous paraffin-embedded tissue blocks, measured at the selected wavelengths biological tissue samples has been performed in these initial exploratory studies.It has been shown that using the results of the analysis of paraffin-embedded tissue blocks one can unambiguously detect the malignant transformation of tissue.The distributions of depolarization Δ(m × n) for healthy tissue are fractal at all studied wavelengths, while corresponding distributions for cancerous tissue loose fractality for medium size features (few hundreds of microns).Adding paraffin to tissue alters its optical properties and increase scattering [25].Our experiments with isotropic scattering phantoms have proved that simple increase in scattering (from single scattering regime to multiple scattering regime) does not destroy the fractality of distributions of Δ(m × n).It suggests that possible difference in paraffin intake by healthy and cancerous tissue can not be the reason for the observed loss of fractal properties of two-dimensional distributions of depolarization Δ(m × n) for cancerous tissue.The most plausible explanation of this effect is related to the structural changes of tissue induced by cancer development.
The obtained results suggest that Mueller-matrix polarimetry can be an effective approach for screening optical anisotropy variations in tissue-like highly scattering media, with a high potential in clinical application for diagnosis of cancerous tissues.Using the thick blocks of excised tissue for the preliminary optical analysis by pathologist may considerably reduce the time and cost of diagnostics.To capitalize on our initial findings the measurements and statistical, correlation and fractal analysis of larger number of paraffin-embedded tissue samples will be undertaken in a future work.

Fig. 1
Fig. 1 Statistical, correlation and fractal characteristics counted based on the experimentally measured depolarization maps for three phantom samples

Table 1
Statistic (Z i ), correlation (K j ), and fractal (D f ) parameters used for characterization of depolarization distributios Δ(m × n) of three phantom samples