Distinct Tumor-Immune Ecologies in Patients with Lung Cancer Predict Progression and Define a Clinical Biomarker of Therapy Response.

Introduction
Despite major advances in elucidating the molecular mechanisms of lung cancer, and improvements in prognosis, diagnosis and treatment, lung cancer remains the leading cause of cancer death worldwide (1). Non-Small Cell Lung Cancer (NSCLC) accounts for about 85% of all cases (2). Disease progression and treatment response in NSCLC vary widely among patients. As targeted molecular therapies, immunotherapies, and combination therapies have become central in managing patients with NSCLC, the ability to identify which patients are appropriate for each therapeutic approach has become even more important. This has accelerated the research to develop biomarkers for identifying patients who would likely respond to therapy. Additionally, biomarkers monitor disease development and help predict disease outcomes and treatment response (3,4). In the case of immunotherapy, existing treatment response biomarkers (PD-L1) are useful, but there are well recognized limitations (3), such as, the spatial and temporal heterogeneity of PD-L1 expression, the dynamic nature of PD-L1 expression, and different methods used to stain and measure PD-L1 expression within tumors. These may account for the challenges in correlating PD-L1 expression with patient response (5–8). Therefore, identifying novel biomarkers to accurately predict the likelihood of response to immunotherapy is crucial in treatment selection and planning for each NSCLC patient.
Multiplexed imaging of tissues, and their analysis, is an emerging and proficient approach aiding early detection, clinical cancer diagnosis, treatment planning, and prognosis (9–13). These images enable the precise determination of spatial distributions of cells and cellular states and the characterization of tumor-immune interactions in situ at the single-cell level. Furthermore, they allow simultaneous detection of various protein biomarkers on the same tissue sample, permitting molecular and immune profiling of NSCLC, while preserving tumor architecture (10,14), and enable the direct quantification of response to a given treatment. However, quantitative image analytic tools must be developed to enable these data to drive predictive biomarkers of response.
The analysis of multiplexed images generally utilizes one of the two state-of-the-art image analysis approaches, i.e., cell segmentation and quadrats. Single-cell data, generated using cell segmentation methods (15), is a fundamental step in many biomedical studies. It quantifies the marker expression of individual cells in a field of view (FoV) and is used to identify cell morphologies, cell subpopulations and neighborhoods, and the stochasticity of marker expression. The approach can be used to extract tumor, stromal, and invasive regions, as well as identify the cells that colocalize in such regions. Dimensionality reduction of the single-cell data can aid with visualization of this high-dimensional data. Single-cell data relies on accurate cell segmentation which can be cumbersome depending on the number of samples and the technical noise of the samples.
The second approach, quadrat analysis, is an important process in image analysis where the high-dimensional image is sequentially tiled, summarized, and classified. Quadrats can be used to identify aggregate marker expression and differential markers across the population of cells in a sample and is typically used to detect interesting regions across high-dimensional images. This does not require any cell-specific information, making it a simple yet effective approach. In a significant amount of high-throughput image analysis research, using tiles for pattern matching and classification is standard (15). As ecologists have long worked with quadrat count data, an additional benefit of tiling the images is that it allows one to leverage many of the spatial analyses originating from ecology, an approach taken in several studies of the tumor microenvironment (TME) (16–19).
Traditional, yet foundational, techniques using multiplexed images often rely on non-spatial analyses which focus on generating various cellular characteristics without taking into account their location or arrangement within the TME. For instance, this involves determining the average biomarker/protein intensity within a set of cells or across an entire image, regardless of the cells’ location, or calculating the cell type proportions in an entire image, without considering their spatial distribution (20).
With recent advances in spatial omics, there has been an increased interest in performing spatial analysis to explore spatial tumor ecologies, i.e., the cellular and quadrat patterns and interactions within the tumor (21–23), how these tumor ecologies are related to patient response to treatment, and how these are impacted by single-agent drugs or immune checkpoint inhibitor (ICI) combination therapies. As tumors are dynamic, rapidly evolving and undergoing ecological transformations, the study of cancer evolution has to incorporate an understanding of the dynamic ecology to explain heterogeneity, behavior, and eco-evolutionary feedbacks (24–26). In this paper, we undertake spatial analyses of images, which preserves and utilizes the spatial relationships between cells and quadrats to understand their interactions with the TME and influence on disease progression (20,27,28).
In this study, our goal is to uncover the evolving spatiotemporal ecologies that drive progression in NSCLC patients and identify predictive biomarkers of response, by leveraging multiplexed images. Here, we analyzed multiplex histological images from pre- and on-treatment (days 15–21) in nine patients with immunotherapy-refractory NSCLC treated with an oral HDAC inhibitor (vorinostat) combined with a PD-1 inhibitor (pembrolizumab) (29). We developed a computational multiplexed-image analysis pipeline using both cell segmentation and quadrats to better understand the spatial and temporal features of multiplexed NSCLC images, and their ability to predict disease progression, as well as identify clinical biomarkers. Performing non-spatial and spatial analyses on the images, and using cell segments and quadrats, has led to novel insights into the TME of patients who can benefit and those who resist immunotherapy, which would have been difficult to obtain using traditional methods alone. Specifically, we defined distinct spatial neighborhoods and associated cellular interactions, and inferred spatial ecologies using only patient-level labels. An important aspect of our study is that combining both cell segment and quadrat approaches allowed us to quantify the unique spatial ecologies in stable disease (SD) and progressive disease (PD) patients, and enabled response predictions based on the abundances and diversity of discovered marker phenotypes and neighborhoods. These approaches, if validated in larger cohorts, may provide improved predictive histological biomarkers from pretreatment samples to drive subsequent treatment decisions. Further, these approaches can be generalized to study multiplexed images from any cancer type.

Materials and Methods

Patient data
We obtained written informed consent from all the patients before enrollment. The trial was conducted in accordance with the provisions of the Declaration of Helsinki, Good Clinical Practice guidelines (as defined by the International Conference on Harmonization), and applicable regulatory requirements. The trial was approved by the central institutional review board (29). Patients were assigned to SD/PD groups based on their response and RECIST criteria. No cell lines were used in this study. All analyses were performed blinded such that experimenters performing data analysis were unaware of the patients.

Multiplexed image processing pipeline for cell segmentation in FoVs
Each 7-stain Vectra multiplexed image is available as a tiff stack with 7 sub tiffs: one tiff for each marker. We have built an image processing pipeline that consists of three parts to analyze each tiff stack (Supplementary Figure 1i):
the pre-processing unit that denoises and cleans each 1008-pixel × 1344-pixel sub tiff of the 7-stain Vectra Tiff stack (Supplementary Figure 1ia). We start with generating the grayscale image of each sub tiff and choose the nuclear image stained with DAPI for cell segmentation (Supplementary Figure 1ib). Since the images are ridden with technical noise, it is crucial to extract true stain signal from noise. For this we apply an Otsu thresholding to the images which involves iterating through all the possible threshold values and calculating a measure of spread for the pixel levels each side of the threshold, i.e., the pixels that either fall in foreground or background. Foreground pixels are regarded as true stain signals while background pixels are noise. Next, the images are tiled into 256-pixel × 256-pixel frames (Supplementary Figure 1ic).

