Differentiation of Pancreatic Ductal Adenocarcinoma and Chronic Pancreatitis using Graph Neural Networks on Histopathology and Collagen Fiber Features

doi:10.21203/rs.3.rs-1951132/v1

Download PDF

Article

Differentiation of Pancreatic Ductal Adenocarcinoma and Chronic Pancreatitis using Graph Neural Networks on Histopathology and Collagen Fiber Features

https://doi.org/10.21203/rs.3.rs-1951132/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Pancreatic ductal adenocarcinoma (PDAC) is one of the most lethal human cancers. However, the symptoms and follow-up radiographic and histopathology assessments of PDAC are similar to chronic pancreatitis (CP) and can lead to misdiagnosis. The need for accurate differentiation of PDAC and CP has become a major topic in pancreatic pathology. These two diseases can present similar histomorphological features, such as excessive deposition of fibrotic stroma in the tissue microenvironment and inflammatory cell infiltration. In this paper, we present a quantitative analysis pipeline empowered by graph neural networks (GNN) capable of automatic detection and differentiation of PDAC and CP in human histological specimens. Modeling histological images as graphs and deploying graph convolutions can enable the capture of histomorphological features at different scales, ranging from nuclear size to the organization of ducts. The analysis pipeline combines image features computed from co-registered hematoxylin and eosin (H&E) images and Second-Harmonic Generation (SHG) microscopy images, with the SHG images enabling the extraction of collagen fiber morphological features. Evaluating the analysis pipeline on a human tissue micro-array dataset consisting of 786 cores and a tissue region dataset consisting of 268 images, it attained 86.4% accuracy with an average area under the curve (AUC) of 0.954 and 88.9% accuracy with an average AUC of 0.957, respectively. Moreover, incorporating topological features of collagen fibers computed from SHG images into the model further increases the classification accuracy on the tissue region dataset to 91.3% with an average AUC of 0.962, suggesting that collagen characteristics are diagnostic features in PDAC and CP detection and differentiation.

Pancreatic cancer is one of the most lethal malignancies with a poor 5-year survival rate of around 5–9% that has remained almost stagnant since the 1960s [1, 2]. Among pancreatic cancers, more than 85% of the cases are pancreatic ductal adenocarcinoma (PDAC), which mostly originates from pancreatic ductal epithelium in the pancreatic head [2]. Chronic pancreatitis (CP), though an entirely different disease with a distinct prognosis compared to PDAC, can present morphological features similar to PDAC ranging from radiographical level to histopathological level [3–7]. Therefore, differentiating PDAC and CP is a challenge in pathology and the consequence of misdiagnosis can be severe due to the rapid progression of PDAC and the high frequency of distant metastases [8–10]. Moreover, multiple attempts at diagnosis with fine-needle aspiration or needle biopsy techniques are frequently required due to the small amounts of tissue recovered with these procedures.

In clinics, computed tomography and ultrasound are often used to detect irregular masses in the pancreas, but histopathological assessment based on tissue sections is required to ensure the accurate differential diagnosis of PDAC and CP [3–5, 11, 12]. CP can present similar histomorphological features such as inflammatory infiltration, dense stroma, angulated glands, cytologic atypia, and tumor-like duct organization that mimic PDAC. The difficulty in distinguishing PDAC from CP is further aggravated by the fact that PDAC may induce CP in the surrounding pancreatic tissue, posing additional challenges to tissue sampling and diagnosis [3–6].

One important hallmark of the tissue microenvironment for both PDAC and CP is a dense desmoplastic stroma, characterized by the increased deposition of fibrillar collagen. In PDAC, the extensive desmoplasia present in the tumor has been shown to have a negative impact on stromal vascularization and immune cell infiltration, thus, it may contribute to resistance to radiotherapy and hinder drug delivery [13–15]. Besides PDAC, evidence in a variety of cancer types has also shown that stromal properties are important factors for disease diagnosis, cancer progression, and tumor response to therapy [16–21]. In differentiating PDAC and CP, the collagen fiber topology at the stromal-epithelial interfaces has been shown to be a statistically significant discriminating feature [22]. In PDAC, ordered periductal stroma characterized by highly aligned and elongated collagen fibers is a negative prognostic factor [18]. Quantification of the size, shape, and patterns of collagen fibers show that these factors might impact tissue stiffness and are associated with an increased risk of cancer progression [23]. Many of the findings are derived from studies using label-free collagen-sensitive imaging modalities such as second-harmonic generation microscopy (SHG) and polarization-based optical microscopy, which can enable the quantification of collagen fibers and desmoplasia [24–28]. With SHG, new insights into tissue properties associated with stroma have been obtained to inform biological understanding, diagnosis, prognostication, drug development, and therapeutic innovations [13–15, 29, 30].

In the last 30 years, the advent of digital slide scanners and whole slide imaging (WSI) has made digital pathology a fruitful research field for image feature analysis with strong potential for clinical applications [31–34]. The potential for WSI-based computational analysis is further increasing due to the recent success of deep neural networks in computer vision and language modeling [35, 36]. The ability to digitize routinely processed tissue on glass slides paved the way for computer algorithms to quantify histological image features. Now, deep learning has demonstrated state-of-the-art performance when applied to many computational histopathology analysis tasks such as cancer detection, tissue segmentation, disease prognostication, and spatial omics analysis [37–39].

In this paper, we present a computational pipeline for the detection and differentiation of histological samples of PDAC and CP, empowered by advances in artificial intelligence and collagen-targeted tissue imaging. The design of our analysis methods is motivated by the observation that discriminative histological features of PDAC and CP span multiple fields and multiple magnifications, from cell morphology to duct organization, and both diseases are characterized by dense stroma with potentially different collagen topography. Thus, a model with the ability to learn image features across different scales and incorporate collagen-based image features would be favored.

