Raman Spectroscopy and Machine Learning Reveals Early Tumor Microenvironmental Changes Induced by Immunotherapy

The emergence of immune checkpoint inhibitor (ICI)–based immunotherapy has revolutionized cancer treatment (1). The antitumoral effect of ICIs is orchestrated by successfully blocking the receptors on cancer cells (e.g., PD-L1) and tumor-reactive T cells (e.g., CTLA4 and PD-1) to limit the interactions that inhibit T-cell activation and unleash an immune response against cancer cells (2, 3). Although immunotherapy has shown incredible responses in certain cancer types such as melanoma (4, 5), renal cell carcinoma (6, 7), and colorectal cancer (8, 9), only a small subset of patients (sometimes as low as 25%) derive clinical benefit from the therapy (3, 10). Therefore, there is an urgent need for predictive biomarkers of response that can identify patients who would benefit from immunotherapy (11, 12). Early determination of nonresponders can aid in exploring alternative treatment strategies to alleviate the toxic immune-related adverse effects of ineffective therapy.

The current clinical metrics for prediction and evaluation of response to immunotherapy are not very effective (13, 14). For example, the use of IHC for PD-L1 expression as a predictive biomarker is highly debated due to the diversity of assays, lack of universal positivity thresholds, and observations of clinical benefit in patients with low levels of expression (14, 15). Increase in tumor-infiltrating lymphocytes has been proposed as an alternate marker, as it was associated with a better response of anti-PD1 therapy in a melanoma study involving a small number of patients (16). A liquid biopsy strategy combining the blood counts parameters, clinical characteristics, and serum lactate dehydrogenase predicted the response of patients without the metastatic disease to anti–PD-1 therapy with about 60% accuracy (17). Several imaging biomarkers are also currently at various stages of development (18, 19). FDG-PET scans of patients treated for melanoma and other cancers provided the evidence of an association between metabolic changes and response to immunotherapy (19). Studies have also leveraged PD-1/PD-L1 and CTLA4–targeting antibodies radiolabeled with ⁸⁹Z for evaluating the tumor uptake of therapeutics using PET imaging; however, such measurements are associated with challenges similar to PD-L1 IHC (19).

Tumor microenvironment (TME) composition and metabolism have been shown to play important roles in the regulation of immune response and disease progression (20, 21). Although some early studies provided evidence of a correlation between collagen density and immunosuppression, the compositional changes of TME in response to immunotherapy have remained largely unexplored (22, 23). A similar gap in knowledge exists due to the lack of investigations into the overall therapy–induced changes in the lipid composition of the tumors, despite emerging evidence of changes in lipid metabolism in individual immune cells due to immune checkpoint inhibition (24). Therefore, development of tools that can provide a holistic picture of the evolution of the TME in response to immunotherapy is highly desirable. Such methods are critical to discovering novel biomarkers of response that not only capture the differences in the composition of individual cellular and extracellular matrix constituents but also leverage their complex interactions to achieve higher sensitivity and specificity.

In this study, we propose Raman spectroscopy, a label-free technique based on the inelastic scattering of photons, to provide an alternate route for evaluating compositional changes in the tumor and its stroma in response to immunotherapy (25). In recent years, Raman spectroscopy has emerged as a valuable optical tool for non-perturbative study and molecular characterization of the biological specimen based on transitions in the molecular vibrational modes. The exquisite chemical specificity of this technique has been leveraged by us and others to study a variety of biological applications across length scales, including cell differentiation, drug uptake in cancer cells, and cancer onset and metastatic progression (26–29). However, there are no reports of Raman spectroscopic characterization of tumor response to immunotherapy.

Using a CT26 murine model of colorectal cancer treated with anti–CTLA4 or anti–PD-L1 antibodies, we sought to determine whether Raman spectroscopy is sensitive to the changes associated with ICI therapy. We used multivariate curve resolution–alternating least squares (MCR-ALS) decomposition of the spectral dataset to uncover pure component spectra that resemble those of the underlying biomarkers. The scores of these key component spectra that resemble lipid, nucleic acids, and collagen as well as their spatial heterogeneity change significantly with the administration of each ICI. Furthermore, leave-one-mouse-out implementation of support vector machines (SVM)–derived classifier displays high accuracy, thereby establishing the feasibility for automated recognition of prospective tumor samples according to their response to immunotherapy. Finally, the feature-ranking capability of bagged random forests was used to identify the spectral markers of response to each ICI therapy, which was corroborated through independent proteomics analysis. To the best of our knowledge, the current study is the first Raman spectroscopic investigation of changes in the tumor composition induced by immunotherapy. The results of our pilot study using a portable clinical Raman system show the substantial promise of optical spectroscopy in guiding personalized treatment decisions for patients.

