Multispectral Imaging Method for Rapid Identification and Analysis of Paraffin-Embedded Pathological Tissues

The study of the interaction between light and biological tissue is of great help in the identification of diseases as well as structural alterations in tissues. In the present study, we have developed a tissue diagnostic technique by using multispectral imaging in the visible spectrum combined with principal component analysis (PCA). We used information from the propagation of light through paraffin-embedded tissues to assess differences in the eye tissues of control mouse embryos compared to mouse embryos whose mothers were deprived of folic acid (FA), a crucial vitamin necessary for the growth and development of the fetus. After acquiring the endmembers from the multispectral images, spectral unmixing was used to identify the abundances of those endmembers in each pixel. For each acquired image, the final analysis was performed by performing a pixel-by-pixel and wavelength-by-wavelength absorbance calculation. Non-negative least squares (NNLS) were used in this research. The abundance maps obtained for the first endmember revealed vascular alterations (vitreous and choroid) in the embryos with maternal FA deficiency. However, the abundance maps obtained for the third endmember showed alterations in the texture of some tissues such as the lens and retina. Results indicated that multispectral imaging applied to paraffin-embedded tissues enhanced tissue visualization. Using this method, first, it can be seen tissue damage location and then decide what kind of biological techniques to apply.


Introduction
Multi-and hyperspectral imaging are non-destructive techniques that integrate both spectroscopic and imaging techniques or computer vision in one system [1,2]. These methods are used in numerous medical applications, especially in disease diagnosis and image-guided surgery [3][4][5]. It is often applied to obtain quantitative information for several biological processes in both healthy and diseased tissues, e.g., in tumor detection and diagnosis, cancer type differentiation and cancer cell metabolism [6][7][8], and cardiovascular conditions [9,10], to quantify metabolic activity in cerebral tissue [11,12] and to measure the changes in the spatial distribution of tissue oxygenation during vascular occlusion [13], among many other applications.
Spectral imaging is based on constructing an image cube which consists of multiple slices of the same scene acquired at a narrow band of the electromagnetic spectrum [2]. The final data acquired by the construction are called "hypercube" since they can be illustrated as containing 3 dimensions with two for spatial coordinates (x, y) collected in the X-Y plane and the other one for spectral information λ represented in the Z-direction [11]. Photon propagation through disordered media such as biological tissues can undergo both scattering (from inhomogeneity of biological structures like extracellular matrix, mitochondria, cell nuclei) and absorption by chromophores such as water, lipids, hemoglobin, or melanin [14]. The relative probability of these processes in a given tissue depends on wavelengths [15].
Folic acid (FA) (vitamin B9) is an essential micronutrient for the growth and development of embryos. Maternal FA deficiency is a serious public health problem [16,17] and has been implicated in a number of fetal abnormalities, especially, but not exclusively, the nervous system [18]. It is known that during the progression of a disease, the scattering and absorption properties of tissue change [2]. For this reason, in the present work, we aimed to evaluate alterations in the absorbance behavior of embryonic eye tissues due to a maternal folic acid-deficient (FAD) diet by using multispectral analysis.
In previous studies, we have observed that maternal FAD diet alters protein expression patterns in the whole embryonic eye [19,20], as well as disorganizes fibers/cell orientation [21]. To observe these alterations, before applying image processing, tissues were stained by antibodies to recognize collagen IV laminin-1 in tissue. However, we know that maternal FAD could alter many molecules and tissues of the eye. Traditionally, to analyze pathological tissues, different biological techniques are usually used, among which those that use antibodies for immunocytochemical methods, as well as in situ hybridization experiments that allow for the localization of RNA or DNA targets in cells and tissues. We describe a new imaging technique based on multispectral image analysis applied to paraffin-embedded tissues for noninvasive biomedical diagnosis. The goal of our study was to detect structural alterations of eye tissues as a diagnostic aid before applying other techniques.

