Introduction

Reprogramming of energy metabolism is an established hallmark of malignant tumors. Tumor cells adjust their metabolism to sustain uncontrolled proliferation and tumor progression, even in the conditions of low nutrient supply and hypoxia. There are two main metabolic pathways for energy production in the form of ATP – glycolysis and oxidative phosphorylation (OXPHOS). For a long time, enhanced glycolysis has been considered as a central feature of tumor metabolism (DeBerardinis and Chandel, 2016). Unlike most normal cells, tumor cells rely on glycolysis not only in hypoxia (anaerobic glycolysis), but also in normoxia (aerobic glycolysis, or the Warburg effect), which provides them with many metabolic intermediates and a high-speed production of ATP for fast growth. The glycolytic shift in the tumors traditionally correlates with negative prognosis (Zhou et al., 2022). At the same time, mitochondrial OXPHOS is as important in malignant cells as glycolysis. Although defective mitochondria and lower OXPHOS capacity are often observed in cancer cells, this is not an absolute phenomenon – many tumors preserve functional mitochondria and normal OXPHOS rate (Gentric et al., 2017). Along with glucose, tumors use glutamine and fatty acids as alternative substrates. Interestingly, tumor cells are capable of switching between different metabolic pathways depending on their own needs and local environment, thus demonstrating a high degree of metabolic plasticity.

Therefore, tumor metabolism is currently viewed as a dynamic, variable system with a diversity of cellular metabolic states. In this context, intratumor heterogeneity means that there are the cells with different metabolic profiles in a tumor concurrently, and intertumor heterogeneity means that tumors of the same type and stage are characterized by the different metabolic strategies (Kim and DeBerardinis, 2019).

Heterogeneity of tumor metabolism can be caused by various factors (Shirmanova et al., 2023). Some of them («intrinsic») are related to tumor cells themselves: these are histogenesis of the tumor, mutation profile, (epi)genetics factors, differentiation state and proliferation activity, to name a few (Pavlova and Thompson, 2016; Sengupta and Pratx, 2016; Seth Nanda et al., 2020; Farhadi et al., 2022). Other reasons («extrinsic») are associated with the nonuniform microenvironment, e.g., local hypoxia, nutrient distribution, heterogeneity of the vasculature network, interaction of tumor cells with the extracellular matrix and stromal cells, etc. (Marusyk and Polyak, 2010; Masson and Ratcliffe, 2014; Yoshida, 2015; Stine et al., 2015). Collectively, these factors lead to variability of tumor metabolism in space and time. Metabolic heterogeneity as a part of general heterogeneity of tumors imposes some difficulties in patients’ diagnosis and treatment and is believed to have a prognostic value.

Although metabolic heterogeneity is a well-recognized feature of tumors, its characterization at the cellular level remains scarce. In part, this is associated with the lack of the highly sensitive and direct techniques for its observation. Clinical imaging modalities, such as metabolic PET with radiolabeled tracers (¹⁸F, ¹¹C) and MRI/MRS, allow to capture inter– and intratumor heterogeneity among patients with low (∼4.5 mm) spatial resolution (Plathow and Weber, 2008). Recently evolved methods of single-cell sequencing deliver comprehensive information about the metabolic landscape of a tumor with a resolution up to 55–100 µm based on the expression of the metabolic genes, but these approaches are rather complex, laborious and not widely available (Evers et al., 2019; Huang et al., 2023).

Interrogation of cellular metabolism is possible with the laser scanning fluorescence microscopy that enables the detection of autofluorescence of the redox cofactors, such as the reduced form of the nicotinamide adenine dinucleotide (phosphate) – NAD(P)H and oxidized flavins (Georgakoudi and Quinn, 2023). Recording of the fluorescence decay parameters using the option of the fluorescence lifetime imaging (FLIM) allows the evaluation of the states of the cofactor molecules attributed to different metabolic pathways. Specifically, the unbound state of NAD(P)H has a short lifetime (∼0.4 ns) and is associated with glycolysis, while the protein-bound state has a long lifetime (∼1.7–3.0 ns) and is associated with OXPHOS (Lakowicz et al., 1992). FLIM of NAD(P)H is currently considered as a promising, label-free technique for the assessment of metabolic heterogeneity of tumors at the (sub)cellular scale. Several studies, including ours, demonstrate the possibilities of NAD(P)H FLIM to visualize metabolic heterogeneity in cultured cells, multicellular structures in vitro, tumor xenografts in vivo and patients’ tumors ex vivo (Shah et al., 2015; Skala et al., 2022; Lukina et al., 2019). However, quantifying the degree of tumor heterogeneity is lacking in most of these works.

