Predicting biochemical recurrence of prostate cancer with artificial intelligence.

Background: The first sign of metastatic prostate cancer after radical prostatectomy is rising PSA levels in the blood, termed biochemical recurrence. The prediction of recurrence relies mainly on the morphological assessment of prostate cancer using the Gleason grading system. However, in this system, within-grade morphological patterns and subtle histopathological features are currently omitted, leaving a significant amount of prognostic potential unexplored.
Methods: To discover additional prognostic information using artificial intelligence, we trained a deep learning system to predict biochemical recurrence from tissue in H&E-stained microarray cores directly. We developed a morphological biomarker using convolutional neural networks leveraging a nested case-control study of 685 patients and validated on an independent cohort of 204 patients. We use concept-based explainability methods to interpret the learned tissue patterns.
Results: The biomarker provides a strong correlation with biochemical recurrence in two sets (n = 182 and n = 204) from separate institutions. Concept-based explanations provided tissue patterns interpretable by pathologists. Conclusions: These results show that the model finds predictive power in the tissue beyond the morphological ISUP grading.

Introduction

Prostate cancer is a common malignancy among men, affecting 1.4 million per year[1]. A significant proportion of these men will receive the primary curative treatment of a prostatectomy. This surgery’s success can partly be judged by the concentration of prostate-specific antigen (PSA) in the blood. While it has a dubious role in prostate cancer screening[2,3], this protein is a valuable biomarker in PCa patients’ follow-up post-prostatectomy. In a successful surgery, the concentration will mostly be undetectable (<0.1 ng/mL) after 4–6 weeks[4].

However, in ~30% of the patients[5-7], PSA will rise again after surgery, called biochemical recurrence, pointing to regrowth of prostate cancer cells. Biochemical recurrence is a prognostic indicator for subsequent progression to clinical metastases and prostate cancer death[8]. Estimating chances of biochemical recurrence could help to better stratify patients for specific adjuvant treatments.

The risk of biochemical recurrence of prostate cancer is currently assessed in clinical practice through a combination of the ISUP grade[9], the PSA value at diagnosis and the TNM staging criteria. In a recent European consensus guideline, these factors were proposed to separate the patients into a low-risk, intermediate-risk and high-risk group[10]. A high ISUP grade independently can, independently of other factors, assign a patient to the intermediate (grade 2/3) or high-risk group (grade 4/5).

Based on the distribution of the Gleason growth patterns[11], which are prognostically predictive morphological patterns of prostate cancer, pathologists assign cancerous tissue obtained via biopsy or prostatectomy into one of five groups. They are commonly referred to as International Society of Urological Pathology (ISUP) grade groups, the ISUP grade, Gleason grade groups, or just grade groups.[9,12-14]. Throughout this paper, we will use the term ISUP grade. The ISUP grade suffers from several well-known limitations. For example, there is substantial disagreement in the grading using the Gleason scheme[15,14]. Furthermore, although the Gleason growth patterns have seen significant updates and additions since their inception in the 1960s, they remain relatively coarse descriptors of tissue morphology. As such, the prognostic potential of more fine-grained morphological features has been underexplored. We hypothesize that artificial intelligence, and more specifically deep learning, has the potential to discover such information and unlock the true prognostic value of morphological assessment of cancer. Specifically, we developed a deep learning system (DLS), trained on H&E-stained histopathological tissue sections, yielding a score for the likelihood of early biochemical recurrence.

Deep learning is a recent new class of machine learning algorithms that encompasses models called neural networks. These networks are optimized using training data; images with labels, such as recurrence information. From the training data, relevant features to predict the labels are automatically inferred. During development, the generalization of these features is tested on separated training data, which is not used for learning. Afterwards, a third independent set of data, the test set, is used to ensure generalization. Since features are inferred, handcrafted feature engineering is not needed anymore to develop machine learning models. Neural networks are the current state-of-the-art in image classification[16].

