The following method was used to construct TMAs. Hematoxylin and eosin slides were reviewed to select tumor foci for each patient. A TMA instrument (3DHISTECH) was used to obtain cylindrical tissue cores from the selected areas of each donor block. Cores were assembled and embedded in the recipient block. Each core was 0.6 mm in diameter and its surface measured 0.282 mm² (2 or 3 high-power fields). The distance from one core to the other was 0.7 or 0.8 mm. 5 μm thick sections were cut from the recipient block to perform immunohistochemistry ( Figure 1 ).

Firstly, the expression of Notch Intracellular Domain (NICD VAL 1744 clone D3B8, dilution 1:100, Cell Signaling Technology), Jag1 (JAG1 clone D4Y1R, dilution 1:100 Cell Signaling Technology) and E-cadherin (clone 36, Ventana Roche, ready-to-use) was studied and only certain staining patterns were considered positive. NICD, Jag1 and E-cadherin staining were considered positive when they showed cytoplasmic and/or nuclear reactivity, cytoplasmic and/or membrane reactivity and membrane reactivity, respectively. Both the intensity and the percentage of stained cells were evaluated. The intensity was assessed as 0 = negative, 1 = weak, 2 = moderate and 3 = strong.

Secondly, the absolute number of CD4+ (clone SP35) and CD8+ (clone SP57) lymphocytes was counted in three different high-power fields and the mean values obtained were used to determine the CD4/CD8 ratio. A CD4/CD8 ratio of 2.0 was considered normal. The assessment of CD3+ (clone 2GLV6), CD44+ (clone SP37) and CyclinD1+ (clone SP4-R) was given as a percentage of positive cells (Ventana Roche, ready-to-use).

Thirdly, the expression of DLL3 (clone SP347, Ventana Roche, ready-to-use) and DLL4 (clone 4A11F8 dilution 1:100 Biorbyt) was studied. Positivity was defined as ≥25% tumor cells, high expression of DLL3/DLL4 was defined as ≥75% tumor cells. The intensity was assessed as 0 = negative, 1 = weak, 2 = moderate and 3 = strong.

Lastly, mismatch repair proteins (MLH1 - clone M1, PMS2 - clone A16-4, MSH2 - clone G219-1129, MSH6 - clone SP93; Ventana Roche, ready-to-use) were studied. Negative expression of one of them was considered proof of microsatellite instability. In fifteen cases we assessed angiogenesis by counting the microvessels at 20x magnification (Nikon, Eclipse E400).

The immunostained sections for CD31 (Ventana, ready-to-use solution) were examined at low power to select the three areas with the highest vascularity (hotspots).

Two pathologists separately assessed each case without any clinical information.

Patients that underwent abdominal Computed Tomography (CT) at the University Hospital of Parma for CRC staging were included in the study. CT scans were performed with different CT scanners and imaging protocols; images were retrieved from Picture Archive and Communication System (PACS) and were subsequently imported into a dedicated software (3D Slicer) for tumor segmentation.

One radiologist (ML) evaluated all CT scans visually and identified the target lesion on portal venous phase. The reader was instructed to draw manually multiple regions of interest (ROI) at different levels by tracing the boundaries of the lesions: subsequently, a dedicated tool (SlicerRadiomics) software interpolated the ROIs to obtain the volume of interest (VOI) which allowed the extraction of 852 radiomic features (RF). The VOI was manually modified by the reader in case of inaccurate segmentation. Image preprocessing based on wavelet decomposition was performed by SlicerRadiomic before feature calculation to generate independent radiomic predictors.

The radiomic dataset included shape, first-order, Gray-Level- Co-occurrence-Matrix (GLCM), Gray-Level-Run-Length-Matrix (GLRLM), Gray-Level-Size-Zone-Matrix (GLSZM), Neighboring-Gray-Tone-Difference-Matrix (NGTDM), Gray-Level-Dependence-Matrix (GLDM).

The chi-square test and Fisher’s exact test were used to perform univariate comparisons between categorical variables.

