Predictors of Outcome after (Chemo)Radiotherapy for Node-Positive Oropharyngeal Cancer: The Role of Functional MRI

Simple Summary We have conducted a prospective study on patients with locally advanced oropharyngeal cancer who are candidates for concomitant radio-chemotherapy; we considered their anamnestic findings, tumor characteristics and evaluated the role of innovative radiological features, particularly the magnetic resonance imaging (MRI) biomarkers. Our aim was to identify those elements correlated with worse tumor control. Diffusion-weighted (DWI) imaging and dynamic contrast–enhanced (DCE) can help identify hypoxic regions in head and neck cancer which are known to be more resistant to the effects of radiation. A better understanding of these factors may help us improve our knowledge on tumor behavior and thus provide a more tailored treatment in patients that respond poorly. Abstract The prognosis of a subset of patients with locally advanced oropharyngeal cancer (LA-OPC) is still poor despite improvements in patient selection and treatment. Identifying specific patient- and tumor-related factors can help to select those patients who need intensified treatment. We aimed to assess the role of historical risk factors and novel magnetic resonance imaging (MRI) biomarkers in predicting outcomes in these patients. Patients diagnosed with LA-OPC were studied with diffusion-weighted imaging (DWI) and dynamic-contrast enhanced MRI at baseline and at the 10th radiotherapy (RT) fraction. Clinical information was collected as well. The endpoint of the study was the development of disease progression, locally or distantly. Of the 97 patients enrolled, 68 were eligible for analysis. Disease progression was recorded in 21 patients (11 had loco-regional progression, 10 developed distant metastases). We found a correlation between N diameter and disease control (p = 0.02); features such as p16 status and extranodal extension only showed a trend towards statistical significance. Among perfusion MRI features, higher median values of Kep both in primary tumor (T, p = 0.016) and lymph node (N, p = 0.003) and lower median values of ve (p = 0.018 in T, p = 0.004 in N) correlated with better disease control. Kep P90 and N diameter were identified by MRMR algorithm as the best predictors of outcome. In conclusion, the association of non-invasive MRI biomarkers and patients and tumor characteristics may help in predicting disease behavior and patient outcomes in order to ensure a more customized treatment.


Introduction
The initial approach in patients with locally advanced head and neck squamous cell carcinoma (HNSCC) is a combination of radiotherapy and concomitant chemotherapy (RT-CT) in order to avoid surgery and to preserve organ functionality. However, the outcome of these patients is still poor, with 40-50% experiencing disease recurrence [1,2]. Patientand tumor-related predictors of response to RT-CT would be helpful for a more efficient selection of candidates to non-surgical approach and to avoid unnecessary toxicity in nonresponders. In the last decades, some major changes have occurred in the epidemiology and management of these patients, especially those with oropharyngeal cancer (OPC), and new prognostic factors have been identified. Currently, HPV-positive OPC is considered a distinct entity that is usually associated with a younger age of onset, distinct clinical features, limited tobacco exposure and a more favorable oncological outcome than HPVnegative OPC [3,4]. Histologically, it is characterized by basaloid, lymphoepithelial and poorly differentiated features [5], and radiologically by cystic-appearing lymph nodes on both magnetic resonance imaging (MRI) and computed tomography (CT).
However, despite the better prognosis, a significant proportion of these patients continues to experience distant disease progression [6,7]. Historical risk factors such as significant smoking history and extranodal extension (ENE) in lymph nodes maintain a significant predictive role of an increased risk of both locoregional (LRF) and distant failure [8][9][10][11].
Response to treatment is also influenced by biological (rather than anatomical) features such as intrinsic tumor radiosensitivity and hypoxia [12,13]. Advances in MR imaging techniques currently allow the estimation of tissue cell density and the localization of hypoxic regions within HNSCC using novel functional biomarkers: diffusion-weighted (DWI) imaging quantifies the water molecule mobility in tissues which is strictly related to cell architecture, while dynamic contrast-enhanced (DCE) MRI allows for the derivation of semi-quantitative hemodynamic maps by the estimation of the passage of blood through vessels [14].
The aim of this study is to prospectively assess by a machine learning (ML) approach the role of both traditional and emerging/novel factors in predicting outcomes in patients with OPC who underwent (chemo)radiotherapy.