Animals and Diet
All mice were housed under temperature-and light-controlled conditions at 22 ± 2 °C on a 12:12 light-dark cycle with free access to food and water. Experiments were carried out in 8-week-old C57/BL/6 J female mice. They were assigned in equal numbers to 3 groups according to the diet: (1) control diet (a standard FA rodent diet, 2 mg/kg diet), (2) FAD diet for 2 weeks (D2 group), and (3) FAD diet for 8 weeks (D8 group). At the conclusion of the experiment, mice were anaesthetized by carbon monoxide and sacrificed by cervical dislocation at 14.5 days of gestation (E14.5). Embryos were removed by cesarean section, placed in cold sterile phosphate-buffered saline, fixed for 48 h in buffered formaldehyde, and decapitated. The heads were then embedded in paraffin following standard procedures. Finally, frontal sections were made with a thickness of 5 µm and placed on slides. Finally, to observe a green autofluorescence without the use of any labels or stains, sections were visualized utilizing a Zeiss Axioplan 2 imaging microscope equipped with a FITC filter. In this study, we studied nine mouse embryos in each group, and one eye per individual was used in the analysis.
The experimental protocol was reviewed and approved by the Animal Experimentation Committee of the Complutense University Madrid (UCM). Mice were maintained in the Animal Facility at the Faculty of Medicine, UCM. Figure 1A shows the experimental setup used for multispectral image acquisition. In general, the system is made up of the following elements:

Description of the Experiment
Monochromatic illumination for multispectral measurements was provided by a tunable light source TLS-25. Multispectral images were obtained in the spectral range 390-705 nm with a resolution of 25 nm. As shown in Fig. 1B, the spectra have a small bandwidth; besides, the intensity level is different for the different wavelengths. The incident beam passes first through an achromatic doublet lens (L1), to collimate the monochromatic light beam. The emerging light passes through tissue samples; then, light emerging from the tissue is collected and focused on two lenses (L2, f = 50 mm; L3, f = 400 mm). The acquisition of multispectral images was done using a 1/4″ Sony CCD progressive scan FireWire CCD Color Camera (DFK21AF04); the CCD chip produces a resolution of 640 × 480 pixels per frame. Finally, CCD camera was operated via a FireWire connection to a desktop PC. The software, IC Capture, was used to capture and display image sequences with a frame rate set at 50 frames per second. Once the sample has been focused on, 50 consecutive frames were taken at a certain selected wavelength 0 .
Taking into account the spectral sensitivity of the camera, the multispectral images were taken from 390 to 705 nm in 15-nm wavelength steps. The 50 frames were averaged pixel-by-pixel to minimize the temporal effects of noise in CCD detectors. Images were acquired in gray scales from 0 to 256 levels using the Matlab Image Acquisition tool. In general, all the images have been taken in darkness. First, for each wavelength, the image has been taken without the sample and then with it. Therefore, for each wavelength, we have divided the image captured without the sample by the image collected with the sample; in this way, the final image recorded is the transmission of the sample. Finally, we have for each of the 27 eyes analyzed (9 eyes in each sample group) a set of multispectral images T Gi (x, y, 0 ) , where G is the group (control, D2, and D8), i the analyzed eye (i = 1,…, 9), x, y the pixel in the image, and 0 is the selected wavelength. The acquired dataset is called spectral data cube (Fig. 1D). The spectral cube combines spatial and spectral information. In total, all captured images form a set of 594 multispectral images. As can be seen, it is a huge amount of information; so, to analyze it, we apply Principal Component Analysis (PCA) method (described below) to reduce the dimensionality of the data.

