Analysis of 3D pathology samples using weakly supervised AI.

RESOURCE AVAILABILITY Lead Contact Further information and requests regarding this manuscript should be sent to and will be fulfilled by the lead contact Faisal Mahmood ( faisalmahmood@bwh.harvard.edu ).

Materials availability

This study did not generate any unique reagents.

Data and code availability

The codebase and instructions are available at https://github.com/mahmoodlab/tripath . The OTLS dataset (core needle biopsy, development cohort) is publicly available from the TCIA ( https://www.cancerimagingarchive.net/collection/pca_bx_3dpathology/ ). DOIs are listed in the key resources table . The internal OCT and microCT prostate cancer datasets reported in this study cannot be deposited in a public repository due institutional policy corresponding data sharing and the large file size of 3D tissue resection images. Requests for access to internal data may be sent to the lead author, all requests will be processed following institutional policy governing data sharing and will require a data user agreement. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS Human Participants

The study involved retrospective analysis of pathology specimens and patients were not directly involved or recruited for the study. Specimens used for OTLS imaging were retrieved from genitourinary biorepository at the University of Washington and specimens used for micoCT imaging were retrieved from the Brigham and Women’s Hospital. IRB approvals were received for retrospective image analysis of excess human tissue at the respective institutions. METHOD DETAILS Patient cohorts The clinical end-point utilized for assessing patients’ risk levels is the duration between prostatectomy and the occurrence of biochemical recurrence (BCR), marked by an elevation in prostate-specific antigen (PSA) levels surpassing a defined threshold. The precise PSA threshold triggering intervention by the treating clinician varies among clinicians and the specific laboratory conducting the PSA test, owing to variations in reference ranges across different assays. To account for this inherent variability, we deemed a patient to have reached BCR based on the date of their most recent PSA test before any intervention by the treating clinician, such as modifications in medical treatment or the initiation of radiotherapy. For both the University of Washington (UW) and Brigham and Women’s Hospital (BWH) cohorts, we identified patients who had at least five years of follow-up post-radical prostatectomy (RP) with Gleason grades of 3+3, 3+4, 4+3, and 4+4 (Gleason group 1 ~ 4). For the UW cohort, archived FFPE prostatectomy specimens were collected from n = 74 patients involved in the Canary TMA case-cohort study[ 109 ] with n = 40 patients who experienced BCR within five years of prostatectomy and n = 34 patients who did not. For the BWH cohort, archived FFPE prostatectomy specimens were collected from n = 64 patients with n = 32 patients who experienced BCR within five years of prostatectomy and n = 32 patients who did not. Upon the microCT volumetric image quality check, n = 19 patients were discarded due to high noise level, resulting in the cohort size of n = 45 . The detailed dataset summary can be found in Table S1 .

Data acquisition

Simulation data

We use the simulated 3D digital phantom datasets to enable comparative analyses between different network architectures and data processing approaches (e.g., 2D vs. 3D) while also ensuring every component of the TriPath pipeline is fully functioning[ 67 , 65 , 66 ]. The phantom datasets hold useful advantages over real-world clinical data in that the data is generated according to a set of pre-specified parameters, and thus all morphological characteristics that ought to be captured by the model are known already. This allows us to efficiently evaluate how well the model can capture these characteristics in the larger data regime. The user can specify different cell types, each with their own distributions for size, color, and eccentricity, as well as multiple shapes generated for each cell (e.g., small nuclei residing inside larger cell membranes). To demonstrate its utility, we design a binary class classification dataset in our study, each class populated with different distributions of normal 3D cells (represented as spheroids) and abnormal 3D cells (represented as more eccentric spheroids). Eccentricity corresponds to a mathematical concept of how stretched-out a spheroid is. Normal cells had eccentricity drawn from a normal distribution N(0.25,0.05 2 ), while the eccentricity of the abnormal cells followed N(0.7,0.05 2 ). The length along the axis of symmetry (the semi-major axis) of the normal cells was drawn from N(20,(10/3) 2 ), while the length along the axis of symmetry of the abnormal cells was drawn from N(28,(10/3) 2 ). The samples of class 1 were populated with 90% normal cells and 10% abnormal cells. The samples of class 2 were populated with 66% normal, 34% abnormal cells. Finally, the thickness of the hollow spheroids was set to 3 pixels. The dimension of each sample was 512 × 1024 × 1024 voxels. For each class, 50 random images were generated, each populated with 500 cells. Examples of the simulated dataset can be found in Figure S1 . In a similar context, we create a simulated survival dataset to test the cohort stratification performance of the networks. To be consistent with the OTLS and microCT dataset task, we define two risk groups defined by distinctive morphological characteristics, each with n = 75 . We use the same data generation specification as the simulated classification dataset for simplicity. The corresponding survival timepoints were synthetically generated as follows[ 110 , 71 ]. First, risk scores generated from N(1.0,0.1 2 ) (class 1) and N(2.8,0.1 2 ) (class 2) were assigned to the samples. Based on these risk scores, the survival times were generated with a Cox-Exponential model, with lower (higher) risk scores likely generating longer (shorter) survival times. Finally, the survival times were censored by generating a cutoff point and censoring all survival times past the cutoff, where the cutoff was chosen to have approximately 30% of the samples censored. MicroCT MicroCT scanning for a series of formalin-fixed and paraffin-embedded (FFPE) cancer tissue blocks was done by Versa 620 X-ray Microscope (Carl Zeiss, Inc., Pleasanton, California, USA). Each unstained FFPE sample is attached to a plastic cassette for patient identification, which needs to be removed to avoid the plastic material absorbing X-ray and distorting the image contrast. We remove the plastic cassette by heating the entire block to partially melt the paraffin, which allows the cassette detachment with a razor blade. The separated FFPE block is mounted vertically to a custom-designed steel sample holder and placed on the stage to minimize the thermal vibration of the sample throughout the scanning. For each sample, two scans were performed at different resolutions. First, a quick scan at the low resolution (22.04 μm /voxel) with large field-of-view was performed to capture the whole paraffin block, followed by zooming into tumor-specific locations within the block to capture morphological details at higher resolution (3.98 μm /voxel) (Scout and Zoom protocol). For the high-resolution scan, a microfocus X-ray source with a tube voltage of 40 kV and filament current of 75 μA (3 Watts) was used. A total of 4,501 projection images, with the sample rotated 0.08 degrees (360 degrees/4,501) per projection, were captured for the entire sample on the 16-bit 3,064 pixels by 1,928 pixels flat panel detector and yielded a stack of 1,300 2D images (the depth dimension). Each projection had 15 frames for averaging and 0.5 seconds of exposure for each frame to improve the signal-to-noise ratio (a total of 7.5 seconds for each projection), with the detector recording the raw grayscale intensities for each voxel. The total scan time amounted to 11.5 hours per sample (0.5 hours for the low-resolution scan and 11 hours for the high-resolution scan), and the field-of-view 5.2mm × 12.8mm × 7.68mm (1,300 × 3,200 × 1,920 voxels). All images’ grayscale intensities were first scaled using the paraffin control block (no tumor), then reconstructed with Zeiss Reconstructor software v16 (Carl Zeiss, Inc., Pleasanton, California, USA). During the scaling process, the density of the air was treated as 10 g/cc 3 with an average intensity of 17,030, and the density of the paraffin material (100% wax) was treated as 27 g/cc 3 with an average density of 44,250. The density values were chosen to elevate the noise to well above the minimum 0 and keep the tissue material well below the possible maximum of 65,535(= 2 16 − 1) intensity value. No additional filtering operations were performed on the data. OTLS For each patient, FFPE tissue blocks were identified that correspond to the six prostate regions targeted by urologists in standard sextant and 12-core biopsy procedures. A roughly 1-mm diameter simulated core-needle biopsy was extracted from each of the six blocks, with each biopsy having a tissue volume of roughly 1 × 1 × 15 mm . The simulated biopsies were first deparaffinized in xylene and ethanol, and then stained with a T&E staining protocol[ 16 ] that serves as a fluorescent analog of hematoxylin and eosin (H&E) staining. Specifically, these biopsies were first washed with 100% ethanol twice for 1 hour each to remove any excess xylene and then treated in 70% ethanol for an hour to partially rehydrate them. Each biopsy was then placed in a 0.5 mL Eppendorf tube and stained for 48 hours in 70% ethanol at pH 4 using a 1:200 dilution of Eosin-Y and a 1:500 dilution of To-PRO-3 Iodide at room temperature with gentle agitation. These biopsies were then dehydrated twice with 100% ethanol for 2 hours. Finally, the biopsies were optically cleared by placing them in ethyl cinnamate for 8 hours. A detailed step-by-step protocol and troubleshooting guide can be found in Bishop et al.[ 77 ]. A custom OTLS microscope[ 24 ] was used to image each biopsy across 2 wavelength channels (with laser wavelengths of 488nm and 638nm). Ethyl cinnamate was utilized as the immersion medium, and a multi-channel laser system was used to provide the illumination. Tissues were imaged at near-Nyquist sampling of approximately 0.44 μm /voxel resolution and the volumetric imaging time was approximately 0.5 minutes/ mm 3 of tissue for each wavelength channel. For efficient computational processing, we downsample the data by a factor of 2× to 0.94 μm/ voxel. Each volumetric image amounted to 320 × 520 × 9,500 voxels. The resulting data is saved as 16-bit unsigned integers. Upon inspection of 3D OTLS images of 444 biopsies (6 per patient, for n = 74 patients) by a pathologist (L.D.T.), 171 biopsies that contained tumors (1 to 5 cancer-containing biopsies per patient) were selected for the study.