The cell segmentation unit: The image tiles are used as input to U-Net (30), a Fully Convolutional Network, for cell segmentation (Supplementary Figure 1id). We further trained U-Net using nuclear segments available from Kaggle’s Data Science Bowl.

Downstream analysis unit: Once the cell segmentation predictions are made per tile (Supplementary Figure 1ie), we stitch the tiles (Supplementary Figure 1if) and perform clustering on the stitched image (Supplementary Figure 1ig). For clustering, we use a Gaussian Mixture model (31,32) which identifies heterogeneous cell types based on presence of stain expression, and does not require the number of clusters to be known in advance.

In the FoVs, cell segments consisted of cells from tumor and stromal regions.

PD patients have high PD-1 and T regulatory cells (FoxP3+) in both the tumor border and tumor regions
Using our multiplexed imaging pipeline for single cells (Supplementary Figure 1i), we extract the cell segments per FoV using a U-Net (30), build a count matrix with cells as rows and markers as columns, and cluster the count matrix to identify heterogeneous cell types. From these cell types, we automatically demarcate tumor-rich regions, across images, that are higher in PanCK expression (red regions in Supplementary Figure 1ii, right panel). The tumor regions are approximated using multiple convex hulls to allow for spatial quantification of tumor and immune cells within and outside the tumor border. We note that for both pre and on treated patients, there is a higher presence of immune cells within the tumor (Supplementary Figure 1iii) and a higher abundance of tumor and immune cells across the tumor border (Supplementary Figure 1iv) in PD than SD patients signaling that there are distinct cellular compositions for PD and SD patients.

Convex-hull approximation
To demarcate tumor regions within an FoV, we draw convex hulls around those clusters (obtained from the above cell segmentation and clustering process) that are high in PanCK. For this we use the scipy.spatial.ConvexHull (33) function that finds the smallest convex set that encloses all the region coordinates, forming a convex polygon. The coordinates for plotting the polygon are obtained from the regionprops (34) function generated during cell segmentation. Once the initial convex hulls are identified, the process is repeated to plot hulls extending both inside and outside the initial hulls using fixed boundary offsets. For example, a fixed offset id is ‘added’ to the coordinates of the initial hulls to plot hulls extending ‘outside’ the initial hulls (Supplementary Figure 1ii).

Quantifying the tumor-immune cell colocalization at the tumor border
Next, to quantify the tumor-immune cell colocalization at the tumor border we cluster the tumor-immune cell counts at the tumor border approximated using convex hulls (see Supplementary Figure 1ii, right panel). Clustering is performed using a Gaussian mixture model (35) (GMM). The input matrix to the GMM consists of cells (as rows) and their marker distribution as features (columns). The rows of the input matrix are ordered based on the cluster assignments and the z-score of the marker expression (columns) is averaged over vectors per cluster representing a cell type (rows). Those markers that have a higher z-score per cluster are identified as the differentially expressed markers for that cluster. The clusters are visualized using a 2D t-SNE plot where each point abstracts an image (Supplementary Figure 2A, B). We obtain distinct clusters showing higher colocalization of PanCK+PD-1+FoxP3+ for PD patients and higher colocalization of PanCK+PD-L1+ along with CD3+CD8+ immune cells for SD patients indicating underlying structural differences between PD and SD patient groups. Further, to better understand the relationship of the clusters to the actual images, we generate an image t-SNE (Supplementary Figure 2C) where we replaced each dot in the t-SNE of Supplementary Figure 2A, B with the corresponding multiplexed images. This image t-SNE was rendered using Mistic (36) and this arrangement of images as an image t-SNE clearly highlights the difference in immune cell abundance across PD and SD patient groups. The image t-SNE also helps in identifying groups of immunologically hot, warm or cold tumors across patients based on immune cell colocalization at the tumor margins and tumor regions.

Multiplexed image processing pipeline for quadrat counts
The quadrat-based analyses are based on marker abundances in the entire quadrat (i.e., those possibly outside nuclear masks), and so used a slightly different approach to remove channel bleed and noise.
Markers were assigned into mutually exclusive groups, i.e., those expressed only on T-cells (CD3, CD8, FoxP3, PD-1), or tumor cells (PanCK, PD-L1). Next, each pixel in the image is assigned to one of those groups based on which marker had the maximum expression. For example, if a pixel has fluorescence for PanCK and FoxP3, it would be assigned to the T-cell group if FoxP3 (or any other T-cell marker) is greater than the expression of PanCK (or any other tumor marker), or tumor group if PanCK is higher than FoxP3. Once assigned to a group, the pixel values in the other group’s channel(s) are set to zero. The result is that a pixel can be positive for more than 1 marker in the same group (i.e., T-cell that is CD3+/CD8+), but not for markers in mutually exclusive groups (i.e. not possible to have CD3+/PanCK+).

Channel bleed is further reduced by subtracting the intensity of the adjacent channel, but only if that other channel belongs to a mutually exclusive group and uses an adjacent fluorophore. For example, if a pixel is positive for FoxP3 (Opal 650) and PanCK (Opal 690), and FoxP3 expression is greater than PanCK, then the pixel is classified as T-cell and the PanCK channel is subtracted from the FoxP3 channel.

After removing channel bleed, each channel is thresholded to determine the number of pixels that are positive for each marker. In the case of the FOV, Otsu thresholding was used to separate foreground (pixels positive for each marker) from background).

Non-spatial analyses
To compare the underlying patterns of tumor evolution in NSCLC between PD and SD patients, we examined the density of markers across the FoVs (FoVdata dataset) using basic ecological analysis, for example, average quadrat counts. We view each FoV as a collection of fixed-size patches, called quadrats, defined as partitions of an image into multicellular regions, creating sequential tiles in order. Each quadrat is 100μm × 100μm in size and expresses the area normalized sum of marker positivity (after thresholding) of all cells present within that quadrat. In Figure 1A, an example 7-marker FoV (with tumor regions in cyan) is shown (i) along with its corresponding cell segments (ii) and quadrats (iii). We plot the marker abundances before and during treatment using: a) average positive pixel counts (PPC) across all the quadrats of PD and SD patients (Figure 1B); and b) using average cell counts from PD and SD patients (Supplementary Figure 3). We further investigate the difference in PD and SD marker densities by looking at the non-metric multidimensional scaling (NMDS) (37) ordination of the positive pixel counts (PPCs) for each FoV. NMDS shows the tumor composition of each treatment category. This is then followed by a confirmatory permutational analysis of variance (PERMANOVA (38)) test on the ecological dissimilarity matrices (computed using PPCs of the quadrats) (Figure 1C). Additionally, using the PPCs, we plot the pairwise marker correlations across SD and PD patients (Figure 1D).

