Haralick texture feature analysis for quantifying radiation response heterogeneity in murine models observed using Raman spectroscopic mapping

Tumour heterogeneity plays a large role in the response of tumour tissues to radiation therapy. Inherent biological, physical, and even dose deposition heterogeneity all play a role in the resultant observed response. We here implement the use of Haralick textural analysis to quantify the observed glycogen production response, as observed via Raman spectroscopic mapping, of tumours irradiated within a murine model. While an array of over 20 Haralick features have been proposed, we here concentrate on five of the most prominent features: homogeneity, local homogeneity, contrast, entropy, and correlation. We show that these Haralick features can be used to quantify the inherent heterogeneity of the Raman spectroscopic maps of tumour response to radiation. Furthermore, our results indicate that Haralick-calculated textural features show a statistically significant dose dependent variation in response heterogeneity, specifically, in glycogen production in tumours irradiated with clinically relevant doses of ionizing radiation. These results indicate that Haralick textural analysis provides a quantitative methodology for understanding the response of murine tumours to radiation therapy. Future work in this area can, for example, utilize the Haralick textural features for understanding the heterogeneity of radiation response as measured by biopsied patient tumour samples, which remains the standard of patient tumour investigation.


Introduction
Radiation therapy is a standard treatment for approximately 50% of all cancer patients [1]. While significant improvements in the technological development of radiation therapy have occurred in the past several decades, a number of challenges in treatment efficacy remain unmet. Among these challenges, optimizing, or personalizing, radiation therapy remains PLOS  difficult due to the considerable inter-and intra-patient heterogeneity of response to radiation [2]. Indeed, heterogeneity of radiation response can exist within individual tumours, and can lead to differential patient response [3][4][5].
However, the precise mechanisms in which tumours establish radioresistance depend on numerous factors. For example, in radiation biology, it is well established that oxygen can provide the cell with a source of reactive species to generate DNA-damaging radicals. Moreover, oxygen may also contribute to the fixation of DNA damage once the initial insult has been established. It is fair to say that the complete mechanisms of radiation resistance and response in tumours is a complex combination of factors, and although differential responses to radiation therapy have been observed in the clinic for decades, the molecular basis of such responses remains an enigma. For a variety of cancers, recent studies have unequivocally highlighted the significant molecular heterogeneity that exists in patients' tumours and in tumour radiation response [6].
Tumour heterogeneity remains a challenge to measure in any scenario. Although a number of assays have been proposed in the literature, no one technique has proven to provide a comprehensive and clinically viable assessment strategy [7]. In previous investigations, we have demonstrated that Raman spectroscopy can offer multiplexed, biologically significant molecular-level information on cellular and tumour radiation response [8]. We have demonstrated that Raman mapping can be used to, for example, provide spatially resolved information on glycogen production in murine models of H460 lung tumours post irradiation [9]. However, the quantification of tumour response heterogeneity is challenging owing to architectural complexity, temporal changes, spatial variation, inherent subpopulations within host that are part of the tumour environment, and potential inaccuracies in data collection, just to name a few.
To overcome the issue of heterogeneity, it has been suggested that the average spectra be used as a representative of the target population, [1,2]. However, such strategies by their very nature lose information on the spatial origin of given biomolecular components, and is true for genomic studies in cancer. Textural analysis attempts to quantitatively describe characteristics of images based on the spatial arrangement of intensity values. While it has been established in pattern recognition [10], and image processing [11], it more recently has been finding application in the biomedical field [12][13][14][15][16]. For instance, textural features of PET scans extracted preand post-treatment from patients with esophageal cancer have been used to differentiate between nonresponders, partial responders, and complete responders [15]. Moreover, the use of PET imaging relies on tumor uptake of the radiotracing compound which could be impacted by the profusion of the microenvironment. In other work, textural measurements (such as heterogeneity, contrast, and energy) were observed to correlated with the fracture toughness of bone tissue [16]. While there is great potential for image analysis to improve our understanding of complex systems like tumours [4,5], this remain an active area of research.
Using the techniques described by Haralick [17], this paper calculates several textural features on Raman spectral maps. Besides using them to gauge tumour heterogeneity, these quantities are employed to assess dose response in murine models of lung tumours irradiated with clinically relevant doses of ionizing radiation. The next section outlines the Raman spectroscopy data, as well as the details on how the spectral maps were generated. Subsequent sections provide the specific definitions for select Haralick features and explore whether certain characteristics of these tumours change in response to radiation dose and time post irradiation.

