3D Mueller matrix mapping of layered distributions of depolarisation degree for analysis of prostate adenoma and carcinoma diffuse tissues

Prostate cancer is the second most common cancer globally in men, and in some countries is now the most diagnosed form of cancer. It is necessary to differentiate between benign and malignant prostate conditions to give accurate diagnoses. We aim to demonstrate the use of a 3D Mueller matrix method to allow quick and easy clinical differentiation between prostate adenoma and carcinoma tissues with different grades and Gleason scores. Histological sections of benign and malignant prostate tumours, obtained by radical prostatectomy, were investigated. We map the degree of depolarisation in the different prostate tumour tissues using a Mueller matrix polarimeter set-up, based on the superposition of a reference laser beam with the interference pattern of the sample in the image plane. The depolarisation distributions can be directly related to the morphology of the biological tissues. The dependences of the magnitude of the 1st to 4th order statistical moments of the depolarisation distribution are determined, which characterise the distributions of the depolarisation values. To determine the diagnostic potential of the method three groups of histological sections of prostate tumour biopsies were formed. The first group contained 36 adenoma tissue samples, while the second contained 36 carcinoma tissue samples of a high grade (grade 4: poorly differentiated—4 + 4 Gleason score), and the third group contained 36 carcinoma tissue samples of a low grade (grade 1: moderately differentiated—3 + 3 Gleason score). Using the calculated values of the statistical moments, tumour tissues are categorised as either adenoma or carcinoma. A high level (> 90%) accuracy of differentiation between adenoma and carcinoma samples was achieved for each group. Differentiation between the high-grade and low-grade carcinoma samples was achieved with an accuracy of 87.5%. The results demonstrate that Mueller matrix mapping of the depolarisation distribution of prostate tumour tissues can accurately differentiate between adenoma and carcinoma, and between different grades of carcinoma. This represents a first step towards the implementation of 3D Mueller matrix mapping for clinical analysis and diagnosis of prostate tumours.