The CT26 colon carcinoma cells were purchased from the ATCC (Product: ATCC CRL-2638, Lot Number: 58494154) and passaged according to established protocols. Cells were Mycoplasma negative when shipped from the third and fourth passages were used for tumor inoculation. Cells were injected into the hind flanks of 6–8-weeks-old Balb/c mice (n = 25) purchased from The Jackson Laboratory. The mice were maintained at the University of Arkansas Central Laboratory Animal Facility with 12 hours light/dark cycles and provided ad libitum access to food and water. All experiments were approved by the Institutional Animal Care and Use Committee at the University of Arkansas (IACUC protocol 18066, 20025). Mice were randomly assigned to one of three treatment groups—mouse IgG isotype (control, n = 10), anti-mouse CTLA4 (αCTLA4, n = 8), and anti-mouse PD-L1 (αPDL1, n = 7).

ICIs used in this study were purchased from BioXcell (Control: MPC-11, αCTLA4: UC10–4F10–11, αPDL1: 10F.9G2). Mice were treated with 3 doses of 100 μg anti–CTLA4 or 200 μg anti–PD-L1, or 100 μg of IgG control dissolved in 100 μL of sterile saline via intraperitoneal injection as done previously (30), when the average tumor volume was approx. 80–120 mm³. Mice were treated on days 1, 4, and 7 with day 1 indicating the first day of treatment (31). Tumors were excised and snap-frozen for Raman spectroscopy measurements and histopathology on day 10.

We performed Raman spectroscopic measurements on the tumors according to previously established protocols. The snap-frozen tumor samples were thawed in PBS and fixed in 10% neutral-buffered formalin. The fixed tumors were washed thoroughly in PBS, flattened, and sandwiched between a quartz coverslip (1 × 1 inch × 0.2 ± 0.05 mm) and an aluminum block with the regular addition of PBS to prevent dehydration during Raman measurements. The thin quartz coverslip is chosen due to the lack of significant fluorescence and Raman signature in the wavenumber region used in this study. The tumors were mapped using a fiber-optic probe mounted on a motorized translational stage (T-LS13M, Zaber Technologies Inc., travel range: 13 mm; ref. 32). The lensed fiber-optic probe (Emvision LLC, probe diameter: 2 mm) comprised of separate optical fibers for laser delivery and collection of backscattered Raman photons to provide an estimated tissue sampling volume of 1 mm³. The probe was interfaced with a custom-built portable clinical Raman system that consists of an 830-nm diode laser (Process Instruments, maximum power: 500 mW), a spectrograph (Holospec f/1.8i, Kaiser Optical Systems), and a thermoelectrically (TE)-cooled CCD camera (PIXIS 400BR, Princeton Instruments; ref. 32). Using the quartz coverslip to flatten the tumor and maintaining a uniform gap between the probe and the coverslip, the laser power was maintained at approximately 20 mW. The spectra were acquired from spatially distinct locations on both sides of the flattened tumor by raster scanning along a 2D rectangular grid of points spaced approximately 1 mm apart using a LabVIEW interface. We acquired a total of approximately 7500 spectra from the 25 tumors in treatment and control groups using 5 seconds (10 accumulations of 0.5 seconds each to prevent CCD saturation) acquisition time for each spectrum.

The analysis of Raman data was performed using in-house code developed in MATLAB 2017b (Mathworks) software. We calibrated the wavenumber axis using acetaminophen and restricted our analysis to the fingerprint region between 600 cm⁻¹ and 1,800 cm⁻¹ in this study. The acquired Raman spectra were pre-processed using fifth-order polynomial-based background subtraction, median filtering, and vector normalization (Euclidean norm of each spectrum is set to unity) to alleviate the contributions of fluorescence, cosmic rays, and laser power variations, respectively (33). The pre-processing steps preserved the size of the spectral dataset as the employed transformations did not involve any spatial averaging or dimensionality reduction steps.

