Data from 19 patients diagnosed with localised prostate cancer and recruited as part of a Human Research Ethics Committee (HREC)-approved study (HREC/15/PMCC125) were included. Patients were scheduled for radical prostatectomy at the Peter MacCallum Cancer Center in Melbourne, Australia, and had mpMRI and PET/CT imaging prior to surgery. Patient clinical, imaging, and pathological details are given in Table 1. The Prostate Imaging-Reporting and Data System (PI-RADS) for the index lesion on mpMRI ranged from 2 to 5, with patient 59 being indeterminate due to imaging artefacts. SUVmax on PSMA PET ranged from 2.43 to 59.40. The Grade Group (GG) of index lesions and other tumour foci show that 10 patients had high-grade disease (defined as an index lesion with GG ≥ 3) and the remaining 9 patients had low-grade disease (defined as an index lesion with GG ≤ 2).

In vivo mpMRI was acquired using two 3 T scanners, the first five patients scanned with a Siemens MAGNETOM Trio (Siemens Healthcare GmbH, Erlangen, Germany) and all other patients scanned with a Siemens MAGNETOM Skyra. The imaging protocol followed guidelines from the European Society of Urogenital Radiology (ESUR) [12] and included T2-weighted (T2w), Diffusion-Weighted Image (DWI) and Dynamic Contrast-Enhanced (DCE) imaging. A surface body coil was used without an endorectal coil to reduce the chance of deformation of the prostate, and patients without contraindications were given Buscopan to reduce peristaltic motion. T2w imaging was obtained using a 2D turbo spin echo sequence using acquisition matrix = 320 × 320, FOV = 160 mm × 160 mm, slice thickness = 3 mm, TE = 89–96 ms, TR = 3500–4830 ms. DWI images were obtained using a 2D spin echo sequence with echo planar readout, with b-values = 50, 400, 800 and 1200 s/mm2, acquisition matrix = 250 × 250, FOV 250 mm × 250 mm, slice thickness = 4 mm. Apparent Diffusion Coefficient (ADC) maps were computed from DWI images using inline software.

Pre-contrast 3D T1-weighted images with variable flip angles (5°, 10°, 15°, 20°, 30°) were acquired. DCE-MRI was performed using a 3D spoiled gradient echo with a time-resolved view sharing sequence for high temporal resolution imaging (TWIST, Siemens Healthineers, Erlangen, Germany). Each patient received a 10-ml bolus injection of contrast agent Dotarem (gadoterate meglumine, Guerbet, USA), followed by a saline flush. Semi-quantitative parameters and pharmacokinetic parameters were computed using Dynamika software (Image Analysis Group, London, UK) [13]. Semi-quantitative parameters included the initial rate of enhancement (IRE), the time to peak enhancement (TTP), the maximum enhancement (ME), the time of contrast agent onset (Tonset), the time of contrast agent washout (Twashout) and the initial rate of washout (IRW). Pharmacokinetic parameters were computed using the Tofts model [14] including Ktrans (the volume transfer constant between blood plasma and extra-vascular extra-cellular space) and Ve (volume of extra-vascular extra-cellular space). The initial area under the gadolinium contrast agent concentration curve for the first 60 s post-injection (iAUGC60) was also computed.

All patients had PET/CT imaging after injection with a 68 Ga PSMA-HBED-CC (PSMA-11) tracer. Five different PET/CT scanners were used, with scanning from the base of the skull or the vertex to the upper thighs. The PET scanning bed steps were acquired starting at the upper thighs to minimise the chance of spatial shifts between PET and CT at the level of the prostate, which can be caused by patient movement, bladder filling or intestinal movements. Further acquisition details are given in Additional file 1: Table S1, including the PET image reconstruction method, Gaussian filter kernel size, the time of bed positions, tracer uptake time and PET and CT image resolution information. PET images were corrected for attenuation using the contemporaneous low-dose non-contrast CT scan and normalised by body weight to obtain PET Standardised Uptake Values (SUV) images for each patient.

After prostatectomy, the prostate specimens were embedded in agarose gel in a custom-made sectioning box for ex vivo MRI scanning, after which the specimen was cut into 5-mm sections and then microtomed at 3 µm to obtain whole-mount haematoxylin and eosin (H&E) stained sections. The most apical and basal histology sections obtained were not included in this study, as standard pathology processing required them to be cut in a parasagittal orientation which could not be easily co-registered with ex vivo MRI. Table 1 details the number of axial histology sections obtained for each patient, and the estimated prostate volume calculated using the ellipsoid method and prostate specimen measurements [11]. Each axial H&E-stained section was annotated for tumour and assigned a GG by an experienced uro-pathologist (CM) and digitised with an Epson Perfection V700 scanner to give images approximately 0.01 mm resolution. Further details of this are given in Reynolds et al. [15].

