Performance and Dimensionality of Pretreatment MRI Radiomics in Rectal Carcinoma Chemoradiotherapy Prediction

1. Introduction

Rectal carcinoma was the third most common cancer and the second most fatal cancer with an estimated 1.4 million new cases worldwide in 2020 [1]. Due to increasing incidence in younger people [2, 3], the American Cancer Society now recommends screening starting at age 45 [4].

Rectal carcinoma prognosis varies by stage, with 5-year survival rates of over 90% for early-stage and ~20% for late-stage disease [5]. The standard treatment for locally advanced rectal carcinoma is neoadjuvant chemoradiotherapy (nCRT) followed by total mesorectal excision, with or without adjuvant chemotherapy [6]. This treatment is highly effective, with low rates of local recurrence [6].

The clinical assessment of neoadjuvant treatment response for rectal carcinoma occurs before surgery, using MRI scans. Clinical complete response (cCR) is determined as the absence of residual disease in the rectum after neoadjuvant treatment, as confirmed by digital rectal examination, endoscopic evaluation (rectosigmoidoscopy) and control MRI examination. Another method of assessing treatment response is tumor regression grading (TRG), which is based on the histopathological examination of tumor tissue postoperative specimen. TRG is graded on a scale of 1 to 5, with TRG1 representing the best response and TRG5 indicating the poorest response, characterized by no observable tumor regression. TRG1 corresponds to pathohistological complete regression (pCR), which is observed in 10-30% of cases [7] and defined as the absence of viable tumor cells after total mesorectal excision (ypT0N0). It has been shown that pCR is associated with improved outcomes, regardless of the initial clinical T and N stages of the disease [8].

TRG is an important prognostic factor for local recurrence, disease-free survival, and overall survival (OS) of rectal carcinoma patients. Patients with TRG 1-2 (good responders) have significantly better outcomes than poor responders with TRG 3-5 [9].

More reliable predictors of nCRT response are needed to enable precision medicine in routine clinical practice, thereby maximizing the effectiveness of existing treatments. Personalized treatment for chemoresistant patients may include second-line chemotherapy, participation in experimental trials with more targeted therapies, more intense radiotherapy, or combination therapy with immunotherapy. Conversely, to improve the quality of life for patients responding well to neoadjuvant treatment, adjustments such as less invasive surgery or a non-operative approach ("watch and wait") may be considered. However, such precision in clinical decision-making is hindered by the limited reliability of current clinical, pathological, and molecular predictors of chemoradioresistance in rectal carcinoma, such as tumor stage, tumor regression grade, tumor markers (carcinoembryonic antigen), circulating tumor DNA, DNA methylation level, and cancer-related inflammatory markers [10, 11].

In our previous work, we compared the proteomic profiles of LARC patients who responded well and poorly to nCRT [12]. We also identified Methylenetetrahydrofolate reductase (MTHFR) 667C and 1298A alleles as low-penetrant risk factors for rectal cancer [13]. We also investigated the predictive value of clinicopathological features in rectal carcinoma [14]. In this study, we enhance our predictive research by integrating radiomics analysis.

Radiomics analysis of MRI scans can provide valuable information for predicting nCRT response in rectal carcinoma, complementary to traditional clinical and molecular methods. Radiomics extracts quantitative data from medical images that are not visible to the naked eye. Most predictive imaging models for rectal carcinoma previously used PET scans [15, 16], while several used MRI. Some of these studies have used post-nCRT MRI to assess treatment response rather than predict it [17-19]. Other studies have used pre-treatment MRI to predict nCRT response, predominantly relying on features computed from unfiltered images [20], with occasional inclusion of wavelet [21] and Laplacian of Gaussian (LOG) filters [22-24]. Image filters are used in MRI radiomics to better emphasize texture or structural information. However, due to the limited use of filters in radiomics analysis of rectal carcinoma, the predictive potential of several commonly used filters, such as logarithmic, exponential, square, square-root, exponential and LBP, has remained unexplored. Additionally, the limited use of image filters in previous predictive radiomic MRI studies has prevented direct evaluation of how radiomics dimensionality affects predictive performance. Notably, a recent study investigating the effect of image filtering on radiomic features recommended conducting radiomic analysis using all available filters [25].