Next, we used MCR-ALS analysis to decompose the spectral matrix into a small set of tissue constituent-like pure component spectra (loadings) and their contributions (scores) to each spectrum in the dataset under a non-negativity constraint on both the resultant matrices (34). The imposition of spectral equal length constraint (the Euclidean norm of each loading is set to unity) allowed us to use the MCR score for each loading as a measure of the abundance of the pure component it resembles and compare the scores across the treatment groups. We plotted the distribution of scores of key MCR loadings that resembled biological constituent spectra using violin plots with embedded box and whisker plots. The few outliers (defined as data points that are 1.5 times the interquartile range away from the bottom or top of the box) beyond the whiskers were not shown for clarity and optimal plotting. The statistical significance of differences in the medians of MCR score distributions across the classes was assessed using a two-sided Wilcoxon rank-sum test statistic with a conventional threshold (P < 0.05 to consider medians different). We also quantified the effect size of the differences across classes using the Wendt formula for rank biserial correlation (35). The MCR score maps were used to determine spatial heterogeneity using the method described previously by Brooks and colleagues (36) for medical images. The proposed index of heterogeneity quantifies spatial heterogeneity based on the smoothness of intensity transitions in the images (Supplementary Fig. S1A–S1D, Supplementary Information). For each 4 × 4-pixel sub-image within an MCR score map, the deviation of intensity transitions from the smoothest possible transition is calculated for all possible pixel pairs and integrated over normalized length to obtain heterogeneity index. The heterogeneity indices obtained for all sub-images are plotted using box and whisker plots to visualize differences in spatial heterogeneity induced by immune checkpoint inhibition therapy.

We used SVM–based classification routine to detect response to immunotherapy and to distinguish between the response of different immunotherapeutic agents. SVM uses constrained quadratic optimization to derive class boundaries in higher dimensional spaces for nonlinear classification problems. We used the LIBSVM library to perform C-SVM classification using a radial basis kernel for nonlinear mapping between the spectra and their class labels (37). The optimal C-SVM parameters (cost and kernel parameter gamma) were found by a grid search over a wide range. To avoid the spectral representation of tumors from the test dataset in the training dataset, we adopted a leave-one-mouse-out approach for our SVM analysis. This involves iteratively excluding all the spectra belonging to each mouse tumor from the training dataset and subjecting them as an independent test dataset to the SVM models trained on the remaining spectral dataset. The leave-one-mouse-out analysis allows objective assessment of the classifiers' ability to probe ICI-induced biomolecular changes beyond the expected biological variability. To avoid skewing the SVM predictions due to varying class sizes across the treatment groups, we performed 20 iterations of SVM training within each instance of the leave-one-mouse-out analysis by choosing a random subset of the training data that included equal representation of both the classes. For each spectrum in the test dataset, the median of predicted class labels obtained by testing it against all the SVM models derived from the 20 iterations was considered. The median predicted class labels thus obtained by testing all the spatially distinct spectra from every mouse allow us to assign an overall class prediction for each animal by thresholding on prediction frequencies. For the binary classification in the study, we assigned an animal to the majority class only if its membership was at least 30% higher than the random chance, that is, over 65%. By extension, the mice were misclassified as the incorrect class only if fewer than 35% spectra are correctly classified. All other mice, where the rate of correct classification lies between 35% and 65% were considered unclassified.

We trained bagged random forest classifiers using the TreeBagger class in MATLAB to quantify the out-of-bag error of the spectral prediction between the control and treatment classes. We trained 100 decision trees in the forest to track the evolution of the out-of-bag error with the increase in the number of trees and obtain measures of permutation predictor importance for all the Raman spectral features (38). We next ranked the features in the decreasing order of their importance to prediction and identified the top 5 non-neighbor features for each comparison to understand the similarities and differences between them.

The snap-frozen tissue was sectioned, and the slides were covered with a solution of 0.5% Oil Red O in propylene glycol and incubated for 45 minutes for lipid staining. Slides were rinsed, and counterstained in hematoxylin, dipped in differentiation solution, and blue buffer. The finished slides were mounted with Fluoromount-G and imaged with a Nikon wild-field microscope at a ×20 magnification. The formalin-fixed tumors were embedded in paraffin and sectioned on to glass slides for histopathology after the acquisition of Raman spectra. The sections were used to perform hematoxylin and eosin staining to observe tumor morphology and Masson's trichrome staining for collagen using standard protocols by the Oncology Histology Core at Johns Hopkins University. The slides were imaged using a Nikon fluorescence microscope using a ×10 objective.

