This retrospective study received approval from the local institutional review board of the Peking University First Hospital Medical Science Research ethics committee (IRB number: 2021-342), and the requirement for written consent was waived.

A total of 2,232 patients suspected of PCa between 2016 and 2020 in our institution were included in the study. The inclusion criteria were as follows: (a) Availability of pretreatment MR images in the Picture Archiving and Communication Systems (PACS). (b) Clinically diagnosed as advanced stage and initial treatment with RT, HT, or a combination of both. (c) Regular follow-up at least every three months in the first year and every six months in the second year. (d) Availability of clinical data. (e) Follow-up period of at least two years with documented biochemical recurrence (BCR+) or non-BCR (BCR-). BCR was defined as any PSA increase greater than 2 ng/mL compared to the PSA nadir value for patients who underwent RT with or without HT. For patients who received only HT, BCR was defined as castrate serum testosterone less than 1.7 nmol/L with either three consecutive PSA increases at least one week apart, resulting in a two-fold increase exceeding the nadir by 50% and a PSA greater than 2 ng/mL, or radiological progression evidenced by the appearance of new lesions (26).

Finally, a total of 131 consecutive patients were included in this study, consisting of 100 BCR- and 31 BCR+ patients (Figure 1). Of these patients, 4 received only RT (including 4 BCR- and 0 BCR+), 75 received RT with HT (74 BCR- and 1 BCR+), and 52 received only HT (22 BCR- and 30 BCR+).

A comprehensive description of acquisition parameters for seven different MR scanners are provided in Table 1. The magnetic field strength of MR scanners included 1.5 T (n = 46, 35.1%) and 3.0 T (n = 85, 64.9%). The radiomics analysis exclusively focused on ADC images, which were obtained with a median b value of 1400 s/mm². Diffusion weighted imaging (DWI) was acquired by using single-shot echo planner imaging (SS-EPI) sequence. The calculation of ADC maps was performed using vendor-specific software associated with each MR scanner (GE: Advantage Workstation, Philips: IntelliSpace Portal, Siemens: syngo.via, UIH: uWS-MR Advanced Postprocess Workstation). Other imaging sequences, including T1-weighted imaging (T1WI), T2-weighted imaging (T2WI) and dynamic contrast enhancement (DCE), were acquired simultaneously but were not subjected to analysis in current study.

The areas of PCa were predicted by priorly trained models (22) on the ADC maps (Figure 2). Initially, the areas corresponding to prostate gland were segmented by artificial intelligence (AI) models. Subsequently, the regions of suspected PCa were segmented in a sequential manner. The three-dimensional volume ROI was utilized for extracting image features required for model development. If multiple PCa foci were segmented, the largest one was automatically taken as the ROI. Prostate Imaging and Reporting Archiving Data System (PI-RADS) v2.1 was used to score lesions of interest by an experienced radiologist (work experience more than 10 years). Among 131 patients, 106 had only one suspicious lesion, 20 had two, 4 had three, and 1 had four. The volume, location, and PI-RADS scores of these lesions are detailed in Table 2. In 121 patients, the largest lesion coincided with the one rated as the highest PI-RADS score. In 10 cases, the largest lesion shared the same PI-RADS score as the second largest lesion.

The collected data were divided into two datasets, namely, the training dataset (n = 93) and the test dataset (n = 38), using a random allocation with a ratio of 7:3 (27). Subsequently, three prediction models were developed using the training dataset: a clinical model, a conventional radiomics model, and a deep-radiomics model. The clinical characteristics, including age, PSA level, PI-RADS score, International Society of Urological Pathology (ISUP) group for biopsy pathology, clinical stage and treatment, remained consistent across all models. However, the image features varied among models, and their specific extraction methods will be elaborated in the subsequent manuscript.

To construct the clinical model, location, volume and ADC value of PCa area segmented by the AI models were selected as image features. Univariable logistic regression and multivariable logistic regression analysis were conducted to identify the significant predictors. The selection of predictors was performed using a stepwise selection approach guided by the Akaike information criterion (AIC). The adequacy of the model was evaluated by calculating the R2 value using the Nagelkerke method. To help comprehensive interpretation of the model, a nomogram graph was generated.

