Polarimetric data-based model for tissue recognition

: We highlight the potential of a predictive optical model method for tissue recognition, based on the statistical analysis of different polarimetric indicators that retrieve complete polarimetric information (selective absorption, retardance and depolarization) of samples. The study is conducted on the experimental Mueller matrices of four biological tissues (bone, tendon, muscle and myotendinous junction) measured from a collection of 157 ex-vivo chicken samples. Moreover, we perform several non-parametric data distribution analyses to build a logistic regression-based algorithm capable to recognize, in a single and dynamic measurement, whether a sample corresponds (or not) to one of the four different tissue categories.


Introduction
Polarimetry comprises a set of non-invasive and non-contact testing techniques that provide information about the optical properties of a sample [1]. Those techniques appear as very useful tools to be applied in different research fields, such as material classification [2], gas sensing [3], astronomy [4], biomedicine [5][6][7][8], and remote sensing [9] but among all, for the purpose of this work, they are an especially useful for analyzing and characterizing biological samples [10][11][12]. Biological tissues are turbid media that strongly scatter light and present certain characteristic degree of anisotropy due to their ultrastructure components (e.g., the dense collagen and elastin matrix with oriented fibers, etc.) [13] so these properties lead to certain polarimetric response (depolarization, retardance, etc.). In several pathologies, morphological changes occur that alter the polarization properties of tissues [14]. Therefore, polarimetry constitutes a powerful and promising optical tool for the study and classification of biological samples [11,15]. Moreover, polarimetry can be combined with other optical techniques such as multispectral imaging [15][16][17][18][19][20][21] or microscopy [22][23][24].
The Mueller-Stokes (M-S) is a widespread mathematical formalism which allows to characterize the state of polarization of light beams (Stokes vector) and the polarimetric properties of samples (Mueller matrices, M). The Mueller matrix, MM, describes the effect of a material media, such as an optical element, or a sample, upon the polarization state of light by generating an output Stokes vector, S out , from an input light state of polarization, S in , for linear processes [25,26]. One interesting property of the (M-S) formalism is that it is based on radiometric measurements and therefore is easy to be implemented to interpret experimental data. In addition, unlike other formalism used to describe polarization, as for instance, the Jones formalism, the M-S is valid to deal with fully polarized light, but also with depolarization content. In of the statistical analysis conducted on polarimetric data (subsection 3.2). In section 4, we describe the construction of an optical probabilistic model suitable for the classification of biological samples, which is based on a logistical regression applied to polarimetric data, and we discuss the validity of developed models in terms of sensitivity and specificity, as well as, their suitability for imaging classification of tissues. Finally, the main conclusions of the study are provided in section 5.

Methods
In this section, we present the mathematical background necessary to build the proposed tissue recognition model. The model is based on a set of metrics derived from experimental Mueller matrices of biological samples. Those observables allow us to obtain a clear interpretation of some physical properties of the samples. We also provide the experimental methodology followed for the completion of the study.

Polarimetric metrics
Recalling the M expression in Eq. (1) we define the diattenuation vector D as the metric that describes the property of an optical element by which the intensity of the exiting beam depends on the polarization state of the incident beam. Analogously, the polarizance vector P describes the polarization capability of a sample, when illuminated by an unpolarized input state [41]. Such magnitudes can be computed and described directly from the combination of M coefficients by using: where m ab (a,b = 0, 1, 2, 3) represents an element of M.
Other polarimetric characteristics of a sample are also encoded in its Mueller matrix in a complex way. This information can be synthesized in a product of three pure Mueller matrices, which are functions of well-defined polarimetric magnitudes for an easier physical interpretation, by employing the so-called Lu-Chipman decomposition [32]: where m 00 is related to the intensity transmittance or reflectance of the sample and those pure matrices are defined as depolarizers (M ∆P ), pure retarders (M R ) and pure diattenuators (M D ). Therefore, the information related to the depolarization, retardance and diattenuation properties of samples is encoded in different scalar polarimetric observables derived from those pure matrices. In particular, the total retardance R, the linear retardance δ and the optical rotation Ψ are obtained from the pure retarder matrix M R in this way [28]: where Tr denotes for the trace of a matrix, and M Ri,j are the coefficients of the pure retarder matrix. The total retardance (R) allows the description of the global behavior of a general retarder (i.e., an elliptical retarder), as it combines the linear retardance and the optical rotation [19].
However, in the context of tissue discrimination, to include individual characteristics of the retarder constituents, as the linear retardance (δ) or the rotation applied to the linear retarder (Ψ ), may be of interest. Thereby, in the present work, all them will be analyzed individually in order to consider all the contributions.
Further description of a sample can be achieved by studying its depolarizing properties. To this aim, it is worth using the indices of polarimetric purity (IPP) [11,12,20]: three invariant depolarizing indicators that contain complete and detailed information of the capability of samples to induce randomness to input polarization states [29]. The nature of the IPPs lies in the statement that any depolarizer response can be synthesized as an incoherent sum of four non-depolarizing components whose relative statistical weights are performed by combinations of IPPs. In this way, we can discriminate different types of depolarizers by only looking at the statistical weight of each pure component.
The IPP observables are obtained from the four eigenvalues λ i (which fulfill the condition λ 0 ≥ λ 1 ≥ λ 2 ≥ λ 3 ) of the covariance matrix H(M), a positive semidefinite Hermitian matrix whose elements are defined through linear combinations of the Mueller matrix [20]. Particularly, the IPPs (real magnitudes P 1 , P 2 and P 3 ) are defined as: In addition, the depolarization index P ∆ , which estimates the overall depolarization of a M, is built from IPPs as follows: These indicators characterize a pure (non-depolarizing) system when P 1 = P 2 = P 3 = P ∆ = 1 is found [42]. Due to eigenvalue restrictions, IPPs are limited by the 0 ≤ P 1 ≤ P 2 ≤ P 3 ≤ 1 inequalities. The ideal depolarizer is described for P 1 = P 2 = P 3 = P ∆ = 0. Moreover, IPPs allow synthetizing the global information provided by P ∆ , allowing to transform a 1-dimensional space into a 3-dimensional information space.