The analysis pipeline includes a deep learning model built on graph neural networks (GNNs) [40, 41] that can be trained with manually classified, coarsely annotated regions of varying sizes, along with a region proposal algorithm that generates candidate regions in unannotated slides. Built on the canonical scheme of WSI data processing, the proposed method extracts local features from image patches but further models the patch-to-patch interactions by constructing graphs from the patch sets and utilizing graph convolutions [41, 42]. GNNs increase the expressivity of the model by letting the information flow between adjacent image patches, thus, could capture histomorphological features that span multiple patches [43–45]. This more closely mimics histopathological examinations conducted by pathologists where relationships between various tissue features are frequently characterized. The concept of using GNNs aligns well with the fact that the major diagnostic criteria of PDAC and CP include features ranging from the cellular level (e.g. variation of epithelial nuclei sizes) to the tissue level (e.g. spatial organization of ducts in comparison to nerves and arteries).

The proposed GNN-based method is evaluated on two types of datasets, one collected from human tissue micro-arrays (TMAs) and the other from human tissue sections, both consisting of PDAC, CP, and normal pancreas tissue samples. The resulting model outperforms the widely used multiple-instance learning (MIL) framework [46, 47], achieving 86.4% accuracy with an average area under the curve (AUC) of 0.954 on the TMA dataset and 88.9% accuracy with an average AUC of 0.957 on the tissue section dataset. Furthermore, we demonstrate that the incorporation of collagen-based features extracted from SHG images leads to higher classification accuracy compared to using brightfield hematoxylin and eosin (H&E) features alone, with classification accuracy increasing from 88.9–91.3%. This result confirms the diagnostic potential of utilizing characteristic collagen topology in differentiating PDAC, CP, and normal pancreas tissue.

2.1 Sample sets and pathology annotations

Human tissue microarrays. 7 H&E stained formalin-fixed paraffin-embedded human pancreas TMAs were purchased from US Biomax (Rockville, MD, USA), the slide IDs used were: BIC14011b, HPan-Ade180Sur-01, HPanAde170Sur-01, HPan-Ade120Sur-01, PA2081b, BBS14011, PA485. The TMAs consist of cores from 1.0 to 1.5 mm in diameter and were reviewed by two US Biomax pathologists. All slides were scanned with an Aperio CS2 (Leica Biosystems, Wetzlar, Germany) at 20× magnification. All PDAC (380 cores), chronic pancreatitis (193 cores), and normal (213 cores) TMA cores from the slides were used for analysis and image classification.

Human tissue sections. 24 H&E stained tissue section slides were obtained from Metro Health Medical Center (Cleveland, Ohio, USA). The tissue sections were obtained from surgical specimens including Whipple procedures and distal pancreatectomies, and from autopsies. The slide set contains tissue sections from 12 cases of PDAC and 12 cases of chronic pancreatitis. Slides were reviewed and annotated by two pathologists (Omid Savari and Agnes G. Loeffler). For the PDAC slides, roughly 6 regions containing PDAC, 6 regions containing chronic pancreatitis, and 3 regions containing histologically normal pancreas were coarsely marked with rectangles; for chronic pancreatitis slides, on average 3 regions containing chronic pancreatitis and 6 regions containing normal-appearing pancreas were coarsely marked with rectangles. Note that since the regions were coarsely annotated, a PDAC region might include chronic pancreatitis sub-regions, and a PDAC or a chronic pancreatitis region might contain sub-regions with normal tissues. The annotations were made using Aperio ImageScope (12.4.3) software (Leica Biosystems, Wetzlar, Germany) on brightfield images scanned by an Aperio CS2 digital pathology slide scanner (Leica Biosystems, Wetzlar, Germany) at 20× magnification. The annotation process produced a total of 70 PDAC regions, 112 chronic pancreatitis regions, and 86 normal regions.

2.2 Second-harmonic generation microscopy imaging of samples

Annotated regions from the tissue section slides were imaged using a custom-built multiphoton laser-scanning microscope (LSM) that has been previously described [48]. In brief, a Tsunami Ti:Sapphire laser (Spectra-Physics, Santa Clara, CA, USA) tuned to 800 nm, with a pulse length of approximately 100 fs was focused on the sample by a 20X/0.75NA air immersion objective (Nikon, Melville, NY, USA). The SHG signal was collected in the forward direction, filtered by a narrow bandpass filter centered at 400 ± 10 nm (ThorLabs, Newton, NJ, USA), onto a H7422-40P photomultiplier tube (Hamamatsu, Hamamatsu, Japan). The scanning process and data acquisition signal were controlled by our custom laser-scanning software, OpenScan. The acquisition workflow utilizes the open-source software QuPath [49], Micro-Manager, and Pycro-Manager [50, 51].

2.3 Graph neural networks

2.3.1 Graph construction

Image registration, pre-processing, and normalization. Brightfield images of individual TMA cores and annotated regions were cropped from the Aperio scans. Brightfield images of annotated regions in the tissue section set were registered with the corresponding SHG images using the automatic H&E to SHG registration algorithm described in [52]. Image pairs that failed to be registered by the algorithm were further registered using the method described in [53] or manually registered by corresponding landmarks in Fiji [54].

Each brightfield image was tiled into patches with a size of 256×256 pixels (without overlap); low-entropy patches (< 0.05) that corresponded to empty backgrounds or hollow intra-gland areas were discarded. The tiles were then normalized using the H&E image normalization algorithm described in [55] with batch processing. Every batch contained 256 patches sampled from all patches, the normalization algorithm then treated the 256 patches as a single input image and performed normalization. Running the algorithm in batch helped stabilize the normalization by evenly sampling the whole dataset and avoiding the generation of extreme pixel values in a single image patch that might contain insufficient pixels for meaningful statistics.

Registered SHG images of annotated regions in the tissue section set were normalized by stretching the pixel values between 5 − 99.5 percentile to the full dynamic range. The images were then tiled into patches with a size of 256 × 256 pixels (without overlap), resulting in SHG image patches with pixel-wise correspondences to the brightfield image patches.

