Maps of cellularity (estimated as number of nuclei per unit area) were calculated for each digitized slide according to the method presented in Figure 2A: HES images were initially converted from the RGB color space into the lightness channel (L) of the Lab color space. An optimum threshold, t_opt, was then used to classify all pixels with L < t_opt as positively stained nuclei, thus creating a binary segmentation mask. The binary mask was subdivided into a regular grid of square subregions, each measuring 500 × 500 pixels (0.115 mm² resolution), to represent the pixels in the final cellularity map. Within each subregion, the ratio of the total area covered by the mask to the mean nucleus cross-sectional area, provided an estimate of nuclear count (the mean ± SE cross-sectional area of nuclei was estimated to be 41.29 ± 0.42 µm² from 1,110 manually contoured nuclei from representative HES histology images including all tissue types). The final cellularity pixel value was thus derived by taking the ratio of the estimated nucleus count to the total area in each subregion.

Our optimization strategy for determining segmentation thresholds is outlined in Figure 2B. A representative region was chosen from each digitized slide measuring 2,000 × 2,000 pixels (0.46 mm × 0.46 mm) and then subdivided into a regular grid of the same dimensions used in cellularity calculation described above. To provide a gold-standard cellularity estimate, a manual nucleus count was performed in each subregion (Figure 2B, i). Using 14 out of the 16 subregions, we then used a Levenberg–Marquardt algorithm to find the luminosity threshold that minimized the sum of absolute differences between the cellularity estimates for each subregion estimated using the above algorithm and the reference gold-standard (Figure 2B, ii). The remaining two subregions provided an unbiased estimate of the uncertainty of the threshold (the cross-validation score, CVS). This process was repeated eight times by leaving-out each combination of neighboring subregions but ensuring a subregion was not re-used. The final threshold used for the digitized slide was chosen as the mean result from cross-validation and the CVS provided an estimate on the uncertainty of the segmentation strategy (a median uncertainty of within 9.42% with a range of 0.98% for the fifth percentile to 42.16% for the 95th percentile across all samples).

Our methodology for extracting luminal fraction (λ) maps from histology is illustrated in Figure 3: images were converted from the RGB into the Lab color space and an unsupervised K-means clustering algorithm (k = 6) was applied on the color channels (A and B). A subset of the segmented clusters was manually selected and converted to binary. Morphological closing was performed using a small disk-shaped structuring element and spatially separate regions were identified using region labeling. Any regions whose size was deemed to be too small/large (e.g., due to background noise) were discarded. Each pixel of the λ-map was computed by calculating the fraction of pixels within each 20 × 20 pixel grid identified as luminal space in these images. This method was implemented in Matlab (Mathworks, Natick, MA, USA).

We simulated the motion of water molecules within the diffusing space derived from histology samples according to the method outlined in Figure 4. Digitized images from HES slides were divided into subregions, each measuring 0.115 mm × 0.115 mm to represent a single pixel in the final sADC map (resolution matched to cellularity maps). Nuclei were segmented within each subregion using the strategy detailed in the preceding section, which provided a two-dimensional set, M(x,y) = M(x), of regions in which the particles were allowed to freely diffuse. The random trajectories of N_p particles, with positions at time t denoted by xit∈R2, were simulated within the segmented diffusing space over N_t time increments using the following algorithm:

t = 0

Uniformly sample xi0∈M

while t < Nt do

for i ∈ Np do

Δxi∼N(0,σ2),σ2=2DfTNt

if xit+Δxi∈M then

xit+1←xit+Δxi

else

xit+1←xit

end if

end for

t = t + 1

end while

N(0,σ2) represents a zero-centered Normal distribution with isotropic covariance σ², D_f is the diffusion coefficient of free particles (set to 3.0 × 10⁻³ mm²/s to represent free water at 37°C) and T is the total time particles are allowed to diffuse. If a particle attempted to go past the edge of the image field-of-view, then the subregion was tessellated to approximate an infinite spatial field. The gradient of the mean-square-displacement curve of all particles over the duration of the simulation time provided an estimate of the simulated-diffusion coefficient of the particles within each sub-region. In our experiment, we used N_p = 5,000 and N_t = 2,630 (see Appendix A in Supplementary Material) with the total diffusion-time, T = 8.67 ms matched to that of our DW-MRI protocol.

