Quantitative Diffusion-Weighted Imaging Analyses to Predict Response to Neoadjuvant Immunotherapy in Patients with Locally Advanced Head and Neck Carcinoma

Simple Summary Immunotherapy may induce early treatment response in head and neck squamous cell carcinoma (HNSCC) for some patients. Routine imaging parameters fail to diagnose these responses; however, magnetic resonance (MR) diffusion-weighted imaging (DWI) may be able to do so. This study sought to correlate DWI parameters with treatment response early after immunotherapy treatment in HNSCC. We analyzed 24 patients with advanced HNSCC with imaging before and after the immunotherapy. We found that rounder tumors that were smaller in diameter before treatment were more likely to respond. A decrease in skewness of the tumor after treatment compared to before treatment, as well as an overall low skewness post-treatment, were linked to better treatment response. Though this study was explorative in nature, these results are promising for the predictive use of MR-DWI in HNSCC treated with immunotherapy. Abstract Background: Neoadjuvant immune checkpoint blockade (ICB) prior to surgery may induce early pathological responses in head and neck squamous cell carcinoma (HNSCC) patients. Routine imaging parameters fail to diagnose these responses early on. Magnetic resonance (MR) diffusion-weighted imaging (DWI) has proven to be useful for detecting HNSCC tumor mass after (chemo)radiation therapy. METHODS: 32 patients with stage II–IV, resectable HNSCC, treated at a phase Ib/IIa IMCISION trial (NCT03003637), were retrospectively analyzed using MR-imaging before and after two doses of single agent nivolumab (anti-PD-1) (n = 6) or nivolumab with ipilimumab (anti-CTLA-4) ICB (n = 26). The primary tumors were delineated pre- and post-treatment. A total of 32 features were derived from the delineation and correlated with the tumor regression percentage in the surgical specimen. Results: MR-DWI data was available for 24 of 32 patients. Smaller baseline tumor diameter (p = 0.01−0.04) and higher sphericity (p = 0.03) were predictive of having a good pathological response to ICB. Post-treatment skewness and the change in skewness between MRIs were negatively correlated with the tumor’s regression (p = 0.04, p = 0.02). Conclusion: Pre-treatment DWI tumor diameter and sphericity may be quantitative biomarkers for the prediction of an early pathological response to ICB. Furthermore, our data indicate that ADC skewness could be a marker for individual response evaluation.


Introduction
Head and neck squamous cell carcinoma (HNSCC) accounts for approximately 4% of worldwide cancer cases. Known risk factors for HNSCC development are the consumption of tobacco and alcohol, and infection with human papillomavirus (HPV) [1] Current curative treatment options for HNSCC patients are definitive radiotherapy with or without cisplatin-based chemotherapy [(C)RT], or primary surgery with or without adjuvant (C)RT. Despite these invasive treatments, advanced stage HNSCC patients have a poor prognosis [2,3]: their 5-year overall survival (OS) after extensive (salvage) surgery and adjuvant (C)RT is only 40-53%, thus underpinning the yet-unmet need for other treatment options [4].
Immune checkpoint blockade (ICB) of programmed cell death protein 1 (PD-1) therapy nearly triples the 2-year overall survival rate compared to investigator's choice in platinumrefractory, recurrent, or metastatic HNSCC patients [5,6]. In a curative setting, pre-surgery, neoadjuvant anti-PD-1 monotherapy alone or combined with cytotoxic T-lymphocyteassociated protein 4 (CTLA-4) ICB induces a major response, with >90% pathological tumor regression in 20-35% of HNSCC patients [7,8]. However, biomarkers predicting a response to neoadjuvant ICB are not available. Within the IMCISION trial, patients who developed a major pathological response (MPR) after neoadjuvant anti-PD-1 monotherapy or concurrent anti-PD-1 and anti-CTLA-4 all remained free of HNSCC at a median followup of two years [8], raising the question of whether extensive and potentially mutilating surgery is necessary after a major response to ICB. De-escalating or delaying surgery in this population, however, requires an accurate and, preferably, minimally invasive method to establish ICB response reliably. Evaluation by computed tomography (CT) and magnetic resonance (MR) imaging, according to the current response evaluation criteria in solid tumors (RECIST 1.1 [9]), unfortunately underestimates the frequency and depth of pathological response after neoadjuvant ICB in HNSCC [7,8,10,11], highlighting the need for other modalities to establish response early upon neoadjuvant treatment.
In IMCISION, 18[F]-fluorodeoxyglucose (FDG)-positron emission tomography (PET)based metabolic response evaluation has shown promise in detecting pathological responses in HNSCC patients treated with neoadjuvant ICB, using the delta of the calculated total lesion glycolysis (TLG) [12]. Although the accuracy at the primary tumor site is 94%, this technique is limited due to the relatively large tumor volume needed to assess TLG and FDG avidity for the influx of immune cells upon ICB treatment, which can lead to false positive results [13].
MR-imaging is a fundamental modality in clinical practice, providing extensive anatomical and functional information regarding soft tissue [14]. Even so, anatomic MRimaging performs poorly in distinguishing post-treatment effects from tumor recurrence after (C)RT [13]. Diffusion-weighted imaging (DWI) is a functional MR technique able to assess random-water (Brownian) motion, which is predominantly hindered by cellularity and, thus, negatively correlated with highly cellular tumor tissue [15,16] DWI has shown value in differentiating between tumor recurrence or benign, post-(chemo)radiotherapy effects in HNSCC [13,17]. Studies have shown the value of pre-treatment functional MR imaging in predicting treatment outcome or HNSCC after (C)RT [18,19].
However, whether DWI can differentiate HNSCC tumor cells from the influx of immune cells in the context of a major response to neoadjuvant ICB remains unknown. In this study, multiple quantitative DWI radiomic parameters were assessed, and their ability to predict pre-treatment or to diagnose early post-treatment of a major pathological response to neoadjuvant ICB was explored in HNSCC patients.