Consequently, to address the clinical need for improved prediction of nCRT response in rectal carcinoma, this study systematically evaluated the predictive potential of all available image filters and dimensionality in radiomic MRI analysis.

2. Results

The predictive model consisted of two radiomics scores: a moderately dimensional score derived from 72 clinicopathological (CP) and 107 radiomic features, and a very highly dimensional score derived from 72 CP and 1862 radiomic features. We also computed a CP score to serve as the benchmark for state-of-the-art performance in predicting nCRT response. The nCRT continuous outcome was defined in the order: clinical complete response (cCR), tumor regression grade 1 (TRG1), TRG2, TRG3, TRG4, and TRG5. For receiver operating characteristic (ROC) analysis, we also categorized this outcome into responders (cCR, TRG1, TRG2) and non-responders (TRG 3-5) [21].

Table 1 provides a comprehensive overview of patient characteristics, disease specifics, treatment details, and outcomes. Notably, direct tumor spread was the primary form of involvement when the mesorectal fascia was affected. Pathological regional lymph nodes were present in 97.3% of cases, and one-third of patients displayed extramural vascular invasion (EMVI). In our cohort, 46.7% of patients were categorized as treatment responders.

We observed that LASSO sometimes selected features that were weakly associated with the outcome. Therefore, we removed all features that did not significantly correlate with the outcome before applying LASSO feature selection (not shown). This pre-selection improved the predictive performance of the scores selected by LASSO. In contrast, removing highly correlated features before LASSO did not affect the final results. The workflow of our analysis is presented in Figure 1.

Table 2 presents features achieving the strongest predictive performance, from each of the three categories: CP features, radiomic features computed from original MR images and radiomic features computed from both original and filtered MR images. Notably, out of the 72 available CP features, only the five presented here significantly associated with the outcome. Univariate predictive analysis of the features in Table 2 using the Pearson correlation test and univariate linear regression showed that individual features from the CP and two radiomic feature groups had comparable associations with the outcome.

We calculated a basic set of 107 radiomics features from the original images, including 14 size/shape, 18 intensity, and 75 second-order texture features. To enhance the depth of our analysis, we applied 21 image transformations using various filters, including LoG with five settings, wavelet with eight subbands, LBP 3D with two settings, and square, square-root, logarithmic, exponential, and gradient filters. This increased the total number of features to 2141. However, we removed 279 features computed on images produced by wavelet decomposition into LHH, HLH, and HHH subbands due to low gray level intensity ranges. This reduced the total number of features to 1862.

To address the high dimensionality of the data and improve the predictive performance of the final models, we pre-selected radiomics features by removing those that did not show significant Pearson correlation with the outcome. This reduced the number of radiomics features from 1862 to 365. We did not pre-select CP features because only five of them were significantly associated with the outcome.

After pre-selection, we used linear LASSO regression followed by multivariate stepwise backward linear regression to select features. From the initial two radiomics sets consisting of 72 CP + 107 radiomics features and 72 CP + 1862 radiomics features, LASSO selected nine and eleven features, respectively (not shown). Multivariable linear regression further refined the models to five features each (Table 3).

To assess the benefits of comprehensive and very high-dimensional tumor characterization, we compared the predictive performance of the two models described above: one derived from a smaller set of 107 features and one derived from a much larger set of 1862 radiomics features (Table 3). Table 3 shows that the very high-dimensional model had only a marginal improvement in predictive capability, with an R-squared value of 0.47, compared to the R-squared value of 0.45 for the model with 107 features. This finding was consistent with the ROC analysis presented in Figure 2.

ROC analysis requires a binary outcome, so we categorized patients into responders (cCR, TRG1, and TRG2) and non-responders (TRG3 and TRG4) (Figure 2). cCR means that no cancer was detected in the primary tumor by imaging and physical examination. TRG ranges from TRG1 and TRG2 (complete response and near complete response, respectively) to TRG3 and TRG4 for moderate and poor response. Figure 1 shows that the score with highest dimensionality computed using 75 CP and all 1862 radiomics features had the best predictive performance.