Deep learning has previously been shown to find visual patterns to predict genetic mutations from morphology, for example, in lymphoma[17] and lung cancer[18]. Additionally, deep learning has been used for feature discovery in colorectal cancer[19] and intrahepatic cholangiocarcinoma[20] using survival data. Although deep learning has been used with biochemical recurrence data on prostate cancer, Leo et al.[21]. assumed manual feature selection beforehand, strongly limiting the extent of new features to be discovered. Yamamoto et al.[22]. used whole slide images and a deep-learning-based encoding of the slides to tackle the slides’ high resolution. They leverage classical regression techniques and support-vector machine models on these encodings. The deep learning model was not directly trained on the outcome, limiting the feature discovery in this work as well.

A common critique of deep learning is its black-box nature of the inferred features[23]. Especially in the medical field, decisions based on these algorithms should be extensively validated and be explainable. Besides making the algorithms’ prediction trustworthy and transparent, from a research perspective, it would be beneficial to visualize the data patterns which the model learned, allowing insight into the inferred features. We can visualize the patterns learned by the network leveraging a new technique called Automatic Concept Explanations (ACE)[24]. ACE clusters patches of the input image using their intermediate inferred features showing common patterns inferred by the network. We were interested in finding these common concepts over a range of images to unravel patterns that the model has identified.

This study aimed to use deep learning to develop a new prognostic biomarker based on tissue morphology for recurrence in patients with prostate cancer treated by radical prostatectomy. As training data, we used a nested case-control study[25]. This study design ensured we could evaluate whether the network learned differentiating patterns independent of Gleason patterns. The prognostic biomarker provides a strong correlation with biochemical recurrence in two sets (n = 182 and n = 204) from separate institutions. Furthermore, the Automatic Concept-based Explanations provided tissue patterns interpretable by our pathologist.

Methods

Cohorts

Two independent cohorts of patients who underwent prostatectomy for clinically localized prostate cancer were used in this study. Patients were treated at either the Johns Hopkins Hospital in Baltimore or New York Langone Medical Centre. Both cohorts were accessed via the Prostate Cancer Biorepository Network[26]. The Johns Hopkins University School of Medicine Institutional Review Board and The New York University School of Medicine Institutional Review Board provided ethical regulatory approval for collection and disbursement of data and materials from the respective institutions. The need for acquiring informed consent was waived by the institutional ethical review boards.

For the development of the novel deep-learning-based biomarker (further referred to as DLS biomarker), we used a nested case-control study of patients from Johns Hopkins. This study consists of 524 matched pairs (724 unique patients) containing four tissue spots per patient. They were sampled from 4860 prostate cancer patients with clinically localized prostate cancer who received radical retropubic prostatectomy between 1993 and 2001. Men were routinely checked after prostatectomy at 3 months and at least yearly thereafter. Surveillance for recurrence was conducted using digital rectal examination and measurement of serum PSA concentration. Patients were followed for outcome until 2005, with a median follow-up of 4.0 years. The outcome was defined as recurrence, based on biochemical recurrence (serum PSA > 0.2 ng/mL on 2 or more occasions after a previously undetectable level after prostatectomy), or events indicating biochemical recurrence before this was measured; local recurrence, systemic metastases, or death from prostate cancer. Controls were paired to cases with recurrence using incidence density sampling[27]. For each case, a control was selected who had not experienced recurrence by the date of the case’s recurrence and was additionally matched based on age at surgery, race, pathologic stage, and Gleason sum in the prostatectomy specimen based on the pathology reports. Given the incidence density sampling of controls, some men were used as controls for multiple cases, and some controls developed recurrence later and became cases for that time period.

The TMA spots were cores (0.6 mm in diameter) from the highest-grade tumour nodule. Random subsamples were taken in quadruplicate for each case. The whole slides were scanned using a Hamamatsu NanoZoomer-XR slide scanner at 0.23 μm/px. TMA core images were extracted using QuPath (v0.2.3[28],). We discarded analysis of cores with <25% tissue. The cores were manually checked (HP) for prostate cancer, excluding 535 cores without clear cancer cells present in the TMA cross-section, resulting in a total of 2343 TMA spots. The nested case-control set was split based on the matched pairs into a development set (268 unique pairs), and a test set (91 pairs); the latter was used for evaluation only. We leveraged cross-validation by subdividing the development into three folds to tune the models on different parts of the development set. We divided paired patient, randomly, keeping into account the distribution of the matched variables. The random assignment was done using the scikit-multilearn package[29], specifically the ‘IterativeStratification’ method in ‘skmultilearn.model_selection’. After splitting the dataset into training and test, we split the training dataset into three folds using the same method for the cross-validation.

