Pancreatic ductal adenocarcinoma (PDAC) is a highly aggressive disease with poor 5-year survival rates, necessitating identification of novel therapeutic targets. Elucidating the biology of the tumor immune microenvironment (TiME) can provide vital insights into mechanisms of tumor progression. In this study, we developed a quantitative image processing platform to analyze sequential multiplexed IHC data from archival PDAC tissue resection specimens. A 27-plex marker panel was employed to simultaneously phenotype cell populations and their functional states, followed by a computational workflow to interrogate the immune contextures of the TiME in search of potential biomarkers. The PDAC TiME reflected a low-immunogenic ecosystem with both high intratumoral and intertumoral heterogeneity. Spatial analysis revealed that the relative distance between IL10+ myelomonocytes, PD-1+ CD4+ T cells, and granzyme B+ CD8+ T cells correlated significantly with survival, from which a spatial proximity signature termed imRS was derived that correlated with PDAC patient survival. Furthermore, spatial enrichment of CD8+ T cells in lymphoid aggregates was also linked to improved survival. Altogether, these findings indicate that the PDAC TiME, generally considered immuno-dormant or immunosuppressive, is a spatially nuanced ecosystem orchestrated by ordered immune hierarchies. This new understanding of spatial complexity may guide novel treatment strategies for PDAC.

Significance:

Quantitative image analysis of PDAC specimens reveals intertumoral and intratumoral heterogeneity of immune populations and identifies spatial immune architectures that are significantly associated with disease prognosis.

Pancreatic cancer is projected to be the second leading cause of cancer-related death worldwide by 2030 (1). Pancreatic ductal adenocarcinoma (PDAC) is the most prevalent type of pancreatic cancer accounting for over 90% of all pancreatic malignancies diagnosed (2). PDAC is highly aggressive and demonstrates a strikingly poor 5-year survival rate of approximately 10% in the United States (3). Surgical resection is generally considered the only curative treatment for localized PDAC, when clinically a viable option (4). However, more than 80% of tumors are nonresectable at the time of diagnosis, and most patients who benefit from surgery may eventually recur (5, 6). As a result, developing novel treatment strategies for PDAC is an urgent yet unmet medical need.

In recent decades, advancements in immunotherapy have led to some modest improvements in PDAC treatment. Currently, there are several clinical trials evaluating anticancer efficacy of immune checkpoint inhibitors and cancer vaccines (7–9). For instance, Jaffee and colleagues evaluated neoadjuvant GVAX with or without low-dose cyclophosphamide in patients with resectable PDAC and reported evidence of proinflammatory immunologic changes within PDAC tumors in response to GVAX (10). Bockorny and colleagues evaluated the synergistic effects of a CXC chemokine receptor 4 antagonist (BL-8040) in combination with PD-1 inhibition (pembrolizumab) and chemotherapy in PDAC, and reported that dual blockade improved clinical outcomes (11). Despite significant survival benefits associated with aforementioned treatments, the overall response rates were modest. Therefore, identifying novel therapy targets and predictive biomarkers will guide new clinical strategies for PDAC treatment. In general, searching for novel biomarkers relies on the notion that an activated tumor immune microenvironment (TiME) favors positive response to immunotherapy. The PDAC TiME has traditionally been considered T cell suppressive, thus presenting considerable barriers to establishing effective biomarkers for patient stratifications (12, 13). As a result, a more nuanced understanding of the biology of TiME is warranted for rational disease control.

Recent developments in high-dimensional multiplexed imaging technologies allow for simultaneous phenotyping of multiple cell populations with their spatial domains preserved, thus facilitating a better understanding of TiMEs. For instance, previous studies have reported the prognostic significances of spatial landscapes for a series of immune populations featuring B cells, T cells, myelomonocytes, and cancer-associated fibroblasts in PDAC TiMEs (14–18). Although these studies shed light on features of distinct immune populations, how they cooperatively impact antitumor immunity remains unknown. Indeed, antitumor immunity requires structured, spatially coordinated interplay between components of the TiME. Schürch and colleagues reported that coupling of distinct cellular neighborhoods alters antitumor behaviors in colorectal cancer (19); similarly, Jackson and colleagues described higher-order interactions between immune phenotypes, namely spatial communities, to be correlated with clinical outcomes for patients with triple-negative breast cancer (20).

To this end, we deployed a multiplexed IHC (mIHC)-based computational pathology framework and applied it to tumor resections from archival single-institution patient cohorts with PDAC consisting of a treatment-naïve group and a presurgically treated group. This framework is a significant expansion of our previous studies (21–23), while also building on studies from other groups (19, 24). The imaging system herein utilized a comprehensive lymphoid and myeloid marker antibody panel to profile leukocyte populations in histopathologically defined regions within PDAC TiMEs. An image processing pipeline was then implemented to capture the multiscale immune architectures. Altogether, such workflows enable assessment of complex TiME biology and relationship to clinical outcomes, with potential applications to informing effective patient stratification and rationally designed clinical trials.

Data acquisition with mIHC

mIHC staining was performed as described previously (25). mIHC-stained mages utilized in the current study included specimens from cohort 1 of the treatment-naïve group and presurgically treated group, stained previously with the functional panel (Supplementary Table S1). A subset of biomarkers from this panel was used to evaluate higher-order cell types with greater expected abundance. In the current study, the treatment-naïve group was used for biomarker discovery, denoted as discovery cohort; presurgically treated group was used for validation, denoted as validation cohort. Patient characteristics for both cohorts were summarized in Supplementary Tables S2 and S3. Images were acquired using an Aperio AT2 scanner (Leica Biosystems) at 20× magnification. Human PDAC specimens were obtained in accordance with the Declaration of Helsinki and were acquired with written informed consent by Institutional Review Boards at Johns Hopkins University (Baltimore, MD) and Oregon Health and Science University (Portland, OR).

Image processing and analysis

Regions were selected and annotated in Aperio ImageScope (Leica Biosystems, RRID: SCR_014311) and saved as an XML. Selected regions of each stain from each tissue were extracted and registered to the final hematoxylin stain in MATLAB version R2018b (The MathWorks, Inc., RRID: SCR_001622) using detectSURFFeatures to identify the strongest matching key points from each of the RGB channels and selecting the channel with the highest matching key points. Region of interest (ROI) areas for each patient are summarized in Supplementary Table S4. A geometric transformation for each mIHC stain to the corresponding hematoxylin was estimated on the basis of similarity of matched key points with outliers excluded using the M-estimator Sample Consensus algorithm (MSAC). The transformation was applied to each RGB channel and merged to create the composite registered RGB image region. Color deconvolution was performed in FIJI (RRID: SCR_002285) by converting the RGB image to CMYK and using the Y channel to extract AEC chromogenic signal. Nuclei segmentation was watershed based and performed in FIJI by first using color deconvolution (H AEC) to separate hematoxylin signal, followed by smoothing, background subtraction and Otsu thresholding. Cell Profiler was used to quantify mean intensity of each segmented cell, along with area and centroid location (Fig. 1A). The resulting feature matrix was used in FCS Express 7 Image Cytometry RUO (De Novo Software) where each tissue region was manually gated for single-cell classification and visually validated with live rendering of the segmentation mask on top of the signal extracted image (Fig. 1B; Supplementary Fig. S1). Phenotyping results as well as the registered pseudocolored multi-channel images were visualized in Fig. 1C.

Figure 1.

A, Overview of mIHC staining and analysis pipeline including ROI selection from a PDAC resection, image registration, nuclei segmentation, color deconvolution, and single-cell measurements. B, Table describing biomarkers used for identifying cell types and functional states (left) alongside a visual of the gating strategy (right). The full gating strategy and description can be found in Supplementary Fig. S1. C, Gated phenotype map shows cells classified from manual gating (left); same region shown as a pseudofluorescent image (middle). White box indicates area of the right panel image; zoomed in region from middle panel image showing cells with biomarkers used for identification. White arrows (left to right) point to a CD3 T cell, CD20 B cell, CD8 T cell, and CD68 myelomonocytic cell. (Schematics created with BioRender.com.)

Figure 1.

A, Overview of mIHC staining and analysis pipeline including ROI selection from a PDAC resection, image registration, nuclei segmentation, color deconvolution, and single-cell measurements. B, Table describing biomarkers used for identifying cell types and functional states (left) alongside a visual of the gating strategy (right). The full gating strategy and description can be found in Supplementary Fig. S1. C, Gated phenotype map shows cells classified from manual gating (left); same region shown as a pseudofluorescent image (middle). White box indicates area of the right panel image; zoomed in region from middle panel image showing cells with biomarkers used for identification. White arrows (left to right) point to a CD3 T cell, CD20 B cell, CD8 T cell, and CD68 myelomonocytic cell. (Schematics created with BioRender.com.)