Volumetric image preprocessing Volume segmentation

We treat the volumetric image as a stack of 2D images and perform tissue segmentation serially on the stack. First, the mean voxel intensity is computed for each image to identify a subset of stacks containing air and images below a user-defined threshold are disregarded before segmentation. Images in the remaining stack are then converted to grayscale color space, median-blurred to suppress edge artifacts, and binarized with modality-specific thresholds. The tissue contours are identified based on the binarized images, and the stack of tissue contours serves as the contour for the volume input. Images with tissue area below a certain threshold are removed to ensure sufficient tissue exists in each image. 3D patching & 2D patching The segmented volume is patched into a set of smaller 2D patches (from a stack of planes) or 3D patches (from a stack of cuboids) to make direct computational processing of the volume feasible. The patch size and the overlap between the patches are chosen to ensure that context is sufficiently covered within each patch and enough patches exist along each dimension. For OTLS, we use 3D patch size of 128 × 128 × 64 voxels (≃ 128 × 128 × 64 μm ). An overlap of 32 voxels along the depth dimension is used to ensure that enough patches exist along the depth dimension, as it is only comprised of 320 voxels. For microCT, we use 128 × 128 × 32 voxels (≃ 512 × 512 × 128 μm ) without any overlap as the size of tissue allows a sufficient number of patches along all dimensions. For 2D patch, we use a non-overlapping patch of 128 × 128 pixels (≃ 128 × 128 μ m for OTLS and 512 × 512 μ m for microCT) for both modalities. For 3D patching, a reference plane is required from which the patching operation along the depth dimension is started. We use the largest plane by tissue area (identified by the tissue contour from the volume segmentation step) as the reference and compute the two-dimensional patch coordinates within the tissue contour. We then perform 3D patching along both directions of the depth dimension starting from the reference plane. The collection of two-dimensional coordinates computed in the reference plane is utilized across the entire volume. Upon completion, we remove 3D patches if more than 50% of the volume (area) constitutes the background to ensure each patch contains sufficient tissue. Additional details on the patch dimensions can be found in Table S2 . After patching, the intensity in each patch is clipped at modality-specific lower and upper thresholds and then normalized to [0,1] for the next feature encoding step. For microCT, the lower threshold is set to 25,000 intensity value and the upper threshold to the top 1% of each tissue volume’s intensity value. For OTLS, the lower threshold is set to 100, and the upper threshold to the top 1% of each tissue volume’s intensity value. For OTLS, we additionally invert the normalized intensity values.

Clinical validation OTLS cohort

We organized a reader study for Gleason grading of the OTLS cohort by recruiting 6 board-certified Genitourinary pathologists (R.C., M.R.D., X.F., P.H., C.V., L.D.T.). For each false-colored OTLS 3D biopsy image, one experienced pathologist (L.D.T) pre-selected the area that contained evidence of cancer. The area was then cropped to a size ranging from 1,024 × 2,000 pixels to 1,024 × 8,000 pixels with a sampling pitch of 1 μm /pixel. In total, two rounds of reader study were conducted. In the first round, 3 slices were taken from the center level along the depth dimension of each cropped biopsy, and ±20 μm from the center, to mimic the standard clinical practice of examining three H&E-stained levels per biopsy (discarding 5 sections or roughly 20 μm between levels). Each pathologist was simultaneously provided the 3 images of the biopsy and was asked to provide the following diagnostic information: primary and secondary Gleason pattern, percentage of Gleason pattern 4, and presence of cribriform. In the second round, which was conducted after 2 months of washout period, all of the slices from the biopsy were shown to each pathologist and was asked to provide the same diagnostic information. The entire reader study was performed with a custom-developed web tool (GroundTruthLab), which allowed pathologists to freely scroll through the image slices within each biopsy. Given the Gleason grade diagnoses, we assess their prognostic value by using linear logistic regression against the binary 5-year BCR status. Specifically, we use a one-hot encoding scheme where a Gleason grade is represented as a 4-dimensional binary vector, i.e. , [1,0,0,0] for grade 3 + 3, [0,1,0,0] for grade 3 + 4, and so on. This allows the different grades to have differential effects on the BCR status. For a patient with multiple cancer-containing biopsies, we take the biopsy with the maximum Gleason grade to represent the patient, per standard practice. We employ 5-fold cross-validation, where the regression model is trained on 80% of the cohort and computes the predicted probability for BCR on the remaining 20% of the cohort. The predictions are aggregated across the folds and the cohort-level AUC is computed. To be consistent with all other experiments, we repeat this procedure 5 times to ensure a specific data split is not favored. We observe that the result is robust to different choices of L2 regularization penalty for the logistic regression.

MicroCT cohort

For each patient, we digitize a H&E slice taken from the same block that was subject to microCT scanning. The H&E slice is scanned as a WSI at 10x magnification (1 μ m/pixel) per standard clinical practice for Gleason grading. The WSI is segmented and patched at 256 × 256 pixels without overlap. A 2D feature encoder (Resnet50) encodes each of the 2D WSI patch in the set as a 1,024-dimensional feature. The set of patches is then processed with an attention-based aggregation module, in the same manner as other TriPath experiments, and yields patient-level predicted risk. We also repeat the same procedure with an ROI (4 × 4 mm ) within each WSI that matches the lateral field of view of the microCT images to minimize potential bias arising from different fields of view between H&E and microCT datasets.

