Background

Despite improvements with BRAF and immune-checkpoint inhibitors (ICI), at least half of patients with melanoma succumb to disease.1 2 Tissue-based biomarker testing requires invasive procedures and inadequately describes response. Meanwhile, interlesion metastatic heterogeneity impacts the survival of ICI and is influenced by the site of metastasis.3–10

Factors intrinsic to tumor cells, such as tumor mutational burden and immunogenicity, among others, vary by anatomic location of metastasis and impact therapy.11–13 In melanoma and non-small cell lung cancer (NSCLC), metastasis-specific patterns impact the treatment outcome of anti-programmed cell death protein-1 (PD-1).3–5 Interlesion heterogeneity of immune cell content has also been described in ovarian and colorectal cancers where non-responding lesions have been associated with immune exclusion.14 15 Particularly myeloid cells are known to have a complex phenotypic distribution across organs16 and impact on antitumor immunity.17 18 Local therapies for oligometastatic progression to extend BRAF or immunotherapy benefit is a clinical strategy in melanoma oncology.19 20 An unmet need remains for non-invasive biomarkers of lesion and organ-specific response to inform treatment selection and the study of tumor resistance.

Multiple previous studies have used radiomics to predict ICI responses in melanoma, however, these have focused on anti-PD-1 or analyzed only single lesions or a predefined set of features to predict patient outcomes.21–25 Few studies have focused on identifying features specific to metastatic sites. We report patient and organ-level response patterns as well as radiomic models across treatment modalities of ipilimumab plus nivolumab, anti-PD-1 monotherapy, and BRAF±MEK inhibitor (MEKi) combinations in advanced melanoma.

Methods

A full description is supplied in the online supplemental methods. Patients with unresectable stage III/IV melanoma who received ICI (anti-PD-1/cytotoxic T-lymphocyte associated protein 4 (nivolumab plus ipilimumab), referred to as I+N cohort, n=91; anti-PD-1, referred to as PD-1 cohort, n=151; 242 total), or BRAF±MEKi targeted therapy (BRAF cohort, n=49) from 2015 to 2020 were identified from UPMC Hillman Cancer Center registry (online supplemental figure S1, table S1). Best overall response and organ-specific tumor responses (adrenal, brain, liver, lymph node (LN), lung, and soft tissue) were evaluated by Response Evaluation Criteria in Solid Tumors (RECIST) V.1.1. Across all cohorts, 1,166 CT and 168 MRI scans at baseline and best response were analyzed. After quality control, scans were segmented by 3D Slicer, with 400 radiomics features extracted,26 including 10 first-order (FOF001-010), 195 second-order (SOF001-195), and 195 volume-adjusted second-order features (SOVF001-195). Machine-learning (ML) models by XGBoost were constructed to predict two classes: disease control (DC, including complete response (CR), partial response (PR), stable disease (SD)) or progressive disease (PD). Models were developed to predict overall response or organ-specific response. Two types of ML models (radiomics features only; radiomic plus clinical features) were developed when appropriate. Statistical analysis was performed using R (V.4.1.2) and Bioconductor (release 3.14), with false discovery rate (FDR) controlled at 0.10, and Benjamini-Hochberg-FDR adjustment for multiple comparisons. An analysis overview is provided in online supplemental figure S2, with the 400 radiomics features described in online supplemental tables S2 and S3.

ResultsPopulation characteristics and overall response

Among 242 patients receiving ICI (I+N and PD-1), 199 (82%) received first-line immunotherapy (online supplemental table S1). The I+N and BRAF cohorts had higher combined incidences of M1C (extrapulmonary metastasis) and M1D (central nervous system metastasis) by American Joint Cancer Classification (AJCC) staging criteria (69% and 70%, respectively), while the PD-1 cohort showed similar distribution across M1A-D. In all cohorts, the majority of patients exhibited high levels of circulating lactate dehydrogenase (LDH>upper limit normal (ULN), the institutional ULN of LDH is 170) and elevated neutrophil-to-lymphocyte ratio (NLR>3.0) at baseline. DC of 53% (8% CR, 25% PR, 20% SD) was observed in all ICI (figure 1A). Cohort review included I+N, PD-1, and BRAF, respectively, with 41%, 60%, and 86% DC (figure 1A). DC varied by line of therapy with I+N (45% first line vs 31% in ≥2 lines) and PD-1 (62% first line vs 41% in ≥2 lines; figure 1B), consistent with the literature.27