Experimental methodology and statistical analysis
In this section, we provide the experimental methodology and the description of the inspected tissue samples (subsection 3.1) as well as the description of the data analysis conducted on Mueller matrix-based polarimetric data (subsection 3.2). The study provided in this section lays the foundation for the subsequent construction of an optical model suitable for the classification of biological samples, that will be discussed in section 4.

Experimental methodology and sample description
In the current study, the experimental Mueller matrices of 157 biological samples are measured at three different wavelengths covering the visible range (625 nm, 530 nm and 470 nm), since different wavelengths are associated with different light penetration capability in tissues [43]. Such measurements are performed in a scattering configuration (capturing diffuse light) by using a complete Mueller imaging polarimeter (described in section 1 of Supplement 1) based on parallel-aligned liquid crystal retarders [11,12,15,25].
The inspected tissues were obtained from 25 different ex-vivo specimens of chicken thighs provided by the same commercial brand. Tissues were taken from unambiguous anatomic locations under the supervision of a pathologist: we disinserted the muscles in order to dissect skeletal muscle, tendon and myotendinous junction as well as bone and worked with the same decomposition conditions. To do so, two samples per each tissue type were obtained (a total of 23 bones, 50 muscles, 50 myotendinous junctions and 34 tendons) and subsequently frozen at −16 • C. Three hours before measuring, the corresponding samples were systematically defrosted in order to equalize the decomposition stage. Also, histological sections from the same regions that were studied were collected and stained using a Masson Trichromic technique, which allows a clear-cut differentiation of all connective tissue (and otherwise mesodermal as bone and muscle). The detailed process and the histological images (Figs. S2 and S3) are provided in the section 2 of Supplement 1.
An example of an image of each one of the four studied tissues is shown in Fig. 1. Soft tissues (muscle, tendon and myotendinous junction) and bone present different properties owing to their physiologic functionality and structure. In the following, we describe a physiological analysis of each soft tissue in order to well-differentiate its structure and biological components. The skeletal muscle ( Fig. 1(a)) is composed by contractile myofibril chains bundled into fascicles sheathed by perimysium (connective tissue which contains type-I collagen fibers [44]) and packed within a collagen-rich epimysium layer, which protects muscle fibers from friction. By comparison, tendons ( Fig. 1(b)) are composed by dense collagen (60-80% type-I collagen) fibers arranged into parallel fascicles by following the same orientation as muscle fibers [45,46] which are completely enclosed by paratenon (a thick fibroadipose layer) and fascia (collagen-rich tissue that covers both muscle and tendon) [47,48]. In turn, myotendinous junction ( Fig. 1(c)) is a variable combination of both previously mentioned soft tissues, which are progressively mixed along the transition muscle-tendon: contractile myofibrils and collagen fibers are bundled into separated fascicles [49] covered by fasciae. Each fascicle contains fibers of single type (either contractile or collagenous) but fascicles of each type are intermingled in varying proportions. Finally, bones ( Fig. 1(d)) are composed by a dense matrix of collagen fibers with a varying arrangement that depends on the type of bony tissue examined. The cortical bony tissue of a long bone like the one here examined is formed by parallel collagen fibers arranged following the main axis of the bone and interrupted by lacunae (spaces containing cells) at more or less regular distances. Around lacunae collagen fibers are distributed concentrically. Mineral deposits are located in the collagen fibers and, particularly, in the spaces between adjacent collagen fibers. Obtained experimental Mueller images correspond to a region of interest (ROI) of 512 × 512 pixels which correlates to an area of 1.1 × 1.1 cm 2 (the images of the experimental Mueller matrices, MMs, as well as MMs of the decomposed parameters according to section 2 are presented on Figs. S4 to S11 on Section 3 of Supplement 1). The polarimetric analysis is performed for a pure tissue region, which exclusively contains the pixels related to the desired type of tissue and from which the background has been mostly removed. By using the measured experimental Mueller images, polarimetric images corresponding to different polarimetric metrics were calculated, according to the equations provided in section 2. To suppress the possible specular reflection effects on camera and to perfectly isolate the desired type of tissue, we selected, from the original 512 × 512 pixels polarimetric images, a specific sub-image with a sub-ROI of 150 × 150 pixels (which does not include the saturated -specular reflection originated-regions) from which we compute the mean value and standard error for the different polarimetric metric sub-images. The obtained standard errors of the mean for all the metrics calculated are within a value ranging between 0.006% and 0.02%, except for the optical rotation Ψ , whose highest standard errors are ranging from 0.3% to 3.4%: those optical rotation errors could be related with the heterogeneity of the tissue fibers orientations through the whole sample. The general low values of the standard errors suggest that the computed mean values are reliable and consistent, probably due to the large number of pixels (150 × 150) measured per sample.