Model architecture

Feature encoder choice Feature encoders serve the purpose of extracting and encoding compressed and representative descriptor h j ∈ ℝ K , j = 1 , ⋯ , J of the patch input x j ∈ ℝ L × D × H × W (3D patch) or x j ∈ ℝ L × H × W (2D patch), where K corresponds to the encoded feature dimension, J denotes the number of patches, L denotes number of input channels, and D , H , W denotes the depth, height, and width dimension respectively. TriPath provides the choice between a range of 2D and 3D pretrained feature encoders based on convolutional neural networks (CNN) or Vision Transformer (ViT)[ 111 ] for transfer learning. For the 3D feature encoders, TriPath provides spatiotemporal CNN[ 108 , 35 ] pretrained on a large collection of human action recognition videos[ 112 ] and video sliding-window transformer (Video SwinViT)[ 113 , 114 ] pretrained on a human action recognition videos or 3D medical imaging dataset[ 46 ]. For the 2D feature encoders, TriPath provides ResNet-50[ 115 ] pretrained on natural images and SwinViT pretrained on a large collection of histopathology images[ 116 ] or natural images[ 117 ]. Due to the scarcity of patient-level labels (clinical endpoints) for 3D pathology datasets and generally larger encoder network size for processing the depth dimension, fine-tuning the pretrained feature encoders or training the encoders from scratch in the volumetric image data domain results in network overfitting and poor generalization performance. To address the domain gap from transfer learning, we apply a fully-connected linear layer to the feature encoder outputs { h j } j = 1 J , parameterized by W enc ∈ ℝ 256 × K and b enc ∈ ℝ 256 , followed by GeLU nonlinearity. This further converts patch feature h j from the feature encoder to a more-compressed and domain-specific feature z j ∈ ℝ 256 conducive to downstream tasks with better generalization performance: (1) z j = G e L U ( W enc h j + b enc ) For our study, we used deep residual CNN (ResNet-50), truncated after the third residual block and pretrained on natural images (ImageNet) for the 2D experiments and spatiotemporal CNN with ResNet-50 backbone[ 108 , 35 ] pretrained on action recognition videos (Kinetics-400[ 112 ]), both of which yield K = 1 , 024 . The choice of the 2D and 3D feature encoders was driven by 1) a further ablation study on feature encoders showing that the spatiotemporal CNN performs consistently well for both OTLS and microCT datasets ( Figure S3 ) and 2) the choice ensuring a fair comparison of network architecture, as both networks are based on deep residual components. Nevertheless, we encourage testing out different feature encoders as other encoders could yield the best performance depending on the task. Feature encoding step As most feature encoders take three-channel RGB inputs, we emulate the setting by replicating channel information. For the dual-channel OTLS data, we replicate the nuclear channel data across the first two channels and set the eosin channel as the third. For the single-channel microCT data, we replicate the data across all three channels. For the feature encoding step, we use a batch size of 500 for 2D patches and 100 for 3D patches. For CNN-based feature encoders, the immediate output of the feature encoder is 3-dimensional for a 2D patch ( K , H ˜ , W ˜ ) and 4-dimensional for a 3D patch ( K , D ˜ , H ˜ , W ˜ ), where D ˜ , H ˜ , W ˜ correspond to the downsampled depth, height, width dimension respectively. The intermediate features get compressed to one-dimensional feature h j ∈ ℝ K with adaptive average-spatial pooling operation and subsequently to z j ∈ ℝ 256 with the fully-connected network. For ViT-based feature encoders, we treat the CLS token output of the ViT as h j and subsequently to z j ∈ ℝ 256 with the fully-connected network. Aggregation module The patching and feature encoding operation results in a collection of 256-dimensional features (also referred to as instances) { z j } j = 1 J , constituting the volume with a single patient-level supervisory label, a setting which is referred to as multiple instance learning (MIL). It is also referred to as weakly-supervised learning, due to the substantial size of the input (number of patches) in comparison to the supervisory label. To this end, we use an attention-based aggregation module[ 48 , 40 ], a lightweight attention network that learns to automatically compute the importance score of each patch feature and aggregates by weighted-averaging the features to form a single volume-level feature. The attention network consists of three sets of parameters V ∈ ℝ 64 × 256 , U ∈ ℝ 64 × 256 , and W ∈ ℝ 1 × 64 . The network assigns an importance score a j ∈ [ 0 , 1 ] to feature z j : (2) a j = exp ( W ( tanh ( V z j ) ⊙ s i g m ( U z j ) ) ) ∑ j ′ = 1 J exp ( W ( tanh ( V z j ′ ) ⊙ s i g m ( U z j ′ ) ) ) . with tanh and sigm denoting hyperbolic tangent and sigmoid function respectively, and ⊙ denoting element-wise multiplication operation. A high score ( a j close to 1) indicates that the corresponding patch is very relevant for sample-level risk prediction, with a low score ( a j close to 0) indicating no prognostic value. Finally, the volume-level feature z volume is computed as (3) z volume = ∑ j = 1 J a j z j ∈ ℝ 256 . Classification module The volume-level feature is fed into the final classification layer parameterized by W cls ∈ ℝ 1 × 256 and bias b cls ∈ ℝ resulting in the probability for the high-risk group p ∈ [ 0 , 1 ] , (4) p = s i g m ( W cls z volume + b cls ) .

Training & Evaluation

We train all the networks for a fixed number of 50 epochs and the initial learning rate of 2×10 −4 with the cosine decay scheduler. We use AdamW optimizer with default parameters of β 1 = 0.9 and β 2 = 0.999 , with weight decay of 5 × 10 −4 . We use a mini-batch size of 1 patient sample, pooling together patches across samples if multiple tissue samples exist for the patient. We also employ a gradient accumulation of 10 training samples for training stability. For each volume, we randomly sample 50% of the patches as a means of data augmentation to prevent overfitting and also inject diversity into training samples. Rather than sampling randomly from the entire volume, we sample 50% of the patches per plane or cuboid to ensure all depths are equally accounted for. In addition, we employ a heavy dropout of p = 0.5 after each fully-connected layer. We also apply conventional data augmentation schemes to patches, such as rotation and intensity jittering on top of random sampling of the patches. We use the binary cross-entropy loss for the loss function.

Integrated gradients interpretability analysis