Spatial analyses with cell segments
Using our multiplexed imaging pipeline for single cells (Supplementary Figure 1i), we extract cell segments and tumor regions per FoV. Once the cell segments are extracted using a U-Net (30), we area-normalize the cells per FoV, build a count matrix with cells as rows and markers as columns, and cluster the count matrix to identify heterogeneous cell types using a Gaussian mixture model (GMM) (39). To understand the spatial organization of tumor and immune cells across PD and SD patients, we investigate each cell in the context of its spatial neighbors. This is done by generating cellular neighborhoods (CNs) where the marker expression of each cell is the average of ten of its nearest spatial neighbors in Euclidean space (Figure 2A-C). A CN can be defined as the minimal set of cell types that are both functionally and spatially similar.

Construction of CNs:
Each cell is represented as a binary vector of size 7 representing the 7 markers (1 = marker present, 0 = marker absent). To calculate the CN for each cell, we average ten (or less) of each cell’s nearest spatial neighbors in Euclidean space. Through this construction, for each cell, we have a vector of size 7 that quantifies its spatial neighborhood. In addition to its co-expressed markers, one cell now jointly describes which of the 7 markers make up its spatial neighborhood. This approach is inspired by a previous CN implementation (40) with the major difference that we do not have a fixed set of cell types at the onset but work with cell types that are organically inferred through the clustering process. Along with the cell segments, we also computationally demarcate tumor-rich regions (approximated using multiple convex hulls; Methods), across images, that are higher in PanCK expression (red regions in Supplementary Figure 1ii, right panel). To compute the colocalization of markers, we use a modified Morisita-Horn co-localization index (16,41).

Weighted Morisita-Horn measure
To quantify the co-localization of tumor and immune cells regions, we use the Morisita-Horn (MH) similarity measure (16,41). For two populations, immune and tumor, with proportions and respectively, the MH measure can be written as:.
We consider a weighted form of the MH measure where weights across cells of the two populations are calculated using inverse-distance weighting. This weighting influences the similarity measure based on the cells that are closer to, than farther away from, the cell in question. The weighted MH measure is given as:

where is the average distance of a tumor cell, , to its neighboring immune cell, . is the total number of immune cells and is the total number of tumor cells.

Spatial analyses with quadrats
In addition to the cell segments, we process image quadrats (refer Methods) by identifying quadrat neighborhoods (QNs) (Figure 2D, E) and building species association networks (42,43) (SANs) and Gaussian copula graphical lasso (GCGL) using the ecocopula package (44) applied to pre- and during-treatment images (Figure 2F).

Description of the GCGL framework:
The GCGL model combines the Gaussian copula framework with the graphical lasso method. It is used to estimate graphical models for multivariate data, particularly when the variables are not necessarily normally distributed but can be transformed to normality using a Gaussian copula. This approach is useful in various fields, including finance, ecology, and genomics, for modeling dependencies between variables.

Gaussian Copula:
The Gaussian copula is a statistical tool that models the covariance structure between multiple random variables, independent of their individual marginal distributions. It works by transforming the variables to a multivariate normal distribution and then estimating the covariance structure using a correlation matrix.

Graphical Lasso:
The graphical lasso is a method for estimating a precision matrix (also known as the inverse covariance matrix) of a Gaussian vector. It does this by applying an L1 penalty (lasso penalty) to the precision matrix, which promotes sparsity and identifies conditional dependencies between variables.
The GCGL combines the strengths of both methods. It uses the Gaussian copula to model the covariance structure between variables (even for non-normally distributed random variables), and then uses the graphical lasso to estimate the precision matrix, which identifies the conditional independence relationships between variables.
In our setting, the variables GCGL operates on are the markers.
Similar to CNs, QNs are defined as the minimal set of functionally and spatially similar quadrats. These are extracted using the Leiden clustering algorithm so as to identify spatial communities.

Inferring underlying networks of SD and PD patients using prior information
To quantify the distinct architectures across PD and SD patients, we build a statistical model to infer the difference in marker interactions, in the form of a network, between these patient categories. This is done by: utilizing the species association networks (42,43) (SANs; Figure 2F) along with the edge weights, and marker expression both at the pre-treatment and on-treatment states for each category; and capturing the ‘difference network’ that potentially drives the patients from the pre-treatment state to the on-treatment state. The inference problem can be stated as a least-squares minimization problem (35) and can be further regularized using prior biological knowledge of marker interactions.

Problem formulation
We denote the change in the ith marker expression in the pre-treatment case at time , , as:

where is the number of markers, , , and is the weight matrix. Equation 1 can be written in matrix form as:

where and * denotes pre or on states. We can solve Equation 2 as a least-squares minimization problem given markers and as steady states. This can be written as:

Further, we can add prior knowledge of biological mechanisms between markers as additional constraints to the least-squares problem (Supplementary Figure 4A). These mechanisms can be captured in pairwise format as a matrix where a 1 indicates positive interaction, 0 indicates no interaction and −1 indicates negative interaction. Equation 3 can be written as follows to yield :

Implementation
We construct weighted count matrices for the four categories: SD (pre, on) and PD (pre, on), by combining the SANs and the edge weights obtained in Figure 2F. Starting with the SAN for each category, a weighted adjacency matrix, , is created. is a square distance matrix representing connections between markers, wherein ‘1’ indicates a marker-marker connection, and ‘0’ otherwise, and each connection is multiplied by weights from the SAN where the weights are bound by −1 (red/negative association) and 1 (blue/positive association). For each , we construct a similarity matrix, , using with where is the Kronecker delta function and is the number of markers, based on the centering operation in kernel PCA37. The weighted count matrix, can be computed by the eigenvalue decomposition37 of . Assuming that each marker follows a Gaussian distribution, we can model these weighted count matrices of cells by markers to follow a multivariate Gaussian distribution, without loss of generality. We compute the first and second order moments to describe each of the four categories. Next, for each category, we randomly draw 10,000 samples creating a random matrix of 10000 rows and 5 markers. This is repeated 100 times and the random matrices are averaged out to create one representative matrix denoting that category. This is done to ensure that we work with random matrices that a) are invariant to biases and shifts while generating them from distance and similarity matrices and b) have the same number of rows in the ‘pre’ and ‘on’ states that obey the observed/empirical distributional moments. Implementing Equation 4 by incorporating the biological prior, gives us the optimized for SD and PD (Supplementary Figure 4B i,ii) respectively. It becomes evident through for SD (Supplementary Figure 4Bi) that there is an increased CD8 interaction with the tumor, a decrease in FoxP3 interaction with tumor and a positive interaction between CD8, PD-L1 and PanCK. For PD (Supplementary Figure 4Bii), the indicates an increase in FoxP3 interaction with tumor, a decrease in CD8 interaction with tumor and a positive interaction between FoxP3, PD-1 and PanCK. Note that these interactions that differentiate SD from PD were previously identified from the image t-SNE analysis (see Supplementary Figure 2) as well as the quadrat analysis using (Figure 2D-F).

