Quantifying Intratumoral Heterogeneity and Immunoarchitecture Generated In-Silico by a Spatial Quantitative Systems Pharmacology Model

Abstract

Spatial heterogeneity is a hallmark of cancer. Tumor heterogeneity can vary with time and location. The tumor microenvironment (TME) encompasses various cell types and their interactions that impart response to therapies. Therefore, a quantitative evaluation of tumor heterogeneity is crucial for the development of effective treatments. Different approaches, such as multiregional sequencing, spatial transcriptomics, analysis of autopsy samples, and longitudinal analysis of biopsy samples, can be used to analyze the intratumoral heterogeneity (ITH) and temporal evolution and to reveal the mechanisms of therapeutic response. However, because of the limitations of these data and the uncertainty associated with the time points of sample collection, having a complete understanding of intratumoral heterogeneity role is challenging. Here, we used a hybrid model that integrates a whole-patient compartmental quantitative-systems-pharmacology (QSP) model with a spatial agent-based model (ABM) describing the TME; we applied four spatial metrics to quantify model-simulated intratumoral heterogeneity and classified the TME immunoarchitecture for representative cases of effective and ineffective anti-PD-1 therapy. The four metrics, adopted from computational digital pathology, included mixing score, average neighbor frequency, Shannon's entropy and area under the curve (AUC) of the G-cross function. A fifth non-spatial metric was used to supplement the analysis, which was the ratio of the number of cancer cells to immune cells. These metrics were utilized to classify the TME as "cold", "compartmentalized" and "mixed", which were related to treatment efficacy. The trends in these metrics for effective and ineffective treatments are in qualitative agreement with the clinical literature, indicating that compartmentalized immunoarchitecture is likely to result in more efficacious treatment outcomes.

Keywords: agent-based Model (ABM); computational digital pathology; immune checkpoint inhibitor; immunoarchitecture; intratumoral heterogeneity; quantitative systems pharmacology (QSP).

Figure 1

(a) spQSP platform schematic. Cellular and molecular interactions in different compartments. Top: QSP model: tumor-draining lymph node (LN) compartment, central (blood) compartment, tumor compartment, peripheral compartment. Mature antigen-presenting cells process tumor antigen in the tumor compartment and these cells are transported via lymphatic vessels to the lymph node where they prime naïve cytotoxic T lymphocytes (CTL) and T regulatory cells (Treg). These processes result in clonal expansion of the T cells, and these cells are trafficked through blood circulation and extravasate into the tumor microenvironment. Middle: spatio-temporal ABM module substitutes a portion of the tumor compartment from the QSP model, where different subtypes of cancer cells, including cancer stem-like cells, cancer progenitor cells, cancer senescent cells and different T cells, including CTL effector, CTL cytotoxic and CTL exhausted, are represented spatially at a single-cell resolution. Spatio-temporal distribution of soluble cytokines IL-2 and IFNγ are described by partial differential equations (PDEs). (b) Algorithm flowchart, coupling between the QSP and ABM modules.

Figure 2

(a) spQSP initial condition within a flattened computational box with dimensions of 10 × 10 × 0.2 mm. All the cancer cell types, all the TCL subtypes and Treg are displayed with blue, dark green and red, respectively. The initial cancer cells are assigned to 25 × 25 × 4 voxels (0.5 × 0.5 × 0.16 mm) containing 10% CSC and 90% cancer progenitor cells. There is no immune cell in the computational box at the beginning. To make results visually comparable with multiplex immunohistochemistry and immunofluorescence, 2-D snapshots using a slice of the 3-D simulation are chosen. At each time step, different metrics are calculated on a fixed center-cropped tile from the 2-D slides; (b) Relative change in the number of cancer cells for high dose (HD), low dose (LD) and no treatment (NT) cases within the whole tumor. For each case, the line with different patterns is the mean of 10 replications, and the shaded band displays the standard deviation.

Figure 3

Snapshots of the spQSP simulation at different time points for the HD case. Blue: cancer cells; green: CD8+ T cells; red: regulatory T cells. (a) Whole slide image (WSI) represents the cell distribution on a 2-D slide in the middle of the computational box. As time progresses, the T cells penetrate into the tumor space between the cancer cells and eradicate them; (b) 3-D visualization in a sample cube with dimensions of 200 × 200 × 200 $\mathsf{\mu }$m located in the center of the computation domain; (c) region of interest (ROI) marked in WSI represents the cell distribution on the 2-D center-cropped tile from WSI with dimensions of 1 × 1 mm. As time proceeds, TME immunoarchitecture from cold becomes mixed and then compartmentalized.

Figure 4

Snapshots of the spQSP simulation at different time points for the LD case. Blue: cancer cells; green: CD8+ T cells; red: regulatory T cells. (a) Whole slide image (WSI) represents the cell distribution on a 2-D slide in the middle of the computational box. As time proceeds, the T cells infiltrate into the tumor space between the cancer cells but cannot kill them; (b) 3-D visualization in a sample cube with dimensions of 200 × 200 × 200 $\mathsf{\mu }$m located in the center of the computation domain; (c) region of interest (ROI) marked in WSI represents the cell distribution on the 2-D center-cropped tile from the WSI with dimensions of 1 × 1 mm. As time proceeds, TME immunoarchitecture from cold becomes and remains mixed.

Figure 5

