Automatic Detection of Prostate Tumor Habitats using Diffusion MRI

A procedure for identification of optimal Apparent Diffusion Coefficient (ADC) thresholds for automatic delineation of prostatic lesions with restricted diffusion at differing risk for cancer was developed. The relationship between the size of the identified Volumes of Interest (VOIs) and Gleason Score (GS) was evaluated. Patients with multiparametric (mp)MRI, acquired prior to radical prostatectomy (RP) (n = 18), mpMRI-ultrasound fused (MRI-US) (n = 21) or template biopsies (n = 139) were analyzed. A search algorithm, spanning ADC thresholds in 50 µm2/s increments, determined VOIs that were matched to RP tumor nodules. Three ADC thresholds for both peripheral zone (PZ) and transition zone (TZ) were identified for estimation of VOIs at low, intermediate, and high risk of prostate cancer. The determined ADC thresholds for low, intermediate and high risk in PZ/TZ were: 900/800; 1100/850; and 1300/1050 µm2/s. The correlation coefficients between the size of the high/intermediate/low risk VOIs and GS in the three cohorts were 0.771/0.778/0.369, 0.561/0.457/0.355 and 0.423/0.441/0.36 (p < 0.05). Low risk VOIs mapped all RP lesions; area under the curve (AUC) for intermediate risk VOIs to discriminate GS6 vs GS ≥ 7 was 0.852; for high risk VOIs to discriminate GS6,7 vs GS ≥ 8 was 0.952. In conclusion, the automatically delineated volumes in the prostate with restricted diffusion were found to strongly correlate with cancer aggressiveness.

The use of multiparametric MRI (mpMRI) is rapidly gaining momentum in the management of prostate cancer because of its improved diagnostic potential. Diffusion Weighted Imaging (DWI) and the associated Apparent Diffusion Coefficient (ADC) maps play a central role in the Prostate Imaging Reporting and Data System (PI-RADS) 1,2 , as well as in computer-aided diagnosis (CAD) systems and in the search for quantitative imaging biomarkers for tumor assessment [3][4][5][6][7] . Mean ADC values or other parameters of the ADC distribution within suspect cancer areas are typically investigated in relation to cancer aggressiveness. However, the size of the prostatic 3D regions of decreased diffusion have not been considered, partly due to the fact that setting a threshold below which to determine these volumes of "low" ADC, a.k.a Volumes of Interest (VOIs), is not trivial. In previous work, ADC = 800 µm 2 /s and 850 µm 2 /s were used as such thresholds in the Peripheral (PZ) and Transition Zone (TZ), correspondingly 8 . These ad-hoc values were selected based on literature and empirical observations [9][10][11][12][13][14] . In a subsequent pilot study of a retrospective patient cohort (n = 139), the size of the VOIs, defined by pixels with ADC < 800/850 µm 2 /s (PZ/TZ) were significantly correlated with biopsy Gleason Score (GS) 15 .
These preliminary results motivated the current study for systematic determination of optimal ADC thresholds, associated with high, intermediate and low risk of cancer. The underlying hypothesis is that the size of the VOIs with reduced diffusion are related to tumor aggressiveness. The purposes of this study are: (i) to develop a procedure for identification of optimal ADC thresholds for automatic delineation of prostatic lesions at differing risk for cancer; and (ii) to evaluate the relationship between the size of these VOIs and GS. The ADC thresholds were determined in prostatectomy samples, using histopathology location and aggressiveness of the tumor nodules as "ground truth". The size of the VOIs, defined by these thresholds, were evaluated for correlation with GS in independent datasets.

