Multiparametric Modelling of Survival in Pancreatic Ductal Adenocarcinoma Using Clinical, Histomorphological, Genetic and Image-Derived Parameters

Rationale: Pancreatic ductal adenocarcinoma (PDAC) remains a tumor entity of exceptionally poor prognosis, and several biomarkers are under current investigation for the prediction of patient prognosis. Many studies focus on promoting newly developed imaging biomarkers without a rigorous comparison to other established parameters. To assess the true value and leverage the potential of all efforts in this field, a multi-parametric evaluation of the available biomarkers for PDAC survival prediction is warranted. Here we present a multiparametric analysis to assess the predictive value of established parameters and the added contribution of newly developed imaging features such as biomarkers for overall PDAC patient survival. Methods: 103 patients with resectable PDAC were retrospectively enrolled. Clinical and histopathological data (age, sex, chemotherapy regimens, tumor size, lymph node status, grading and resection status), morpho-molecular and genetic data (tumor morphology, molecular subtype, tp53, kras, smad4 and p16 genetics), image-derived features and the combination of all parameters were tested for their prognostic strength based on the concordance index (CI) of multivariate Cox proportional hazards survival modelling after unsupervised machine learning preprocessing. Results: The average CIs of the out-of-sample data were: 0.63 for the clinical and histopathological features, 0.53 for the morpho-molecular and genetic features, 0.65 for the imaging features and 0.65 for the combined model including all parameters. Conclusions: Imaging-derived features represent an independent survival predictor in PDAC and enable the multiparametric, machine learning-assisted modelling of postoperative overall survival with a high performance compared to clinical and morpho-molecular/genetic parameters. We propose that future studies systematically include imaging-derived features to benchmark their additive value when evaluating biomarker-based model performance.


Introduction
Pancreatic ductal adenocarcinoma (PDAC), despite its relative rarity, remains among the deadliest tumor entities in the developed world. For instance, PDAC is the 4th leading cause of cancer-related death while representing only 3% of newly diagnosed cancer cases in the United States [1]. A large body of research into therapeutic targets, and newly introduced therapy regimens, have hitherto not been able to improve overall PDAC patient prognosis beyond a five-year survival of about 9% [2]. Even in the setting of tumor resection and adjuvant chemotherapy [3], patients often develop therapy-resistant tumor recurrence or metachronous metastatic disease. The key driver of such therapy escape phenomena is likely to be tumor heterogeneity, both on the genetic level, with the four driver genes kras, smad4, tp53 and p16 [4,5] having been shown to lead to distinct survival outcomes [6,7]. On the molecular/histo-morphologic level, quasi-mesenchymal and epithelial tumors demonstrate a distinct chemotherapy response [8,9]. Hence, despite the overall poor prognosis in PDAC, significant differences in individual patient survival are noted. Therefore, outcome prediction and patient stratification based on available parameters represent important goals in clinical patient care. The current gold standard for tumor assessment, biopsy and histopathology, entail a major risk of misclassifying tumors due to undersampling and sampling errors. Moreover, non-invasive biomarkers have yet to reach clinical maturity. Recently however, evidence has emerged that the machine learning-based analysis of pre-therapeutic imaging can provide a decision guidance based on the non-invasive derivation of quantitative, whole-tumor characteristics [10,11]. Image pattern analysis and machine learning approaches (often termed radiomics) have yielded encouraging results in the prediction of molecular PDAC subtypes and of patient survival from magnetic resonance or computed tomography imaging [12,13]. Much of recent research, driven by the increase in attention towards the field of artificial intelligence, has concentrated exclusively on proving the superiority of radiomics to human observers or other diagnostic and prognostic parameters [14][15][16]. However, in a disease as complex as PDAC, the integration of clinical information, invasive biomarkers and imaging-derived data may provide a more comprehensive prognostic assessment of the tumor than any single modality [17]. Therefore, a fair comparison elucidating the differential contribution of each parameter type and their combined value, is warranted.
In this study, we present a multimodal, data-driven workflow including clinical information, histo-morphologic/genetic parameters and computed tomography-derived imaging biomarkers for modelling overall patient survival in the post-operative setting.