Close modal

Characterization of spatial heterogeneity

A previously introduced spatial form of Shannon entropy (ESP) was implemented to quantify the spatial heterogeneity (22). ESP is formally defined as:
formula
where | $d_i^{{\rm{int}}}$ | denotes the average Euclidean distance between all cell centroids of type i; | $d_i^{{\rm{ext}}}$ | represents the average distance between all cells of set i and cells of all other types; pi is the percentage of type i within the core.

Evaluation of spatial correlations

To evaluate the spatial correlation between each pair of immune cell phenotypes, we utilized the spatial G(r) function (Gcross) to model the clustering pattern. Gcross computed the likelihood of seeing at least one cell of phenotype j around a cell of phenotype i, denoted as Gcross(i-j) within a series of radii. Here, we chose the maximum evaluable radius as 50 μmol/L, to ensure clustering patterns were fully captured. Gcross values were estimated at every 0.05 μmol/L to maximize continuity, then stratified into four bins based on the evaluated radius: (i) 10–20 μmol/L, (ii) 20–30 μmol/L, (iii) 30–40 μmol/L, and (iv) 40–50 μmol/L. The bin selection was based upon previous studies that evaluated direct cell–cell interactions and enables a plausible profiling of clustering (26–28). G(r) was computed using “Gcross” function from R package “spatstat” (29).

Enrichment score

ROIs with at least 10 CD4+ T cells and 10 myelomonocytes were subjected to analysis. The enrichment score was computed for both CD4+ T cells and myelomonocytes on a per-ROI basis. For example, for myelomonocytes, the procedure is as follows: For each evaluable ROI, we computed the distance from all myelomonocyte centroids to CD4+ T-cell centroids. Next, myelomonocytes located within 50 μmol/L of at least one CD4+ T cell were identified and defined as neighbors of those CD4+ T cells. Then, for every protein x ∈ X = {PD-L1, PD-1, IL-10, GZMB, Ki67, EOMES, ICOS}, the following procedures were conducted: (i) the ratio of neighbors positive for x was computed, denoted as R0 (ii); repeatedly randomize the location of myelomonocytes. For the ith randomization, the ratio was also computed and denoted as Ri, where i ∈ {1, 2, 3, …., n} (iii); Enrichment score was formally defined as the proportion of R0 greater than Ri. The protocol was also implemented in reverse with CD4+ T cells as neighbors of myelomonocytes. Here, distances were computed using “dist2” function from R package “flexclust”; randomization was implemented using “runifpoint” function from R package “spatstat”.

IL10+ myelomonocyte RiskScore

ROIs with at least one IL10+ myelomonocyte, at least one PD-1+ CD4+ T cell, and at least one granzyme B-positive (GZMB+) CD8+ T cell were subjected to analysis. For every IL10+ myelomonocyte, the distances to its nearest PD-1+ CD4+ T cell and GZMB+ CD8+ T cell were computed and denoted as d1 and d2, respectively. The IL10+ myelomonocyte RiskScore (imRS) is formally defined as:
formula
RiskScore, also termed as imRS in this study, reflects the balance of IL10+ myelomonocyte between high immunosuppressive risk and low immunosuppressive risk. Close proximity to PD-1+ CD4+ T cells will result in d1 decreasing and d2 increasing, therefore resulting in low imRS.

CD8+ T-cell–B-cell network morphometrics

ROIs with at least 100 total CD8+ T cells and 100 total B cells were subjected to analysis. For each evaluable ROI, cell dataset were truncated to include CD8+ T cells and B cells only. A HDBSCAN algorithm was first applied to identify clusters. For each cluster, Voronoi tessellation was then generated. Here, each cell was treated as a node. For each pair of Voronoi polygons that share at least one border, the associated nodes were connected to formulate a link. The set of nodes and links formed the cluster network. To facilitate the computation of cluster morphometrics, the α-shape and convex hull of each cluster was also derived. In this study, four morphometrics were gauged:

Cluster area: measures the area of convex hull.

Circularity: measures the roundness of the object. Here, circularity is defined as:
formula
where Aconc is the area of the α-shape, Pconc is the perimeter of the concave hull.
Eccentricity: measures the degree of the object deviation from being circular. Here, we computed the fitting ellipse by first computing the covariance matrix from the cluster point clouds. Assuming the point clouds follow Chi-squared distribution | ${\rm{Q}}\sim{{\rm{x}}}^2( k )$ |⁠, the eigenvectors can be calculated from the covariance matrix yielding the orientations. Thus, the major and minor axis length can be computed as:
formula
respectively, where λ1 and λ2 are eigenvalues of the covariance matrix. According to the definition, the ellipse represents the contour that covers 95% of point clouds. Thus, the eccentricity is computed as:
formula
Convexity: measures the curvature of the object. Here, convexity is defined as:
formula
where Aconv is the area of convex hull.

HDBSCAN algorithm was implemented using “hdbscan” function from R package “dbscan”; neighboring Voronoi polygons were identified using “voronoi_adjacentcy” function from R package caramellar; α-shape was computed using the “ashape” function from R package “alphahull”; concave hull was computed using the “chull” function from R package “grDevice”; covariance matrices from the cluster point clouds were computed using “cov.wt” function from R package “stats”; eigenvectors were computed using “eigen” function from base R package; χ2 distribution at 95% confidence interval (CI) was computed using “qchisq” function from R package “stats”.

Statistical analysis

Two-sided Wilcoxon rank-sum test was performed for pairwise comparisons using “wilcox.test” function from R package “stats”. FDR adjustment for P values were performed for multiple comparisons using “p.adjust” from R package “stats” and adjusted P < 0.05 was considered significant. Pearson correlation coefficients were computed to assess the linear correlation between two sets of data using “cor” function from R package “stats”; χ2 test was performed to assess the correlation between categorical variables using “chisq.test” function from R package “stats”; heatmaps and hierarchical clustering were generated using “Heatmap” function from R package “ComplexHeatmap” (RRID: SCR_017270).

Data availability

