The distribution of tumor-infiltrating lymphocytes (TILs) within the tumor microenvironment provides strong prognostic value, which is increasingly important with the arrival of new immunotherapy modalities. Both visual and image analysis–based assays are developed to assess the immune contexture of the tumors. We propose an automated method based on grid subsampling of microscopy image analysis data to extract the tumor-stroma interface zone (IZ) of controlled width. The IZ is a ranking of tissue areas by their distance to the tumor edge, which is determined by a set of explicit rules. TIL density profiles across the IZ are used to compute a set of novel immunogradient indicators that reflect TIL gradient towards the tumor. We applied this method on CD8 immunohistochemistry images of surgically excised hormone receptor–positive breast and colorectal cancers to predict overall patient survival. In both cohorts, the immunogradient indicators enabled strong and independent prognostic stratification, outperforming clinical and pathologic variables. Patients with breast cancer with low immunogradient levels had a prominent decrease in survival probability 5 years after surgery. Our study provides proof of concept that data-driven, automated, operator-independent IZ sampling enables spatial immune response measurement in the tumor-host interaction frontline for prediction of disease outcomes.

One of the most prominent discoveries in modern genetics—analysis of the whole human genome sequence—has made it possible to determine specific genetic mutations associated with cancer. This determination has led to the definition of cancer as a cell disease caused by genetic mutations. Although genetic mechanisms explain many aspects of tumor progression, many hallmarks of cancer, including host immune and inflammatory response, angiogenesis, and metabolic disarrangements, evolve in the context of the tumor microenvironment (TME). In particular, a major driving force of local tumor–host tissue interactions are inflammatory and immune cell infiltrates. Solid tumors are infiltrated by both innate immunity cells (natural killer cells, macrophages, neutrophils, and phagocytes) and adaptive immunity cells (T lymphocytes, B lymphocytes, and dendritic cells). Recently, discoveries of the mechanisms by which cancer cells inhibit host antitumor immune response and new immunomodulating therapies have shifted focus toward a search for antitumor immune response components and biomarkers.

^{,}^{2}

^{4}

^{5}

^{,}Tumor-infiltrating lymphocytes (TILs) and their distributions within TME compartments have been reported as potential prognostic and predictive biomarkers in various types of cancer. Studies in clinical and experimental settings have revealed the prognostic role of CD3 Because TILs are represented by several subsets of T and B cells with complex interactions and roles, their assessment in the TME requires taking both functional and spatial aspects into account to understand their roles as major components of antitumor response. A comprehensive study of colorectal cancer (CRC) immunome by Galon et al, based on digital image analysis (DIA) of immunohistochemistry (IHC) slides, revealed that the densities of CD3 Recent studies have found prognostic and predictive value of high TIL infiltration in triple-negative and human epidermal growth factor receptor 2–positive breast cancer (BC). Naturally, intratumoral T cells are considered an important cornerstone in the emerging concept of Immunogram—a comprehensive assessment of a patient's antitumor response to guide immunotherapies in various cancers. The success of these studies has spawned into a search for novel image analytics methods, also benefiting from multiple immunofluorescence techniques, to develop combinatorial image biomarkers for colorectal, prostate, and renal cancer.

^{+}, CD4^{+}, CD8^{+}, and FOXP3^{+}TILs in many types of solid human tumors, such as melanoma, colorectal, breast, lung, bladder, prostate, renal, and hepatocellular carcinomas.^{5}

^{,}^{, }^{8}

^{, }^{9}

^{, }^{10}

^{11}

^{12}

^{+}and CD8^{+}TILs in the core tumor and the invasive margin (IM) correlate with the outcome of the disease. This discovery led to a clinically validated Immunoscore indicator, which was superior to the conventional TNM staging system.^{13}

^{14}

^{,}^{15}

- Denkert C.
- von Minckwitz G.
- Darb-Esfahani S.
- Lederer B.
- Heppner B.I.
- Weber K.E.
- Budczies J.
- Huober J.
- Klauschen F.
- Furlanetto J.
- Schmitt W.D.
- Blohmer J.U.
- Karn T.
- Pfitzner B.M.
- Kummel S.
- Engels K.
- Schneeweiss A.
- Hartmann A.
- Noske A.
- Fasching P.A.
- Jackisch C.
- van Mackelenbergh M.
- Sinn P.
- Schem C.
- Hanusch C.
- Untch M.
- Loibl S.

Tumour-infiltrating lymphocytes and prognosis in different subtypes of breast cancer: a pooled analysis of 3771 patients treated with neoadjuvant therapy.

*Lancet Oncol.*2018; 19: 40-50

^{16}

^{, }^{17}

^{, }^{18}

^{19}

^{, }^{20}

^{, }At the core of any method that aims to obtain a better understanding of tumor immune contexture is the task of quantifying the individual immune cell subtypes and their location relative to the tumor cells. The enumeration of immune cells in the TME implies not only accurate detection of the tumor and stroma regions but also a need for a clear and reproducible delineation of the IM. Regardless of whether visual or DIA methods for outlining the IM are applied, the definitions of the IM remain rather ambiguous. An early description of the IM configuration was proposed by Jass et al in 1986, who studied histomorphologic prognostic indicators in rectal carcinoma and defined two different configurations of the IM: expansive (or pushing) and infiltrative. A pushing IM is identified visually when a clear delineation of the tumor and host tissue is possible during examination of the histologic slide. Tumors with an infiltrative IM configuration have a relatively irregular growth pattern in which it is difficult to delineate host tissue from tumor cell aggregates. Many other studies used the IM definition by Mlecnik et al: a 1-mm-wide area around the border separating the host tissue from malignant glands. However, this definition does not provide an explicit explanation of the border; it actually requires an expert's judgment to manually draw it. This remains a source of bias because it leads to interobserver and intraobserver variance in tumors with irregular and highly infiltrative growth patterns, and it surely decreases the capacity of analysis, even if other analysis steps are automated. The only way to achieve a faster manual IM annotation is to simplify the IM shapes, which means leaving out the finer structural details. For instance, in the Immunoscore border, the stromal pathways within the core tumor may not be included as part of the TME. Consequently, the informative power and clinical utility of TILs and other TME-context assays may be underachieved. Recently, Harder et al applied a

^{8}

^{,}^{22}

- Hendry S.
- Salgado R.
- Gevaert T.
- Russell P.A.
- John T.
- Thapa B.
- et al.

Assessing tumor-infiltrating lymphocytes in solid tumors: a practical review for pathologists and proposal for a standardized method from the international immunooncology biomarkers working group: part 1: assessing the host immune response, TILs in invasive breast carcinoma and ductal carcinoma in situ, metastatic tumor deposits and areas for further research.

*Adv Anat Pathol.*2017; 24: 235-251

^{,}^{23}

- Hendry S.
- Salgado R.
- Gevaert T.
- Russell P.A.
- John T.
- Thapa B.
- et al.

Assessing tumor-infiltrating lymphocytes in solid tumors: a practical review for pathologists and proposal for a standardized method from the international immuno-oncology biomarkers working group: part 2: TILs in melanoma, gastrointestinal tract carcinomas, non-small cell lung carcinoma and mesothelioma, endometrial and ovarian carcinomas, squamous cell carcinoma of the head and neck, genitourinary carcinomas, and primary brain tumors.

*Adv Anat Pathol.*2017; 24: 311-335

^{8}

^{24}

^{24}

^{,}^{25}

^{26}

- Mlecnik B.
- Bindea G.
- Kirilovsky A.
- Angell H.K.
- Obenauf A.C.
- Tosolini M.
- Church S.E.
- Maby P.
- Vasaturo A.
- Angelova M.
- Fredriksen T.
- Mauger S.
- Waldner M.
- Berger A.
- Speicher M.R.
- Pages F.
- Valge-Archer V.
- Galon J.

The tumor microenvironment and Immunoscore are critical determinants of dissemination to distant metastasis.

*Sci Transl Med.*2016; 8: 327ra26

^{20}

*tissue phenomics*approach to search for image-based biomarkers in their study of prostate cancer recurrence prediction. In particular, their DIA step used morphologic operations to automatically delineate tumor gland and stroma areas and subsequently sample the tumor border as a region reaching equally far to both tumor and stroma regions. Another recent study proposed a data-driven method to discover immune contexture biomarkers; however, it included a manual step of tumor area delineation.In this study, we present a novel set of immunogradient indicators based on a new method of automated grid-based extraction of the tumor-stroma interface zone (IZ). The method first identifies the tumor edge (TE) using a set of explicit rules based on IHC DIA data. Subsequently, the IZ is extracted and ranked by the distance around the TE to allow computing TIL density profiles across the IZ. The indicators, reflecting TIL density gravitation toward the tumor, were independent prognostic factors of better overall survival (OS) in patients with hormone receptor–positive BC and CRC