Patient Population and Treatment Characteristics
We evaluated 97 consecutive patients treated within a non-randomized, prospective, single-institution trial funded by the Italian Association for Cancer Research (AIRC, project No.17028), and specific informed consent was obtained before enrollment.
Inclusion criteria to be met were: (a) age older than 18 years; (b) histologicallyconfirmed squamous cell carcinoma of the oropharynx; (c) locally advanced tumor stages III and IV according to 8th edition of the AJCC Cancer Staging Manual [15]; (d) definitive treatment with concomitant RT-CT. Exclusion criteria considered were: any contraindication for MR examination due to previous allergic reaction to intravenous contrast material administration or renal disease; patients previously treated with surgery, chemotherapy or radiotherapy to the head and the neck. Disease staging required MRI examination and CT of the neck and chest or PET-CT for distant staging. Baseline characteristics, including smoking pack-years and alcohol consumption, were collected as well. Patients received intensity modulated radiation therapy (IMRT) with a simultaneous integrated boost technique in 35 fractions over seven weeks. Prescribed doses were 70 Gy to macroscopic primary (T) and lymph node disease (N), and 63 Gy and 58.1 Gy to areas at high risk and at low risk of subclinical disease, respectively [16]. Concomitant chemotherapy consisted of cisplatin, three times weekly (100 mg/m 2 for three cycles every 21 days) or weekly (40 mg/m 2 for 6 cycles) [17].

MRI Protocols
MRIs were done on a 1.5 T scanner (Optima™ MR450w, GE Healthcare, Milwaukee, WI, USA) with a head and neck phased-array coil. All patients underwent three serial studies: at baseline, at the 10th fraction of RT and eight weeks after RT. MRI follow-up examinations were performed every six months for the first two years, and once per year afterwards. The pretreatment MRI protocol included coronal and axial T2-weighted images (field of view, 26-28 cm; acquisition matrix, 288 × 256; slice thickness, 4 mm), intravoxel incoherent motion diffusion-weighted imaging (IVIM-DWI), and dynamic-contrast enhanced MRI (DCE-MRI) sequences. IVIM-DWI was performed using nine b values (b = 0, 25, 50, 75, 100, 150, 300, 500, and 800 s/mm 2 , field of view 26 × 28 cm; acquisition matrix, 128 × 128; slice thickness, 4 mm; scanning time, 6 min 13 s). DCE-MRI was performed using a 3D fast-spoiled gradient echo sequence (flip angle, 30 • ; field of view, 28 cm; acquisition matrix, 128 × 128; slice thickness, 4 mm; spacing between slices, 2 mm; temporal resolution, 5 s; scanning time, 5 min). Three pre-bolus phases were acquired, followed by contrast-arrival phases after the intravenous administration of gadoliniumbased contrast agent, at a rate of 3 mL/s (60 total dynamic phases). To reduce the use of contrast medium, only IVIM-DWI was performed during treatment.

DCE-MRI and DWI Analysis
3D Slicer Software (version 4.6.2) was used for the lesion visualization and segmentation [18]. Commercial software (GenIQ General, GE Advanced Workstation, Palo Alto, CA, USA) was used to derive the quantitative maps from the DCE-MRI on the basis of a two-compartment pharmacokinetic model and automatic population-based selection of the arterial input function [19]. Three perfusion parameters were calculated at the single voxel level: K trans , defined as the transfer constant between plasma and the extravascular extracellular space (EES), and K ep , defined as the transfer constant between EES and plasma and v e , which represents the fractional EES volume. The baseline contours of the lesions, both T and the largest N, were performed on T2-weigheted images and, after rigid propagation, were transferred on the corresponding perfusion maps to perform quantitative analyses. The medians, percentiles (P) P10, P25, P75, and P90, skewness, kurtosis, energy, and entropy values were calculated from the voxel-based distribution of perfusion parameters within the entire lesion. The bin size used for each patient to extract the data was 0.2 min −1 , 0.4 min −1 and 0.025 for K trans , K ep and v e , respectively. The lesions at baseline and at the 10th fraction were also outlined on DWI at b = 800 s/mm 2 to extract the diffusion coefficients from the signal intensity curve at increasing b values, by means of home-made scripts developed in MATLAB (Release 2020b, MathWorks Inc., Natick, MA, USA). To reduce the instability of the diffusion-weighted signal within single voxels and increase the robustness of the quantifications, the median of the signal from the entire lesion at each b value was extracted and used for the fitting process. The conventional ADC was derived from data at b values of 0, 500, and 800 s/mm 2 , while the perfusion-free diffusion coefficient Dt was derived from data at b values of 300, 500, and 800 s/mm 2 , by a mono-exponential function [20]. The Levenberg-Marquardt algorithm was used to perform the fits. The perfusion fraction f was derived from an asymptotic estimate through an extrapolation of the signal intensity S0extr to b = 0 s/mm 2 from the above calculation of Dt: f = [(S0meas − S0extr)/S0meas] × 100. S0meas indicates the measured signal intensity at b value of 0 s/mm 2 [21].

