GRASP: GRAph-Structured Pyramidal Whole Slide Image Representation

With recent progress in deep learning, deep features, i.e. high-level embeddings from a deep network, have advanced past handcrafted features and are considered the most robust sources for image representation. Pre-trained networks such as DenseNet Huang et al. (2017), ResNet He et al. (2016), or Swin Liu et al. (2021) draw their features from millions of non-medical and non-histopathological images, where they cannot necessarily produce high-level embeddings for complex images, especially rare cancers Wang et al. (2023); Riasatian et al. (2021); Ciga et al. (2022); Wang et al. (2022). In this context, the use of Variational Autoencoders (VAEs) has been evaluated in Chen et al. (2022a), where the authors show that DenseNet pre-trained on ImageNet performs better for extracting semantic features from WSIs than VAEs. However, domain-specific vision encoders such as KimiaNet Riasatian et al. (2021), CTransPath Wang et al. (2022), PLIP Huang et al. (2023), UNI Chen et al. (2023), Virchow Vorontsov et al. (2023), etc. were developed and trained on large sets of histopathology images (patches) outperforming pre-trained models on ImageNet across various tasks.

MIL Approaches: Several domains of deep learning have been explored in an attempt to effectively address the task of classification in digital pathology. Models such as AB-MIL Ilse et al. (2018), CLAM Lu et al. (2021), and Trans-MIL Shao et al. (2021) have utilized MIL with promising results. Such approaches have generally focused only on instance-level feature extraction and have not yet explored modeling global, long-range interactions within and across different magnifications. In Waqas et al. (2024), a detailed overview of different MIL methods has been provided.
Graph-based Approaches: To incorporate contextual information and long-range interactions, models such as PatchGCN Chen et al. (2021) and DGCN Zheng et al. (2021) have been designed with a graph structure that can capture and learn context-aware features from interactions across the WSI. These models represent WSIs as graphs where the nodes are usually embeddings and edges are defined based on clustering or neighborhood node similarity, which in turn adds new hyperparameters and increases inference time. The similarity between nodes can be measured in terms of spatial or latent space, leading to the construction of different graphs for each WSI.
Multi-Magnification Approaches: Multiple efforts have been made to incorporate multi-magnification information in the context of gigapixel histopathology subtyping tasks. Models such as HiGT Guo et al. (2023), ZoomMIL Thandiackal et al. (2022), CSMIL Deng et al. (2023), $H^{2}$-MIL Hou et al. (2022), and DSMIL Li et al. (2021) address this by aggregating contextual tissue information using features from multiple magnifications in WSIs. DSMIL concatenates embeddings from different magnifications together by duplicating lower-magnification features, making the model biased toward lower magnifications and unable to look into inter-magnification information. On the other hand, ZoomMIL aggregates information from $5x$ to $20x$ in a fixed hierarchy with no interaction in the opposite direction. Chen et. al explore this in the context of vision transformers with their Hierarchical Image Pyramid Transformer (HIPT) Chen et al. (2022b). Their architecture incorporates regions of size $256\times 256$ and $4096\times 4096$ pixels to leverage the natural hierarchical structure of WSIs. $H^{2}$-MIL also adopts a graph-based approach, where it pools the nodes in each magnification using an Iterative Hierarchical Pooling module. Our proposed model, on the other hand, is designed to dynamically aggregate information within and across different magnifications without using traditional pooling layers in its intra-magnification interactions. It also employs a similar mechanism to zoom-in and zoom-out through its inter-magnification interactions, from lower to higher magnifications and vice versa.

Contrary to Multiple Instance Learning (MIL) approaches, which use a bag of instances to represent a given WSI, GRASP benefits from a graph-based, multi-magnification structure to objectively represent connections between different instances across and within different magnifications. To build a graph and learn a graph-based function $\mathcal{F}$ that predicts slide-level labels with no knowledge of patch labels, the following formulation is adopted.