Results
Patients. An Institutional Review Board (IRB) at the University of Miami approved a protocol for retrospective review of prostate mpMRI. The IRB waived the need for informed consent. All research was performed in accordance with the relevant guidelines/regulations associated with the protocol. Imaging datasets from three cohorts of patients were analyzed: (i) eighteen patients, who underwent radical prostatectomy (RP); (ii) 21 patients, who underwent MRI/Ultrasound fused (MRI-US) biopsy; and (iii) 139 patients with template transrectal ultrasound-guided (TRUS) biopsies. The inclusion criteria was that mpMRI was acquired on the 3 T Discovery MR750 magnet (GE, Waukesha, WI). Here and throughout the manuscript, these datasets are referred as prostatectomy, target biopsies, and template biopsy datasets. The clinical characteristics of the patients are presented in Table 1. mpMRI were acquired between May 2012 and November 2016. Surgery/biopsies were carried out (mean ± SD) 74 ± 45/48 ± 56 days after imaging.
Prostatectomy dataset. Mapping tumor lesions from RP on MRI. RP specimens were handled in accordance with well-established procedures 16 . All specimens were reviewed by an urologic pathologist (ONK) who was blinded to the imaging findings. Cancer was mapped on hematoxylin and eosin (H&E) stained slides (histopathology ROIs). Each tumor nodule was individually staged and graded according to the latest recommendations 17,18 . The annotated slides were scanned and prostate quadrants were "stitched" into a pseudo whole-mount RP sample (Figs 1a and 2a). The locations of the histopathology Regions of Interest (ROIs) from radical prostatectomy were mapped onto MRI (rpROI) in MIM (Cleveland, OH). Out of a total of 53 H&E nodules (Table 1), three regions could not be mapped on MRI. The rpROI volumes were measured in MIM and summarized in Supplementary Table S1, together with their location (PZ/TZ) and GS. In Fig. 3, the box plot of the rpROI volumes versus GS is presented (ρ = 0.674, p < 0.001).
ADC thresholds for identifying prostate volumes at differing risk for cancer. To identify and assign a level of risk of prostate cancer to various regions within the prostate, an automated software based on ADC was developed. The algorithm searched through ADC values to determine VOIs that matched the location of rpROIs. The prostate was divided in 9 quadrants (Fig. 4). The algorithm then determined ADC thresholds that maximized the sum of r, the Pearson Correlation coefficients between the nine pairs of fractions of rpROIs and ADC subvolumes in each quadrant, and ρ, Spearman's ρ between VOIs and GSs. Three thresholds: D HR , D IR, D LR (D HR < D IR < D LR ) were determined for both PZ and TZ, representing cut-points defining volumes related to areas at high (VOI HR ),      Fig. 1b. Based on the histopathology review, there is large inter-and intra-lesional heterogeneity. For instance, the anterior lesions in Fig. 1a vary in GS across the slices (GS = 6 in the apex and GS = 7 in mid and base). Note the excellent correlation between the red color intensities (habitats at high-risk) in Fig. 1b with the higher microscopic tumor grade. There are also areas of low (yellow color) and intermediate (blue color) risk in these lesions. VOI HR , VOI IR , and VOI LR depict also the inter-and intra-lesional heterogeneity in the example in Fig. 2b. The regions in yellow depict areas of low risk and in general appear to be larger than the tumor on RP. This is illustrated on Fig. 5, where the relationship between rpROIs and VOI HR / VOI IR /VOI LR are shown (Pearson correlation coefficient 0.86, 0.93 and 0.91, respectively). While all volumes were significantly correlated with the rpROIs (p < 0.0001), the high and intermediate risk volumes are smaller than the rpROIs. This is consistent with the notion that VOI HR and VOI IR are related to the more aggressive tumor habitats. The low-risk volumes, by design, are supposed to map all tumor areas, including GS6. These volumes appear to be larger by a factor of 1.4, that is consistent with reported RP tissue-shrinkage factor (1.22-1.5) 19 .
Biopsy datasets. Twenty-eight target lesions were identified in the targeted biopsy dataset and 139 biopsies were analyzed in the template biopsy datasets (Table 1). Biopsy GS and PI-RADS were recorded for the locations of biopsies in the targeted biopsy dataset (tbROI). VOI HR , VOI IR and VOI LR were determined using D HR , D IR and D LR cut-points. Subvolumes, containing the tbROIs, were paired with the corresponding GS. For the patents in the template biopsy dataset only the highest GS from the biopsy session (6 to 18 biopsies) was available. The largest subvolume of VOI HR , VOI IR and VOI LR was matched with the patients's GS.
In Fig. 6 the associations of the sizes of VOI HR and VOI IR in all three datasets with GS are shown. All correlations, including VOI LR (data not shown) were statistically significant.
As expected, the highest correlations were in the prostatectomy dataset, as the tumors were most comprehensively annotated on histological slides and the thresholds were identified to maximize exactly these correlations. The correlations were decreasing with the reduction in precision of localization of the tumor lesions (from targeted to template biopsy datasets). The correlations remained very similar when the relative tumor volume (ratio between VOI HR /VOI IR /VOI LR and prostate volume) was used (Fig. 7).
Sensitivity, specificity, negative and positive predictive values and AUC for the VOIs for discriminating (i) GS6 versus 7+; and (ii) GS6,7 versus 8+ lesions are presented in Table 3 as well as the results for PI-RADS in the targeted biopsy dataset. AUC curves are shown in Fig. 8.   Figure 4. Prostate quadrants. The prostate and the Peripheral Zone (PZ) were manually contoured on the T2weighted MRI in MIM. Transition Zone (TZ) contour was estimated algebraically by subtracting PZ from the prostate. The prostate is automatically segmented in three sections: apex, mid and base. 3 segments: PZ and two for TZ are automatically generated for each section, resulting in a 9-element 3D grid. Representative T2weighted axial images in apex, mid and base with the corresponding grid are shown.