The data supporting the findings of this study (single-cell and patient data) are available online at Zenodo (RRID: SCR_004129, https://zenodo.org/record/6416102#.Yky33jfMIrY). Raw images are available from corresponding author upon reasonable requests.

Code availability

The codes for computational methods are made publicly available at GitHub: https://github.com/popellab/PDAC-mIHC-processing-pipeline. And an online-executable capsule is made available at Code Ocean (https://codeocean.com/capsule/2525016/tree/v1) and the link is provided in the GitHub repository.

Patient stratifications

Acknowledging the fact that PDAC is highly aggressive and has poor prognosis, we stratified treatment-naïve patients into two groups reflecting long-term survivors, that is, above median, and short-term survivors, that is, below median [based on overall survival (OS) days; threshold = 619 days], for downstream evaluations. A merit of such approach is that survival differences can be maximized therefore the subtle spatial heterogeneity can be scaled to facilitate the prognostic marker discovery. Identifying differences first and then correlating to survival may weaken the significance due to misclassification of certain patients. While different splitting strategies have been employed by previous studies for biomarker discovery (25, 30, 31), there is no consensus on the optimal split approach. In this study, 50% was selected as the cut-off point produced two groups of comparable size, also reducing the risk of overfitting for the downstream analysis. The resulting long-term group consisted of 22 patients with a median OS = 832 days, whereas the short-term group consisted of 23 patients with a median OS = 313 days (Fig. 2A).Using 50% as the cut-off point produced two groups of comparable size, reducing the risk of sample limitations in downstream analysis. In addition, this stratification also exhibits statistical significance between survival groups at the same level (log-rank test P < 0.001), as opposed to alternative stratification strategies (Supplementary Fig. S2).

Figure 2.

Characterization of the PDAC immune ecosystem. A, A discovery cohort was split into short-term and long-term groups by OS using 50% as the cut-off threshold (long-rank test). B, Density of each immune phenotype was compared between long-term and short-term groups; CD45+ other: CD45+ population other than identified immune phenotypes (B cell, myelomonocyte, CD4+ T, and CD8+ T); ns, not significant (Wilcoxon rank-sum test). C, Immune cell densities of all patients were sorted in ascending order and the associated compositions at both patient and ROI level were visualized. D, Correlation analysis identified strong negative correlations between B-cell and myelomonocyte density, and positive correlation between myelomonocyte and other immune phenotypes (Pearson coefficient). E, Correlation heatmaps for all immune phenotype pairs showed similar patterns between short- and long-term groups (Pearson coefficient). F, Immune-infiltration levels were quantified into low, medium, and high categories and compared within each category; ns (Wilcoxon rank-sum test); the distribution among all categories was also compared between short- and long-term groups (χ2 test). G, Exemplary ROI with high and low spatial Shannon entropy. H, Spatial Shannon entropy was computed for each ROI and then averaged to the patient level for comparison between short- and long-term groups; ns (Wilcoxon rank-sum test).

Figure 2.

Characterization of the PDAC immune ecosystem. A, A discovery cohort was split into short-term and long-term groups by OS using 50% as the cut-off threshold (long-rank test). B, Density of each immune phenotype was compared between long-term and short-term groups; CD45+ other: CD45+ population other than identified immune phenotypes (B cell, myelomonocyte, CD4+ T, and CD8+ T); ns, not significant (Wilcoxon rank-sum test). C, Immune cell densities of all patients were sorted in ascending order and the associated compositions at both patient and ROI level were visualized. D, Correlation analysis identified strong negative correlations between B-cell and myelomonocyte density, and positive correlation between myelomonocyte and other immune phenotypes (Pearson coefficient). E, Correlation heatmaps for all immune phenotype pairs showed similar patterns between short- and long-term groups (Pearson coefficient). F, Immune-infiltration levels were quantified into low, medium, and high categories and compared within each category; ns (Wilcoxon rank-sum test); the distribution among all categories was also compared between short- and long-term groups (χ2 test). G, Exemplary ROI with high and low spatial Shannon entropy. H, Spatial Shannon entropy was computed for each ROI and then averaged to the patient level for comparison between short- and long-term groups; ns (Wilcoxon rank-sum test).

Close modal

Quantitative image analysis reveals poor immunogenicity in PDAC TiME

We applied automated image analysis to profile the immune populations for all patients. Specifically, we quantified the cell densities per ROI for every immune cell phenotype and averaged to the patient level. Region properties associated with each ROI are summarized in Supplementary Table S4. Average densities of each cell type were compared between survival groups, but there were no significant differences found (Wilcoxon rank-sum test P > 0.05; Fig. 2B). Nonetheless, all immune cell phenotype densities showed moderate to large variability, thus we further sought to analyze whether such variations were driven by different patterns reflecting immune subset compositions. We first sorted leukocyte densities for each patient from low to high, and then calculated the percentage of each leukocyte subset at the per-patient and per-ROI level. We then correlated overall immune cell densities with each immune cell phenotype. Previous studies reported that the compositions of immune subsets were highly correlated with overall immune cell abundances in triple-negative breast cancer (24). Specifically, CD4+ T cells were enriched whereas myelomonocytes were sparse in patients with high-leukocyte infiltration. We observed the same significance for CD4+ T cells but not myelomonocytes in this study for both short- and long-term groups (Supplementary Fig. S2A–S2F). However, correlations were identified between immune cell populations densities (Fig. 2C); for instance, myelomonocyte densities were positively correlated with other CD45+ density, but were negatively correlated with B-cell densities (Fig. 2D) in both short- and long-term groups. In fact, the density correlation heatmaps revealed that the two groups shared similar correlation profiles. That is, whichever cell correlation pair that was significant in the short-term group, was also significant in the long-term group, with the only difference reflected in significance level (Fig. 2E). Next, we stratified all patients into tertiles reflecting low-, medium-, and high-leukocyte infiltration. We then decomposed each tertile based on the patient's survival group. Results demonstrated no correlation between immune infiltration levels and survival, as short- and long-term survivors were evenly distributed across all categories (χ2P = 0.9149; Fig. 2F). The immune infiltration profiles were further gauged by evaluating compositional heterogeneity. Here, a spatial form of Shannon entropy (Esp) was implemented. Briefly, the metric attributes the sources of spatial entropy to the distance between cells and the diversity of cell species and their abundance, that is, close proximity between cells of different types and enrichment of unique cell phenotypes both lead to the increase of Esp (Fig. 2G; Materials and Methods). Previously we reported that high Esp, that is, high spatial immune heterogeneity, is associated with response to immune therapy in the context of hepatocellular carcinoma (32). The results here however show no significant difference between short-term and long-term groups (Fig. 2H). Univariate analysis on the aforementioned first-order features included immune population densities and Shannon entropy (Supplementary Fig. S3A). We first aggregated features from core level to patient level by taking the mean. Univariate Cox regression analysis revealed that none of the features significantly associated with survival (P > 0.05) with four features conferred HRs = 1, suggesting these features provide zero risk reduction. Similarly, none of the features was prognostic in multivariate Cox regression analysis (Supplementary Fig. S3B). It is noteworthy that Shannon entropy yields a relatively higher HR (mean = 4.21); however, the 95% CI extended over a wider range, suggesting less clinical interest.

Together, these findings may potentially indicate that characterization of TiME at low resolution is not able to make distinctions among patients and it requires higher-order, more complex characterizations to establish potential biomarkers in separating short- versus long-term survivors, though validations on additional cohort is required to eliminate overfitting issues due to small sample size.

Spatial correlation analysis of PDAC delineates orchestrated architecture of immune cell phenotypes

To further model the spatial organization of the immune landscape in PDAC, we applied spatial G(r)-cross function (Gcross) to evaluate cell-cell colocalizations. At this point, CD45+ other cells were not included because their phenotype and functional states were not determined. Gcross was applied to each pair of immune cell phenotypes with both directions (Fig. 3A; Materials and Methods). For instance, Gcross (CD4+–myelomonocyte) evaluates the likelihood of observing CD4+ T cells colocalized with at least one myelomonocyte within a series of radii. High Gcross (CD4+–myelomonocyte) values indicate clustering of myelomonocytes around CD4+ T cells, reflected by steep climbing curves along with the radius gradient (Fig. 3B). Because Gcross is not a symmetrical metric, cell type pairs were evaluated in both directions, where CD4+ T cells or myelomonocytes, were evaluated as the focal cell type. Gcross values were computed with a radius gradient and then binned into four categories for every 10 μmol/L (Materials and Methods) and compared between survival groups within each category. Comparisons identified three pairs that were significant in each category spanning all distances (FDR-adjusted Wilcoxon rank-sum test P < 0.05; Fig. 3C; Supplementary Table S5). Here, we focused on two-way clustering involving CD4+ T cells and myelomonocytes, indicating a spatial mutual signaling axis, and one-way clustering from B cells to CD8+ T cells, highlighting lymphoid aggregate (LA)-like structures. Indeed, the prognostic significance of both patterns has been reported in broader contexts (33–35). Specifically, Gcross values were significantly higher in long-term survivors for both directions between CD4+ T cells and myelomonocytes, while only Gcross (B cell–CD8+ T) was significantly higher in short-term survivors. These results indicate that CD4+ T cells and myelomonocytes exhibited a higher level of mutual colocalization in long-term survivors, while short-term survivors tended to reflect increased density of CD8+ T cells in areas containing B cells (Fig. 3D).

Figure 3.

Spatial correlation analysis identifies topological features associated with prognosis. A, schematic of spatial G(r) function (Gcross). B, Exemplary Gcross curves computed for ROIs associated with short-term (red) and long-term (blue) groups. C, Comparisons of Gcross values between short- and long-term groups at four distance intervals identifies three spatial clustering pairs; ****, P < 0.0001 (Wilcoxon rank-sum test). D, Aforementioned spatial clustering pairs were summarized into a two-way clustering involving myelomonocytes and CD4+ T cells, and a one-way clustering involving B cells to CD8+ T cells. E, Functional marker expressions across immune phenotypes. F, Schematic of computation of enrichment score. G, Enrichment score computed for each functional marker on both myelomonocytes and CD4+ T cells. H, Schematic of imRS, derived from the relative distances from IL10+ myelomonocytes (mean of five nearest neighbors, if any) to PD-1+ CD4+ T cells and GZMB+ CD8+ T cells (mean of five nearest neighbors, if any). I, imRS was computed at per-cell basis and compared between short- and long-term groups; ****, P < 0.0001 (Wilcoxon rank-sum test). J, Exemplar topological graph of IL10+ myelomonocytes, PD-1+ CD4+ T cells, and GZMB+ CD8+ T cells in short- and long-term groups. (Schematics created with BioRender.com.)

Figure 3.

Spatial correlation analysis identifies topological features associated with prognosis. A, schematic of spatial G(r) function (Gcross). B, Exemplary Gcross curves computed for ROIs associated with short-term (red) and long-term (blue) groups. C, Comparisons of Gcross values between short- and long-term groups at four distance intervals identifies three spatial clustering pairs; ****, P < 0.0001 (Wilcoxon rank-sum test). D, Aforementioned spatial clustering pairs were summarized into a two-way clustering involving myelomonocytes and CD4+ T cells, and a one-way clustering involving B cells to CD8+ T cells. E, Functional marker expressions across immune phenotypes. F, Schematic of computation of enrichment score. G, Enrichment score computed for each functional marker on both myelomonocytes and CD4+ T cells. H, Schematic of imRS, derived from the relative distances from IL10+ myelomonocytes (mean of five nearest neighbors, if any) to PD-1+ CD4+ T cells and GZMB+ CD8+ T cells (mean of five nearest neighbors, if any). I, imRS was computed at per-cell basis and compared between short- and long-term groups; ****, P < 0.0001 (Wilcoxon rank-sum test). J, Exemplar topological graph of IL10+ myelomonocytes, PD-1+ CD4+ T cells, and GZMB+ CD8+ T cells in short- and long-term groups. (Schematics created with BioRender.com.)

Close modal

Spatial proximity signature of CD4+ T cell, CD8+ T cell, and myelomonocyte subpopulations associated with PDAC survival

To further decipher the underlying biology of the association between survival and mutual clustering of CD4+ T cells and myelomonocytes, we studied each pair of the relevant cell types and presence of key mediators they express (e.g., PD-L1, PD-1, IL10, granzyme B, Ki67, EOMES, ICOS) to evaluate their influence on the clustering. We first quantified the composition for each immune cell phenotype within the cohort and revealed similar patterns across cell types: a majority of myelomonocytes, B cells, CD4+ and CD8+ T cells, were IL10 positive, with GZMB positivity in a minority of these subsets (Fig. 3E). On the basis of this, we investigated the enrichment score to account for such potential bias caused by the skewness of positive-protein distribution. For each ROI, we repeatedly randomized the geolocations of cell of type X positive for protein α and recorded the ratio of positive cells within 50 μmol/L of all cells of type Y (neighbors). The enrichment score SEn of the ROI is formally defined as the ratio of the times that the randomization ratio is smaller than the ratio derived from the real point pattern (Fig. 3F; Materials and Methods). The benefit of this approach lies in the fact that the proportion of cells positive for each protein within the ROI is kept consistent throughout the randomization, resulting in the global proportion bias becoming a recurrent constant factor. Under such conditions, high SEn indicates that more neighbors are observed than at random. Given that the global proportion always has the same weight, such close proximity is likely driven by an underlying biological process rather than bias or coincidence. We computed SEn for all ROIs and found that PD-1 positivity had the highest mean for CD4+ T cells (median SEn = 1, mean SEn = 0.65) whereas IL10 had the highest mean for myelomonocytes (median SEn = 1, mean SEn = 0.70) compared with other markers; the conclusion held under different simulation rounds (Supplementary Table S6), indicated that PD-1+ CD4+ T cells and IL10+ myelomonocytes were key contributing factors to the clustering (Fig. 3G).

Previously, Phillips and colleagues developed SpatialScore and reported relationships between spatial proximity of immune subpopulations and response to immunotherapy in cutaneous T-cell lymphoma (36). Of note, the concept of SpatialScore is to characterize the balance between CD4+ T-cell effector and immunosuppressive activity. Here, we propose a modification to represent the relative distance of PD-1+ CD4+ T cell to IL10+ myelomonocyte and GZMB+ CD8+ T cell. To recapitulate this proxy, we included GZMB+CD8+ T cells given that IL10+ myelomonocyte have been reported to inhibit CD8+ T-cell cytotoxicity by direct and indirect mechanisms (immunosuppressive activity; refs. 37, 38). In this study, the prediction is that closer proximity between CD4+ T cells and myelomonocytes will correlate with better survival (inhibited immunosuppressive activity). In practice, we first identified the mean distance of each IL10+ myelomonocyte (evaluable interaction pairs n = 13,451) to its k-nearest PD-1+ CD4+ T cell and k-nearest GZMB+ CD8+ T cell for each ROI (evaluable n = 141), then computed the relative distance ratio to generate the IL10+ imRS (Fig. 3H; Materials and Methods). Here, we set the threshold k = 5 (or equals to the total number of candidates if not abundant) to reduce the sensitivity to outlier or misclassified cells. Quantitative results revealed that the imRS was significantly higher in short-term survivors (mean imRS = 0.44) than long-term survivors (mean imRS = 0.23). It is also noteworthy that the distribution of imRS from long- and short-term groups indicated a distinct IL10+ myelomonocyte behavior, as they were clearly separated by the 0.5 threshold (Fig. 3I). Visualizations further corroborated such finding: we linked IL10+ myelomonocyte to their nearest PD-1+ CD4+ T cell and GZMB+ CD8+ T cell on a per-cell basis and visualized the proximity of IL10+ myelomonocyte adjacent to PD-1+ CD4+ T cells, while remaining distal to GZMB+ CD8+ T cells in long-term survivors (Fig. 3J). Given that the sample size is relatively modest, we sought to test whether such signal was biased that originated from a single ROI or patient. Therefore, we further performed exclusion analysis on the imRS data. We iteratively excluded data from a patient (evaluable N = 31). For each iteration, imRS between short- and long-term survivors was compared using Wilcoxon rank-sum test and P values were recorded. Results showed that the statistical significances retained over all iterations under both conditions, therefore confirmed the robustness of imRS (Supplementary Table S7). Collectively, we demonstrated that the spatial proximity signature of IL10+ myelomonocyteto PD-1+ CD4+ T cell and GZMB+ CD8+ T cell correlates with patient survival in PDAC.

Compositions of B-cell and CD8+ T-cell clustering associated with PDAC survival

Next, we sought to further interrogate the second significant pattern identified previously, that is, one-way colocalization between CD8+ T cells and B cells. Visual inspections indicated that the clustering patterns resembled LAs (Fig. 4A). Here, we quantified the CD8+ T–B-cell aggregates by applying the HDBSCAN algorithm to identify aggregates (n = 125). For all cells within each identified aggregate, a network was constructed by first generating the Voronoi tessellation and then connecting all neighboring cells (Materials and Methods). For each network, we measured the following metrics: convex hull area, CD8+ T-cell density, B-cell density, circularity, eccentricity, and convexity. While results revealed that cluster morphometrics did not associate with survival (Wilcox rank-sum test P > 0.05), CD8+ T-cell density per aggregate was significantly higher in long-term survivors (Fig. 4B). In addition, we aggregated the density from per-LA basis to patient level and fitted a univariate Cox regression model. Results revealed that the density feature significantly contributed to the risk reduction (Wald test P = 0.03) with HR = 0.997 (95% CI, 0.9938–0.9997). Correlation analysis using Spearman method demonstrated same pattern that the density positively correlated with OS (ρ = 0.3658, P = 0.35; Supplementary Fig. S3C).

Figure 4.

Spatial correlation analysis revealed prognostic significance of LAs. A, First-order properties (cell densities) and morphometrics were computed and compared between short- and long-term groups; **, P < 0.005; ns, not significant (Wilcoxon rank-sum test). B, Exemplar circular (criteria: eccentricity < 0.8 and convexity > 0.8 and circularity > 0.5) and elongated (criteria: eccentricity > 0.8 or convexity < 0.3 or circularity < 0.3) LA. C, Hierarchically clustered heatmap of marker expression compositions in LAs, for example, PD1._CD8T represents PD-1+ CD8+ T cells. Bar plot shows the B-cell and CD8+ T-cell densities for each LA, respectively. D, Comparison of cluster distribution in different survival groups (χ2 test).

Figure 4.

Spatial correlation analysis revealed prognostic significance of LAs. A, First-order properties (cell densities) and morphometrics were computed and compared between short- and long-term groups; **, P < 0.005; ns, not significant (Wilcoxon rank-sum test). B, Exemplar circular (criteria: eccentricity < 0.8 and convexity > 0.8 and circularity > 0.5) and elongated (criteria: eccentricity > 0.8 or convexity < 0.3 or circularity < 0.3) LA. C, Hierarchically clustered heatmap of marker expression compositions in LAs, for example, PD1._CD8T represents PD-1+ CD8+ T cells. Bar plot shows the B-cell and CD8+ T-cell densities for each LA, respectively. D, Comparison of cluster distribution in different survival groups (χ2 test).

Close modal

To further characterize the CD8+ T-cell–B-cell aggregates, we computed the percentage of CD8+ T cells and B cells that are positive for each functional marker and then applied hierarchical clustering to group networks with similar compositions (Fig. 4C). Results indicated that the clustering algorithm identified three major clusters: in cluster 1, all networks were dominated by IL10+ B cells and CD8+ T cells; in cluster 2, an accumulation of EOMES positive cells was identified and the majority of aggregations associated with short-term survivor samples (30/36) from 11 patients; in cluster 3, diversified dominancy was observed. Of note, the distribution of clusters indicated a significant association to survival (χ2 test P = 0.00018; Fig. 4D).

Validation of the spatial TiME architectures with a presurgically treated PDAC cohort

In this study, we selected an independent cohort (evaluable N = 12) that had received presurgically treatment to validate our previous findings. The cohort composed of 6 patients treated with (chemo)radiation and chemotherapy (class 1), 5 patients treated with chemotherapy only (class 2), and 1 patient treated with radiation or chemotherapy only (class 3), prior to resection surgery (Supplementary Table S3). We first characterized the differences in survival and TiME compositions between classes. Class 3 was removed from this analysis since it only contains 1 patient. Importantly, we observed trend toward better prognosis in class 2 patients (mean OS days = 809.4) compared with class 1 (mean OS days = 531.8) but such survival advantage is not significant (log-rank test P = 0.07; Supplementary Fig. S4A). TiME composition analysis revealed no significant differences in immune populations (Wilcoxon rank-sum test P > 0.05). However, it is noteworthy that class 2 patients tended to have global higher numbers of CD4+ T cells and CD8+ T cells, suggesting that chemotherapy may be more effective in activating the antitumor immunity (Supplementary Fig. S4B).

Among 12 patients, 3 had marked response/minimal residual disease (grade 1), 3 had moderate response (grade 2), and 3 had poor or no response (grade 3). In this study, we define patients with grade 1 and 2 as responders (R) and patients with grade 3 as nonresponders (NR).While we recognize that this presurgically treated validation cohort does not represent the discovery cohort in terms of clinical scenario, due to limitations in data availability, the neoadjuvant cohort provided the best option for further investigation. Importantly, we demonstrated that the presurgical treatment-naïve long-term and presurgically treated validation groups both had significantly longer OS than the treatment-naïve short-term group, with no statistical difference in OS between the treatment-naïve long-term and presurgically treated groups (Fig. 5A). Therefore, we sought to determine whether the presurgically-treated cohort also resembled the long-term survival group in terms of the spatial features described above. We first computed the imRS for all IL10+ myelomonocytes (evaluable interaction pairs n = 102) using the same ROI selection criteria and protocol (evaluable ROI n = 16). Results revealed that the mean imRS in the validation cohort was significantly lower than both short-term (Wilcoxon rank-sum test P < 0.001) and long-term survivors (Wilcoxon rank-sum test P < 0.05), validating our hypothesis that the imRS is indeed associated with improved OS (Fig. 5B). Again, visualization recapitulated the same proximity pattern observed in long-term survivors (Fig. 5C). However, it is noteworthy that presurgically treated group that bear significantly higher imRS seems to not benefit from further improved survival (log-rank test P = 0.11), likely due to other changes in TiME altered by the therapy and warrants further analysis to confirm.

Figure 5.

Validation of potential spatial biomarkers with an independent cohort (presurgically treated). A, Validation of imRS by comparing of OS between short-term group, long-term group, and the validation cohort (log-rank test). B, Comparisons of imRS between the short-term group, long-term group, and validation cohort. *, P < 0.05; ****, P < 0.0001 (Wilcoxon rank-sum test). C, Exemplar topological graph of IL10+ myelomonocytes, PD-1+ CD4+ T cells, and GZMB+ CD8+ T cells in the validation cohort. D, Immunodetection of CD8 and CD20 featuring LAs of B cells and CD8+ T cells. E, Validation of immune cell density by comparing B-cell and CD8+ T-cell densities between short-term and long-term groups, and the validation cohort; ns, not significant; **, P < 0.01 (Wilcoxon rank-sum test). F, Hierarchically clustered heatmap of marker expression compositions in LAs. Bar plot show the B-cell and CD8+ T-cell densities for each LA, respectively. G, imRS was computed at per-cell basis and compared between R and NR groups; ****, P < 0.0001 (Wilcoxon rank-sum test). H, Comparison of CD8+ T-cell density in LAs between R and NR groups.

Figure 5.

Validation of potential spatial biomarkers with an independent cohort (presurgically treated). A, Validation of imRS by comparing of OS between short-term group, long-term group, and the validation cohort (log-rank test). B, Comparisons of imRS between the short-term group, long-term group, and validation cohort. *, P < 0.05; ****, P < 0.0001 (Wilcoxon rank-sum test). C, Exemplar topological graph of IL10+ myelomonocytes, PD-1+ CD4+ T cells, and GZMB+ CD8+ T cells in the validation cohort. D, Immunodetection of CD8 and CD20 featuring LAs of B cells and CD8+ T cells. E, Validation of immune cell density by comparing B-cell and CD8+ T-cell densities between short-term and long-term groups, and the validation cohort; ns, not significant; **, P < 0.01 (Wilcoxon rank-sum test). F, Hierarchically clustered heatmap of marker expression compositions in LAs. Bar plot show the B-cell and CD8+ T-cell densities for each LA, respectively. G, imRS was computed at per-cell basis and compared between R and NR groups; ****, P < 0.0001 (Wilcoxon rank-sum test). H, Comparison of CD8+ T-cell density in LAs between R and NR groups.

Close modal

Next, we identified CD8+ T cell–B cell aggregates (n = 31) from the validation cohort and compositions were quantified (Fig. 5D). Similar to the discovery cohort, results revealed no significant difference in terms of B-cell density per cluster; however, CD8+ T-cell density was significantly elevated in the validation cohort compared with short-term survivors (Fig. 5E). Of note, hierarchical clustering demonstrated that the majority of aggregates in the validation cohort were EOMES dominant, which was contrary to the findings from the discovery cohort (Fig. 5F).

In addition, we explored the prognostic values of imRS and CD8+ T-cell density in LAs in predicting presurgical treatment efficacy. We observed that imRS in Rs group (mean = 0.179) is significantly lower compared with NRs (mean = 0.303, Wilcoxon rank-sum test P = 4.098 × 10−7; Fig. 5G), However, we did not observe the same significance of CD8+ T-cell density in Rs (Wilcoxon rank-sum test P = 0.4945; Fig. 5H). This is likely related to the limited sample size because LAs were identified in only 6 of 12 patients. Specifically, 15 LAs were detected in 3 Rs and 16 in 3 NRs, hence limiting the statistical power of the comparison. It is worth noting, the mean CD8+ T-cell density in Rs is 297.4 cells · mm−2, which is much higher compared with the density in NRs of 176.4 cells · mm−2. Univariate Cox regression model indicated that R group with higher density trended toward lower risk (HR = 0.981; CI = 0.94–1.02; Wald test P > 0.05). Correlation analysis revealed the same trend that CD8+ T-cell density positively correlated with OS compared with the discovery cohort (ρ = 0.6, P = 0.35; Supplementary Fig. S5).

Collectively, the validations indicated that the stratification power of imRS was in agreement with both the presurgically treated cohort and long-term survival group, hence suggesting prognostic values. However, the LA-associated features were not recapitulated in the validation cohort and responders to presurgical treatment, demonstrating that these patterns shall be interpreted with caution and further validations on extended dataset are required.

Collectively, these results indicated that pretreatment differences in immune architectures within TiME likely possess prognostic values in patient with pancreatic adenocarcinoma. We observe that PD-1+ CD-4+ T cells appear to present adjacent to IL10+ myelomonocytes in long-term survivors but adjacent to GZMB+ CD8+ T cells in short-term survivors. Though the signaling insights remains unknown, we proposed two hypothesis models that likely reflect the underlying biology (Fig. 6). Previous research provides evidence of the immunosuppressive role of IL10 by downregulating the cytotoxicity of CD8+ T cells (38). Meanwhile, activated CD4+ T cells, signatured by expression of PD-1, are able to secrete IFNγ that can inhibit the production of IL10 (39–41). Taken together, we hypothesized that in the long-term survivors, the immunosuppression of IL10+ myelomonocytes is impeded by the IFNγ production on activated CD4+ T cells; whereas in short-term survivors, IL10+ myelomonocytes may inhibit the cytotoxic CD8+ T cell to promote tumor progression. We also reported that the CD8+ T-cell densities in lymphocytes aggregates are significantly elevated in long-term survivors. Therefore, distinct functional states of immune populations, coupled with their spatial topology, can likely to predict the survival of patients with PDAC.

Figure 6.

Proposed TiME landscape modeling the differences between patients with long-term survival and short-term survival PDAC with regard to immune population activities. In long-term survivors, the immunosuppression activity of IL10+ myelomonocytes were hypothesized to be inhibited by their adjacent PD-1+CD-4+ T cells through IFNγ production; in short-term survivors, IL10+ myelomonocytes were hypothesized to suppress their adjacent GZMB+ CD8+ T cells by direct inhibition on cytotoxicity; in addition, low density of CD8+ T, rather than B cells, in LAs was also observed in short-term survivors. (Schematics created with BioRender.com.)

Figure 6.

Proposed TiME landscape modeling the differences between patients with long-term survival and short-term survival PDAC with regard to immune population activities. In long-term survivors, the immunosuppression activity of IL10+ myelomonocytes were hypothesized to be inhibited by their adjacent PD-1+CD-4+ T cells through IFNγ production; in short-term survivors, IL10+ myelomonocytes were hypothesized to suppress their adjacent GZMB+ CD8+ T cells by direct inhibition on cytotoxicity; in addition, low density of CD8+ T, rather than B cells, in LAs was also observed in short-term survivors. (Schematics created with BioRender.com.)

Close modal

In this study, we investigated features of immune contexture and the spatial landscape of archival PDAC specimens from 45 treatment-naïve PDAC surgical resections, using an mIHC pipeline followed by quantitative spatial characterizations using a computational image processing workflow. The imaging approach and computational pipeline enable simultaneous profiling of multiple leukocyte populations and quantitative assessment of their spatial architectures to identify potential prognostic biomarkers.

Using sequential mIHC on 45 archival pathologic formalin-fixed paraffin-embedded (FFPE) samples with a panel of 27 antibodies enabled identification of cell lineages and their functional states from at least three ROIs per patient and ensuring tumor-immune heterogeneity was captured. The image processing workflow was then applied to extract a single-cell database of 4,026,079 cells featuring CD4+ T cells, CD8+ T cells, CD20+ B cells, myelomonocytes, and neoplastic tumor cells. Their functional states were also determined, reflecting expression of PD-L1, PD-1, IL10, GZMB, Ki67, EOMES, and ICOS. The methodology described in this study represents a multi-scale analysis of the tumor TiME ecosystem spanning from single-cell level properties to spatial clustering patterns. In addition to the PDAC specimens evaluated herein, elements of this platform have been used for biomarker discovery in the context of triple-negative breast cancer, muscle-invasive bladder cancer, and hepatocellular carcinoma, and revealed efficacy in predicting response to cancer treatment, thus providing a general digital pathology framework for deciphering complex spatial biology (22, 23, 32). Given that FFPE specimens are ubiquitously available in laboratories that conduct diagnostic clinical tasks, the merit of the methodology also lies in the minimal materials required for analysis, enabling broad applicability using archival tissue samples. Considering the current advancements in single-cell multiplex proteomics, the framework can also be easily adapted to high-dimensional imaging systems for detailed immune phenotyping that further accelerates biomarker studies with increased profiling bandwidth.

The major goal of this study was to discern subtle biological differences between patients with poor and improved survival for the establishment of prognostic biomarkers. Although PDAC is an aggressive cancer with poor survival, we split the patients into long-term and short-term groups with near equivalent sizes (42). It is important to note that long-term survivors are significant only at statistical level and do not necessarily possess clinical meaning since the median OS is 832 days. This limitation in searching for predictive biomarkers could apply to PDAC generally. Nevertheless, we reasoned that the distinction in OS between long- and short-term survivors was driven by underlying biological behaviors, considering the strong statistical difference.

Here, we focused on the immune populations within the PDAC TiME. To start, we quantified the first-order properties of each immune phenotype. Specifically, we observed that short-term and long-term survivors tended to have similar compositions of immune populations. Although the immune infiltration levels varied across tissue regions, highly immune-infiltrated regions existed in TiME from both long-term and short-term groups. By computing spatial Shannon entropy, we observed that the TiME from both groups were equally heterogeneous. Taken together, these results support the notion that the high intertumoral and intratumoral heterogeneity is a hallmark of PDAC independently of survival status (43). Previous studies reported that immune components within TiME were highly coordinated to orchestrate antitumor immunity, thus we further sought to study the TiME as an ecosystem (25). Using spatial G(r) function, we identified two pairs of immune phenotypes that exhibited spatial clustering and associated with survival: CD4+ T cell–myelomonocyte and CD8+ T cell–B cell. CD4+ T cell and myelomonocyte pairs were a two-way clustering, that is, cells accumulated around cells of different phenotypes. To interrogate the clustering pattern, we developed enrichment scores and discovered that IL10-expressing myelomonocytes and PD-1–expressing CD4+ T cells were key determinants of the spatial dependence. In addition, we proposed an imRS and found that long-term survivors associated with decreased distance between IL10+ myelomonocytes and PD-1+ CD4+ T cells; whereas decreased distance between IL10+ myelomonocytes and GZMB+ CD8+ T cells instead associated with short-term survivors. imRS revealed the balance of IL10+ myelomonocytes between high and low risk of immunosuppression for patient prognosis. Previous research provides evidence of the immunosuppressive role of IL10 by directly and indirectly impacting cytotoxicity of CD8+ T cells (38). Moreover, activated CD4+ T cells, indicated by expression of PD-1, can secrete IFNγ and thereby inhibit neighboring production of IL10 (39–41). Taken together, we hypothesized that in the long-term survivors, immunosuppression of IL10+ myelomonocytes is impeded by the IFNγ production from activated CD4+ T cells; whereas in short-term survivors, IL10+ myelomonocytes in turn inhibit the cytotoxic properties of CD8+ T cells and thereby foster tumor progression. It is worth mentioning that some studies report contradictory findings: Wang and colleagues found that instead of being immunosuppressive, IL10 could enhance antitumor immunity by hampering suppressive CD4+ T cells, thus indicating an opposite signaling axis as compared to the aforementioned hypothesis (44). Therefore, this finding should be interpreted with caution, as we were not able to also evaluate FOXP3 expression to differentiate suppressive T-regulatory CD4+ T cells from Th cells. Second, we examined the one-way clustering pattern of CD8+ T cell–B cell, also termed as LAs in this study. We found that the high CD8+ T-cell density in LA was associated with long-term survival. Similarly, Gunderson and colleagues reported a distinct type of T- and B-cell aggregates, namely early-stage tertiary lymphoid structure, featuring high CD8+ T-cell infiltration that associated with improved survival in PDAC (45).

There are important limitations in our study. First, the discovery dataset (treatment-naïve) is relatively small. However, it is noteworthy that we sampled multiple regions from each patient, sampling on average 89,469 cells per patient. The approach also preserved the intratumoral heterogeneity, thus entailing a generalized platform less sensitive to sampling bias. In addition, the biology of IL10 in PDAC is still poorly understood and limits us to corroborating the connection of PD-1+ CD4+ T cell – IL10+ myelomonocyte – GZMB+ CD8+ T-cell signaling axis to overall patient survival. Although we did validate our findings with an independent cohort, the clinical parameters of the validation cohort were not identical to the discovery cohort, thus the hypothesis requires further validation on larger external datasets.

In conclusion, we developed a multiscale quantitative assay combining sequential IHC imaging and spatial analysis to study the TiME in patients with PDAC in search of potential prognostic biomarkers. While the platform still requires intensive validation on external cohort, it might complement current standard pathology staging and serve as a research tool for quantitative immuno-oncology domain. The proposed platform also reveals a broad applicability and represents a novel application in the field of translational medicine, as well as prospects in initialization and parameterization of computational models (46–52).

C.B. Betts starting working at a new position at Akoya Biosciences during the article review process. E.M. Jaffee reports grants from Lustgarten Foundation and grants and other support from BMS during the conduct of the study; personal fees from Abmeta, Genocea, Achilles, DragonFly, Carta, and Nextcure; other support from Parker Institute; and grants and other support from Genentech outside the submitted work. L.M. Coussens reports nonfinancial support and other support from Cell Signaling Technology, Syndax Pharmaceuticals Inc, and Hibercell, Inc.; nonfinancial support from ZellBio, Inc.; other support from Prospect Creek Foundation, Lustgarten Foundation for Pancreatic Cancer Research, Susan G. Komen Foundation, National Foundation for Cancer Research, Carisma Therapeutics Inc., CytomX Therapeutics, Inc., Kineta Inc., Alkermes, Inc., PDX Pharmaceuticals, Inc., AstraZeneca Partner of Choice Network, Genenta Sciences, Pio Therapeutics Pty Ltd., (P30) Koch Institute for Integrated Cancer Research, Massachusetts Institute of Technology, Bloomberg-Kimmel Institute for Cancer Immunotherapy, Sidney Kimmel Comprehensive Cancer Center at Johns Hopkins, advisory board (honorarium), (P30) Dana-Farber/Harvard Cancer Center, (P30) University of California, San Diego Moores Cancer Center, (P30) The Jackson Laboratory Cancer Center, (P01) Columbia University Medical Center, Prostate P01, (P50) MDACC GI SPORE, Cancer Research Institute, Lustgarten Foundation for Pancreatic Cancer Research, Therapeutics Working Group, American Association for Cancer Research (AACR), NIH/NCI-Frederick National Laboratory Advisory Committee, AACR senior editor for Cancer Immunology Research, AACR scientific editor for Cancer Discovery, Editorial Board member for Cancer Cell, Acerta Pharma, LLC, Innate Pharma, Pharmacyclics, Inc.; Steering Committee for PCYC-1137-CA (NCT02436668), Verseau Therapeutics, Inc., Zymeworks, Inc., Starr Cancer Consortium, and The V Foundation for Cancer Research; and grants from Susan G Komen Foundation, Komen Scholar, AbbVie Inc., Shasqi, Inc., and Cell Signaling Technology outside the submitted work. A.S. Popel reports grants from NIH during the conduct of the study; grants from AstraZeneca, Boehringer Ingelheim, and CytomX Therapeutics; and personal fees from AsclepiX Therapeutics and CytomX Therapeutics outside the submitted work. No disclosures were reported by the other authors.

H. Mi: Conceptualization, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. S. Sivagnanam: Resources, data curation, software, investigation, visualization, methodology, writing–review and editing. C.B. Betts: Writing–review and editing. S.M. Liudahl: Writing–review and editing. E.M. Jaffee: Writing–review and editing. L.M. Coussens: Resources, supervision, funding acquisition, project administration, writing–review and editing. A.S. Popel: Resources, supervision, funding acquisition, project administration, writing–review and editing.

This research was supported by NIH grants U01CA212007 and R01CA138264 to A.S. Popel. L.M. Coussens acknowledges funding from the NIH (1U01CA224012, U2CCA233280), the Knight Cancer Institute, and the OHSU-Brenden-Colson Center for Pancreatic Care. Parts of the schematic figures were produced from BioRender (https://biorender.com).

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).

1.
Rahib
L
,
Smith
BD
,
Aizenberg
R
,
Rosenzweig
AB
,
Fleshman
JM
,
Matrisian
LM
.
Projecting cancer incidence and deaths to 2030: the unexpected burden of thyroid, liver, and pancreas cancers in the United States
.
Cancer Res
2014
;
74
:
2913
21
.
2.
Kleeff
J
,
Korc
M
,
Apte
M
,
La Vecchia
C
,
Johnson
CD
,
Biankin
AV
, et al
.
Pancreatic cancer
.
Nat Rev Dis Primers
2016
;
2
:
16022
.
3.
Mizrahi
JD
,
Surana
R
,
Valle
JW
,
Shroff
RT
.
Pancreatic cancer
.
Lancet
2020
;
395
:
2008
20
.
4.
Wolfgang
CL
,
Herman
JM
,
Laheru
DA
,
Klein
AP
,
Erdek
MA
,
Fishman
EK
, et al
.
Recent progress in pancreatic cancer
.
CA Cancer J Clin
2013
;
63
:
318
48
.
5.
Neoptolemos
JP
.
Adjuvant treatment of pancreatic cancer
.
Eur J Cancer
2011
;
47
:
S378
80
.
6.
Jain
T
,
Dudeja
V
.
The war against pancreatic cancer in 2020—advances on all fronts
.
Nat Rev Gastroenterol Hepatol
2021
;
18
:
99
100
.
7.
Ostios-Garcia
L
,
Villamayor
J
,
Garcia-Lorenzo
E
,
Vinal
D
,
Feliu
J
.
Understanding the immune response and the current landscape of immunotherapy in pancreatic cancer
.
World J Gastroenterol
2021
;
27
:
6775
93
.
8.
Wu
AA
,
Jaffee
E
,
Lee
V
.
Current status of immunotherapies for treating pancreatic cancer
.
Curr Oncol Rep
2019
;
21
:
60
.
9.
Hosein
AN
,
Dougan
SK
,
Aguirre
AJ
,
Maitra
A
.
Translational advances in pancreatic ductal adenocarcinoma therapy
.
Nat Cancer
2022
;
3
:
272
86
.
10.
Lutz
ER
,
Wu
AA
,
Bigelow
E
,
Sharma
R
,
Mo
G
,
Soares
K
, et al
.
Immunotherapy converts nonimmunogenic pancreatic tumors into immunogenic foci of immune regulation
.
Cancer Immunol Res
2014
;
2
:
616
31
.
11.
Bockorny
B
,
Semenisty
V
,
Macarulla
T
,
Borazanci
E
,
Wolpin
BM
,
Stemmer
SM
, et al
.
BL-8040, a CXCR4 antagonist, in combination with pembrolizumab and chemotherapy for pancreatic cancer: the COMBAT trial
.
Nat Med
2020
;
26
:
878
85
.
12.
Karamitopoulou
E
.
Tumour microenvironment of pancreatic cancer: immune landscape is dictated by molecular and histopathological features
.
Br J Cancer
2019
;
121
:
5
14
.
13.
Siret
C
,
Collignon
A
,
Silvy
F
,
Robert
S
,
Cheyrol
T
,
André
P
, et al
.
Deciphering the crosstalk between myeloid-derived suppressor cells and regulatory T cells in pancreatic ductal adenocarcinoma
.
Front Immunol
2020
;
10
:
3070
.
14.
Castino
GF
,
Cortese
N
,
Capretti
G
,
Serio
S
,
Di Caro
G
,
Mineri
R
, et al
.
Spatial distribution of B cells predicts prognosis in human pancreatic adenocarcinoma
.
Oncoimmunology
2016
;
5
:
e1085147
.
15.
Carstens
JL
,
De Sampaio
PC
,
Yang
D
,
Barua
S
,
Wang
H
,
Rao
A
, et al
.
Spatial computation of intratumoral T cells correlates with survival of patients with pancreatic cancer
.
Nat Commun
2017
;
8
:
15095
.
16.
Wang
Y
,
Liang
Y
,
Xu
H
,
Zhang
X
,
Mao
T
,
Cui
J
, et al
.
Single-cell analysis of pancreatic ductal adenocarcinoma identifies a novel fibroblast subtype associated with poor prognosis but better immunotherapy response
.
Cell Discov
2021
;
7
:
36
.
17.
Väyrynen
SA
,
Zhang
J
,
Yuan
C
,
Väyrynen
JP
,
Costa
AD
,
Williams
H
, et al
.
Composition, spatial characteristics, and prognostic significance of myeloid cell infiltration in pancreatic cancer
.
Clin Cancer Res
2021
;
27
:
1069
81
.
18.
Pereira
BA
,
Vennin
C
,
Papanicolaou
M
,
Chambers
CR
,
Herrmann
D
,
Morton
JP
, et al
.
CAF subpopulations: a new reservoir of stromal targets in pancreatic cancer
.
Trends Cancer
2019
;
5
:
724
41
.
19.
Schürch
CM
,
Bhate
SS
,
Barlow
GL
,
Phillips
DJ
,
Noti
L
,
Zlobec
I
, et al
.
Coordinated cellular neighborhoods orchestrate antitumoral immunity at the colorectal cancer invasive front
.
Cell
2020
;
182
:
1341
59
.
20.
Jackson
HW
,
Fischer
JR
,
Zanotelli
VR
,
Ali
HR
,
Mechera
R
,
Soysal
SD
, et al
.
The single-cell pathology landscape of breast cancer
.
Nature
2020
;
578
:
615
20
.
21.
Gong
C
,
Anders
RA
,
Zhu
Q
,
Taube
JM
,
Green
B
,
Cheng
W
, et al
.
Quantitative characterization of CD8+ T cell clustering and spatial heterogeneity in solid tumors
.
Front Oncol
2019
;
8
:
649
.
22.
Mi
H
,
Gong
C
,
Sulam
J
,
Fertig
EJ
,
Szalay
AS
,
Jaffee
EM
, et al
.
Digital pathology analysis quantifies spatial heterogeneity of CD3, CD4, CD8, CD20, and FOXP3 immune markers in triple-negative breast cancer
.
Front Physiol
2020
;
11
:
583333
.
23.
Mi
H
,
Bivalacqua
TJ
,
Kates
M
,
Seiler
R
,
Black
PC
,
Popel
AS
, et al
.
Predictive models of response to neoadjuvant chemotherapy in muscle-invasive bladder cancer using nuclear morphology and tissue architecture
.
Cell Rep Med
2021
;
2
:
100382
.
24.
Keren
L
,
Bosse
M
,
Marquez
D
,
Angoshtari
R
,
Jain
S
,
Varma
S
, et al
.
A structured tumor-immune microenvironment in triple negative breast cancer revealed by multiplexed ion beam imaging
.
Cell
2018
;
174
:
1373
87
.
25.
Liudahl
SM
,
Betts
CB
,
Sivagnanam
S
,
Morales-Oyarvide
V
,
Da Silva
A
,
Yuan
C
, et al
.
Leukocyte heterogeneity in pancreatic ductal adenocarcinoma: phenotypic and spatial features associated with clinical outcome
.
Cancer Discov
2021
;
11
:
2014
31
.
26.
Giraldo
NA
,
Nguyen
P
,
Engle
EL
,
Kaunitz
GJ
,
Cottrell
TR
,
Berry
S
, et al
.
Multidimensional, quantitative assessment of PD-1/PD-L1 expression in patients with Merkel cell carcinoma and association with response to pembrolizumab
.
J Immunother Cancer
2018
;
6
:
99
.
27.
Berry
S
,
Giraldo
NA
,
Green
BF
,
Cottrell
TR
,
Stein
JE
,
Engle
EL
, et al
.
Analysis of multispectral imaging with the AstroPath platform informs efficacy of PD-1 blockade
.
Science
2021
;
372
:
eaba2609
.
28.
Failmezger
H
,
Muralidhar
S
,
Rullan
A
,
de Andrea
CE
,
Sahai
E
,
Yuan
Y
.
Topological Tumor Graphs: a graph-based spatial model to infer stromal recruitment for immunosuppression in melanoma histology
.
Cancer Res
2020
;
80
:
1199
209
.
29.
Baddeley
A
,
Turner
R
.
Spatstat: an R package for analyzing spatial point patterns
.
J Stat Softw
2005
;
12
:
1
42
.
30.
Verlato
G
,
Marrelli
D
,
Accordini
S
,
Bencivenga
M
,
Di Leo
A
,
Marchet
A
, et al
.
Short-term and long-term risk factors in gastric cancer
.
World J Gastroenterol
2015
;
21
:
6434
43
.
31.
Wang
Y
,
Wang
YG
,
Hu
C
,
Li
M
,
Fan
Y
,
Otter
N
, et al
.
Cell graph neural networks enable the precise prediction of patient survival in gastric cancer
.
NPJ Precis Oncol
2022
;
6
:
45
.
32.
Mi
H
,
Ho
WJ
,
Yarchoan
M
,
Popel
AS
.
Multi-scale spatial analysis of the tumor microenvironment reveals features of cabozantinib and nivolumab efficacy in hepatocellular carcinoma
.
Front Immunol
2022
;
13
:
892250
.
33.
Biswas
SK
,
Mantovani
A
.
Macrophage plasticity and interaction with lymphocyte subsets: cancer as a paradigm
.
Nat Immunol
2010
;
11
:
889
96
.
34.
Rodriguez
AB
,
Engelhard
VH
.
Insights into tumor-associated tertiary lymphoid structures: novel targets for antitumor immunity and cancer immunotherapy: tertiary lymphoid structures in cancer
.
Cancer Immunol Res
2020
;
8
:
1338
45
.
35.
Binnewies
M
,
Roberts
EW
,
Kersten
K
,
Chan
V
,
Fearon
DF
,
Merad
M
, et al
.
Understanding the tumor immune microenvironment (TIME) for effective therapy
.
Nat Med
2018
;
24
:
541
50
.
36.
Phillips
D
,
Matusiak
M
,
Gutierrez
BR
,
Bhate
SS
,
Barlow
GL
,
Jiang
S
, et al
.
Immune cell topography predicts response to PD-1 blockade in cutaneous T cell lymphoma
.
Nat Commun
2021
;
12
:
6726
.
37.
Smith
LK
,
Boukhaled
GM
,
Condotta
SA
,
Mazouz
S
,
Guthmiller
JJ
,
Vijay
R
, et al
.
Interleukin-10 directly inhibits CD8+ T cell function by enhancing N-glycan branching to decrease antigen sensitivity
.
Immunity
2018
;
48
:
299
312
.
38.
Biswas
PS
,
Pedicord
V
,
Ploss
A
,
Menet
E
,
Leiner
I
,
Pamer
EG
.
Pathogen-specific CD8 T cell responses are directly inhibited by IL-10
.
J Immunol
2007
;
179
:
4520
8
.
39.
Hu
X
,
Paik
PK
,
Chen
J
,
Yarilina
A
,
Kockeritz
L
,
Lu
TT
, et al
.
IFN-γ suppresses IL-10 production and synergizes with TLR2 by regulating GSK3 and CREB/AP-1 proteins
.
Immunity
2006
;
24
:
563
74
.
40.
Herrero
C
,
Hu
X
,
Li
WP
,
Samuels
S
,
Sharif
MN
,
Kotenko
S
, et al
.
Reprogramming of IL-10 activity and signaling by IFN-γ
.
J Immunol
2003
;
171
:
5034
41
.
41.
Keir
ME
,
Butte
MJ
,
Freeman
GJ
,
Sharpe
AH
.
PD-1 and its ligands in tolerance and immunity
.
Annu Rev Immunol
2008
;
26
:
677
704
.
42.
Siegel
RL
,
Miller
KD
,
Jemal
A
.
Cancer statistics, 2020
.
CA Cancer J Clin
2020
;
70
:
7
30
.
43.
Wandmacher
AM
,
Mehdorn
A-S
,
Sebens
S
.
The heterogeneity of the tumor microenvironment as essential determinant of development, progression and therapy response of pancreatic cancer
.
Cancers
2021
;
13
:
4932
.
44.
Wang
L
,
Liu
J-Q
,
Talebian
F
,
Liu
Z
,
Yu
L
,
Bai
X-F
.
IL-10 enhances CTL-mediated tumor rejection by inhibiting highly suppressive CD4+ T cells and promoting CTL persistence in a murine model of plasmacytoma
.
Oncoimmunology
2015
;
4
:
e1014232
.
45.
Gunderson
AJ
,
Rajamanickam
V
,
Bui
C
,
Bernard
B
,
Pucilowska
J
,
Ballesteros-Merino
C
, et al
.
Germinal center reactions in tertiary lymphoid structures associate with neoantigen burden, humoral immunity and long-term survivorship in pancreatic cancer
.
Oncoimmunology
2021
;
10
:
1900635
.
46.
Ma
H
,
Wang
H
,
Sové
RJ
,
Wang
J
,
Giragossian
C
,
Popel
AS
.
Combination therapy with T cell engager and PD-L1 blockade enhances the antitumor potency of T cells as predicted by a QSP model
.
J Immunother Cancer
2020
;
8
:
e001141
.
47.
Sové
RJ
,
Jafarnejad
M
,
Zhao
C
,
Wang
H
,
Ma
H
,
Popel
AS
.
QSP-IO: a quantitative systems pharmacology toolbox for mechanistic multiscale modeling for immuno-oncology applications
.
CPT Pharmacometrics Syst Pharmacol
2020
;
9
:
484
97
.
48.
Wang
H
,
Ma
H
,
Sové
RJ
,
Emens
LA
,
Popel
AS
.
Quantitative systems pharmacology model predictions for efficacy of atezolizumab and nab-paclitaxel in triple-negative breast cancer
.
J Immunother Cancer
2021
;
9
:
e002100
.
49.
Gong
C
,
Ruiz-Martinez
A
,
Kimko
H
,
Popel
AS
.
A Spatial quantitative systems pharmacology platform spQSP-IO for simulations of tumor—immune interactions and effects of checkpoint inhibitor immunotherapy
.
Cancers
2021
;
13
:
3751
.
50.
Zhang
S
,
Gong
C
,
Ruiz-Martinez
A
,
Wang
H
,
Davis-Marcisak
E
,
Deshpande
A
, et al
.
Integrating single cell sequencing with a spatial quantitative systems pharmacology model spQSP for personalized prediction of triple-negative breast cancer immunotherapy response
.
Immunoinformatics
2021
;
1
:
100002
.
51.
Ruiz-Martinez
A
,
Gong
C
,
Wang
H
,
Sov´e
RJ
,
Mi
H
,
Kimko
H
, et al
.
Simulations of tumor growth and response to immunotherapy by coupling a spatial agent-based model with a whole-patient quantitative systems pharmacology model
.
PLoS Comput Biol
2022
;
18
:
e1010254
.
52.
Wang
H
,
Zhao
C
,
Santa-Maria
CA
,
Emens
LA
,
Popel
AS
.
Dynamics of tumor-associated macrophages in a quantitative systems pharmacology model of immunotherapy in triple-negative breast cancer
.
iScience
2022
;
25
:
104702
.