## Materials and Methods

### Study Population and Tumor Characteristics

This study was performed in BC and CRC cohorts. The BC cohort consisted of 102 patients diagnosed with early hormone receptor–positive invasive ductal breast carcinoma as described previously. Briefly, the patients were women aged 27 to 87 years (median age, 59 years) who underwent surgery in 2007 to 2009 at the National Cancer Institute of Lithuania and were investigated at the Lithuanian National Center of Pathology. The CRC specimens were obtained from 101 patients diagnosed with stage I to III of primary invasive adenocarcinoma and treated by surgical excision without preoperative therapy in 2010 at Vilnius University Hospital Santaros Clinics and investigated at the National Center of Pathology. The CRC cohort included 60 women and 41 men 45 to 89 years old (median age, 70 years). The clinicopathologic and follow-up characteristics of the BC and CRC patient cohorts are summarized in Table 1.

^{28}

^{,}^{29}

Table 1Patient and Tumor Clinicopathologic Characteristics

Characteristic | Breast cancer cohort | Colorectal cancer cohort |
---|---|---|

Patients, n | 102 | 101 |

Age, years | ||

Median | 59 | 70 |

Range | 27–87 | 45–89 |

Sex, n (%) | ||

Female | 102 (100) | 60 (59.0) |

Male | 0 | 41 (41.0) |

Follow-up, months | ||

Median | 112 | 66 |

Range | 17–121 | 2–75 |

Died, n (%) | ||

5-year follow-up | 8 (7.8) | 29 (28.7) |

10-year follow-up | 22 (21.6) | Not available |

Histological grade, n (%) | ||

G1 | 23 (22.5) | 5 (4.9) |

G2 | 48 (47.1) | 85 (84.2) |

G3 | 31 (30.4) | 11 (10.9) |

Tumor invasion stage (pT), n (%) | ||

T1 | 56 (54.9) | 5 (4.9) |

T2 | 46 (45.1) | 19 (18.8) |

T3 | 0 | 62 (61.4) |

T4 | 0 | 15 (14.9) |

Lymph node metastasis status (pN), n (%) | ||

N0 | 55 (53.9) | 57 (56.4) |

N1 | 35 (34.3) | 24 (23.8) |

N2 | 9 (8.8) | 19 (18.8) |

N3 | 3 (2.9) | 1 (1.0) |

The patient studies were approved by and performed in accordance with the guidelines stated by the Lithuanian Bioethics Committee (protocol number 40, April 26, 2007, update September 12, 2017, for the BC patient cohort; protocol numbers L-13-03/1 and L-13-03/2 for CRC patient cohort). Informed written consent was collected from all patients prior to inclusion in the study.

### IHC

Formalin-fixed, paraffin-embedded tissue was used for the study; tumor tissue samples that contained the maximum content of invasive tumor tissue (one formalin-fixed, paraffin-embedded block per patient) were selected. The FFPE tissue sections were cut at 3-μm thickness and mounted on positively charged slides. IHC was performed by Roche Ventana BenchMark ULTRA automated slide stainer (Ventana Medical Systems, Tucson, AZ). Antibodies against cytotoxic T-cell marker CD8 (clone C8/144 B, Dako, Glostrup, Denmark; antibody dilution 1:400) was used followed by use of an ultraView Universal DAB Detection kit (Ventana Medical Systems). The sections were counterstained with Mayer hematoxylin. Positive staining controls were performed in parallel using paraffin-embedded human tonsil tissue.

### Image Acquisition and Analysis

The IHC slides were digitized at ×20 magnification (0.5-μm resolution) using a ScanScope XT Slide Scanner (Leica Aperio Technologies, Vista, CA). DIA of the whole slide images was performed using HALO version 2.2.1870 (Indica Labs, Corrales, NM). Cancer-specific artificial intelligence tissue classifiers were trained using HALO AI module to segment tumor, stroma, and background (consisting of necrosis, artifacts, and glass) and tumor, stroma, lymphoid follicles, necrosis, and background (consisting of glass and artifacts) in BC and CRC, respectively. Subsequently, the Multiplex IHC algorithm version 1.2 in HALO was used to detect and extract coordinates of cytoplasmic CD8

^{+}cells. For quality assurance, all image analysis results were reviewed by a pathologist (A.Laurinavicius): three whole slide images were excluded because of the need to train a specific tissue classifier for the mucinous type of tumor histology. Examples of IHC and DIA analysis output images are presented in Figure 1, A and B .### TE Detection

The automated edge extraction uses the pixel-wise classification of the tissue by DIA (Figure 1, A and B). A hexagonal grid (hexagon side, 65 μm) is overlaid in a random position to subsample the DIA results (Figure 1C), as described in previous studies. Inside each hexagon, a set

^{28}

^{,}^{29}

*V*of data variables (IHC-positive and IHC-negative cells counts and tissue class areas) is collected. For robustness to varying hexagon size, the method relies on area fractions inside the grid elements instead of absolute areas. Regardless of how many classes were identified by the classifier, only three area fractions are of interest here: tumor area fraction*t*, stroma area fraction*s*, and background area fraction*b*combined from the remaining classes extracted by DIA. Figure 1D illustrates area fractions (*t*,*s*, and*b*) as red, green, and blue, respectively.For any grid hexagon, $h$, denote the lookup of data variable $v\in V\phantom{\rule{0.25em}{0ex}}$as $v\left(h\right)$, for example, tumor area fraction as $t\left(h\right)$ or simply $t$ when $h$ is unambiguous. To identify hexagons with abrupt changes, the change in any $v\in V\phantom{\rule{0.25em}{0ex}}$across the hexagonal grid is calculated using derivatives along the three hexagonal axes; ${d}_{x}^{v}\phantom{\rule{0.25em}{0ex}}$denotes derivative of variable

*v*along axis*x*(at implicit hexagon*h*); for example, ${d}_{x}^{t}$ is the change of tumor area fraction along*x*. Total change is given by the norm of the gradient $\left|{\nabla}_{v}\left(h\right)\right|$. Thus, the total change in tumor area fraction used for interface extraction is as follows (again omitting*h*):$\left|{\nabla}_{t}\right|=\sqrt{{\left({d}_{x}^{t}\right)}^{2}+{\left({d}_{y}^{t}\right)}^{2}+{\left({d}_{z}^{t}\right)}^{2}}.$

(1)

The numerical implementation of the hexagonal gradient is similar to the hexagonal gradients as calculated according to previous studies. However, whereas those methods used a linear combination among derivatives to optimize for computation speed, here all derivatives contribute to the gradient to maximize the level of detail extracted. In Figure 1E, $\left|{\nabla}_{t}\right|$ is overlaid in the red color channel and illustrates how it captures

^{30}

^{,}^{31}

*all*changes of tumor.For the extraction of the TE, only hexagons where the tumor area changes because of neighboring stroma regions are of interest. Thus, $\left|{\nabla}_{t}\right|\phantom{\rule{0.25em}{0ex}}$is separated into a tumor-stroma part and a tumor-background part by weighting it with the relative changes in stroma and background area fractions. For ${d}_{x}^{t}$ the relative changes in

*s*and*b*are as follows:${s}^{{d}_{x}^{t}}={d}_{x}^{t}\phantom{\rule{0.25em}{0ex}}\frac{\left|{d}_{x}^{s}\right|}{{\sum}_{i\ne t}\left|{d}_{x}^{i}\right|}={d}_{x}^{t}\phantom{\rule{0.25em}{0ex}}\frac{\left|{d}_{x}^{s}\right|}{\left|{d}_{x}^{s}\right|+\left|{d}_{x}^{b}\right|}\phantom{\rule{0.25em}{0ex}}\text{and}\phantom{\rule{0.25em}{0ex}}{b}^{{d}_{x}^{t}}={d}_{x}^{t}\phantom{\rule{0.25em}{0ex}}\frac{\left|{d}_{x}^{b}\right|}{{\sum}_{i\ne t}\left|{d}_{x}^{i}\right|}={d}_{x}^{t}\phantom{\rule{0.25em}{0ex}}\frac{\left|{d}_{x}^{b}\right|}{\left|{d}_{x}^{s}\right|+\left|{d}_{x}^{b}\right|}$