Patient Population and Trail Treatment Details
This retrospective study used the data from the 32 patients treated within the nonrandomized, open-label phase Ib/IIa IMCISION trial carried out at the Netherlands Cancer Institute (NKI) between February 2017 and October 2019 (NCT03003637). The IMCISION trial included adult patients with primary or recurrent, advanced stage (T2-T4, N0-N3b, M0) HNSCC of the oral cavity, oropharynx, hypopharynx, or larynx, eligible for curative (salvage) surgery. Staging was reported according to the 8th edition of the American Joint Committee on Cancer (AJCC) staging manual. All patients had a World Health Organization (WHO) Performance Score (PS) of 0 or 1. Critical exclusion criteria were the presence of distant metastases, a medical history of autoimmune disease, immunodeficiency, hepatitis B or C virus infection, or the use of immunosuppressive medication prior to treatment. Both HPV-negative and HPV-positive HNSCC patients were allowed to enter the trial. Complete in-and exclusion criteria, methods, and main results if the IMCISION trial have previously been published [8].
As part of the phase Ib trail safety run-in, the first 6 patients were treated with nivolumab monotherapy (240 mg) in weeks 1 and 3 (NIVO MONO), followed by surgery in week 5-6. Upon establishing feasibility, defined as the absence of surgical delay due to immune-related toxicity beyond week 6, the following 6-phase Ib patients were treated with a combination of nivolumab (240 mg) and ipilimumab (1 mg/kg) in week 1 and nivolumab (240 mg) in week 3 (COMBO), prior to surgery in week 5-6. When the COMBO regimen proved feasible, phase IIa opened and included 20 extra patients treated with the COMBO regimen prior to surgery in week 5-6. If indicated according to treatment guidelines, adjuvant (C)RT was performed. The institutional review board of the NKI reviewed and approved the IMCISION trial. All patients provided written informed consent prior to enrolment in IMCISION, which covered the use and presentation of the patient data enclosed in this manuscript [8].