Methods and theory
Biological samples. For determination of the type of tumour, native samples of histological sections of the examined prostate tissue were made. This study was conducted in accordance with the principles of the Declaration of Helsinki, and in compliance with the International Conference on Harmonisation-Good Clinical Practice and local regulatory requirements. Ethical approval was obtained from the Ethics Committee of the Bureau of Forensic Medicine of the Chernivtsi National University and the Bukovinian State Medical University (Chernivtsi, Ukraine), and written informed consent was obtained from all subjects prior to study initiation. The type of prostate tumour was determined by an independent assessment of stained histological samples (Fig. 1). The histological analysis was conducted by the following procedure: 1. Fixation of prostate tissue with formalin (40% formaldehyde aqueous solution); 2. Washing of samples in running water for 24 h; 3. Dehydration with alcohols with increasing concentration (70-100%) within 48 h; 4. Fixing the material in a xylene-paraffin mixture for 1-2 h, at a temperature of 52°-56°, and cutting out a block with a sample enclosed in it; 5. Making histological sections on a standard freezing microtome; 6. Staining of histological sections with haematoxylin and eosin ( Fig. 1 shows real colour microscopic images); 7. Microscopic examination of images of the obtained samples with differentiation of their structure by grades and determination of the position of the prostate tumor sample according to the Gleason Pattern scale.
Three representative groups of histological biopsy sections, obtained from radical prostatectomy, were formed: • Group 1 consisted of n = 36 adenoma samples; • Group 2 consisted of n = 36 moderately differentiated (3 + 3 on Gleason's pattern scale) carcinoma samples; • Group 3 consisted of n = 36 poorly differentiated (4 + 4 on Gleason's Pattern scale) carcinoma samples. Table 1 presents the optical and geometric parameters of the samples of native histological sections of prostate tumour biopsies from each of the groups. The geometric thickness of the histological sections was determined according to the standard values of the freezing microtome scale. The extinction coefficient was determined by measuring the attenuation of the illuminating beam intensity as it passed through the sample. This was achieved using an integral light-scattering sphere following a well-established method 30,31 . The measurement of the integral  www.nature.com/scientificreports/ degree of depolarisation ( ) was performed using a standard Mueller matrix polarimeter in accordance with similar measurements previously underatken 23 . By comparing the optical and geometric parameters of the different histological section groups, we observe that the depolarisation degree is high for all samples. Additionally, there is little difference in the depolarisation values for the different groups. As discussed below, this significantly limits the diagnostic potential of traditional 2D Mueller matrix polarimetry.
Theoretical background. The diagnostic efficiency of detection of oncological conditions by applying the technique of 2D Mueller matrix polarimetry to the polycrystalline component of optically thin (non-depolarising) histological sections of biological tissues of various human organs has previously been demonstrated [32][33][34][35] . These previous studies are devoted to the search for, and subsequent diagnostic application of, a set of diagnostically relevant relationships between Mueller matrix images ( F ik x,y ) and maps of linear and circular birefringence ( LB x.y , CB x,y ) and dichroism ( LD x.y , CD x.y ): As a general rule, samples of biological tissues of human organs are diffuse and thus strongly depolarize optical radiation (as is the case with our prostate tissue samples). This significantly limits the potential for differential diagnostics by traditional 2D Mueller matrix polarimetry. In this situation, almost all off-diagonal ( F ik,i =k ) elements of the Mueller matrix of the diffuse biological layer are significantly reduced 5,17,21 . However, the diagonal elements ( F ik,i=k ) are still diagnostically relevant. The superposition of the diagonal elements determines the overall value of depolarisation for the optical radiation 11,14,15 .
The value of the parameter is an integral equivalent of the overall optical properties of the object under investigation.
In all other cases, the value of the integral depolarisation parameter is determined by a combination of two components. Firstly, the 'local' component resulting from the formation of an orthogonal component of the laser radiation amplitude (i.e. a change in the initial state of polarisation) due to the optical anisotropy of the biological layer. We will call this the "A-mechanism". Secondly, the 'diffuse' component resulting from the statistical averaging of the polarisation state due to the superposition of laser waves, scattered in the volume of the biological layer, with different states of polarisation. We will call this the "B-mechanism".
For optically thin ( τ ≤ 0.01 ), singly-scattering, optically anisotropic biological layers, the A-mechanism is dominant. The influence of the optical anisotropy (linear and circular birefringence and dichroism) is seen through a slight decrease in the values of the diagonal elements F 22;33;44 < 1 of the Mueller matrix 23 . Therefore, in comparison with an optically isotropic layer, there is an increase in the value of depolarisation ( � > 0 ). The specific value of this parameter is interrelated with the specificity of the polycrystalline structure of biological tissues 10, 32,35 .
As the optical thickness increases ( τ > 0.01 ), the multiplicity of the light scattering (i.e. the number of scattering events each photon experiences on average) in the volume of the histological sections increases (2) F ik,i =k → 0; F ik,i=k = 0; www.nature.com/scientificreports/ correspondingly. The value of the integral depolarisation parameter is determined by the superposition of the effects of the A-and B-mechanisms. For small values of the attenuation coefficient ( τ ≤ 0.5 ), the contribution of the two optical mechanisms is comparable. Averaging the polarisation states of the scattered (differently polarized) wavefronts leads to a further increase in the degree of integral depolarisation. There is hence a decrease in the dependence of integral depolarisation on the specific morphological structure of biological tissue. For diffuse (optically thick, τ > 1.5 ) biological layers, the polarising effects of the individual features of the polycrystalline structure are effectively smoothed by the B-mechanism. Therefore, the effectiveness of differential 2D diagnostics is extremely low or unsatisfactory in this limit. The 2D polarisation diagnostics of frozen prostate tissue samples is limited by the impossibility of obtaining optically thin ( h ∼ 10−15 µm ) histological sections with geometric consistency between samples. To obtain consistent samples, a geometric thickness h ∼ 40 µm was required. Hence, the B-mechanism is dominant for all samples, with the diffuse scattering resulting in both high and comparable values of the integral depolarisation for all the types of prostate tissues investigated (Table 1) ) into the plane of the polarisation-inhomogeneous image of the object. The angle between the illuminating and reference beams was selected as = 4 • . As a result, an interference pattern is formed, the period of which is 10 µm . According to the Nyquist-Shannon sampling theorem, this periodicity ensures reliable recording of a single fringe by 5 pixels of the digital camera (15, see Fig. 2).
The error in setting the transmission axes of the polarisers 8,11,13,14,36 was �γ = 0.2 • . The phase shift error of the quarter-wave plates 13,37 was �δ = 0.5 • . The error in setting the fast axis of the quarter-wave plates was �ϕ = 0.2 • . The resulting error in the linear polarisation states did not exceed 0.0021, and that in the circular polarisation states did not exceed 0.0056 38 .
Before carrying out measurements on histological sections of the prostate, metrological certification of the experimental setup was conducted using model objects (air, linear polariser, quarter-and half-wave plates). From 50 measurements for each model object, the errors in the determination of the diagonal elements of the Mueller matrix were determined: for F 22;33 -0.25% and for F 44 -0.5%.