To construct the conventional radiomics model, the image features were extracted from ROIs on the ADC maps using the Pyradiomics package in Python (https://pyradiomics.readthedocs.io/en/latest/changes.html). To account for potential variations across different scanners, preprocessing of the ADC maps was performed by applying image normalization to mitigate confounding effects. Three types of images were analyzed: “Original Images” representing unfiltered images, “LoG Images” obtained by applying the Laplacian of Gaussian filter, and “Wavelet Images” generated through a three-dimensional wavelet transformation using the PyWavelet package in the x, y, and z directions. The ROIs were then preprocessed to ensure a consistent size. A total of 14 shape features, 18 first-order statistical features, and 24 texture features were extracted from the images. The shape features were exclusively extracted from the “Original Images,” while the first-order statistical and texture features were extracted from all three types of images. Therefore, a total of 14 shape features, 216 first-order statistical features, and 840 texture features were obtained. The mathematical expressions and semantic meanings of the features extracted can be found at https://pyradiomics.readthedocs.io/en/latest/. After feature extraction, several additional steps were performed. Z score normalization was applied to rescale the extracted features, and Pearson correlation coefficients (PCCs) were calculated to identify highly correlated features. Features with a PCC value exceeding 0.99 were removed to mitigate multicollinearity. Analysis of variance (ANOVA) was then used to select features for the final radiomics model. As the classifier, eXtreme gradient boosting (XGBoost) algorithm was employed. And SHAP (SHapley Additive exPlanations) values were used to interpret the results obtained from the XGBoost models.

To construct the deep-radiomics model, the image features were extracted by employing a deep learning algorithm. The construction process involved several steps. First, the intensities of ADC maps were preprocessed by normalizing them. Second, ROIs were resampled to ensure a consistent voxel size. Third, a pre-trained deep learning model was utilized to extract features from ROIs, leveraging the power of learned representations. Specifically, the ROIs were input into the convolutional layers of the MedicalNet architecture, which had been initialized with pre-trained weights to extract discriminative features. A comprehensive description of the network can be found in previous work (https://github.com/Tencent/MedicalNet). The resulting channel feature maps were then subjected to feature dimension reduction by filtering with the maximum value. This process yielded a set of 2048 one-dimensional features. After extracting deep features, the construction of the deep-radiomics model followed a similar procedure to the conventional radiomics model.

The parameters used for feature extraction, feature selection, internal validation, and model building are presented in the Supplementary Material.

Multiple metrics were employed to evaluate of the constructed models performance using the test dataset. First, the receiver operating curve (ROC) analysis with the use of the area under the curve (AUC) was utilized as a widely accepted metric.

Because of the imbalanced distribution of the BCR- and BCR+ cases, the precision-recall (PR) curve was employed. The AUC of the PR and 95% bootstrap confidence interval (BCa) were calculated. If the BCa does not include zero, it indicates that the difference in PR AUC between the models is statistically significant.

To gain insights into the clinical utility of the models, decision curve analysis (DCA) was employed. The net reclassification index (NRI) and integrated discrimination improvement (IDI) were utilized to quantify the enhancement in prediction accuracy achieved by the radiomics and deep-radiomics models compared to the clinical model. A positive NRI indicates improved risk classification, while a negative NRI suggests misclassification. A higher IDI value indicates enhanced discrimination between individuals with and without the event of interest.

IBM SPSS® 20.0 (www.ibm.com), MedCalc 20.014 (www.medcalc.org) and R 3.5.1 (www.r-project.org) were used for statistical analysis.

Descriptive statistics were used to summarize the data, with mean (standard deviation) reported of continuous variables that followed a normal distribution and median [Q1, Q3] of continuous variables that did not follow a normal distribution. Categorical variables were reported as frequencies (percentage%). A Mann-Whitney U test or chi-square test was used to assess differences in clinical characteristics both between the training and test cohorts and between the BCR- and BCR+ groups. The DeLong test was used to compare the difference between the AUCs of ROC. The level of statistical significance was set at P < 0.05.

Out of 131 cases, 62 (47.3%) were classified as locally advanced PCa, while 20 (15.3%) exhibited lymph node metastasis, and 49 (37.4%) presented with distant metastasis. There was no significant difference in any feature between the training and test datasets (Table 3). BCR rates were 24.7% (23/93) in the training cohort and 21.1% (8/38) in the test cohort, without significant difference (P = 0.824).

Compared the BCR+ and BCR- groups (Table 4), no significant difference was found in age, PI-RADS score, PCa location, PCa ADC value, and ISUP group (all P > 0.05). There were significant differences in PSA level, PCa volume, clinical stage and treatment method between the BCR- and BCR+ groups (all P < 0.05).

The clinical model was constructed using stepwise multivariable logistic regression. PSA and stage were the final predictor variables included in the model after selection process (Table 5, Figure 3). The conventional radiomics model incorporated PSA, stage, treatment, and four image features extracted from ADC maps as the final predictor variables (Figure 4). Similarly, the deep-radiomics model included PSA, stage, treatment, and three image features extracted using a deep learning algorithm as the final predictor variables (Figure 5).

The AUC of the ROC was 0.717 (95% CI: 0.492, 0.941) for the clinical model, 0.771 (95% CI: 0.607, 0.935) for the conventional radiomics model, and 0.954 (95% CI: 0.892, 1.000) for the deep-radiomics model in the test dataset (Figure 6A). The deep-radiomics model exhibited a significantly higher AUC than the clinical model (P = 0.033) and the conventional radiomics model (P = 0.026). However, there was no significant difference between the AUC of the clinical model and the conventional radiomics model (P = 0.570).

The AUC of the PR curve was 0.385 (95% CI: 0.191, 0.696) for the clinical model, 0.373 (95% CI: 0.209, 0.665) for the conventional radiomics model, and 0.805 (95% CI: 0.451, 0.985) for the deep-radiomics model (Figure 6B). The PR AUC difference between the clinical model and the deep-radiomics model was 0.420 (95% BCa: 0.292, 0.651), indicating a statistically significant distinction. Similarly, the PR AUC difference between the conventional radiomics model and the deep-radiomics model was 0.432 (95% BCa: 0.327, 0.638), demonstrating statistical significance. However, the PR AUC difference between the clinical model and the conventional radiomics model was 0.011 (95% BCa: -0.070, 0.216), revealing no statistical significance.

The DCA curve demonstrated superior performance of the deep-radiomics model compared to the conventional radiomics model and the clinical model across all risk thresholds. Additionally, the DCA curve of the conventional radiomics model outperformed the clinical model across all risk thresholds (Figure 6C).

Compared to the clinical model, the conventional radiomics model demonstrated an NRI of 0.359 (95% CI: 0.099, 0.973), and the deep-radiomics model showed a higher NRI of 0.508 (95% CI: 0.309, 1.175). Compared to the conventional radiomics model, the deep-radiomics model showed an NRI of 0.149 (95% CI: 0.032, 0.300).

When using the clinical model as reference, the conventional radiomics model exhibited an IDI of 0.516 (95% CI: 0.294, 0.666), and the deep-radiomics model showed a higher IDI of 0.679 (95% CI: 0.431, 0.795). When comparing the conventional radiomics model to the deep-radiomics model, the IDI was 0.164 (95% CI: -0.045, 0.363) without statistical difference.