Patch feature extraction. A feature vector with a shape of 512×1 was then extracted from every brightfield patch using the ImageNet pretrained ResNet18 [56]. For each SHG patch, collagen fiber-specific features such as width, length, and angle were extracted using CurveAlign and CT-FIRE [57]. Fiber segments, alignment coefficients, numbers of fibers, circular mean and variance of fiber angles, mean and variance of fiber lengths, mean and variance of fiber straightness, and mean and variance of fiber widths were finally computed and formed a feature vector with a shape of 10×1 for each patch. The features were normalized to a range from 0∼1.

Constructing graphs from images. The image of a TMA core or an annotated region was treated as a single graph, and the patches tiled from the image were treated as nodes of the graph. Each patch (node) was connected with the nearest neighboring patches (if they existed) with edge weights equal to the cosine similarity between the feature vectors of the two nodes. The information about a tile (a node feature) consisted of the corresponding 512×1 brightfield feature vector in the case of the brightfield graph or the corresponding 10×1 SHG feature vector in the case of the SHG graph. A brightfield-SHG graph was constructed by concatenating the brightfield and SHG feature vectors as the node features (with a vector size of 522×1).

2.3.2 Graph convolutional network architecture

Feature pre-processing layer. The node features were first processed by a fully connected neural network layer followed by a parameterized ReLU (PReLU) [58]. For the brightfield graph and brightfield-SHG graph, the fully connected layer produced an output vector with a size of 512×1. For the SHG graph, the fully connected layer produced an output vector with a size of 32×1. The processed features were then fed to the following graph convolution layers.

Graph convolution network (GCN). Let V be the set of nodes and ${\mathbf{h}}_{v}$ where $v \in \mathcal{V}$ be the node feature vector processed using the procedures described above. The updated note feature ${\mathbf{h}}_{v}^{{\prime }}$ after a 1-hop graph convolution operation [41] is given by

$${\mathbf{h}}_{v}^{{\prime }}=\text{P}\text{R}\text{e}\text{L}\text{U}\left(\frac{1}{‖\mathcal{N}\left(v\right)‖}\sum _{u\in \mathcal{N}\left(v\right)}{\Theta }{\mathbf{h}}_{u}^{ }+\mathbf{b}\right), \forall v\in \mathcal{V}$$

where ${\Theta }$ is a weight matrix and $\text{b}$ is a bias vector, and $u \in \mathcal{N}\left(v\right)$ denotes the neighbors of node $v$ (including $v$ itself). Three 1-hop graph convolution layers are stacked and the update ${\mathbf{h}}_{v}^{{\prime }}$ of each layer is passed to the next graph convolution layer. The resulting embedding ${\mathbf{h}}_{o}$ then represents information up to 3-hop neighbors from the node. Besides, another type of graph convolution operation was evaluated in the design. This alternative model used a single-layer 3-hop simplified graph convolution (SGC) operator that directly aggregates information in the 3-hop neighborhood of the nodes without non-linearity between each hop [42]. This operation was more memory and computationally efficient and empirically showed comparable performance to the standard graph convolution. ${\mathbf{h}}_{o}$ was concatenated with ${\mathbf{h}}_{v}$ to form a new node embedding as output.

Readout layer. The output node embeddings of each graph were averaged to produce a graph embedding. The graph embedding was scored by a linear classification head followed by a sigmoid function to produce a graph prediction with the same shape as the graph label representation using one-hot encoding. The overall analysis pipeline including the network training procedure is shown in Fig. 1. The objective function with respect to a dataset containing $N$ image-label pairs can be written as:

$$\{\theta ,\mathbf{W}\}=\text{a}\text{r}\text{g}\underset{\theta ,\mathbf{W}}{\text{min}}\frac{1}{N}\sum _{i=0}^{N-1}\mathcal{L}(\mathbf{W}f\left({\text{G}}_{i};\theta \right),{\mathbf{y}}_{i})$$

where L is the loss function, θ is the set of weights parameterizing the graph convolution network f described above, W is the weight matrix of a linear layer, and ${\text{G}}_{i}$ and ${\mathbf{y}}_{i}$ are a graph constructed from an image and its class label, respectively.

Class coding and loss function. PDAC can induce CP in surrounding areas and areas with a mixture of the two are also often observed. However, when both CP features and PDAC features are present in the same region, the sample will be labeled as PDAC. Therefore, we used distributed coding instead of one-hot coding for class labeling to account for the potential co-existence of both classes in a sample. Specifically, a normal region or core is labeled as [0,0], a CP region or core is labeled as [0,1], and a PDAC region or core is labeled as [1, 1]. The final class label is obtained using the following rule: a prediction is regarded as normal when both digits are 0, CP when the first digit is 0 and the second digit is 1, and PDAC when the first digit is 1 (regardless of the second digit). The loss function is multi-label binary cross-entropy loss. The optimal thresholds for the two digits are obtained by finding the threshold leading to the best true-positive to false-positive ratio in the validation set for each digit.

2.3.3 Network training specifics

The networks were optimized with stochastic gradient descent with the Adam solver and the initial learning rate was set to 0.0001 [59]. The learning rate was scheduled with cosine annealing without warm restart [60]. The batch size was 1. For the TMA dataset, the cores were randomly separated into 4 equal-size chunks and a 4-fold cross validation was carried out. For the tissue section dataset, the slides were first randomly separated into 4 equal-size sets, and the annotated areas were then cropped from the slides to form 4 sets for a 4-fold cross validation. Cross-validations were repeated 3 times for each dataset, each time with a new random split.

2.4 Region proposal in whole slide images

A region proposal pipeline was used to identify candidate regions in WSIs of tissue sections. This pipeline took coarsely annotated target regions as inputs and searched for similar regions in WSIs based on the distribution of patches in the annotated regions. The overall process can be seen as a task of slide region retrieval for some annotated regions as queries.