Statistical treatment
The object of the study is a data matrix composed of 157 tissues cases, analysed through 27 predictor variables (9 polarimetric metrics measured at three wavelength channels; 625 nm, 530 nm, and 470 nm) and one classifier variable with four categories (the type of biological tissues: muscle, tendon, myotendinous junction and bone). The 9 selected polarimetric metrics, calculated from the experimental Mueller matrix images according to section 2, are: P ∆ , P 1 , P 2 , P 3 , P, R, D, δ and Ψ . Under this scenario, the whole data matrix corresponds to 157 x (27 + 1) items. It is important to note that there are not conditional relations when experimentally measuring the 157 cases, so they can be considered as independent observations. In the following, two statistical studies we applied on measured polarimetric data are described: (1) an exploratory univariant statistical analysis of the polarimetric data distributions profile; and (2) a factor analysis on data matrix with principal components extraction, based on the results obtained with the exploratory analysis, to obtain a more adequate basis to construct a logistic model capable to discriminate between the four studied tissue categories. These two studies are provided in subsections 3.2.1 and 3.2.2, respectively. Note that the different statistical tests, following discussed, applied on experimental data were conducted by using the SPSS software.

Exploratory univariant statistical analysis on polarimetric data distributions
The exploratory univariant statistical analysis applied on the data matrix above-mentioned was performed at a global level and according to classificatory variable groups (the type of tissues). Particularly, data distribution is analyzed by studying the central tendency, the dispersion and asymmetry behavior, outlier's identification and graphical differences between samples sorted by type of tissue characterization from Boxplot [50,51], as well as by checking normal distribution fit assumptions (assess of normality charts and contrast tests Kolmogorov-Smirnov and Shapiro-Wilk [52][53][54][55]). The conducted analysis proved that, in general, neither these data distributions nor their transformations (according to Tuckey's ladder of powers [56]) fit the normal distribution behavior. Particularly, primary statistics show that an asymmetric data cloud cannot be fitted into a normal distribution. In this context, we observe that the 27 metrics present asymmetry: skewness values different from 0 (associated with normal distribution) indicate the mean position with respect to the median. The closest skewness value to zero is -0.111 (linear retardance, δ, distribution at 625 nm) meanwhile the highest one, 2.257, is achieved by polarizance, P, measured at 530 nm. Only 7 from the 27 measured polarimetric distributions can be associated to normal distributions, with skewness values close to zero (they present values lower than |0.840|), and they are mostly found at 625 nm. The remaining 20 polarimetric distributions present higher skewness values, being the polarimetric distributions based on the 470 nm channel those with the clearest asymmetric tendency.
Additionally, regarding the Boxplot analysis [50,51], we evaluate the total percentage of outliers (mild and extreme) to show that there is no relevant measurement bias: among the different polarimetric distributions analyzed, the outlier values range from 0% to 5.09% of the total distribution data. Because we have not done any normality assumption on measured data and because the subsequent statistical analysis is based on robust statistical techniques, we conduct a conservative approach and those detected outlier values are not removed from the study.
In Fig. 2, we present, as a representative example to illustrate the small number of outliers associated to the observed data distributions, the Boxplot for the Index of Purity P 2 distributions measured at 625 nm illumination channel and corresponding to the four studied tissues. The few observed outliers are represented in Fig. 2 by circles (mild) and stars (extremes). What is more, the asymmetric tendency of the distributions is suggested by the relative position of the median, which is not equidistant from the first and third quartiles (they are not placed at the middle of boxes). In addition, the tissue discrimination capability of the metrics is also suggested by the boxplot analysis as medians do not fit within the boxes of other tissues. For instance, Fig. 2 illustrates the potential of P 2 to discriminate muscle among the remaining tissues (see dashed red line). Similar tendency is obtained when analyzing 530 nm and 470 nm illumination wavelength (represented in Supplement 1, Figs. S12 and S13, respectively). These results can be extrapolated to the remaining M-metrics studied in this work. The normality analysis above conducted on polarimetric distributions proved that subsequent data analysis must be done by following non-parametric statistics. Under this scenario, the different polarimetric distributions, sorted by type of tissue, were studied by performing a non-parametric homogeneity multiple contrast analysis: the Kruskal-Wallis test [57,58]. This test allows us to determine if there are significant statistical differences between two or more data distributions. Therefore, in our case, the outputs are connected with the tissue discriminatory potential of the studied metrics: p-values lower than 0.05 confirm, with a 95% of confidence, the capability of the polarimetric indicator to discriminate the origin of the data between different pair of tissues.
Concerning the Kruskal-Wallis test results (whole data is presented in Table S1 of Supplement 1) we confirm that P ∆ , IPPs (P 1 , P 2 and P 3 ) and polarizance P, are the most sensitive and interesting polarimetric indicators for the studied tissue characterization, while the optical rotation Ψ is not capable to differentiate between any pair of tissues. Particularly, polarizance P seems to be the only metric capable to discriminate between myotendinous junction and bone pair. With regards to retardance indicators, they also provide some predictive potential. From this, we conclude that the main polarimetric characteristics of samples (selective absorption, retardance and depolarization) present certain signatures in the studied tissues, and therefore, the ensemble of polarimetric indicators selected seems to be adequate to be used in optical models to discriminate between the studied tissues.
The univariant analysis results above-discussed, particularly the non-normal tendency of studied polarimetric distribution, manifest the impossibility to implement our solution based on techniques such as discriminant analysis: this requires, in addition to the assumption of equality of original metrics, sorted by type of tissue, variance-covariance matrices, that these metrics follow a multivariant normal distribution behavior. As a consequence, our proposal considers an elaborated model based on a more robust technique, the logistic regression, which does not require these assumptions. At the same time, a dimension reduction of the original matrix in pursuit of non-correlated predictors (that minimize or delete the multicollinearity that generates instability on models) is aimed. Particularly, this solution considers: (1) A factor analysis with principal components extraction [59][60][61][62][63][64][65][66][67][68][69][70][71][72][73] that provides a set of independent predictor variables (defined as principal components) that maintains the information contained in the original variables and describes an acceptable proportion of the variance or inertia of point clouds.
(2) The implementation of multivariant predictive models, based on binary logistic regression [74] techniques with the extracted principal components as parameters, for the type of biological tissue classification, as well as the independent predictor and implemented model associated to ROC curves [75][76][77][78] draft for their predictive capacity evaluation, comparison and optimization.
The factor analysis is described in subsection 3.2.2, whereas the implementation of the statistical models and their evaluation are described in section 4.