Principal Component Analysis Method
Our multispectral images have been analyzed using a multivariate statistical method known as principal component analysis (PCA). It is the best, in the sense of least-square error, linear dimension reduction technique [22,23]. It reduces the dimensions of a dataset with minimal loss of information. In essence, PCA seeks to reduce the dimension of the data by performing an orthogonal transformation to the basis of correlation eigenvectors and projects onto the subspace spanned by those eigenvectors corresponding to the largest eigenvalues [24]. In the case of multispectral imaging, PCA re-expresses the original spectra bands by combining their orthogonal linear components with the largest variance. It computes an orthonormal basis for a series of images, transforming the spectra cube to a "noise-less" structure with a lower dimension [14].
We name F to a set of frames that can be written as F = F 1 , F 2 , ..., F , ..., F N , where N is the number of frames chosen registered with a determined spatial frequency. If the detector array camera has R rows and C columns, it is possible to organize the M = R × C signals of the pixels as a two-dimensional column vector; we can so characterize the position of the pixels, considering its row r and its column c, as k = (c − 1)R + r , where r is from 1 to R and c is from 1 to C. Therefore, k takes values from 1 to M, where M is the total number of pixels in the array. Then, a set of N frames can be written as an F matrix of M rows and N columns; the matrix F is defined as: where p M ( ) is the signal of the pixel M in the space λ. Linear combinations of the original frames form eigenvectors E of the covariance matrix S of the original frames. These For all this, the expansion of the PCs is obtained as a M × N matrix, with the following expression: Each PC, established by columns, is given by: In this way, each PC is an image that is properly combined to produce the original set of frames. Using the eigenvalues to form a new basis, any pixel vector can be expressed as a linear combination of the eigenvectors: This leads us to the following equation: The element e , is the weight of the PC in the frame F .

Calculation of the Absorbance or Optical Density
For a given wavelength 0 , the absorbance of a given sample can be defined as: where T Gi x, y, 0 is the transmission of the sample at a defined wavelength; it can be defined as follows: In this way, the absorbance defined in (7) will be written as: where j is the total absorption coefficient of a given component j present in the sample, j (x, y) is the concentration of j at the image point (x, y), and L is the thickness of the sample that, in our case, is 5 µm.
If now we apply a decomposition of principal components onto absorbance data, where the average value of each absorbance image ⟨A Gi ( 0 )⟩ x,y was removed to obtain images with zero mean, so we would have: This decomposition interprets each of the images as a mixture of original images uncorrelated between them. These images are the PCs. Since the different images were taken at different wavelengths, the eigenvectors e are those that carry the trend in 0 and explain how we must blend the principal components to obtain the initial images, while the PCs carry the spatial dependence. Then, if we compare (9) with (10), we see that eigenvectors are related to the absorption coefficients j and the PCs to the spatial distribution of said component. The eigenvalues represent the total variance of each PC. Figure 2 shows the original images located in a configuration space, where each direction symbolizes a certain λ. In bold, there are three axes associated with λ1, λ2, and λ3 taken, which are therefore associated with the directions (1, 0…0), (0, 1…0), and (0, 0, 1… 0), respectively. Next to these axes appear the absorbance images taken for each of these wavelengths. The images have a non-zero covariance since they show great similarity between them. The decomposition of the CPs can be seen as a search for new directions in this space (eigenvectors), such that when the original images are mixed, taking the components of the eigenvectors as coefficients of the linear combination, new images are obtained (PCs) that show a null correlation between them. These new principal directions (eigenvectors) appear as dotted lines in the figure. The eigenvector and PC associated with that direction are represented next to them. Lastly, the eigenvalue would represent the variance of the PC associated with that direction.
The analysis of the PCs involves decomposing the absorbance images according to the following decomposition: where x is a spatial coordinate within the image, 0 is the wavelength, e 0 is the eigenvector, and PC (x) is the image corresponding to the PC α .
This decomposition makes all Pcs have a zero mean as well as can contain both positive and negative values. Likewise, eigenvectors are difficult to interpret directly as spectra, since they also have positive and negative values that reflect values above or below the mean absorbance, A( ) . For this reason, the need has arisen in the hyperspectral image analysis environment to carry out a decomposition analogous to that described by Eq. (11), but in the following way: where E k ( ) is called "endmember k" generally defined as a spectrally unique, idealized, and pure signature spectrum weighted by their correspondent fractions, or abundances Ab k (x) . The pure spectral signature signifies the complete reflectance of a pixel exclusively occupied by single surface material. These pure signatures can decompose the hyperspectral scene into abundance fractions by means of spectral unmixing algorithms [25]. Many endmember extraction algorithms have been developed including pixel purity index [26], N-FINDR algorithm [27], and the orthogonal subspace projection (OSP) algorithm [28].
In the current study, Eq. (12) involves solving not a single system of equations, but one for each point x in the image, with the condition that the abundances must be numbers greater than zero. This non-negative least square (NNLS) algorithm is subject to all the solutions being greater than or equal to zero. Mathematically, obtaining the abundances is equivalent to solving the system Ux = b , for each point x, where b is the regression vector of the least squares problem; hence, x would be calculated as x = U −1 b . However, the U-matrix, in our case the collection of endmembers k, is generally not a square matrix (the number of wavelengths is usually greater than the number of spectra considered relevant); hence, it only makes sense to talk about solving the system approximately by using pseudo-inverses. That is, we use the formula U T × U × x = U T × best with the following solution x = U T U −1 U T best , which is not exact since best is an approximation of b [29]. Furthermore, the solution is restricted to considering only values in which "x" takes The initial problem of calculating the abundance is to find the spectra to be considered endmember. In this study, we used the structure of the eigenvectors associated with the PCs to obtain endmembers. As can be seen in Fig. 3, the pixel values of the absorbance of three images of a control eye at three wavelengths 1 = 435nm, 2 = 540nm, 3 = 660nm have been represented in a "scatter plot" diagram. The directions associated with the eigenvectors of that "scatter plot" have also been represented; their length is proportional to data propagation in these directions. The standard deviation of this dispersion is given by √ Γ , where Γ is the eigenvalue associated with the direction α; note that eigenvectors have a unit module. Hence, if we form the vector E ( ) = A( ) + 2 √ Γ e ( ) , we are forming a spectrum that would be represented as a pixel located at the ends of the red arrows shown in Fig. 3. These new vectors continue to be related to the eigenvectors and, at the same time, represent spectra located on the "outside" of the point cloud that forms the scatter plot. The new vectors also satisfy the conditions to be considered endmembers. After obtaining the endmembers E , the NNLS algorithm was used to estimate the corresponding abundances Ab(x) for each of them as a positive number that we can directly relate to the concentration of the said spectrum in the image. Given the orthogonality of the eigenvectors from which the endmembers come, they will reflect structures in the image as disjoint as possible and ordered in order of importance in the data [31]. Table 1 shows the eigenvalues and percentages of the variance associated with each component and for each sample group (control, D2, and D8). As can be seen, in all groups studied, there are five eigenvalues, in order of importance, that accounted for greater than 99% of the total variance. On the other hand, Fig. 4 shows the mean value of the eigenvectors, from #1 to #5, together with its errors. From this figure, it can be deduced that each eigenvector is approximately the same for all eyes that comprise the same group since the error bars are small. For this reason, their mean values will be used as a representative eigenvector base of the entire sample group.