Assessing the reliability of CNs distinguishing SD and PD patients
To show the importance/reliability of the characteristic cellular neighborhoods (CNs) defining PD and SD patients, we set up a supervised learning framework to compare the reliability of the SD CNs (PanCK+PD-L1+ with CD3+CD8+) with that of the PD CNs (PanCK+FoxP3+PD-1+) for each patient. Specifically:
We consider two subgroups of data per patient: the SD CNs group (PanCK+PD-L1+ with CD3+CD8+) and the PD CNs group (PanCK+FoxP3+PD-1+)

Calculate the geometric mean per group for each patient. The geometric mean is an average that is calculated by multiplying the spatial contribution of all the markers in each group and then taking the nth root, where ‘n’ is the total number of entries in that group. This geometric mean approach is particularly useful when dealing with data that involves multiple averages (like the CN implementation described in Materials and Methods section), and percentages (spatial contribution for each marker is normalized based on area).

Compute the ratio of geometric means (RoGM) for each patient:
RoGM = geometric_mean(SD_CN_group) / geometric_mean(PD_CN_group)

The ratio indicates how much larger or smaller the geometric mean is in the SD marker group compared with the PD marker group.

Goal is to show that the RoGM can separate the two response groups clearly.
Higher RoGM indicates higher presence of the SD CN group (PanCK+PD-L1+ with CD3+CD8+)