(2)

The rationale is that if the amount of background area changes very little across some hexagons, any change in tumor area can be interpreted as being

*caused*by change in stroma area and vice versa. Note that ${}_{s}{}^{\phantom{\rule{0.25em}{0ex}}}{d}_{x}^{t}+{}_{b}{}^{\phantom{\rule{0.25em}{0ex}}}{d}_{x}^{t}={d}_{x}^{t}$ ensures that no information is lost or added but merely separated. The separation weights are similar along*y*and*z*, and thus total change $\left|{\nabla}_{t}\right|$ can be separated by*s*and*b*:$\left|{\nabla}_{t}\right|=\phantom{\rule{0.25em}{0ex}}\left|{\nabla}_{t}^{s}+{\nabla}_{t}^{b}\right|$

(3)

Figure 1F illustrates the separation of tumor changes into tumor-stroma and tumor-background changes on the hexagonal grid using

as green and

as blue.

$\left|{\nabla}_{t}^{s}\right|=\phantom{\rule{0.25em}{0ex}}\sqrt{{\left({}_{s}{}^{\phantom{\rule{0.25em}{0ex}}}{d}_{x}^{t}\right)}^{2}+{\left({}_{s}{}^{\phantom{\rule{0.25em}{0ex}}}{d}_{y}^{t}\right)}^{2}+{\left({}_{s}{}^{\phantom{\rule{0.25em}{0ex}}}{d}_{z}^{t}\right)}^{2}}$

(4)

as green and

$\left|{\nabla}_{t}^{b}\right|=\phantom{\rule{0.25em}{0ex}}\sqrt{{\left({}_{b}{}^{\phantom{\rule{0.25em}{0ex}}}{d}_{x}^{t}\right)}^{2}+{\left({}_{b}{}^{\phantom{\rule{0.25em}{0ex}}}{d}_{y}^{t}\right)}^{2}+{\left({}_{b}{}^{\phantom{\rule{0.25em}{0ex}}}{d}_{z}^{t}\right)}^{2}}$

(5)

as blue.

To consistently classify the hexagons into tumor, stroma, background, TE, and tumor-background interface, the normalized tumor gradient, denoted ${\nabla}_{t}\in \left[0;1\right]$

_{,}is used as a*probability*for whether a hexagon should be part of the TE. The complete set of steps for hexagon classification is as follows:i) Hexagons with abrupt changes in tumor area fraction are identified by testing if ${\nabla}_{t}>0.5$. For all hexagons for which this is the case, the maximum of ${\nabla}_{t}^{s}$ and ${\nabla}_{t}^{b}$ will subsequently determine whether the hexagon is part of the TE or the tumor-background interface, respectively.

ii) Hexagons not deemed part of the TE by step 1 are grouped into tumor, stroma, or background by the maximum of area fractions of

*t*,*s*, and*b*because t+ s+b=1, where t,s, and b∈[0;1].iii) TE invasive regions are identified as hexagons classified as tumor or stroma in step 2 and where the tumor and stroma areas are of similar amount. Specifically, if $\left|t-s\right|<0.25$ the hexagon is deemed part of the invasive areas of TE.

Figure 1G shows the result of hexagon classification.

### Tumor-Stroma IZ Extraction around the TE

Using a simple hexagonal distance transform, each hexagon's shortest distance to the extracted TE was computed. Figure 1H shows the hexagonal distances using a random color for each distance value. To identify the tumor and stroma aspects of the IZ,

*ranks*are established. For any hexagon*h*,$rank\left(h\right)\underset{\xaf}{\underset{\xaf}{\text{def}}}\{\begin{array}{c}if\phantom{\rule{0.25em}{0ex}}class\left(h\right)=tumour,\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}rank=distance\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\\ if\phantom{\rule{0.25em}{0ex}}class\left(h\right)=stroma,\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}rank=\phantom{\rule{0.25em}{0ex}}-distance\phantom{\rule{0.25em}{0ex}}\\ if\phantom{\rule{0.25em}{0ex}}class\left(h\right)=TE,\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}rank=0\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\\ otherwise,\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}ignore\phantom{\rule{0.25em}{0ex}}hexagon\phantom{\rule{0.25em}{0ex}}h\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\end{array}$

(6)

The extracted TE has rank 0; inside the tumor, the rank is simply the positive distance from TE, and in stroma, it is the negative distance. The remaining tissue classes (background and tumor-background) are not included in further analyses.

The IZ, now consisting of TE with adjacent tumor and stroma tissue, can be defined for different choices of width. Here, IZ

_{w}denotes IZ consisting of rank interval $\left[-\frac{w}{2};\frac{w}{2}\right]$, where [ ] is rounding to nearest integer toward 0. Similarly, it is possible to define different widths of central TE by ranks. In this article only, TE = TE_{1}consisting of ${r}_{0}$ and TE_{3}consisting of ranks [−1; 1] are relevant. Unless otherwise mentioned, TE refers to TE_{1}. Figure 1I shows IZ_{9}; a nine-hexagon-wide IZ with tumor and stroma aspects labeled as red and green hexagons; their color intensity reflects distance to the yellow TE (with an implicit width of 1).### Computing Immunogradient Indicators in the IZ

Information about cell densities within and across the IZ can be computed by summarizing the hexagonal data values for each rank into

*rank quantities*(Table 2). The rank quantities form a collective interface profile that reveals how cell densities (and other features) vary inside and across the IZ. Examples of CD8^{+}cell density profiles are given for three tumors in Figure 2. An additional set of cases from both the BC and CRC cohorts illustrating the extracted IZ and immunogradient indicators in various tumor growth patterns and CD8^{+}density profiles are included in Supplemental Figure S1, Supplemental Figure S10, Supplemental Figure S11, Supplemental Figure S12, Supplemental Figure S2, Supplemental Figure S3, Supplemental Figure S4, Supplemental Figure S5, Supplemental Figure S6, Supplemental Figure S7, Supplemental Figure S8, Supplemental Figure S9.Table 2Summary of Grid Data Variables, Rank Quantities, and Profile Indicators

Value | Formula | Statistical notation |
---|---|---|

Grid data variables | ||

Positive cell counts | Stain No. | CD8 |

Areas | $Area=\frac{No.\phantom{\rule{0.25em}{0ex}}of\phantom{\rule{0.25em}{0ex}}pixels\phantom{\rule{0.25em}{0ex}}of\phantom{\rule{0.25em}{0ex}}tissue\phantom{\rule{0.25em}{0ex}}class}{pixels/\mu m}$ | T, S, and B |

Area fractions | $\frac{area\phantom{\rule{0.25em}{0ex}}of\phantom{\rule{0.25em}{0ex}}class}{hexagon\phantom{\rule{0.25em}{0ex}}area}$ | t, s, and b |

Cell density | $\frac{No.\phantom{\rule{0.25em}{0ex}}of\phantom{\rule{0.25em}{0ex}}positive\phantom{\rule{0.25em}{0ex}}stains}{T+S}$ | CD8_d |

Ranks and rank quantities | ||

${r}_{i}$ | Any rank in IZ of width w | ${r}_{i}\phantom{\rule{0.25em}{0ex}}\in I{Z}_{9}$ are $[{r}_{-4},{r}_{-3},\dots ,\phantom{\rule{0.25em}{0ex}}{r}_{3},{r}_{4}]$ |

$q\left({r}_{i}\right)$ | A statistic of hexagons in ${r}_{i}$ | Simple $q$‘s are mean and sd |

Mean of CD8 density in ${r}_{2}$ | CD8_mean (r_{2}) | |

Mean CD8 in tumor aspect | CD8_mean_T | |

Profile indicators | ||

CM | $CM\left(q\right)\phantom{\rule{0.25em}{0ex}}=\phantom{\rule{0.25em}{0ex}}\frac{{\sum}_{{r}_{i}}{r}_{i}\phantom{\rule{0.25em}{0ex}}q\left({r}_{i}\right)}{{\sum}_{{r}_{i}}q\left({r}_{i}\right)}$ | CD8_CM_mean |

ID | $ID\left(q\right)\phantom{\rule{0.25em}{0ex}}=\phantom{\rule{0.25em}{0ex}}\raisebox{1ex}{$q\left({r}_{-1}\right)$}\!\left/ \!\raisebox{-1ex}{$q\left({r}_{1}\right)$}\right.$ | CD8_ID_mean |