The Kaplan-Meier method was used to estimate the mean and median time for progression free survival (PFS) followed by a Cox regression analysis to evaluate the relationship between survival and covariates in a multivariable framework. The model was evaluated by making use of model diagnostics. This included checking for the overall goodness of fit, model adherence to key assumptions, influential observations and nonlinearity. The variables considered in the Cox regression were KRAS, type of chemotherapy protocol, site of primary tumor, NICD, CD44, Jag1, CD3, DLL4 expression; only NICD expression resulted statistically significant and was maintained in the final model. The regression coefficients were reported as hazard ratios (HRs). The 95% confidence intervals (CIs) were also estimated from the analysis.

The commercial package IBM-SPSS v.28 and the open-source statistical system Jamovi version 2.3.0, which is based on the widely used open-source system R, were used to perform survival analysis. A p-value less than 0.05 was considered statistically significant (p < 0.05).

Classification models were developed to predict time to disease progression. With this aim in mind, PFS at 9 months was used to stratify patients in two groups, namely short and long-term survivals. We considered PFS at 9 months as a target variable because a comprehensive meta-analysis has recently showed that PFS ranges between 7 and 10.8 months for CRC patients treated with bevacizumab (18). Therefore, we acquired the central value of that interval from the meta-analysis to further stratify the prognosis of our patients according to the integrated profile. Three models were developed: radiomic (R), clinical/Notch signaling (C/N) and the comprehensive integrated model (I). In the R model, we removed redundant highly correlated features by calculating their Spearman Rho correlation coefficient: RFs with a coefficient greater than 0.99 were excluded from the successive analyses. Subsequently, feature standardization by z-score was applied. In the C/N model, the same variables considered in the Cox regression were added as predictors ( Table 1 ). In both R and C/N models, a L2 penalized logistic regression algorithm was implemented for features selection and validation. Most predictive features were selected by means of a wrapper approach, i.e. the sequential forward feature selection algorithm with 20 Monte-Carlo cross-validation (MCCV) splits. We chose the area under the receiver operating characteristic (ROC) curve (ROC-AUC) as performance metric. We iterated the selection process 50 times, using 50 different random states and, subsequently, we selected the features that had higher frequency of occurrence for both R and C/N models. Through 5000 MCCV splits (train:0.7, test:0.3), different numbers of clinical/genomic and radiomic features were selected and used respectively for R and C/N model training and validation. Likelihood ratio test was applied to verify if the addition of another feature significantly improved the model performance. The selected C/N and R features were used together to build the I model.

For all models, the ROC-AUC and accuracy scores for each MCCV were calculated and averaged over these iterations. Mean ROC curve and mean learning curve were also plotted. Recall and precision metrics were also calculated

Survival analysis was performed in radiomic dataset. Kaplan–Meier survival curves PFS for two risk groups were calculated and then compared using log-rank test. The risk groups were assessed by using continuous RFs, previously selected by the machine-learning model. Risk groups based on RFs were developed using ROC analysis to determine the cutoff value of each RF for optimal stratification into two classes: Youden index was chosen as optimal threshold. Subsequently, we combined the selected features in a single variable and we performed Kaplan–Meier analysis again. Finally, we calculated the probability to predict longer-term class in each risk group of combined features. Probabilities derived from the R model were averaged over MCCV splits.

Machine-learning model, analysis and plots were performed by means of Python v. 3.8.5; scikit-learn and MLextend machine learning libraries were used for features selection and model development.

A total of 111 patients have been included in the analysis. The cohort is shown in Table 1 . A positive association was found in univariate analysis between NICD and Jag1 expression (p < 0.001; Table 2 ). No significant association was found for the other analyzed markers and KRAS mutation (data not shown). All main clinical characteristics were comparable among the subgroups of patients (data not shown). Specifically, no significant associations of NICD and Jag1 immunostaining scores with age, baseline CEA levels, number of metastatic sites and subsequent chemotherapy were observed.