Figure 1

Overall and organ-specific response in ICI (I+N, PD-1) and BRAF cohorts. (A) Overall response to therapy by treatment. For each cohort, the outer circle shows the percentage of patients who experienced CR, PR, SD, or PD. The inner circle shows DC% (including CR, PR, and SD) and PD%. (B) Patients stratified by previous exposure to immunotherapy. Color represents CR, PR, SD, and PD same as in (A). The number above each bar shows the percentage of patients who experienced DC in each subset. (C) Heatmap showing the organ-specific response (adrenal, brain, liver, LN, lung, soft tissue, on the row) and in the context of overall response per patient (on the column). n=291 patients are shown in (A), (B), and (C), with 91 from I+N, 151 from PD-1, and 49 from BRAF. (D) Comparison of interorgan heterogeneity in patients with mixed response versus those with uniform progression versus those with uniform disease control. The y-axis represents the SD of weighted RECIST scores across all metastases in a patient. Each data point represents one patient. n=111 ICI patients who had at least two metastasis sites are shown. (E) Intraorgan heterogeneity by organ site comparing lesions 01 and 02. The y-axis represents the individual lesion’s tumor size change in percentage. Each data point represents one lesion. Lines connect lesion 01/02 from the same metastasis site in the same patient. n=168 sites from ICI patients who had two lesions per site are shown. Wilcoxon rank-sum test was used in (D), Wilcoxon signed-rank test was used in (E). FDR-adjusted p values are shown in (D) and (E). FDR was controlled at 0.10. All tests are two-sided. Denotations: **p<0.01, *p<0.05, + p<0.10. CR, complete response; DC, disease control; FDR, false discovery rate; ICI, immune-checkpoint inhibitors; I+N, ipilimumab+nivolumab; Imtx, immunotherapy; LN, lymph node; PD, progressive disease; PD-1, programmed cell death protein-1; PR, partial response; SD, stable disease.

Mixed organ responses are associated with a higher risk of progression at the patient level

Patients with mixed response between metastatic sites had a significantly higher likelihood of PD, suggesting heterogeneity in organ-specific response determines overall response (p<0.0001 in all cohorts, two-sided Fisher’s exact test; figure 1C). Organ-specific response was defined by weighted RECIST (referred to as, RECISTweighted), which accounts for baseline size and lesion size changes at best response compared with baseline, particularly when multiple lesions were identified in a single organ (see online supplemental methods). Patients were categorized into three groups: uniform PD (all organs showing progression), mixed response (some organs were DC while others showed PD), and uniform DC (all organs showing DC). Patients with uniform PD exhibited significantly higher variation in RECISTweighted compared with those with uniform DC (FDR-adjusted p=0.0062; two-sided Wilcoxon test; figure 1D). However, patients with mixed response showed similar variation in RECISTweighted versus uniform PD (p=0.27), but significantly higher variation than uniform DC (p=0.021; figure 1D). To address potential confounding by higher means, we compared the SD of absolute and original RECISTweighted scores within each group, showing that mixed directionality significantly contributed to response heterogeneity (p=0.000051; online supplemental figure S3A), with similar findings in patients with uniform PD (p=0.014) and a more homogeneous pattern in uniform DC (p=0.70). Those findings suggested that among patients with some levels of organ response, those with higher organ response heterogeneity are more likely to progress.

Disease control and resistance to ICI or BRAF differ between and within metastatic sites

Tumor growth or reduction across cohorts was stratified by the best response at the patient level comparing DC/PD. Among all ICI experiencing DC, liver (n=46) and lung (n=92) metastases experienced greatest reduction in tumor (−66%±8% and −63%±5%, respectively; mean±SEM; online supplemental figure S3B; I+N and PD-1 each are shown in online supplemental figure S3C and D). For BRAF, liver metastases experienced the greatest reduction (−58%±8%) (online supplemental figure S3E). Additional tumor volume metrics are provided in online supplemental table S4.