Firstly, patches with a size of 256×256 were extracted from the training WSIs of the tissue section set without overlap, and the patch feature vectors were extracted by a ResNet18 pre-trained on ImageNet. A k-means clusterer $\mathcal{K}$ was trained on the feature vectors using mini-batch k-means [61]. Secondly, the same patch extraction and feature extraction processes were applied to the annotated regions of PDAC and chronic pancreatitis, and the feature vectors of the patches were assigned to the closest centroid of $\mathcal{K}$. The distribution of centroid assignments was computed and represented in a vector $\mathbf{q}$. The idea of learning a dictionary of pattern basis is frequently used in fine-grained image classification [62]. Thirdly, a sliding window of size 8 patches × 8 patches was used to scan through the WSIs, and the distribution of patch feature vector assignments regarding the centroids of $\mathcal{K}$ within the window was computed per sliding step. The distribution of the assignments was denoted as ${\mathbf{p}}_{ij}$, where $i,j$ are the coordinate of the window center. The cosine similarity between ${\mathbf{p}}_{ij}$ and $\mathbf{q}$ was computed and placed at $(i,j)$. A similarity map was produced after the window slides through the entire WSI. Lastly, a multi-threshold Otsu’s method was applied to the similarity map to separate the similarity values into 3 groups corresponding to background, low-similarity, and high-similarity areas [63]. The high similarity regions were processed by a compact watershed algorithm that grouped local high-similarity areas into superpixels [64]. The bounding boxes of the superpixels were used as the region proposals. The region proposal algorithm is summarized in Algorithm 1 (Supplemental material).

The region proposals generated from testing slides were converted into graphs using the procedure described above and scored by the trained GNN for predictions. The overall workflow is illustrated in Fig. 2.

3.1 Classification of TMA dataset

Baselines. The proposed GNN-based model was compared to three other baseline methods, including a simple model with multiple-layer perceptron (MLP) followed by element-wise mean-pooling, max-pooling MIL, and a recent attention-based MIL model [46]. The max-pooling and mean-pooling MIL used the standard assumption of MIL by treating the patches of a TMA core as a bag of instances. The embeddings of the patches of a core were element-wise pooled using mean-pooling (element-wise averaging) to form a bag embedding, and a linear classification (fully-connected layer) then scored the bag embedding to reach a core-level prediction. In max-pooling MIL, a linear classification head was used to score each patch embedding, and the resulting scores were averaged to produce the final score. The attention-based MIL used attention-based pooling on the embeddings of the patches to form a bag embedding by computing a weighted sum of the patch features.

Evaluation metrics. Average accuracy, average F1 score, and area under the curve (AUC) were used to compare the classification performance. The mean scores were computed from a 4-fold cross validation following the scheme described in the Materials and Methods section 2. The classification performance of the models on brightfield images of the TMA dataset is summarized in Fig. 3.

3.2 Classification of tissue section dataset

ROI classification performance. The same evaluation process was conducted on the brightfield images of the tissue section dataset with results summarized in Fig. 3. The results show that GNN-based models that consider the interactions between patches outperform patch-wise models including MLP and MIL. We give a brief explanation of the difference in model performance in the discussion section.

Region proposal and whole slide scoring. The described WSI region proposal method was evaluated on the whole tissue dataset. Region proposals were generated and scored following the scheme described in the Materials and Methods section 2. Representative results are shown in Fig. 4. The effectiveness of the region proposal method was also quantitatively evaluated. The purpose of the region proposal step is to extract candidate regions from unannotated WSIs that cover positive regions as much as possible while keeping the total number of regions needed to be scored by the downstream classification model small.

Figure 5 plots the recall against the ratio of the area in the WSI covered by region proposals. A region proposal is considered true positive if it overlaps with any of the annotations (PDAC and chronic pancreatitis). False negatives are defined as the total area of the annotated regions that are not covered by any region proposals divided by the average size of region proposals. Intuitively, decreasing the threshold that selects for high similarity areas results in more region proposals; a threshold of 0 leads to a set of region proposals that covers the whole WSI with a recall of 1. The ratio of area covered by the region proposal and the corresponding threshold returned by multi-Otsu’s method is marked with a purple line in the figure. A set of good region proposals should have a high recall of disease-positive region coverage while keeping the overall tissue coverage relatively small. Overall, the region proposals achieve an average recall of 0.909 while covering 40.1% of the tissue regions.

3.3 Effects of incorporating collagen features

Classification using only collagen features extracted from SHG. The same GNN model was used on the SHG images collected from the tissue section dataset. The node features of the graph were collagen-based features extracted from the SHG images using CurveAlign/CT-FIRE. The classification results are summarized in Fig. 6, showing that collagen-based features alone are diagnostic and the trained model is significantly better than a random classifier.

Classification using combined features from brightfield and SHG. To further demonstrate the diagnostic value of collagen-based features, brightfield-SHG graphs were constructed by concatenating the brightfield feature vector and the SHG feature vector as node features to train the GNN model. The SHG image and brightfield H&E image of an annotated region were first co-registered with pixel-wise correspondences such that the tiled patches of the two remained co-registered. The resulting concatenated feature vector contained brightfield features and SHG features of the same tissue area. Representative co-registered brightfield H&E images and SHG images of the three classes are shown in Fig. 7.

The evaluation results suggest that combining the brightfield H&E features with collagen-based features leads to improved performance for both the complex GNN-based model and the simple MLP (Fig. 6). Specifically, the GNN model using combined features achieves an accuracy of 91.3% while the accuracy using brightfield features alone is 88.9% (Fig. 6).

3.4 Feature Interpretability

ResNet extracted H&E features. The node features of H&E patches were directly computed from ImageNet pre-trained ResNet18, as suggested in [65]; this could be an efficient and practical way to extract features from WSI patches. To validate that the computed features are discriminative, the top 16 patches closest to each k-means centroids are shown in Fig. 8 in the case of k = 7. The feature vectors of patches are also embedded in 2 dimensions using t-SNE [66] and labeled with k-means clustering results in the figure. The choice of k = 7 was experimentally determined. Our pathologist validated the following: the patches cluster of class 1 corresponds to blank areas; most of the patches in class 2 have adipose tissue; patches in class 3 and class 4 contain stroma and collagen without any obviously malignant cells; patches in class 5 mostly contain chronic inflammation, sometimes in a collagenous stroma and sometimes admixed with acinar tissue; class 6 contains benign tissue, including acini, islets of Langerhans, and ducts; patches in class 7 contain cancer-appearing duct parts with disorganized collagen, and some contain stroma with acute inflammation or chronic inflammation.