In this paper, we present the results of assessment of cell-level metabolic heterogeneity of colorectal cancer using FLIM-microscopy of NAD(P)H. The distributions of the fluorescence decay parameters have been analyzed in the four colorectal cancer cell lines (HT29, HCT116, CaCo2, CT26), in the mouse tumor models in vivo generated from these cell lines and in the post-operative samples of patients’ colon tumors. Quantification of heterogeneity has been performed with the dispersion and the bimodality index (BI) of the decay parameter a₁ (contribution of free NAD(P)H). Correlative analysis has been conducted to find associations of the heterogeneity degree with clinical prognosis.

Results

1. Metabolic heterogeneity assessment in colorectal cancer cell lines

First, using FLIM of NAD(P)H, cellular metabolism was assessed in monolayer cultures of different colorectal cancer cell lines (Figure 1A). Typical values of NAD(P)H fluorescence lifetimes were registered for all cell lines – short τ₁ ∼0.39 ns and long τ₂ ∼2.30–2.90 ns (Table S1). Due to specifics of cellular metabolism, the relative contributions of the free (a₁) and bound (a₂) forms of NAD(P)H and, consequently, the mean lifetimes τ_m varied between the cell lines. The a₁ value decreased and τ_m increased in the following order: CT26 (∼86%, 0.57 ns) > HCT116 (∼84%, 0.71 ns) > HT29 (∼80%, 0.80 ns) > CaCo2 (∼73%, 0.94 ns) (Figure 1B, Table S1). A high NAD(P)H-a₁ (low τ_m) is an indicator of glycolytic shift, which is typical for cells with intense proliferation in a monolayer culture, like CT26 and HCT116. CaCo2 cells with the lowest a₁ value had more OXPHOS-shifted metabolism, which can be associated with the expression of morphological and functional characteristics of mature enterocytes of normal small intestine by this cell line (Lea, 2015). In all cell lines inter-cellular variations of the a₁ parameter were minor (<3%).

Figure 1.

To quantify the metabolic heterogeneity, the BI was calculated for distributions of both a₁ (Figure 1C) and τ_m values (Figure S1). For all cell lines, except CaCo2, the BI of the a₁ distribution did not exceed 1.1 (the threshold of bimodality), justifying the uniformity of cell metabolism in a culture, which is consistent with the general view on standard cell lines as homogenous populations. CaCo2 cell culture appeared bimodal in a₁ parameter (BI-a₁=1.12), which is rational because this cell line is known to be heterogeneous and contain cells of a colon carcinoma and differentiated cells with properties of mature enterocytes (Lea, 2015). The BI-τ_m was, however, rather high (>1.0) in all cell lines except HCT116 (0.84). In the further experiments on tumors and tissue samples the BI-a₁ was used as more relevant.

Additionally, the dispersion of NAD(P)H-a₁ (D-a₁) was evaluated for each cell line to describe the extent of distribution of the data around the median (Table 1). The D-a₁ value varied from 1.67% in CT26 cell culture to 3.5% in CaCo2.

Table 1.

2. Metabolic heterogeneity in mouse tumor models in vivo

Next, FLIM of NAD(P)H was performed in vivo on mouse tumor models, obtained from the colorectal cell lines shown above (Figure 2A). All the tumors were verified by histopathological analysis (Figure 2B).

Figure 2.

Among the tumors, CT26 had the highest (∼83%), and CaCo2 had, on average, the lowest (∼76%) contribution of free NAD(P)H a₁, similar to the cultured cells (Figure 2C, Table S1).

Most of the individual tumors showed larger (≥3%) inter-cellular variations of NAD(P)H-a₁ than corresponding cell lines, where the deviation from the median was ≤3%. The dispersion D-a₁ was in the range of 2.6–4.02%, and for each specific tumor type was higher than in the corresponding cell line (Table 1). At that, intertumor differences across each type of tumor were insignificant.