To assess intraorgan heterogeneity of response to ICI, tumor volume changes are reported in patients with multiple lesions in the same organ. Among patients with multiple same-organ target lesions, brain metastases demonstrated the highest variability in intraorgan response comparing RECISTweighted of lesions 01 versus 02 (FDR-adjusted p=0.084, two-sided Wilcoxon test for paired samples) (figure 1E). For all organ sites, the largest lesion at baseline was designated as 01 and second as 02. In addition, the absolute differences in the two lesions were compared between overall response groups. Lung metastases had greater intraorgan absolute differences in DC/PD (FDR-adjusted p=0.074; two-sided Wilcoxon test; online supplemental figure S4A). Two cases involving liver and lung metastases (red arrows, online supplemental figure S4A) that experienced the largest differences in lesion 01 versus 02 tumor size change after ICI treatment are highlighted in online supplemental figure S4B and C. No substantial interlesion heterogeneity was found across metastatic sites for BRAF (online supplemental figure S4D and E). Adding clinical variables, we observed that non-cutaneous melanoma showed significantly higher RECISTweighted than cutaneous melanoma in liver metastases from the ICI cohort (online supplemental figure S5).

Radiomic signature differentiates disease control rate at the patient and organ levels after BRAF-targeted and immunotherapy

To identify radiomic features associated with DC or PD, we performed comparisons between DC/PD at the patient level, or within each organ. Starting with I+N, we detected 39 differential radiomic features at FDR 0.10 (figure 2A). Of these, we found that the DC/PD differences were mostly driven by lung compared with other organ sites (online supplemental figure S6) however this was not shared by PD-1 or BRAF.

Figure 2

Radiomic features differentiate overall response or organ-specific response DC versus PD. (A) 39 features that distinguish overall response in patients who received I+N at FDR-adjusted p<0.10. Patients are clustered on the column and features are clustered on the row with dendrograms shown. The horizontal annotation bar on top of the heatmap indicates the overall response PD and DC. Feature names (eg, IMC_1_variance) are shown on the right side of the heatmap, which correspond to the feature IDs (eg, SOF194, SOF038). n=82 patients from I+N cohort are shown. (B) Overlapping or unique features across patient cohorts. The DC versus PD differences of the 39 features from (A) (I+N) are shown in patients who received PD-1 monotherapy or BRAF-targeted therapy. Features are shown in the same order as on the heatmap from (A). (C) 14 features that distinguish organ-specific response in one or more cohorts (left to right) each organ site is shown in I+N or PD-1: lung, LN, liver, and soft tissue. For (A) and (C), full feature IDs and names are described in online supplemental tables S2 and S3. Wilcoxon rank-sum test was used in (A), (B), and (C). All tests are two-sided. DC, disease control; FDR, false discovery rate; I+N, ipilimumab+nivolumab; IMC, informational measure of correlation; LN, lymph node; PD, progressive disease; PD-1, programmed cell death protein-1; FOF, first-order feature; SOF, second-order feature; SOVF, volume-adjusted second-order feature.

Across the three cohorts, the I+N features were partially detected comparing DC with PD in the PD-1 cohort, but absent in BRAF (figure 2B). On comparing DC/PD within each organ, we found that I+N and PD-1 showed non-overlapping features associated with organ level response (figure 2C; nominal p<0.01). In lung metastases, seven radiomic features of “correlation average”, “correlation range”, “informational measure of correlation [IMC]_1 range”, and “IMC_1 variance” are associated with outcome in I+N, which were not detected in PD-1 (figure 2C, left). At LN and liver metastases, I+N cohort did not return significant features, while the PD-1 cohort showed differences in “IMC_2 average” (figure 2C, middle). In soft tissue metastases, two features “minimum” and “percentile_1” which are first-order radiomics features distinguish DC from PD organs in I+N cohort, and no differences were detected in PD-1 (figure 2C, right). Collectively, our results suggested that binary group comparisons distinguish radiomics features in I+N cohort with organ heterogeneity, and organ-specific features were specific to I+N or PD-1 cohort each.

Radiomics features predict response at the patient level and integrating clinical features improves the performance of DC versus PD classification

To predict the overall response across all metastatic sites per individual, we evaluated 221 ICI patients with high-quality images and performed radiomics modeling (86 I+N and 135 PD-1). BRAF cohort was not included due to sample size. MRIs were only used for brain metastases, while CT scans were used for other organ metastases. We constructed machine-learning (ML) models classifying DC/PD within I+N and PD-1 each using XGBoost by 80/20 training/test with 10-fold cross-validation (CV) for feature selection and model optimization. Training/test set samples were randomly split prior to any processing steps, ensuring that samples in the test set were blind from the training set. In the training set, after proper scaling and normalization, features that are highly correlated (Spearman’s correlation ρ>0.80), show high collinearity, or low variance were removed before model construction. This step reduced the number of radiomic features from 400 to approximately 10–30 total, which were then used for feature selection and hyperparameter tuning via 10-fold CV in the training set. The performance of the final optimized model on unseen data is reported in the test set.