Representative MR images of nCRT responders and non-responders are shown in Figure 3. The computed values of features (Figure 3) illustrate the classification of images based on nCRT response. In the case of these four patients, original NGTDM complexity and 72 CP + 1862 radiomics score were consistently higher in non-responder patients, while square-root 90th percentile and exponential gray level nonuniformity were not consistent. This association of increasing feature values with increasing TRG (reflecting lower response to nCRT) was in line with the data presented in Table 1, because these features yielded AUC values above 0.5. The texture features of tumors are subvisual, making it difficult to visually distinguish between images of tumors from responder and non-responder patients (Figure 3).

3. Discussion

In this study, we developed two novel MRI-based radiomics models to predict the response of LARC nCRT. We also examined the relationship between the depth of radiomics analysis and predictive performance.

This study highlights the complementary predictive power of clinicopathological and radiomics features. The association of mucinous tumor differentiation with poor treatment response was consistent with the findings of Simha et al. [26]. Similarly, MRI-derived EMVI was associated with poor responses in this study, as in previous findings [27]. Such EMVI-associated chemoradioresistance may be explained by tumor hypoxia, prompting the investigation of dose escalation with adaptive MRI-guided radiotherapy for this group of patients [28]. The identification of basophil counts as a predictive parameter in LARC was also consistent with our previous study [14] and findings in advanced gastric cancer [29].

In line with most previous predictive studies in rectal carcinoma that used radiomics analysis of MRI data, we investigated the T2W sequence. This sequence is considered the most suitable for predictive purposes because it can capture predictive clues well, is less susceptible to artifacts, and has good reproducibility. [30].

With inclusion of all available image filters, this study reports that the features selected for predictive scores derived from the original images and those generated using previously unutilized exponential and square root filters. This finding is consistent with a previous study that investigated the impact of preprocessing filters on the predictive performance of radiomics analyses across seven radiomics datasets [31], whereby square root, exponential, and wavelet filters achieved the best predictive performance [31]. The exponential filter can enhance image boundaries and contrast, while the square root filter highlights subtle textural patterns that might otherwise be difficult to detect and has noise reduction properties.

Omics methodologies like radiomics exploit the advantages of comprehensive approaches involving thousands of features. However, using many features risks false discoveries, where statistical significance may arise by chance. One goal of this study was to assess the predictive benefits of the deepest radiomics approach. To achieve this, we compared predictive models based on large and small numbers of radiomics features. The model trained on 1862 features had a moderate predictive advantage over the model trained on 107 features (AUC 0.90 vs. 0.86). This suggests that using a very large number of features can capture additional predictively relevant information, even if a moderately dimensional analysis with only 107 features already delivers very good predictive performance. Although the improvement from AUC 0.86 to 0.90 may seem incremental, further improvements become increasingly difficult as predictive performance nears the ideal AUC of 1.0. Therefore, even small improvements in predictive performance by 0.04 AUC can have substantial clinical value, reducing the number of incorrectly classified patients by 5%.

In this study, we developed radiomics signatures with 3-5 features after dimensionality reduction, while previous studies used up to 30 features [18, 19, 32] A common guideline is to use one feature for every ten patients [33, 34]. As expected, our predictive scores, which integrated both clinicopathological and radiomics features, outperformed the benchmark predictive score based only on clinicopathological features. This finding aligns with the previous report that radiomics features add value when combined with qualitative MRI features [20] while most studies did not provide such comparison. However, some studies have reported clinicopathological scores outperforming models that incorporate radiomics features [35].

While chemoradiotherapy responsiveness is a surrogate endpoint that provides valuable insights into treatment efficacy and short-term outcomes, it is important to recognize its limitations because it serves as a measure for the true clinical outcomes of interest, like disease recurrence or long-term survival [36]. Not all patients who respond to nCRT will experience long-term disease control and conversely, not all non-responders will have poor outcomes. Furthermore, although our patient cohort was highly homogeneous and largely exceeded the required sample size, its overall size remains a limitation. Although we performed thorough internal validation using cross-validation and bootstrapping, further studies in larger external patient groups are needed to confirm the clinical validity of the reported associations between the computed scores and chemoradiotherapy outcomes. Additionally, while the computational analysis technique is entirely objective, there is still some residual subjectivity involved in tumor VOI selection for radiomics analysis. The retrospective design of the predictive model is another limitation. Also, the feature selection using LASSO is susceptible to multicollinearity and lacks the statistical test for significance. In our study, we found that LASSO effectively handled multicollinearity, as removing inter-correlated features before the LASSO selection did not impact the results. However, the removal of features with no significant correlation with the outcome yielded favorable results, indicating that through feature pre-processing, we successfully addressed a major limitation in the feature selection process. Other studies have employed similar multi-step approaches to feature selection, involving initial steps like the Wilcoxon rank-sum test, Spearman correlation analysis, followed by LASSO, and multivariate logistic regression analysis [32, 37].