Determination of the layered depolarisation distributions.
When obtaining the polarisation interference images, one can determine the Mueller matrix elements by the following procedure: 1. Forming six distinct polarisation states of both the illuminating and reference beams: four linear polarisations at angles of 0°, 90°, 45°, and 135°, and two circular polarisations-right ( ⊗ ) and left ( ⊕ ) circular polarisation respectively. 2. Collection of six images for each partial interference pattern passed through the polariser-analyser [14] with sequential orientation of its transmission axis at angles 0 = 0 • ; 90 = 90 • . We then have a series of twelve digital photo of the original field with the applied reference wave. www.nature.com/scientificreports/ 3. For each partial interference distribution, we perform a two-dimensional discrete Fourier transform ( DFT(υ, ν)) on the image. The two-dimensional DFT(υ, ν) of a two-dimensional array I 0;90 x,y (i.e. the image) f(x,y)is a function of two discrete variables coordinates x,y the DFT(υ, ν) function is defined by by 39 : The results of this transformation should contain three peaks, one central (main) peak and two additional side peaks. The Fourier transform acts like a low-pass filter. It removes interference fringes, which enables extraction of the complex representation of the real field due to the object. Also, since the extracted part has limited size, it acts like a low-pass filter for the object field too. 5. Either of the additional side peaks (in complex representation) can be used to create a new Fourier spectrum by first extracting the peak and then placing it into centre of a newly generated spectrum DFT 0 where n is the birefringence; is the wavelength; and h is the sample thickness of the object field, separated by an arbitrary step of �θ. 8. In each phase plane θ k , the corresponding sets of parameters of the Stokes vector of the object field of the biological layer are calculated: The obtained relations (9) form the basis of the method for determining layered distributions of the values of the Mueller matrix elements F ik x, y, θ k . The elements are obtained from the following relationships:    (10) Here, K is the total number of pixels of the CCD-camera. These parameters characterize the mean value ( Z 1 ), dispersion ( Z 2 ), skewness ( Z 3 ) and kurtosis ( Z 4 ) of the distributions � x, y, θ k .
To determine the statistical significance of a representative sampling of the number of samples by the crossvalidation method 39 , the standard deviation σ 2 of each of the calculated values of the central statistical moments Z i=1;2;3;4 (n) (15)(16)(17)(18), which characterise the distribution of the degree of depolarisation (4), was determined. The specified number (36 for each group) of samples provided the level σ 2 ≤ 0.025 . This standard deviation corresponds to a confidence interval p < 0.05 , which demonstrates the statistical reliability of the 3D Mueller matrix mapping method for layered depolarisation maps within a representative sample.
Diagnostic method. Our proposed optical procedure, used herein for the differential diagnosis of the prostate tissue samples, comprised the following steps: www.nature.com/scientificreports/ 1. We determined the optimal phase plane for further analysis and diagnostics. A "macro" step of discrete phase scanning was chosen-�θ max k = 0.25 rad . Algorithmically, a series of layered coordinate distributions of the depolarisation degree � x, y, θ k corresponding to each �θ max k = 0.25 rad was reconstructed. 2. For each phase section θ max k of the object field of scattered radiation, the set of central statistical moments of the first to fourth orders Z i=1;2;3;4 , which characterize the distributions � x, y, θ k , was calculated.