CM, center of mass; ID, immunodrop; IZ, interface zone.

To express the properties of directional variance (gradient) from stroma to tumor aspect, profile indicators can be calculated. Several indicators of cell density variance within and across the IZ were tested for statistical power to predict OS of the patients. Two cell density indicators were found to be the most powerful:

### Center of Mass

Center of mass (CM) was calculated as follows:

where ${r}_{i}$ indexes all ranks in the IZ, for example, ${r}_{i}\in \phantom{\rule{0.25em}{0ex}}\left[-4;4\right]$ for IZ

$CM\left(q\right)\phantom{\rule{0.25em}{0ex}}=\phantom{\rule{0.25em}{0ex}}\frac{{\sum}_{{r}_{i}}{r}_{i}\phantom{\rule{0.25em}{0ex}}q\left({r}_{i}\right)}{{\sum}_{{r}_{i}}q\left({r}_{i}\right)},$

(7)

where ${r}_{i}$ indexes all ranks in the IZ, for example, ${r}_{i}\in \phantom{\rule{0.25em}{0ex}}\left[-4;4\right]$ for IZ

_{9}, and $q\left({r}_{i}\right)$ denotes choice of rank quantity, for example, mean of CD8^{+}cell density. The CM calculates a coordinate along the horizontal axis under which one would have to support the profile for it to remain in perfect horizontal balance. Conversely, it can be understood as a measure of which part of the interface a biomarker gravitates (or has a positive gradient) toward.### ID

Immunodrop (ID) was calculated as follows:

$ID\left(q\right)\phantom{\rule{0.25em}{0ex}}=\phantom{\rule{0.25em}{0ex}}\frac{q\left({r}_{-1}\right)}{q\left({r}_{1}\right)}$

(8)

This indicator represents abrupt change in the cell density in the immediate vicinity of the TE and is calculated as the ratio of the quantities in ranks −1 and 1. Because an abrupt decrease in CD8

^{+}cell density from rank −1 to rank 1 was frequently observed, the indicator was named*immunodrop*.All data values

*v*, rank quantities*q*, and profile indicators are summarized in Table 2. The total computation time for the image analysis of the whole slide images and the extraction of the TE and IZ for computation of the final indicators was less than 2 hours per case. Summary statistics of CM, ID, and IZ densities are presented in Table 3 for BC and CRC using IZ width of nine and three hexagons, respectively.Table 3Summary Statistics of CD8 Density Indicators

Immunogradient indicators (CD8^{+} cell density, cells/mm^{2}) | BC | CRC | ||||
---|---|---|---|---|---|---|

Mean | Median | SD | Mean | Median | SD | |

CD8_CM_mean | −1.06 | −1.18 | 0.81 | −0.38 | −0.37 | 0.15 |

CD8_CM_sd | −1.16 | −1.25 | 0.54 | −0.30 | −0.29 | 0.12 |

CD8_ID_mean | 4.53 | 3.23 | 5.33 | 4.94 | 3.61 | 3.61 |

CD8_mean_S | 208.52 | 154.10 | 198.42 | 251.67 | 210.81 | 216.35 |

CD8_sd_S | 365.48 | 315.74 | 234.20 | 313.28 | 283.19 | 174.72 |

CD8_mean_TE | 147.64 | 109.12 | 164.55 | 178.01 | 123.30 | 194.66 |

CD8_sd_TE | 257.86 | 212.61 | 159.14 | 244.58 | 213.62 | 166.54 |

CD8_mean_T | 71.13 | 39.93 | 84.60 | 93.74 | 49.11 | 137.14 |

CD8_sd_T | 116.37 | 95.91 | 87.66 | 122.60 | 102.73 | 91.06 |

The interface zone (IZ) and tumor edge (TE) settings were IZ

_{9}(nine hexagons wide) and TE_{3}(3 hexagons wide) for the breast cancer (BC) cohort and IZ_{3}and TE_{1}for the colorectal cancer (CRC) cohort. CD8_mean and CD8_sd are summarized in the stroma aspect (S), TE, and tumor (T) aspects, respectively.*n*= 102 BC;*n*= 101 CRC.BC, breast cancer; CD8_CM_mean, center of mass for CD8 density by mean in ranks; CD8_CM_sd, center of mass for CD8 by variance in ranks; CD8_ID_mean, immunodrop of mean density; CRC, colorectal cancer.

### Statistical Analysis

After initial analyses, one CRC and five BC samples with tumor areas in the whole slide image <4.5 mm

^{2}were found too small to extract sufficient IZ data; therefore, these cases were excluded from further analyses. The distributions of CD8^{+}cell counts and densities revealed left asymmetry; therefore, logarithm-transformed values were used for parametric statistics (for the sake of readability, the prefix log is not used in the text or graphs). Summary statistics and distribution analyses were performed with significance tests based on one-way analysis of variance followed by Bonferroni post hoc tests for pairwise comparisons, and a two-sided Welch*t*-test was used to compare homogeneity of variances. Factor analysis was performed using the factoring method of principal component analysis; factors were retained based on the threshold of the eigenvalue of 1. General orthomax rotation of the initial factors was performed. The software tool Cutoff Finder version 2.1 (Charité University, Berlin, Germany) was used to determine a cutoff value for each indicator to test univariate and multiple predictions of OS. The OS distributions for the patients were estimated using the Kaplan-Meier method. The log-rank test was used to assess the statistical significance of differences between the stratified groups.Regression analysis of the OS data was performed using a Cox proportional hazards regression model obtained by a stepwise likelihood ratio test procedure to analyze the independent prognostic significance of TIL density indicators in the context of clinicopathologic variables. For both cohorts, the best representable subset of variables was found using leave-one-out cross-validation. To investigate the combined power of CD8 density and CD8 gradient variables, a second model was obtained by replacing the immunogradient CD8 variables with the aggregated IZ CD8

^{+}cell response factor; see models 1 and 2 in Table 4.Table 4Multivariate Cox Proportional Hazards Regression Analysis with the CD8

^{+}Cell Density (CD8) Indicators and Their Factor Scores in IZ_{9}with TE_{3}for the Cohort of Patients with BC and in IZ_{3}with TE_{1}for the Cohort of Patients with CRCBC | CRC | ||||
---|---|---|---|---|---|

Model 1 (LR = 18.20, P = 0.0004) | HR (95% CI) | P value | Model 1 (LR = 22.54, P < 0.0001) | HR (95% CI) | P value |

pN status (N0 vs N1–3) | 2.61 (1.03–6.62) | 0.0428 | Age group by median | 2.04 (1.15–4.09) | 0.0453 |

CD8_CM_mean | 0.33 (0.13–0.84) | 0.0207 | pT stage (T1–2 vs T3–4) | 6.22 (1.48–26.21) | 0.0128 |

CD8_mean_T | 0.05 (0.01–0.28) | 0.0005 | CD8_CM_mean | 0.39 (0.20–0.77) | 0.0071 |

Model 2 (LR = 15.36, P = 0.0005) | Model 2 (LR = 21.58, P < 0.0001) | ||||

pN status (N0 vs N1–3) | 2.40 (1.01–5.72) | 0.0488 | Age group by median age | 2.24 (1.12–4.50) | 0.0235 |

Aggregated IZ CD8^{+} cell response factor | 0.23 (0.1–0.53) | 0.0007 | pT stage (T1–2 vs T3–4) and aggregated IZ CD8^{+} cell response factor | 5.25 (1.22–22.61) and 0.41 (0.19–0.87) | 0.0260 and 0.0196 |

*n*= 102 BC;

*n*= 101 CRC.

BC, breast cancer; CD8_CM_mean, center of mass for mean density; CD8_CM_sd, center of mass for density variance, CD8_ID_mean, immunodrop of mean CD8 density; CRC, colorectal cancer; HR, hazard ratio; IZ, interface zone (numbers indicate number of hexagons); LR, likelihood ratio; S, stroma; T, tumor; TE, tumor edge.