To validate the DLS biomarker on a fully independent external set, we used the cohort from New York Langone Medical Centre. This external validation cohort consists of 204 patients with localized prostate cancer treated with radical prostatectomy between 2001 and 2003. Patients were followed for outcome until 2019, with a median follow-up of 5 years. Biochemical recurrence was defined as either a single PSA measurement of ≥0.4 ng/m or PSA level of ≥0.2 ng/ml followed by increasing PSA values in subsequent follow-up. Cores were sampled from the largest tumour focus or any higher-grade focus (>3 mm). Subsamples were taken in quadruplicate for each case. Images were scanned using a Leica Aperio AT2 slide scanner at 0.25 μ/px.

Model details

For developing the convolutional neural networks (CNNs) we used PyTorch[30]. As an architecture, we used ResNet50-D[31] pretrained on ImageNet from PyTorch Image Models[32]. We used the Lookahead optimizer[33] with RAdam[34], with a learning rate of 2e-4 and mini-batch size of 16 images. We used weight decay (7e-3), and a drop-out layer (p = 0.15) before the final fully-connected layer. We used EfficientNet-style[35] dropping of residual connections (p = 0.3) as implemented in PyTorch Image Models. We used Bayesian Optimization to find the optimal values (See Supplementary Notes 1 for details about the searchspace).

We resized the TMAs to 1.0 mu/pixel spacing and cropped to 768 × 768 pixels. Extensive data augmentations were used to promote generalization. The transformations were: flipping, rotations, warping, random crop, HSV colour augmentations, jpeg compression, elastic transformations, Gaussian blurring, contrast alterations, gamma alterations, brightness alterations, embossing, sharpening, Gaussian noise and cutout[36]. Augmentations were implemented using albumentations[37] and fast.ai[38].

TMA spots from cases experiencing recurrence were assigned a value of 0–4, depending on the year on which the first event, either biochemical recurrence, metastases, or prostate cancer-related death, was recorded, with 0 meaning recurrence within a year, four meaning after 4+ years. TMA spots from cases without an event were also assigned the label 4.

We validated the model on the development validation fold each epoch with a moving average of the weights from five subsequent epochs. We used the concordance index as a metric to decide which model performed the best.

As the final prediction at the patient level, the TMA spot with the highest score was used. The final DLS consists of an ensemble of 15 convolutional neural networks. Using cross-validation as described above, 15 networks were trained for each fold, of which the five best performing were used for the DLS. See Fig. 1 for a graphical overview of the methods, further details can be found in the Supplementary Methods.

Fig. 1

Overview of the methods summarizing the biomarker development and the Automatic Concept Explanations (ACE) process.

Cores were extracted from TMA slides and used to train a neural network to predict the years to biochemical recurrence. On the nested case-control test set, a matched analysis was performed. For ACE, patches were generated from the cores, inferenced through the network and clustered based on their intermediate features.

Statistics and reproducibility

For primary analysis of the nested case-control study, odds ratios (OR) and 95% confidence intervals (CI) were calculated using conditional logistic regression, following Dluzniewski et al.[39]. Due to the study design, calculating hazard ratios using a Cox proportional hazard regression is not appropriate. For the primary analysis, the continuous DLS marker was given as the only variable. For a secondary analysis, we added the non-matched variables PSA, positive surgical margins, and a binned indicator variable for year of surgery. Since matching was done on Gleason sum, and our goal was to identify patterns beyond currently used Gleason patterns, we corrected for the residual differences of the ISUP grade between cases and control (see Table 1). A correction was performed by adding a continuous covariate since, due to the small differences, an indicator covariate did not converge. Analysis was done using the lifelines Python package (v. 0.25.10)[40] with Python (v. 3.7.8). P-values were calculated as a Wald test per single parameter. Since the DLS predicts the time-to-recurrence, high values indicate a low probability of recurrence. We multiplied the DLS output by −1 to make the analysis more interpretable. For three patients (1 from the Johns Hopkins cohort and 2 from the New York Langone cohort), PSA values were missing and were therefore replaced by the median.