FFPE tissue scrolls were deparaffinized by solubilization in xylene. Extracted proteins were then reduced, alkylated, and digested using filter-aided sample preparation with sequencing grade–modified porcine trypsin (Promega; ref. 39). Tryptic peptides were labeled using tandem mass tag isobaric labeling reagents (Thermo) following the manufacturer's instructions and combined into one 11-plex sample group. The labeled peptide multiplex was separated into 46 fractions on a 100 × 1.0 mm Acquity BEH C18 column (Waters) using an UltiMate 3000 UHPLC system (Thermo) with a 50 minutes gradient from 99:1 to 60:40 buffer A:B ratio under basic pH conditions, and then consolidated into 24 super-fractions. Each super-fraction was then further separated by reverse phase XSelect CSH C18 2.5 μm resin (Waters) on an in-line 150 × 0.075 mm column using an UltiMate 3000 RSLCnano system (Thermo). Peptides were eluted using a 90 minutes gradient from 98:2 to 60:40 buffer A:B ratio. Eluted peptides were ionized by electrospray (2.2 kV) followed by mass spectrometric analysis on an Orbitrap Eclipse Tribrid mass spectrometer (Thermo) using multi-notch MS3 parameters. MS data were acquired using the FTMS analyzer in top-speed profile mode at a resolution of 120,000 over a range of 375 to 1,500 m/z. Following CID activation with normalized collision energy of 35.0, MS/MS data were acquired using the ion trap analyzer in centroid mode and normal mass range. Using synchronous precursor selection, up to 10 MS/MS precursors were selected for HCD activation with normalized collision energy of 65.0, followed by acquisition of MS3 reporter ion data using the FTMS analyzer in profile mode at a resolution of 50,000 over a range of 100–500 m/z.

Protein TMT MS3 reporter ion intensity values are assessed for quality using our in-house ProteiNorm app, a user-friendly tool for a systematic evaluation of normalization methods, imputation of missing values and comparisons of different differential abundance (40). Popular normalization methods are evaluated, including log₂ normalization (Log₂), median normalization (Median), mean normalization (Mean), variance stabilizing normalization, quantile normalization (Quantile), cyclic loess normalization (Cyclic Loess), global robust linear regression normalization, and global intensity normalization (Global Intensity). The individual performance of each method can be evaluated by comparing the following metrices: total intensity, pooled intragroup coefficient of variation, pooled intragroup median absolute deviation, pooled intragroup estimate of variance, intragroup correlation, sample correlation heatmap (Pearson), and log₂-ratio distributions. The normalized data were used to perform statistical analysis using Linear Models for Microarray Data (limma) with empirical Bayes (eBayes) smoothing to the standard errors (41). We performed limma differential abundance analysis using a paired sample design to evaluate differences between injured and naïve samples. Proteins with an FDR adjusted P < 0.05 and a fold change >2 are considered significant. Significant proteins were used to identify important protein networks and pathways.

The data generated in this study are available upon request from the corresponding author.