The TE extraction and computation of hexagonal data variables were implemented in C++ (g++ 7.3.8), using libtiff version 5.2.4 (

*https://www.libtiff.org*) and Boost version 1.67 (*https://www.boost.org*). SAS software version 9.4 (SAS Institute Inc., Cary, NC) was used to perform statistical analyses. Plots were produced using R software version 3.4.4 (R Foundation for Statistical Computing, Vienna, Austria).Data tables with raw data and TE and immunogradient results are available for both the BC and CRC cohorts, including all subsampled hexagonal grids, extracted hexagonal variables, pathologic variables, and survival data. The tables can be downloaded from Vilnius University research storage (

*https://www.midas.lt/action/resources/790c85e2-f6e5-4db3-9d55-69bb6527b6a3*, last accessed January 30, 2020).## Results

### Summary Statistics of CD8^{+} Cell Density Indicators in the IZ

A total of 102 BC and 101 CRC samples were used to explore the prognostic power of the CD8

^{+}cell densities in the stroma aspect, TE, and tumor aspect of the IZ and the density profiles across it. The density profile indicators included CM and ID. IZ widths of three, five, seven, and nine hexagons, indicated as IZ_{3}, IZ_{5}, IZ_{7}, and IZ_{9}, respectively, were tested with TE width of both 1 (TE_{1}) and 3 (TE_{3}). The most significant prognostic stratifications were achieved in IZ_{9}with TE_{3}for BC and in IZ_{3}with TE_{1}for the CRC samples. The relevant summary statistics of the indicators are presented in Table 3.The distribution of CD8

^{+}cells within the IZ was similar in both cancers: the densities were highest and most dispersed in the stroma aspect, less abundant and dispersed within the TE (both TE_{1}and TE_{3}), and lowest and less dispersed in the tumor aspect of the IZ (*P*< 0.05) (Table 3). Importantly, in both patient cohorts, a trend of IZ CD8^{+}cell density decrease toward the tumor aspect was observed (Figure 3).### Factor Analysis of the IZ CD8^{+} Cell Density Indicators

Two orthogonally independent factors of IZ CD8

^{+}cell density indicators were extracted for both the BC and CRC cohort (Figure 4). Altogether, the two factors explained 90.92% and 95.64% of the variance in the data set in BC and CRC, respectively, and the factor patterns were similar between the cohorts. Factor 1 was characterized by strong positive loadings of the variables that reflected the level of CD8^{+}density within all IZ aspects and was therefore interpreted as the IZ CD8^{+}density level factor, also referred to as CD8^{+}density factor. Factor 2 was characterized by strong positive loadings of CM for both means and SDs of the CD8^{+}densities and strong negative loading of the ID indicator and thus represents the CD8^{+}density change across the IZ. Higher values of factor 2 scores reflect increasing CD8^{+}densities toward the tumor aspect of the interface; therefore, factor 2 can be interpreted as the IZ CD8^{+}density gradient, more specifically IZ CD8^{+}density gradient toward tumor, also referred to as CD8^{+}gradient factor. Because both factors may represent important (although independent) properties of the immune response, the sum of factor 1 and factor 2 scores were computed as aggregated IZ CD8^{+}cell response factor to test their combined prognostic power.### Prognostic Value of the IZ CD8^{+} Cell Density and Gradient Indicators

Patient OS probability stratifications based on univariate Kaplan-Meier analyses are presented in Figure 5, and prognostic models obtained by multiple Cox proportional hazards regression analyses are presented in Table 4. Two models were investigated for both cohorts: model 1 used the set of IZ CD8

^{+}cell density indicators along with the clinicopathologic variables, and model 2 replaced the IZ CD8^{+}cell density indicators with the aggregated IZ CD8^{+}factor score. The models in both patient cohorts were essentially similar, but the indicators derived from a wider IZ_{9}with TE_{3}were more informative for the patients in the BC cohort, whereas a narrow IZ_{3}with TE_{1}provided better prognostic stratification of the patients in the CRC cohort. In the BC cohort, better OS was predicted by the CM for the CD8^{+}mean density and mean of the CD8^{+}density in the tumor aspect of the IZ (model 1) as well as the aggregated IZ CD8^{+}cell response factor (model 2) in the context of worse OS predicted by pN status. For the patients in the CRC cohort, model 1 included the CM for the CD8^{+}mean density indicator to predict better OS; in model 2, the aggregated IZ CD8^{+}cell response factor predicted better OS in the context of worse OS predicted by the patient age and pT status. Importantly, the models in both BC and CRC revealed strong independent prognostic value of the IZ CD8^{+}density indicators measuring the density gradient toward the tumor aspect alone or in combination with absolute density of CD8^{+}infiltrate in the IZ expressed as the aggregated IZ CD8^{+}cell response factor score.## Discussion

Our study found that automated detection of the tumor-stroma IZ in digital microscopy images of tumor tissues with extraction of immune cell density profiles provides novel indicators to assess antitumor immune response in the very frontline of tumor-host interaction. Although the IZ in our method is similar to the concept of IM commonly used in pathology and defined in previous DIA studies,

^{22}

- Hendry S.
- Salgado R.
- Gevaert T.
- Russell P.A.
- John T.
- Thapa B.
- et al.

Assessing tumor-infiltrating lymphocytes in solid tumors: a practical review for pathologists and proposal for a standardized method from the international immunooncology biomarkers working group: part 1: assessing the host immune response, TILs in invasive breast carcinoma and ductal carcinoma in situ, metastatic tumor deposits and areas for further research.

*Adv Anat Pathol.*2017; 24: 235-251

^{,}^{23}

- Hendry S.
- Salgado R.
- Gevaert T.
- Russell P.A.
- John T.
- Thapa B.
- et al.

Assessing tumor-infiltrating lymphocytes in solid tumors: a practical review for pathologists and proposal for a standardized method from the international immuno-oncology biomarkers working group: part 2: TILs in melanoma, gastrointestinal tract carcinomas, non-small cell lung carcinoma and mesothelioma, endometrial and ovarian carcinomas, squamous cell carcinoma of the head and neck, genitourinary carcinomas, and primary brain tumors.

*Adv Anat Pathol.*2017; 24: 311-335

^{,}^{33}

- Galon J.
- Costes A.
- Sanchez-Cabo F.
- Kirilovsky A.
- Mlecnik B.
- Lagorce-Pages C.
- Tosolini M.
- Camus M.
- Berger A.
- Wind P.
- Zinzindohoue F.
- Bruneval P.
- Cugnenc P.H.
- Trajanoski Z.
- Fridman W.H.
- Pages F.

Type, density, and location of immune cells within human colorectal tumors predict clinical outcome.

*Science.*2006; 313: 1960-1964

^{, }^{34}

- Lechner A.
- Schlosser H.
- Rothschild S.I.
- Thelen M.
- Reuter S.
- Zentis P.
- Shimabukuro-Vornhagen A.
- Theurich S.
- Wennhold K.
- Garcia-Marquez M.
- Tharun L.
- Quaas A.
- Schauss A.
- Isensee J.
- Hucho T.
- Huebbers C.
- von Bergwelt-Baildon M.
- Beutner D.

Characterization of tumor-associated T-lymphocyte subsets and immune checkpoint molecules in head and neck squamous cell carcinoma.

*Oncotarget.*2017; 8: 44418-44433

^{, }^{35}

- Bordry N.
- Broggi M.A.S.
- de Jonge K.
- Schaeuble K.
- Gannon P.O.
- Foukas P.G.
- Danenberg E.
- Romano E.
- Baumgaertner P.
- Fankhauser M.
- Wald N.
- Cagnon L.
- Abed-Maillard S.
- Hajjami H.M.
- Murray T.
- Ioannidou K.
- Letovanec I.
- Yan P.
- Michielin O.
- Matter M.
- Swartz M.A.
- Speiser D.E.

Lymphatic vessel density is associated with CD8(+) T cell infiltration and immunosuppressive factors in human melanoma.

*Oncoimmunology.*2018; 7: e1462878

^{, }essential differences should be emphasized: i) the IZ detection is completely data driven and independent of arbitrary visual evaluation; ii) the transition from stroma to tumor is achieved by spatial ranks within the IZ, which enables computation of immune cell density indicators across the IZ to reflect the stroma-to-tumor gradient of the infiltrates; and iii) the grid subsampling allows consideration of the variance of the indicators within each rank of the IZ. The set of immunogradient indicators provides strong and independent prognostic stratification in the BC and CRC patient cohorts. Importantly, the proposed method is flexible and can be adapted for various assays of inflammatory and immune response in tumor and nontumor pathology.Our method uses systematic grid-based subsampling with explicit data-driven rules to extract and rank the tumor-stroma transition area from digital images of microscopy slides; in essence, the probability of interfaceness is computed. In comparison, current DIA platforms, for instance, Indica Labs HALO and QuPath (Queen's University, Belfast, Northern Ireland), only provide the option to select a fixed-width margin along a hand-drawn line around the region of interest. In addition, the Immunoscore procedure presumably depends on a manual annotation of IM, and most immune profiling studies were based on manual delineation of IM from the core tumor. A manual IM annotation step is time-consuming and may cause bias because of inevitable intraobserver and interobserver variability. For example, manual IM delineation discrepancies are very likely in tumors with irregular growth patterns, including tumor budding features.

^{37}

^{,}^{38}

^{12}

^{,}^{13}

^{,}^{33}

- Galon J.
- Costes A.
- Sanchez-Cabo F.
- Kirilovsky A.
- Mlecnik B.
- Lagorce-Pages C.
- Tosolini M.
- Camus M.
- Berger A.
- Wind P.
- Zinzindohoue F.
- Bruneval P.
- Cugnenc P.H.
- Trajanoski Z.
- Fridman W.H.
- Pages F.

Type, density, and location of immune cells within human colorectal tumors predict clinical outcome.

*Science.*2006; 313: 1960-1964

^{,}^{36}

^{39}

- Kather J.N.
- Suarez-Carmona M.
- Charoentong P.
- Weis C.A.
- Hirsch D.
- Bankhead P.
- Horning M.
- Ferber D.
- Kel I.
- Herpel E.
- Schott S.
- Zornig I.
- Utikal J.
- Marx A.
- Gaiser T.
- Brenner H.
- Chang-Claude J.
- Hoffmeister M.
- Jager D.
- Halama N.

Topography of cancer-associated immune cells in human solid tumors.

*Elife.*2018; 7: e36967

^{, }^{40}

- Steele K.E.
- Tan T.H.
- Korn R.
- Dacosta K.
- Brown C.
- Kuziora M.
- Zimmermann J.
- Laffin B.
- Widmaier M.
- Rognoni L.
- Cardenes R.
- Schneider K.
- Boutrin A.
- Martin P.
- Zha J.
- Wiestler T.

Measuring multiple parameters of CD8+ tumor-infiltrating lymphocytes in human cancers by image analysis.

*J Immunother Cancer.*2018; 6: 20

^{, }^{41}

^{42}

^{,}In contrast, our approach is entirely data driven and independent of arbitrary visual evaluation. Instead of providing only a crude high-level outline between tumor and stroma regions or seeking to delicately delineate every tumor cell conglomerate, our method evaluates grid-subsampled tissue areas to represent the tumor-host interface and its spatial vicinity to the tumor or host aspect of the interface. Consequently, it is less dependent on a specific tumor growth pattern and captures the IZ in tumors with pushing and/or infiltrative IM, even in cases with very irregular infiltrating growth patterns and budding phenomena.Besides detecting the TE, the proposed method provides the important advantage of adjustable width of the IZ and spatial ranking the stroma and tumor aspects within the IZ. This method enables computation of IZ profile indicators, in particular immune cell density indicators that reflect directional properties (gradient) of the infiltrates. The flexibility in IZ ranking and compartmentalization of tissue aspects are of particular interest because the choice of IM width is not consistent throughout individual studies. The international immuno-oncology biomarkers working group defined IM as a region of 1-mm width centered on the border between tumor and immune cells as proposed by other investigators. Other studies reported prognostic significance of high TIL density for different widths of the IM region. For example, Lechner et al used an IM that ranged from 50 μm within the tumor to 300 μm outside the tumor border in head and neck squamous cell carcinoma; Bordry et al defined the IM as a tumor region of approximately 400-μm width between the tumor and the reticular dermis in primary melanoma. Hermitte defined the IM as a region with a span of 360 μm into the healthy tissue and 360 μm into the tumor in CRC. Another study found a survival benefit for patients with BC and high TIL densities within the 0- to 10-μm distance from a manually outlined tumor border. These findings suggest that only TIL in close contact with the tumor cells deliver direct and potent antitumor effects.Type, density, and location of immune cells within human colorectal tumors predict clinical outcome. Consequently, TIL densities measured in narrow IM regions can provide more precise prognostic markers compared with those based on wider IM or global tumor or stroma sampling. Our grid-based technique allows flexibility in analysis settings to adjust the resolution of tissue area sampling and IZ sampling width. We tested various analysis settings in our study with similar results. The best prognostic models were achieved with relatively wide IZ in BC and narrow IZ in CRC, which may be explained by the need to sample larger IZ in hormone receptor–positive BC with relatively sparse TILs to obtain sufficient data for the immunogradient indicators.

^{22}

- Hendry S.
- Salgado R.
- Gevaert T.
- Russell P.A.
- John T.
- Thapa B.
- et al.

Assessing tumor-infiltrating lymphocytes in solid tumors: a practical review for pathologists and proposal for a standardized method from the international immunooncology biomarkers working group: part 1: assessing the host immune response, TILs in invasive breast carcinoma and ductal carcinoma in situ, metastatic tumor deposits and areas for further research.

*Adv Anat Pathol.*2017; 24: 235-251

^{,}^{23}

- Hendry S.
- Salgado R.
- Gevaert T.
- Russell P.A.
- John T.
- Thapa B.
- et al.

Assessing tumor-infiltrating lymphocytes in solid tumors: a practical review for pathologists and proposal for a standardized method from the international immuno-oncology biomarkers working group: part 2: TILs in melanoma, gastrointestinal tract carcinomas, non-small cell lung carcinoma and mesothelioma, endometrial and ovarian carcinomas, squamous cell carcinoma of the head and neck, genitourinary carcinomas, and primary brain tumors.

*Adv Anat Pathol.*2017; 24: 311-335

^{,}^{33}

- Galon J.
- Costes A.
- Sanchez-Cabo F.
- Kirilovsky A.
- Mlecnik B.
- Lagorce-Pages C.
- Tosolini M.
- Camus M.
- Berger A.
- Wind P.
- Zinzindohoue F.
- Bruneval P.
- Cugnenc P.H.
- Trajanoski Z.
- Fridman W.H.
- Pages F.

Type, density, and location of immune cells within human colorectal tumors predict clinical outcome.

*Science.*2006; 313: 1960-1964

^{34}

- Lechner A.
- Schlosser H.
- Rothschild S.I.
- Thelen M.
- Reuter S.
- Zentis P.
- Shimabukuro-Vornhagen A.
- Theurich S.
- Wennhold K.
- Garcia-Marquez M.
- Tharun L.
- Quaas A.
- Schauss A.
- Isensee J.
- Hucho T.
- Huebbers C.
- von Bergwelt-Baildon M.
- Beutner D.

Characterization of tumor-associated T-lymphocyte subsets and immune checkpoint molecules in head and neck squamous cell carcinoma.

*Oncotarget.*2017; 8: 44418-44433

^{35}

- Bordry N.
- Broggi M.A.S.
- de Jonge K.
- Schaeuble K.
- Gannon P.O.
- Foukas P.G.
- Danenberg E.
- Romano E.
- Baumgaertner P.
- Fankhauser M.
- Wald N.
- Cagnon L.
- Abed-Maillard S.
- Hajjami H.M.
- Murray T.
- Ioannidou K.
- Letovanec I.
- Yan P.
- Michielin O.
- Matter M.
- Swartz M.A.
- Speiser D.E.

Lymphatic vessel density is associated with CD8(+) T cell infiltration and immunosuppressive factors in human melanoma.

*Oncoimmunology.*2018; 7: e1462878

^{36}

^{44}

^{33}

- Galon J.
- Costes A.
- Sanchez-Cabo F.
- Kirilovsky A.
- Mlecnik B.
- Lagorce-Pages C.
- Tosolini M.
- Camus M.
- Berger A.
- Wind P.
- Zinzindohoue F.
- Bruneval P.
- Cugnenc P.H.
- Trajanoski Z.
- Fridman W.H.
- Pages F.

*Science.*2006; 313: 1960-1964

Several recent studies applied common DIA morphologic operations to automatically delineate tumor cell clusters and stroma. Harder et al performed comprehensive screening of immune contexture signatures to predict disease recurrence in 90 patients with prostate cancer. They used morphologic operations to smoothen and dilate the tumor regions into the stroma. Subsequently, they defined the tumor border as a region that reaches equally far to the inside and outside of the tumor region (total width of 112.5 μm) as well as two TME zones inside the stroma, effectively extracting three transition zones. A large set of density and distance indicators were computed for CD3 investigated spatial aspects of CD8

^{20}

^{+}, CD8^{+}, CD34^{+}, CD68^{+}, and CD163^{+}cells*within*each transition zone. Their best prognostic indicators expressed the ratio of different biomarkers within the transition zones (best candidate being the CD8/CD34 ratio); however, none of their selected significant indicators in prostate cancer represented a biomarker density profile across the transition zones (gradient) properties. Li et al^{45}

^{+}cell infiltration into tumor cell clusters in 28 patients with triple-negative BC. Although they performed manual delineation of the IM, their study also included an automated tumor cell boundary detection based on pan-cytokeratin fluorescence images. This method allowed them to measure CD8^{+}pixel ratios across the boundary to explore potential factors of TIL motility by mathematical modeling, and they concluded that chemokines rather than dense collagen fibers prevent TILs from entering the tumor cell clusters. In both studies, the tumor-stroma delineation within the TME goes beyond the classic IM concept, which is similar to our definition of the IZ. However, the IZ is based on grid-based subsampling of the relevant tissue classes and is less dependent on the individual tumor growth pattern and precision in segmenting individual tumor cells; also, it is not forced to be a ribbon of a fixed width in a two-dimensional image of an actual three-dimensional interface.We explored the prognostic value of the immune cell density profiles in the IZ in the BC and CRC cohorts, which led to a set of immunogradient indicators that enabled strong and independent prognostic stratification in both patient cohorts. We tested various analysis settings to generate the most informative prognostic models for each of the cohorts and found that a relatively narrow IZ generated more powerful prognostic models for CRC, whereas BC benefited from broader IZ data collection. Remarkably, in both the BC and CRC cohorts, it was found that the indicators revealing a gravitation of CD8

^{+}cells toward the tumor aspect within the IZ were independent predictors of better OS. In particular, the CM indicator, which cumulatively reflects the gradient of CD8^{+}cells toward the tumor within the IZ, was strongly associated with better OS in both patient cohorts. Similarly, the ID indicator, which reflects an abrupt decrease in CD8^{+}cell density in the tumor aspect of the IZ in the nearest vicinity of the TE, was associated with worse OS. By definition, the ID is independent of the width of the IZ because it always represents the narrowest (IZ_{3}) sampling. Therefore, this indicator may perform better in cases of low tumor content or core biopsy samples in which tumor aspect ranking can be limited.The indicators that reflect the spatial shift of the IZ CD8

^{+}cells toward the tumor revealed added prognostic value to the density of IZ CD8^{+}cells per se. TIL densities and their ratios in the global or remaining tumor and stroma areas were less informative (data not shown). Both IZ CD8^{+}density level and IZ CD8^{+}density gradient factor scores, although orthogonally independent, were significantly associated in univariate analyses with better OS in both the BC and CRC patients (Tables 4 and 5, and Figure 5, A–F). Furthermore, the aggregated IZ CD8^{+}cell response factor, obtained by the sum of the two factor scores, was an independent prognostic marker of OS in both patient cohorts. This finding indicates that combining absolute density and gradient aspects of the immune cell infiltrates across the IZ may provide added prognostic value of the image-based biomarker of antitumor immune response. Indeed, the aggregated IZ CD8^{+}factor allowed dichotomization of the patients with CRC into prognostic groups with 5-year OS probability of 82% and 58% (Figure 5H), which is similar to the patients with CRC dichotomized by Immunoscore based on CD3 and CD8 IHC as reported by Pages et al.^{46}

Table 5Univariate Kaplan-Meier Analysis with the CD8

^{+}Cell Density (CD8) Indicators and Their Factor Scores in IZ_{9}with TE_{3}for the Cohort of Patients with BC and in IZ_{3}with TE_{1}for the Cohort of Patients with CRC,IZ indicators and clinicopathologic variables | BC | CRC | ||
---|---|---|---|---|

HR (95% CI) | P value | HR (95% CI) | P value | |

CD8_CM_mean | 0.28 (0.11–0.71) | 0.0044 | 0.33 (0.16–0.65) | 0.0008 |

CD8_CM_sd | 0.39 (0.13–1.14) | 0.0740 | 0.41 (0.20–0.81) | 0.0077 |

CD8_mean_S | 0.43 (0.15–1.27) | 0.1200 | 0.11 (0.01–0.80) | 0.0080 |

CD8_sd_S | 0.40 (0.14–1.18) | 0.0860 | 0.41 (0.20–0.86) | 0.0140 |

CD8_mean_TE | 0.38 (0.16–0.91) | 0.0250 | 0.27 (0.09–0.76) | 0.0076 |

CD8_sd_TE | 0.34 (0.12–1.01) | 0.0420 | 0.44 (0.21–0.93) | 0.0260 |

CD8_mean_T | 0.21 (0.09–0.53) | 0.0002 | 0.35 (0.18–0.69) | 0.0014 |

CD8_sd_T | 0.26 (0.10–0.66) | 0.0024 | 0.39 (0.20–0.76) | 0.0044 |

CD8^{+} density factor | 0.31 (0.11–0.84) | 0.0140 | 0.43 (0.20–0.96) | 0.0330 |

CD8^{+} gradient factor | 0.34 (0.14–0.84) | 0.0140 | 0.35 (0.18–0.70) | 0.0018 |

Aggregated IZ CD8 factor | 0.23 (0.10–0.55) | 0.0003 | 0.32 (0.15–0.67) | 0.0015 |

G stage (G1–2 vs G3) | 1.40 (0.59–3.33) | 0.4523 | 0.95 (0.85–1.06) | 0.3379 |

T stage (T1–2 vs T3–4) | 1.19 (0.52–2.75) | 0.6786 | 6.47 (1.55–27.10) | 0.0106 |

N status (N0 vs N1–3) | 2.30 (0.97–5.49) | 0.0598 | 1.01 (1.00–1.02) | 0.0214 |

Age group by median age | 2.63 (1.07–6.46) | 0.0347 | 1.86 (0.93–3.71) | 0.0793 |

Sex | Not applicable | 1.56 (0.79–3.05) | 0.1974 |

*n*= 102 BC;

*n*= 101 CRC.

BC, breast cancer; CD8_CM_mean, center of mass for mean density; CD8_CM_sd, center of mass for density variance, CD8_ID_mean, immunodrop of mean CD8 density; CRC, colorectal cancer; HR, hazard ratio; IZ, interface zone (numbers indicate number of hexagons); S, stroma; T, tumor; TE, tumor edge.

Remarkably, a good prognostic stratification was achieved in early hormone receptor–positive BC, which is a relatively indolent disease without an established role for immune response. Furthermore, a prominent time-dependent effect was discovered in the prognostic stratification of the patients with BC (Figure 5): >92% of the patients were alive after 5 years from surgery in both groups. However, the OS curves diverged sharply at this infliction point to end at 87% and 55% OS probability after 10 years in the patients with high and low aggregated IZ CD8

^{+}cell response factor scores, respectively. Of note, this time-dependent effect is revealed only with the gradient-type indicators (Figure 5, C, E, and G) but not in the stratification by the CD8^{+}cell mean density on tumor aspect of the IZ (Figure 5A). This finding emphasizes the potential of the immunogradient indicators, obtained from a surgical excision sample and stained for a single CD8^{+}IHC, to predict long-term prognosis of patients with this relatively well-managed disease. Furthermore, it implies that antitumor immune response properties, which potentially determine the long-term outcome of the disease, are encoded in the TME and captured by the immunogradient assay. A biologic interpretation of this phenomenon leads to the hypothesis that the individual immune response properties are established in the early stages of BC and are critical to reach full recovery without remaining dormant cancer cells under standard therapies. Subsequently, this finding raises the question of whether patients with hormone receptor–positive BC with low immunogradient scores could benefit from immunotherapy at some point.In this study, we present proof of concept and a detailed description of a new method that provides prognostic value in two independent but relatively small patient cohorts. Large-scale studies are needed to further optimize the analysis settings and to establish clinical utility of the immunogradient indicators for various cancer types and prognostic tasks. In addition, direct, in-depth comparisons to existing methods should be performed in appropriately designed studies with all methods applied to the same patient cohorts.

In conclusion, our study demonstrates the potential of high-capacity, automated, and data-driven sampling of tumor-stroma interfaces in microscopy images to extract novel, clinically useful indicators of antitumor immune response. The IZ immunogradient enables strong and independent prognostic stratification of patients with BC and CRC, with prominent divergence of OS probability in patients with hormone receptor—positive BC 5 years after surgery.

## Acknowledgments

A.R. and A.L. constructed the interface zone extraction rules; A.R., D.Z., and R.A. implemented the interface zone extraction rules; A.R., D.Z., A.N., R.A., and A.L. performed statistical analyses; A.R., in collaboration with D.Z., A.N., and A.L., drafted essential parts of the manuscript; A.R., D.Z., A.N., V.O., T.P., and A.L. reviewed and edited the manuscript; all authors participated in conception and design of the study, data analysis, and critical revision and approval of the final manuscript.

## Supplemental Data

## References

- The cancer genome.
*Nature.*2009; 458: 719-724 - Unraveling the genetics of cancer: genome sequencing and beyond.
*Annu Rev Genomics Hum Genet.*2011; 12: 407-430 - Hallmarks of cancer: the next generation.
*Cell.*2011; 144: 646-674 - The tumor microenvironment and its role in promoting tumor growth.
*Oncogene.*2008; 27: 5904-5912 - Implications of the tumor immune microenvironment for staging and therapeutics.
*Mod Pathol.*2018; 31: 214-234 - State-of-the-art of profiling immune contexture in the era of multiplexed staining and digital analysis to study paraffin tumor tissues.
*Cancers (Basel).*2019; 11: E247 - Inflammation and cancer.
*Nature.*2002; 420: 860-867 - The immune contexture in human tumours: impact on clinical outcome.
*Nat Rev Cancer.*2012; 12: 298-306 - Increased p53 mutation load in noncancerous colon tissue from ulcerative colitis: a cancer-prone chronic inflammatory disease.
*Cancer Res.*2000; 60: 3333-3337 - T lymphocyte subsets in cancer immunity: friends or foes.
*J Leukoc Biol.*2019; 105: 243-255 - Adoptive cell transfer as personalized immunotherapy for human cancer.
*Science.*2015; 348: 62-68 - Cancer classification using the Immunoscore: a worldwide task force.
*J Transl Med.*2012; 10: 205 - Towards the introduction of the 'Immunoscore' in the classification of malignant tumours.
*J Pathol.*2014; 232: 199-209 - Prognostic value of tumor-infiltrating lymphocyte density assessed using a standardized method based on molecular subtypes and adjuvant chemotherapy in invasive breast cancer.
*Ann Surg Oncol.*2018; 25: 937-946 - Tumour-infiltrating lymphocytes and prognosis in different subtypes of breast cancer: a pooled analysis of 3771 patients treated with neoadjuvant therapy.
*Lancet Oncol.*2018; 19: 40-50 - The cancer immunogram as a framework for personalized immunotherapy in urothelial cancer.
*Eur Urol.*2019; 75: 435-444 - An immunogram for the cancer-immunity cycle: towards personalized immunotherapy of lung cancer.
*J Thorac Oncol.*2017; 12: 791-803 - Cancer immunology: the "cancer immunogram".
*Science.*2016; 352: 658-660 - Prognostic value of the neo-immunoscore in renal cell carcinoma.
*Front Oncol.*2019; 9: 439 - Tissue phenomics for prognostic biomarker discovery in low- and intermediate-risk prostate cancer.
*Sci Rep.*2018; 8: 4470 - Caie PD: automated analysis of lymphocytic infiltration, tumor budding, and their spatial relationship improves prognostic accuracy in colorectal cancer.
*Cancer Immunol Res.*2019; 7: 609-620 - Assessing tumor-infiltrating lymphocytes in solid tumors: a practical review for pathologists and proposal for a standardized method from the international immunooncology biomarkers working group: part 1: assessing the host immune response, TILs in invasive breast carcinoma and ductal carcinoma in situ, metastatic tumor deposits and areas for further research.
*Adv Anat Pathol.*2017; 24: 235-251 - Assessing tumor-infiltrating lymphocytes in solid tumors: a practical review for pathologists and proposal for a standardized method from the international immuno-oncology biomarkers working group: part 2: TILs in melanoma, gastrointestinal tract carcinomas, non-small cell lung carcinoma and mesothelioma, endometrial and ovarian carcinomas, squamous cell carcinoma of the head and neck, genitourinary carcinomas, and primary brain tumors.
*Adv Anat Pathol.*2017; 24: 311-335 - The grading of rectal cancer: historical perspectives and a multivariate analysis of 447 cases.
*Histopathology.*1986; 10: 437-459 - A new prognostic classification of rectal cancer.
*Lancet.*1987; 1: 1303-1306 - The tumor microenvironment and Immunoscore are critical determinants of dissemination to distant metastasis.
*Sci Transl Med.*2016; 8: 327ra26 - Data-driven discovery of immune contexture biomarkers.
*Front Oncol.*2018; 8: 627 - Ki67/SATB1 ratio is an independent prognostic factor of overall survival in patients with early hormone receptor-positive invasive ductal breast carcinoma.
*Oncotarget.*2015; 6: 41134-41145 - Immunohistochemistry profiles of breast ductal carcinoma: factor analysis of digital image analysis data.
*Diagn Pathol.*2012; 7: 27 - Tri-directional gradient operators for hexagonal image processing.
*J Vis Commun Image Represent.*2016; 38: 614-626 - Edge detection in a hexagonal-image processing framework.
*Image Vis Comput.*2001; 19: 1071-1081 - Cutoff Finder: a comprehensive and straightforward Web application enabling rapid biomarker cutoff optimization.
*PLoS One.*2012; 7: e51862 - Type, density, and location of immune cells within human colorectal tumors predict clinical outcome.
*Science.*2006; 313: 1960-1964 - Characterization of tumor-associated T-lymphocyte subsets and immune checkpoint molecules in head and neck squamous cell carcinoma.
*Oncotarget.*2017; 8: 44418-44433 - Lymphatic vessel density is associated with CD8(+) T cell infiltration and immunosuppressive factors in human melanoma.
*Oncoimmunology.*2018; 7: e1462878 - Biomarkers immune monitoring technology primer: Immunoscore(R) Colon.
*J Immunother Cancer.*2016; 4: 57 - Quantifying immune cell distribution in the tumor microenvironment using HALO® spatial analysis tools.Indica Labs, 2016
- Open source software for digital pathology image analysis.
*Sci Rep.*2017; 7: 16878 - Topography of cancer-associated immune cells in human solid tumors.
*Elife.*2018; 7: e36967 - Measuring multiple parameters of CD8+ tumor-infiltrating lymphocytes in human cancers by image analysis.
*J Immunother Cancer.*2018; 6: 20 - Computer-assisted stereology and automated image analysis for quantification of tumor infiltrating lymphocytes in colon cancer.
*Diagn Pathol.*2017; 12: 65 - Tumour budding at invasive margins and outcome in colorectal cancer.
*Colorectal Dis.*2008; 10: 41-47 - Automatic evaluation of tumor budding in immunohistochemically stained colorectal carcinomas and correlation to clinical outcome.
*Diagn Pathol.*2018; 13: 64 - Detailed resolution analysis reveals spatial T cell heterogeneity in the invasive margin of colorectal cancer liver metastases associated with improved survival.
*Oncoimmunology.*2017; 6: e1286436 - Infiltration of CD8(+) T cells into tumor cell clusters in triple-negative breast cancer.
*Proc Natl Acad Sci U S A.*2019; 116: 3678-3687 - International validation of the consensus Immunoscore for the classification of colon cancer: a prognostic and accuracy study.
*Lancet.*2018; 391: 2128-2139

## Article Info

### Publication History

Published online: March 17, 2020

Accepted:
January 28,
2020

### Footnotes

Supported by European Social Fund project 09.3.3-LMT-K-712-01-0139 under grant agreement with the Research Council of Lithuania.

Disclosures: Vilnius University and A.L., A.R., D.Z., A.N., and R.A. submitted a Lithuanian patent (LT2019 509) for the invented immunogradient indicators.

### Identification

### Copyright

© 2020 American Society for Investigative Pathology. Published by Elsevier Inc.

### User License

Creative Commons Attribution – NonCommercial – NoDerivs (CC BY-NC-ND 4.0) | How you can reuse

Elsevier's open access license policy

Creative Commons Attribution – NonCommercial – NoDerivs (CC BY-NC-ND 4.0)

## Permitted

### For non-commercial purposes:

- Read, print & download
- Redistribute or republish the final article
- Text & data mine
- Translate the article (private use only, not for distribution)
- Reuse portions or extracts from the article in other works

## Not Permitted

- Sell or re-use for commercial purposes
- Distribute translations or adaptations of the article

Elsevier's open access license policy