For a given WSI, $W_{r}\in\mathbb{R}^{N\times M\times 3}$ with label $\mathcal{Y}$, three sets of $m$ patches, $\{p_{i}\in\mathbb{R}^{n\times n\times 3}:\forall i\in[1,...,m]\}$, $\{p^{\prime}_{i}\in\mathbb{R}^{n\times n\times 3}:\forall i\in[1,...,m]\}$, and $\{p^{\prime\prime}_{i}\in\mathbb{R}^{n\times n\times 3}:\forall i\in[1,...,m]\}$ are extracted for each magnification of $\textbf{M}_{1}=5x$, $\textbf{M}_{2}=10x$, and $\textbf{M}_{3}=20x$, respectively. It is important to note that $p^{\prime\prime}_{i}$ is the high-resolution window located at the center of $p^{\prime}_{i}$, and $p^{\prime}_{i}$ is the high-resolution window located at the center of $p_{i}$.  These patches provide $3m$ patches in total that are then fed into an encoder $\phi$ to encode extracted patches into a lower dimension space as follows:

where $h_{i}$ is the feature vector corresponding to the patch $p_{i}$. Correspondingly, $h^{\prime}_{i}$ represents $p^{\prime}_{i}$, and so $h^{\prime\prime}_{i}$ does $p^{\prime\prime}_{i}$. Using all the feature vectors for each $W_{r}$, graph $\mathbb{G}_{r}$ is constructed using the transformation $\Gamma$:

Eventually, classifier $\mathcal{C}$ is applied on top of graph convolutional layers $\mathcal{G}$ to build the graph-based function $\mathcal{F}$ to predict slide-level label $\mathcal{Y}$ as follows:

We start by extracting multi-magnification patches as described earlier. Then, for any $i$, we use the same encoder to encode $p_{i}$, $p^{\prime}_{i}$, and $p^{\prime\prime}_{i}$ into features $h_{i}$, $h^{\prime}_{i}$, and $h^{\prime\prime}_{i}$ respectively. Having instances features, we use the transformation $\Gamma$ to build $\mathbb{G}_{r}$ as introduced in Eq. 2.

The mechanism of connecting every two nodes in $\mathbb{G}_{r}$ through $\Gamma$ is premised upon an intuition of the pyramidal nature of WSIs as well as the way in which a conventional light microscope works when one switches from one magnification to another. When using a microscope, increasing magnification preserves the size of the image yet increases resolution by showing the central window of the lower magnification. This is the exact procedure we use to extract our patches in three magnifications. Therefore, for any $i$, $h_{i}$, $h^{\prime}_{i}$, and $h^{\prime\prime}_{i}$ are connected to each other via undirected edges, where this connection represents inter-magnification information contained in the features. On the other hand, for any $i$, all $h_{i}$’s contain information in $\textbf{M}_{1}$, such that they are connected to each other, forming a fully connected graph at $\textbf{M}_{1}$ magnification to represent intra-magnification information. Similarly, all $h^{\prime}_{i}$s are connected to each other and also all the $h^{\prime\prime}_{i}$s to represent intra-magnification information contained in $\textbf{M}_{2}$ and $\textbf{M}_{3}$, respectively.

Figure 2 shows a small example of such a graph for $\mathbb{G}_{r}=(V_{r},E_{r})|_{m=4}$ where blue, red, and green nodes each form a fully connected graph of size $m$; the inter-magnification relationship can also be seen via the edges between blue & red nodes as well as the red & green nodes. So far, each WSI, $W_{r}$, has been represented by a fixed graph $\mathbb{G}_{r}$ with $3m$ nodes and $\frac{(3m+1)m}{2}$ edges. These graphs are thus deployed to train the GCNs and predict the label $\mathcal{Y}$ at the output.

Following Eq. 3, we are defining $\mathcal{G}$ which includes three GCN layers. The intuition behind using three layers is that as a pathologist begins to look for a tumor in a given WSI, they use an initial magnification to find the region of interest; Once found, they consult with other magnifications, which may require zooming in and out back and forth, to confirm their final decision. Therefore, as shown in Figure 2, all nodes in the graph interact with one another in a hierarchical fashion through the GCN layers. Consequently, each node gradually gathers information from all other nodes, and therefore, if there are any important messages carried by some nodes, it is guaranteed to be broadcast to all other nodes, which is the equivalent of zoom-in and zoom-out mechanism. This dynamic and hierarchical structure imposes theoretical properties on the model which is going to be discussed here. Following the graph convolutional layer introduced in Kipf & Welling (2016), the graph nodes are updated as follows:

where $b^{(l)}$ is bias; $h_{i}^{(l+1)}$ is the node feature’s update of the graph at $(l+1)$-th step at $\textbf{M}_{1}$; $\mathcal{N}(i)$ is the set of neighbors of node $i$, $c_{ji}=\sqrt{|\mathcal{N}(j)||\mathcal{N}(i)|}$, where given the symmetry of the graph, all $c_{ji}$s are equal; and $\alpha(.)$ is the activation function, which is ReLU in our implementation. $h_{i}^{\prime(l+1)}$ and $h_{i}^{\prime\prime(l+1)}$ expressions follow the same logic as $h_{i}^{(l+1)}$ in terms of the parameters mentioned above.  After the last graph convolutional layer, where the intra-magnification convergence happens, the graph is passed through an average readout module to pool the three-node graph mean embedding and then is fed into the two-layer classifier $\mathcal{C}$ to predict $\mathcal{Y}$.

Based on the idea of capturing the information within and across magnifications, we now show that node features converge in each magnification in the graph to one node Having this, GRASP essentially encodes a graph of size $3m$ nodes to only 3 nodes. We interpret this as the model learning the information contained in the magnification through interaction with other magnifications without a need for traditional pooling layers.

Proof: Please see Section 7.3 (Theoretical Analysis).

Note: It is worth mentioning that the assumptions made in Theorem 1 are minimal. In the case of $L_{2}$-bounded weights, it has been further explained in Wu et al. (2023); Cai & Wang (2020). The assumption that graph node features at $l=0$ are bounded is also minimal as we use a frozen features extractor, which ideally should not generate nondefinitive values for a finite feature vector.
With that, we show

which essentially means that $\left\|h_{i}^{(3)}-h_{j}^{(3)}\right\|_{2}$ is upper bounded inversely with the reciprocal of $(m+1)^{3}$ . Thus, as $m$ increases, the upper-bound gets tighter and eventually leads to $\lim_{m\to\infty}||h_{i}^{(3)}-h_{j}^{(3)}||_{2}=0$. As a result, we can conclude the following corollary.

This is necessarily equivalent to pooling the nodes in magnification level, yet with a completely new approach than the traditional pooling layer, without further imposing computational load on the network for pooling. Taking into account the fact that $h^{*}$, $h^{\prime*}$, and $h^{\prime\prime*}$ are not necessarily equal, our model is fusing node features in each magnification while it consults with other magnifications, and draw the conclusion via averaging nodes across three magnifications at the end of the convolutional layers by means of the readout module. This means that the final embedding of the graph is $\dfrac{h^{*}+h^{\prime*}+h^{\prime\prime*}}{3}$. We believe that this process helps the model reduce variance and uncertainty in making predictions as $m$ grows. To support this claim, we provide empirical evidence detailed in Section 7.5.1 (Monte Carlo Test).

The structure of the graph has been designed in such a way that it does not get stuck in the bottleneck of over-smoothing, a common issue in deep GCNs Cai & Wang (2020). Our intuition is that nodes in $\textbf{M}_{1}$, $\textbf{M}_{2}$, and $\textbf{M}_{3}$ are interacting via message passing, and the flow of inter-magnification information helps the model to keep its balance and continue the process of learning. Nevertheless, by increasing the number of GCN layers for the model to four or higher, the so-called over-smoothing problem will take place, which can possibly deteriorate the model’s performance. On the other hand, less than three layers of GCNs might not be able to fully capture the inter-magnification interactions. This leaves us with three layers of GCNs, which is equal to the graph’s diameter. In addition to this theoretical description, we empirically support our claim in Section 7.5.5 (Graph Depth).

We utilize three datasets: Esophageal Carcinoma (ESCA) from The Cancer Genome Atlas (TCGA), which includes 135 WSIs across two subtypes, and Ovarian Carcinoma and Bladder Cancer, where Ovarian Carcinoma consists of 948 WSIs with five histotypes, while Bladder Cancer contains 262 WSIs with two histotypes. These datasets were curated using HistoQC Janowczyk et al. (2019). A detailed breakdown of each dataset is available in Table 3.

To compare with state-of-the-art approaches, we repeated these experiments with the same cross-validation folds and random seeds to have a fair comparison; the choice of ten random seeds is to capture statistical significance and reliability. For evaluating the models, we adopt Balanced Accuracy and F1 Score since these metrics show how reliable a model performs on imbalanced data, and more importantly on clinical applications. We compare our proposed model, GRASP, with models using different approaches to have a broad spectrum of evaluation. These models include Ab-MIL Ilse et al. (2018), Trans-MIL Shao et al. (2021), CLAM-SB Lu et al. (2021), and CLAM-MB Lu et al. (2021) from the attention/transformer-based family; ZoomMIL (2021) Thandiackal et al. (2022), H2MIL (2022) Hou et al. (2022), and HiGT (2023) Guo et al. (2023) from multi-magnification approaches since they have a hierarchical structure and are compatible with our patch extraction paradigm; and PatchGCN: latent & spatial Chen et al. (2021) and DGCN: latent & spatial Zheng et al. (2021) from graph-based learning approaches.