To evaluate the sensitivity of Raman spectroscopy to the biological response to immunotherapy, we used murine tumors treated with anti–CTLA4 (αCTLA4, n = 8), anti–PD-L1 (αPDL1, n = 7) or IgG control (CTRL, n = 10; Fig. 1A). Ex vivo Raman mapping of tumors excised 3 days after completion of treatment provided a spectral dataset comprised of 7585 spectra from spatially distinct points (∼300 on an average) on the 25 tumors (Fig. 1B). The mean spectra after background subtraction for each treatment group along with the standard deviation are shown as shaded error bar plots in Fig. 2A. The background Raman signals (×2 for clarity) obtained from the quartz-aluminum plate sandwich without any tissue samples are shown in Fig. 2B. The prominent Raman peaks associated with the biological constituents of the tumors can be observed at 849 cm⁻¹ (C–O–C skeletal mode of polysaccharides), 1,260 cm⁻¹ (amide III of proteins and CH₂ in-plane deformation of lipids), 1,301 cm⁻¹ (CH₃/CH₂ twisting or bending modes of lipids and collagen), 1,448 cm⁻¹ (CH₂ bending modes in lipids and collagen), and 1,657 cm⁻¹ (amide I of proteins and C=C stretching in lipids). However, gross visual inspection does not reveal any striking differences between the average spectra across the three-treatment groups. Therefore, we used MCR-ALS to decompose the spectral dataset into pure constituent-like component loadings and their contribution (scores) in each spectrum (34). A five-component MCR-ALS decomposition provided loadings that harbored features of lipids (MC1), nucleic acids (MC2), and collagen (MC3) as shown in Fig. 2C. The remaining two loadings show weak collagen features (MC4) and features that stem from formalin fixation (MC5). The minimal spectral interference due to formalin fixation of tissues has been shown previously to not have significant impact on spectral quality (42). Furthermore, the presence of all spectral features related to formalin in a single MCR component facilitates their digital removal if necessary. The important features for each component loading are shown in the figure and their band assignments are tabulated in Supplementary Table S1 (Supplementary Information).

To evaluate changes in the biomolecular composition of tumors in response to immunotherapy, we compared the MCR scores of the pure component loadings for spectra in the three treatment groups. As shown in Fig. 2D–G, we performed three comparisons for each component to study the effect of ICIs (CTRL vs. αCTLA4, CTRL vs. αPDL1) on the distribution of their scores and to elucidate possible differences between the two therapeutics (αCTLA4 vs. αPDL1). The violin plots (with embedded box and whisker plots) in Fig. 2D show that there is a slight but significant increase (r = 0.06, P < 0.05) in the scores of lipid-like MC1 component in the αCTLA4 treatment group compared with the CTRL. In contrast, there is a much more significant decrease (r = −0.32, P < 0.05) in the lipid scores due to αPDL1 treatment. The scores of the nucleic acid-like MC2 component decreased significantly (P < 0.05) for both the treatment groups in comparison with the control group CTRL with the effect being more pronounced (r = −0.31) for the αPDL1 group (Fig. 2E). In contrast with MC1 and MC2 scores that showed larger effect sizes in the αPDL1 group, the scores of collagen-like MC3 component decreased significantly (P < 0.05) and with a larger effect size (r = −0.26) for treatment with αCTLA4, whereas they did not change significantly with αPDL1 treatment (Fig. 2F). The variation in the effect sizes (Fig. 2G) shows the differential perturbation of the TME composition upon treatment with αCTLA4 and αPDL1.

The differences in the trend and effect sizes of MCR score distributions across the three classes motivated us to explore the possibility of discriminating the treatment classes from control in the subspace formed by the pure component-like MCR loadings. We used the median MCR scores for each animal to illustrate their relative positions in the 3D space defined by MC1, MC2, and MC3 loadings (Fig. 2H). We did not find clear class boundaries between the three classes in this subspace, consistent with the subtle differences observed in the univariate comparison of the MCR scores for each component loading. We also examined the 2D projection (Fig. 2I) of the 3D data on to the MC1-MC2 plane due to the relative clustering of mice in the αPDL1 treatment group and the observed higher effect sizes for MC1 and MC2 scores in the CTRL versus αPDL1 comparison (Fig. 2G). We found a weak clustering of the CTRL and αPDL1 mice in this plane with a couple of CTRL mice residing in the relatively tighter αPDL1 cluster, consistent with the observations in univariate analysis of MCR scores. These results along with the absence of a nonoverlapping cluster for αCTLA4 group indicate that the MCR scores alone may be insufficient for assessment of response to ICIs.

Next, we sought to quantify the spatial heterogeneity of the MCR scores and assess their variation with immunotherapy. We used a metric (referred to as heterogeneity index in the ensuing discussion) developed for measuring spatial heterogeneity in medical images based on the distance-dependent deviation of the spatial distribution from the smoothest intensity gradation (36). The calculated heterogeneity indices showed significant differences for each component with αCTLA4 and αPDL1 treatment (Supplementary Fig. S1B–S1D, Supplementary Information). The spatial heterogeneity of scores of lipid-like and nucleic acid-like MCR loadings decreased with ICI therapy for both the groups. However, for the collagen-like component, we observed a decrease in the spatial heterogeneity of MCR scores for the αCTLA4 group and an increase for the αPDL1 group. These subtle differences in spatial differences show that Raman spectroscopy can not only capture variations in the compositional constituents of the tumors but also shed light on their spatial distribution in response to therapy.