The difference (�Z
between the values of each of the statistical moments for the different phase planes was calculated. 4. The phase interval �θ * = θ max j+1 − θ max j was determined, within which the monotonic growth of the value

5.
Within the limits of �θ * , a new series of values �Z i = Z i θ min q+1 − Z i θ min q was calculated with a discrete phase scanning "micro" step �θ min q = 0.05 rad. 6. The optimal phase plane θ * was determined, in which �Z i (θ * ) = max . The specified analytical procedure for one sample of prostate tissue took less than 7 min. A more accurate phase scanning algorithm can be chosen, for example with a step of 0.01 rad. In this case the processing time increased to 11 min. 7. In the plane θ * , the mean Z * i=1;2;3;4 and standard deviations σ �Z * i were determined within the representative samplings of histological sections from group 1, group 2 and group 3. 8. To differentiate benign and malignant tumours, for each of the statistical moments Z i=1;2;3;4 , the sensitivity ( Se 12 = a 12 a 12 +b 12 100% ; Se 13 = a 13 a 13 +b 13 100% ), specificity ( Sp 12 = c 12 c 12 +d 12 100% ; Sp 13 = c 13 c 13 +d 13 100% ) and balanced accuracy ( Ac 12 = The position of the optimal phase plane θ * , determined by the indicated methods for each of the three groups separately, coincided with an accuracy of 0.02 rad. The operating characteristics of the method 40 (sensitivity Se , specificity Sp , accuracy Ac ) remained practically unchanged. Figure 3 shows the phase dependences of the integral magnitude (averaged over all the pixels of the CCD-camera) of the degree of depolarisation �(θ k ) of laser radiation passing through prostate tissues from each of the three groups defined previously. The degree of depolarisation was experimentally determined using the 3D Mueller matrix mapping method described above. www.nature.com/scientificreports/ Analysis of the data obtained (see Fig. 3) revealed that, despite the overall depolarisation through the whole samples being more or less the same (Table 1), there are significant differences observed within the volume of the tissue samples. The overall ranges of �(θ k ) by histological biopsy sections of all samples span (15%−19%)for(θ k = 0.25 rad) ≤ � ≤ (78%−81%)for(θ k = 1.5 rad) . The rate at which the dependencies �(θ k ) change is very different though. The most rapidly increasing depolarising ability is seen for samples of histological sections of adenoma biopsy (see Fig. 3, dark curve), while the least rapid is for samples of native sections of grade 4 carcinoma (see Fig. 3, red curve), with grade 1 carcinoma (see Fig. 3, blue line).

Results
For small values of the phase ( θ k ≤ 0.5 rad ), the dominant contribution to the depolarisation formation scenario is the "A-mechanism". In general, soft biological tissues (including the prostate) have a relatively low level of optical anisotropy. The magnitude of the circular and linear birefringence n does not exceed 1.5 × 10 −3 [5][6][7][8] . Therefore, for small values of the phase, the integral depolarisation is low for histological sections of all types of prostate tumours. With increasing θ k , the contribution of the "B-mechanism" to the depolarising ability of the prostate samples increases. The contribution of the "B-mechanism" characterizes the measure of diffraction expansion of partial waves on optical inhomogeneities of biological tissue. The most pronounced effect of this mechanism manifests on small-scale (well differentiated) structures of the prostate tissue. Therefore, the most rapidly growing value of �(θ k ) is found for samples of histological sections of prostate adenoma (see Fig. 3, dark curve), the morphological structure of which is at the smallest scale. Conversely, the values of �(θ k ) for histological sections of biopsy of poorly differentiated prostate carcinoma increase most slowly (see Fig. 3, red curve).
For the largest values of phase shifts ( θ k ≥ 1.5 rad ), the diffraction "B-mechanism", which characterizes the multiple superposition of laser waves scattered in the volume of the biological layer, becomes dominant. The degree of depolarisation then reaches a maximum level for all types of samples. The maximum differences between the values of integral depolarisation by samples of histological sections of biopsy of prostate tumours are realized for a certain "intermediate" range of phase shifts ( 0.75 rad < �θ * < 1.05 rad ). The maximum difference is found at θ * = 0.85 rad , so we choose this value for further investigations.
Figures 4 and 5 present exemplar depolarisation maps obtained for samples of prostate adenoma (see Fig. 4), and carcinoma with a 3 + 3 Gleason score (see Fig. 5). However, no obvious differences or relation to the tissue structure are immediately visible. For carcinoma with a higher (4 + 4) Gleason score, more obvious differences are visible. By analysing the results obtained, we can see that there is a complex and individual topographic structure of the depolarisation maps � x, y, θ k = 0.85 rad of native histological sections of prostate tumours obtained during radical prostatectomy (Fig. 6).
We can then assess the depolarising ability of histological biopsy sections of benign and malignant prostate tumours within the framework of a statistical approach (15)(16)(17)(18) to the analysis of the distributions � x, y, θ k . The largest mean and maximum range of coordinate fluctuations of the magnitude of the degree of depolarisation � x, y, θ k = 0.85 rad are found for histological sections of benign adenoma biopsies. The smallest mean and minimum range of coordinate fluctuations of the magnitude of the degree of depolarisation � x, y, θ k = 0.85 rad are found for histological sections of poorly differentiated carcinoma biopsies. The experimentally revealed differences in the series of maps x, y in the phase section θ * = 0.85 rad can be associated with the different degrees of differentiation of the structures within the prostate tumour tissues. For well differentiated adenoma tissue, diffraction effects are most pronounced. Therefore, the average value and fluctuations of the degree of depolarisation are maximal. For poorly differentiated (4 + 4) carcinoma tissue, the diffraction broadening of laser waves is smaller. Therefore, both the mean and variance of the individual coordinate values � x, y, θ k = 0.85 rad decrease.
The results of quantitative statistical analysis and comparison of the distributions � x, y, θ k for the different groups of tissues are presented in Table 2. The data obtained shows that statistical moments of higher orders, which characterize the skewness ( Z 3 (θ k ) ) and kurtosis ( Z 4 (θ k ) ) of the distributions of the degree of   Table 2. Intergroup differences in the magnitude of the central statistical moments of the 1st-4th orders Z i=1;2;3;4 , which characterize the distributions of the magnitude of the degree of depolarisation � x, y, θ k in the layered phase sections corresponding to θ k = 0.5 rad , θ k = 0.85 rad , and θ k = 1.5 rad.

Adenoma -Carcinoma (3+3)
A  x,y, θ k ↓ . Table 3 shows the changes in the magnitude of the sensitivity Se , specificity Sp , and balanced accuracy Ac of the diagnostic power of the 3D Mueller matrix mapping method for native histological sections of benign and malignant prostate tumours.

Conclusions
The method of 3D Mueller matrix mapping of diffuse biological layers with digital holographic reconstruction of layered depolarisation maps (1-18) is analytically substantiated. This method, built on the platform of polarisation interferometry, was experimentally tested for the task of express (less than 15 min) differential diagnosis of prostate tumours obtained during radical prostatectomy. A representative sample of diffuse samples of and variously differentiated (Gleason scores 3 + 3 and 4 + 4) carcinoma biopsies was studied. Using the method of phase scanning of the object field of histological sections, the optimal cross section ( θ k = 0.85 rad ) was determined. In this cross-section, the maximum differences ( Z i=1;2;3;4 → max ) between the values of the central statistical moments Z i=1;2;3;4 , which characterize the depolarisation maps of samples of benign and malignant prostate tumours, were realized. The most sensitive statistical parameters were established to be the skewness and kurtosis of the coordinate distributions � x, y, θ k . The operational characteristics (sensitivity Se , specificity Sp and balanced accuracy Ac ) have been determined, which demonstrate the diagnostic power of the 3D Mueller matrix mapping of diffuse biological layers with digital holographic reconstruction of layered depolarisation maps. High accuracy ( 90.3% ≤ Ac 12;13 (θ * = 0.85 rad) ≤ 93.1% ) of differentiation of benign and malignant samples of native histological sections of biopsy of prostate tumours was achieved. In addition, the possibility of diagnosing samples of variously differentiated malignant tumours with accuracy Ac 23 (θ * = 0.85 rad) = 87.5% has been demonstrated. Table 3. Selectivity, specificity, and balanced accuracy of the diagnostic power of the 3D Mueller matrix mapping method of native histological sections of benign and malignant prostate tumours in the layered phase sections corresponding to θ k = 0.5 rad , θ k = 0.85 rad , and θ k = 1.5 rad.