Mr Imaging Acquisition
The MR-examinations were acquired on a 1.5 or 3.0T MRI scanner (Philips Healthcare) at baseline (week 0, 17 days (±10) pre-treatment) and after neoadjuvant ICB (post-treatment) at the end of week 4, 3 days (±0.6) before surgery, using a 16-channel head and neck coil combined with neck attachment. In 3 patients, the baseline MR imaging obtained in the referring institution was used. An overview of the varying examination parameters can be found in Appendix A.
The MR imaging protocol included 2D T1-weighted imaging (T1W) prior to contrast administration, a 2D T1W after contrast injection (T1Wc), and a diffusion-weighted imaging (DWI) sequence. Depending on the center and protocol, either a T2-weighted sequence (T2W), a T2W with short tau inversion recovery (STIR) image, or a T2W Spectral Attenuated Inversion Recovery (SPAIR) was included.
The T1W and T1Wc images were generated using a turbo spin echo (TSE), and had a slice thickness of 3 to 4 mm, an echo time (TE) of 10 ms, and a repetition time (TR) of 538 to 778 ms. Contrast was administered with a standard dose of 15 mL gadolinium (0.5 mM Dotarem, Guerbet, France), followed by a 20 mL saline flush.
DWI was acquired with b-values of either 0, 200, and 1000 s/mm 2 , or b-values of 0, 200, 400, 600, 800 s/mm 2 , depending on the protocol. The DWI had a slice thickness of 3 to 4 mm, a TE range between 67 and 80 ms, and a TR range between 3658 and 5255 ms. The apparent diffusion coefficient (ADC) map was automatically calculated on a pixel-by-pixel basis by the Philips Healthcare-provided MRI software using 3 b-values: 0, 200, and either 800 or 1000 s/mm 2 , depending on availability (Appendix A).