Table 1

Baseline characteristics of test set and development set from the John Hopkins Hospital, prostate cancer recurrence cases and controls, men who underwent radical prostatectomy for clinically localized disease between 1993 and 2001.

	Development set			Test set
	Recurrence cases	No events cases	P	Recurrence cases	Controls*	P
N	368	135		91	91
Age, mean (SD)	58.9 (6.2)	59.3 (6.3)	p = 0.540	58.4 (6.1)	58.3 (6.3)	Matched
preop. PSA (ng/mL), mean (SD)	12.3 (10.0)	10.1 (7.5)	p = 0.010	12.3 (10.8)	10.5 (7.7)	p = 0.195
Race, n (%)			p = 0.599			Matched
White	327 (88.9)	120 (88.9)		72 (79.1)	75 (82.4)
Black or African–American	32 (8.7)	14 (10.4)		12 (13.2)	10 (11.0)
Other	9 (2.4)	1 (0.7)		7 (7.7)	6 (6.6)
Pathological stage			p = 0.107			Matched
pT2	43 (11.7)	25 (18.5)		20 (22.0)	19 (20.9)
pT3a	199 (54.1)	63 (46.7)		50 (54.9)	51 (56.0)
pT3b or N1	126 (34.2)	47 (34.8)		21 (23.1)	21 (23.1)
Gleason sum prostatectomy (%)			p = 0.179			Matched
6	38 (10.3)	25 (18.5)		20 (22.0)	23 (25.3)
7	233 (63.3)	76 (56.3)		51 (56.0)	50 (54.9)
8+	97 (26.4)	34 (25.2)		20 (22.0)	18 (19.8)
ISUP grade, n (%)			p = 0.002			p = 0.851
1	38 (10.3)	25 (18.5)		20 (22.0)	23 (25.3)
2	140 (38.0)	61 (45.2)		35 (38.5)	38 (41.8)
3	93 (25.3)	15 (11.1)		16 (17.6)	12 (13.2)
4	49 (13.3)	21 (15.6)		13 (14.3)	10 (11.0)
5	48 (13.0)	13 (9.6)		7 (7.7)	8 (8.8)
Positive surgical margins	140 (38.1)	24 (17.8)	p < 0.001	36 (39.6)	20 (22.0)	p = 0.016
Mean year of surgery	1997.0 (2.3)	1995.5 (2.3)	p < 0.001	1997 (2.3)	1995 (2.1)	p < 0.001

*Due to the nested case-control nature, some controls could have a biochemical recurrence, but always later than their matched case.

For primary analysis of the New York Langone cohort, we calculated hazard ratios (HR) using a Cox proportional hazards regression. We report a secondary multivariable analysis including indicator variables for relevant clinical covariates, Gleason sum, pathological stage, and surgical margin status. We tested the proportional hazards assumption as satisfactory (every p-value > 0.01) using the Pearson correlation between the residuals and the rank of follow-up time. Kaplan–Meier plots were generated for the New York Langone cohort. Due to the nested case-control design for the Johns Hopkins set, this set could not be visualized in a Kaplan–Meier plot.

Automatic concept explanations

To generate concepts, we picked the best performing single CNN from the DLS based on its validation set fold. We used a combination of the methods of Yeh et al., 2020[41] and Ghorbani et al., 2019[24].

We tiled the TMA images into 256 × 256 patches within the tissue, discarding patches with >50% whitespace. These patches were padded to the original input shape of the CNN (768 × 768 pixels). The latent space of layer 42 of 50 was saved for each tile. Afterwards, we used PCA (50 components) to lower the dimensionality and then performed k-means (k = 15) to cluster the latent spaces.