The integrated gradient (IG) method[ 73 , 74 ] assesses the relationship between an input to a network and the corresponding prediction. In our study, the input and the prediction correspond to the set of instance features { h j } j = 1 J and the probability for high-risk group p , respectively. The IG method assigns an IG score for each input ( h j in our case), signifying the strength of each input’s influence on the prediction, with the sign of the score indicating the direction of influence. In the context of prognosis, this can directly be translated as positive IG values increasing the risk (unfavorable prognosis) and negative IG values decreasing the risk (favorable prognosis). The IG values close to 0 have no prognostic influence. Denoting F as the sequence of the fully-connected layer for feature encoder, attention aggregation, and classification modules, i.e. , p = F ( { h j } j = 1 J ) , and also M as the total number of IG interpolation steps, h j , k as the k th element of the feature h j , we can compute the IG score for h j as (5) I G ( h j ) = ∑ k = 1 K h j , k × 1 M ∑ m = 1 M ∂ F ( { m M ⋅ h j } j = 1 J ) ∂ h j , k , assuming zero feature baseline. Once all the IG scores are computed for the patch features of a patient, we normalize the negative IG values to [−1,0] and the positive IG values to ( 0 , 1 ] , to ensure the sign and the influence of a patch does not get flipped. Cross-modal analysis The OTLS and microCT dataset characteristics are vastly different and thus necessary adjustments are required for fair cross-modal experiments. To this end, we downsample the OTLS dataset by a factor of 4 (from 1 μm /voxel to 4 μm /voxel) and use only the nuclear channel to match with 4 μm /voxel single-channel microCT dataset, producing the converted OTLS dataset. Each of these adjustments results in information loss and contributes to the drop in test AUC for OTLS dataset ( i.e. , trained and tested on single-channel 4 μm /voxel OTLS data vs. dual-channel 1 μm /voxel OTLS data). No adjustments were made to the microCT dataset. Five models trained in a cross-validation setting on one cohort are applied to the other cohort to compute the test AUC.

Visualization

False-coloring the raw input For the OTLS dataset, we use a false-coloring module[ 80 ] which utilizes the physics model (Beer-Lambert law absorption of light) of the dual-channel information of the raw OTLS data to render hematoxylin and eosin appearance. Integrated gradients heatmap To generate fine-grained 3D IG heatmaps, we use 3D cuboid patches with 75% overlap in 2D plane direction and 50% overlap along the depth dimension, to reduce blocky effects. To compute an IG score for a given region, the raw IG scores (prior to normalization) of all the patches covering the region are accumulated and divided by the number of overlapping patches. These IG scores are then normalized in the manner described in the previous section. A coolwarm colormap, with red and blue colors indicating positive and negative IG values respectively, is then applied to the normalized IG scores, which is then overlaid on the raw volumetric image with a transparency value of 0.4. The IG heatmaps are shown in Figure 2E , Figure 3E , Figure S5 , Figure S6 .

QUANTIFICATION AND STATISTICAL ANALYSIS

Training and testing on both the microCT and OTLS cohort were performed with 5-fold cross-validation stratified by the BCR status. Upon completion of training for all five folds, the predicted probabilities from each fold are aggregated and cohort-level AUC is computed. We repeated the experiments with five different random splits of the training and testing data, such that the result is not biased towards a certain data split. For the comparison of high and low-risk survival curves, we used the log-rank test. The stratification of the high and low-risk groups was always performed at 50 percentile of the predicted probability for the high-risk group, unless specified otherwise. For each baseline, we displayed the survival curve for the data split that results in the model with the highest cohort-level AUC (among five runs). Since we deal with a small number of samples and treat the risk stratification task as that of a binary classification problem, not as a survival task (via Cox regression), we primarily rely on AUC to assess TriPath performance, with the survival curve as an auxiliary metric. We used the two-sided unpaired t-test to assess the statistical significance of two individual groups. Spearman’s correlation coefficient r was used to assess the correlation between two quantities. Differences in the compared group were considered statistically significant when P values were smaller than 0.05 ( P > 0.05, not significant; * P ≤ 0.05, ** P ≤ 0.01, *** P ≤ 0.001, **** P ≤ 0.0001).

Computational hardware and software