Lower RoGM indicates higher presence of the PD CN group ((PanCK+FoxP3+PD-1+)

Figure 3A plots the RoGM per patient. By using the SD and PD markers in this supervised manner, we can show that these biomarker combinations can indeed separate the two patient groups: Note the clear separation of the green dots (RoGM for SD patients) and orange dots (RoGM for PD patients).
We also generated the Kaplan Meier (KM) survival curves using Python’s lifelines package (RRID:SCR_024899) for the above RoGM analysis. In Figure 3B, we show the KM analysis using the response groups (SD and PD (i)), and the KM analysis using RoGM (ii). For Panel B, we binarize the RoGM vector at a decision boundary = 0.5 i.e., we set RoGM values > 0.5 = 0, and RoGM values < 0.5 = 1. Comparing both Panels, we see that our RoGM metric (ratio of the SD CN markers to the PD CN markers) successfully differentiates the SD and PD groups without overlap.

RoGM applied to a non-spatial setting:
We repeat the RoGM and variations of the RoGM in a non-spatial setting i.e. the input to the RoGM analysis are the markers from FoVs without the position information (Supplementary Figure 5). We look at the ratios of the SD and PD CN groups for each patient in the non-spatial setting (Supplementary Figure 5B). Additionally, we look at the ratio of individual markers for each patient (Supplementary Figure 5C-H). Note that in these panels B-H, there is no clear decision boundary (or separator line) that can be drawn to separate the SD from the PD patients. It is only when we analyze the markers collectively in their spatial setting that we obtain a clear separation between SD and PD patients (Figure 3A). With these experiments, along with our clustering and prediction experiments, we can authenticate the importance of spatially analyzing the images and further validate the presence of specific SD and PD spatial neighborhoods.

Support Vector Machine (SVM) used to make predictions with single cell data
We implemented the SVM classifier using Python’s (RRID:SCR_008394) scikit learn’s (45) (sklearn; (RRID:SCR_002577)) SVM, SVC classes. The classifier uses a linear kernel. Classification statistics (accuracy, precision, recall, and F1-score) were calculated using scikit learn’s metric package. The LOOCV was implemented using the LeaveOneOut class provided in sklearn’s model_selection suite.

Neural Network Species Distribution Model (NN SDM) for predicting disease progression using quadrats
We use a feed forward neural network (NN) to identify species (marker) distribution interactions (SDM) to predict disease progression in patients (Supplementary Figure 6). The inputs to the NN SDM are quadrats generated from the FoVdata images. Quadrat dimensions for the NN SDM are 50μm × 50μm. The activation function used is Leaky Rectified Linear Unit (46) (ReLU). The loss function used is cross balanced loss (47), which was minimized using the AdamW (48). Due to the unbalanced nature of the data, the Synthetic Minority Over-sampling Technique (SMOTE) (49) was used to oversample the minority class, thereby balancing the training data. To analyze feature importance (Figure 4D), we use integrated gradients (50) available through PyTorch’s model interpretability library, Captum (51).

Training classifiers for predicting disease progression
Using FoVdata, we trained and tested classifiers based on support vector machines (SVMs), and neural network species distribution networks (NN SDMs) to predict disease progression of patients. The response of the patients was known a priori. Using the Leave-one-out cross-validation (LOOCV) setup, at the patient level, the classifiers were trained on the training set and predictions of disease progression were made on the test set. To generate these predictions, we used images from 8 patients (5 SD and 3 PD) from the 9-patient cohort since one patient (classified as PD) had only a single image taken before treatment.
We used both SVMs and NN SDMs as these are popular and well established supervised learning techniques to build prediction models (52). While SVMs are powerful for finding optimal decision boundaries, especially in high-dimensional spaces, and with non-linear separation boundaries, NN SDMs are considered non-parametric, making no assumption on the distribution of data and the structure of the true data generator. The similarity of biomarkers obtained by both the SVMs and NN SDMs show that the results obtained are not artefacts of the specific prediction algorithm used but are prominent underlying patterns in the data.

Approaches to combine cell segments and quadrats
We describe three approaches to combine the insights derived from cell segments and quadrats.
In the first crossover experiment, we created a weighted cell marker matrix where each cell’s marker expression is multiplied by the average of its corresponding quadrat progression probabilities (as obtained from the NN SDM). This weighted marker matrix is then used to train the SVM. The motivation to use such a weighted approach was to primarily handle unbalanced sample sizes (5 SD and 3 PD patients in FoVdata), as well as unbalanced regions of disease progression within each FoV. We observe that this approach leads to an improvement in the overall classification accuracy at the cell/quadrat level per patient, especially for PD patients (Supplementary Table 1), and that these results can be further augmented by additional patient samples.
For the second approach, we describe a crossover experiment between the cell segmentation and quadrat approaches to further reinforce the distinct architectures existing between PD and SD patients (Figure 5A,B). For this, we use images from only the pre-treatment cases and redefine the quadrat into both finer and coarser scales, and cluster the quadrats. In the finer scale, each quadrat is represented as a cell with its neighborhoods (Figure 5Ai, ii), and in the coarser scale, each quadrat is a fixed-size patch with a collection of cells and their neighborhoods (Figure 5B).
In our third approach to quantify those markers that potentially drive the patients from the pre-treatment state to the on-treatment state, we build a statistical model to infer the difference in marker interactions, in the form of a network, between these patient categories. This is done by utilizing the SANs (Figure 2F) from the FoVs, across response groups and across pre- and during-treatment states. The inference problem can be stated as a least-squares minimization problem37 and can be further regularized using prior biological knowledge of marker interactions (Supplementary Figure 4).

Results

Study design and data source
The patient cohort used in this study was part of a randomized Phase I trial (NCT02638090) (29). They were treated with an oral HDAC inhibitor (vorinostat) combined with a PD-1 inhibitor (pembrolizumab). Eligible patients had histologically confirmed advanced or metastatic (Stage 4) NSCLC with progression, and they were also refractory or relapsed on immunotherapy. Core biopsies were collected from all patients both pre- and on-treatment (days 15–21). All pre-treatment biopsies were obtained immediately before treatment and are not archival tissues. Patient responses were classified into two groups: those who remained in stable disease (SD) after 24 weeks; and those who had progressive disease (PD) before 24 weeks (29).
We used 92 PerkinElmer Vectra Field of Views (FoVs) (hereby referred to as FoVdata) from nine patients from the above cohort. Of the nine patients, five qualified as SD and four as PD. There are 58 images from the SD patients and 34 images from the PD patients. Using the seven-color multiplex fluorescence assay by PerkinElmer-Opal-Kit, each image is stained for seven markers: DAPI (nuclei); CD3 (T cells); CD8 (effector T cells); FoxP3 (regulatory T cells); PD-1 (inhibitory receptor on immune cells); PD-L1 (immune checkpoint marker); and Pan-cytokeratin (PanCK, epithelium). We had a median of 10 images per patient in FoVdata, split evenly across pre- and on-treatment, except for one patient for whom we had access to only a single image taken pre-treatment. We analyzed the FoVdata to identify spatial ecologies and to make disease progression predictions (Figure 6). Both stromal and tumor regions were imaged in the FoVdata. We summarize FoVdata in Table 1. An exemplar FoV with its corresponding cell segments and quadrats is shown in Figure 1A, along with its individual markers and composite marker image (Supplementary Figure 7). The various metastatic sites from where the FoVs were chosen are provided in Supplementary Table 2.

Non-spatial analyses indicate PD and SD patients have significantly different ecologies
Through our non-spatial analysis using positive pixel counts (PPCs) derived from FoVdata quadrats, it is evident that PD tumors demonstrate increased marker levels throughout (pre- and during-treatment) (Figure 1B). We overlay the NMDS (non-metric multidimensional scaling) projection of the average PPCs across all FoV quadrats with the PERMANOVA test (see Methods) to depict that the centroids, indicating the spread of each group, are different; this reinforces that PD and SD have significantly different pre-treatment ecologies (p-value = 0.007). Since each groups’ centroid shifted while on treatment, this also indicates that treatment itself alters the tumor ecology, as shown by the arrows in Figure 1C. In Figure 1D, we plot the pairwise marker correlations across SD and PD patients, computed based on PPCs, before and during treatment, and again note that marker correlations vary across response groups. For instance, in PD patients (both pre- and during-treatment), the correlations between FoxP3 with PanCK, and PD-1 with PanCK, are noticeably higher (blue) than in SD patients (red). These significant differences in composition derived from basic non-spatial analyses, where neither the pixel’s positional information nor their spatial distribution was considered, suggest that the pre-treatment ecology may be sufficiently different in PD and SD to be predictive of response.
Our finding that a higher presence of FoxP3+PD-1+ with PanCK cells (Figure 1B, D) and a reduced presence of CD8+CD3+ with FoxP3+ cells distinguishes PD patients is consistent with the findings of other multiplexed immunofluorescence (mIF) studies showing that immune suppression via FoxP3 Tregs is a sign of poor prognosis (53–55).
To further gauge the prognostic significance of tumor-immune heterogeneity, and to gain a detailed understanding of the presence and evolution of this spatial heterogeneity in the pre- and during-treatment phases, we spatially analyzed these images by partitioning pixels into cells or quadrats. Through this approach, the cell/quadrat locations are preserved to identify specific immune neighborhoods and their roles in enabling disease prediction.

Spatial analyses highlight distinct spatial communities in SD and PD patients
We further study the spatial heterogeneity of tumors across response groups, by investigating the differences pertaining to the spatial organization of the tumor with respect to its microenvironment in the FoVdata cohort using cell segmentation and quadrat approaches. By clustering cells from all FoVs along with their spatial neighbors, we obtain twelve distinct cellular neighborhoods (CNs) (Figure 2A, Supplementary Figure 8A, Methods). An example of the spatial orientation of the CNs mapped onto a FoV is shown in Figure 2B. The t-SNE rendering of these CNs demonstrates that there are characteristic cellular neighborhoods defining PD and SD patients (Figure 2Ci, ii) such as PanCK+FoxP3+PD-1+ for PD patients (clusters 6 and 11) and PanCK+PD-L1+ with CD3+CD8+ for SD patients (clusters 4, 5, 7 and 9). It is evident that prior to treatment, PD patients show a higher number of cells expressing FoxP3+PD-1+ and PanCK+ (Supplementary Figure 8Ai). On the other hand, there is always an increased number of cells expressing CD3+, CD8+, PD-L1+ and PanCK+, both before and during treatment in SD patients. Further, we set up a supervised learning framework to compare the reliability of the SD CNs (PanCK+PD-L1+ with CD3+CD8+) with that of the PD CNs (PanCK+FoxP3+PD-1+) for each patient using a ratio of geometric means framework, and successfully validated using a Kaplan Meier survival curve analysis (Figure 3; Supplementary Figure 5; Methods). These findings are consistent with prior research, where FoxP3+ cells were known to foster an immunosuppressive niche that promoted tumor growth. One study that analyzed images from 261 NSCLC patients revealed that strong interactions between FoxP3+ and neighboring PanCK cells in the tumor border and tumor centers were associated with adverse patient prognosis (56). While increased FoxP3 cells in tumor regions was shown as an independent predictor of worse overall survival in NSCLC, increased CD8+ among PanCK cells was significantly associated with better survival (55,57).
These interactions have been further reinforced visually using an image t-SNE generated using Mistic38 (Supplementary Figure 2). Further, the tumor boundary analyses using convex hulls indicates that: a) for both pre- and on-treatment patients, there is a higher presence of immune cells (primarily, FoxP3) within the tumor of PD patients, and more so within PD patients during treatment (Supplementary Figure 1iii); b) a higher abundance of tumor and immune cells at the tumor boundary (Supplementary Figure 1iv) in PD than SD patients; and c) the potential to identify groups of immunologically hot, warm or cold tumors across patients based on immune cell colocalization at the tumor margins and tumor regions (Supplementary Figure 2). Together, these results conform with previous studies (55–58) and support that there are distinct cellular compositions, both within the tumor regions and at the tumor boundaries, for PD and SD patients.
Next, we cluster the quadrats generated from FoVdata using a Gaussian mixture model (GMM) (39) and illustrate the inferred quadrat neighborhood (QN) clusters on t-SNE renderings (Figure 2D) and the marker contributions per cluster (Figure 2E). It is again evident that the QNs capture the distinct ecologies defining PD and SD patients such as PanCK+FoxP3+PD-1 (clusters 6 and 11) for PD patients, and PanCK+PD-L1 with CD3+CD8 (clusters 4 and 5) for SD patients.
To further explore the differences spatial structure between PD and SD patients, we next construct species association networks (42,43,59) (SANs) (Figure 2F) using the Gaussian copula graphical lasso (GCGL) (44) method to identify marker interactions. SANs provide a modelling framework to extract direct spatial associations based on spatial co-localization, while removing indirect effects. This is a popular approach in quantitative ecology used by landscape ecologists (60,61). The fundamental modeling steps in SANs are data preparation, model fitting and assessment, and prediction. The quadrats are combined and fit with a Gaussian copula graphical lasso (GCGL) using the ecocopula package (44) applied to pre- and during-treatment images (refer Methods for additional details). Gaussian copula can directly model covariance and strength between markers, revealing direct associations while removing indirect associations. We leverage this aspect of the GCGL to create a single association network with all direct links between markers and the strength of association increases from red (negative/dispersed) to blue (positive/clustered). Spatial communities (i.e. groups of cells that cluster together in space) are shown by the larger polygons capturing sets of markers that are clustered using the Leiden algorithm (34). By collectively examining each marker’s importance, response, and interactions across quadrats within the community associations, we again observe distinct cellular communities distinguishing PD and SD patients: the pre-treatment association of FoxP3 with PD-1 and PanCK in PD patients; and the on-treatment association of CD8 with PD-L1 and PanCK in SD patients.
Additionally, we inspect the response heterogeneity across CNs and QNs (Supplementary Figure 8). We note that in general, FoxP3+PD-1+PanCK+ is higher in PD CNs and QNs (refer orange outlined neighborhoods) while CD3+CD8+PD-L1+PanCK+ (refer green outlined neighborhoods) is higher in SD CNs and QNs. Through these clustering analyses, we demonstrate that PD tumors are characterized by an immune-suppressive environment (due to presence of FoxP3+ and PD-1+ cells) prior to treatment. PD-1+ immune cells may have other checkpoint receptors such as LAG3, TIM3 or TIGIT that would indicate a deeper immune cell exhaustion state (62) further reducing the response to treatment. On the other hand, SD tumors show the presence of higher immune cell infiltration indicating durable clinical benefit (56,57).
In a similar NSCLC study that involved 8 patient FoVs obtained pre-treatment, it was shown that responders (complete responders (CR), partial responders (PR)) associated with a higher spatial proximity of CD8+FoxP3+PD-1+ to tumor cells, whereas nonresponders (SD and PD) exhibited a low spatial proximity of CD8+FoxP3+PD-1+ (58). Since our study analyzed SD and PD separately, the result that CD8+FoxP3+PD-1+ is a potential predictive biomarker (58) conforms with our findings when considering the marker combinations for SD and PD taken together. A study on 57 NSCLC cases showed that a higher spatial proximity of CD8+ with PanCK correlated with better patient survival (11). Yet another study that spatially analyzed 120 NSCLC images from patients with stage I-III disease found that a greater degree of interaction between tumor cells and FoxP3 was associated with worse survival outcomes (57). These studies cumulatively indicate that the proximity of FoxP3 to tumor cells may hinder the immune system’s ability to fight cancer. These results resonate with our observation that a higher presence of FoxP3+ PD-1+ with PanCK was prevalent in PD patients while CD8+ with PanCK had a higher presence in SD patients.