Factor analysis: principal components extraction
The goal of this work is to implement multivariant probabilistic models, particularly a binary logistic regression, with the capability of determining whether a certain sample corresponds to any of the chicken tissue categories studied in this work: muscle, bone, tendon and connective tissue. Because the application of logistic regression on data does not necessarily converge when using a non-orthogonal basis, it is convenient to obtain a set of uncorrelated (orthogonal) predictors for the regression formulation which, in turn, define a reduced-dimension space. To do so, when performing factor analysis, we choose the principal components extraction as it always provides a solution in which the output factors compound a set of independent (uncorrelated) variables obtained by linear combinations of the original indicators. The extracted principal components maintain the same information of the original variables while keeping the maximum possible variance: they constitute the ideal option to implement on predictive models. It is important to remark that, because of the uncorrelation of these new components, the multicollinearity problem in a regression model is avoided.
Before conducting principal components analysis (PCA) on our particular polarimetric data matrix, we have previously computed the Bartlett's test of Sphericity [79] and the Kaiser-Meyer-Olkin Measure of Sampling Adequacy (KMO) [80] to ensure the suitability of factor analysis for the particular studied polarimetric data. Particularly, Bartlett's test of Sphericity is based on the null hypothesis according to which the population correlation matrix is equal to the identity matrix: in such case, the factor model is not appropriate to treat the input data [79]. On the other hand, KMO allows evaluating the degree in which each variable is predictable from the remaining ones: obtained values larger than 0.5 indicate that at least one common factor underlying the observed variables exists [80]. The output of both tests shows the suitability of applying factor analysis (PCA) to our studied data matrix. Bartlett's test clearly rejects the null hypothesis (p-value 0.000). Likewise, Kaiser-Meyer-Olkin outputs a value of 0.650, larger than the typical cut-off, meaning that there's a common subjacent factor for the studied variables. Because a correlation between the variables exists, the factor model is an appropriate methodology for treating the data. Hence, considering the output of both Bartlett's and KMO tests, they all suggest that PCA may suit our particular data set when being applied.
Once the Principal Components have been calculated by using the SPSS software on the data matrix, we have computed the scree plot [71] (Fig. 3) which allows us to select a reasonable number of components to be used for the predictive models. In particular, we decided to hold all components preceding the sedimentation zone: 10 components explain more than 90% of the original metrics variance. The eigenvalues of those components, together with the data variance (in %) and the cumulative variance (in %) are shown in Table S2 of the Supplement 1. In this way, PCA reduces the 27-dimension space of M-metrics into a new 10-dimension space by identifying the original variables set underlying dimensions. The extracted factors cannot be directly measured but they enable the data structuration around a reduced number of variability axis. As the new factors are mathematical constructions, they do not necessarily have to be connected to physical properties, but they can be further interpreted by analyzing their dependence with polarimetric metrics. From a geometrical point of view, the principal components correspond to axes perpendicular to each other that better fit with the point cloud that sets the data matrix.