Because the MCR scores plot did not provide clear class boundaries between the control tumors and tumors treated with αCTLA4 and αPDL1, we used SVMs to build a supervised classification routine for predicting response to immunotherapy. We used a leave-one-mouse-out protocol (described previously in the Materials and Methods section) to iteratively use the spectra from each individual mouse as a test dataset to the binary SVM models trained on the spectra from the remaining mice in the dataset. On the basis of the spectral prediction frequencies obtained for each class, the test mice are either classified correctly, misclassified, or unclassified according to predetermined thresholds that are detailed in the Materials and Methods section. We performed three binary comparisons to test the ability of SVM-derived classifier to predict the response to therapy using individual ICIs versus CTRL, and to distinguish between the tumors treated with αCTLA4 and αPDL1 (Fig. 3). The binary SVM classification between CTRL and αCTLA4 yielded accurate predictions for mice belonging to both the groups except two mice in the αCTLA4 group that were unclassified (Fig. 3A) whereas that between CTRL and αPDL1 groups resulted in accurate prediction of all the mice in both the classes (Fig. 3B). Finally, the comparison between αCTLA4 and αPDL1 groups to evaluate inter-therapeutic differences between the responses provided an accurate prediction of all the mice treated with αPDL1 and most of the mice treated with αCTLA4 (Fig. 3C). However, one of the eight mice in the αCTLA4 cohort was misclassified as αPDL1, whereas another was unclassified.

The availability of spectral datasets from the CT26 tumors treated with two different ICIs allowed us to use random forest classification analysis for uncovering the spectral markers of their response. We leveraged the ability of random forests to rank the predictors in the order of their contribution to the classification tasks and determined the five most important Raman predictors (spectral features) for response to each therapeutic (Fig. 4). Because we aimed to compare the differential composition of the tumors treated, respectively, with αCTLA4 and αPDL1 to CTRL, we chose the binary classification tasks'—CTRL versus αCTLA4 and CTRL versus αPDL1 for this analysis. We first tested the out-of-bag classification error for binary classification tasks using random forest training on the entire spectral dataset with bootstrapped aggregation (bagging). The out-of-bag classification error for both the comparisons decreased to less than 2% upon the inclusion of fewer than 50 trees in the forest (Fig. 4A). We also observed that the out-of-bag classification error for both the comparisons followed a similar trajectory with the inclusion of additional trees. Next, we plotted that the predictor importance estimates provided by the random forests for each comparison with uncover putative spectral markers of response to immunotherapy. We found that the primary contributors to the accuracy of random forest prediction of tumor response to αCTLA4 treatment were the Raman spectral features at ca. 981 cm⁻¹, 1,065 cm⁻¹, 1,386 cm⁻¹, 1,574 cm⁻¹, and 1,652 cm⁻¹ (Fig. 4B). The same analysis for the tumors in the αPDL1 group revealed spectral features around 1,323 cm⁻¹, 1,574 cm⁻¹, 1,596 cm⁻¹, 1,652 cm⁻¹, and 1,763 cm⁻¹ as the key markers of therapeutic response (Fig. 4C). We observed that although the two sets of spectral markers had two common features around 1,574 cm⁻¹ and 1,652 cm⁻¹, the remaining features were unique to the ICIs used for the treatment (Fig. 4D). The band assignments for the identified predictors in both the treatment classes are tabulated in Supplementary Table S2 (Supplementary Information).