Collagen-based features. In a previous study where collagen-based features were used for a simple binary linear classifier on PDAC and CP, the importance of each collagen-based feature was evaluated separately [22]. Here we show an alternative way to quantify the importance of the input feature in the setting of the GNN model. We removed one of the input features (for width, angle, straightness, length, both the mean and variance were removed) and trained the same GNN to perform classification, and observed the performance percentage drop, as shown in Fig. 9. A more significant performance drop after the removal suggests that the removed features are more important.

Importance of graph models in pancreas histomorphology. The results show that modeling WSI regions as graphs is beneficial for building a computational classifier for the differential diagnosis of PDAC and CP. A potential reason for the accuracy improvement is that the diagnostic features of PDAC and CP span multiple patches and are difficult to capture using patch-wise representations. In the histomorphology of pancreatic cancer, tumor regions are defined as malignant epithelium arranged in atypical ducts, clusters of cells or individual cells infiltrating into the surrounding tissue. The diagnostic criteria range from nuclear to architectural, including large variations in epithelial nuclear size, incomplete epithelium around ducts, angular and branching ducts, disorganized duct distribution, and ducts appearing in architecturally aberrant locations. Many of these features overlap with those of CP, including the increasing deposition of fibrillar collagen and desmoplastic stroma in the tissue microenvironment. These issues make obtaining expressive features from the image data in both patch-wise methods and MIL-based methods a challenge. Even though many MIL-based methods use bag-of-patch representations and learn a joint representation of all image patches, they could still perform incompetently as the spatial relations of the patches are not modeled and the features to aggregate in MIL operators are restricted to the field-of-view of a single patch [46, 47].

By contrast, GNNs model the patch-to-patch interactions by deploying message passing between adjacent patches and updating the node features to encode information from the neighboring patches. The updated node features can eventually incorporate features up to k hops away, determined by the number of graph convolution layers, forming representations that correspond to a larger field of view. For example, the GNN could potentially learn the features of duct organization or an incomplete duct by propagating the information of a tile with disrupted epithelium to the adjacent tiles that form a duct object 10. Therefore, compared to MIL-based models, a GNN-based model could capture the spatial arrangement and relationship of malignant cells to each other, to the stroma, and to other tissue elements. It is interesting that this kind of carry-over or comparison from field to field begins to approximate the manner in which pathologists evaluate or “read” slides in the clinical setting.

Construction of graphs from histopathological images. One important step in using graph models in WSI analysis is the construction of graphs from images. The proposed method follows a nearest neighbor scheme where each patch is connected to its non-empty adjacent patches. Though easy to implement, this method does not exploit the nested structure of tissue components and redundant texture-like information that is often spatially continuous. More complex approaches such as multiscale quadtree approximation [67] and constructing refined graphs based on the results of patch-wise spatial agglomeration [68, 69] could potentially remove redundancies and result in graphs that encode the tissue structure and relations of components in more efficient ways.

Model the interaction between histopathology modalities. Examination of H&E stained slides is the standard clinical histopathologic approach to assess surgically excised pancreatic tissue. Special modalities such as histochemistry or immunohistochemistry are seldom used to differentiate between benign (CP) and malignant tissue. However, the research applications of multimodal histopathology have been rapidly advancing with promising results that could benefit clinical histopathology.

Given our prior knowledge of the importance of fibrillar collagen in the tissue microenvironment of PDAC and CP [22], we incorporated collagen-based features extracted from SHG images, superimposed on the original H&E images, into the analysis. Even with the collagen fiber features alone, the GNN model achieves classification accuracy significantly better than a random classifier, showing that the features are discriminative for differentiating PDAC from CP. In the multimodal analysis pipeline, the feature vector pairs of registered image patches from both imaging modalities are simply concatenated to form the input to the model. After combining the collagen-based features from SHG images with the features from the corresponding H&E images, the downstream GNN achieves better accuracy in classification, as expected. To fully leverage the strength of GNNs, one could treat the features of another modality also as graph nodes and build a graph, while also connecting the nodes of two modalities that are paired spatially. Note that the inter-modality edge weights can be learned separately from the weights that connect the nodes spatially, forming an architecture. This could grant GNNs the flexibility to let information flow and interact between modalities. Another future direction is to incorporate information from patient demographics and radiographics that could provide more systematic views of the patient’s condition.

Out-of-distribution samples. The quality of tissue sampling is critical to both visual inspection and computational analysis in histopathology. For PDAC and CP detection, low-quality tissue sampling or slide preparation could lead to inadequate tissue regions that fail to support an accurate diagnosis. For example, if the pancreas lesion is not sampled to get adequate epithelium in fine-needle aspiration, the cytopathologist would not be able to make a diagnosis of adenocarcinoma. In computational analysis, tissue sampling could affect the distribution of training samples. Thus, if the tissue components included in the training set are not exclusive, the resulted classifier can encounter out-of-distribution tissue regions and give over-confident false-positive results.

In our experimental setting, the pancreas tissue sections are annotated coarsely with PDAC, CP, or normal regions. However, when applying the classifier to all candidate regions extracted from the entire tissue on a slide, the classifier predicts some false-positive regions (PDAC or CP) with high confidence. These were areas of smooth muscle (duodenal muscularis propria), or tissue from another organ entirely (mostly from autopsy cases). These were not annotated by the pathologist because they were intuitively not “regions of interest.” For the classifier, however, these were clearly out-of-distribution samples that had not been previously labeled and were therefore erroneously predicted as malignant (Fig. 11). This could be related to a known artifact of discriminative classifiers including deep neural networks [70]. The discriminative distribution density saturates at 1 even for the case where the input value to the last Sigmoid layer is arbitrary large which could be caused by an out-of-distribution sample that is very far away from any of the correct class distributions. Therefore, to facilitate the real-world applications of deep learning models in assisting clinical histopathology, it is important to incorporate uncertainty estimation into the classifier [71], or construct generative classifiers where out-of-distribution inputs could be detected [72].