Statistics
The endpoint of the study is the development of disease progression at the primary site, neck or distantly, alone or in combination at a follow up time of two years. Patients experiencing second primary tumors or death due to intercurrent disease who precluded a two year minimum follow up time were excluded.
Since an integrated statistical approach is suggested to improve data interpretability and prediction accuracy [22], data analyses were approached with both conventional statistical methods (to infer the relationships between variables) and ML algorithms (to at best design the prediction model). Regarding the former, univariate analyses on disease progression were performed considering various patient-, tumor-and treatment-related variables (p16 status, smoking habit, alcohol abuse, T subsite, T size, lymph node diameter, matted lymph nodes, cystic lymph nodes, ENE). Groups were compared with the chi-squared test or the Mann-Whitney rank test when appropriate. Statistical significance was claimed for p values < 0.05.

Machine Learning Analysis
Before the model building, among the above mentioned ones, the parameters with the best classification performance were selected using the minimum redundancy maximum relevance (MRMR) algorithm. The selection of the most informative predictors is a separate and mandatory step of the ML analysis pipeline to avoid the model overfitting and reduce to the minimum the subset of variables to make the best predictions, without loss of information [23]. The MRMR algorithm ranks both categorical and continuous parameters for classification and it finds the optimal set that was mutually and maximally dissimilar, based on the mutual information of variables [24]. The most common ML algorithms were quickly trained on our dataset to compare their performances and find the most appropriate [23]. Details are available in the Supplementary Materials. The model classification ability was evaluated in terms of accuracy, sensitivity, specificity, positive predictive value (PPV) and negative predictive value (NPV) after a stratified five-fold cross-validation to prevent overfitting and improve the possibility of generalizing the models. The prediction accuracies between different models were compared using the mid-p-value McNemar test. The ML analysis was carried out using the MATLAB environment.

Patient Characteristics
Out of the 97 enrolled patients, 68 (70.1%) were considered eligible. Eight patients without lymph node involvement were excluded. Twenty-one patients were not considered evaluable for disease response due to death due to intercurrent disease (N = 17) or to disease (N = 4) within 24 months from treatment end. There were 53 male patients (77.9%) and the median age was 61 years (IQR: 55-69 years). Selected patient and tumor characteristics are summarized in Table 1. Median follow-up was 33.2 months (IQR: 26.3-46.9 months). The majority of tumors (N = 53, 77.9%) were p16 positive. The distribution of primary subsites within the oropharynx was significantly different in p16 positive and negative patients (Table 1, p = 0.012). Similarly, smoking status and alcohol consumption were also distributed differently according to p16 status (p = 0.041 and p < 0.001 for smoking and drinking, respectively). P16 negative disease presented with a slightly smaller lymph node disease (median 2.5 cm vs. 1.6 in p16 negative patients, p = 0.06). Five patients refused chemotherapy and were treated with RT alone. At a median follow up of 9.7 months (IQR: 8.5-15.5 moths), disease progression was observed in 21 patients (30.8%). Interestingly, all the failures had been observed within two years from treatment end, and none afterwards. Regarding the pattern of disease progression, 10 patients (47.6%) developed distant metastases (lung metastasis in eight patients, bone in two patients); 11 patients (52.4%) had loco-regional disease progression (eight patients in lymph nodes, one in primary disease, three patients in both). Predictors for disease control following univariate analysis are shown in Table 2. Among all covariates, a significant correlation with disease control was found only for N diameter (p = 0.02). P16 status and ENE had a marginal impact on disease control (p = 0.134 and 0.073, respectively).