All volumetric images were processed on AMD ® Ryzen multicore CPUs (central processing units) and a total of 6 NVIDIA GeForce RTX 3090 GPUs (graphics processing units) using our custom, publicly available TriPath package processing pipeline implemented in Python (version 3.10.9). TriPath uses pillow (version 9.2.0) and opencv-python (version 4.6.0) for image processing. All deep learning implementations were performed with PyTorch (2.0.1). The spatiotemporal CNN and Swin transformer feature encoders were adapted from pytorchvideo (version 0.1.5). 3D Visualization was accomplished via napari (version 0.4.16). Plots were generated in Python using matplotlib (version 3.5.2) and numpy (version 1.22.4) was used for vectorized numerical computation. Other Python libraries used to support data analysis include pandas (version 1.4.3), scipy (version 1.9.0), torchvision (version 0.15.2), and timm (version 0.9.7). The scientific computing library scikit-learn (version 1.0.2) was used to compute various classification metrics and estimate the AUC ROC. The survival analysis was performed with lifelines (0.26.0). The interactive demo website was developed using THREE.js (version 0.152.2) and jQuery (version 3.6.0). ADDITIONAL RESOURCES Description: URL Selected OTLS and microCT raw images and heatmaps can be visualized in our interactive demo website ( https://mamba-demo.github.io/demo/ ). This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Materials availability

This study did not generate any unique reagents.

EXPERIMENTAL MODEL AND SUBJECT DETAILS Human Participants

The study involved retrospective analysis of pathology specimens and patients were not directly involved or recruited for the study. Specimens used for OTLS imaging were retrieved from genitourinary biorepository at the University of Washington and specimens used for micoCT imaging were retrieved from the Brigham and Women’s Hospital. IRB approvals were received for retrospective image analysis of excess human tissue at the respective institutions.

METHOD DETAILS Patient cohorts The clinical end-point utilized for assessing patients’ risk levels is the duration between prostatectomy and the occurrence of biochemical recurrence (BCR), marked by an elevation in prostate-specific antigen (PSA) levels surpassing a defined threshold. The precise PSA threshold triggering intervention by the treating clinician varies among clinicians and the specific laboratory conducting the PSA test, owing to variations in reference ranges across different assays. To account for this inherent variability, we deemed a patient to have reached BCR based on the date of their most recent PSA test before any intervention by the treating clinician, such as modifications in medical treatment or the initiation of radiotherapy. For both the University of Washington (UW) and Brigham and Women’s Hospital (BWH) cohorts, we identified patients who had at least five years of follow-up post-radical prostatectomy (RP) with Gleason grades of 3+3, 3+4, 4+3, and 4+4 (Gleason group 1 ~ 4). For the UW cohort, archived FFPE prostatectomy specimens were collected from n = 74 patients involved in the Canary TMA case-cohort study[ 109 ] with n = 40 patients who experienced BCR within five years of prostatectomy and n = 34 patients who did not. For the BWH cohort, archived FFPE prostatectomy specimens were collected from n = 64 patients with n = 32 patients who experienced BCR within five years of prostatectomy and n = 32 patients who did not. Upon the microCT volumetric image quality check, n = 19 patients were discarded due to high noise level, resulting in the cohort size of n = 45 . The detailed dataset summary can be found in Table S1 .

Data acquisition

Simulation data

We use the simulated 3D digital phantom datasets to enable comparative analyses between different network architectures and data processing approaches (e.g., 2D vs. 3D) while also ensuring every component of the TriPath pipeline is fully functioning[ 67 , 65 , 66 ]. The phantom datasets hold useful advantages over real-world clinical data in that the data is generated according to a set of pre-specified parameters, and thus all morphological characteristics that ought to be captured by the model are known already. This allows us to efficiently evaluate how well the model can capture these characteristics in the larger data regime. The user can specify different cell types, each with their own distributions for size, color, and eccentricity, as well as multiple shapes generated for each cell (e.g., small nuclei residing inside larger cell membranes). To demonstrate its utility, we design a binary class classification dataset in our study, each class populated with different distributions of normal 3D cells (represented as spheroids) and abnormal 3D cells (represented as more eccentric spheroids). Eccentricity corresponds to a mathematical concept of how stretched-out a spheroid is. Normal cells had eccentricity drawn from a normal distribution N(0.25,0.05 2 ), while the eccentricity of the abnormal cells followed N(0.7,0.05 2 ). The length along the axis of symmetry (the semi-major axis) of the normal cells was drawn from N(20,(10/3) 2 ), while the length along the axis of symmetry of the abnormal cells was drawn from N(28,(10/3) 2 ). The samples of class 1 were populated with 90% normal cells and 10% abnormal cells. The samples of class 2 were populated with 66% normal, 34% abnormal cells. Finally, the thickness of the hollow spheroids was set to 3 pixels. The dimension of each sample was 512 × 1024 × 1024 voxels. For each class, 50 random images were generated, each populated with 500 cells. Examples of the simulated dataset can be found in Figure S1 . In a similar context, we create a simulated survival dataset to test the cohort stratification performance of the networks. To be consistent with the OTLS and microCT dataset task, we define two risk groups defined by distinctive morphological characteristics, each with n = 75 . We use the same data generation specification as the simulated classification dataset for simplicity. The corresponding survival timepoints were synthetically generated as follows[ 110 , 71 ]. First, risk scores generated from N(1.0,0.1 2 ) (class 1) and N(2.8,0.1 2 ) (class 2) were assigned to the samples. Based on these risk scores, the survival times were generated with a Cox-Exponential model, with lower (higher) risk scores likely generating longer (shorter) survival times. Finally, the survival times were censored by generating a cutoff point and censoring all survival times past the cutoff, where the cutoff was chosen to have approximately 30% of the samples censored. MicroCT MicroCT scanning for a series of formalin-fixed and paraffin-embedded (FFPE) cancer tissue blocks was done by Versa 620 X-ray Microscope (Carl Zeiss, Inc., Pleasanton, California, USA). Each unstained FFPE sample is attached to a plastic cassette for patient identification, which needs to be removed to avoid the plastic material absorbing X-ray and distorting the image contrast. We remove the plastic cassette by heating the entire block to partially melt the paraffin, which allows the cassette detachment with a razor blade. The separated FFPE block is mounted vertically to a custom-designed steel sample holder and placed on the stage to minimize the thermal vibration of the sample throughout the scanning. For each sample, two scans were performed at different resolutions. First, a quick scan at the low resolution (22.04 μm /voxel) with large field-of-view was performed to capture the whole paraffin block, followed by zooming into tumor-specific locations within the block to capture morphological details at higher resolution (3.98 μm /voxel) (Scout and Zoom protocol). For the high-resolution scan, a microfocus X-ray source with a tube voltage of 40 kV and filament current of 75 μA (3 Watts) was used. A total of 4,501 projection images, with the sample rotated 0.08 degrees (360 degrees/4,501) per projection, were captured for the entire sample on the 16-bit 3,064 pixels by 1,928 pixels flat panel detector and yielded a stack of 1,300 2D images (the depth dimension). Each projection had 15 frames for averaging and 0.5 seconds of exposure for each frame to improve the signal-to-noise ratio (a total of 7.5 seconds for each projection), with the detector recording the raw grayscale intensities for each voxel. The total scan time amounted to 11.5 hours per sample (0.5 hours for the low-resolution scan and 11 hours for the high-resolution scan), and the field-of-view 5.2mm × 12.8mm × 7.68mm (1,300 × 3,200 × 1,920 voxels). All images’ grayscale intensities were first scaled using the paraffin control block (no tumor), then reconstructed with Zeiss Reconstructor software v16 (Carl Zeiss, Inc., Pleasanton, California, USA). During the scaling process, the density of the air was treated as 10 g/cc 3 with an average intensity of 17,030, and the density of the paraffin material (100% wax) was treated as 27 g/cc 3 with an average density of 44,250. The density values were chosen to elevate the noise to well above the minimum 0 and keep the tissue material well below the possible maximum of 65,535(= 2 16 − 1) intensity value. No additional filtering operations were performed on the data. OTLS For each patient, FFPE tissue blocks were identified that correspond to the six prostate regions targeted by urologists in standard sextant and 12-core biopsy procedures. A roughly 1-mm diameter simulated core-needle biopsy was extracted from each of the six blocks, with each biopsy having a tissue volume of roughly 1 × 1 × 15 mm . The simulated biopsies were first deparaffinized in xylene and ethanol, and then stained with a T&E staining protocol[ 16 ] that serves as a fluorescent analog of hematoxylin and eosin (H&E) staining. Specifically, these biopsies were first washed with 100% ethanol twice for 1 hour each to remove any excess xylene and then treated in 70% ethanol for an hour to partially rehydrate them. Each biopsy was then placed in a 0.5 mL Eppendorf tube and stained for 48 hours in 70% ethanol at pH 4 using a 1:200 dilution of Eosin-Y and a 1:500 dilution of To-PRO-3 Iodide at room temperature with gentle agitation. These biopsies were then dehydrated twice with 100% ethanol for 2 hours. Finally, the biopsies were optically cleared by placing them in ethyl cinnamate for 8 hours. A detailed step-by-step protocol and troubleshooting guide can be found in Bishop et al.[ 77 ]. A custom OTLS microscope[ 24 ] was used to image each biopsy across 2 wavelength channels (with laser wavelengths of 488nm and 638nm). Ethyl cinnamate was utilized as the immersion medium, and a multi-channel laser system was used to provide the illumination. Tissues were imaged at near-Nyquist sampling of approximately 0.44 μm /voxel resolution and the volumetric imaging time was approximately 0.5 minutes/ mm 3 of tissue for each wavelength channel. For efficient computational processing, we downsample the data by a factor of 2× to 0.94 μm/ voxel. Each volumetric image amounted to 320 × 520 × 9,500 voxels. The resulting data is saved as 16-bit unsigned integers. Upon inspection of 3D OTLS images of 444 biopsies (6 per patient, for n = 74 patients) by a pathologist (L.D.T.), 171 biopsies that contained tumors (1 to 5 cancer-containing biopsies per patient) were selected for the study.

Volumetric image preprocessing Volume segmentation

We treat the volumetric image as a stack of 2D images and perform tissue segmentation serially on the stack. First, the mean voxel intensity is computed for each image to identify a subset of stacks containing air and images below a user-defined threshold are disregarded before segmentation. Images in the remaining stack are then converted to grayscale color space, median-blurred to suppress edge artifacts, and binarized with modality-specific thresholds. The tissue contours are identified based on the binarized images, and the stack of tissue contours serves as the contour for the volume input. Images with tissue area below a certain threshold are removed to ensure sufficient tissue exists in each image. 3D patching & 2D patching The segmented volume is patched into a set of smaller 2D patches (from a stack of planes) or 3D patches (from a stack of cuboids) to make direct computational processing of the volume feasible. The patch size and the overlap between the patches are chosen to ensure that context is sufficiently covered within each patch and enough patches exist along each dimension. For OTLS, we use 3D patch size of 128 × 128 × 64 voxels (≃ 128 × 128 × 64 μm ). An overlap of 32 voxels along the depth dimension is used to ensure that enough patches exist along the depth dimension, as it is only comprised of 320 voxels. For microCT, we use 128 × 128 × 32 voxels (≃ 512 × 512 × 128 μm ) without any overlap as the size of tissue allows a sufficient number of patches along all dimensions. For 2D patch, we use a non-overlapping patch of 128 × 128 pixels (≃ 128 × 128 μ m for OTLS and 512 × 512 μ m for microCT) for both modalities. For 3D patching, a reference plane is required from which the patching operation along the depth dimension is started. We use the largest plane by tissue area (identified by the tissue contour from the volume segmentation step) as the reference and compute the two-dimensional patch coordinates within the tissue contour. We then perform 3D patching along both directions of the depth dimension starting from the reference plane. The collection of two-dimensional coordinates computed in the reference plane is utilized across the entire volume. Upon completion, we remove 3D patches if more than 50% of the volume (area) constitutes the background to ensure each patch contains sufficient tissue. Additional details on the patch dimensions can be found in Table S2 . After patching, the intensity in each patch is clipped at modality-specific lower and upper thresholds and then normalized to [0,1] for the next feature encoding step. For microCT, the lower threshold is set to 25,000 intensity value and the upper threshold to the top 1% of each tissue volume’s intensity value. For OTLS, the lower threshold is set to 100, and the upper threshold to the top 1% of each tissue volume’s intensity value. For OTLS, we additionally invert the normalized intensity values.

Clinical validation OTLS cohort

We organized a reader study for Gleason grading of the OTLS cohort by recruiting 6 board-certified Genitourinary pathologists (R.C., M.R.D., X.F., P.H., C.V., L.D.T.). For each false-colored OTLS 3D biopsy image, one experienced pathologist (L.D.T) pre-selected the area that contained evidence of cancer. The area was then cropped to a size ranging from 1,024 × 2,000 pixels to 1,024 × 8,000 pixels with a sampling pitch of 1 μm /pixel. In total, two rounds of reader study were conducted. In the first round, 3 slices were taken from the center level along the depth dimension of each cropped biopsy, and ±20 μm from the center, to mimic the standard clinical practice of examining three H&E-stained levels per biopsy (discarding 5 sections or roughly 20 μm between levels). Each pathologist was simultaneously provided the 3 images of the biopsy and was asked to provide the following diagnostic information: primary and secondary Gleason pattern, percentage of Gleason pattern 4, and presence of cribriform. In the second round, which was conducted after 2 months of washout period, all of the slices from the biopsy were shown to each pathologist and was asked to provide the same diagnostic information. The entire reader study was performed with a custom-developed web tool (GroundTruthLab), which allowed pathologists to freely scroll through the image slices within each biopsy. Given the Gleason grade diagnoses, we assess their prognostic value by using linear logistic regression against the binary 5-year BCR status. Specifically, we use a one-hot encoding scheme where a Gleason grade is represented as a 4-dimensional binary vector, i.e. , [1,0,0,0] for grade 3 + 3, [0,1,0,0] for grade 3 + 4, and so on. This allows the different grades to have differential effects on the BCR status. For a patient with multiple cancer-containing biopsies, we take the biopsy with the maximum Gleason grade to represent the patient, per standard practice. We employ 5-fold cross-validation, where the regression model is trained on 80% of the cohort and computes the predicted probability for BCR on the remaining 20% of the cohort. The predictions are aggregated across the folds and the cohort-level AUC is computed. To be consistent with all other experiments, we repeat this procedure 5 times to ensure a specific data split is not favored. We observe that the result is robust to different choices of L2 regularization penalty for the logistic regression.

MicroCT cohort

For each patient, we digitize a H&E slice taken from the same block that was subject to microCT scanning. The H&E slice is scanned as a WSI at 10x magnification (1 μ m/pixel) per standard clinical practice for Gleason grading. The WSI is segmented and patched at 256 × 256 pixels without overlap. A 2D feature encoder (Resnet50) encodes each of the 2D WSI patch in the set as a 1,024-dimensional feature. The set of patches is then processed with an attention-based aggregation module, in the same manner as other TriPath experiments, and yields patient-level predicted risk. We also repeat the same procedure with an ROI (4 × 4 mm ) within each WSI that matches the lateral field of view of the microCT images to minimize potential bias arising from different fields of view between H&E and microCT datasets.

Model architecture

Feature encoder choice Feature encoders serve the purpose of extracting and encoding compressed and representative descriptor h j ∈ ℝ K , j = 1 , ⋯ , J of the patch input x j ∈ ℝ L × D × H × W (3D patch) or x j ∈ ℝ L × H × W (2D patch), where K corresponds to the encoded feature dimension, J denotes the number of patches, L denotes number of input channels, and D , H , W denotes the depth, height, and width dimension respectively. TriPath provides the choice between a range of 2D and 3D pretrained feature encoders based on convolutional neural networks (CNN) or Vision Transformer (ViT)[ 111 ] for transfer learning. For the 3D feature encoders, TriPath provides spatiotemporal CNN[ 108 , 35 ] pretrained on a large collection of human action recognition videos[ 112 ] and video sliding-window transformer (Video SwinViT)[ 113 , 114 ] pretrained on a human action recognition videos or 3D medical imaging dataset[ 46 ]. For the 2D feature encoders, TriPath provides ResNet-50[ 115 ] pretrained on natural images and SwinViT pretrained on a large collection of histopathology images[ 116 ] or natural images[ 117 ]. Due to the scarcity of patient-level labels (clinical endpoints) for 3D pathology datasets and generally larger encoder network size for processing the depth dimension, fine-tuning the pretrained feature encoders or training the encoders from scratch in the volumetric image data domain results in network overfitting and poor generalization performance. To address the domain gap from transfer learning, we apply a fully-connected linear layer to the feature encoder outputs { h j } j = 1 J , parameterized by W enc ∈ ℝ 256 × K and b enc ∈ ℝ 256 , followed by GeLU nonlinearity. This further converts patch feature h j from the feature encoder to a more-compressed and domain-specific feature z j ∈ ℝ 256 conducive to downstream tasks with better generalization performance: (1) z j = G e L U ( W enc h j + b enc ) For our study, we used deep residual CNN (ResNet-50), truncated after the third residual block and pretrained on natural images (ImageNet) for the 2D experiments and spatiotemporal CNN with ResNet-50 backbone[ 108 , 35 ] pretrained on action recognition videos (Kinetics-400[ 112 ]), both of which yield K = 1 , 024 . The choice of the 2D and 3D feature encoders was driven by 1) a further ablation study on feature encoders showing that the spatiotemporal CNN performs consistently well for both OTLS and microCT datasets ( Figure S3 ) and 2) the choice ensuring a fair comparison of network architecture, as both networks are based on deep residual components. Nevertheless, we encourage testing out different feature encoders as other encoders could yield the best performance depending on the task. Feature encoding step As most feature encoders take three-channel RGB inputs, we emulate the setting by replicating channel information. For the dual-channel OTLS data, we replicate the nuclear channel data across the first two channels and set the eosin channel as the third. For the single-channel microCT data, we replicate the data across all three channels. For the feature encoding step, we use a batch size of 500 for 2D patches and 100 for 3D patches. For CNN-based feature encoders, the immediate output of the feature encoder is 3-dimensional for a 2D patch ( K , H ˜ , W ˜ ) and 4-dimensional for a 3D patch ( K , D ˜ , H ˜ , W ˜ ), where D ˜ , H ˜ , W ˜ correspond to the downsampled depth, height, width dimension respectively. The intermediate features get compressed to one-dimensional feature h j ∈ ℝ K with adaptive average-spatial pooling operation and subsequently to z j ∈ ℝ 256 with the fully-connected network. For ViT-based feature encoders, we treat the CLS token output of the ViT as h j and subsequently to z j ∈ ℝ 256 with the fully-connected network. Aggregation module The patching and feature encoding operation results in a collection of 256-dimensional features (also referred to as instances) { z j } j = 1 J , constituting the volume with a single patient-level supervisory label, a setting which is referred to as multiple instance learning (MIL). It is also referred to as weakly-supervised learning, due to the substantial size of the input (number of patches) in comparison to the supervisory label. To this end, we use an attention-based aggregation module[ 48 , 40 ], a lightweight attention network that learns to automatically compute the importance score of each patch feature and aggregates by weighted-averaging the features to form a single volume-level feature. The attention network consists of three sets of parameters V ∈ ℝ 64 × 256 , U ∈ ℝ 64 × 256 , and W ∈ ℝ 1 × 64 . The network assigns an importance score a j ∈ [ 0 , 1 ] to feature z j : (2) a j = exp ( W ( tanh ( V z j ) ⊙ s i g m ( U z j ) ) ) ∑ j ′ = 1 J exp ( W ( tanh ( V z j ′ ) ⊙ s i g m ( U z j ′ ) ) ) . with tanh and sigm denoting hyperbolic tangent and sigmoid function respectively, and ⊙ denoting element-wise multiplication operation. A high score ( a j close to 1) indicates that the corresponding patch is very relevant for sample-level risk prediction, with a low score ( a j close to 0) indicating no prognostic value. Finally, the volume-level feature z volume is computed as (3) z volume = ∑ j = 1 J a j z j ∈ ℝ 256 . Classification module The volume-level feature is fed into the final classification layer parameterized by W cls ∈ ℝ 1 × 256 and bias b cls ∈ ℝ resulting in the probability for the high-risk group p ∈ [ 0 , 1 ] , (4) p = s i g m ( W cls z volume + b cls ) .