Fundamentally distinct ecologies across PD and SD patients enable disease progression prediction and clinical biomarker identification
For the FoVdata cohort, we trained a support vector machine (SVM) using the single cell data to predict response to therapy at the patient level. Using the leave one out cross validation (LOOCV) approach, the SVM achieves an overall mean accuracy of classification of 87.5%, with class-specific mean accuracy of 100% (SD) and 66.7% (PD) (Supplementary Table 3
row1). To identify the marker importance per patient and therefore for each response group, we assess the weights obtained from the SVM. These weights represent the vector coordinates which are orthogonal to the SVM hyperplane, and their direction indicates the predicted class. The absolute size of the coefficients in relation to each other can then be used to determine marker importance for the prediction task. The marker importance per response group, generated by the SVM, indicates that CD3+CD8+PD-L1 are important markers for SD patients and FoxP3+PD-1 for PD patients, respectively (Figure 4A). These predictive markers for SD and PD are in agreement with the markers identified using spatial analyses of cells and quadrats in the previous section, as well as studies that show that FoxP3 infiltration into tumors is a contributor to poor prognosis due to its immunosuppressive role, while CD3, CD8 infiltration tends to mitigate this effect and is associated with better survival (55–58).
A neural network (NN) species distribution model (SDM) (see Methods, Supplementary Table 1) trained on quadrat counts has an overall mean accuracy of 75% with class-specific mean accuracy of 60% (SD) and 100% (PD) (Supplementary Table 3
row2). The distribution of quadrat probabilities of PD classification, generated through the NN SDM (Figure 4B), could be used to assess risk of progression, as on average patients with stable disease (SD) had lower quadrat probabilities than PD patients. This distribution for each patient tested under the LOOCV method is compared to the distributions of both SD and PD patients generated during training, using the Bhattacharyya distance (BD) (63,64) as the distance/similarity metric. The test distribution is assigned to either SD or PD based on the closest BD distance, i.e. the more similar distribution of progression probabilities (Figure 4C). A feature-importance analysis (50,51) reveals that PD-1, FoxP3, CD8 and PanCK had the greatest impact on the NN SDM’s ability to accurately predict treatment response, suggesting these markers play the most important roles in determining whether a patient responds to treatment (Figure 4D). Positive values indicate the marker is important in predicting PD, while negative values indicate importance in predicting SD. Further, the NN SDM predictions allow us to create a risk map that identifies areas of the tumor associated with different probabilities of overall patient progression while on treatment (Figure 4B).
In addition to the above, we perform three crossover experiments involving cell segments and quadrats. In the first crossover experiment, we implement a weighted SVM. For this, each cell’s marker expression is weighted by the average of the corresponding quadrat progression probabilities (as obtained from the NN SDM), and this updated matrix is used to train the SVM. With this approach, while the predictions at the patient level remain the same (Supplementary Table 3
rows1 and 3), the predictions at the number of cells/quadrats correctly classified show an improvement in the overall classification accuracy (up to 88.5%), especially for PD patients which improved from 66.7% to 69.39%, while also increasing the overall precision (Supplementary Table 1
row3). For the second crossover experiment, we cluster the cells across all quadrats from the pre-treatment cases (finer scale setting; see Methods), project the cells onto a 2D t-SNE, and depict the NN SDM disease predictions (generated using the coarse scale setting) to obtain a modest overlap and clearer distinction between PD (orange) and SD (green) regions (Figure 5Ai, ii). In the third crossover experiment, we use the coarser scale setting, where quadrats are sequential tiles of the images, and utilize the NN SDM predictions to color the quadrats for each FoV generating a risk map (Figure 5B). These risk maps are arranged in a 2D layout using Mistic (36) where the 2D co-ordinates are generated by clustering the cellular marker densities for each FoV. Their orientation suggests that spatial measurements of the tumor ecology, and thus the interactions between markers, can be used to make accurate predictions of how a patient will respond to therapy. We observe that most of the PD images are colored in yellow/orange/red, correctly identifying these images, and therefore the patients, of having a higher chance to progress. Similarly, SD images have green/blue hues denoting a lower chance to progress. Additionally, we verified these differences using a statistical network model (Methods) that confirmed the presence of ecological diversity in pre-treatment biopsies from PD and SD patients.
The predictions from both the SVM and NN SDM are suggestive that distinct ecologies exist between PD and SD patients. This is consistent with prior studies that found infiltration of intra-tumoral CD3+ and CD8+ T cells in NSCLC is associated with better survival outcome (11,56,57,65) and that infiltration of FoxP3 cells is associated with a negative prognosis (55–58) due to FoxP3’s role of suppressing anti-tumor immunity and promoting tumor progression. Further, in NSCLC patients, spatial proximity of tumor and FoxP3 cells has been found to be associated with diminished survival59. Additionally, since both the SVM and NN SDM leveraged the same predictive ecologies, it implies that these ecologies are inherently present in the data and are not artifacts of the prediction algorithms themselves.