MRI Analysis and Prediction Models Driven by Machine Learning
The summary statistics of DCE-MRI and DWI parameters are reported in Tables 3-5.  Among radiological features, higher median values of K ep and lower median values of v e at both the primary tumor (p = 0.016 and p = 0.018 for K ep and v e , respectively) and the lymph nodes (p = 0.003 and p = 0.004 for K ep and v e , respectively) were associated with better disease control. Diffusion parameters did not correlate with clinical outcome, with the exception of ∆ADC (%), which showed a trend towards significance (p = 0.073). The MRMR algorithm identified K ep P90 and the diameter of N as best predictors among all the categorical and continuous parameters included in the analyses ( Figure S1). The performance of different families of ML algorithms is reported in Figure S2 (Tables S1 and S2). As illustrated in Figure S3, the best model is the one which incorporates nodal K ep P90 along with baseline N diameter. All other models lacking DCE-parameters were statistically inferior to the one including K ep P90 in terms of accuracy and/or specificity and/or sensitivity ( Figure S3). Table 5. Summary statistics of the diffusion parameters and their variations at the 10th fraction in primary tumor (T) and lymph node (N), for patients with and without disease control.

Discussion
In the present study, we found that baseline (perfusion) MRI features have an inde-

Discussion
In the present study, we found that baseline (perfusion) MRI features have an independent predictive role of disease outcome at two years at conventional statistics. Moreover, with a ML approach, the model incorporating selected DCE-MRI parameters along with clinical ones was the best predictor of outcome in terms of accuracy, specificity or sensitivity. Therefore, perfusion MRI features add to the clinical ones in predicting oncologic outcome after (chemo)radiotherapy for OPC. Recently, there has been an increasing interest in baseline or pre-treatment quantitative imaging features [25,26] with the purpose of customizing and further refining treatment strategies [27]. Similarly, to other recent investigations [28][29][30], we found several DCE-based parameters at both the primary tumor and the nodes to be correlated with outcome. In particular, higher median values of K ep and lower median values of v e were found to be associated with better tumor control rates. This is consistent with the fact that hypoxic tumors usually need more aggressive treatments [31] and are at higher risk of treatment failure [28]. In our dataset, the Decision Tree classification learner provided the best accuracies with respect to the other common ML algorithms. The histogram analysis of perfusion maps and the calculation of percentiles seemed to improve the capability to identify relationships between DCE-MRI parameters and outcome. Indeed, the model incorporating the parameter K ep P90, along with the nodal diameter, was the best predictor for disease progression among those tested ( Figure S3), providing a good accuracy, a high specificity and a fair sensitivity. These results are consistent with the current literature investigating prognostic functional imaging in patients with squamous cell carcinoma of the head and neck. In the review of Bos et al. [30], higher K ep values were associated to a higher probability of both overall survival and loco-regional control. Moreover, K ep was also identified as an independent predictor of disease free survival. Based on the pharmacokinetic model [19], the parameter K ep is positively associated to K trans but compared to K trans , it should be less influenced by blood volume/flow and more influenced by vessel permeability and the characteristics of the extravascular and extracellular compartments [32]. In our study, patients with disease control showed only a trend towards higher K trans (median) values while v e values were significantly lower in both primary tumors and lymph nodes. Being derived from the ratio K trans / v e , K ep was able to amplify these differences, showing a superiority compared to the other DCE-derived parameters. Concerning DWI features, we assumed that the perfusion-free diffusion coefficient Dt could be more appropriate to quantify the thermal diffusivity of water molecules in tissue and, consequently, to better indirectly evaluate the tissue cellular density and its early radiation-induced modifications compared to ADC. However, neither ADC nor Dt were significantly different between patients with disease control and progression, except for the nodal percentage variation of ADC during treatment. This parameter showed a trend towards significance, suggesting that a smaller increase in ADC was indicative of a higher risk of progression in accordance with previous investigations [33,34]. With regard to patient and tumor characteristics, it is a matter of fact that survival in patients with HNSCC is influenced by the presence of nodal metastases [35][36][37]. In HPV-related disease, there is a proven association between the infection and neck lymph-node disease [38]. Our finding is in line with the well-known clinical evidence; advanced N-disease staging is a well-established risk factor for disease progression [9]. Patients with high burden neck disease, particularly the presence of a matted node, have a higher risk of developing distant metastasis [39,40]. Lymph node size is also tightly correlated with regional recurrence [41]. Among morphological characteristics, we found a significant association between the nodal diameter and disease control, with a higher risk of progression for patients with larger pathologic lymph nodes. Concerning HPV status and ENE, we found a trend towards worse disease control in patients with a higher proportion of HPV-negative OPC and with ENE, though the difference was not statistically significant. Both HPV negative status and ENE are well-established factors of poor prognosis [42,43]. In the retrospective analysis of Ang et al. [3], patients with HPV-related disease showed a better three-year overall survival compared to HPV-negative patients, though factors such as heavy smoking history and advanced T and N stage had an independent detrimental effect on outcome. Despite a large number of HPV-positive patients, the vast majority of patients were smokers, thus limiting the benefits of HPV status. Similarly, even if the presence of ENE has been acknowledged as a poor predictor of outcome after radiotherapy for head and neck squamous cell carcinomas [10] and along with positive surgical margin, it represents an indication for chemo-radiotherapy in the post-operative setting [44]; its role in patients with HPV-positive disease is less clear. In the systematic review by Benchetrit et al. [45], radiological and pathological ENE in HPV-positive OPC leads to higher distant recurrence and worse overall survival but had a limited impact on locoregional recurrence. The present study has some limitations. In order to limit the impact of confounding factors on outcome (i.e., death from intercurrent causes), we selected only patients with a 2-yr minimum follow up. This strengthened the potential relationships between the various covariates and the oncologic outcome, though it reduced the sample size and thus the power of the study. Moreover, we focused only on a particular subset of head and neck cancer, oropharynx, mostly of which were actually virus-related. Therefore, generalization of the present findings to other head and neck subsites may not be appropriate. The accuracy of v e and K ep estimations may have been affected by the scan time of the DCE-MRI sequence. Indeed, a scan duration ≥10 min is recommended to allow the contrast agent to reach the equilibrium after leaking into the EES and leaking back to the vascular space [46]. However, such a long scan time is difficult to be used in clinical practice, particularly in head and neck cancer patients, due to their limited tolerance; moreover, the variability in DCE-MRI acquisition protocols and image post-processing in the literature has prevented, up to now, a harmonization between results from different institutions, limiting the replicability of our data. Furthermore, we did not include analyses on ADC/Dt maps to reduce the instability of the diffusion-weighted signal within single voxels; in future investigations we could explore the potential of a histogram-based approach to diffusion maps for enhancing the performance of the classification model.

Conclusions
Perfusional features at baseline MRI such as K ep and v e are predictive of tumor response at two years in patients affected by locally advanced oropharynx cancer and, when incorporated in a model along with clinical factors, help to build the predictive model with the highest performance in predicting clinical outcome. Therefore, DCE-MRI findings can potentially provide additional and unique information to offer more personalized treatment. Further studies on a larger population are needed to confirm our observations. Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers14102477/s1, Materials and Method: Details on model building; Figure S1: Bar plot of predictor importance scores after selection by MRMR algortithm; Figure S2: Bar plot of the accuracies obtained by different families of ML algorithms; Figure S3: Bar plot of predictive performances of various models obtained from different combination of the most relevant predictors; Table S1: List of predictors included in the different models; Table S2: Comparison between predictions of the models shown in Figure S3.
Author Contributions: All listed co-authors have made a substantial and intellectual contribution to the work; they were all involved with the conception, manuscript writing and final approval of the manuscript. Specific individual contributions are listed as follows: Conceptualization, methodology, validation, formal analysis, investigation, writing-original draft preparation, supervision: P.D., S.M. and G.S.; resources, data curation, project administration, funding acquisition: P.D., A.F., L.M., G.S.,