Introduction

Raman imaging has significant potential to provide information that not only reflects the (relative) population of various chemical species but also their spatial distribution on biological samples1,2,3,4,5,6,7,8. Its label-free and noninvasive nature, combined with the ability to perform measurements without any pretreatment, such as staining, allows pathological specimens to be preserved for subsequent histological and biomolecular analysis.9,10. Each pixel on a two-dimensional Raman image possesses a Raman spectrum, the intensity of Raman scattering by molecules in the pixel region for each Raman shift. Therefore, a Raman image has a data structure spanned by the 2D spatial position and Raman shift. Because the number of measured Raman shifts reaches several hundreds or even thousands, Raman images contain much chemical information.

Recent advances in Raman microscopy have enabled us to access rich information in subcellular, in vitro11,12, and ex vivo tissue imaging13,14,15. Because the spectral differences of biological samples are often subtle, chemometrics aids their use in many aspects16,17,18, involving classification models such as k-nearest neighbor and discriminant analysis19,20,21,22,23,24, support vector machines25,26,27,28,29,30,31,32,33,34,35, neural networks17,36,37,38,39,40,41,42,43,44,45,46,47,48 , and fuzzy clustering49,50,51,52,53,54. Regression models in terms of Raman signals have succeeded in predicting pH55. However, these machine learning methods focus on chemical information and rarely consider spatial information; that is, they treat Raman images as ensembles of Raman spectra, neglecting the spatial structure of how the spectra assemble to each other in space over the field of view (FOV)14,56,57. Chemical species in cells and tissues exhibit spatially-heterogeneous structures and patterns. Such spatial structure is considered versatile, especially in pathology, observing spatial contrasts in cytoplasm, nuclei of cells, and tissues in stained images58. Using the spatial information of Raman signals in addition to its chemical information can also accelerate measurements through selective illumination on the fly59,60,61,62. In addition, most pathological inspections of biological samples are performed in relation to morphological characteristics. For instance, nonalcoholic fatty liver disease (NAFLD), characterized by abnormal lipid accumulation in more than 5% of the liver without significant alcohol consumption, is classified pathologically into nonalcoholic fatty liver (NAFL) and nonalcoholic steatohepatitis (NASH) through its microscopic inspection of stained images by pathologists in terms of spatial structure such as hepatocellular inflammation, ballooning, and fibrosis63,64. Here, NAFL has generally a stable and non-progressive course, with the presence of lipid droplets in the hepatocytes. In contrast, NAFL followed by NASH is associated with a poor prognosis and the risk of advancing to liver cirrhosis or liver cancer65,66,67. Recent evidence suggests a more complex and diverse progression of NAFLD68,69. Research on rapid progression to advanced fibrosis and NASH with cirrhosis has also progressed68,70,71, and recently developed histological and staging criteria have improved the clinical diagnosis of NAFLD63,72,73,74. However, prediction for future fibrosis progression and rapid progressors among NAFL patients remains challenging68,70. While other hyperspectral imaging techniques, such as near-infrared spectroscopy48, provide powerful capabilities for biological imaging, spontaneous Raman scattering microscopy offers distinct advantages essential for this study: non-invasive spectral collection from water-rich tissues, broad spectral coverage to capture chemo-spatial diversity, and high-sensitivity detection of biomolecules including cytochromes and vitamin A, facilitated by resonance Raman effects. Furthermore, the non-necessity of sample staining is another advantage compared to fluorescence-based hyperspectral imaging75.

In this paper, we present a method based on information theory to classify Raman images of biological samples in terms of both chemical (Raman) information and spatial information of chemical signals by incorporating an auxiliary variable in the classification. The auxiliary variable quantifies the spatial heterogeneity of the spectra, that is, how much Raman spectra reflecting the microscopic chemical states at neighboring pixels are different from each other. We analyze liver tissues dissected from the NAFLD rat model annotated with pathological classifications such as normal, NAFL, and NASH, respectively. We address how the addition of spatial information on top of chemical information further provides a clue to differentiate tissue state, especially at NASH samples. In what follows, we refer to the spectral information combined with the spatial position information as chemo-spatial information.

Results

Spatial heterogeneity in a Raman image

Fig. 1
figure 1

Procedure for computing the descriptor of spatial heterogeneity \(p_{V|R}(v|r)\). The first step is to calculate the spectral distances \(v\) between a pixel r and its surrounding pixels \(s \in S_r\) located within the circle centered at r with a radius \(\lambda _0\). The second step is to compute the probability density function of \(v\), \(p_{V|R}(v|r)\), from the histogram of \(\{v(I_{r},I_{s}) \mid s\in S_{r}\}\). \(p_{V|R}(v|r)\) quantifies how the Raman spectra at pixel r differs from those at the surrounding pixels.

Full size image

Spatial heterogeneity is defined at each pixel based on the Raman spectral distances between the pixel and its surrounding pixels as follows (see Fig. 1): We define a set of the surrounding pixels around the pixel r as the pixels in the circle centered at r with a radius of \(\lambda _0\). We denote the set of the surrounding pixels of r by \(S_r\). In this paper, as a typical spatial scale, we employed \(\lambda _0=20\) µm, ca. 4 pixels, corresponding to being slightly larger than the typical size of hepatocytes so that our chemo-spatial information is relevant to the spatial diversity in Raman signals in the scale of hepatocyte.