The BI-a₁ values in mouse tumors were generally higher than in monolayer cultures (Figure 2D, Table 1). In 4 of 12 tumors the BI-a₁ was ≥1.1, indicating the bimodal distribution of the a₁ parameter, i.e., the presence of two subpopulations of cells with different metabolism. Seven tumors had the BI-a₁ in the range of 0.7–1.0, which suggests that, while the distribution of the estimated parameter was unimodal, it was either wide or asymmetric, thus, also indicating some degree of heterogeneity. In 1 of 12 tumors the BI-a₁ was small (0.23), suggesting uniformity of cells’ metabolism. Given the genetic identity of cells in standard cell lines, we can assume that the nonuniform microenvironment in the tumors was a major source of their variable metabolism.

To ensure that the observed heterogeneity resulted from the diversity of cancer cells but not from the contribution of stromal cells, that cannot be reliably identified on NAD(P)H FLIM images, immunohistochemical staining for EpCam (epithelial marker) and vimentin (mesenchymal marker) was performed. It was found that in tumor xenografts the fraction of cells stained with vimentin was less than 10% (Figure S2).

3. Metabolic heterogeneity in colorectal cancer samples from patients

NAD(P)H FLIM images were collected from 21 postoperative samples of patients’ colorectal adenocarcinomas, among which were the tumors of the I–IV stages, poorly and highly differentiated. The representative FLIM and histological images are presented in Figure 3A and B correspondingly. Patients’ tumors showed fluorescence lifetimes of τ₁ ∼0.45 ns and τ₂ ∼1.80–3.20 ns (Table S3), which were comparable with the values in human tumor xenografts (Table S1). The parameter NAD(P)H-a₁ was in the range of ∼62–80% and NAD(P)H-τ_m was in the range of 0.80–1.20 ns, indicating that patients’ tumors generally had larger variability of metabolic statuses than cancer cells in vitro and in xenografts in vivo.

Figure 3.

A high degree of inter– and intratumor variability of cellular metabolism was detected in patients’ tumors (Figure 3C). Half of the tumors (10 of 21) showed deviation of NAD(P)H-a₁ from the median <10% across the cells, and for the rest (11 of 21) the variations were 10–25%. The dispersion D-a₁ varied significantly among the samples, from 2.19% to 11.99%.

According to the heterogeneity assessment using the bimodality index, 9 of 21 tumors were metabolically highly heterogeneous (BI-a₁≥1.1) (Figure 3D). In 11 tumors the BI-a₁ value was in the range of 0.59–1.0, which indicated the presence of metabolically different cells but not clearly separated into two clusters. Only one tumor sample (p16) was metabolically homogeneous (BI-a₁=0.24).

Notably, the bimodality index showed no correlation with dispersion. That is, among the samples there were those with bimodal distribution, but small dispersion of NAD(P)H-a₁ in a cell population (e.g., p12, p13, p21), and vice versa (e.g., p3, p7, p8).

Therefore, using NAD(P)H FLIM we have observed and quantified metabolic heterogeneity of patients’ colorectal tumors. Unlike tumor xenografts obtained from the cell lines that are thought to be genetically stable and identical, patients’ tumors are genetically diverse, which could also contribute to their metabolic heterogeneity, in addition to microenvironmental factors.

4. Interrelation between metabolic heterogeneity and clinicopathological characteristics of tumors

The relationships between the metrics of cellular heterogeneity BI-τ_m, BI-a₁ and the D-a₁ NAD(P)H with the clinical parameters of the tumor, such as the stage according to the TNM system and the grade (G), were analyzed. Due to the small sample size, all tumors were divided into two groups by each parameter: T1+T2 and T3+T4; N0M0 and all the others with metastases; G1+G2 and G3.

The results showed that the high-grade tumors displayed a higher value of dispersion of NAD(P)H-a₁ than the low-grade ones, indicating their higher heterogeneity (Figure 4). R point-biserial correlation coefficient between the grade and the dispersion a₁ was –0.49 (p-val = 0.03). Thus, the dispersion of free NAD(P)H fraction a₁ allowed differentiating between the low-grade and high-grade tumors with a high accuracy (p-val = 0.01).

Figure 4.

Tumors of the advanced stages T3 and T4 were characterized by reduced dispersion of D-a₁ and increased bimodality index BI-a₁ compared with early stages T1 and T2 (Figure 4). The stage of the tumor T has a tendency to be associated with D-a₁ (p-val = 0.06); the correlation between T stage and D-a₁ was –0.43 (p-val = 0.04).

Importantly, if metastases were present (N and M were different from 0 in TNM), then those primary tumors showed higher values of the bimodality index BI-a₁ compared with tumors for which metastases were absent (Figure 4). The NM parameter had a significant trend of association with BI-a₁ (p-val = 0.07), although the dispersion of a₁ for metastatic tumors was low.