The training CV area under curve (AUC) and test AUC are shown in figure 3. XGBoost models consisting of optimized radiomic signatures classified DC from PD across I+N (AUC=0.71, CV; 0.85, test set; figure 3A, online supplemental figure S7A), PD-1 (AUC=0.65, CV; 0.71, test set; figure 3B, online supplemental figure S7B). Integration of clinical variables (age, body mass index (BMI), sex, AJCC stage, melanoma subtype, baseline eosinophil count, LDH, and NLR) to the radiomic signature achieved an improved predictive model for I+N (AUC=0.71, CV; 0.89, test set; figure 3C), PD-1 (AUC=0.70, CV; 0.90, test set; figure 3D). We found that BMI, albumin, age, and eosinophil count were important predictors for I+N (VarImp≥20; figure 3E). These four, plus NLR and AJCC M1c stage, also ranked as the top variables for PD-1 model (figure 3F). While our models incorporated volume-adjusted second-order radiomics features (SOVF001-195; see online supplemental methods), we also explored simpler models based on lesion volume to assess their predictive capability. However, the volume-only patient models showed poor performance compared with radiomic or radiomic+clinical models (online supplemental figure S7C and D). The lesion volume model’s performance was evaluated using receiver operating characteristic (ROC)-AUC analysis. Variable importance was not computed for the volume-only model as it was not applicable to single-predictor models. ROC-AUC provides a direct comparison of the overall prediction capabilities of radiomic versus volume-only models. Taken together, these results demonstrate that integrating clinical variables with radiomic features significantly enhances the predictive performance for classifying DC versus progression.

Figure 3

Radiomics models predict overall response DC/PD in ICI cohorts. For each cohort, models were optimized in the training set with 10-fold CV, and the final performance was reported on unseen data in the test set. We show both the training set 10-fold CV ROC curve as well as the test set ROC curve. AUC, sensitivity (Sens), specificity (Spec), precision (Prec), and F-score (F1) were reported. (A) Model of radiomic features only in I+N cohort. (B) Model of radiomic features only in the PD-1 cohort. (C) Model of radiomic features and clinical variables in I+N cohort. (D) Model of radiomic features and clinical variables in PD-1 cohort. For I+N models in (A) and (C): n=67 and 15 patients in training/test set (80%/20% split), respectively (total is 82). 400 radiomic features were reduced to 17 prior to model training. For PD-1 models in (B) and (D): n=104 and 25 patients in training/test set (80%/20% split), respectively (total is 129). 400 radiomic features were reduced to 23 prior to model training. (E) Variable importance (VarImp) of the features from I+N model in (C). (F) Variable importance (VarImp) of the features from the PD-1 model in D). Features with VarImp>1 are shown in (E) and (F); red vertical dashed line indicates VarImp=20; features with VarImp≥20 are generally considered important in predicting outcome. Color indicates whether a feature is greater in overall response PD (blue) or DC (gold). Clinical variables are bolded. The AUC p value shown at the top left corner of each ROC panel in (A-D) was computed using function roc.area from R package verification (V.1.42), which implements a two-sided Wilcoxon rank-sum test. ASM, angular s moment; AUC, area under the curve; BMI, body mass index; CV, cross-validation; DC, disease control; FPR, false positive rate; ICI, immune-checkpoint inhibitors; I+N, ipilimumab+nivolumab; IMC, informational measure of correlation; NLR, neutrophil-to-lymphocyte ratio; PD, progressive disease; PD-1, programmed cell death protein-1; ROC, receiver operating characteristic; TPR, true positive rate.

Radiomics features predictive of organ-specific response vary by site of metastasis in ICI cohorts