Predictive model construction: results and discussion
The selected principal components encode the polarimetric information so we can write each component as the linear combination of the 27 studied polarimetric metrics weighed by some constants provided by the so-called component score coefficient matrix, shown in Table S3 of Supplement 1. The fact that the component score coefficient matrix allows expressing the different polarimetric metrics as a linear combination of the principal components, or vice versa, is determinant because this situation defines the basis for building a probabilistic model. In a general way, the principal components (written asˆ︁ C) are obtained by multiplying the M-metrics 27×1-dimension vector,ˆ︁ P, by the transposed component score matrix,ˆ︂ CS (a 27×10-dimension data matrix in Table S3):ˆ︁ whereˆ︁ C elements correspond to C i (i = 1, . . . , 10) andˆ︁ P is written as follows,︁ where R, G and B indicate the wavelength channel (625 nm, 530 nm and 470 nm, respectively) of the polarimetric indicator and the order of the metric is conserved along the vector.
To illustrate the correspondence between the principal components and the physical variables, we represent a two-dimensional space plot where the 27 indicators are scattered as a function of the two main principal components, C 1 and C 2 . The obtained results are shown in Fig. 4.   Fig. 4. The plot showing principal components C 1 against C 2 represents the correlation coefficients between the 27 polarimetric indicators and the two first extracted principal components. The notation of the polarimetric indicators names is composed of the word that represents the measured parameter: P1, P2, P3 and PA (IPPs and P ∆ , respectively), D and P (diattenuation and polarizance), R and Delta (global and linear δ retardance) and Phi (optical rotation Ψ ), followed by R, G or B (corresponding to red, green and blue measured wavelength, respectively).
Note that the metrics are clustered in four main areas (highlighted by black circles) and more interestingly, those clusters mainly coincide with their physical characteristics regardless its measured wavelength: retardance related indicators (total and linear retardance, R and δ respectively), depolarization (IPPs and P ∆ ), dichroism (diattenuation and polarizance, D and P, respectively) and optical rotation Ψ . By comparing the weight or influence of the two first principal components when analyzing polarimetric indicators we find that IPPs and P ∆ (scattering, depolarization) are mostly described by C 1 as they are clustered with weights around 0.5 and 1.0 meanwhile the weights of C 2 are between 0 and 0.5. When looking at how retardance indicators are grouped, we find that the presence of the C 2 principal component (between 0.5 and 1.0) is higher than the influence on the C 1 component (for the last one, its weight value is close to 0). For diattenuation and polarizance, they are both influenced by the two main principal components but they are especially represented by the C 1 component as most dispersion is observed in such a component. Finally, optical rotation principal components C 1 and C 2 values are both near 0 so they do not provide enough information to describe Ψ (the two first principal components are not able to explain optical rotation so we expect to find weights different from zero when analyzing the remaining components).
Pointing out some interesting facts: plotting principal components against each other provides us a very intuitive way to understand in which way those components are connected with actual physical information of samples (the 27 polarimetric indicators), and which one of those polarimetric indicators (or collection of indicators) carry more data variability in the original raw data matrix. Interestingly, the fact that physical variables of the same type (for instance retardance) are grouped in well-defined areas of the C 1 -C 2 space (or of other combinations of 2-dimensional C i spaces) indicates that different origin in data variability (orthogonal dimensions) is well connected with actual different physical structures present in the studied tissue samples.
The discussion of the principal component's space characteristics (such as its discriminatory potential, sensitivity and specificity) and the predictive model construction of each type of tissue and its predictive features are described in subsections 4.1 and 4.2, respectively.