We also note that the BI-τ_m had no significant associations with clinicopathological parameters of the tumors (Figure S3).

Discussion

The altered energy metabolism is known to support tumor progression because tumor cells are critically in need of ATP for uncontrolled proliferation and growth. It has been established that tumor cells explore different metabolic programs and utilize multiple fuels even within one tumor, thus leading to metabolic heterogeneity. Until recently, direct observation of cell-level heterogeneity of tumor metabolism was a challenging task, but became possible with the evolution of advanced optical microscopic techniques, such as two-photon fluorescence lifetime microscopy, FLIM. In the present research, using FLIM of metabolic cofactor NAD(P)H, which possesses endogenous fluorescence, we compared metabolic features of colorectal cancer cells across in vitro and in vivo models and patients’ samples. For the first time, the heterogeneity of cellular metabolism quantified on the basis of FLIM data was correlated with clinical characteristics of the tumors.

A lot of studies, including the works of our group, demonstrate the possibilities of FLIM for assessment of the intra– and intertumor heterogeneity of metabolism in different models. Using NAD(P)H FLIM, metabolic heterogeneity was revealed in colorectal cancer cell cultures obtained from the patients’ tumors, whereas the cells of the standard cell lines were uniform in their metabolism (Shirshin et al., 2022). J. Chacko and K. Eliceiri observed intercellular heterogeneity of metabolism in the MCF10A human breast epithelial cell line (Chacko and Eliceiri, 2019). Several studies show metabolic heterogeneity of 3D multicellular structures. For example, tumor spheroids obtained from the cervical cancer cell line HeLa (Lukina et al., 2018) or murine colorectal carcinoma cell line CT26 (Shirmanova et al., 2021) displayed the differences in metabolism between the outer and the inner cell layers with the outer layers being more glycolytic (higher NAD(P)H-a₁). The spheroids generated from the patients’ derived glioblastoma cultures did not show metabolic zonality, but generally had greater intercellular variations of NAD(P)H lifetime parameters than the spheroids from standard line U373 MG (Yuzhakova et al., 2023). A spectrum of works by M. Skala’ group demonstrates an opportunity of the optical metabolic imaging for heterogeneity assessment in patient tumor-derived organoids on the example of a breast, pancreatic (Sharick et al., 2020; Walsh et al., 2014; Walsh et al., 2017), head and neck (Shah et al., 2017), colorectal cancers (Skala et al., 2022).

In the research in vivo, the heterogeneous structure of tumor and its microenvironment was shown on PyMT mammary tumors (Burkel et al., 2022). Analysis of the optical metabolic parameters revealed heterogeneity of the B78 mouse melanoma model, caused by different microenvironments (Heaton et al., 2023). The authors suggested that tumor heterogeneity could be induced by such factors as pH, nutrient and oxygen availability, or activation of immune cells. In the in vivo study on HeLa tumor xenografts, metabolically different cells were detected in the cellular and collagen-rich areas (Shirmanova et al., 2018). Recently, we demonstrated a high variability of cellular metabolic statuses in CT26 tumors of large size compared with small tumors, which correlated with heterogeneous oxygen distribution (Parshina et al., 2022). A significant intertumor heterogeneity of metabolism was observed in patient-derived glioblastoma xenografts, which differed them from standard U87 glioma (Yuzhakova et al., 2022).

In our study, we compared metabolic FLIM parameters of colorectal cell lines, mouse tumor models and patients’ colorectal tumors. As expected, intercellular heterogeneity of metabolism was the highest in patients’ tumor samples as followed from a wide and often bimodal distribution of NAD(P)H-a₁. This result is consistent with the study by J. Auman and H. McLeod using genome-wide gene expression data; the authors concluded that colorectal cancer cell lines lack the molecular heterogeneity of clinical colorectal tumors (Auman and McLeod, 2010).

Cell-level metabolic heterogeneity may be present in the cell population initially, as discussed above, or develop during the treatment. The heterogeneous metabolism after treatment was detected in tumor organoids (Walsh et al., 2016; Walsh et al., 2014; Gillette et al., 2021) and in animal models (Heaster et al., 2019). In any case, the presence of metabolically distinct cells led to the varying drug responses.