We constructed radiomics models to predict response at metastasis level for I+N or PD-1. To maximize the number of samples for model construction, we included all samples for training with leave-one-out cross-validation (LOOCV). Features were selected based on a bootstrapping strategy using 80% of the samples, repeated 100 times. In each bootstrapped set, we compared features between DC/PD within each organ and selected the top 10 features that consistently passed a lenient p value threshold (nominal p<0.20) across bootstrapped sets (see online supplemental methods). The performance of the final model was reported based on LOOCV from the training set, recognizing that the generalization of the models will require independent validation. From I+N, the best performance of organ-specific models was observed in LN (AUC=0.85; figure 4A, online supplemental figure S8A). For PD-1, the best performance of organ-specific models was observed in liver (0.94; figure 4B, online supplemental figure S8B). In volume-only organ-specific models, I+N lung and PD-1 soft tissue models achieved AUCs of 0.69 and 0.72, respectively, while the rest showed AUCs ranging from 0.50 to 0.61 (online supplemental figure S8C and D).

Figure 4

Radiomics models predict organ-specific response DC/PD in immune-checkpoint inhibitors cohorts. For each cohort, leave-one-out cross-validation was applied to all samples to generate the ROC curve. AUC, sensitivity (Sens), and specificity (Spec) were reported. 400 radiomic features were reduced to 10 prior to model construction. (A) Model of radiomic features in I+N cohort (left to right): lung (34), LN (37), liver (21), soft tissue (32), and brain (20). Numbers in paratheses indicate the number of metastases for model construction per organ. (B) Model of radiomic features in PD-1 cohort (left to right): lung (54), LN (52), liver (22), soft tissue (50). Brain models of the PD-1 cohort were not constructed considering the small sample size. (C) Heatmap summarizing radiomic features used by overall response or organ-specific response models to predict DC/PD. Models are shown on the row, and features are shown on the column. Same radiomic features at different gray levels were collapsed as one entry for visualization purposes. Features of variable importance (VarImp)≥20 from each model were included in these comparisons. Asterisks highlight top shared features high in DC (IMC_2_Avg, kurtosis) or high in PD (skewness). Blue arrow indicates features shared by I+N organ models (correlation_Avg). The AUC p value shown at the top left corner of each ROC panel in (A) and (B) was computed using function roc.area from R package verification (V.1.42), which implements a two-sided Wilcoxon rank-sum test. ASM, angular s moment; AUC, area under curve; CV, cross-validation; DBV, divided by volume (indicating this is a volume-independent second-order feature); DC, disease control; FPR, false positive rate; IMC, informational measure of correlation; I+N, ipilimumab+nivolumab; LN, lymph node; PD, progressive disease; PD-1, programmed cell death protein-1; ROC, receiver operating characteristic; TPR, true positive rate.

The radiomic features of each model were compared across cohorts or metastatic sites to assess features common or unique to each organ site. The same features at different gray levels were collapsed as one. Across all models, “kurtosis” and “IMC_2 average” are the most shared DC-related features, whereas the most shared PD-related feature is “skewness” (figure 4C, asterisks). For the patient-level models, I+N radiomics only and radiomics+clinical models shared~60% of features, while PD-1 radiomics only and radiomics+clinical models shared~30% of features (figure 4C, upper panels). Looking into models constructed at each organ site, the majority of the features are unique to each organ, with overlapping ones including “kurtosis” associated with DC, among others (figure 4C, middle panels). When comparing models for I+N, “correlation average” was a shared feature across the lung and LN models (figure 4C, blue arrow). Taking into account different organ sites, we sought to build a pan-organ model that includes all metastases in predicting organ response DC/PD, with organ site as a covariate (lung, LN, etc). We confirmed that individual patients’ organ metastases were either all in training or test set to prevent data leaking. The pan-organ models reached an AUC of 0.63 (CV) and 0.75 (test set) in I+N, 0.68 (CV) and 0.79 (test set) in PD-1 (online supplemental figure S9A-D). Volume-only pan-organ models returned AUCs of <0.60 in I+N and PD-1 cohorts (online supplemental figure S9E and F). Online supplemental table S5 describes the full list of feature comparisons between models.