Postprocessing and Feature Extraction
All available MR-DWI scans were delineated by physician researcher (H.v.d.H.) under supervision of an experienced H&N radiologist (J.C.). Regions of interest (ROIs) were manually placed on the entire tumor volume using 3DSlicer software (Version 4.10.2; http://www.slicer.org [20]). ROIs were placed on the solid tumoral components according to high signal intensity on DW MR-images acquired by the highest available b-value (either 800 or 1000 s/mm 2 ), combined with corresponding low intensity on the ADC map [13]. Window and level values were kept consistent between patients. Large cystic or necrotic areas were excluded in order to focus on the viable tumor cells. T1Wc and STIR or SPAIR were used for anatomical correlation and tumor location. Researchers were blinded for treatment outcome. In Figure 1, an example of a delineated tumor is shown on DWI and ADC maps.
Signal intensity normalization was not needed for calculation of ADC map data, as this had already been derived from the DWI data. Image resampling was set to isotropic voxels of 2.0 mm. Bin width for the quantification of texture images was set to 5. Due to the variability in examination parameters and the small sample size, only first order and shape features were extracted using open-source package PyRadiomics 2.2.0 (Amsterdam, the Netherlands) [21]. Thus, a total of 32 radiomic features per patient was calculated from the ADC map data.

Outcome Assessment
The primary outcome of pathological response to neoadjuvant immunotherapy was assessed by an experienced head and neck pathologist (L.S.), using histological examination of the H&E-stained surgically resected specimen versus baseline biopsy specimen [8].
In the baseline biopsies and resection specimens, the percentage of viable tumor cells was calculated by measuring the area of viable tumor cells divided by the entire tumor bed area (defined by necrosis, fibrosis, keratinous debris, scarring, and immunoreaction sites).
Since previously (C)RT-treated, salvaged patients had already low viable tumor cell count at baseline biopsies, the change in viable tumor cell percentage from baseline biopsy to the post-treatment resection specimen was calculated in all patients, and, henceforth, was called the tumor regression percentage [22].
Patients with ≤10% residual viable tumor cells and 90-100% of tumor regression percentage were defined as major pathological responders (MPRs). Patients with ≤50% residual viable tumor cells and 50-89% regression percentage had a partial pathological response (PPR), and patients with any percentage of residual viable tumor cells, but <50% tumor regression percentage, had no pathological response (NPR) [23].
In contrast to patients with a PPR and NPR, patients with an MPR after neoadjuvant immunotherapy in IMCISION were characterized by an excellent clinical outcome, without recurrent disease, at a median follow-up of 2 years [8]. Consequently, patients with a primary tumor MPR were considered responders, whereas patients with PPR or NPR were considered non-responders. Two patients were not available for pathological evaluation as they did not undergo curative surgery [8]. These patients were classified according to their response as defined by the radiological RECIST-criteria on the post-treatment MRI, 10 to 16 days after the last ICB infusion, combined with clinical follow-up.

Statistical Analysis
Analyses were conducted using R (R Core Team (2022), R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria, https://www.R-project.org/, version 4.1.2). As the radiomics features are extracted from MR-imaging of different scanners, the stability of the radiomic features between the two predominately used MR-machines was assessed using the one-way ANOVA test. A p-value < 0.05 was considered significant for the instability of the feature. . DWI and STIR imaging of 2 patients before and after immunotherapy treatment with different treatment outcomes. Pt39, with a primary cT3N0 HNSCC of the oral cavity, is depicted before ((a1) b1000; (a2) ADC-map; (a3) STIR) and after ((b1) b1000; (b2) ADC-map; (b3) S = TIR) immunotherapy; this patient had MPR to treatment. Pt37 with a primary cT3N0 HNSCC of the oral cavity is depicted before ((c1) b1000; (c2) ADC-map; (c3) STIR) and after ((d1) b1000; (d2) ADC-map; (d3) STIR) immunotherapy; this patient had NPR to treatment.
Signal intensity normalization was not needed for calculation of ADC map data, as this had already been derived from the DWI data. Image resampling was set to isotropic voxels of 2.0 mm. Bin width for the quantification of texture images was set to 5. Due to the variability in examination parameters and the small sample size, only first order and shape features were extracted using open-source package PyRadiomics 2.2.0 (Amsterdam, the Netherlands) [21]. Thus, a total of 32 radiomic features per patient was calculated from the ADC map data.
Because of the small sample (<30 patients), and the fact that the data are not normally distributed, the means could not be compared using a Student's t-test. A Wilcoxon signed rank test would not provide effect sizes for the data. Regression methods were used to analyze potential associations between the radiomics features and the pathological response.
Regression analyses were applied for two classifications of pathological outcome; the continuous rate of pathological tumor regression (0-100%) was analyzed using linear regression, and the two defined responders (MPR) and non-responders groups (NPR + PPR) were analyzed using logistic regression. Considering the small number of patients and most clinical variables being of the categorical type, only 2 covariates (age and sex) were included in the multivariate analyses.
A separate regression analysis was performed separately for each feature at the baseline and post-treatment time points, and for the delta of each feature between the time points. The delta features were calculated only when both baseline and post-treatment data were available per patient. The significance level for the regression analyses was set at p < 0.05. Due to the explorative nature of the study, we chose not to correct for multiple testing.
Finally, we explored the predictive value of the pathological regression groups based on the evolution of all radiomics features together. To do so, logistic regression was used, and regularization penalization was optimized using elastic net. To minimize overfitting, cross-validation was applied by splitting the dataset on several folds. Considering the small sample size, a 3-fold cross-validation was used. The performance of the model for detecting response was assessed for each receiver operating characteristic (ROC) and area under the ROC curve (AUC), combined with specificity and sensitivity values.

Results
Two of the thirty-two patients enrolled in IMCISION had no DW MR imaging available, and were excluded. Six patients were excluded due to large artefacts or poor DWI quality at the tumor location on pre-and post-treatment imaging. Thus, twenty-four eligible patients with pre-treatment, post-treatment, or both time points' DWIs remained for analysis (total of 44 MRI-examinations) ( Figure 2). The mean age was 62 (±12) years, and 67% (n = 16) were male. All patients had oral cavity or oropharyngeal tumors. One patient had an HPVpositive HNSCC; all other tumors were HPV-negative. Seven patients were included with recurrent or residual HNSCC after previous (C) RT  Two of the thirty-two patients enrolled in IMCISION had no DW MR imaging available, and were excluded. Six patients were excluded due to large artefacts or poor DWI quality at the tumor location on pre-and post-treatment imaging. Thus, twenty-four eligible patients with pre-treatment, post-treatment, or both time points' DWIs remained for analysis (total of 44 MRI-examinations) ( Figure 2). The mean age was 62 (±12) years, and 67% (n = 16) were male. All patients had oral cavity or oropharyngeal tumors. One patient had an HPV-positive HNSCC; all other tumors were HPV-negative. Seven patients were included with recurrent or residual HNSCC after previous (C)RT (n = 3) or previous surgical treatment, with or without adjuvant RT (n = 4). On average, patients were diagnosed with a T3 tumor (45.8%), N0 stage (54.2%), and with disease stage III-IV (66.7%). Of the 24 patients, 7 had MPR to treatment, 15 patients had NPR, and 2 patients had PPR. See Table  1 for all baseline characteristics.  Baseline Characteristics (n = 24 Patients)

All
Responders Non-  From the ADC-map data, 32 radiomic features were extracted. After analyzing the stability of the data throughout the different scanners, 10 features were deemed unstable (p < 0.05) and were excluded from further analyses. (Appendix B) The remaining 22 stable features were analyzed for possible associations with the outcome. The mean values of these remaining features at the different time points are shown in Table 2.
Based on baseline imaging, responding patients had significantly smaller tumor diameters prior to therapy, measured from 3D and 2D dimensions on varying planes, compared to patients who did not respond (p = 0.04). Furthermore, these patients had an overall more spherical (round) tumors (p = 0.03) and lower entropy values (p = 0.05) than non-responders ( Figure 3, Table 3). Entropy specifies the randomness in the image values, where lower entropy means more homogeneous tissue [21]. These variables, except for tumor sphericity, were likewise associated with the continuous tumor pathological regression percentage ( Table 4). All findings were unaffected after correcting for sex and age (Tables 3 and 4). Table 2. Mean value of responding and non-responding tumor features at pre-and post-treatment timing, and the calculated delta features using the patients with data available at both time points. Pre-tr: pre-treatment, On-tr: post-treatment. *: Significant in analyses in Table 3.

Responders
Non Responders Upon post-treatment imaging after immunotherapy and prior to surgery, sphericity and entropy were no longer correlated with response. The post-treatment analysis did yield significant difference in three of the diameter measurements: the 3D diameter (p = 0.05), the diameter in the sagittal plane (row) (p = 0.04), and the major axis length (p = 0.04) between responders and non-responders (Table 3). Similar results were seen for the tumor pathological regression percentage. In addition, a new, significant difference was seen for skewness of the ADC, as negative skewness was associated with higher pathological response percentage in the univariate analyses (p = 0.04) as well as in the multivariate analyses (p = 0.05). (Table 4) Table 3. Features analyses of responder groups at pre-treatment, post-treatment, and delta time points.   In the multivariate analysis of the delta between the features collected from pre-and post-treatment DWI, ADC skewness was the only significant feature that was correlated to tumor regression percentage upon ICB (p = 0.02) ( Table 4). However, skewness was not correlated to either of the two response groups (p = 0.07) ( Table 3).

Pre-Treatment
As the degree of pathological tumor regression is based on changes in viable tumor cells in the tumor bed between the baseline biopsy and the surgical specimen, we would expect the delta of the features post-treatment to depict this regression most accurately. (A) Sphericity or tumor roundness is significantly higher for responders at pre-treatment; this dissipates at post-treatment imaging. (B) Significantly lower levels of entropy can be seen for responders compared to non-responders, but only at pre-treatment imaging. (C) The maximum diameter measured using 3D dimensions is significantly lower at pre-treatment imaging for responders compared to non-responders. (D) No significance of ADC skewness can be seen between the two response groups at pre-or post-treatment. A near-significant result (p = 0.066) is observed for the delta of ADC skewness.
In the multivariate analysis of the delta between the features collected from pre-and post-treatment DWI, ADC skewness was the only significant feature that was correlated to tumor regression percentage upon ICB (p = 0.02) ( Table 4). However, skewness was not correlated to either of the two response groups (p = 0.07) ( Table 3).
As the degree of pathological tumor regression is based on changes in viable tumor cells in the tumor bed between the baseline biopsy and the surgical specimen, we would expect the delta of the features post-treatment to depict this regression most accurately.
However, except for ADC skewness, no significant associations between the delta features were observed. Nonetheless, combining all delta features together in a model did yield a considerable AUC of 0.846 (0.716-0.977) in predicting responders from all non-responders (Appendix C).

Discussion
This study aimed to explore the ability of MR diffusion-based imaging parameters, extracted from MR-imaging acquired at baseline and after immunotherapy treatment (shortly prior to surgery), to predict or detect major pathological responses to neoadjuvant ICB early on in patients with resectable HNSCC.
Several differences were seen at the baseline between the pathological response groups in our data. Responding tumors were characterized by a significantly smaller diameter, greater sphericity, and lower entropy values.
Higher sphericity (tumor roundness) has previously been correlated to improved progression-free survival (PFS) and local control after (C)RT in HNSCC [24][25][26]. Sphericity is associated with an expansive pattern of tumoral growth often seen in HPV-positive HNSCC, in contrast to the infiltrative growth pattern more common in HPV-negative tumors [24,27,28]. As baseline DWI was unavailable for the only HPV-positive patient in this dataset, the higher sphericity observed in the ICB-responsive population was fully reflective of HPV-negative HNSCC. In this cohort, necrotic area(s) were excluded in the delineation, creating a more complex and, thus, less spherical shape in necrotic tumors. This is illustrated in Figure 4a,b. This may result in possible confounding of the tumor sphericity, as this feature may now also be negatively associated with the baseline presence of necrotic area(s) within the tumor. As poor perfusion may impede the delivery of intravenous drugs [29], necrotic (non-spherical) tumors may have been overrepresented in the nonresponding group. In addition, the inclusion of patients with recurrent or residual disease after prior RT may also have influenced the sphericity, due to previous treatment effects on the tumoral area resulting in regional necrosis or fibrosis, though our data showed no clear trend for sphericity.
On the contrary, the highest entropy values of all evaluated tumors were discernable for the previously treated cancers. Entropy is a statistical measure of the randomness in the image values used to characterize the texture of the input image [21]. A lower level of entropy is defined by lower levels of chaos and randomness, and resembles more homogenous tissue, whereas higher entropy is linked to heterogeneity of the tissue. Though heterogeneous tumors have also been linked to higher chances of local failure, our data may have been confounded by the inclusion of previously irradiated tumors [25,30]. Finally, the diameter parameter, measured in different planes and dimensions, was significantly different between the two response groups. Correlations between pre-treatment smaller tumor volume and increased locoregional control or PFS have often been described for (C)RT treated patients [31][32][33][34]. Initially, it was hypothesized that a significant change in volume would occur based on the applied delineation method of low ADC combined with high b1000/800 areas, which should result in a decline in tumor size at post-treatment as a result of a decrease in tumor cellularity density upon response [13]. The maximum diameter can even be fairly reliably used for tumor volume estimation in anatomical MRIs [35]. Yet, in contrast to the diameter, no significance was seen for tumor volume throughout our dataset, though a lower baseline p-value (p = 0.1) was discernable using the tumor regression percentage (Table 4). When considering the substantial standard deviation of the volume parameter, as shown in Table 2, it becomes apparent that the range of the volume parameter is too broad for this small dataset to be significant, though an overall smaller tumor volume can be seen in responders compared to non-responders. The diameter appears to be less variable between patients, and is, therefore, significant within this cohort.
At post-treatment imaging, most significant correlations dissipated, with the exception of the newly found parameter skewness of the ADC and the maximum diameter in 3D, the sagittal plane and major axis. Intriguingly, a recent study describes a moderate inverse correlation between ADC skewness and programmed death-ligand 1 (PD-L1) expression scores [36] Though this correlation alone is not strong enough to predict PD-L1 expression in a clinical routine, it does indicate the ADC-map may depict more complex histopathology than just the viable tumor cells. Moreover, this correlation is especially interesting since nivolumab functions as a PD-1 inhibitor. A significantly lower post-treatment skewness or a decrease in the skewness during treatment, as seen with the delta analyses, may have some usefulness in treatment monitoring in the future [8]. However, within this dataset, significance was only observed for the continuous tumor regression percentage, and not as clearly in the two defined response groups. The maximum diameter, measured in several dimensions and planes, was also significant post-treatment, but less so than at baseline. This could be due the challenges of delineating post-treatment, as is shown by increased variation between delineation post-treatment, and as illustrated in Figures 1 and 4.  Figure 4a and b depict pre-immunotherapy imaging of a patient with a residual tumor (rT3N0) after previous CRT with a necrotic area (n) within the tumor region (t) on T1Wc (a), and as shown on ADC-map (b). On the T1Wc imaging, a possible infiltrate (i) could be included in the tumor area, marked by t(+i), as this is difficult to discern. This patient did not respond to treatment. Figure 4c and d depict post-treatment imaging of a cT4N1 primary tumor of the oral cavity, with possible immune infiltrate (i) surrounding the tumor (t), or, conceivably, within the tumor area (t(+i)) on T1wc (c) and ADC-map (d). This patient had a partial pathological response to treatment.
At post-treatment imaging, most significant correlations dissipated, with the exception of the newly found parameter skewness of the ADC and the maximum diameter in 3D, the sagittal plane and major axis. Intriguingly, a recent study describes a moderate inverse correlation between ADC skewness and programmed death-ligand 1 (PD-L1) expression scores [36] Though this correlation alone is not strong enough to predict PD-L1 expression in a clinical routine, it does indicate the ADC-map may depict more complex histopathology than just the viable tumor cells. Moreover, this correlation is especially  Figure 4a and b depict pre-immunotherapy imaging of a patient with a residual tumor (rT3N0) after previous CRT with a necrotic area (n) within the tumor region (t) on T1Wc (a), and as shown on ADC-map (b). On the T1Wc imaging, a possible infiltrate (i) could be included in the tumor area, marked by t(+i), as this is difficult to discern. This patient did not respond to treatment. Figure 4c and d depict post-treatment imaging of a cT4N1 primary tumor of the oral cavity, with possible immune infiltrate (i) surrounding the tumor (t), or, conceivably, within the tumor area (t(+i)) on T1wc (c) and ADC-map (d). This patient had a partial pathological response to treatment. We constructed a model that yielded a high AUC, underlining the potential value of these parameters when assessed in concert, rather than individually. However, this model was trained and tested on the same limited dataset, and, therefore, was likely overfitted. Without external validation, its value remains uncertain. Recent work by Corino et al., however, promisingly illustrated a higher AUC of a CT radiomics model compared to a clinical model for predicting 10-month OS to nivolumab treatment of HNSSC-patients [37]. This highlights the potential for models based on MR-radiomics that could additionally provide extensive data on soft tumoral tissue.
This retrospective analysis had several limitations, such as the small sample size, making the study explorative in nature. Within this limited dataset, previously (C)RTtreated tumors were also included. These are radiologically different from primary tumors, which may have diluted the results. Furthermore, the MR-imaging was not all acquired according to the same protocol, on the same MR-machine, or in the same treatment center. While stability between the two main MRI types and field strength was established, and instable features excluded, minor variations may still have limited the interpretability of the results. In addition, some patients were excluded based on artefacts, indicating that a possible MRI model will not be applicable for everyone. No correlation between molecular pathological biomarkers and imaging markers was explored due to the limited sample size.
Neoadjuvant immunotherapy trials employing combined aPD1 and aCTLA4 immunecheckpoint blockade prior to extensive curative surgery and radiation therapy provided promising results, namely a MPR rate of 30% upon immunotherapy at the time of surgery, at the primary tumor site [7,8]. This study examined a novel cohort of immunotherapy-treated HNSCC, and attempted to explore the possible value of DWI for response monitoring after immunotherapy in HNSCC. In view of future clinical trials aiming at de-escalation of surgery in these patients, biomarkers to identify these patients with a favorable response upon immunotherapy prior to surgery are needed. RECIST 1.1 [9] criteria used for (chemo) radiation-type treatments do not currently suffice to monitor treatment response in immunotherapy, suggesting a different imaging approach for immunotherapy response monitoring [7,8,10,11]. Though more research in larger datasets is required to definitively link DWI parameters to immunotherapy outcomes, this study might provide guidance for the direction of such research.

Conclusions
In conclusion, while this study has certain limitations, our analysis of DWI features identified a significant association between baseline tumor diameter and sphericity and pathological response after neoadjuvant nivolumab and ipilimumab in HNSCC. In addition, ADC skewness may play an important role in pathological response evaluation, either as a stand-alone post-treatment value or calculated as the difference between baseline and post-treatment imaging.  Informed Consent Statement: Informed consent was already acquired from all patients prior to enrolment in the IMCISION-trail (NCT03003637), which covered the use and presentation of patient data enclosed in this manuscript.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest: C.L.Z. is linked to Investigator-initiated clinical trials in collaboration with
Bristol Myers Squibb. This funder had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results. All other authors have no conflicts of interest to declare.