An assessment of metabolic heterogeneity in patients’ tumors seems valuable from a clinical point of view. Previously, using NAD(P)H FLIM of postoperative samples, we observed a high degree of inter– and intratumor heterogeneity of T3 stage colorectal tumors in comparison with normal colon samples (Lukina et al., 2019). Unfortunately, most of the tumor types are inaccessible for microscopic investigation in patients using clinically approved techniques. The exception is the skin tumors that can be inspected with clinical systems, like MPTflex or MPTcompact (JenLab, Germany). At the same time, endoscopic FLIM-microscopy has been developing actively since the last ten years, which opens the prospects for in vivo examination of a wide spectrum of patients’ tumor types (Sun et al., 2013; Farhadi et al., 2022). In principle, an assessment of tumor metabolism in biopsy samples is also possible, which widens the potential clinical applications of FLIM.

Although metabolic heterogeneity at the cellular level has been reported for different models and patient-derived cells, the comparisons between different samples and studies are hard to make as most of them lack any quantification of the heterogeneity.

Different approaches have been proposed to quantify the metabolic differences identified by FLIM. For example, in the papers (Shah et al., 2015; Sharick et al., 2019) a heterogeneity index and its modified form the weighted heterogeneity index (wH-index) were used for quantitative analysis of cellular heterogeneity. The wH-index is based on the Gaussian distribution models and is a modified form of the Shannon diversity index. Using this index the authors described the variations of a combination of parameters, such as NAD(P)H fluorescence lifetime, FAD fluorescence lifetime and optical redox ratio, defined as the OMI-index. The methods developed in Ref. (Heaster et al., 2019) established the combination of optical metabolic imaging variables and spatial statistical analysis (spatial proximity analysis, spatial clustering, multivariate spatial analysis, spatial principal components analysis) to quantify the spatial heterogeneity of tumor cell metabolism. Recently, we have proposed a new quantitative criterion – the bimodality index (BI) – to accurately discriminate between metabolically diverse cellular subpopulations on the basis of NAD(P)H FLIM data (Shirshin et al., 2022). The BI provides dimensionless estimation on the inherent heterogeneity of a sample by checking the hypothesis about approximation of a fluorescence decay parameter distribution by two Gaussians. Using the BI, the metabolic heterogeneity has been identified in standard and primary cancer cells cultures after chemotherapy with 5-fluorouracil. Here, we used the dispersion of NAD(P)H-a₁ and BI-a₁ for quantifying metabolic heterogeneity of patient’s tumors and compared it with in vitro and in vivo models.

Metabolic heterogeneity of colorectal cancer is also discussed in the literature. The focus of many of these studies is on the molecular classification of tumors by the analysis of their metabolic features. For example, Zhang et al. showed that colorectal tumors could be classified into three distinct metabolism-relevant subtypes and developed a metabolism-related signature consisting of 27 metabolic genes, which were expressed differentially among the three subtypes and correlated with patients’ overall survival (Zhang et al., 2020). Varshavi et al. performed metabolic characterization of colorectal cancer cells depending on KRAS mutation status using ¹H NMR spectra of the metabolites. They revealed that some mutations led to an increase of glucose consumption and lactate release, while others, on the contrary, decreased it (Varshavi et al., 2020). Liu et al. found a correlation between intratumor metabolic heterogeneity parameters of ¹⁸F-FDG PET/CT and KRAS mutation status in colorectal cancer – KRAS mutant tumors had more ¹⁸F-FDG uptake and heterogeneity than wild-type KRAS (Liu, Wang et al., 2022). In a different study they showed that intratumor metabolic heterogeneity assessed from ¹⁸F-FDG PET/CT is an important prognostic factor for progression-free survival and overall survival in patients with colorectal cancer (Liu, Xiang et al., 2022). The value of intratumoral metabolic heterogeneity in ¹⁸F-FDG PET/CT for prediction of recurrence in patients with locally advanced colorectal cancer was also demonstrated by Han et al. (Han et al., 2018). In the study by Zhang et al. intratumoral metabolic heterogeneity derived from ¹⁸F-FDG PET/CT was higher in colorectal tumors with a high microsatellite instability (MSI) – an important prognostic biomarker, and predicted MSI in stage I–III colorectal cancer patients preoperatively (Zhang et al., 2023). Lin et al. identified two subtypes of colorectal cancer using a metabolic risk score based on genes that were mostly involved in lipid metabolism pathways; this criterion was applied for survival prediction – patients with a higher metabolic risk score had worse prognosis (Lin et al., 2021). Therefore, these clinical studies consider intratumor metabolic heterogeneity as a useful prognostic factor. Our results on patients’ colorectal tumors revealed some associations between cell-level heterogeneity and tumor differentiation grade, a stage-independent prognostic factor in colorectal cancer where the lower grade correlates with better the prognosis. The high-grade (G3) tumors had significantly higher heterogeneity, quantified by the dispersion parameter D NAD(P)H-a₁, compared with low-grade ones.