Animal care and experiments were carried out in accordance with Norwegian and EU guidelines for care and use of laboratory animals and were approved by the Norwegian National Animal Research Authority and the Norwegian Food Safety Authority (FOTS application 3823). The colony of TRAMP mice used as model organism in all animal experiments, were genetically modified from C57BL/6 mice (Jackson Labs, USA) and established in-house (NTNU, Norway). Genotyping was performed by PCR. The animals were kept in a standardized environment and monitored for general health status and body weight for the duration of the experiments. Male TRAMP (n = 9) and control (n = 6) mice from the same genetic background (C57BL/6) were imaged using MRI every 4 weeks from 8 weeks of age, and terminated at 28–30 weeks or when visual inspection of images indicated unacceptable tumor burden. Mouse body weights were recorded before each MRI session. In this study, results are presented from the final scanning time point only, where parameters from MRI and histology were compared.

MRI was performed on a 7T scanner (Biospec 70/20 Avance III, Bruker Biospin MRI, Ettlingen, Germany) with a volume resonator (86 mm diameter) for RF transmission and a phased array mouse heart surface coil for reception. Mice were anesthetized (≈2% isoflurane in medical air with 36% O₂) for the duration of the MRI scan and positioned on the scanner bed in a prone position. Breathing motion in the pelvic region was reduced by firmly securing the mouse to the scanner bed with adhesive tape across its lower back. The respiration rate was monitored (SA Instruments, USA), and the body temperature was maintained at 37°C by circulating warm water through the bed. The following MR-imaging sequences were used: T2-weighted (T2W) images were acquired in the axial plane with an isotropic in-plane resolution of 0.1 mm and slice thickness of 0.33 mm using a RARE spin echo sequence (TE = 36 ms, TR = 5.5 s). Diffusion-weighted (DW-MRI) were acquired over the same region of the mouse as T2W images to allow for image registration; a fat-suppressed Stejskal–Tanner prepared multi-shot EPI sequence was used with the following parameters: TE = 28.5 ms; TR = 3 s; b-values = 0, 100, 200, 400, 800 s/mm² acquired along three orthogonal directions; averages = 4; matrix size = 128 × 128; slice thickness = 0.9 mm; in-plane resolution 0.2 mm × 0.2 mm; number of EPI segments = 4. ADC maps were calculated in Matlab (MathWorks, Natick, MA, USA) by voxel-wise fitting of the signal (S) averaged over all gradient directions using a monoexponential model for all b-values according to: S(b)=S(0)⋅e−b⋅ADC where S(0) is the signal intensity where b = 0 s/mm². Regions of interest (ROIs) were drawn by DKH on ADC maps around regions corresponding to those drawn on histology. T2W and high-resolution HES images were consulted to ensure corresponding regions were chosen on ADC and cellularity maps.

Upon sacrifice, the genitourinary (GU) tract (prostate, seminal vesicles, emptied bladder) was excised, weighed, and fixed in formalin (10%) for at least 48 h. Samples were embedded such that the sectioning plane was aligned with the MRI images, as described in Ref. (34), and MR-images were used to identify cancerous areas from which histology samples were acquired. Formalin-fixed paraffin embedded samples were sectioned (4 µm slice thickness) and stained with hematoxylin (ChemiTeknik AS, Norway), erythrosine B (Sigma-Aldrich, Norway), and saffron (ChemiTeknik AS, Norway) using an automatic slide stainer (Sakura Tissue-Tek^© Prisma™). Whole HES-stained slides were digitized into jpeg format using a Hamamatsu NanoZoomer XR (Hamamatsu, Japan) scanner (40× magnification) such that each pixel in the high-resolution digitized image had a square resolution of 0.23 µm. Anatomopathological regions of interest were highlighted on digital slides by a pathologist (Daniel N. Rodrigues) using the tools provided by NanoZoomer Digital Pathology (Hamamatsu, Japan). Benign areas were defined for analysis on basis of predominant tissue type and included: fat, smooth muscle, and benign glands. As the distinction between intra-glandular overgrowth of luminal cells, i.e., prostatic intra-epithelial neoplasm (PIN), and invasive gland-forming malignancies, i.e., adenocarcinomas, occasionally depends on markers that highlight micro-anatomical boundaries, PIN/adenocarcinoma was considered as a single category. Overtly invasive tumors with a solid growth pattern and complete absence of gland formation were defined as poorly differentiated carcinomas.