Comparing Tables 1 and 2, all the models performed better with KimiaNet embeddings than with Swin embeddings, which is mostly because KimiaNet has domain knowledge and can provide more contextual features than Swin. Furthermore, GRASP, H2MIL, ZoomMIL, Ab-MIL, and CLAMs showcase robust generalization and effective performance even when utilizing features from different backbones, especially with GRASP being the most robust model.

Although Deng et al. (2023) has used attention score distribution to show their model is reliable across different magnifications, we want to step further and adopt a similar logic as first introduced in Selvaraju et al. (2017) to define the concept of energy of gradients for graph nodes (Please see Graph-Based Visualization 7.6 in Appendix). Therefore, for the first time in the field, we show that an AI model such as GRASP can learn the concept of magnification and behave according to the subtype and slide characteristics. To this end, we formulate an experiment to obtain a sense of each magnification’s influence on the model which leads to Figure 3. The main takeaway of this experiment is that depending on the subtype, the distribution of referenced magnifications by GRASP is different. From a pathological point of view, this finding fits our knowledge of the biological properties of each subtype. As an example, we conducted a case study on the bladder dataset, where the micropapillary subtype is known to be diagnosed generally in lower magnification owing to its morphological properties and the structure of micropapillary tumors, whereas UCC needs to be examined in higher magnifications due to its cell- and texture-dependent structure.

On the Ovarian dataset, for the subtype ENOC, endometrioids can often be recognized and identified at low power as they tend to have characteristic glandular architecture occupying contiguous, large areas. Low-grade serous carcinomas (LGSC) can be very difficult at low power due to the necessity of confirming low-grade cytology at high power. For CCOC, clear cell carcinomas have characteristic low-power architectural patterns but can also require high-power examinations to exclude high-grade serous carcinoma with clear cell features, meaning that important information is distributed on all magnifications. The other subtypes, MUC and HGSC, may either show pathognomic architectural features at low power or require high-power examination on a case-by-case basis. According to Figure 3, GRASP collects the information from all three magnifications for MUC and CCOC.

Furthermore, to examine whether GRASP understands the biological meaning of the data, i.e., differentiating between tumor vs non-tumor regions, we conduct a visualization experiment to plot the pixel-level heatmap of patches in multiple magnifications as depicted in Figure 4.

Here, we design five experiments to evaluate our proposed model. Firstly, a Monte Carlo test on graph size, i.e., the number of nodes to investigate the impact of the number of nodes on the model’s performance (Section Monte Carlo Test 7.5.1). Secondly, analyzing model performance on individual or pairs of magnifications to study the effectiveness of multi-magnification representation (Section Magnification Test 7.5.2). Thirdly, we investigate the effect of different graph convolution types (Section Graph Convolutions 7.5.3). In the fourth experiment, we study how different models perform when all the patches from a WSI are used (Section Patch Number 7.5.4). Lastly, we empirically show that the graph depth of $d=3$ is the appropriate choice for our design (Section Graph Depth 7.5.5).

A total of 1,133,388 patches of size $1000\times 1000$ pixels for the Ovarian dataset, 602,874 patches of size $224\times 224$ from the ESCA dataset, and 313,191 patches of size $1000\times 1000$ pixels for the Bladder dataset are extracted in multi-magnification setting (approximately 2 TB of Gigapixel WSIs). Patches being extracted such that they do not overlap at $\textbf{M}_{3}$ while overlapping at $\textbf{M}_{2}$ and $\textbf{M}_{1}$ is inevitable. From each magnification, $m\leq 400$ patches (note that we use a large field of view, meaning this number eventually covers much of tissue regions) have been extracted per slide in both Ovarian and Bladder datasets, as it’s been shown in Chen et al. (2022a); Wang et al. (2023); Rasoolijaberi et al. (2022) that a subset of patches is enough to represent WSIs.
For patch-level feature extraction, we utilized two backbones: KimiaNet and Swin_base. Given that KimiaNet was trained on TCGA data in a supervised fashion, we intentionally refrained from extracting features from the ESCA dataset using this backbone to ensure an unbiased and leakage-free comparison. Conversely, we employed Swin, pre-trained on ImageNet, to extract features from all three datasets.
For each cancer dataset, we trained our proposed method in a 3-fold cross-validation and repeated the experiments ten times with different random seeds, where random seeds were randomly generated, to ensure a rigorous comparison. In order to prevent data leakage in our cross-validation splits, we split the slides based on their patients, since some patients have more than one slide, meaning that slides were split in a way that all slides from the same patient remain in the same set.