Discriminatory potential of the principal components space
In this section, we are interested in studying the the performance of 10-principal components space as tissue classifiers. With this aim, the Receiver Operating Characteristic (ROC) curve analysis [75][76][77][78] is carried out: it describes the performance of a classifier (an algorithm or a particular variable) when classifying measures into two categories by plotting the true positive rate (TPR), or sensitivity, against the false positive rate (FPR), or 1-specificity, for multiple threshold values of the classifier. For each classifier and measure, the combination of what the sample actually is (real value, positive or negative) and how it is classified by the model (prediction; positive or negative) results in four possible outcomes [76,77]. The mathematical relations connecting these four criteria define the sensitivity and specificity, which are the base of ROC curves, and are given by: where TP, P, TN and N denote for true positives, total positives, true negatives and total negatives, respectively. Another way to define sensitivity is the ability to correctly identify positive samples while specificity results in the opposite case, correctly identifying negative samples [78]. By using the SPSS software, we computed the ROC curve for all the 10-principal components and for the four different studied tissues (a specific dichotomic predictive model is constructed for each studied tissue, so four models are constructed). As an example, the ROC curve of the principal component C 1 for each biological tissue is shown in Fig. 5. Note that the different threshold are all the possible values in C 1 . However, all principal components present statistical significance in ROC curves for one or more tissues, which bolsters the interest in keeping them in the process.
Usually, to compare the performance of classifiers (in our case, to estimate which component provides better sensitivity-specificity values for a particular tissue), the area under the curve (AUC) is calculated, with values ranging from 0, when the variable has no predictive capability, to 1, with 100% both sensitivity and specificity. The AUC values of all principal components when discriminating each biological tissue against the remaining three are summarized in Table 1.
With regards to the exposed information, we see how most of the principal components provide certain discriminative potential, with values larger than 0.5. The components with the largest AUC for a particular tissue are highlighted in Table 1. Notice that C 1 present significant  discrimination values in all tissues, C 2 and C 3 mostly provide discriminatory information for tendon, C 4 and C 5 focus on muscle and tendon discrimination, C 6 on tendon and bone, C 7 on bone, C 8 keeps the information about muscle, tendon and bone and, finally, C 9 and C 10 provide discrimination for the myotendinous junction.
Interestingly, it has been noticed that the principal components predictive potential (given by the large area under the curve, AUC, values exposed in Table 1) is different for distinct tissues. Such behavior was expected as every component is written as a linear combination of the polarimetric indicators (rows in Table S3 in Supplement 1), and the metric values strongly depend on the composition characteristics of each biological tissue. Therefore, considering the structural differences between them, some polarimetric indicators and, consequently, some principal components, are expected to have a higher capacity to discriminate a particular tissue between the rest.