Materials and methods
All mouse protocols, tumour growth parameters, irradiations, tumour sectioning, and Raman analysis were performed as described previously [9]. Briefly, H460 tumours were implanted and grown in murine models. Four mice were irradiated to 5 Gy, four mice were irradiated to 15 Gy, and a sham group (four mice) were left as controls, for a total of 12 mice. Tumours were harvested 3 days post irradiation and sectioned into slices using established protocols [9]. In terms of Raman analysis of resulting tumour sections, each Raman spectrum was collected utilizing a Renishaw Raman microscope (Renishaw Inc, Chicago, USA) operating at 785 nm laser excitation, a 100x dry objective, and a 10s acquisition time. Raman maps were collected on up to 4 separate tumour slices for each of the 12 mice in the study and were acquired on a 15 × 15 μm grid. The slices were selected randomly, as described in our previous work [9] and shown by way of illustrative example in Fig 1. In all cases 100 × 100 or 200 × 200 micron regions were studied.
Post-Raman acquisition spectral processing included cosmic ray removal, spectral smoothing, baseline subtraction, and volume normalization as described previously [8]. The only minor deviation here is within the baseline correction algorithm where we here use the baseline package for spectra data available in R [18] using the second derivative constrained weighted regression [19] (method = "als") with the second derivative constraint, lambda, equal to 4. Our complete data set, X n×p , contain n = 6581 spectra comprised of 74 Raman maps with a grid size of either 8 × 8 (15 μm × 15 μm per pixel) or 14 × 14 pixels (20 μm × 20 μm per pixel). That is, each map is stored in either 64 or 196 rows in X n×p , where Principal Component Analysis (PCA) [20][21][22] is a commonly used technique for data reduction that allows the user to project a high-dimensional data set into a new set of variables in a lower dimensional space. PCA is an orthogonal transformation which ensures that each variable, henceforth referred to as principal component, is uncorrelated and ordered such that the first principal component explains as much of the variability in the data as possible, the second component has the second largest possible variance, and so on. Applying PCA to our data results in a score matrix Y n×q with q � p where each row represents the corresponding Raman spectrum transformed into the new PC coordinate system. Herein, we concentrate on the first principle component which explains approximately 55 percent of the variability in our data. In previous work we have established that PC component 1 in the irradiated H460 tumours used in the present study corresponds to glycogen production post irradiation [8] (see S1 To generate our 2D Raman images (see Fig 2 for example and S2-S4 Figs in the supplementary data for the complete collection), we populate the grid matrices described above with PC component 1 scores (i.e. glycogen production scores) scaled between 0 and 1. In essence, the gray scale image provides a visual representation of tumour glycogen levels. As an example, Fig 3 represents how the pixel values for Map 27 relates to the Raman spectra. Spectra having high and low glycogen production scores are plotting in red and blue, respectively. We remark that the spectral interpretation in certain areas having a large discrepancy between red and blue spectra (eg. 490, 850, 1050 cm −1 ) correspond to spectral peaks characterized by glycogen [23].
The original paper by Haralick [17] introduced a total of 14 measurements; however, many of these values are highly correlated rendering all but five useful in practice [24,25]. These features rely on a gray level cooccurrence matrix (GLCM), denoted here by M, which contains the textural information of an image. Assuming we have a gray-scale image composed of N g gray levels, the ijth element of GLCM, denoted M(i, j), counts the number of times a pixel with the value of i is neighbouring a pixel with the value of j. Neighbours are defined by a user-specified distance d and angle θ. Finally, the GLCM is normalized so that the total sum of the elements sum to 1. We shall denote the elements of this N g × N g normalized GLCM by m d,θ (i, j), where we can write m d,θ (i, j) as (for the example case of θ = 0˚) [17]: Following the recommendations in [25], we focus on homogeneity (H), contrast (Con), correlation (Cor), entropy (E) and local homogeneity (LH). These texture features are calculated as follows: Con ¼ Cor ¼  [27] provide some helpful comments on how interpret the Haralick texture features. Briefly, homogeneity will be close to 1 when only a few dominant gray-tones are present. Contrast, will be equal to 0 for a constant image and become larger as the pixel intensities in a local neighbourhood become more disparate. Correlation ranges from -1 to 1 and reflects how correlated a pixel is to its neighbours. Entropy reflects the amount of randomness in intensity of an image and will increase with the images local complexity. Local homogeneity measures the similarity of pixels and is larger when there are minimal local changes.

Results
For each of the 2D gray-scale Raman images, the Haralick textural features given in Eqs (2)-(6) were calculated. Since the correlation between pixels is likely to decrease with distance, we use q = 1 in accordance with [28]. To avoid any directional dependency, we averaged the values across four angles θ = 0, 45, 90, and 135˚ [29].
For illustration, we have provided the Haralick features corresponding to the maps displayed in Fig 2. Images that appear more monochromatic (eg. Map 10) have higher values of H and LH than images having more pixel variation (eg. Map 34). Although Map 10 and 39 are similar in terms of homogeneity, the neighbouring pixels in image 10 are closer in gray-tone and consequently yield a lower value for contrast. Map 39 and 67 have a higher measure of gray-tone linear dependencies and consequently produce larger correlation values. On the contrary, the lack of any discernible linear pattern in Map 28, for example, generates a correlation closer to zero. Finally, images with larger entropy values exhibit more randomness in pixel intensity. Representative results are summarized in Table 1.
The Haralick features for the complete set of 74 gray-scale images are summarized in the box-and-whisker plots in Fig 4. Each of the twelve mice have 5-8 associated Raman maps; the legend provides a key to the mouse identification number. For each pair of box plots within each plot, a comparison for the difference of means is made using a Kruskal-Wallis test with a Among the 74 images considered, a significant difference between non-irradiated and irradiated groups was uncovered. For instance, the group of non-radiated mice has a significantly different mean contrast value than mice treated with 5 Gy and 15 Gy; p con 0 5 ¼ 0:016, p con 0 15 ¼ 0:021. Furthermore, there was no statistically significant difference between mice irritated with 5 Gy and 15 Gy; eg p con 5 15 ¼ 0:869. This trend, which is prevalent among all Haralick features considered, indicates that the textural complexity of the images extracted from irradiated mice tends to be greater than those extracted from non-radiated mice. As these images were generated from the score values linked to glycogen, this work suggests that the affect of radiation on glycogen production is not uniform across the tumour.

Discussion
The Haralick feature analysis provides a new methodology for shedding light onto the evolution of textural features within Raman maps of glycogen production within murine-irradiated tumours. For example, in this study we quantitiavely show a statistically significant decrease in the homogeneity of glycogen production (Fig 4a) as a function of irradiation dose. Put another way, glycogen production is heterogeneously distributed throughout the tumour and, furthermore, the extent of heterogeneity exhibits a radiation dose dependence. In a separate study we show that glycogen production is negatively correlated with tumour regression post radiation [30]. Furthermore, the current and past hypoxic state of local tumour morphology is correlated with glycogen production, thus affecting the spatial extent of tumour regression and glycogen production [30]. While we tackle the radiobiological implications of glycogen production in a separate work [30], it is clear that the Haralick features calculated here are (i) able to quantitate the extent of textural variation as a function of radiation dose, and (ii) correlate well with expected radiobiological trends in our murine models. Furthermore, it is apparent that glycogen heterogeneity exists not only intra tumour, but also between murine tumours. This variation is inherently interesting from the point of view of future assays based on tailoring treatment based on individual sample response and points towards future experiments dedicated to understanding this inter-murine variability.
Within our current approach, only a single principal component (PC) is considered. Although this component describes a large percentage (54.92%) of the total variability within the Raman dataset, the inclusion of additional PCs may be advantageous in terms of understanding variability of Haralick indices for PCs that may harbour additional biochemical information relating to radiation response. In order to analyze multiple PCs simultaneously, one could investigate an extension of Haralick textural features that uses the colours discriminators to describe colour images. In the context of this paper, rather than using a single PC score to create a gray-toned image, a pseudo-coloured Raman map could be generated using the first three PC scores to represent the red, green, and blue components in an RGB colour model; an example is provided in Fig 5. Owing to the fact that these principal components explain 79.34% of the total variability in our data, these images present a richer representation of the underlying spectroscopic information and could provide useful insight on the radiationinduced differences in succeeding PCs that may be linked to other biological interpretations. In future work, we will explore if and how the colour-indexed Haralick features are affected by this extension. Haralick features provide possibilities for future implementation in the long term strategy of personalized radiotherapy. For example, it is possible that Haralick analysis on biopsy-based Raman maps can provide information on response heterogeneity that, along with orthogonal biological information, can be used in the assessment of overall response. Other in vivo image modalities provide information regarding tumor architecture as well as spatial context relative to normal adjacent or host tissue within an organ system. In some cases such as PET, indirect assessments of phenotype (e.g. glucose uptake) of a single biomarker can help infer some biological phenotype (e.g. metabolically active). While improvements may allow for potentially dual or even triple biomarker uptake analysis, these techniques are different than Raman spectroscopy. Raman spectroscopy provides a global picture of biochemical signatures across multiple biomolecules in a single sample and does not rely on uptake of a tracer. As such, the information collected from imaging and Raman spectroscopy are complementary in this regard rather than overlapping. Future Raman spectroscopic point-of-care probes under research development may in the future alleviate the need for biopsy samples. Further research is required to continue to develop this long term goal.
On a final note, we recognize that an array of texture analysis methods exist in the literature. We here do not presume that Haralick features will necessarily outperform other methods. Rather, we here demonstrate that textural features provide an excellent platform to tackle the issue of response heterogeneity. While this well-known model serves as a reasonable first choice in our exploratory analysis, alternative techniques such as the ones described in [31]  could also prove useful. Additional work is required to establish the feasibility of other textural feature analysis techniques.

Conclusion
The Haralick texture features display biologically relevant trends for dose response in tumours extracted from our murine model. We here have demonstrated the ability of Haralick textural features to quantitatively characterize the evolution of textural components (e.g. contrast, homogeneity, entropy) within Raman maps of radiation-induced glycogen production in murine tumours. Our results indicate that Haralick feature calculations provide a new, quantitative assessment within a radiobiological context and may, in future studies, help quantify the extent of radiation response. In turn, such quantification may be valuable in the long term goal of personalized radiation therapy response monitoring.