The main goal of the current study was to determine whether Raman spectroscopy was sensitive to the compositional changes in the tumors in response to treatment using ICIs. The second goal was to test the specificity of these label-free measurements in distinguishing between the treatment responses mediated by distinct biological pathways. Because of their widespread use in preclinical studies of immunotherapy, we used the CT26 murine colorectal tumor xenografts in these studies to compare their responses with treatment with three doses each of anti–CTLA4 and anti–PD-L1 antibodies (43, 44). Unlike other forms of cancer treatment, immunotherapy does not result in an immediate and predictable reduction in tumor size; indeed, an initial increase in volume mediated by immune cell infiltration has often been observed (45). Here, we did not observe any profound changes in tumor volumes in response to αCTLA4 and αPDL1 treatment (Supplementary Fig. S2, Supplementary Information). In addition, histopathological assessment of the tumors did not reveal any significant differences in the levels of tumor necrosis, fibrosis, or lipids between the ICI-treated and control tumors (Fig. 5A–I). Yet, Raman spectroscopic information, as evinced by machine learning, reveals latent biomolecular changes in tumor composition in response to immunotherapy—in advance of their morphological manifestations. Specifically, the MCR-ALS analysis provided promising evidence of spectral differences that can be tied to the compositional constituents of the tumors. The subtle but statistically significant differences observed in the scores of MCR loadings resembling lipids, nucleic acids, and collagen showed that the response to αCTLA4 and αPDL1 manifests in changes in the composition of the TME. These findings are in line with prior observations that clearly demonstrated that vibrational spectroscopy is sensitive to molecular changes that precede morphological manifestations in dynamically changing TMEs (27, 46). On the basis of these findings, we sought to determine whether these changes correspond to proteomic alterations via quantitative mass spectrometry (Fig. 5J and K). Of the more than 6,600 proteins identified, about 136 proteins were found to be significantly different in treated tumors relative to controls (P < 0.05 and log₂ fold change of >2). A subset of these differentially expressed proteins (Fig. 5J and K) is known to regulate lipid metabolism and extracellular matrix composition, whereas others are known to control transcriptome dynamics or are associated with response to ICI therapy, thereby corroborating our Raman spectroscopic measurements.

The principal observations of this study are also consistent with a rapidly emerging body of work of the role of metabolism and TME in modulating immune responses. The differences in the scores of lipid-like components between the treatment classes can be attributed to the differential lipid metabolism observed in the cellular and extracellular matrix components of the TME in response to immunotherapy (24, 47, 48). For example, Lin and colleagues (49) recently observed that the tissue-resident memory T cells rely on fatty acid oxidation for survival and compete with gastric adenocarcinoma cells for lipid uptake. Immunotherapy with PD-L1 blockade enabled the memory T cells to outcompete cancer cells and contribute to antitumor immune response (49). Similarly, another study in regulatory T cells showed that the inhibition of a member of fatty acid–binding proteins (FABP5) that facilitate uptake and intracellular lipid trafficking resulted in mitochondrial changes, impaired lipid metabolism, and reduced proliferation, thus providing evidence that lipid metabolism is critical for the formation of an immunosuppressive microenvironment by regulatory T cells (50). Lipid metabolism, particularly fatty acid oxidation, is also being actively investigated in other cells that are known to play a role in immunotherapy such as tumor-associated macrophages, natural killer cells, and myeloid-derived suppressor cells (24, 48, 51–53).

Several studies have shown that the fibrotic TME protects tumors from immune surveillance and promotes resistance to therapy (22). A study in pancreatic cancer showed that the inhibition of focal adhesion kinase reduced fibrosis in the TME and improved T-cell infiltration into tumors (23). Reduction of fibrosis through CXCR4 blockade in preclinical models of metastatic breast cancer alleviated immunosuppression and significantly improved their response to immunotherapy (54). Similarly, hypoxia (which correlates with high collagen density) is known to suppress the antitumor immune response, particularly due to PD-L1 upregulation, in a variety of cancers (55). It is worth noting, though, that although most investigations have focused on the study of immunosuppressive effects of increased collagen density, they do not evaluate the changes in the collagen content in sensitive tumors in response to ICI therapy. Our results show that there is a significant decrease in the collagen-like MCR scores due to αCTLA4 treatment. However, a similar decrease in the scores of the same component due to αPDL1 was not statistically significant at the chosen threshold. We hypothesize that the observed decrease in the collagen content of the tumors after immunotherapy may be a result of extracellular matrix remodeling due to immune cell infiltration and tumor cell kill due to immune checkpoint blockade. The results of our MCR-ALS analysis and the differences in the ICI-induced spatial heterogeneity of the component scores should provide important starting points for detailed investigations to understand the underlying biological processes.