Training & Evaluation

We train all the networks for a fixed number of 50 epochs and the initial learning rate of 2×10 −4 with the cosine decay scheduler. We use AdamW optimizer with default parameters of β 1 = 0.9 and β 2 = 0.999 , with weight decay of 5 × 10 −4 . We use a mini-batch size of 1 patient sample, pooling together patches across samples if multiple tissue samples exist for the patient. We also employ a gradient accumulation of 10 training samples for training stability. For each volume, we randomly sample 50% of the patches as a means of data augmentation to prevent overfitting and also inject diversity into training samples. Rather than sampling randomly from the entire volume, we sample 50% of the patches per plane or cuboid to ensure all depths are equally accounted for. In addition, we employ a heavy dropout of p = 0.5 after each fully-connected layer. We also apply conventional data augmentation schemes to patches, such as rotation and intensity jittering on top of random sampling of the patches. We use the binary cross-entropy loss for the loss function.

Integrated gradients interpretability analysis

The integrated gradient (IG) method[ 73 , 74 ] assesses the relationship between an input to a network and the corresponding prediction. In our study, the input and the prediction correspond to the set of instance features { h j } j = 1 J and the probability for high-risk group p , respectively. The IG method assigns an IG score for each input ( h j in our case), signifying the strength of each input’s influence on the prediction, with the sign of the score indicating the direction of influence. In the context of prognosis, this can directly be translated as positive IG values increasing the risk (unfavorable prognosis) and negative IG values decreasing the risk (favorable prognosis). The IG values close to 0 have no prognostic influence. Denoting F as the sequence of the fully-connected layer for feature encoder, attention aggregation, and classification modules, i.e. , p = F ( { h j } j = 1 J ) , and also M as the total number of IG interpolation steps, h j , k as the k th element of the feature h j , we can compute the IG score for h j as (5) I G ( h j ) = ∑ k = 1 K h j , k × 1 M ∑ m = 1 M ∂ F ( { m M ⋅ h j } j = 1 J ) ∂ h j , k , assuming zero feature baseline. Once all the IG scores are computed for the patch features of a patient, we normalize the negative IG values to [−1,0] and the positive IG values to ( 0 , 1 ] , to ensure the sign and the influence of a patch does not get flipped. Cross-modal analysis The OTLS and microCT dataset characteristics are vastly different and thus necessary adjustments are required for fair cross-modal experiments. To this end, we downsample the OTLS dataset by a factor of 4 (from 1 μm /voxel to 4 μm /voxel) and use only the nuclear channel to match with 4 μm /voxel single-channel microCT dataset, producing the converted OTLS dataset. Each of these adjustments results in information loss and contributes to the drop in test AUC for OTLS dataset ( i.e. , trained and tested on single-channel 4 μm /voxel OTLS data vs. dual-channel 1 μm /voxel OTLS data). No adjustments were made to the microCT dataset. Five models trained in a cross-validation setting on one cohort are applied to the other cohort to compute the test AUC.