In contrast to Yeh et al. and Ghorbani et al., we did not sort the concepts on completeness of the explanations or importance for prediction of individual samples. We sorted the concepts to find interesting new patterns related to recurrence across images by ranking the concepts based on the DLS score of the TMA spot from which they originated.

For each concept, 25 examples were randomly picked and visually inspected by a pathologist (JvI), with a special interest in uropathology, blinded to the case characteristics and prediction of the network.

Table 2

Baseline characteristics of the cohort from New York Langone hospital, prostate cancer recurrence cases and controls, men who underwent radical prostatectomy between 2001 and 2003.

	Recurrence cases	Controls	P
N	38	166
preop. PSA (ng/mL), mean (SD)	11.6 (11.5)	6.7 (3.9)	p = 0.014
Age, mean (SD)	61.7 (8.9)	60.3 (6.6)	p = 0.359
Race, n (%)			p = 0.401
African–American	2 (5.3)	4 (2.4)
Asian	2 (5.3)	3 (1.8)
Caucasian	33 (86.8)	144 (86.7)
Not reported	0 (0)	2 (1.2)
Other	1 (2.6)	13 (7.8)
Pathological stage, n (%)			p < 0.001
pT2a	0 (0)	12 (7.2)
pT2b	3 (7.9)	5 (3.0)
pT2c	16 (42.1)	114 (68.7)
pT3a	10 (26.3)	27 (16.3)
pT3b	9 (23.7)	8 (4.8)
ISUP grade, n (%)			p < 0.001
1	3 (7.9)	67 (40.4)
2	13 (34.2)	76 (45.8)
3	6 (15.8)	13 (7.8)
4	5 (13.2)	3 (1.8)
5	11 (28.9)	7 (4.2)
Surgical Margins, n (%)			p = 0.060
Focal	10 (26.3)	20 (12.0)
Free of tumour	27 (71.1)	144 (86.7)
Widespread	1 (2.6)	2 (1.2)

Table 3

Conditional logistic regression analyses of the Johns Hopkins test set.

Covariate	Matched analysisJohns Hopkins (OR)a	Multivariate analysisJohns Hopkins (OR)
Biomarker	3.28 (CI 1.73–6.23; p < 0.005)	3.32 (CI 1.63–6.77; p = 0.001)
preop. PSA (ng/mL)		1.04 (CI 0.99–1.10; p = 0.10)
Surgical margins (pos)		1.69 (CI 0.69–4.18; p = 0.25)
ISUP grade (cont.)b		1.34 (CI 0.64–2.82; p = 0.44)
Mean year of surgery
1992–1994 (n = 75)		1.0
1994–1997 (n = 55)		3.35 (CI 1.13–9.91; p = 0.03)
1997–2001 (n = 52)		8.22 (CI 2.38–28.37; p = 0.0009)

aCases and controls were matched on age at surgery, race, pathologic stage, and Gleason sum in the prostatectomy specimen.

b The ISUP grade covariate was added to correct for the residual differences left after matching cases with controls on prostatectomy Gleason sum.

Table 4

Cox proportional hazard analyses of New York Langone external validation cohort.

Covariate	Univariate analysisNYU (HR)	Multivariate analysisNYU (HR)
Biomarker	4.79 (CI 2.09–10.96; p = 0.0002)	3.02 (CI 1.10–8.29; p = 0.03)
preop. PSA (ng/mL)		1.07 (CI 1.02–1.12; p = 0.004)
ISUP grade
1		1.0
2		2.64 (CI 0.73–9.58; p = 0.14)
3		8.74 (CI 2.16–35.30; p = 0.00)
4		12.78 (CI 2.82–57.91; p = 0.00)
5		9.60 (CI 2.32–39.69; p = 0.00)
Pathological stage
pT2a + b		1.0
pT2c		1.02 (CI 0.27–3.80; p = 0.98)
pT3a		1.26 (CI 0.28–5.67; p = 0.77)
pT3b		2.77 (CI 0.66–11.62; p = 0.16)
Surgical margins
Free		1.0
Focal		2.13 (CI 0.76–5.96; p = 0.15)
Widespread		0.20 (CI 0.01–3.39; p = 0.27)