Figure 6A demonstrates the relationship between mean estimates of MR-derived ADC and histology-derived cellularity (in units of 10³ cells/mm²) within the regions of interest defined on both modalities. There is a significant inverse relationship between ADC and tissue cellularity that is well described by our model (equation (5)); we estimate the diffusion coefficient of the fluidic compartment of our model to be D₀ = 2.18 × 10⁻³ mm²/s (95% CI: 1.90, 2.55), the average cross-sectional area of the non-diffusing compartment to be An¯=68.17 µm² (59.52, 71.39) and the tortuosity exponent to be τ = 0.599 (0.597, 0.600). In addition to our model fitting approach, we performed a Leave-One-Animal-Out (LOAO) cross validation of the data that revealed this model was able to predict cellularity based on MR-derived ADC estimates to within a median deviation of 34.5% (normalized-root-mean-square-error = 0.18). This indicates that in this particular animal model, longitudinal measurements of ADC could provide a useful surrogate measurement of cellularity for monitoring disease progression and/or treatment response. We also observe a clear separation between normal and cancerous prostate tissue by both parameters, and ROI clustering of different cancer grades, providing evidence that ADC and cellularity serve as biomarkers for disease aggressiveness in prostate cancer.

Figure 6B illustrates the significant positive correlation (r² = 0.49) observed between average ADC and fractional space occupied by prostate lumina within histology (as a percentage). A linear fit to the slope provided a gradient of 16.3 (95% CI: 8.05, 24.55) with an intercept of −3.33 (95% CI: −12.94, 6.29). LOAO cross-validation analysis revealed this model was able to predict luminal space based on MR-derived ADC estimates to within a median deviation of 29.9% (normalized-root-mean-square-error = 0.20). This demonstrates that ADC may also provide a suitable biomarker for monitoring changes occurring to the glandular structure of the prostate during the onset of disease and/or during treatment. Luminal fraction also provided clear separation between poorly differentiated and well differentiated adenocarcinoma, indicating its utility for assessing tumor grade in the TRAMP model.

A plot of the average MR-derived ADC against sADC revealed a significant linear trend (r² = 0.55) between both parameters (Figure 6C). Linear fitting provided a gradient of 0.77 (95% CI: 0.42, 1.13) and intercept 1.29 (95% CI: 0.88, 1.70) indicating a positive bias in sADC compared to conventional ADC measurements. We attribute this bias to be a consequence of the assumption of free diffusion (diffusion coefficient of 3.0 × 10⁻³ mm²/s for free water at 37^∘C) when simulating diffusion in histology images, ignoring the presence of any cell membranes, organelles, and/or macromolecules in the extracellular space. These data demonstrate that sADC may serve as a useful surrogate for MR-derived ADC when performing pixel-wise analysis of quantitative histology maps; spatial registration between MRI and histology is known to be problematic owing to a decrease in the size of histopathologic specimens following fixation and slide preparation.

Using sADC as a surrogate for in vivo MR-measurements of ADC allowed us to consider pixel-wise relationships between all three histologically derived parameters of prostate tissue morphology. Scatter-plots of these parameters (Figure 7) suggest that there are two principal linear compartments in these data; by performing a 2-component linear mixture model (LMM), we derived maps of the probability of each pixel on histology belonging to each of these two compartments. One compartment (C1) was observed to correspond to luminal epithelium, smooth muscle, and stroma, while the other corresponded to the presence of luminal space (C2). The linear relationships between sADC and cellularity in both of these compartments were revealed to be similar (linear gradient of −7.0 × 10⁻⁶ s/mm⁴ for C1 versus –8.2 × 10⁻⁶ s/mm⁴ for C2).