To tackle the data imbalance problem, for all models in the study, we deployed a weighted cross entropy loss. A learning rate of $0.001$ and a weight decay of $0.01$ for Adam optimizer have been adopted, and in case competing models were not converging, learning rate of $0.0001$ resolved the problem. Models were trained for 100, 50, and 10 epochs for the Ovarian, ESCA, and Bladder datasets, respectively. Specific to GRASP, the first two layers are of size 256 and the last layer output is of size 128. For all training and testing, the GPU hardware used was either a GeForce GTX 3090 Ti-24 GB (Nvidia), Quadro RTX 5000-16 GB (Nvidia), RTX 6000-48 GB (Nvidia) based on availability. Deep Graph Library (DGL), PyTorch, NumPy, SciPy, PyGeometric, and Scikit-Learn libraries have been used to perform the experiments.

Here, we prove Theorem 1 for any $h_{i}^{(3)}$ and $h_{j}^{(3)}$, and conclude the case for $h_{i}^{\prime(3)}$s and $h_{i}^{\prime\prime(3)}$s similarly. To start with, we demonstrate Lemma 1.

Description: given Eq. 17, every two arbitrary nodes $h_{i}^{(3)}$ and $h_{j}^{(3)}$ in one level of magnification are converging to each other. This means that all nodes are converging to the same value, which we name $h^{*}$. Thus, $\forall i\in[1,...,m]$, $\lim_{m\to\infty}\left\|h_{i}^{(3)}-h^{*}\right\|_{2}=0$ or equivalently $h_{i}^{(3)}\longrightarrow h^{*}$. Using the same logic as above, one can conclude $h_{i}^{\prime(3)}\longrightarrow h^{\prime*}$ and $h_{i}^{\prime\prime(3)}\longrightarrow h^{\prime\prime*}$. Since each of $h^{*}$, $h^{\prime*}$, and $h^{\prime\prime*}$ are a function of $m$, increasing $m$ would result in them being a better estimation/representation for the intra-magnification information.

Let’s call the output of the classifier, $S$, where the logit for correctly classified slides/graphs is $S_{c}$. To visualize the importance of magnifications, we compute the magnitude of the gradient of a graph with respect to its node features at $l=0$ ($h_{i}^{(0)}$; for the sake of brevity, we drop the superscript $(0)$ form $h_{i}^{(0)}$ and show it as $h_{i}$), which is $\left|\frac{\partial S_{c}}{\partial h_{i}}\right|$ for the magnification $\textbf{M}_{1}$. $\left|\frac{\partial S_{c}}{\partial h_{i}}\right|$ is a vector of size $1024\times 1$, and we have $m$ nodes giving result to $m$ such vectors. Likewise, we can define $\left|\frac{\partial S_{c}}{\partial h_{i}^{\prime}}\right|$ and $\left|\frac{\partial S_{c}}{\partial h_{i}^{\prime\prime}}\right|$ for $\textbf{M}_{2}$ and $\textbf{M}_{3}$, respectively. Arranging these absolute gradients for each magnification in a matrix of size $m\times 1024$ as follows,

As such, putting matrices in 20, 21, and 22 together in a matrix gives us the overall heatmap for the graph, of the size $3m\times 1024$, to compare the influence of each magnification:

This is the heatmap depicted in Figure 9 as a graph-based heatmap, which shows how model focuses on different magnifications.

Having the gradient for each node in the graph, we develop the concept of energy of gradients to find out which magnification(s) play a more important role in GRASP’s final decision. To do so, we start by defining $\mathcal{E}_{\mathbf{M}_{1}}$ as follows,

similarly, for $\mathbf{M}_{2}$ and $\mathbf{M}_{3}$,

Having these energies, the energy contribution of each magnification is calculated based on their relative share in the whole energy spent in the graph:

Accordingly, the importance of each magnification can be quantified for further investigations. More samples are provides in Figure 3.