Discussion
The use of mpMRI is rapidly gaining momentum in the management of prostate cancer. Currently, PI-RADS are the standard of care for identification of areas for targeted biopsy. The five-score system does not tap into the wealth of quantitative imaging information contained in the multiple sequences of mpMRI, nor does it elucidate intralesional spatial heterogeneity. In contrast, the proposed approach, based only on diffusion MRI, maps prostate tumor heterogeneity by assigning each pixel at three levels of risk for prostate cancer. This is the first study which identified ADC thresholds for volumes of reduced diffusion and associated the automatically identified volumes with cancer grade. The size of the automatically delineated volumes in the prostate with restricted diffusion strongly correlate with cancer aggressiveness. The simplicity of the approach and strength of the determined relationships are quite revealing. The association of ADC volumes and GS (Figs 6 and 7), analogous to the correlation between rpROIs and GS (Fig. 3), is quite novel.
In comparison, when compared to high and intermediate ADC volumes in discriminating between GS6,7 and GS8+ lesions, PI-RADS assessments for sensitivity, specificity, NTV and PPV were lower (Table 3). In addition, PI-RADS contours do not depict the 3D tumor volume. While this may not be an objective of the radiologists/ PI-RADS, areas of most aggressive disease need to be consistently and objectively identified in 3D for assigning biopsy targets or dose-boost areas in radiotherapy of the prostate.
This study has several limitations. First, this is a retrospective study with all inherent limitations of the design. Second, ADC thresholds were derived using a limited set of 18 RPs. The small sample size may present a problem for the generalizability of the identified thresholds, especially since the clinical characteristics in Table 1 indicate differences between the three datasets. For instance, based on PSA and T-stage, the patients in the prostatectomy and target biopsy datasets appear to be at higher risk relative to the patients in the template biopsy datasets. Based on Gleason Score, however, the patients in the prostatectomy dataset are of lower risk: roughly 50% of the mapped 50 lesions were GS6. Gleason Score is the strongest, most reproducible predictor of clinical outcome and in this sense can be concluded that the prostatectomy patients are a lower risk cohort. Thus, the ADC thresholds, determined on this dataset warrant high sensitivity in detection of prostate cancer lesions. Third, only patients with proven cancer were included in the study as the diagnostic potential of the approach (i.e. cancer vs no cancer) was not a study objective. The application of the thresholds also requires contours of PZ and TZ. The use of a prostate atlas will decrease significantly the need of manual contouring 20 . Finally, the generalizability of the determined thresholds to other MRI sequences, vendors, magnetic field strengths and coils (endorectal vs body) should be investigated.
The ADC maps allow for: (i) better definition of biopsy targets for the identification of high grade disease and potentially the acquisition of tissue more representative of aggresiveness; and (ii) better definition of tumor volumes for any focal type of therapy, such as focused radiotherapy dose escalation. The identified ADC thresholds were recently integrated in Habitat Risk Score (HRS), a pixel by pixel ten-scale system that combines the quantitative Dynamic Contrast-Enhanced (DCE-)MRI and ADC 21 . HRS maps were used prospectively to guide radiotherapy (RT) boost volumes in a randomized Phase II clinical trial, comparing two methods of increasing dose to the mpMRI-defined tumor habitat region(s). In conclusion, ADC thresholds for determining volumes at risk in the prostate are presented. The values can be utilized to generate 3D volumes or serve as guide for prostate evaluation. Mapping tumor lesions from radical prostatectomy on MRI. RP specimens were handled in accordance with a well-established procedure 16 . The radical prostatectomy (RP) specimens were first fixed overnight in ambient formalin without injection. The seminal vesicles were amputated and the prostate weight was recorded without seminal vesicles 16 . The prostates were inked for microscopic margin assessment and submitted entirely for histologic examination in regular size cassettes. Apex and base were cut at 5-7 mm into the prostate and submitted as perpendicular sections to the margin. The rest of the prostate was cut at 3-4 mm intervals perpendicular to the urethra and submitted as quadrants. All specimens were reviewed by a urologic pathologist (ONK) who was blinded to the imaging findings. Cancer was mapped on hematoxylin and eosin (H&E) stained slides (histopathology ROIs). Tumor nodules were considered spatially separate if they were at least 3 mm apart in a plane of section or at least 4 mm on consecutive sections 22 . Each tumor nodule was individually staged and graded according to the latest recommendations 17,18 . The prostate was then cut in quadrants, and tumor nodules contoured and graded by a urologic pathologist (ONK) 17,18 . The annotated slides were scanned and prostate quadrants were "stitched" into a pseudo whole-mount RP sample (Figs 1a and 2a). Patient's T2w and ADC were uploaded to a commercial image software package (MIM, Cleveland, OH). MIM is an image processing and analysis platform incorporating DICOM input/output, and includes fusion and contouring tools as well as an interface for user-defined algorithms. The goal was to 'transfer' the locations of the histopathology Regions of Interest (ROIs) from radical prostatectomy (rpROI) onto MRI. Correlation of MRI images with histopathologic whole-mount sections is challenging because of prostate deformation following prostatectomy and during fixation and processing 23 . The issues of one-to-one correspondence of MR images and histopathologic whole mount sections are well understood and acknowledged 23,24 . Following the approach, described in Wibmer et al. 5 and using anatomical landmarks (urethra, ejaculatory ducts, hyperplastic nodules), freehand ROIs were drawn on the ADC maps and denoted rpROI. ROIs were drawn on all slices encompassing the lesion. Pathologist's contours were then transferred on MRI using rpROIs 5,21 . The rpROIs volumes were estimated in MIM.