Potential of distinct ecologies as companion biomarkers for PD-L1 status
Through two tests, we further gauge the disease progression predictions using PD-L1 status alone, which is a recommended biomarker for NSCLC, particularly when deciding whether a patient is likely to benefit from immunotherapy treatments. In a first test, we train a SVM on the total number of cells (from cell segmentation) and a NN SDM on the total number of PPCs (from quadrats), using only the PD-L1 expression, from the pre-treatment images in FoVdata to predict patient response categories, i.e., SD versus PD (Supplementary Table 3
rows 4–6). We show that using PD-L1 alone for response prediction has a predictive power of only 63% across both single cells and quadrats for FoVdata (Supplementary Table 3).
In our second test, we utilize the clinically assigned PD-L1 status for each patient. The PD-L1 status is defined as the percentage of tumor cells that express PD-L1, and it is shown that tumors expressing high amounts of PD-L1 (50% or more) may respond particularly well to immunotherapy75. We re-render the FoVdata t-SNE in Figure 5Ai where cells per patient are colored according to the patient’s PD-L1 status (Figure 5C). There are four clinically-assigned PD-L1 status categories for our nine-patient cohort: 0: No PD-L1; 1–49%: Low PD-L1; >50%: High PD-L1; Unknown: unable to assess PD-L1. Additionally, each patient is clinically categorized into either SD (green) or PD (orange) (Figure 5A). We observe that within each PD-L1 status category, there is a mix of both response categories. For instance, the No PD-L1 category has one PD and one SD patient, while the high PD-L1 category has four SD and one PD patient. Further, there was one patient in the unknown category that was assigned to be a PD patient. Note that for this patient, we only had one pre-treatment image and therefore the patient’s data was not used in the prediction experiments and is not depicted in Figure 5C. This shows that there were 3 misclassified patients out of 9, indicating that using PD-L1 status alone, the best patient response prediction we would achieve is 66.7%. In comparison, our ML predictions (Figure 5Aii) using PD-L1 in combination with other markers along with their spatial context, can predict patients with an accuracy between 75%−88%. This immediately highlights the clinical challenge of using PD-L1 status as a lone diagnostic biomarker as well as strengthens the utility of our marker combinations as a companion biomarker for PD-L1 status to predict treatment response.