Experimental Section
The study was designed as a retrospective cohort study. The study was approved by the institutional ethics review board (Ethics Committee of the Technical University of Munich Faculty of Medicine) (Protocol Number 180/17S; date of approval: 9 May 2017), and the requirement for individual written consent was waived. All procedures were carried out in accordance with pertinent laws and regulations, as well as the Helsinki declaration. The STrengthening the Reporting of OBservational studies in Epidemiology (STROBE) checklist [18] and patient inclusion flowchart are included in the supplemental material (Table S1 and Figure S1). In brief, 177 patients with a confirmed diagnosis of PDAC were considered for inclusion. Patients without baseline computed tomography (CT) imaging, incomplete imaging, insufficient image quality (due to motion artifacts or significant beam hardening from adjacent foreign matter), patients who died sooner than six weeks after surgery or those with incomplete clinical data were excluded. A total of 103 patients were analyzed. All patients underwent tumor resection. The follow-up interval began on the 10 October 2006 and ended on the 14 April 2019. Clinical data were obtained using the hospital's information system, and survival information was obtained from the hospital's information system and from the national cancer registry. The genetic and histopathological analyses were carried out by the Department of Pathology of the Technical University of Munich and were available prior to the study commencement. Image-derived parameters data were collected during the analysis.
The following clinical data was obtained for all patients: sex, age at diagnosis, tumor site (head, body or tail of the pancreas), type of adjuvant chemotherapy received (gemcitabine-based vs. no chemotherapy), time to progression, type of progressive disease (local recurrence, metastasis, both), first-line chemotherapy after progression (gemcitabine, FOLFIRINOX, none/best supportive care) and overall survival time including censoring.
CT images were exported from the hospital picture archive and processed as previously described [13,21] and according to current best practices. Segmentation was performed fully automatically by in-house developed software and manually corrected by an expert observer. Image-derived features were extracted using PyRadiomics version 2.2.0 [22] with standard settings. 1409 image features were extracted. Features with a variance of less than 1e-5 were excluded, yielding 1384 radiomic features which were included in the analysis. Feature values were normalized to the (0, 1) interval. A detailed description of the radiomic extraction process can be found in the supplement (File S1). The PyRadiomics extraction settings file (.yml) can also be found in the supplement (File S2).
Due to the very large number of features compared to the cohort size, an unsupervised machine learning-based dimensionality reduction and feature selection were performed as follows: linear principle component analysis (PCA) was employed for radiomic features with the aim of capturing 99% of the parameter variance. A univariate pre-selection was performed on all features using the log-rank-test statistic with a cut-off of p < 0.1.
The Cox Proportional Hazards model was fit to each parameter group separately at first to assess the model concordance index (CI) after asserting that the proportional hazards assumptions were met. Parameters which did not meet the proportional hazards assumptions were excluded from the final analysis of the combined features. For testing the predictive performance, a five-fold hold-out cross-validation was employed in all analyses, and the concordance indexes across the five folds were averaged. 95% confidence intervals were calculated using 100-fold bootstrap resampling. All statistical analyses were performed in Python v. 3.7.6. using the packages lifelines and scikit-learn. A two-sided significance level of p < 0.05 was chosen.