4. Materials and Methods

We utilized a retrospective cohort of patients with locally advanced rectal carcinoma (LARC) who have undergone diagnostic MRI and neoadjuvant chemoradiotherapy. Pre-treatment magnetic resonance imaging (MRI) scans were subjected to radiomics feature extraction, encompassing morphological and textural aspects of the tumor. LASSO was then employed to train and validate the predictive model, employing statistical and cross-validation techniques to select the predictively most valuable features, thus ensuring robustness and generalizability.

Ethics Approval Statement

The study received approval from the Ethics Committee of the Institute for Oncology and Radiology of Serbia (Approval No. 2211-01 from 11.06.2020) and Ethics Committee of the Faculty of Medicine, University of Belgrade (Approval No. 1322/XII-17 from 03.12.2020). It adheres to The Code of Ethics of the World Medical Association (Declaration of Helsinki), as published in the British Medical Journal (July 18, 1964) and its 7th revised edition in 2013. All patients signed an informed consent.

Patients

A total of 75 patients diagnosed with LARC were enrolled in this study. The patient cohort was selected from the Institute for Oncology and Radiology of Serbia, spanning the period between June 2020 and January 2022. Inclusion criteria required patients to have histopatologically confirmed adenocarcinoma of the rectum, with the tumor located within up to 12 cm from the anal verge, as determined by rigid proctoscopy. LARC was defined as encompassing T3-T4N0 stages or any T stage with positive lymph nodes (N+). Pretreatment assessment included abdominal and pelvic MRI scans, as well as computed tomography (CT) scans or chest X-rays. All patients underwent long-course nCRT. Radiotherapy (RT) was administered using the volumetric modulated arc therapy-simultaneous integrated boost technique (VMAT-SIB). The prescribed dose to the mesorectum and pelvic lymph nodes was 45 Gy, administered in daily fractions of 1.8 Gy. Additionally, a simultaneous integrated boost (SIB) was administered to the macroscopic disease region with a 2 cm margin, totaling 54 Gy, delivered in daily fractions of 2.16 Gy. Concomitant chemotherapy was initiated on the first day of radiotherapy and continued during the first and fifth weeks of the treatment regimen. The chemotherapy regimen consisted of 5-fluorouracil (5-FU) at a dose of 350 mg/m2 on the first day of the first and fifth weeks of radiotherapy, along with leucovorin (25 mg/m2 daily) administered during the five consecutive days of the first and fifth weeks of radiotherapy.

The evaluation of tumor response was conducted eight weeks following the completion of nCRT and included pelvic MRI scans, rigid proctoscopy and digital rectal examination. For patients achieving complete clinical response (cCR) and initially distant tumor location, immediate radical surgery was not recommended. Instead, they were enrolled in a stringent follow-up program ("watch and wait" approach). Patients with cCR who were candidates for sphincter preservation surgery underwent surgical resection within a window of eight to twelve weeks following the completion of nCRT. For patients exhibiting a partial response (PR), surgery was conducted approximately twelve to fifteen weeks after the completion of nCRT.

Patients' responses to treatment were categorized using the Mandard classification system based on the pathohistological TRG observed in postoperative specimens. Responders included patients with cCR who did not require surgery, as well as those with TRG 1 and TRG 2 postoperative categories. Non-responders encompassed patients classified as TRG 3-5.

MRI

Initial magnetic resonance imaging (MRI) data of the pelvic region were available for 71 out of the total 75 patients included in this study. All MRI examinations featured 1.5-Tesla 3D T2-weighted (T2W) contrast sequences, which were utilized for precise tumor delineation. It's noteworthy that the MRI procedures were performed at a different institution, presenting both a limitation and an advantage to this study. While this difference in the origin of MRI data is a limitation, it simultaneously offers an advantage by increasing the generalizability of the study results, as they are not contingent on specific MRI acquisition protocols.