Discussion
Analyzing multiplexed images of patient tissues has the potential to become a routine clinical procedure for early detection, clinical cancer diagnosis, treatment planning and prognosis for many cancer types including NSCLC. We analyze multiplexed images to quantify the spatial ecology of tumor and immune cell types to understand if these ecologies can predict disease progression dynamics under therapy. Through two cornerstone image processing techniques (cell segmentation and image tiling to generate quadrat counts) combined with spatial statistics, machine learning, and deep learning, we perform in-depth analyses of multiplexed images based on the frequency, phenotype, and spatial distribution of immune and tumor cells within the immune landscape of NSCLC. Multiplexed field of views (FoVs) were obtained from pre- and during-treatment (days 15–21) core biopsies in nine patients with advanced/immunotherapy-refractory NSCLC treated with an oral HDAC inhibitor (vorinostat) combined with a PD-1 inhibitor (pembrolizumab). We characterize cellular composition and spatial neighborhoods that predict distinct tissue ecologies in patients whose NSCLC has progressive disease or stable disease following treatment. We have identified potential biomarkers for treatment response, e.g., PD-1+PanCK+FoxP3. We demonstrate that PD tumors are characterized by a highly suppressed immune response (due to presence of FoxP3+ and PD-1+ immune cells) prior to treatment, whilst SD tumors showed a higher colocalization of PD-L1+CD3+CD8+ with tumor cells. Together, these findings indicate that there are distinct cellular compositions, both within the tumor regions and at the tumor boundaries, for PD versus SD patients. Furthermore, we use these distinct biomarkers to train support vector machines (SVMs) and neural network species distribution models (NN SDM) to predict disease progression.
Previous studies have shown that HDACi increases anti-PD-1 response by enhancing tumor immunogenicity, in part through upregulation of major histocompatibility complex (MHC) and T-cell chemokine expression (66). As HDACi alone is not effective in NSCLC (67), we hypothesize the TME associations described here are primarily associated with response to anti-PD-1. In this study, we have shown that both single cell and quadrat analysis methods can be synergized to predict disease progression, as well as visualize these predictions at the patient level through risk maps (Figure 5B). Further, these methods are generalizable across various cancer types as well as imaging modalities. We identify that disease progression predictions made by SVMs using PD-L1 alone have an overall mean accuracy of 63% (when using both cell segments and quadrats (at the patient level, Supplementary Table 3; at the cell/quadrat level per patient, Supplementary Table 1)). On the other hand, through our approach of using multiple markers along with PD-L1, the overall mean accuracy using SVM trained on cells is 88%, the NN SDM trained on quadrats is 75%, and the weighted SVM combining both cells and quadrats is 88.5% (Supplementary Table 1). These accuracies can be further improved in our ongoing follow up study with larger patient cohorts.
The predictions obtained from single cells are generally better than and comparable to those from the quadrats. This is mainly due to the marker resolution, i.e., in single cells, each marker is analyzed at the cellular level whereas in the quadrats, each marker represents a group of cells present within that quadrat. For research questions related to cellular co-localization and its impact on cellular processes (such as, cell signaling), one can rely on cellular analysis approaches. For studies seeking to summarize spatial interactions of cells (such as, interaction networks), quadrat approaches are appropriate. Nevertheless, it is important to note that the SVM and NN SDM predictions in this study, generated by utilizing all the markers along with their spatial context, are notably higher than those obtained by using PD-L1 alone (Supplementary Table 3).
Through these findings, our study reinforces that there are technical and clinical challenges in using PD-L1 as a standard diagnostic biomarker (5–8). Thus, there is a need for complementary diagnostic biomarkers along with PD-L1 for different cancer types and treatment (29). From both our single-cell and quadrat analyses, we find the tumor tissue immune architecture in pre-treatment biopsies is different in responding and non-responding tumors. Herein, we have shown that these pre-treatment ecologies of NSCLC patients are distinct, so much so that their spatial organization (CNs or QNs) can be used to predict response to treatment more accurately than PD-L1 alone. In PD patients, pre-treatment biopsies demonstrated exhausted T cells (PD-1) and immune suppressive Tregs (FoxP3), along with decreased CD8 T-cell infiltration. However, in SD patients, we observe increased CD8 infiltration in tumor regions. These findings are in alignment with studies that show PD tumors to exhibit more immune suppressed environments, and SD tumors to have higher numbers of infiltrated immune cells (68). Additionally, we verified these differences using a statistical model (Methods) that confirmed the ecological changes in pre-treatment biopsies from PD and SD patients.
The current study does have its limitations. One limitation is that it uses retrospective clinical trial data. In addition, it has a small sample size (FoVdata: n = 9 patients, 92 FoVs) with stable and progressive disease cases, but no complete or partial response cases. While selection bias may be present in FoVdata, we minimized this risk by selecting FoVs that had a tumor to stromal compartment ratio close to 0.5. Additionally: a) While we assumed PanCK+ to be indicative of tumor cells, not all PanCK+ are tumor cells, as non-tumor epithelial cells can also contribute to PanCK+ staining. Therefore, the presence of tumor cells in biopsies is determined by H&E staining. We note that the marker combinations we identified to differentiate responders versus non-responders are mostly immune signatures i.e., SD patients are identified using CD3+CD8+PD-L1+PanCK+, and PD patients using FoxP3+PD-1+PanCK+, where in both cases PanCK is a common denominator; b) The tumor biopsies used in this study were obtained from primary or different metastatic sites (Supplementary Table 2). Biopsy collection is typically based on sampling sites where the procedural risk to the patient is low. However, given the small numbers of individual sites, it would be difficult to assign spatial differences to specific sites, and the presence of different metastatic sites can be a confounding factor in the interpretation of results. Further, the patients enrolled in this Phase 1b trial had all failed previous immunotherapy. While this represents a similar patient population in that respect, previous immunotherapy treatments may impact the TME of patients in the current study. Also, our potentially significant clinical observation - that the pre-treatment tissue ecology is predictive of response or resistance to a treatment - is in the context of combination therapy, i.e., with HDAC inhibitor combined with a PD-1 inhibitor. In another ongoing study, we compare these predictive ecologies with a monotherapy cohort.
Our major motivations to embark the current study on this limited cohort were: a) to glean insights into the presence and effect of spatiotemporal ecologies in shaping treatment resistance as very little is known about the impact of drugs on evolving ecologies in NSCLC; b) to demonstrate the technical feasibility of our combined image analyses methodology; c) to scale our current computational pipeline from smaller datasets to utilizing larger datasets (both with increased number of patient samples and increased image sizes), and to extract meaningful biomarkers while bridging these scales; and d) as a computationally-driven study following the promising clinical observations inferred from FoVdata (29) where such predictive biomarkers were hypothesized to exist.
Acknowledging these limitations, we are currently expanding this study to validate the biomarkers from this initial study in a larger cohort of patients, including those from the ongoing randomized Phase II clinical trial (NCT02638090). This larger cohort has patients with a wider range of responses, including stable and progressive disease, and partial responders. Additionally, this trial has two arms, one in which anti-PD-1 naïve NSCLC patients are treated with either anti-PD-1 pembrolizumab monotherapy (Arm A, n=39), or with pembrolizumab + HDACi vorinostat (Arm B, n=39; similar to FoVdata in this study). Given the complexity of HDACi effect on the TME as well as uncertainty of eventual approval of HDACi treatment of NSCLC patients in combination with PD-1 inhibitors, our future work will be aimed at better understanding how the biomarkers might vary in a monotherapy versus combination therapy setting as well as across patient response categories. This will help interpret the variations in marker interaction and tumor architecture unique to each therapy. We have currently also undertaken mouse studies to mechanistically demonstrate our findings sine hypothesis testing studies are difficult to perform in these human specimens. Further, using Mistic and our statistical network analysis to study varying networks across treatment time points, we will investigate marker changes that characterize immune-cool tumors from immune-hot tumors. We will also develop mechanistic mathematical models to simulate response to treatment, which will aid in clarifying why the markers are predictive.
The greatest novelty and impact of this study are the computational methods and the pipeline, combining both cell segments and quadrats, which point to a role of spatial neighborhood (both, CNs and QNs) compositions in immunotherapy response. Specifically, the associative aspects such as co-occurrence of FoxP3+PD-1+PanCK+ in CNs and QNs of PD patients are unique and suggestive of new strategies for predicting immunotherapy response. In summary, we use single cells and quadrats combined with spatial statistics, machine learning, and deep learning, to perform in-depth analyses of multiplexed images (FoVs) based on the frequency, phenotype, and spatial distribution of immune and tumor cells within the immune landscape of advanced/immunotherapy-refractory NSCLC. We demonstrate that PD tumors are characterized by a highly suppressed immune response (due to presence of FoxP3+ and PD-1+ immune cells) prior to treatment, whilst SD tumors showed a higher colocalization of PD-L1+CD3+CD8+ with tumor cells. This is consistent with previous literature that in NSCLC patients, spatial proximity of tumor and FoxP3 cells results in exhausted immune cells and has been associated with diminished survival (57), while infiltration of CD3+ and CD8+ T cells along with anti-PD-L1 treatment is associated with better overall survival (65,68). These, therefore, validate our cell- and quadrat-based approaches, and further indicate that our dual approaches are applicable beyond NSCLC (including solid tumors such as ovarian cancer (69), endometrial cancer (70), and glioblastoma (71)) as well as the specific subset of markers discussed here, and will be even more relevant for use with new markers where the underlying biology is not clear.

Supplementary Material
1234567891011