(a) Change in mixing score over time within the ROI. Mixing score for LD remains constant while it decreases for HD; (b) time-averaged value of mixing score in the format of a bar plot. LD has a higher time-averaged mixing score; (c) the ratio of the number of cancer cells to immune cells within the ROI. The reduction in the ratio for HD is more significant; (d) bivariant classification TME based on mixing score and ratio of cancer cells to immune cells (green: compartmentalized; red: mixed; purple: cold). Lower mixing score corresponds to the HD case and compartmentalized immunoarchitecture.

Figure 6

(a) Change in average neighbor frequency over time within the ROI. Reduction in average neighbor frequency for HD is more significant; (b) time-averaged value of average neighbor frequency in the format of a bar plot. LD has higher average neighbor frequency; (c) bivariant classification TME based on average neighbor frequency and ratio of cancer cells to immune cells (green: compartmentalized; red: mixed; purple: cold). Lower average neighbor frequency belongs to the HD case and compartmentalized immunoarchitecture.

Figure 6

(a) Change in average neighbor frequency over time within the ROI. Reduction in average neighbor frequency for HD is more significant; (b) time-averaged value of average neighbor frequency in the format of a bar plot. LD has higher average neighbor frequency; (c) bivariant classification TME based on average neighbor frequency and ratio of cancer cells to immune cells (green: compartmentalized; red: mixed; purple: cold). Lower average neighbor frequency belongs to the HD case and compartmentalized immunoarchitecture.

Figure 7

(a) Change in Shannon’s entropy over time within the ROI. The rise in Shannon’s entropy for HD is more significant; (b) time-averaged value of Shannon’s entropy in the format of a bar plot. HD has higher time-averaged value; (c) bivariant classification of TME based on Shannon’s entropy and ratio of cancer cells to immune cells (green: compartmentalized; red: mixed; purple: cold). Higher Shannon’s entropy corresponds to the HD case and has compartmentalized immunoarchitecture.

Figure 8

(a) Change in area under the curve (AUC) of G-cross function over time within the ROI. AUC of G-cross function decreases for the HD case, while it is not sensitive to the change in LD morphology; (b) time-averaged value of AUC of G-cross function in the format of a bar plot. LD has higher time-averaged value for the metric; (c) bivariant classification TME based on AUC of G-cross function and ratio of cancer cells to immune cells (green: compartmentalized; red: mixed; purple: cold). Lower AUC of G-cross function corresponds to the HD case and compartmentalized immunoarchitecture.

Figure 9

Snapshots of the spQSP simulation at different time points for an HD case with higher vascular density to extravasate immune cells (HDV). Blue: cancer cells; green: CD8+ T cells; red: regulatory T cells. (a) Whole slide image (WSI) represents the cell distribution on a 2-D slide in the middle of the computational box. As time proceeds, the larger number of T cells penetrates the space between cancer cells and kills them; (b) 3-D visualization in a sample cube with dimensions 200 × 200 × 200 $\mathsf{\mu }$m located in the center of the computation domain; (c) region of interest (ROI) marked in WSI represents the cell distribution on the 2-D center-cropped tile from WSI with dimensions of 1 × 1 mm. As time proceeds, TME immunoarchitecture changes from cold to become more mixed and then compartmentalized.

Figure 10

(a) Relative change in the number of cancer cells for HD, LD, NT and HDV with higher tumor vascular density (HDV) cases within the whole tumor. For each case, the solid line with different patterns is the mean of 10 replications, and the band displays the standard deviation; (b) the ratio of the number of cancer cells to immune cells within the ROI for HDV. The ratio varies between 2.6 and 0.2 over time. Bivariant classification of TME based on different spatial metrics and ratio of cancer cells to immune cells for HDV (green: compartmentalized; red: mixed; purple: cold). The following correspond to the compartmentalized immunoarchitecture: (c) lower mixing score; (d) lower average neighbor frequency; (e) higher Shannon’s entropy; (f) lower AUC of G-cross function.

Figure 11

(a) Relative change in the total number of cancer cells for 5 responder (R) cases and 5 non-responder (NR) cases. For each case, the solid line with different patterns is the mean of 10 replications, and the band displays the standard deviation (SD). Among R cases, the final number of cancer cells is smaller in R1, R3 and R5 with respect to R2 and R4. NR3, NR5, NR2, NR4 and NR1 have the largest final number of cancer cells among NR cases. (b) Cell distribution in the center-cropped ROI for R cases before treatment (BT). Blue: cancer cells; green: CD8+ T cells; red: regulatory T cells. In R1, R3 and R5, the ratio of cancer cells to immune cells is smaller. (c) Cell distribution in the center-cropped ROI for cases after treatment (AT). R1, R3 and R5 demonstrate more obvious compartmentalized patterns because they show better response to the therapy. (d) Cell distribution in the center-cropped ROI for NR cases before treatment (BT). NR1 and NR4 have more mixed structure before treatment. (e) Cell distribution in the center-cropped ROI for NR cases after treatment (AT). The tumor-immune patterns do not change significantly for NR cases.

Figure 12

The change percentage in different spatial metrics for R and NR cases versus reduction percent of the ratio of cancer cells to immune cells. M1: mixing score; M2: average neighbor frequency; M3: AUC of the G-cross function; M4: Shannon’s entropy. Green: compartmentalized; red: mixed; purple: cold. A change of less than 1% happens when the TME remains cold during treatment. The change between 1% and 10% is linked to the mixed TME, while a change above 10% in the spatial metrics happens when the compartmentalized patterns are observed after treatment. All the R cases have compartmentalized patterns, while all the NR cases have either cold or mixed immunoarchitectural patterns.