Postprocessing

The initial step in our analysis involved the import of MRI Digital Imaging and Communications in Medicine (DICOM) files into 3D Slicer (version 5.3.0) and the subsequent generation of Multiple Resolution Bitmap (MRB files), streamlining the data for further examination [38]. For all 71 patients under consideration, precise segmentation of rectal carcinoma tumor volumes, which refers to the process of delineating regions of interest (ROI), was carried out. Sequences exhibiting any artifacts were excluded from the analysis to ensure the highest possible image quality.

Tumor volume delineation was a collaborative effort involving both a radiation oncologist and reference to the initial radiologist reports, which provided critical guidance in ensuring accuracy and consistency. The radiation oncologist (M.M.) with six years of expertise in oncologic imaging manually segmented the volume of interest (VOI) on all 71 T2-weighted images using open-source 3D Slicer software. The segmentation of the tumor was performed on each image slice where the tumor was visible. Figure 1 shows the MRI sagittal plane before and after segmentation.

Sample size calculation

The prospective sample size calculation was based on a pilot experiment involving 36 patients. It required 40 patients with 18 positive cases for alpha = 0.05, beta = 0.20, and area under the ROC curve AUC=0.75 (Medcalc 14.8.1; MedCalc Software Ltd., Ostend, Belgium). The actually obtained AUCs for the two calculated scores including clinicopathological and radiomics features were 0.86 and 0.90, with a final sample size of 71 patients.

MRI normalization

Z-score normalization was applied to MR images as instructed by the Parameter file listed below. Prior to radiomics feature extraction in 3D, all images and segmentations were resampled and interpolated to obtain isotropic voxels of 1 mm³.

Feature extraction

The computational analysis of rectal MRI images was conducted using the open-source Pyradiomics plugin, integrated into 3D Slicer (version 5.3.0), named "Radiomics"[39]. PyRadiomics is compliant with Imaging Biomarker Standardization Initiative (IBSI) standards. The Pyradiomics software was customized through the parameter file provided below, instructing the generation of all available image transformations and feature types, resulting in computation of a total of 2141 features per image.

The width of discretisation bins for feature extraction was adjusted separately for each image filter in a pilot analysis of 40 MRI, targeting a bincount range between 30-130. Feature extraction was performed on the original images (107 features) as well as after applying eight types of built-in pre-processing filters: wavelet (5x93=465 features), square (93 features), square root (93 features), logarithm (93 features), gradient (93 features), exponential (93 features), Laplacian of Gaussian (LoG, 5x93=465 features), Local Binary Pattern (LBP 2D, 89 features; LBP 3D 3x92=276). Wavelet filtering, yields 8 subbands, all possible combinations of applying either a High or a Low pass filter in each of the three dimensions.

Among the basic set of 107 original features, 18 were first-order intensity and 14 shape, while the remaining 75 were second-order texture features belonging to six distinct classes: gray-level co-occurrence matrix (GLCM, 24 features), gray-level run length matrix (GLRLM, 16 features), gray-level size zone matrix (GLSZM, 16 features), gray-level dependence matrix (GLDM, 14 features), and neighboring gray-tone difference matrix (NGTDM, 5 features). We adhered to the common assumption of dimensionality: moderate: 10-100 dimensions, high dimensionality: 100-1000 dimensions and very high dimensionality: 1000+ dimensions. Detailed descriptions of the extracted radiomics features can be found in the PyRadiomics documentation available at https://pyradiomics.readthedocs.io/en/latest/features.

Normalization of the calculated feature values

The computed variable values ranged between 6.2E10 end -3.810E6. By employing z-score normalization, we were able to bring values of all variables to a similar scale, with values ranging between -6.60 and 8.31, thus enabling comparison among features.

Evaluation of predictive performance