Results
A total of 103 patients were included in the study. Their clinical characteristics are summarized in Table 1. Of the 103 patients included, 51 patients experienced a progressive disease during the observation period, and of these 14 (27%) experienced local recurrence, 29 (57%) liver metastasis and 8 (16%) both. Of those patients, 22 (43%) received first-line gemcitabine, 13 (25%) received first-line FOLFIRINOX and 16 (31%) received best supportive care. The median progression-free survival was 7.3 months. The disease progression (experienced vs. did not experience) or type of treatment after progression were not significantly associated with overall survival and are thus not included in the final analysis.
The Cox proportional hazards method was used for modelling the overall survival using each of the parameter categories. The median overall survival was 16.6 months, and 19/103 patients (18.4%) were censored in the final survival analysis.
The modelling of parameter group 1, i.e., the clinical and standard histopathological data (age, sex, chemotherapy, pT, pN, G, R) achieved a CI of 0.63 (95% conf. int. 0.60-0.66). The nodal status (pN), tumor grading (G), tumor size (pT) and chemotherapy regimen were significantly associated with overall survival (Table 2, Figure 1). The Cox proportional hazards method was used for modelling the overall survival using each of the parameter categories. The median overall survival was 16.6 months, and 19/103 patients (18.4%) were censored in the final survival analysis.
The modelling of parameter group 1, i.e., the clinical and standard histopathological data (age, sex, chemotherapy, pT, pN, G, R) achieved a CI of 0.63 (95% conf. int. 0.60-0.66). The nodal status (pN), tumor grading (G), tumor size (pT) and chemotherapy regimen were significantly associated with overall survival (Table 2, Figure 1).  In parameter group 2, i.e., the morpho-molecular and genetic data (tumor morphology, molecular subtype, and pt53, kras, smad4 and p16 genetics), the tumor morphology did not meet the proportional hazard assumptions and was therefore excluded from the combined analysis. None of the parameters were significantly associated with overall survival (Table 3, Figure 2). The model resulted in a CI of 0.53 (95% conf. int. 0.47-0.59).  In parameter group 2, i.e., the morpho-molecular and genetic data (tumor morphology, molecular subtype, and pt53, kras, smad4 and p16 genetics), the tumor morphology did not meet the proportional hazard assumptions and was therefore excluded from the combined analysis. None of the parameters were significantly associated with overall survival (Table 3, Figure 2). The model resulted in a CI of 0.53 (95% conf. int. 0.47-0.59).  The principal component analysis of the imaging features, i.e., parameter group 3, resulted in 71 unique image-derived feature groups representing 99% of the total parameter variance. Their modelling was preceded by a univariate pre-selection, after which eight image feature groups were retained for inclusion in the multivariate model. Of these, four were significantly associated with overall survival (Table 4, Figure 3), three negatively and one positively. The multivariate model achieved a CI of 0.65 (95% conf. int. 0.62-0.69).  Finally, to assess the predictive performance of the combined feature set and model feature interactions, the features from parameter groups 1 through 3 which were found to obey proportional hazards assumptions were jointly included in a Cox proportional hazards model. In total, seven features, namely three image-derived parameter groups, the tumor subtype, tp53 mutation status, The principal component analysis of the imaging features, i.e., parameter group 3, resulted in 71 unique image-derived feature groups representing 99% of the total parameter variance. Their modelling was preceded by a univariate pre-selection, after which eight image feature groups were retained for inclusion in the multivariate model. Of these, four were significantly associated with overall survival (Table 4, Figure 3), three negatively and one positively. The multivariate model achieved a CI of 0.65 (95% conf. int. 0.62-0.69).  The principal component analysis of the imaging features, i.e., parameter group 3, resulted in 71 unique image-derived feature groups representing 99% of the total parameter variance. Their modelling was preceded by a univariate pre-selection, after which eight image feature groups were retained for inclusion in the multivariate model. Of these, four were significantly associated with overall survival (Table 4, Figure 3), three negatively and one positively. The multivariate model achieved a CI of 0.65 (95% conf. int. 0.62-0.69).  Finally, to assess the predictive performance of the combined feature set and model feature interactions, the features from parameter groups 1 through 3 which were found to obey proportional hazards assumptions were jointly included in a Cox proportional hazards model. In total, seven features, namely three image-derived parameter groups, the tumor subtype, tp53 mutation status, Finally, to assess the predictive performance of the combined feature set and model feature interactions, the features from parameter groups 1 through 3 which were found to obey proportional hazards assumptions were jointly included in a Cox proportional hazards model. In total, seven features, namely three image-derived parameter groups, the tumor subtype, tp53 mutation status, tumor size (pT) and nodal status (pN), were found to significantly and independently predict overall survival (Table 5, Figure 4). The final combined model achieved a CI of 0.65 (95% conf. int. 0.60-0.69). tumor size (pT) and nodal status (pN), were found to significantly and independently predict overall survival (Table 5, Figure 4). The final combined model achieved a CI of 0.65 (95% conf. int. 0.60-0.69).