Results
The error bars most often represent the standard deviation of a dataset. We have observed that within the same individuals of the two groups D2 and D8, the standard deviations were large which indicates that there are a lot of variances in the observed data, that is, there is wide variation among individuals of the same group. As can be seen, the eigenvectors are quite similar but with a slight difference in the spectral behavior between the control and the FAD groups.
Finally, we represent in Fig. 5 the mean absorbance ⟨A Gi ( 0 )⟩ x,y for all groups. As can be seen, the wavelength most absorbed was around 420 nm.  In our case, we have found five relevant endmembers (Fig. 6A). After the endmembers have been found, NNLS algorithm was used to obtain abundance maps of all groups. The abundance map identifies the proportion of each endmember present in the spectra of each pixel. To analyze these maps, first, we study the probability density function (pdf) of the abundance maps of each endmember; the results appear in Fig. 6B. In this figure, it can be seen how the highest abundances are for endmember #1 (first component) and endmember #3 (third component).
Images of the abundance maps are detailed in Fig. 7. Abundance maps are images of the abundance values shown  Figure 7A, D, and G represent the abundance maps of the first endmember for control, D2, and D8 individuals, respectively.
In Fig. 7A (control), it can be seen how the spectrum is mainly concentrated in the vitreous and the choroid, indicating especially the areas of blood vessels (to distinguish the different structures of the eye, see Fig. 1C). However, these points of high abundance decrease in the abundance maps of both groups D2 and D8, especially in the vitreous region. The areas where there is the greatest difference are 450-500 nm, a peak at 550 nm, and a drop from 600 nm. Of all chromophores that are present in the eye and can exhibit this type of behavior is maybe the hemoglobin (Fig. 8). As can be seen, both spectra are similar.
On the other hand, the abundance maps obtained from the third endmember ( Fig. 7 (central column)) seem associated with both structures the retina and lens, since there is the presence of high abundances greater than 0.2.