Methods
Determination of ADC thresholds for identifying prostate volumes at differing risk for cancer. An automated algorithm for identification of optimal ADC thresholds for automatic delineation of prostatic lesions at differing risk for cancer was developed. Because of the differences in the signal intensities in PZ and TZ (Fig. 9), the algorithm determines zone-specific ADC thresholds. The gist of the method is iterating through ADC values (between 600 and 1300 µm 2 at a step of 50 µm 2 ); for each tested ADC threshold, the volumes, containing pixels with lower ADCs are evaluated in terms of overlap with rpROIs. The algorithm can be described as follows: Let N be the number of patients in the dataset. If there are I m radical prostatectomy (RP) lesions for the m th patient (m = 1, …, N), then rpROI i m denotes the i th lesion (i = 1, …, I m ) with Gleason Score GS i m . Three thresholds (D HR , D IR and D LR , where D HR < D IR < D LR ) are sought to identify the location and volumes at high (VOI HR ), intermediate (VOI IR ) and low risk (VOI LR ) areas for cancer. The thresholds in ADC intensities are set separately for the peripheral zone (PZ) and transition zone (TZ). The routine developed in Java initially identifies D HR by spanning ADC values between 700-1400 µm 2 /s in PZ and 600-1100 µm 2 /s in TZ by increments of 50 µm 2 /s. Pixels with intensities under 400 µm 2 /s were considered artefacts and not taken in account.
Let t be a pair of ADC thresholds in PZ and TZ; t = [t PZ , t TZ ]. The volume V t m in the m th patient for the threshold t consists of voxels in PZ with ADC values < t PZ and pixels in TZ with ADC values < t TZ . V t m is composed of multiple disconnected volumes. For consistency, volumes less than β cc are removed using the "cleaning" utility in MIM. Assuming that there are J m continuous volumes for the threshold t in the m th patient, let V tj m be the j th volume (j = 1, …, J m ). Each volume and rpROI were matched spatially using the 9 quadrants technique (Fig. 9).    Once thresholds D HR , D IR, D LR (D HR < D IR < D LR ) were identified, areas at high (VOI HR ), intermediate (VOI IR ) and low risk (VOI LR ) for cancer could be identified. The procedure was repeated for all patients in the prostatectomy dataset, resulting in VOIs that were paired with the GS of the corresponding rpROIs, and Spearman's ρ between VOIs and GS was calculated.
Biopsy datasets. The determined D HR , D IR and D LR were evaluated in the target and template biopsy datasets.
Both cohorts of patients were identified from the Radiation Oncology Prostate database.
Patients in the target biopsy dataset participated in a contemporary radiotherapy clinical trial and underwent MRI-US biopsies as per protocol. MRI-US biopsies were carried out on Uronav (Invivo corp., Gainesville, FL). Prostate and suspicions lesions were contoured by an experienced fellowship-trained radiologist (FM) using  were co-registered using deformable fusion. Biopsy targets on MRI were transferred onto real-time ultrasound. According to the clinical trial protocol, only the designated target areas, usually 1 or 2, were biopsied. A maximum of three biopsy cores were acquired per target. Limited number of cores are acquired in order to decrease the risk for bleeding, as the patients are planned to start radiotherapy within a month from the biopsy session. The biopsy tissue was reviewed by an urologic pathologist (ONK). The locations of biopsies with GS ≥ 6 (tbROI) were transferred to MRI in MIM, using the biopsy needle tracks in Uronav. Biopsy GS and PI-RADS were recorded for each tbROI. VOI HR , VOI IR and VOI LR were determined using D HR , D IR and D LR cut-points. Subvolumes, containing the tbROIs, were paired with the corresponding GS. As 1 to 3 biopsies were acquired per target, the highest GS is used for ADC volume correlations. VOI HR , VOI IR , and VOI LR are calculated as sets of pixels with ADC < D HR , D IR, and D LR, respectively. VOI HR , VOI IR , and VOI LR consist typically of multiple discontinuous volumes across the prostate. Small volumes (<0.1 cc) are "cleaned" and the remaining continuous volumes are sorted by size and labeled. For instance, let VOI HR consist of three subvolumes. These volumes can be visualized in MIM and the one containing the biopsy needle is selected to be paired with the biopsy GS.
For the patents in the template biopsy dataset only the highest GS from the biopsy session (6 to 18 biopsies) was available. The largest subvolume of VOI HR , VOI IR and VOI LR was matched with the patients's GS.

Statistical analysis.
For each of the three cohorts, Spearman's ρ was calculated between: (i) the size of the VOIs of "low" ADC, and (ii) the GS of the corresponding RP lesion/biopsy. Sensitivity, specificity, negative and positive predicative values for GS were calculated. For comparison, the same assessments were carried out for PI-RADS in the targeted biopsy dataset (PI-RADS were not available for the majority of the patients in the prostatectomy and template biopsy datasets).
The area under the curve (AUC) of the receiver operating characteristics (ROC) analysis was computed using a logistic regression model. Statistical analysis was performed using a statistics software package (R, http:// www.R-project.org/). All tests were two-sided. Significance was set at p-value < 0.05.