A Novel Framework for Quantification of Immune Cell Density and Characterization of Tumor-Immune Spatial Relationships in Tumor Microenvironment

Understanding the impact of tumor biology on the composition of nearby cells often requires characterizing the impact of biologically distinct tumor regions. Biomarkers have been developed to label biologically distinct tumor regions, but challenges arise because of differences in the spatial extent and distribution of differentially labeled regions. In this work, we present a framework for systematically investigating the impact of distinct tumor regions on cells near the tumor borders, accounting their cross spatial distributions. We apply the framework to multiplex immunohistochemistry (mIHC) studies of pancreatic cancer and show its efficacy in demonstrating how biologically different tumor regions impact the immune response in the tumor microenvironment. Furthermore, we show that the proposed framework can be extended to largescale whole slide image analysis.


Introduction
The analysis of spatial relationships among tumor infiltrating immune cell types relative to cancer biomarker expression can be used to develop more effective approaches for immunotherapeutic intervention in a wide range of solid tumors, including pancreatic ductal adenocarcinoma (PDAC) [1,2]. The use of high-throughput slide scanning and innovative machine learning-based algorithms to relate cancer biomarker expression to the inflammatory tumor microenvironment (TME), using either state-ofthe-art brightfield, multiplexed immunohistochemical (mIHC) or multiplexed immunofluorescent (mIF) protocols can provide automated, unbiased, and reproducible analyses across hundreds of specimens. In our work, we seek to develop methods that characterize the effects of differential tumor regions on immune cells. We focus on the analysis of biologically relevant subclasses of T cells and macrophages, using mIHC, to generate data on immune cell quantity and spatial relationships to tumor cell, relative to pan-cytokeratin, a marker of tumor epithelial cells, and to keratin 17 (K17), a biomarker of the most biologically aggressive form of PDAC [3,4,5]. The mIHC-based analysis has provided unique insight into the spatial relationships among cells within the complex TME, including infiltrating immune cells, cancer cells, and stromal cells [6]. An earlier work has reported a deep learning approach and a multiparameter quantification strategy to enable comprehensive phenotyping of immune complexity [7]. Building upon the deep learning approach, we now report a modified method to carry out a comparative analysis between immune cells and biomarkerpositive tumor cells in order to further develop an understanding of the TME. The primary innovation of our work is the introduction of a novel, image-based scoring system of cancer biomarkers. The scoring system computes a metric, called "influence score", by a spatial analysis of deep learning-based segmentation and classification results, to measure and compare the influence of tumor regions on different types to immune cells. The proposed metric is carefully designed with biological rationale. We demonstrate how the influence score can be applied to multiplex IHC PDAC data to analyze the effect of K17 on different types of immune cells.

Methodology
Brightfield mIHC of tissue sections was performed using an Ultra Autostainer (Roche/Ventana, Oro Valley, AZ). Antibodies for CD4, CD8, CD16, CD163, pancytokeratin, and K17 were provided by Roche/Ventana. mIHC staining was performed using horseradish peroxidase (HRP) and alkaline phosphatase (AP)-based protocols with different colored chromogens (Red: CD4, Purple: CD8, Yellow: CD16, Green: CD163, Teal: pancytokeratin, and Brown: K17) to enable multispectral imaging of diverse immune cell populations within the cancer microenvironment. A deep learning analysis workflow was used to detect and classify cells in whole slide images (WSIs) of mIHC stained tissue specimens. The WSIs were generated at a resolution of 0.346 μm per pixel using an Olympus VS120 digital microscope (Olympus, Tokyo, Japan). The analysis workflow is illustrated in Fig. 1. The workflow consists of two deep learning models combined in an ensemble configuration [7]. Briefly, pixel-wise predictions from a color auto-encoder (ColorAE) model [8] and a U-Net model [9] are combined to create multi-class masks. A labeled pixel in a mask indicates to which type of cell (e.g., CD4, CD8 lymphocyte) the pixel belongs. The masks are further analyzed as described in the rest of this section. The notations used in the following sections are described in Table S1 (Supplementary Material).

Cell Detection and Classification
The ensemble of ColorAE and U-Net was developed in [7] for detection and classification of cells in mIHC images. Each model is trained separately, and predictions from each model are combined in the inference phase to create multi-class masks. The predictions can be combined in four ways (union, intersection, Union anchor AE, Union anchor U-Net). We used the Union anchor U-Net configuration, which we experimentally validated to be the best performing configuration for our dataset.

Expressing WSIs in terms of K17
We partition each mask into a 20 x 20 grid. We divide the mask into patches of 29 x 29 pixels (~ 10 x 10 microns). We filter the mask patches keeping only those in the regions of interest, and furthermore, excluding those intersecting with the marked non-tumor regions. Using the set of mask patches, we compute a K17 score, which quantitatively characterizes a WSI in terms of tumor percentage, as follows: In equation (1)