Another important observation in our research is that bimodality in NAD(P)H-a₁ distribution in primary tumors was linked to the presence of metastases. Several studies report on the correlations of colon tumor heterogeneity with metastases status. For example, in Ref. (Rajput et al., 2017), analysis of the mutational spectra of stage II and III colon cancer samples using next-generation sequencing approach revealed that there was a strong correlation between higher MATH (the Mutant-Allele Tumor Heterogeneity) score and a risk of metastases. Oh et al. observed a significant difference in progression-free survival according to heterogeneity assessment from targeted deep sequencing data, suggesting that tumor heterogeneity might be a determining factor for recurrence or metastasis in colorectal cancer patients with curative resection (Oh et al., 2019). Joung et al., using the whole-exome sequencing, determined that heterogeneity of primary colorectal tumors can predict the potential for liver metastasis (Joung et al., 2017). Similar results were obtained for different types of tumors, including breast (Yang et al., 2017; Chung et al., 2017) and cervical cancers (Kidd and Grigsby, 2008; Pinho et al., 2020).

In the context of metabolic heterogeneity assay using NAD(P)H FLIM, some limitations should be mentioned. If the investigation of patients’ tumor metabolism is performed on ex vivo tissue samples, one should be careful to work with freshly excised tissue or store the sample in the appropriate conditions. As we found previously, FLIM parameters of autofluorescence change very quickly (within 15 min) on air, but can be preserved in 10% BSA on ice during 3 h (Lukina et al., 2019). A limitation of our study is that microscopic FLIM images of the tissues were processed manually to obtain information about each individual cell, which made the analysis time-consuming. The methods for segmentation of tissue images have been continuously developed in digital histopathology. As for FLIM, there are approaches to automatic segmentation of FLIM images of cell cultures on the basis of machine learning. So, in spite of the complex and heterogeneous tissue architecture, there are expectations that the task of identification of the individual cells in FLIM images of tissues will be solved in the nearest future. Finally, a small size of the microscopic field of view (typically less 300 µm) limits the inspected area within the tumor sample. Taking into account a high spatial heterogeneity of clinical tumors, there is no confidence whether the inspected area is sufficiently representative.

Owing to the high (sub)cellular resolution (∼200–500 nm), FLIM-microscopy provides unique information about metabolic heterogeneity, unavailable with any other methods, like ¹⁸F-FDG PET or spatial transcriptomics. Therefore, FLIM-microscopy of endogenous cofactors is not only a powerful research technique, which is capable of improving our understanding of cellular metabolic diversity, but also the tool with a great potential for clinical translation to predict disease outcome.

Conclusions

It is now evident that metabolic processes in cancer are highly variable, which make the tumor metabolism extremely heterogeneous. Metabolic heterogeneity of tumors complicates treatment efforts and is thought to be a negative prognostic factor. Due to the lack of methods for direct observation and quantification of metabolic heterogeneity at cellular level, it has been poorly characterized so far. Our assessments of cell-level heterogeneity from FLIM-microscopy of NAD(P)H clearly show that heterogeneous metabolic landscape of a patient’ tumor is hardly reproducible in in vitro and in vivo models, which underlie the importance of such investigations on clinical material. Although the present research included a limited number of patients (n=21), the obtained results showing the associations between metabolic heterogeneity and clinical features of tumors allow us to consider it as a potential prognostic marker. We plan to continue this study to collect more tumor samples and estimate the correlations between their metabolic heterogeneity and follow-up clinical outcomes.

Materials and methods

Cell Cultures