Discussion
This paper presents the findings of a multispectral analysis, based on the measured spectral difference between normal embryonic eye tissues and tissues obtained from embryos with a maternal FA deficiency diet. Results of the first abundance maps suggest that eyes from D2 and D8 groups exhibited alterations in the vascular region of both vitreous and choroid. It is known that FA deficiency affects directly red blood cell metabolism and oxygen balance that they can carry to the tissues [32]. In fact, low oxygen content in blood, due to a deficit of FA, triggers a whole series of adaptive responses by the body, especially, reduction of cell proliferation and change in the size of the lumen of the blood vessels, as well as changes in the shape of red blood cells (they tend to deform tending to a spherical shape) [33,34]. Therefore, those alterations should be observed in those areas where it is normal to see blood vessels, such as vitreous and choroid. These areas are precisely those that appear signaled in the images of differences in spectral behavior between the control and the groups D2 and D8. In summary, we could link the changes detected to an imbalance in the levels of hemoglobin/oxyhemoglobin due to maternal FA deficiency. Vascular pattern alteration observed in the vitreous body and choroid could be triggered by the hypoxia induced by FA deficiency. In recent research [35], we have observed that the size and density of choroidal vessels (which lies under the retinal pigment epithelium) and hyaloid blood vessels (located within the vitreous) were increased in various embryos from both groups D2 and D8. Specially, we found large choroidal vessel caliber than those of the control group; this abnormality was the cause of retinal detachments in those embryos. In conclusion, our method was able to detect changes in blood vessels only by evaluating paraffin-embedded tissues (unstained) by using multispectral analysis.
On the other hand, the abundance maps obtained from the third endmember indicate some changes in the texture within the cornea, lens, retina, and eyelids. Maybe they are associated with microstructural alterations like a change in the spatial organization of the extracellular matrix as well as small changes in the fiber orientation. All these texture alterations mentioned were found in previous studies realized by our research group [20,21,36]. Furthermore, a strong green autofluorescence can be observed in tissues from D2 and D8 groups compared to the control. The autofluorescence in the tissue was studied without the use of any exogenous fluorophores. In another recent study [37], we also observed alterations at the lens, retina, vitreous, and choroid levels, using polarized light applied to these tissues.
The present study clearly shows that the optical characterization of tissues by multispectral imaging can be used efficiently to discriminate pathological tissues embedded in paraffin since the findings in both studies would seem to concur. It is known that absorption peaks at a specific wavelength can be related to a particular molecule. These peaks could be used as a fingerprint of the molecules' response to light, providing valuable information that can be becoming useful for diagnosis purposes [38]. Our method reveals important information about structural and/or biochemical changes in the different components of ocular structures, which indicates that this novel technique can serve as a useful tool for studying cellular or molecular alterations in tissues before the applications of molecular labeling.

Conclusions and Future Research
Multi-and hyperspectral imagery are well-known biomedical imaging techniques; both have in common that light-tissue interaction is used as a source of information since the acquired data contain spectral information about optical properties of the tissue such as absorption and scattering. We describe a novel method that combines multispectral imaging in the visible spectrum, data postprocessing, and reconstruction methods, which allows estimating tissue optical density for each point in the image at multiple wavelengths of light. Statistical analysis was computed using the PCA method, where the original multispectral dataset was transformed via PCA from its wavelength space into a new dimensional space. The main advantage of this technique is that the sample preparation is eliminated, that is, non-deparaffinized tissues should be used.
The method has been tested on normal and pathological eye tissues and showed reasonably good results. However, the authors are aware of the fact that more types of tissues are needed for testing the method adequately. Furthermore, new more extended studies should be carried out in the future to evaluate the kind of devices and parameters that can be used to further elucidate these findings before this technique is recognized as an applicable method of nondeparaffinized tissue assessment.