Predictive model construction
In this section, the determination of the predictive model for tissue recognition is described. The constructed model is based on a probabilistic function achieved by fitting the experimental data (array of principal components values for all chicken samples) to a sigmoid function, more specifically, to a logistic function. In other words, we perform a non-linear data model transformation into a linear one by means of a logistic regression fit of data. The logistic function with the principal components as variables and with a curve maximum value equal to 1 is written as follows, where C i are the i principal components, β i the weights of the components and β 0 a constant value. The function in Eq. (13) is bounded between 0 and 1 and can be interpreted as the probability of a given outcome.
We have chosen a logistic function because it is valid for non-parametric data and it does not require the relation between the predictors and the probability of a target outcome to be linear. Moreover, despite other fitting methods, as it is the case of those based on the ordinary least square regression (OLS), the logistic function neither requires the residuals to be normally distributed and exhibit constant variance. Particularly, we have previously verified the non-dependency of the model predictors (principal components) and their connection with the dichotomic dependent variable. Moreover, by performing the PCA, the multicollinearity between predictors have been removed so stable predictive models and convergent iteration processes are achieved. Under this scenario, four logistic regression functions are designed: one model for each studied biological type of tissue. As seen in Eq. (13) the obtained probability function depends on some of the principal components factors which, in turn, depend on the polarimetric indicator values measured from the image sample. The logistic regression fit on the experimental principal components (obtained from the experimental polarimetric indicators according to the relations provided in the component score coefficient matrix in Table S3 in Supplement 1) is conducted by using the SPSS software. By applying a stepwise regression approach (backward elimination) and using the Wald estimator throughout multiple steps, the principal components are removed until only the most significant ones remain [74].
This routine gave us the values for the weights in Eq. (13) and, accordingly, the obtained probabilistic functions for the four studied tissues are exhibit below: To illustrate the goodness-of-fit for the four regressions, the Hosmer-Lemeshow significance [79,80] and the R 2 of Nagelkerke [81,82] indicators are provided in Table S4 of Supplement 1.
Afterwards, to study the efficiency of the probabilistic models in terms of sensitivity and specificity, the associated ROC curves were also computed for each one of the four models. They are provided in Fig. 6, in which their corresponding AUC are also indicated. We see how the AUC are significantly larger than those we obtained when representing the principal components (see Fig. 5 as an example), so we see how the discriminatory potential of the probabilistic functions clearly overcome the obtained with the principal components by themselves. We want to highlight that the AUC values obtained for muscle, tendon and bone tissues models, are significantly high (0.92, 0.95 and 0.89 respectively), providing the efficiency of these models for tissue discrimination. In turn, the model providing lower discriminatory potential is that associated to the myotendinous junction (AUC=0.79). This situation agrees with the fact that only C 1 , C 9 and C 10 components present discriminatory potential for the myotendinous junction tissue (see Table 1), and the two latter, represent a very small data variation of the whole data matrix (2.350% and 1.885% of the variance (Table S2), respectively). Note that models in Eqs. (14)(15)(16)(17) provide as output a real number between 0 and 1, associated to the probability of a given tissue to belong to a particular tissue category. To build a dichotomic model determining if a tissue belongs or not to a particular category, a threshold (specific probability value) must be set in a way that values above/below such threshold are associated with Yes/No answers (i.e., is the unknown sample a bone? Yes-No; is it a tendon? Yes-No, etc.). However, each possible threshold (probability), as can be seen from ROC curves in Fig. 6, will lead to a particular pair of sensitivity and specificity values for the constructed models. In our work, the criterion used to set such threshold is by using the Youden's Index [83] of the ROC curves in Fig. 6. According to this criterion, the optimal cut-off for our models corresponds to the farthest point of the probabilistic function in the ROC curves (blue curve in Fig. 6), from the ROC diagonal (red line in Fig. 6) along Y axis, which is calculated as the maximum value, d max , of d = sensitivity + specificity -1.
The computed Youden's index [83] for each model, and its corresponding sensitivity and specificity values, are shown in Table 2. Note that the models for the muscle, tendon and bone present high values of specificity and sensitivity (values higher than 80% in all the cases). In turn, the worst obtained result is that for the myotendinous junction tissue specificity, which descends to a value of 71% (in agreement with the limited potential of the myotendinous model discussed in the previous analysis). Therefore, the sensitivity and specificity obtained results for the constructed predictive models (data in Table 2) quantify the efficiency of those models to categorize the studied organic tissues. What is more, the current study highlights the potential of predictive models based on polarimetric data to discriminate between animal tissues, and the methods provided could be applied in the future, for instance, for medical applications.
Following this same idea, to show the potential of the method to be applied to imaging techniques to discriminate between tissues, we encoded an algorithm that computes the probabilistic function for each pixel of the measured sample. In particular, for an arbitrary chicken tissue, the Mueller matrix image is measured at the three studied wavelengths (625 nm, 530 nm, 470 nm), from which the 27 polarimetric images are calculated. From this polarimetric image database, the probability image of the arbitrary tissue to be categorized as a particular category (muscle, bone, tendon and myotendinous junction) is calculated according to the logistic probability functions in Eqs. (14)(15)(16)(17). The algorithm outputs which is the probability of the analyzed pixel to be recognized as a particular tissue, so the probability image is constructed, with values ranging between 0 and 1.
As an example, in Fig. 7 we provide the probability function images, corresponding to the four predictive models when analyzing an arbitrary sample of a chicken tendon. Figure 7 presents the intensity image of the tendon (M00, Fig. 7(a)), and the output images of the probabilistic models for the recognition of muscles ( Fig. 7(b)), tendons (Fig. 7(c)), myotendinous junctions ( Fig. 7(d)) and bones (Fig. 7(e)). Probability images are given in grayscale, in which white is associated to the maximum probability value (1) and black to the minimum (0). By comparison, the highest probabilistic values (closer to 1) are clearly obtained by the tendon predictive model (Fig. 7(c)), recognizing almost all the pixels along the tendon tissue image, and thus, perfectly classifying the sample as a tendon. In turn, the probability images for muscle ( Fig. 7(b)), myotendinous junction ( Fig. 7(d)), and bone ( Fig. 7(e)) are based, despite some clustered high output pixels, on significantly smaller probabilistic values than those present in tendon recognition image (Fig. 7(c)). Therefore, those models do not assign the tendon sample as muscle, myotendinous junction or bone, respectively. Among these three probability images, the larger values are obtained when using the myotendinous junction model, providing that this model prediction is the less sensitive and specific one, in agreement with previous discussions on model efficiency.
Analogously, we provide a second example, for the study of an arbitrary muscle sample. Results are provided in Fig. 8, for the intensity image (M00, Fig. 8(a)), and the four predictive model output images (muscle, tendon, myotendinous junction and bone recognition, Figs. 8(b-e), respectively). The obtained predictive results demonstrate the high capability of the musclerecognition model: the highest overall values for the output probability image correspond to the muscle recognition model (Fig. 8(b)), perfectly classifying the muscle tissue. When analyzing the probability distribution for the three remaining models, a good discriminating potential is also demonstrated: the sample's output shows low probability values for classifying the muscle as a tendon (Fig. 8(c)), a myotendinous junction (Fig. 8(d)) or a bone (Fig. 8(e)). Once again, and as expected by the previous analysis (sensitivity and specificity model values in Table 2), the worst predictive efficiency is related to the myotendinous junction model, leading to larger probability values than the tendon and bone models.
Note that the tendon and muscle study cases (results in Figs. 7 and 8) highlight the strong qualitative and quantitative discriminative potential of the models to be applied for tissue recognition imaging techniques. For the sake of completeness, the remaining tissue study cases, bone and myotendinous junction recognition test, can be found in section 6 of Supplement 1. These examples (Fig. 7, Fig. 8 and Fig. S14 and Fig. S15 of Supplement 1) point in the direction that the proposed models are suitable for tissue classification and recognition, and the corresponding sensitivity and specificity outcomes will be limited by data provided in Table 2. , tendon (c), myotendinous junction (d) and bone (e) for chicken muscle measurements. The gray level bars, placed to the right of the corresponding probability function images, defines the probability of the pixel to be recognized as a particular tissue, in a range between one (white) or zero (black).