The human colorectal carcinoma cell lines HТ29, HCT116, CaCo2 and murine colon carcinoma cell line CT26 were used in the study. Cells were grown in Dulbecco’s modified Eagle’s medium (DMEM) (Gibco, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (FBS) (Gibco, Carlsbad, CA, USA), 2 mM glutamine (Gibco, Carlsbad, CA, USA), 10 mg/mL penicillin (Gibco, Carlsbad, CA, USA), 10 mg/mL streptomycin (Gibco, Carlsbad, CA, USA). The cells were incubated at 37 °C, 5% CO2, and 80% relative humidity and passaged three times a week. The passaging of cells was carried out at a confluence of 70– 80% with trypsin-EDTA (Thermo Fisher Scientific, Waltham, MA, USA).

For fluorescence microscopic studies the cells were seeded in 35 mm glass-bottomed FluoroDishes (Ibidi GmbH, Gräfelfing, Germany) in the amount of 5х10⁵ cells in 2 mL of DMEM and incubated for 24 h (37°C, 5% CO2). Before FLIM, DMEM was replaced with the FluoroBrite DMEM (Thermo Fisher Scientific, Waltham, MA, USA) for fluorescence imaging.

Tumor Models

In vivo experiments were performed on female nude and Balb/C mice weighing ∼20–22 g. All animal protocols were approved by the Local Ethical Committee of the Privolzhsky Research Medical University (Nizhny Novgorod, Russia). To generate tumors, the suspensions of cancer cells were injected subcutaneously into the tight in the following amount: 1.5×10⁶ of НТ29 cells in 100 µL PBS, 1.0×10⁶ of НCT116 cells in 100 µL PBS, 5.5×10⁶ of CaCo2 cells in 100 µL PBS. CT26 cells (2.0×10⁵ cells in 20 µL PBS) were inoculated intracutaneously in the ear. The tumor volume was measured using a caliper, and calculated using the formula v=a×b×1/2b, where a is the length and b is the width of the tumor. In vivo studies were done on 21^st day of growth for HT29 tumors (244.8±36.6 mm³), on the day 23 for HCT116 (533.9±42.0 mm³), on the day 68 for CaCo2 (397.8±65.2 mm³), and on the day 14 for CT26 (60.9±5.6 mm³).

Mice were anesthetized by an intramuscular injection of Zoletil (40 mg/kg, Virbac SA, France) and 2% Rometar (10 mg/kg, Spofa, Czech Republic) for intravital microscopy. The skin over the tumor was opened and covered with a coverslip.

Patient samples

Twenty-one surgical samples of patients’ colorectal tumors were obtained in the Nizhny Novgorod Regional Oncological Center (Nizhny Novgorod, Russia) during the tumor resection. All the patients gave informed written consent prior to the enrollment in the study.

All the patients had a histological verification of colorectal adenocarcinoma, the stage definition according to the TNM system, and the definition of the grade. There were tumors of the I–IV stage of low (G1), moderate (G2) and high (G3) grade. The tumors were localized in different sites of the large intestine (caecum, colon, rectum). 19 of 21 patients did not receive any anti-cancer therapy before the surgery, 2 patients (7 and 8) have been pretreated with radiotherapy or chemotherapy. The detailed information about patients’ tumors is presented in Table 2, data about their clinicopathological characteristics is presented in Table S2.

Table 2.

Immediately after surgical excision, tumor samples, 0.5–1 cm³ in size, were wrapped in gauze soaked in a solution of 10% BSA (bovine serum albumin), placed in sterile Petri dish on ice, transferred to the laboratory within 30 min and examined on the laser scanning microscope immediately. This storage protocol allows to preserve autofluorescence lifetime parameters unchanged for at least 3 hours (Lukina et al., 2019).

FLIM-microscopy

FLIM of NAD(P)H was performed using the laser scanning microscope LSM 880 (Carl Zeiss, Germany) equipped with a TCSPC-based FLIM module (Becker & Hickl GmbH, Germany). The Ti:Sa femtosecond laser MaiTai HP (Spectra-Physics Inc., USA, repetition rate 80 MHz, pulse duration 120 fs) was used for two-photon excitation of NAD(P)H at a wavelength of 750 nm. Fluorescence signal was registered in the range 450–490 nm. The laser power applied to the samples was ∼6 mW. FLIM images were obtained using the water-immersion objective C-Apochromat W Korr 40×/1.3 (Carl Zeiss, Germany). Image collection time was 60 s. In total, 5–10 images were obtained from each sample.

FLIM images were processed in the SPCImage 8.5 software (Becker & Hickl GmbH, Germany). NAD(P)H fluorescence was analyzed in the cytoplasm of individual cells, in total 50-200 cells in each sample. Fluorescence decay curves were fitted by a bi-exponential model with a goodness of fit χ² 0.8–1.2. The following fluorescence decay parameters were analyzed: short component corresponding to the free form of NAD(P)H (τ₁), long component corresponding to the protein-bound NAD(P)H (τ₂), their relative contributions (а₁ и а₂, correspondingly, a₁+a₂=100%), and the mean fluorescence lifetime

Histopathology and Immunohistochemistry (IHC)

Formalin-fixed tumor samples were embedded in paraffin according to standard procedure and cut parallel to the optical plane. The sequential sections 7 μm thick were stained with hematoxylin and eosin, sections 4 µm thick were used for immunohistochemical staining. Sections were stained with primary antibodies to EpCAM (ThermoFisher Scientific, mouse, cat. #14-9326-82, dilution 1:700) and vimentin (Abcam, ab137321, rabbit, dilution 1:700), according to the manufacturer’s protocol. Secondary antibodies goat anti-mouse (Alexa Fluor 488) and goat anti-rabbit (Alexa Fluor 555) were used. Tissue sections were imaged using EVOS M7000 Imaging System (Thermo Fisher Scientific Inc., Waltham, MA USA) with LED cubes GFP (ex.470/22 nm, em.510/42 nm for Alexa 488) and RFP (ex.531/40 nm, em. 593/40 nm for Alexa 555) at x20 magnification.

Statistical Analysis

The obtained data were checked for the normal distribution using the Kolmogorov–Smirnov’s criterion. The data with a normal distribution were presented as mean ± standard deviation (SD). The data with an abnormal distribution were presented as a median and 25% and 75% quantiles (Q₁ and Q₃). The ANOVA test for comparison data with a normal distribution and the Kruskal-Wallis’s test for comparison data with an abnormal distribution were used, with p-val<0.05 being considered statistically significant. Statistical data processing was performed in IBM SPSS Statistics 26.0, R-studio, Python.

Bimodality Index and Dispersity calculation

A continuous measure known as the «bimodality index» is utilized to gauge the degree of conformity of a set of univariate data to a two-component mixture model. The score is larger if the two components are balanced in size or if the separation between the two modes is larger. The BI≥1.1 is considered as a cutoff to reliably define a bimodal distribution in the data (Wang et al., 2009; Shirshin et al., 2022).

BI was calculated according to the equation:

where μ – the mean of each Gaussian, σ – the standard deviation of each Gaussian and n – the number of measurements of each Gaussian.

Dispersion was calculated as interquartile range: 75% quantile – 25% quantile (Q₃ – Q₁).

The parameters τ_m and a₁ obtained by FLIM microscopy were used to calculate the bimodality indices (BI) and the dispersion for cultured cells, mouse tumors and patients’ samples.

The resulting BI for τ_m and a₁ and dispersion of a₁ were used to search for the relationships between them and the pathophysiological parameters of tumors: the TNM stage and the grade (G).

The fluorescence decay parameters were normalized using the z-normalization method to level out the numerical difference.

A random forest decision tree model was used to evaluate the importance of fluorescence decay parameters on clinicopathological characteristics. The outcome of the model was used for SHAP analysis. Finally, SHAP analysis allowed us to determine the most significant parameters that were investigated by the classic statistical methods (Mann–Whitney).

SHAP (SHapley Additive exPlanations) analysis is a method used to explain the output of any machine learning model. SHAP analysis portrays the relative impact of variables on the output of the model, and it can give insights into the most significant factors and the influence of variables on the result. SHAP values are a way to explain the contribution of each feature to the model’s output. They measure how much each feature contributes to the model’s prediction and can help identify which features are most important for the model and how they affect the outcome.

The T parameter from TNM was modified in a way that classes 1+2 became group 0, and 3+4 became group 1. The modification of the parameter G resulted in two classes: classes 1+2 became group 0 and class 3 became group 1. Parameters N and M were modified as NM, in which the class was distinguished: no metastases – 0, metastases – 1. Since the clinical parameters, in a modified form, are binary, a point-biserial correlation coefficient was chosen to assess the correlation.

Acknowledgements

The authors are grateful to Dr. Anton Plekhanov (PRMU) for helpful discussions and valuable suggestions.

Funding

This work was supported by the Russian Science Foundation (Project No. 23-15-00294).

Competing interests

The authors declare no conflict of interest.

Supplementary

Figure S1.

Figure S2.

Figure S3.

Table S1.

Table S2.

Table S3.