Compared with patients who had NICD and/or Jag1 low tumors, patients whose pre-treatment tumors expressed high NICD and Jag1 levels showed poor RECIST 1.1 categories with higher rates of stable disease (SD) or progressive disease (PD) as best response, and lower frequencies of complete response (CR) or partial response (PR); p = 0.002 ( Table 3 ). Associations between NICD and Jag1 and therapy response were further evaluated using PFS and Kaplan-Meier and Cox proportional hazard modeling. Median PFS was significantly shorter in patients whose tumors expressed high NICD and Jag1 (6.43 months vs 11.53 months for negative cases; p = 0.001, Figure 2 ). Cox regression following univariate analysis confirmed NICD as the only independent predictor for PFS (HR = 1.820 [1.165 – 2.844]; p = 0.009).

Quite surprisingly, 5 patients with high NICD tumors showed long PFS. Each case was evaluated for the following features: inflammation, staging, grading and microvascular density. The last one was the only noteworthy characteristic. For this reason, we assessed the microvascular density according to Chalkey’s methods: microvessels were counted manually for each hotspot at 20x magnification (high power field) and expressed as MVD score. This assessment was carried out for 15 patients based on response to therapy: 5 were non responder (NICD/Jag1 positive), 10 were responder (5 NICD/Jag1 positive and 5 NICD/Jag1 negative). Those with an MVD score ≥5 (CD31-high, NICD/Jag1 positive) were associated with significantly poorer survival. Low CD31 was seen in all 10 responder patients (both 5 NICD/Jag1 positive and NICD/Jag1 negative) and associated with a better prognosis.

The retrospectively collected 111 CRC cases were decreased due to inclusion criteria that comprised the availability of (i) CT data and (ii) PFS information. Thus, the ensuing results based on the multiomic approach refer to a restricted population of 76 subjects. Regarding feature preprocessing, the Spearman correlation matrix for RFs is reported in Figure S1 . Redundant features were removed, thereby reducing the number of RFs by about 33.6%.

In the R model, the best performance in differentiating short and long survival was obtained by selecting two RFs: Strength (NGTDM) and Skewness (first order) with a ROC-AUC of 0.709 and an accuracy of 0.671. Likelihood ratio showed that the performance of the model would not significantly improve by adding further RFs ( Table 4 ).

In the C/N Model, the most predictive features were Jag1 and NICD, whereas the addition of a third feature was not significantly relevant: the ROC-AUC was higher as compared to the R model alone, with a ROC-AUC of 0.743, while accuracy was slightly lower (0.649) ( Table 4 ).

The I Model included the previously selected R and C/N features, that had a negative effect on the likelihood of predicting long survivors, as shown by odds ratios ( Table S1 ). The I model yielded the highest ROC-AUC (0.823) and accuracy (0.751) values. The mean ROC curves are displayed in Figure 3 .

Learning curves show the ROC-AUC score ( Figure S2A ) and accuracy ( Figure S2B ) as a function of the number of training samples. We plotted performance scores obtained by predictions for both training (blue line) and validation (green line) datasets and averaged over all iterations of the MCCV. For each model, we also calculated and averaged over 5000 MCCV splits the recall metric representing the true predictions of longer survival class: 0.794 [95% CI : 0.790, 0.797] for R model, 0.660 [95% CI: 0.654, 0.666] for C model and 0.767 [95% CI : 0.763, 0.770] for I model. In addition, we calculated precision metrics (i.e. positive predictive value), representing the fraction of true positive cases among the total positive predicted instances: 0.642 [95% CI : 0.641, 0.645] for R model, 0.673 [95% CI: 0.669, 0.676] for C model and 0.751 [95% CI : 0.749, 0.754] for I model.

Kaplan Meier curves of PFS ( Figure S3 ) showed significantly different risk strata for Strength, whereas none for skewness. The combined RF (strength-skewness) created 3 risk groups which significantly stratified in PFS curve ( Figure 4 ). Probabilities of longer-term class prediction are listed in Table 5 .