In summary, we developed a GNN-based analysis pipeline for detecting and differentiating PDAC and CP in histological images. To our knowledge, this is the first DL and GNN-based computational workflow for the differentiation of PDAC and CP in histopathology. Evaluation results indicate that modeling histological images as graphs significantly improves classification accuracy. Meanwhile, conventional analysis models such as MIL can be overly simple for cases where diagnostic features span multiple patches. Furthermore, the incorporation of collagen fiber-based features increased the accuracy of the classifier, highlighting once again the potential diagnostic value of collagen characteristics and the potential clinical benefits of utilizing advanced imaging modalities.

ACKNOWLEDGEMENTS

We acknowledge members of the Eliceiri and Loeffler labs and Dr. Yin Li for useful input and discussion.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

AUTHOR CONTRIBUTIONS

Study concept and design (B.L., M.S.N., A.G.L. and K.W.E.). Analysis pipeline development, imaging experiment, and data quantification (B.L. and M.S.N.). Sample collection and pathological review (O.S. and A.G.L.). Manuscript organizing and writing (B.L., M.S.N., A.G.L. and K.W.E.).

FUNDING

We acknowledge funding from the Morgridge Institute for Research, the Semiconductor Research Corporation, and NIH grants U54CA268069, P41GM135019 and R01CA238191.

ETHICS APPROVAL AND CONSENT TO PARTICIPATE

Tissue sections were used under MetroHealth IRB CR00000311.

DATA AVAILABILITY STATEMENT

Access to the image datasets used in this study and the code is available at https://github.com/uw-loci/gnn-pccp.