Evaluation of the predictive performance for the demographic, MRI, clinicopathological and the radiomics features was approached by the univariate linear regression and receiver operating characteristic (ROC) analysis, both executed using IBM SPSS software package version 28 by IBM Corporation in Armonk, NY, United States. Our study focused on the prediction of the response to chemoradiotherapy in rectal carcinoma as the primary endpoint. In the case of univariate linear regression, we worked with continuous values for independent variables, while the dependent endpoint was also expressed in continuous numerical values, ranging from 0 to 4. Specifically, we assigned values as follows: cCR = 0, TRG1 = 1, TRG2 = 2, TRG3 = 3 and TRG4 = 4. Linear regression assesses whether predictor independent variables explain the dependent variable (endpoint or outcome).

ROC analysis requires categorized dependent endpoint variable. The area under the rate of change curve (AUC) was used to evaluate the discriminatory efficiency of the prognostic features in relation to a binary endpoint. Discrimination refers to the ability of the prognostic features to classify patients based on their actual metastasis occurrence. AUC values range from 0.0 to 1.0, where 0.5 indicates random chance discrimination, while values of 0.0 or 1.0 indicate perfect discrimination. Furthermore, it is of note that the AUC not only quantifies the level of discrimination but also offers directional information. For instance, AUC values of 0.2 and 0.8 suggest the same level of discriminatory efficiency but in opposite directions. We set a significance threshold at P ≤ 0.05 to determine the statistical significance of our findings.

Model selection

Predictive models were constructed using the features selected by LASSO regression (Stata/MP 16, StataCorp, College Station TX, United States) followed by stepwise multivariate linear regression (IBM SPSS v28). LASSO was performed by the inclusion of radiomics features together with all 72 clinicopathological features.

Although LASSO handles multicollinearity and removes features unrelated to the outcome, we performed feature preselection, prior to LASSO and multivariable regression. Radiomics features which did not significantly correlate with outcome by Pearson coefficient were discarded. Clinicopathological features were exempt from the above selection because only five features were significantly associated with the outcome. Dimensionality was further reduced by use of LASSO regression (Stata/MP 17, StataCorp, College Station, TX, United States), followed by backwards stepwise multivariable binary regression (IBM SPSS v28). LASSO regression is a machine learning regularization method that selects covariates and estimates their coefficients using a tuning parameter λ in cross-validation. As λ increases, LASSO eliminates less important variables by shrinking their coefficients to zero, providing robustness against outliers and improving generalization performance. This approach allows for the selection of features with the most robust prognostic performance. Multivariable stepwise backwards linear regression proceeds to remove the least significant variables, one at a time. The removal of variables continues until the p-value of the next feature to be removed falls below 0.05, satisfying the stopping rule.

Features remaining after the applied selection methods were used for computation of the Radiomics Scores, based on the selected features and the coefficients provided by the regression methods, using the formula: variable1*coefficient1 + variable2*coefficient2 + variable3*coefficient3.

Validation

To mitigate the over-optimistic bias in the univariate and multivariable logistic regression analysis, the bootstrap technique was employed, generating 5,000 random resamples of the data, allowing for more robust estimation [40]. The original confidence intervals (95% CIs) and p-values were adjusted using this resampling approach, ensuring more reliable results (IBM SPSS v28).

In addition, split-sample cross-validation was utilized as a validation method to determine the optimal penalty coefficient λ in the LASSO regression analysis conducted in Stata/MP 17. This technique involved dividing the data into two subsets: a training set and a validation set. The training set was used to fit the LASSO regression model and select the most informative features, while the validation set was used to evaluate the performance of the model.

5. Conclusions

We developed very highly dimensional and moderately dimensional predictive models for neoadjuvant chemoradiotherapy (nCRT) response in rectal carcinoma, using clinical parameters and pretreatment MRI radiomics features. These models outperformed the benchmark model relying on clinicopathological parameters alone. We also showed that using all available image filters and high dimensionality improves predictive performance. The combined clinicopathological and radiomics scores obtained in this study provide a foundation for future broader multimodal predictive approaches that also incorporate genomics [14] and proteomics [12] data. Improved prediction of nCRT response might be clinically significant because it may enable personalized treatment decisions such as delaying surgery or using less aggressive post-operative therapies to minimize toxicity in predicted nCRT responders. For predicted non-responders, more intensive therapies with shorter follow-up periods may reduce recurrence risk.