Conclusions
Four predictive models based on the measure of diverse polarimetric metrics derived from the experimental MM of four ex-vivo chicken tissues (bone, tendon, muscle and myotendinous junction tissue) have been designed for tissue recognition and provided in this manuscript. The predictive models start by studying the distribution of the experimental data and their different statistical origins. Obtained results proved that we deal with non-parametric data with significant potential to discriminate between studied tissues. Afterwards, with the idea of constructing a robust non-parametric predictive probabilistic model, we applied a Principal Component Analysis on data, which allowed us to reduce the dimension of the data space (from 27-dimension to 10-dimensions) but dealing with maxima information as well as to obtain an orthogonal basis of parameters. By analyzing the connection between the principal components basis and the original polarimetric metrics, we realized that different variability directions of data in the principal components space can be roughly associated by different physical origins of the tissue structures (dichroism, retardance and depolarization).
Subsequently, based on the computed principal components basis, we constructed a logistic regression on data, leading to four probabilistic models (a model for each one of the four studied tissues) able to categorize if an arbitrary tissue belongs to a certain tissue category. The sensitivity and specificity values reached for each one of the four models are respectively: 85.3% and 93.5% (Tendon-model), 86.0% and 88.8% (Muscle-model), 82.6% and 80.6% (Bone-model) and 82.0% and 71.0% (Myotendinous junction-model), proving the potential of the method.
Finally, to highlight the suitability of the provided methods to be applied in biological samples imaging recognition techniques, we computed an algorithm, based on the constructed probabilistic models, that provides the probability, between 0 to 1, of a given pixel to be categorized as a particular tissue, and this is done for all the pixels of a given tissue image. In short, the procedure outputs four different probability images, corresponding to the probability of an arbitrary tissue to be muscle, bone, tendon or myotendinous junction tissue. The method provided complementary visual interpretation for tissue recognition, and a satisfactory categorization of the analyzed tissue image, according to the sensitivities and specificities stated above.
The proposed non-invasive methods discussed in this work could be applied, by conducting the required statistical data feeding, in multiple biomedical scenarios, for example, for the early diagnosis of some pathologies.