Hidalgo M. Pancreatic cancer. New England Journal of Medicine. 2010;362(17):1605–1617.
Rawla P, Sunkara T, Gaduputi V. Epidemiology of pancreatic cancer: global trends, etiology and risk factors. World journal of oncology. 2019;10(1):10.
Bellizzi AM, Frankel WL. Pancreatic pathology: a practical review. Laboratory Medicine. 2009;40(7):417–426.
Mostafa ME, Erbarut-Seven I, Pehlivanoglu B, Adsay V. Pathologic classification of” pancreatic cancers”: current concepts and challenges. Chinese clinical oncology. 2017;6(6):59–59.
Mihaljevic A, Esposito I, Friess H, Kleeff J. Molecular biology, models, and histopathology of chronic pancreatitis and pancreatic cancer. European Surgery. 2009;41(6):250–267.
Kl¨oppel G. Chronic pancreatitis, pseudotumors and other tumor-like lesions. Modern pathology. 2007;20(1):S113–S131.
Yepuri N, Naous R, Richards C, Dhir M, Jain A. Poorly differentiated signet ring cell carcinoma of pancreas masquerading as chronic pancreatitis. Journal of Surgical Case Reports. 2018;2018(8):rjy218.
Dal Molin M, Zhang M, De Wilde RF, Ottenhof NA, Rezaee N, Wolfgang CL, et al. Very long-term survival following resection for pancreatic cancer is not explained by commonly mutated genes: results of whole-exome sequencing analysis. Clinical Cancer Research. 2015;21(8):1944–1950.
Strobel O, Neoptolemos J, J¨ager D, Bu¨chler MW. Optimizing the outcomes of pancreatic cancer surgery. Nature reviews Clinical oncology. 2019;16(1):11–26.
Yamamoto T, Yagi S, Kinoshita H, Sakamoto Y, Okada K, Uryuhara K, et al. Long-term survival after resection of pancreatic cancer: a single-center retrospective analysis. World journal of gastroenterology: WJG. 2015;21(1):262.
Esposito I, Hruban RH, Verbeke C, Terris B, Zamboni G, Scarpa A, et al. Guidelines on the histopathology of chronic pancreatitis. Recommendations from the working group for the international consensus guidelines for chronic pancreatitis in collaboration with the International Association of Pancreatology, the American Pancreatic Association, the Japan Pancreas Society, and the European Pancreatic Club. Pancreatology. 2020;20(4):586–593.
Dimastromatteo J, Brentnall T, Kelly KA. Imaging in pancreatic disease. Nature reviews Gastroenterology & hepatology. 2017;14(2):97–109.
Sarantis P, Koustas E, Papadimitropoulou A, Papavassiliou AG, Karamouzis MV. Pancreatic ductal adenocarcinoma: Treatment hurdles, tumor microenvironment and immunotherapy. World journal of gastrointestinal oncology. 2020;12(2):173.
Truong LH, Pauklin S. Pancreatic Cancer Microenvironment and Cellular Composition: Current Understandings and Therapeutic Approaches. Cancers. 2021;13(19):5028.
Hosein AN, Brekken RA, Maitra A. Pancreatic cancer stroma: an update on therapeutic targeting strategies. Nature Reviews Gastroenterology & Hepatology. 2020;17(8):487–505.
Conklin MW, Eickhoff JC, Riching KM, Pehlke CA, Eliceiri KW, Provenzano PP, et al. Aligned collagen is a prognostic signature for survival in human breast carcinoma. The American journal of pathology. 2011;178(3):1221–1232.
Nadiarnykh O, LaComb RB, Brewer MA, Campagnola PJ. Alterations of the extracellular matrix in ovarian cancer studied by Second Harmonic Generation imaging microscopy. BMC cancer. 2010;10(1):1–14.
Drifka CR, Loeffler AG, Mathewson K, Keikhosravi A, Eickhoff JC, Liu Y, et al. Highly aligned stromal collagen is a negative prognostic factor following pancreatic ductal adenocarcinoma resection. Oncotarget. 2016;7(46):76197.
Hanley CJ, Noble F, Ward M, Bullock M, Drifka C, Mellone M, et al. A subset of myofibroblastic cancer-associated fibroblasts regulate collagen fiber elongation, which is prognostic in multiple cancers. Oncotarget. 2016;7(5):6159.
Best SL, Liu Y, Keikhosravi A, Drifka CR, Woo KM, Mehta GS, et al. Collagen organization of renal cell carcinoma differs between low and high grade tumors. BMC cancer. 2019;19(1):1–8.
Zunder SM, Gelderblom H, Tollenaar RA, Mesker WE. The significance of stromal collagen organization in cancer tissue: An in-depth discussion of literature. Critical Reviews in Oncology/Hematology. 2020;151:102907.
Drifka CR, Tod J, Loeffler AG, Liu Y, Thomas GJ, Eliceiri KW, et al. Periductal stromal collagen topology of pancreatic ductal adenocarcinoma differs from that of normal and chronic pancreatitis. Modern Pathology. 2015;28(11):1470–1480.
Grizzi F, Fiorino S, Qehajaj D, Fornelli A, Russo C, De Biase D, et al. Computer-aided assessment of the extra-cellular matrix during pancreatic carcinogenesis: a pilot study. Journal of translational medicine. 2019;17(1):1–9.
Han M, Giese G, Bille JF. Second harmonic generation imaging of collagen fibrils in cornea and sclera. Optics express. 2005;13(15):5791–5797.
Strupler M, Pena AM, Hernest M, Tharaux PL, Martin JL, Beaurepaire E, et al. Second harmonic imaging and scoring of collagen in fibrotic tissues. Optics express. 2007;15(7):4054–4065.
Keikhosravi A, Bredfeldt JS, Sagar AK, Eliceiri KW. Second-harmonic generation imaging of cancer. Methods in cell biology. 2014;123:531–546.
Chen X, Nadiarynkh O, Plotnikov S, Campagnola PJ. Second harmonic generation microscopy for quantitative analysis of collagen fibrillar structure. Nature protocols. 2012;7(4):654–669.
Bredfeldt JS, Liu Y, Pehlke CA, Conklin MW, Szulczewski JM, Inman DR, et al. Computational segmentation of collagen fibers from secondharmonic generation images of breast cancer. Journal of biomedical optics. 2014;19(1):016007.
Valkenburg KC, De Groot AE, Pienta KJ. Targeting the tumour stroma to improve cancer therapy. Nature reviews Clinical oncology. 2018;15(6):366– 381.
Zhang J, Liu J. Tumor stroma as targets for cancer therapy. Pharmacology & therapeutics. 2013;137(2):200–215.
Pantanowitz L, Valenstein PN, Evans AJ, Kaplan KJ, Pfeifer JD, Wilbur DC, et al. Review of the current state of whole slide imaging in pathology. Journal of pathology informatics. 2011;2.
Ghaznavi F, Evans A, Madabhushi A, Feldman M. Digital imaging in pathology: whole-slide imaging and beyond. Annual Review of Pathology: Mechanisms of Disease. 2013;8:331–359.
Farahani N, Parwani AV, Pantanowitz L, et al. Whole slide imaging in pathology: advantages, limitations, and emerging perspectives. Pathol Lab Med Int. 2015;7(23-33):4321.
Huss R, Coupland SE. Software-assisted decision support in digital histopathology. The Journal of Pathology. 2020;250(5):685–692.
Srinidhi CL, Ciga O, Martel AL. Deep neural network models for computational histopathology: A survey. Medical Image Analysis. 2021;67:101813.
Ahmedt-Aristizabal D, Armin MA, Denman S, Fookes C, Petersson L. A survey on graph-based deep learning for computational histopathology. Computerized Medical Imaging and Graphics. 2021;p. 102027.
Dimitriou N, Arandjelovi´c O, Caie PD. Deep learning for whole slide image analysis: an overview. Frontiers in medicine. 2019;p. 264.
Fu Y, Jung AW, Torne RV, Gonzalez S, V¨ohringer H, Shmatko A, et al. Pan-cancer computational histopathology reveals mutations, tumor composition and prognosis. Nature Cancer. 2020;1(8):800–810.
Van der Laak J, Litjens G, Ciompi F. Deep learning in histopathology: the path to the clinic. Nature medicine. 2021;27(5):775–784.
Scarselli F, Gori M, Tsoi AC, Hagenbuchner M, Monfardini G. The graph neural network model. IEEE transactions on neural networks. 2008;20(1):61–80.
Kipf TN, Welling M. Semi-supervised classification with graph convolutional networks. arXiv preprint arXiv:160902907. 2016;.
Wu F, Souza A, Zhang T, Fifty C, Yu T, Weinberger K. Simplifying graph convolutional networks. In: International conference on machine learning. PMLR; 2019. p. 6861–6871.
Aygu¨ne¸s B, Aksoy S, Cinbi¸s RG, K¨osemehmeto˘glu K, Onder S,¨ Uner A.¨ Graph convolutional networks for region of interest classification in breast histopathology. In: Medical Imaging 2020: Digital Pathology. vol. 11320. International Society for Optics and Photonics; 2020. p. 113200K.
Adnan M, Kalra S, Tizhoosh HR. Representation learning of histopathology images using graph neural networks. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops; 2020. p. 988–989.
Pati P, Jaume G, Foncubierta-Rodr´ıguez A, Feroce F, Anniciello AM, Scognamiglio G, et al. Hierarchical graph representations in digital pathology. Medical image analysis. 2022;75:102264.
Ilse M, Tomczak J, Welling M. Attention-based deep multiple instance learning. In: International conference on machine learning. PMLR; 2018. p. 2127–2136.
Li B, Li Y, Eliceiri KW. Dual-stream multiple instance learning network for whole slide image classification with self-supervised contrastive learning. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition; 2021. p. 14318–14328.
Keikhosravi A, Li B, Liu Y, Conklin MW, Loeffler AG, Eliceiri KW. Non-disruptive collagen characterization in clinical histopathology using cross-modality image synthesis. Communications biology. 2020;3(1):1–12.
Bankhead P, Loughrey MB, Fernandez JA, Dombrowski Y, McArt DG, Dunne PD, et al. QuPath: Open source software for digital pathology image analysis. Scientific reports. 2017;7(1):1–7.
Edelstein AD, Tsuchida MA, Amodaj N, Pinkard H, Vale RD, Stuurman N. Advanced methods of microscope control using µManager software. Journal of biological methods. 2014;1(2).
Pinkard H, Stuurman N, Ivanov IE, Anthony NM, Ouyang W, Li B, et al. Pycro-Manager: open-source software for customized and reproducible microscope control. Nature methods. 2021;18(3):226–228.
Keikhosravi A, Li B, Liu Y, Eliceiri KW. Intensity-based registration of bright-field and second-harmonic generation images of histopathology tissue sections. Biomedical Optics Express. 2020;11(1):160–173.
Pielawski N, Wetzer E, Ofverstedt J, Lu J, W¨ahlby C, Lindblad J, et al.¨ CoMIR: Contrastive multimodal image representation for registration. Advances in neural information processing systems. 2020;33:18433–18444.
Schindelin J, Arganda-Carreras I, Frise E, Kaynig V, Longair M, Pietzsch T, et al. Fiji: an open-source platform for biological-image analysis. Nature methods. 2012;9(7):676–682.
Macenko M, Niethammer M, Marron JS, Borland D, Woosley JT, Guan X, et al. A method for normalizing histology slides for quantitative analysis. In: 2009 IEEE international symposium on biomedical imaging: from nano to macro. IEEE; 2009. p. 1107–1110.
He K, Zhang X, Ren S, Sun J. Deep residual learning for image recognition. In: Proceedings of the IEEE conference on computer vision and pattern recognition; 2016. p. 770–778.
Liu Y, Keikhosravi A, Pehlke CA, Bredfeldt JS, Dutson M, Liu H, et al. Fibrillar collagen quantification with curvelet transform based computational methods. Frontiers in bioengineering and biotechnology. 2020;8:198.
He K, Zhang X, Ren S, Sun J. Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. In: Proceedings of the IEEE international conference on computer vision; 2015. p. 1026–1034.
Kingma DP, Ba J. Adam: A method for stochastic optimization. arXiv preprint arXiv:14126980. 2014;.
Loshchilov I, Hutter F. Sgdr: Stochastic gradient descent with warm restarts. arXiv preprint arXiv:160803983. 2016;.
Sculley D. Web-scale k-means clustering. In: Proceedings of the 19th international conference on World wide web; 2010. p. 1177–1178.
Huang Z, Li Y. Interpretable and accurate fine-grained recognition via region grouping. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition; 2020. p. 8662–8672.
Liao PS, Chen TS, Chung PC, et al. A fast algorithm for multilevel thresholding. J Inf Sci Eng. 2001;17(5):713–727.
Neubert P, Protzel P. Compact watershed and preemptive slic: On improving trade-offs of superpixel segmentation algorithms. In: 2014 22nd international conference on pattern recognition. IEEE; 2014. p. 996–1001.
Lu MY, Williamson DF, Chen TY, Chen RJ, Barbieri M, Mahmood F. Data-efficient and weakly supervised computational pathology on whole-slide images. Nature Biomedical Engineering. 2021;5(6):555–570.
Van der Maaten L, Hinton G. Visualizing data using t-SNE. Journal of machine learning research. 2008;9(11).
Lafert´e JM, P´erez P, Heitz F. Discrete Markov image modeling and inference on the quadtree. IEEE Transactions on image processing. 2000;9(3):390–404.
Sharma H, Alekseychuk A, Leskovsky P, Hellwich O, Anand RS, Zerbe N, et al. Determining similarity in histological images using graph-theoretic description and matching methods for content-based image retrieval in medical diagnostics. Diagnostic pathology. 2012;7(1):1–20.
Zheng Y, Jiang Z, Shi J, Xie F, Zhang H, Luo W, et al. Encoding histopathology whole slide images with location-aware graphs for diagnostically relevant regions retrieval. Medical image analysis. 2022;76:102308.
Hein M, Andriushchenko M, Bitterwolf J. Why relu networks yield high-confidence predictions far away from the training data and how to mitigate the problem. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition; 2019. p. 41–50.
Leibig C, Allken V, Ayhan MS, Berens P, Wahl S. Leveraging uncertainty information from deep neural networks for disease detection. Scientific reports. 2017;7(1):1–14.
Ren J, Liu PJ, Fertig E, Snoek J, Poplin R, Depristo M, et al. Likelihood ratios for out-of-distribution detection. Advances in Neural Information Processing Systems. 2019;32.

There is NO conflict of interest to disclose.

Algorithm1.pdf
Algorithm 1

Download PDF

Version 1

posted

You are reading this latest preprint version

Differentiation of Pancreatic Ductal Adenocarcinoma and Chronic Pancreatitis using Graph Neural Networks on Histopathology and Collagen Fiber Features

Status:

Version 1

Abstract

Figures

1 Introduction

2 Materials And Methods

2.1 Sample sets and pathology annotations

2.2 Second-harmonic generation microscopy imaging of samples

2.3 Graph neural networks

2.3.1 Graph construction

2.3.2 Graph convolutional network architecture

2.3.3 Network training specifics

2.4 Region proposal in whole slide images

3 Results

3.1 Classification of TMA dataset

3.2 Classification of tissue section dataset

3.3 Effects of incorporating collagen features

3.4 Feature Interpretability

4 Discussion And Future Work

Declarations

ACKNOWLEDGEMENTS

CONFLICT OF INTEREST

AUTHOR CONTRIBUTIONS

FUNDING

DATA AVAILABILITY STATEMENT

References

Additional Declarations

Supplementary Files

Status:

Version 1