Visualization

False-coloring the raw input For the OTLS dataset, we use a false-coloring module[ 80 ] which utilizes the physics model (Beer-Lambert law absorption of light) of the dual-channel information of the raw OTLS data to render hematoxylin and eosin appearance. Integrated gradients heatmap To generate fine-grained 3D IG heatmaps, we use 3D cuboid patches with 75% overlap in 2D plane direction and 50% overlap along the depth dimension, to reduce blocky effects. To compute an IG score for a given region, the raw IG scores (prior to normalization) of all the patches covering the region are accumulated and divided by the number of overlapping patches. These IG scores are then normalized in the manner described in the previous section. A coolwarm colormap, with red and blue colors indicating positive and negative IG values respectively, is then applied to the normalized IG scores, which is then overlaid on the raw volumetric image with a transparency value of 0.4. The IG heatmaps are shown in Figure 2E , Figure 3E , Figure S5 , Figure S6 .

Supplementary Material 1 Figure S1: 3D phantom datasets and analysis with TriPath, related to Figure 1 (A) Examples of single-channel 3D phantom data samples for the binary classification task ( n = 100 ), false-colored for different cell types. Samples from the first class are dominated by normal cells (blue), while samples from the second class are dominated by abnormal cells with large eccentricity (red). (B) Binary classification task AUC for TriPath trained and tested on a random plane from each volume (random plane), the targeted plane that contains both cell types (targeted plane) from each volume, all planes, and cuboids within the whole volume (whole volume planes and cuboids). *** P ≤ 0.001 and **** P ≤ 0.0001. (C) Principal component feature space plot for the sample-level attention-aggregated volume features for the whole volume 3D approach, with the colors indicating ground truth labels. The good separation between the two classes supports the observed high AUC performance.(D) Kaplan-Meier survival analysis stratified at 50 percentile by TriPath-predicted risk on survival prediction phantom dataset ( n = 150 ) for 2D targeted single plane and whole volume 3D approaches. Error bars indicate one standard deviation from the mean, over five different experiments. Further details of the simulation dataset are described in the STAR Methods . 2 Figure S2: Additional metrics for high & low-risk patient classification task, related to Figure 2 , 3 Balanced accuracy and F1 score for high & low-risk patient classification task for (A) simulation (B) OTLS (C) microCT dataset. For OTLS and microCT cohorts, we also display performance metrics for the clinical baseline based on histologic examination of the prostatectomy specimen. In all three datasets and metrics, the 3D treatment of the whole volume and 3D patching is superior to the 2D plane-based alternatives. Error bars indicate one standard deviation from the mean, over five different experiments. 3 Figure S3: Comparison between different feature encoders in OTLS and microCT cohort, related to Figure 2 , 3 TriPath uses transfer learning to extract representative and compressed features from 2D patches and 3D patches. Consequently, the encoder architecture (e.g., CNN, ViT) and the pre-training dataset (e.g., radiology, videos, images) affect the downstream performance. (A) In the OTLS cohort, five encoders of different architectures and pretraining datasets were considered. Resnet-2D and SwinViT (2D) are defined by averaging 2D features of each level within each 3D patch. (B) The 2D slice with the largest tissue area and slices at ±20 μm from the level were selected for feature extraction with 2D encoders. (C) Even with a different feature encoder, we observe increased AUC as larger tissue volume is used for analyses. (D-F) The same set of experiments for the microCT cohort. Results demonstrate that different feature encoders and pretraining datasets lead to varying performance levels. CNNs and ViTs trained on natural images or videos lead to better performance than those pretrained on domain-specific datasets (radiology and histology). We attribute the low performance of the radiology-pretrained SwinViT (3D) to a large difference in image resolution (3D radiology: 1 ∼ 2 mm /voxel vs. 3D pathology: 1 ∼ 4 μm /voxel). This suggests that a feature encoder pretrained on datasets from the same data domain is warranted, which we leave for future work. For the main analyses, we use spatiotemporal CNN[ 108 ] ( i.e. , ResNet-(2+1)D) for 3D analysis as it provides consistently good performance across both cohorts. For 2D analysis, we use Resnet-2D since it shares the same residual network backbone as Resnet-(2+1)D, thereby allowing for a fair comparison between 2D and 3D tasks. Error bars indicate one standard deviation from the mean, over five different experiments. 4 Figure S4: Integrated gradient analysis for open-top light-sheet microscopy (OTLS) dataset, related to Figure 2 (A) Patches from the high IG cluster exhibit infiltrative carcinoma that resembles predominantly poorly-differentiated glands (Gleason pattern 4), exhibiting cribriform architecture. Patches from the middle IG cluster exhibit infiltrative carcinoma that resembles mixtures of Gleason patterns 3 and 4. Patches from the low IG cluster predominantly exhibit large, benign glands, with occasional corpora amylacaea. (B) Scatter plot of the normalized IG patch scores averaged within each sample as a function of predicted risk (the predicted probability for the high-risk group). (C) The scatter plot of the proportion of the number of high, middle, and low IG group patches in each sample as a function of predicted risk, which shows that a sample with a higher predicted risk profile has a larger (smaller) fraction of high (low) IG patches. (D) Kaplan-Meier curve for the cohort stratified (50%) by the ratio of the number of patches in the high and low IG group. The good stratification suggests that the extent to which prognostic morphologies manifest in each sample is also important. The scale bar is 100 μm . 5 Figure S5: Examples of integrated gradient (IG) heatmaps for open-top light-sheet microscopy (OTLS) cohort, related to Figure 2 The integrated gradient scores are assigned to each patch with high IG (low IG) patches indicating that patch contributes to an unfavorable (favorable) prognosis. (A) High IG areas in the high-risk sample contain cancerous glands that resemble poorly differentiated tumor glands (Gleason patterns 4). (B) In the low-risk sample, the high IG areas are those with cancerous glands that are smaller, more tortuous, and more closely resemble Gleason pattern 4, as well as regions with a cellular stroma. All scale bars are 200 μm . The heatmaps can also be visualized in our interactive demo. 6 Figure S6: Integrated gradient (IG) heatmaps for the microcomputed tomography (microCT) cohort, related to Figure 3 The integrated gradient scores are assigned to each patch with high IG (low IG) patch indicating that patch contributes to unfavorable (favorable) prognosis. (A) In this high-risk sample, high IG values are localized in areas with the smallest and densest cancerous glands, especially when they are in or adjacent to the capsule of the prostate, as well as dense stroma that resembles the prostate capsule. (B) Similar to the high-risk case, high IG regions in this low-risk sample correspond to areas with small, dense cancerous glands and dense stroma. The juxtaposition of these two morphologies has particularly high IG values. All scale bars are 500 μm . The heatmaps can also be visualized in our interactive demo. 7 Figure S7: Integrated gradient analysis for microcomputed tomography (microCT) dataset, related to Figure 3 (A) The high IG cluster consists of patches with infiltrative carcinoma that most closely resembles Gleason pattern 4; however, the lower resolution and lack of H&E staining make definitive grading infeasible by visual inspection of the microCT images alone. In the middle IG cluster, most patches contain infiltrating carcinoma that resembles Gleason patterns 3 and 4. The low IG cluster consists mostly of patches containing benign prostatic tissue with occasional foci of infiltrative carcinoma that resembles Gleason pattern 3. (B) Scatter plot of the normalized IG patch scores averaged within each sample as a function of predicted risk (the predicted probability for the high-risk group). (C) The scatter plot of the proportion of the number of high, middle, and low IG group patches in each sample as a function of predicted risk, which shows that the sample with higher predicted risk has a larger (smaller) fraction of high (low) IG patches. (D) Kaplan-Meier curve for the cohort stratified (50%) by the ratio of the number of patches in the high and low IG group. The good stratification suggests that the extent to which prognostic morphologies manifest in each sample are also important. The scale bar is 250 μm . 8 Figure S8: Clinical validation of 3D pathology for OTLS cohort, related to Figure 4 TriPath’s performance is further validated in OTLS cohort with the second part of the reader study. (A) The web interface for the reader study, where a pathologist could scroll through images of OTLS biopsy (three and all slices for the first and second round, respectively). (B) After 2 months of washout period from the first part, each pathologist was shown all slices of the OTLS biopsy image to provide per-biopsy diagnoses. (C) Cohort-level ( n = 50 ) BCR status prediction AUCs are shown based on 6 pathologists’ diagnoses of all image slices (individual and consensus) and the diagnosis from standard post-operative histopathology of the whole prostatectomy specimen. We employ two versions of TriPath: whole volume 2D slices (blue), where 2D patches are generated from the 2D slices across the whole volume, and whole volume 3D (orange), where 3D patches are generated from the whole volume (our 3D pathology baseline). (D) Quadratic weighted kappa to assess agreement between pair of pathologists. The median kappa value of 0.662 is slightly lower than that of 2D reader study (median kappa: 0.677). While the pathologists consensus performance increases compared to that of the first reader study ( All slices AUC: 0.799 vs. three slices AUC: 0.744), we observe that whole volume 3D TriPath still outperforms all clinical baselines. Combined with the fact that the median Kappa value did not change significantly, the results suggest that it is nontrivial for humans to process a huge stack of 2D slices (100-fold increase in number of slices per biopsy). In addition, that whole volume 3D TriPath outperforms both the pathologist baselines and whole volume 2D slices TriPath, both of which use the entire volume and rely on interpreting 2D morphology, suggests the importance of encoding 3D morphology. 9