Co-registration of PET/CT with mpMRI and ground truth histology was carried out using our established framework [10, 15], which utilised ex vivo MRI to account for tissue deformation and shrinkage after prostatectomy. In brief, the PET and CT images were first qualitatively inspected to determine whether there were any spatial shifts between them due to the different timing between the scans; however, none of the PET or CT images required manual correction. Then, PET images were rigidly registered with in vivo 3D T2w MRI in 3D Slicer software [16] by using the contemporaneous CT, which provided anatomical information and higher resolution than the PET image. The computed transformation between CT and in vivo 3D T2w MRI was then applied to co-register PET with in vivo MRI. The 2D T2w, ADC maps from DWI and DCE MRI parameter maps were rigidly registered with the reference 3D T2w MRI in 3D Slicer. All co-registered in vivo MRI and PET/CT data were re-sampled to isotropic 0.8 mm voxels to match the 3D T2w MRI resolution, and deformable image registration was applied to co-register with ex vivo MRI and ground truth histology data. Figure 1 shows an example co-registered dataset with PET, in vivo T2w MRI and histology.

Co-registered imaging data for patient 26 showing a histology and T2w MRI with tumour annotated by a pathologist, and b–d axial, coronal and sagittal (respectively) T2w MRI aligned with PET showing tumour voxels defined from histology outlined in white and benign voxels defined as being 3.3 mm beyond the tumour voxel boundary outlined in blue

Correlation analysis was carried out at a voxel-wise level using the co-registered images to investigate the relationship between signals on mpMRI and PET SUV values. For this and subsequent analyses, DCE MRI parameters Tonset, Twashout and IRW were excluded as the Tonset parameter was inconsistent across the dataset and challenging to reproduce, while Twashout and IRW contained many zero value pixels as the contrast agent had not washed out from the entire prostate during the imaging timeframe. To account for an estimated 3.3 mm average registration uncertainty between histology and in vivo mpMRI and PET/CT computed in our prior study [15], benign voxels were defined as all voxels within the prostate contour which were at least 3.3 mm away from the tumour boundary (see Fig. 1). Additionally, small tumour foci with an area below that for a tumour with an average 5 mm diameter were excluded, corresponding to the upper bounds of the registration uncertainty as well as studies indicating the minimum tumour size that can be identified on mpMRI [17, 18].

Kolmogorov–Smirnov tests were performed to determine whether PET SUV and mpMRI voxel values for tumour and benign tissue were normally distributed, and whether the benign and tumour voxel values for each imaging parameter exhibited the same distribution. Spearman correlation coefficients were computed to assess the degree of correlation between the PET SUV and the mpMRI parameter tumour voxel values. Bonferroni correction was applied to the correlation cut-off p-values to assess its significance.

Voxel-wise radiomic features were extracted from co-registered mpMRI and PET imaging data using the PyRadiomics Library (v.3.0.1) [19] in the Python programming environment (v3.8). Features were extracted from the PET, T2w MRI, ADC maps from DWI and DCE MRI parameters TTP, Ktrans and iAUGC60. These DCE MRI parameters were chosen based on prior studies which indicated they were the most predictive for tumour [20,21,22]. Radiomic features included first-order statistics, shape-based and texture features based on Grey Level Co-occurrence Matrix (GLCM) (number of features, n = 24), Grey Level Run Length Matrix (GLRLM) (n = 16), Grey Level Size Zone Matrix (GLSZM) (n = 16), Neighbouring Grey Tone Difference Matrix (NGTDM) (n = 5), Grey Level Dependence Matrix (GLDM) (n = 14). Radiomic features were extracted from the original images and after applying filters including wavelet transform, gradient magnitude, Laplacian of Gaussian, 2D local binary pattern (LBP) and 3D LBP. A kernel size of 9 voxels in each direction (7.2 mm in each direction) was chosen, as it was closest to the average tumour radius in the dataset when assumed to be spherical.

Patient-based clinical features including age and PSA level, and PET-specific features were added to the model, including SUVmax, SUVmean, radioactive tracer uptake time and tracer radioactivity. Feature reduction was performed to prevent overfitting and reduce model training cost by applying the following strategies: (i) reject all highly correlated features, (ii) retain the top 10% of the features based on the ANOVA test and (iii) retain the top 50 features based on the mean decrease in random forest Gini impurity.

Machine learning classifiers were trained to predict tumour location using a radial basis function kernel with the Python scikit-learn library (v0.24.2) [23]. Two classifiers were used to compare performance, a Random Forest Classifier (RFC) and a Support Vector Classifier (SVC). A fivefold cross-validation scheme was used to optimise the parameters of each classifier via a halving successive grid search procedure. Data imbalance, caused by the number of benign voxels vastly surpassing the number of tumour voxels, was addressed through data augmentation by flipping images 180 degrees about the axial, sagittal and coronal axes. This increased the number of tumour voxels fourfold.

Classifiers were trained to predict tumour location using (1) PET images alone, (2) mpMRI data alone (T2w, ADC, TTP, Ktrans and iAUGC60) and (3) mpMRI and PET images combined. For the best performing tumour detection model, a second model was trained to predict the grade of each predicted tumour voxel, by classifying it as either high grade (defined as GG ≥ 3) or low grade (defined as GG ≤ 2). The performance of the tumour grade prediction model was assessed using balanced accuracy and weighted F1 score. Balanced accuracy was computed via the arithmetic mean of the sensitivity and specificity, measuring the combined performance of the tumour detection and tumour grading models. The weighted F1 score was calculated by taking the average F1 score from the high-grade and low-grade grading performance, which were each weighted by the proportion of voxels with the corresponding grade in the patient.