Furthermore, the results of our binary SVM-based supervised classification routine underscore the promise of an automated route for identification of response to immunotherapy. The leave-one-mouse-out protocol allows us to completely avoid the representation of spectra from test mice in the training dataset. Furthermore, the impressive accuracy of the SVM comparison between the αCTLA4 and αPDL1 groups shows that the two immunotherapy strategies result in distinct spectral responses, consistent with the significantly different mechanisms underlying their antitumor immune responses (3). This observation was also in agreement with our random forest analysis that revealed two separate sets of spectral markers with little overlap that most strongly contribute to the prediction of response to each ICI. Although the peak assignments for these predictors (Supplementary Table S2, Supplementary Information) do not readily allow the identification of independent molecular constituents, they show that multivariate analysis of Raman spectra can distinguish between treatment responses guided by distinct underlying biological mechanisms.

Taken together, our results demonstrated that molecular spectroscopic imaging is capable of noninvasive assessment of tumor response to immunotherapy. The multivariate data analysis of spectra allowed us to overcome the limitations posed by the low throughput of the method, particularly in a challenging application traditionally characterized by subtle molecular changes. However, there are some limitations in the present study that we plan to address as part of our future work involving larger animal cohorts. First, in this study, we have evaluated tumor response to monotherapy to gain a better understanding of the microenvironmental changes associated with each ICI. However, it will be important to also evaluate Raman spectral features associated with combination ICI therapy. Second, it will be important to determine whether Raman measurements made before commencing treatment are predictive of eventual treatment response, especially in combination ICI therapy. Finally, we plan to determine whether there are spectral differences between different lipid and collagen types that will allow us to better elucidate the sources of the Raman scattering signal. Nevertheless, we expect that the promising results obtained in this study coupled with the ability to readily adapt the developed experimental and computational frameworks will catalyze further optical spectroscopic investigations of immunotherapy. Future studies aimed at the prediction of tumor responses to immunotherapy using single or multiple agents and their combination with chemotherapy and radiotherapy will truly pave the way for clinical translation of Raman spectroscopy for treatment planning and monitoring applications.

S.K. Paidi reports grants from Society for Laboratory Automation and Screening during the conduct of the study. C.M. Quick reports personal fees from Allergan, Inc. outside the submitted work. N. Rajaram reports grants from National Institutes of Health and Winthrop P. Rockefeller Cancer Institute, and other support from Winthrop P. Rockefeller Cancer Institute during the conduct of the study. I. Barman reports grants from National Cancer Institute, National Institute of Biomedical Imaging and Bioengineering, and National Institute of General Medical Sciences during the conduct of the study. No disclosures were reported by the other authors.

S.K. Paidi: Data curation, formal analysis, investigation, visualization, methodology, writing–original draft. J. Rodriguez Troncoso: Data curation, formal analysis, investigation, visualization, methodology, writing–original draft. P. Raj: Data curation, formal analysis, investigation, writing–review and editing. P. Monterroso Diaz: Formal analysis, investigation, methodology. J.D. Ivers: Formal analysis, investigation, methodology. D.E. Lee: Formal analysis, investigation, methodology. N.L. Avaritt: Formal analysis, investigation, visualization, methodology, writing–review and editing. A.J. Gies: Formal analysis, investigation, methodology. C.M. Quick: Formal analysis, investigation, methodology, writing–review and editing. S.D. Byrum: Formal analysis, investigation, methodology. A.J. Tackett: Resources, supervision, funding acquisition, investigation, writing–review and editing. N. Rajaram: Conceptualization, resources, data curation, supervision, funding acquisition, investigation, methodology, writing–original draft. I. Barman: Conceptualization, resources, data curation, supervision, funding acquisition, investigation, methodology, writing–original draft.

The work was supported by the SLAS Graduate Education Fellowship grant (to S.K. Paidi). N. Rajaram acknowledges support from the National Cancer Institute (R01CA238025) and the Winthrop P. Rockefeller Cancer Institute (Team Science Award). I. Barman acknowledges support from the National Cancer Institute (R01CA238025), the National Institute of Biomedical Imaging and Bioengineering (2-P41-EB015871–31), and the National Institute of General Medical Sciences (DP2GM128198). A.J. Tackett acknowledges support from the National Institute of General Medical Sciences (R24GM137786 and P20GM121293), the National Cancer Institute (R01CA236209), the Arkansas INBRE (P20GM103429), and the Winthrop P. Rockefeller Cancer Institute (Team Science Award). The schematic in Fig. 1 was partially created with BioRender.com.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.