Discussion
We here present the results of a comprehensive multiparametric analysis for the prediction of overall patient survival in a group of resected PDAC patients. Our results underscore the potential of image-derived biomarkers compared to both conventional histopathological and novel molecular, morphologic and genetic phenotyping in the context of patient risk stratification in PDAC.
In our previous work, we reported on survival prediction based on imaging biomarkers alone and developed models for the stratification of molecular phenotypes, therapy response and patient survival from computed tomography and diffusion-weighted magnetic resonance images [12,13,21]. Our current work expands on these efforts by including broadly available clinical parameters and

Discussion
We here present the results of a comprehensive multiparametric analysis for the prediction of overall patient survival in a group of resected PDAC patients. Our results underscore the potential of image-derived biomarkers compared to both conventional histopathological and novel molecular, morphologic and genetic phenotyping in the context of patient risk stratification in PDAC.
In our previous work, we reported on survival prediction based on imaging biomarkers alone and developed models for the stratification of molecular phenotypes, therapy response and patient survival from computed tomography and diffusion-weighted magnetic resonance images [12,13,21]. Our current work expands on these efforts by including broadly available clinical parameters and unestablished, novel biomarkers in multivariate survival models. Interestingly, both conventional clinical parameters, such as tumor stage, nodal status and adjuvant chemotherapy, and the imaging-derived biomarkers introduced in this study outperform morpho-molecular and genetic parameters alone. We believe this to be an effect of the parameter distribution in this cohort and the overall cohort size. Both previous work and our own multivariate analysis have shown genetic parameters [20] and molecular tumor subtypes [8,19] to significantly predict survival. Furthermore, we believe sampling errors of the histo-morphological and genetic workup to reduce the representativeness of invasive biomarkers, and we thus highlight the necessity for a whole tumor analysis. Despite the fact that imaging-derived biomarkers do not outperform other parameters (0.63 versus 0.65 CI), they are justified due to their non-invasive nature, which applies especially to the palliative and neo-adjuvant settings.
The recent study by Zhang et al. [17] showed a similar survival model concordance index of around 0.65 on validation data by applying convolutional neural networks (CNN) to CT imaging data. The non-inferiority of our approach compared to the more advanced CNN algorithm might be a result of the limited sample size, a common problem in current imaging studies [23].
To our knowledge, this is the first study to systematically evaluate a broad range of currently available PDAC risk biomarkers. The sample size studied, and the absence of an external validation cohort, are a limitation of our study but also an immediate consequence of the high financial and personnel cost of the extensive morpho-molecular and genetic workup. Considering the lack of re-imbursement for the detailed tissue analyses presented, biomarkers derived from routine clinical data-including imaging-are a viable option in the setting of limited healthcare resources. Nevertheless, our results should be tested prospectively to ascertain their external validity.

Conclusions
In summary, this study presents a further step towards the validation of multi-parametric analyses in the context of clinical patient care. Two conclusions arise from our results: A side-by-side direct comparison of whole tumor tissue and imaging-data derived features will be required to adequately compare the true predictive performance of imaging against morpho-molecular tumor constituents. Furthermore, for multiparametric data integration-as demonstrated in our work-to reach clinical maturity, multi-centric cooperation will be required to attain sufficient sample sizes.