In an attempt to explore the potential mechanisms behind the associations of radiomic features with patient outcomes, we evaluated whether specific patient metadata, such as age, sex, BMI, melanoma subtype, or AJCC stage, show different patterns or associations with the important radiomic features from the models. These analyses revealed that the first-order feature, skewness (FOF009, online supplemental table S2), was consistently elevated in tumors of advanced AJCC stage, comparing M1a/b versus M1c/d (p=1.49E-07 by Wilcoxon rank-sum test; p=4.06E-06 by univariable logistic regression). The significance remains after adjusting for cohort (I+N, PD-1, or BRAF) and response (PD or DC) as covariates in a multivariable logistic regression model (p=2.88E-05). Other features did not show consistent patterns in association with the clinical variables investigated.

Discussion

Using CT and MRI images from patients with melanoma treated with anti-PD-1, I+N, and BRAF therapy, we have described organ-specific patterns of DC/PD and have built clinically informed radiomic models that may identify sites of likely progression. We observed that hepatic metastases experienced a significant reduction in size following ICI in DC but also the greatest increase in tumor volume in PD. Consistent with the literature that response and survival are attenuated in patients with liver metastases,28 29 these data emphasize the liver as a primary driver of patient outcome in the context of ICI. In contrast, no variability in BRAF was observed. Few groups have produced models capable of predicting response at the metastatic lesion level, opting instead to average features from all lesions, or select a single representative lesion to analyze.30 The ability to anticipate organ-specific resistance could potentially have clinical utility. Augmentation of systemic treatment regimens with localized therapy based on organ-specific resistance modeling could reserve treatments for subsequent therapy.

Using radiomic features derived from tumors at each metastatic site, we generated predictive models at both the patient and individual organ levels. With caution based on our sample size, we observed that our models performed well compared with published models, especially when including clinical variables.23 We generally observed non-overlapping radiomic features driving model performance, potentially suggesting unique biology driving the antitumor immunity effect from each treatment, within each organ. Several radiomic features we identified have been associated with outcomes in cancer previously. Kurtosis, a measure describing the “peaked-ness” of the voxel intensity distribution of mass, is higher in homogenous lung nodules relative to heterogenous ones31 and associated with good prognosis in metastatic melanoma.32 Consistent with this, a higher baseline kurtosis distinguishes DC from PD in our study at both the patient and organ levels (LN). The same pattern was observed in IMC_2, a gray-level co-occurrence matrix (GLCM) statistic that quantifies the complexity of texture and has been reported to be associated with clinical benefit in NSCLC33 and gastroenteropancreatic neuroendocrine tumors.34 Other features identified (eg, entropy) describe tumor spatial heterogeneity associated with inferior clinical outcome.35–37 It may be worth noting that radiomic heterogeneity has been described to predict outcome to a similar degree to that of tumor size across cancer types.30 38

We emphasize the need for further validation in larger sample sets. Based on post-treatment “real-world” identification, we cannot control for lines of treatment, number and distribution of metastases, or other clinical factors. Our data were aggregated from patients treated within a large health system and we therefore cannot account for differences that might be due to variation between radiology machines in generating scan images. We also acknowledge that the statistical power of our models is less robust given our sample sizes, though only a few studies have described impressively larger melanoma sample sizes. Given the observation that our model performed better on the test set compared with CV on the train/tune set, this discrepancy may be influenced by the small sample size in the test set. ROC-AUC, a commonly used metric for evaluating model performance, can be sensitive to variations in data distribution, especially when working with a limited sample size. Future studies with larger, independent validation cohorts will be necessary to validate our findings and to ensure the generalizability of the predictive models. Furthermore, while we observed promising performance in organ-specific models, these results should be interpreted with caution given the limitations of both sample size and CV methodology. While LOOCV helps mitigate the risk of overfitting, broader validation will be critical to ensure the generalizability of these results to diverse patient populations and treatment settings. Lastly, we believe the medical interpretation of radiomic features remains early in development, and we do not ascribe specific medical meaning to a specific feature. Rather, we would suggest that these features should be studied in larger cohorts to better understand their clinical relevance. In this study, we used a texture-based approach derived from standardized GLCMs and first-order statistics. Future work could incorporate advanced radiomic analytical techniques, such as deep learning models for automatic high-dimensional feature extraction,39 or multimethod ensemble radiomic pipelines for feature integration,40 to capture complex tumor characteristics enhance both mechanistic insights and prediction accuracy.41

In summary, we describe patterns of treatment response and resistance as well as clinically informed radiomic predictive models that can identify individual sites of treatment-refractory lesions. With the validation of this work, future translational investigation of resistance and potentially clinical trials could be enhanced through the integration of these approaches.