Quantification of Tumor-Immune Cell Spatial Relationship
We propose a scoring system (illustrated in Fig. 3.) to quantify the influence of tumor microenvironment on the distribution of different categories of immune cells. For this purpose, we introduce a novel metric called tumor/stroma influence score denoted by by calculated by following equation: The equation represents the approximate count of each immune cell (numerator) normalized by total influence area (denominator). In the equation, represents immune cells (CD4, CD8, CD16 and CD163) and represents the tumor boundary explained later in this section. To calculate the influence area, we define tumor nest area using 29x29 pixels mask containing only tumor cell denoted as . We call a tumor nest mask ( ), if the percentage of tumor cell area within that mask area is more than 35% (decided by expert pathologists). We further subclassify into = 17+ mask if more than 40% area of the tumor cells contain K17+ (decided by expert pathologists) else it is considered to be = 17− . The = 17+ and = 17− are visually illustrated in Fig. 3(B). We accumulate the tumor nest masks, define tumor nests with K17+ and K17-regions and get a rectangular boundary containing two distinct parts denoted by marker + and marker − . The tumor nest boundaries with + (brown colored lines) and − (teal-colored lines) are depicted in Fig. 3(C). We are representing a tumor nest having two distinct regions with boundaries having two distinct parts so that we can compute the approximate relationship between tumor and immune cells in an experimentally achievable time. Next, we define a partitioned influence area having K17+ and K17-influenced region. To calculate the influence area, we consider the region that falls outside the tumor nest and within the vicinity of 100μm tumor microenvironment from a tumor nest boundary.

Fig. 3. Mechanism of Distance Based Tumor-Immune Spatial Relationships Calculation
A conceptual diagram is shown in Fig. 2(B). For this purpose, we make use of the 29x29 pixels mask that is not classified as tumor nest mask denoted by, = − . This represents the discretized non tumor nest space within the tumor region. Considering center point of a single mask ( ) with index ( ) as ( ) (where, ∶ [1, | |] ∈ ) and boundary , we classify each as + or − influenced mask. We denote a single mask ( ) as K17+ influenced mask denoted by + ( ) if the euclidean distance between center point of that mask and the nearest + boundary is less than the nearest − boundary and it in 100μm distance from the + . Otherwise, if the distance between center point and the nearest − boundary is less than the nearest + boundary and it in 100μm distance from t e − then it is considered K17-influenced mask denoted by This method partitions the influence space into two regions: + influenced region (brown colored squares) and − influenced region (teal colored squares) depicted in Fig. 3(D). In a similar manner, considering ( , ) (where, ∶ [1,29 × 29] ∈ ) as a single pixel with index ( ) of immune cells in a single mask, ( ) with index ( ), we can get K17+ and K17-influenced pixels which is expressed by the following equation: Based on equation ( (5) With same technique we can calculate the total influenced immune cell area expressed by following equation: We normalize the influenced cell area by average immune cell area to get an approximate influenced cell count . We also use this influenced cell count for evaluation purpose using random cell distribution (see Section 5). Our expert pathologists approximated the average lymphocyte (CD4, CD8) to be a circle with an approximate radius of 4 microns and macrophage (CD16, CD163) with a radius of 8 microns. Therefore, the area of the average lymphocyte is 2 = 3.1416 × 4 2 = 50 2 or 144 2 (approx.), and the area of the average macrophage is 201 2 or 580 2 (approx.) for our calculation. Finally, we plug in equation (5)  and and thus get our influence score for each immune cell.

Dataset and Model Training
The training and validation datasets for the deep learning models were generated from 23 WSIs. Expert pathologists manually put labeled dots into the vicinity of each cell over 1000x1000 pixel tiles extracted from these 23 WSIs. and assigning a dot's label to all t e pixels in t e overlapping superpixel. e subdivided the tiles into 16 400400 patches with an overlap of 200 pixels in each dimension to generate ground truth. Figure 4 shows the annotation workflow. The training and validation datasets generated through this process are described in Table 1.

Model Validation and Experimental Setup
For quantitative evaluation of our detection and classification model we used F1 score illustrated in Table 2. We also performed qualitative analysis (See Fig. S1. in the supplementary material). Based on both quantitative and qualitative results, we found Union Anchor U-Net is the best model to perform detection and classification.

Results and Conclusions
We have implemented our framework on 40 WSIs to find a pattern between immune cells and tumor biomarkers. Fig. 6 sho s our distance analysis result using 25μm instead of 100μm for 40 S s. or evaluation e ave used 25 omogeneous poisson point process [12,13,14] simulations to scatter each type of cell on the environment space randomly. This simulation addresses our null hypothesis that there is a relationship between the tumor boundary and cell distribution. In the figure, the x-axis represents K17 score, and y-axis represents Tumor/Stroma Influence score. Each group of four bars in the graph corresponds with a single WSI. In the figure, we have only shown results for CD163, where we can observe that the real Tumor/Stroma Influence score largely deviates from the simulation confidence band and there is a noisy linear trend with respect to Influence Score. For an ideal case (i.e., if the cells are randomly distributed) we can observe that K17+ and K17-simulation confidence are approximately equivalent to each other. We have calculated the mean absolute deviation (MAD) between 25 simulations mean and real scenarios for both K17+ and K17-, which is tabulated in Table 3. It is clarified from the table that CD8 and CD163 have higher deviation from the randomness compared to CD4 and CD16, implicating CD8 and CD163 have spatial patterns.  Though this initial noisy pattern and correlation does not conclude the full story of tumor-immune relationship, this preliminary result is inspiring for us to conduct more investigative research on a larger set of data and subdivide into groups to fully understand the internal mechanism of Tumor-Immune relationship. From a technical perspective, we also intend to use the spatial results as a prior for developing a more advanced deep learning-based framework.