In Fig. 1, for two given center pixels r and \(r'\), their Raman spectra at a given wavenumber \(w\) are denoted as \(I_r(w)\) and \(I_{r'}(w)\), respectively (blue box in Fig. 1). In this example, the set of Raman spectra near each of the pixel r and \(r'\) is collectively depicted in the orange box (i.e., the spectra at the surrounding pixels \(s \in S_r\) and \(s' \in S_{r'}\)). One can see that the diversity of Raman spectra is apparently different between the neighborhoods of the two center pixels r and \(r'\). It is clearly impossible to capture this only through the Raman signals of the central pixel. To quantify such diverse characteristics of the chemical composition in space, we first introduce a spectral distance \(v\) between pixel r and one of its surrounding pixels \(s \in S_r\), defined by

$$\begin{aligned} v\left( I_{r},I_{s}\right) = \sum _{w \in W}|I_{r}\left( w \right) - I_{s}\left( w\right) | \end{aligned}$$
(1)

where W denotes the set of all measured wavenumbers.

Not a single Raman spectrum at pixel s near the center pixel r, but the set of surrounding Raman spectra characterizes the spatial heterogeneity at the given center pixel through their spectral distances. That is, for any given (center) pixel r, a set of different \(v\) exists. This provides us with the conditional probability distribution of the spectral distance, \(p(V=v|R=r)\), as a descriptor to quantify the spatial heterogeneity for a given Raman image. Here, according to the standard notation in statistics, capital letters V and R denote stochastic variables, while the small characters \(v\) and r do some of their actual values. In this paper, we briefly represent \(p_{V|R}\left( v| r\right)\) and when we refer to the function with any value of \(v\) and r we use \(p_{V|R}\) (see details of how to compute the conditional probability distribution in Supplementary Information (SI)). The spatial heterogeneity dictated by \(p_{V|R}(v|r)\) reflects the spatial structure of the Raman image within the spatial scale defined by the radius \(\lambda _0\). Figure 1 exemplifies the probability distribution functions (PDFs), \(p_{V|R}(v|r)\) and \(p_{V|R}(v|r')\), respectively. The former has a larger variance with a longer tail, reflecting its heterogeneous chemo-spatial information near the pixel r, while the latter has a smaller variance with a sharper peak, reflecting its homogeneous local chemo-spatial information in the vicinity of the pixel \(r'\). Note that random shuffle of the spatial positions of Raman spectra results in completely different \(p_{V|R}\) profile from original one, while it does not affect any result of commonly used clustering analyses (e.g.,14,17,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54) because they do not take into account any spatial information in the images (see SI Fig. S1).

Using the information bottleneck method76, we performed clustering based on the probability distribution function \(p_{V|R}\) as a descriptor. In the clustering, the total number of clusters was set to 4 (We confirmed that the conclusions below do not change for the number of clusters 3–6). The details of the clustering method are provided in the Methods section below.

Chemo-spatial clustering of Raman images of liver tissues from the NAFLD model

We applied our chemo-spatial clustering method to Raman images of NAFLD-model liver tissues. As recently reported77, a total of 48 rats were housed with different diets: standard diet (SD, 16 rats), high-fat diet (HFD, 60% lard, 16 rats) and high-fat high-cholesterol diet (HFHC, 60% lard, 1.25% cholesterol, 0.5% cholic acid, 16 rats). Four rats from each group were sacrificed at 2, 4, 8, and 16 weeks (after starting designated feeding) for sample collection. Raman images of sliced liver tissues excised from each sacrificed rat were taken. Raman images of NAFLD liver tissues were preprocessed using the well-known protocol as before49. The detailed protocol is provided in the Methods section below.

Figure 2a shows the conditional probability distribution of the spectral distance, denoted by \(p_{V|C}(v|c)\), which represents the chemical heterogeneity of pixels belonging to the cluster c. (Here, C denotes the stochastic variable that represents the cluster index.) The index of clusters is taken from 1 to 4 in decreasing order of its chemical heterogeneity, so that the larger the index is, the lesser the Raman spectral difference near the pixels (i.e., more homogeneous).

Fig. 2
figure 2

Analysis of the clustering based on chemo-spatial information. a Representative spectral distance distribution of each cluster, \(p_{V|C}\left( v|c\right)\), characterizing local chemical heterogeneity in the surrounding region centered on each pixel. Local heterogeneity increases as green\(\rightarrow\)blue\(\rightarrow\)yellow\(\rightarrow\)orange cluster because their spectral distances increase on average. b Spatial distribution of clusters (cluster map) on each rat’s liver tissue image based on chemo-spatial information by using information bottleneck method. Colors are the same as (a). It should be noted that local heterogeneity tends to decrease in the set of tissues pathologically identified as NASH, and some of NAFLs.

Full size image

Then, how is each cluster distributed over space depending on the state of NAFLD? What can be the additional information we may acquire compared to a simple clustering based solely on Raman spectra? Fig. 2b shows the spatial distributions of the clusters, hereinafter referred to as ‘cluster maps’, for each of the liver tissue samples from rats, classified by the dietary condition SD, HFD, and HFHC. The cluster colors on the cluster maps are the same as those in Fig. 2a. In the figure, for example, the label HFD8w3 means a rat #3 out of the four individuals which were fed a high-fat diet (HFD) for 8 weeks. The pathologically identified tissues states, Normal, NAFL and NASH, are indicated by thin black, thick pink, and thick blue frames, respectively. As the overall trend of chemo-spatial analysis for NAFLD, the worst disease tissues NASH exhibit chemically spatial homogeneous patterns (see green contiguous regions, which are the most chemically homogeneous cluster). Interestingly, not all, but some NAFL tissues from the HFHC diet model also exhibit similar chemical homogeneity to NASH (e.g., HFHC4w2 and HFHC8w2). In terms of dietary model categories, Cluster 3 (blue) is distributed over the tissues of the SD and HFD images, in which more chemically spatial heterogeneous Clusters 2 (yellow) and 1 (orange) are scattered sparsely. In turn, HFHC images seem to be divided into two classes roughly in terms of the abundance of chemical spatial homogeneous Cluster 4, not necessarily pathologically identified NAFL and NASH classes. That is, HFHC4w1, HFHC4w2, HFHC8w2 and most HFHC16w images (except HFHC16w4) have highly populated Cluster 4 compared to the other HFHC images. It should be noted that such a distinct pattern of chemo-spatial heterogeneity across dietary and pathologically identified states cannot be captured by standard clustering based solely on Raman spectral features free from spatial information [see further detailed analysis in the Supplemental Materials (Fig. S1)].

Fig. 3
figure 3

Analysis of clustering based solely on chemical information using the fuzzy clustering method77. a The average spectra for the same clusters. The intensity of vitamin A (1585 cm\(^{-1}\)) decreases while that of lipids (2800-3000 cm\(^{-1}\)) increases with the cluster index increases from 1′ to 7′. b Spatial distribution of clusters (cluster map) on each rat’s liver tissue image.

Full size image

Then, what chemical species actually affect chemical heterogeneity in space? To address this question, let us compare the current results with the cluster analysis based solely on Raman spectral information existing in the same field of view (FOV) without referring to the spatial location14,77. Hereinafter, we call this ‘chemical clustering’ while the former clustering taking into account chemo-spatial heterogeneity in Raman signals is called ‘chemo-spatial clustering.’

Figure 3a presents the results of chemical clustering analysis in which seven different spectral clusters are identified14,77. Note that the resolution in cluster difference is higher in chemical clustering than in chemo-spatial clustering. Fat accumulation increases from Cluster 1′ to Cluster 7′ (see the high-wavenumber region of the lipid 2790-3045 cm\(^{-1}\)) while vitamin A decreases (1595 cm\(^{-1}\)). Figure 3b shows the spatial distribution of the clusters in each image. These indicate the localization of vitamin A in normal tissues, while some NAFLs with HFHC diets (termed NAFL-\(\beta\)77) and NASH seem to exhibit a wide range of fat accumulation. It should be noted that the size of cluster (i.e., population) in such standard chemical clustering analysis is exactly preserved, even when the spatial information of Raman signals is completely spoiled by random shuffling among pixels (see Fig. S1). There seem to be partially similar patterns as Fig. 2b, e.g., spatial patterns of Cluster 7′ (gray) (\(\sim\) lipid-rich and less vitamin A and cytochrome) in the chemical clustering and Cluster 4 (green) (\(\sim\) most homogeneous) in the chemo-spatial clustering (Fig. 2b) seem to be similar in the HFHD dietary models.

Fig. 4
figure 4

a Joint probability statistics taken over all liver tissue Raman images. The joint probability value at cth row and \(c'\)th column in the heatmaps means the ratio of the number of pixels simultaneously belonging to the cth cluster in the chemo-spatial clustering and the \(c'\)th cluster using the chemical clustering to the number of all pixels in all images applied in the clusterings. The cluster maps (first row: chemo-spatial clustering, second row: chemical clustering) of HFHC4w2 and HFHC8w2 illustrates the chosen two points classified as \((c,c')=(4,7')\); counting the total number of such pixels over all images and dividing it by the total number of pixels in all images results in \(p_{C,C'}(4,7')\) for all images. b Statistics taken over the liver tissue Raman images only for the same pathological states with the same diet. The marginal probability distributions \(P_C(c)(=\Sigma _{c'} p_{C,C'}(c,c'))\) and \(P_{C'}(c')(=\Sigma _{c} p_{C,C'}(c,c'))\) are attached to each figure vertically and horizontally in (a) and (b) .

Full size image

To quantify how chemo-spatial clustering and the chemical clustering results are related to each other, we examine the joint probability distribution \(p_{C,C'}(c,c')\) of the two types of clustering, shown in Figs. 2b and 3b, that is, the probability of a pixel belonging to the cth cluster in chemo-spatial clustering and the \(c'\)th cluster using chemical clustering simultaneously. As seen in the cluster maps (first row: chemo-spatial clustering, second row: chemical clustering) of HFHC4w2 (NAFL) and also HFHC8w2 (NAFL) in Fig. 4a, for example, the chosen two pixels both belong to Cluster 4 in chemo-spatial clustering (i.e., \(C=4\)) and Cluster \(7'\) in chemical clustering (i.e., \(C'=7'\)) simultaneously, and are counted in the total number of such pixels \(N(C=4,C'=7')\). For example, the joint probability distribution for a Raman image is computed by \(p_{C,C'}(c,c')=N(C=c,C'=c')/N_\text {pixel}\). The leftmost figure in Fig. 4a shows the heatmap of the joint probability distribution \(p_{C,C'}(c,c')\) whose statistics are taken on all Raman images of the liver tissue in order to glimpse the overall relationship between the two clustering schemes. The horizontal and vertical axes correspond to the cluster indices of chemical clustering free from spatial structure and those of chemo-spatial clustering, respectively. Cluster 3 has the highest probability in chemo-spatial clustering, and Cluster \(4'\) has the largest in chemical clustering (see the marginal probability distributions \(P_C(c)\) and \(P_{C'}(c')\) in Fig. 4a). Very roughly, relatively high joint probabilities appear to align diagonally on the heat map as the overall trend: the more the cluster index increases in chemo-spatial clustering [i.e., more spatially homogeneous (Fig. 2a)], the more the cluster index in chemical clustering increases [i.e., more lipids and less vitamin A and cytochrome contents (Fig. 3a)]. This might coincide with the fact that Figs. 2b and 3b look similar as overall trend.

Fig. 5
figure 5

Detailed comparison between the chemo-spatial and chemical clustering maps dependent on HFD and HFHC diet models with pathological classification. Each figure accompanies the information a–f: a the tissue sample label of diet, weeks and index of the individual rat such as HFHC16w4. b The heatmap of the joint probability distribution \(p_{C,C'}(c,c')\) whose statistics are taken over each liver tissue labeled by (a). White to red color grade corresponds to the joint probability value being 0 to 1. c The probability mass function \(p_{C'}(c')\) of chemical clustering whose total number of clusters is seven. The color indicates the index of cluster. Purple to gray bars indicate the densities of Cluster 1′ to 7′, having approximately low to high lipid content and high to low Vitamin A. d The probability mass function \(p_C(c)\) of chemo-spatial clustering whose total number of clusters is four. Orange to green bars indicate the densities of Cluster 1 to 4 in the descent order of chemical heterogeneity in space. e Cluster map of chemo-spatial clustering. f Cluster map of chemical clustering ignoring any spatial information in the samples. Note in pathologically-identified NASH liver tissues that \(p_{C,C'}(c,c')\) heat maps of HFHC16w1, HFHC16w2 and HFHC16w3 exhibit the similar distribution as well as the two corresponding cluster maps, whilst HFHC16w4 and HFHC4w4 differs from the typical HFHC16w model.

Full size image

However, when the joint probabilities are taken over the liver tissue ensemble only for the same pathological states with the same diet, the pattern of the joint probabilities becomes more diversified (Fig. 4b). For example, (diet, pathological state) = (HFD, NAFL) corresponds to the set of liver tissues of rats that have been fed a high-fat diet for a relatively long period (4, 8, 16 weeks) and are identified as pathologically NAFL. Likewise, (HFHC, NAFL) corresponds to the set of those fed a high-fat, high-cholesterol diet (e.g. 2, 4, 8 weeks) and are identified pathologically as NAFL. As an overall trend as (SD, Normal) \(\rightarrow\) (HFD, Normal) \(\rightarrow\) (HFD, NAFL) \(\rightarrow\) (HFHC, NAFL) \(\rightarrow\) (HFHC, NASH), the most populous cluster C shifts from Cluster 3 to 4 in chemo-spatial clustering and the most populous cluster \(C'\) shifts from Cluster \(4'\) to \(7'\) in chemical clustering (see \(P_{C'}(c')\) and \(P_C(c)\) attached to each figure vertically and horizontally in Fig. 4b). This indicates that as NAFLD progresses, lipids increase but vitamin A and cytochrome decrease in liver tissues, while their local chemical environment becomes more homogeneous within the typical spatial scale of a hepatocyte.

The first striking consequence in Fig. 4b is that in the same disease state pathologically identified as NAFL with different diet models (HFD and HFHC), the joint probability distribution \(p_{C,C'}(c,c')\) of (HFHC, NAFL) looks closer to that of (HFHC, NASH) than that of (HFD, NAFL), indicating that the liver tissues of (HFHC, NAFL) resemble those of (HFHC, NASH). This indicates that before the morphological signatures of NASH emerges, Raman imaging has some potential to infer the closeness to NASH (i.e., ill prognosis) among the NAFL not only chemically but also chemo-spatially. Likewise, when (SD, Normal) and (HFD, Normal) are compared to each other for the same pathological state ‘Normal,’ this joint probability distributions are not necessarily similar. For both, the most populous cluster in chemo-spatial clustering is Cluster 3, which may indicate that the heterogeneity of local chemical environments does not change largely but chemical ingredients change (the most populous cluster changes from Cluster \(4'\) to Cluster \(5'\)).

Fig. 6
figure 6

The Raman signature in a region in which the cluster map is uniform in chemical clustering but rather heterogeneous in chemo-spatial clustering. a Cluster map of HFHC16w4 generated by our chemo-spatial clustering. b Cluster map of HFHC16w4 by chemical clustering. c Raman Spectra measured within the region indicated by the white circle in (a) and (b). The thick black line represents the spectrum at the center pixel, and the other lines are the spectra in the surrounding pixels. The different colors used for depicting the Raman spectra are aimed at attempting to differentiate each of the spectra (not indicating different clusters).

Full size image

Let us look deeper into the image-wise analysis especially for the liver tissue Raman images pathologically identified as NAFL and NASH, because NAFL is a “driver” to transition to NASH. In Fig. 5, juxtaposing the two different types of cluster maps (chemo-spatial and chemical clusterings), the set of the heat map of the joint probability distribution \(p_{C,C'}(c,c')\) as well as its marginals \(P_C(c)\) and \(P_{C'}(c')\) attached vertically and horizontally for each of these images is presented. Comparison of the heat map in each image of rats fed the HFHC diet for 16 weeks (HFHC16w)—expected as the most severe condition with the state of NASH — with the overall (HFHC, NASH) heat map (Fig. 4b) tells that the images of HFHC16w1, HFHC16w2, and HFHC16w3 follow the general characteristic of the HFHC-NASH heatmap. In chemo-spatial clustering, the spectra of HFHC16w1, HFHC16w2, and HFHC16w3 belong predominantly to Cluster 4 and simultaneously to Cluster \(7'\), showing that chemical environment characterized with more lipids and less vitamin A and cytochromes spreads over a large spatial region homogeneously.

The second striking consequence of chemo-spatial clustering is as follows: Among the HFHC16w rat model, the fourth individual (represented as HFHC16w4) deviates from the most populous averaged characteristics of chemo-spatial clustering: Unlike the other HFHC16w models, many spectra belong to Clusters 1 through 3 rather than dominantly to Cluster 4. This is opposite to the situation of chemical clustering in that the Raman spectra of HFHC16w4 belong mainly to Cluster \(7'\) just as the other HFHC16w model. This coincides with our visual inspection of the two types of cluster maps of chemical and chemo-spatial clusterings. Even though Raman spectra indicate populous lipid with fewer vitamin A and cytochromes (from the dominance of Cluster \(7'\)), the pattern of how chemical ingredients are distributed in tissue differs from the averaged Raman signature. Note that the distinction between HFHC16w4 and the other HFHC16w model is not manifested using chemical clustering because it does not take into account the spatial structure of Raman signals, and visual inspection of the cluster map for chemical clustering appears not to show significant differences between the HFHC16w model accompanying the populous Clusters \(6'\) and \(7'\).

Finally, we address why spatial heterogeneity increases within regions where similar spectra are distributed in HFHC16w4. To consider the origin of the increased spatial heterogeneity, the cluster maps of HFHC16w4 in chemo-spatial (Fig. 6a) and chemical clusterings (Fig. 6b) are shown with the radius of 4 pixels indicated by the white circles (\(\sim\) the typical size of hepatocyte) composed of similar spectra (Cluster 7′) but relatively heterogeneous in space. The Raman spectra within the white circle zone, which were used to calculate the spatial heterogeneity, are shown in Fig. 6c (c.f., the corresponding Raman spectra in a region of HFHC16w1, shown in Fig. S2, where both clustering methods present a uniform distribution, mainly composed of a single cluster).

Fig. 7
figure 7

The relationship between chemo-spatial heterogeneity and its Raman intensities of hemoglobin-related, vitamin A and lipids across HFHC16w family. White circles indicate the same position in Fig. 6a and b in HFHC16w4, Fig. S2a and S2b in HFHC16w1 as references. Black circles and arrows in HFHC16w2 and HFHC16w3 points regions where heterogeneous ingredients expected to be present. a Heat maps of Raman shifts in HFHC16w images related to Cytochrome, Hemoglobin, Vitamin A and Lipids. Here, to compare the Raman intensities over the different HFHC16w images, focusing on the relatively constant value at 1658 cm\(^{-1}\), reflecting the stretching of C=O at Amide-I78, we standardized the spectra by calculating the ratio of the intensities at each Raman shift divided by the corresponding peaks at 1658 cm\(^{-1}\) as a reference. The use of Amide-I as a reference peak is considered appropriate, since the parent protein of Amide-I is considered to distribute throughout the biological sample. The range of color bar is specified in each Raman shift and by the values from 0.005 quantile to 0.995 quantile of the data set. b Cluster maps of the HFHC16w family generated by chemo-spatial clustering and chemical clustering.

Full size image

The comparison of Fig. 6c and Fig. S2c suggests that there may be three distinct features in the Raman spectra of the HFHC16w4 model. The first characteristic is the relatively high intensities with fluctuations in the blood component (hemoglobin). In Fig. 6c, several peaks are assigned to specific vibration modes of heme b and c, or both, such as 750, 1127, 1360, 1550, 1585, and 1635 cm\(^{-1}\)79,80,81,82. In these experimental conditions, heme b can be mainly attributed to hemoglobin in blood and heme c to cytochrome c. Among the peaks attributed to heme above, 1360 and 1550 cm\(^{-1}\) are generally more intense in heme b than in heme c, and 1635 cm\(^{-1}\) are assigned only to heme b (oxygenated). In Fig. 6c, the intensities of these three peaks exhibit a correlation with each other in the spectra. In this study, we were unable to clarify the origin of hemoglobin, but the most plausible scenario is the blood residues in the tissue from the step of its excision and slicing. This study used freshly resected liver samples: The influence of red blood cells that escape from the resected area cannot be completely ruled out. This interpretation is also supported by the result on the presence of a higher amount of oxygenated hemoglobin (1635 cm\(^{-1}\)) in HFHC16w4 compared to HFHC16w1.

The second is the presence of relatively abundant vitamin A. Spectra with peaks at 1015, 1160, 1200, 1275, and 1595 cm\(^{-1}\) attributed to vitamin A exist within the local region. As seen in the HFHC16w1 model, few spectra reflecting vitamin A were measured. Vitamin A is known to be localized in fat droplets under normal conditions and its decrease is known to occur in several types of disease, including fatty liver83. Thus, the presence of abundant vitamin A suggests that the pathology of HFHC16w4 is still in progress, compared to HFHC16w1. The third feature is the large fluctuation in the Raman intensity in the fat region. The regime indicated by the white circle in HFHC16w4 —where chemical clustering provides almost a single cluster, although chemo-spatial clustering indicates spatially heterogeneous— has higher fat fluctuations than that in HFHC16w1.

To confirm whether the above discussions hold for all pixels in the HFHC16w family, Fig. 7 presents the relationship between chemo-spatial heterogeneity and its Raman intensities of hemoglobin-related, vitamin A, and lipids across that family. Fig. 7a presents the heat maps of Raman intensities at 1130 cm\(^{-1}\), 1360 cm\(^{-1}\), 1550 cm\(^{-1}\), 1635 cm\(^{-1}\) (all hemoglobin-related), and vitamin A at 1595 cm\(^{-1}\), and lipids at 2856 cm\(^{-1}\), respectively, standardized by Amide-I Raman shift, with the cluster maps of the HFHC16w family generated by chemo-spatial clustering and chemical clustering in Fig. 7b. One can see that not only in the white circle region, but also in the other regions on the heat maps of HFHC16w4, hemoglobin-related ingredients are more observed than in the other HFHC16w family (that is, redish pixels are more abundant in across the hemoglobin-related HFHC16w4 image). For example, in HFHC16w2 and HFHC16w3, hemoglobin ingredients at 1550 cm\(^{-1}\) are sparsely distributed (indicated by black circles and black arrows in Fig. 7a). Chemical clustering (with total number of clusters seven) cannot discriminate such sparsely distributed minor components. Specifically, chemo-spatial clustering identifies the hemoglobin-related minor components within the right black circle in HFHC16w3 as belonging to Cluster 2 (yellow) (i.e., more heterogeneous than its neighbor’s Cluster 4 (green)) while chemical clustering does not. Likewise, chemo-spatial clustering differentiates the minor components within the black circles in HFHC16w2 as belonging to Cluster 2 (yellow) (i.e., more heterogeneous than its neighbor’s Clusters 3 (light blue) and 4 (green)) while chemical clustering does not differentiate them from the neighbors, and assigns them as belonging to Cluster \(6'\) in similar to the cytochrome rich (1130 cm\(^{-1}\)) regions.

In order to further quantify predictive performance of the two types of clustering schemes and address the limitation of chemo-spatial clustering, we synthesized a set of Raman images in which chemical contamination is artificially added to the Raman images of HFHC8w2 and SD16w3 as representatives of relatively homogeneous and heterogeneous chemical environments, respectively, and evaluated to what degree chemo-spatial clustering and chemical clustering can predict chemical contamination. Detailed results of the analysis are presented in Supporting Information, Section 4. In the analysis, cholesterol (resp. polystyrene) was randomly introduced into 10% of the pixels of the entire image with several mixing ratios, e.g., mixing ratio 20:80 means that the original Raman spectrum located in the chosen pixel is contaminated 20% by cholesterol (resp. polystyrene) (Fig. S3). The deviations between the clustering results of the two clustering schemes before and after chemical contamination are shown in Figures S4 through S7, as well as the corresponding analysis of the F-measures with respect to the prediction of which pixels are chemically contaminated using the change of the cluster indices of each pixel (Fig. S8). For example, we found that the F-measures increase as the mixing ratio of the contamination increases in both clusterings, and that chemo-spatial clustering yields higher F-measures in a more chemically homogeneous environment, although there is no such systematic trend and the F-measures are relatively low in chemical clustering. However, for chemo-spatial clustering, contamination is harder to locate in more heterogeneous environments. Overall, chemo-spatial clustering is complementary to conventional chemical clustering, which takes into account spatial information of Raman signals across a sample and is more sensitive to chemical contamination.

Conclusion and outlook

We have developed an information-theoretic clustering approach taking into account chemo-spatial information without compromising the structure of hyperspectral images. We introduced the spatial heterogeneity measure to quantify local chemical heterogeneity in space, to which the clustering is performed. In the NAFLD disease model of rats, we demonstrated that the increase and decrease in certain chemical species and the homogenization of the local chemical environment occurred simultaneously. Some tissues in NAFL also exhibited similar surrounding chemical environments as seen in NASH. This feature in NAFL was not apparent from the morphology-based observations. Furthermore, this chemo-spatial clustering revealed in some of the NASH tissues the existence of spatially heterogeneous blood deposits, which could not be easily characterized by chemical clustering.

In June 2023, NAFLD and NASH were renamed metabolic dysfunction-associated steatotic liver disease (MASLD) and metabolic dysfunction-associated steatohepatitis (MASH), respectively84. This change in nomenclature involves a review of diagnostic criteria with a focus on metabolic abnormalities and a comprehensive understanding of the disease, including those who drink small amounts of alcohol. However, MASLD/MASH and NAFLD/NASH are conceptually similar, so the findings of this study, derived from the NAFLD/NASH model, may apply to the updated concept.

Since histopathologies utilize tissue morphological information in stained images, it is natural to reference the spatial structure of Raman images in addition to chemical compositions. In this study, we set the spatial heterogeneity range at a radius \(\lambda _0\) of 20 µm corresponding to the typical size of hepatocytes. Increasing the radius reflects a larger chemical environment in the space, while the difference among each pixel in the spatial heterogeneity measure is expected to decrease because of larger overlap of regimes centered at the individual pixels. Conversely, decreasing the radius may induce a larger fluctuation irrelevant to the subcellular and cellular structure. One of the future subjects to be studied should be the dependency of classification performance on the size of local regions within which chemical heterogeneity in physical space is characterized.

In this paper, we used the information bottleneck method to group the spatial heterogeneity distribution \(p_{V|R}\) at each pixel. This method prevents the whole data set R from being divided into each group more than the experimental error. For example, the cluster C should capture less detailed information about the target data V than the pixel data R does, because the cluster representation C only retains a simplified or “coarse-grained version” of the details of the pixels in R. Due to inherent experimental variation, the degree of information captured by R about V can fluctuate between different experimental instances. Meanwhile, as the number of clusters increases, C tends to capture more information about V. Therefore, if the information retained by C in the current set-up exceeds that retained by R in some (hypothetical) experimental instance, then it would be advisable to reduce the number of clusters because current C captures too much detailed information about V (see also85).

The proposed method may be useful for detecting foreign bodies with specific geometries present in collected tissue samples, e.g. asbestos or microplastics. For example, it is known that various types of asbestos can be discriminated by Raman microscopy even in biological tissues. However, the detection of their Raman spectra requires a relatively higher excitation energy (such as a laser exposure time of several hundred seconds), thus the comprehensive detection of their distribution in tissues by spontaneous Raman imaging has been technically limited86. One of the advantages of this method is the capability to increase the detection sensitivity of specific targets by detecting differences in chemical compositions and spatial distribution. This may make it possible to detect such foreign bodies in biological tissue even under the condition that the Raman spectrum of the target is weak or its pattern is unknown.

Finally, we address the difference of our information-theoretic approach from the deep learning (DL) approach of Raman image analysis. When the Raman images are employed as the input data into the DL, in principle, any chemical and chemo-spatial information would be taken into account in much more abstract, black-box fashion if the spatial information is also versatile to differentiate the disease states of the sample using the Raman images. However, this is faced on the issue of information leakage16, in which DL can use irrelevant information to the sample, such as (unnoticeable) experimental conditions involving inhomogeneous irradiation profiles87. Despite this issue, when appropriate preprocessings of Raman images in question can be undergone, DL could be another means to address both chemical and chemo-spatial information in differentiating the sample states.

Methods

Information bottleneck method

We use the information bottleneck method76 to perform clustering according to spatial heterogeneity \(p_{V|R}\). The information bottleneck method computes the probability that a pixel r belongs to a class (or a cluster) c, denoted by \(p_{C|R}(c|r)\), by solving the following variational problem:

$$\begin{aligned} \operatorname*{minimize}_{p_{C|R}} \left( I(C;R)-\beta I(C;V) \right) . \end{aligned}$$
(2)

Here, the capitals C and R denote the stochastic variables that represent classes and pixels, respectively. I(C; R) is the mutual information between the classes and pixel data, measuring how much information from the original data remains in the clustering result. In turn, I(C; V) is the mutual information between the classes and the spatial heterogeneity data, indicating how much the spatial heterogeneity information is retained in the clustering result. Again, the capital V represents the stochastic variable of spatial heterogeneity. The objective functional of the above variational problem represents a trade-off between maximizing data compression (i.e., the smaller I(C; R) more compresses or unifies the pixel data) and maximizing information preservation of the spatial heterogeneity (i.e., the larger I(C; V) more preserves the details of spatial heterogeneity, requiring more number of clusters), where the parameter \(\beta ~(>0)\) controls the tradeoff between the two quantities.

A formal solution of Eq. (2) is given by76

$$\begin{aligned} p_{C|R}(c|r)=\frac{p_C(c)}{Z(r,\beta )} \textrm{e}^{-\beta D_{\textrm{KL}}\left[ p_{V|R}(v|r)\Vert p_{V|C}(v|c)\right] }, \end{aligned}$$
(3)

where \(Z(r,\beta )\) is the normalizing constant and the symbol \(D_{\textrm{KL}}\) denotes the Kullback–Leibler divergence measuring the “distance” between two probability distributions. In summary, as \(D_{\textrm{KL}}\) decreases — as the spatial heterogeneity profile at r (that is, \(p_{V|R}(v|r)\)) and the representative profile of class c (that is, \(p_{V|C}(v|c)\)) becomes closer, the probability \(p_{C|R}(c|r)\) increases more. In other words, the information bottleneck method classifies each pixel in the Raman image based on the similarity of the spatial heterogeneity characterized by \(p_{V|R}\).

In this study, the hyperparameter \(\beta\), which controls the trade-off between data compression and preservation, was set to 850 to ensure that the resulting clustering corresponds to a hard clustering—that is, the cluster attribution probability for each pixel, \(p_{C|R}(c|r)\), is close to one for a single cluster and close to zero for all others. This parameter choice was also adopted in Ref. 49. The total number of clusters was set at 4 (We confirmed that the conclusions of this paper do not change for the number of clusters 3–6).

Model animals, sample preparation, and measurement condition

This study was performed and described in accordance with the ARRIVE (Animal Research: Reporting of In Vivo Experiments) guidelines 2.0 (https://arriveguidelines.org/). Animal models were prepared according to the following established protocols with minor modifications10,88,89,90,91,92. Eight-week-old male Slc: Sprague Dawley rats (Shimizu Laboratory Supplies, Kyoto, Japan) were housed in 12-hour light/dark cycles with free access to food and water. As recently reported77, the 48 rats were divided into three groups receiving different diets: standard diet (SD), high-fat diet (HFD, 60% lard), and high-fat high-cholesterol diet (HFHC, 60% lard, 1.25% cholesterol, 0.5% cholic acid). All diets were purchased from Oriental Yeast Co., Ltd. (Tokyo, Japan). Four rats from each group were sacrificed at 2, 4, 8, and 16 weeks (after starting designated feeding) for sample collection. The weight (mean ± standard deviation) of rats at the sample collection were, in the unit of gram (g), 409.0±21.9 (SD), 430.2±27.7 (HFD), 401±27.9 (HFHC) for 4 weeks feeding, 472.5±35.5 (SD), 510.3±42.1 (HFD), 529.5±49.0 (HFHC) for 8 weeks feeding, and 563.3±33.3 (SD), 681.3±47.5 (HFD), 633.8±43.1 (HFHC) for 16 weeks feeding, respectively. In the sample collection process, rats were placed under deep general anesthesia by intraperitoneal injection of a combination of three anesthetic agents (0.1 mg/kg medetomidine, 3.0 mg/kg midazolam, and 5.0 mg/kg butorphanol), then their liver was surgically excised. Euthanasia of the rats was conducted by the gentle exsanguination after the liver extraction. All procedures were approved and are in accordance with the guidelines of the Committee for Animal Research of the Kyoto Prefectural University of Medicine (M2022-229/238). Liver tissues excised from each sacrificed rats were sliced to approximately 1 mm thickness, and immediately transferred to 4\(^{\circ }\)C Krebs-Henseleit buffer (KHB), and used for Raman measurements within 2 hours. In the measurements, the tissue slices were placed in glass bottom plates with fully immersion of KHB. The remainder of each liver organ was fixed with 10% formalin buffer solution (Wako Pure Chemical Industires, Ltd., Osaka, Japan), then paraffinized and cut to a thickness of 4 mm with standard protocols for each step. The slices were treated for staining with hematoxylin and eosin (HE) after deparaffinization, then classified into normal liver tissue, NAFL, NASH based on NAFLD Activity Score (NAS)63 and Brunt’s criteria73, focusing on the comprehensive cellular and tissue architecture including the extent of steatosis, lobular inflammation, hepatocellular ballooning, fibrosis staging. The classification process was performed by a pathologist and a hepatologist familiar with liver histopathology. The results of liver histopathology revealed the 26 normal livers, 15 NAFL, 6 NASH (and 1 fibrosis as an unexpected result) as the prepared liver conditions from the 48 rats. The liver histology was mostly the same in the same dietary condition but slightly different between the conditions. Especially, the diet of HFD and HFHC for 16 weeks provided NAFL and NASH, respectively, and the time points 2, 4 and 8 weeks were necessary to obtain the liver condition in progress to NAFL/NASH. Thus, the liver model reflecting the stepwise progress from normal liver to NAFL/NASH was prepared. Raman measurements of thick slices of the liver were performed with a commercial confocal Raman microscope (Raman-11, Nanophoton, Osaka, Japan) with a 532-nm excitation light source and a 20\(\times\)/0.75 dry-type objective lens (Olympus / Evacuation, Tokyo, Japan). Raman scattering of liver tissues was excited by point illumination of excitation light (68 mW · µm\(^{-2}\)), and detected by a CCD camera (Pixis 400BR, Teledyne Princeton Instrument, NJ, USA) at -70 \(^{\circ }\) C through a 50 µm confocal pinhole. Raman images with an area size of 95 µm \(\times\) 345 µm (20 \(\times\) 70 pixels) were recorded with a scan step of 5 µm and an exposure time of 1 s. Raman images of NAFLD liver tissues were preprocessed using the well-known protocol as before49.