The evolution of peripheral immune cell abundance and signaling over time, as well as how these immune cells interact with the tumor, may impact a cancer patient’s response to therapy. By developing an ecological population model, we provide evidence of a dynamic predator–prey-like relationship between circulating immune cell abundance and tumor size in patients that respond to immunotherapy. This relationship is not found either in patients that are nonresponsive to immunotherapy or during chemotherapy. Single-cell RNA sequencing of serial peripheral blood samples from patients shows that the strength of tumor–immune cell interactions is reflected in T cells’ interferon activation and differentiation early in treatment. Thus, circulating immune cell dynamics reflect a tumor’s response to immunotherapy.


The extent to which immune cell phenotypes in the peripheral blood reflect within-tumor immune activity prior to and early in cancer therapy is unclear. To address this question, we studied the population dynamics of tumor and immune cells, and immune phenotypic changes, using clinical tumor and immune cell measurements and single-cell genomic analyses. These samples were serially obtained from a cohort of advanced gastrointestinal cancer patients enrolled in a trial with chemotherapy and immunotherapy. Using an ecological population model, fitted to clinical tumor burden and immune cell abundance data from each patient, we find evidence of a strong tumor-circulating immune cell interaction in responder patients but not in those patients that progress on treatment. Upon initiation of therapy, immune cell abundance increased rapidly in responsive patients, and once the peak level is reached tumor burden decreases, similar to models of predator–prey interactions; these dynamic patterns were absent in nonresponder patients. To interrogate phenotype dynamics of circulating immune cells, we performed single-cell RNA sequencing at serial time points during treatment. These data show that peripheral immune cell phenotypes were linked to the increased strength of patients’ tumor–immune cell interaction, including increased cytotoxic differentiation and strong activation of interferon signaling in peripheral T cells in responder patients. Joint modeling of clinical and genomic data highlights the interactions between tumor and immune cell populations and reveals how variation in patient responsiveness can be explained by differences in peripheral immune cell signaling and differentiation soon after the initiation of immunotherapy.
Immune checkpoint inhibitors can treat a wide range of cancers by targeting immune inhibitory pathways that cancer cells frequently coopt to avoid recognition and to regulate immune proliferation, survival, and effector functions (111). However, clinical response varies substantially, with ∼40% of patients currently experiencing no objective benefit (12, 13). Numerous studies have investigated the role of tumor or tumor-associated immune cell phenotypes in response to immunotherapy (1419). Patient responsiveness has been associated with increased tumor cell mutational load and antigen production (20, 21), and also with greater tumor-associated immune cell infiltration (22), signal production (14), and cross-talk (23). However, the consensus is that these markers are weakly associated with patient response (24). Furthermore, obtaining tumor tissue samples is challenging, especially if a tumor’s immunosuppressive phenotypes evolve over time.
Disease can regulate host immune cell abundance and signaling (2529). Recently, it has been suggested that the frequency of specific peripheral blood immune cells can provide a noninvasive pretreatment indicator of immunotherapy responsiveness, at least in melanoma cancer patients (30). As peripheral blood is easily accessible for serial analysis compared to tumor biopsies, a key question is whether circulating immune cells can serve as a surrogate measurement of a tumor’s interaction with the host immune cells and reflect response to therapy early in the course of treatment. If true, simple blood tests could be developed to guide patient-specific clinical management decisions following the initiation of immunotherapy.
To address these questions, we have measured the strength of patients’ tumor–immune cell interactions, using a data driven ecological mathematical model of the concurrent dynamics of tumor and immune cell abundance. The strength of patients’ tumor–immune cell interactions was then related to immune cell phenotypes experimentally measured using single-cell RNA sequencing (scRNAseq). Fitting the tumor–immune cell interaction model to clinical tumor burden and immune abundance data revealed a consistently increased ability of responders’ immune cells to increase in abundance and indicated that improved tumor cell attack drove decreased tumor burden. The increase in circulating immune cell abundance is concordant with a bolstered antitumor interferon (IFN) signaling state of circulating immune cells and differentiation of T cells to more cytotoxic states, as measured by scRNAseq. This combination of mathematical modeling and genomic analyses suggests that peripheral blood immune cell phenotypes reflect cancer–immune cell interactions and can reliably reveal patient responsiveness to immunotherapy.


Overview of Trial and Patient Cohort.

Patients with advanced gastrointestinal (GI) cancers (colorectal, gastroesophageal, pancreatic, and biliary) were enrolled in a single institution phase I trial (NCT02268825) of a modified FOLFOX6 (mFOLOFX6) chemotherapy regimen followed by a combination of chemotherapy and anti–PD-1 immunotherapy (pembrolizumab) (Fig. 1A). Patient response was assessed according to the RECIST (Response Evaluation Criteria in Solid Tumors) 1.1 guidelines, with responders showing complete/partial response (CR/PR) or stable disease (SD) and nonresponders exhibiting progressive disease (PD) (SI Appendix, Tables S1 and S2). Confirming our classification, 89% of responders survived more than 18 mo after completion of treatment compared to only 26% of nonresponders (Fig. 1B). As reported previously, the tumor’s PD-L1 expression was not strongly predictive of patient response (24). Single-cell phenotypic insights (Fig. 1 C and D) were linked to immune cell function by 1) mathematically modeling patients’ time courses of tumor burden and immune abundance, 2) fitting this model to the clinical data, 3) analyzing temporal changes in the growth rate of the tumor and immune cells, and 4) relating patient-specific model predictions to scRNAseq peripheral immune cell phenotype (Fig. 1E).
Fig. 1.
Overview of the clinical trial treatment strategy, patients’ classification, immune single-cell analysis pipeline, and tumor–immune interaction modeling. (A) Advanced GI patients received mFOLFOX6 chemotherapy at the beginning of the trial for two 14-d cycles. From cycles 3 through 12, they received both mFOLFOX6 and anti–PD-1 immunotherapy. At baseline (cycle 1 = C1), cycle 3 (C3), and cycle 5 (C5) blood was collected and PBMCs were isolated and frozen. (B) Overall survival of responders and nonresponders. (C) PBMC analyses using single-cell RNA sequencing and flow cytometry validation. (D) Flowchart of patient sample selection criteria, showing how patient samples were filtered and analyzed. (E) Mathematical modeling flowchart, depicting how (i) clinical tumor burden data were synthesized and linked to concurrent measurements of PBMC abundance and (ii) how a dynamic model of tumor–immune cell interactions, fitted to this data, allow inference of key biological processes (e.g., the ability of immune cells to kill tumor cells).

Patient-Specific Immune Function Linked to Immunotherapy Success.

Time courses of tumor burden and immune abundance (peripheral blood mononuclear cells, PBMCs) were constructed for each patient (Fig. 1 E, i). Lymphocyte and monocyte abundance was strongly positively correlated with total immune abundance (SI Appendix, Fig. S13), indicating a tight coupling of their population dynamics and motivating the modeling of total immune counts. Tumor burden was measured by combining information from cancer-specific antigen biomarkers and RECIST 1.1 measurements of tumor size, using a Gaussian process latent variable model (SI Appendix). The changes in patients’ tumor burden and immune cell abundance during the trial were described mathematically by a dynamic model of cancer–immune cell interactions (Fig. 1 E, ii). In ecology, interactions between species, where the survival of one depends on attack by another, can be described using predator–prey equations. An adaptation of this ecological theory allowed us to describe the interactions between populations of tumor and immune cells within individual patients. We estimated the strength of this interaction by statistically matching the changing frequency of immune cells and tumor size to model predictions. In the model, the tumor cells (T) are attacked by immune cells (I) and tumor cells induce increase immune cell recruitment. Chemotherapy (C) kills both tumor and immune cells, while PD-1i immunotherapy (P) impacts immune proliferation, recruitment, and cytotoxic tumor activity (Fig. 2A).
Fig. 2.
Patients’ immune cell function in attacking cancer cells and regulating tumor growth measured using a data-driven tumor–immune cell interaction model. (A) Schematic of the mathematical model describing the strength of tumor–immune cell interactions and how their abundances change within a given patient over time. Blue arrows indicate recruitment (triangle tip) and attack interactions (circle tip) between cell types. Green arrows show how immunotherapy influences these interactions and immune population growth. Red arrows indicate chemotherapy effects. Curved arrows indicate intrinsic growth and density dependence within cell types. (B) Statistically fitting the model to clinical data allows an accurate description of observed tumor burden and PBMC abundance across patients and over time. The dashed black line indicates 1:1 model–data correspondence. (C) Histogram showing that responder patients consistently have immune cells with a higher ability to attack cancer cells. (D) Comparison of the speed of growth or decline of the tumor and immune cell populations during the trial, as measured by the RGR of each component between observations. The distinct burst of immune activation in responders (Left) and subsequent tumor decline was negligible in nonresponders (Right). Solid lines show mean trajectories and shaded regions signify model uncertainty intervals. The vertical dashed line indicates the start of immunotherapy and the horizontal gray dashed line shows stable population size). (E) Tumor–immune interaction model predictions of the ability of the immune cells of responders and nonresponders to regulate the growth of the tumor during the trial.
Changes in tumor and immune cell abundance over time were accurately described by statistically fitting the mathematical model to the clinical data, using a Bayesian hierarchical approach (Fig. 2B). This analysis captured the biological differences between tumor and immune populations of responders and nonresponders and the substantial variation between patients within these response categories. Key biological rates that were estimated included 1) how effectively immune cells attack the tumor and 2) the impact of chemotherapy on tumor and immune populations. This identified the consistently improved ability of responder patients’ immune cells to attack the tumor, compared to nonresponders (Fig. 2C).
The timing of most rapid growth/decline of tumor and immune populations was determined by analyzing the populations’ relative growth rates (RGR = speed of population change; positive = growth, negative = decline) (Fig. 2 D and E). The response dynamics were not dependent on the patient’s cancer tissue type. The tumor burden of the responders declined more rapidly during the chemotherapy phase and continued to decline (negative RGR) over time (Fig. 2D). The exception is a time window around day 100 when the immune population was still increasing but the chemotherapy effect was generally decreased; once immune abundance reached a critical level, the tumor began to shrink once again and tumor burden remained substantially below the pretreatment level for the duration of the trial. Interestingly, responders’ PBMCs were also initially less abundant and more sensitive to chemotherapy (more negative RGR) (SI Appendix, Fig. S14). However, their immune cell abundance was boosted following the addition of immunotherapy (Fig. 2D; spike in PBMCs RGR around days 48 to 100). Their immune abundance then stabilized at this level or even increased gradually during the rest of the trial (overall positive RGR).
In contrast to responsive patients, the tumor burden nonresponsive patients declined very little during the preimmunotherapy chemotherapy phase and only marginally in the first weeks of immunotherapy (Fig. 2D). Subsequently, tumor growth accelerated, and the tumor burden returned to the pretreatment level within just 80 to 150 d. Further, nonresponders exhibited a continual decline in immune cell number (negative RGR over most of the trial) and did not experience the immunotherapy-induced boost in immune population growth following the addition of immunotherapy or benefit from immunotherapy. Model analysis showed that prior to immunotherapy the responders’ immune populations less effectively regulated tumor growth (Fig. 2E). However, after immunotherapy induced the growth spike in the responders’ immune population, they became more effective at regulating tumor growth. In contrast, the ability of nonresponders’ immune cells to regulate tumor growth declined continually during the trial and very little benefit of immunotherapy was detected.

Immune Cell Populations Identified Using scRNAseq Profiles.

To understand how phenotype changes of circulating immune cells related to the population dynamics and cell interactions (detailed above), we analyzed phenotypes of PBMCs isolated at three time points during the trial (Fig. 1 A and C). Samples at cycle 1 (C1) provide the baseline before treatment, cycle 3 (C3) reflects treatment with only mFOLFOX6 chemotherapy, and cycle 5 (C5) reflects treatment with both chemotherapy and anti–PD-1 immunotherapy. A total of 13 patients (responder n = 7, nonresponder n = 6) were analyzed by scRNAseq (Fig. 1C). The transcriptional profile of 70,781 immune cells was obtained, revealing a diverse set of 35 cell types. All major PBMC lineages were identified using canonical gene expression markers and analysis of a uniform manifold approximation and projection (UMAP) (Fig. 3 and SI Appendix, Figs. S1–S3 and Table S3).
Fig. 3.
Validated classification of immune cell types, T cells, and monocyte subtypes and identification of the major phenotypic variation within these populations. (A) UMAP of the scRNAseq data of all patient’s PBMCs across analyzed time points. Major PBMC types are labeled (RBC, red blood cells; pDC, plasmacytoid dendritic cells). (B) The agreement between our predicted clusters and public classifications of cell types annotated in two published datasets. (Top) Machine learning prediction: the distribution of immune cells in public datasets predicted to our annotation clusters by Random Forest learner using our predicted clusters as a training set. (Bottom) Shared marker genes: the number of shared genes between public datasets and our predicted clusters (SI Appendix). NKT, natural killer T cells; DCs, dendritic cells). (C) UMAP identification of CD4+ and CD8+ T cell subclusters (TFH, follicular helper) and monocyte subtypes. (D) UMAP representing phenotypic gradients of CD4+ differentiation (top of left subplot: lowest score at right and highest to the left), CD8+ cytotoxic differentiation (bottom of left subplot: lowest score toward the top right and highest at the bottom), and monocyte IFN activation.
The cell-type annotations were validated by comparing our transcriptional profiles and corresponding annotations with published studies of PBMCs (31) and tumor-infiltrating immune cells (32). We found that 96.5% of T cells from the PBMC database and 94.1% of T cells from the tumor-infiltrating dataset were correctly predicted using a machine-learning classifier trained using our annotations (Fig. 3B and SI Appendix, Fig. S3). A similarly high agreement was found between our annotations and published annotations when examining cell-type-specific marker genes and comparing the cell-type connections (Fig. 3B and SI Appendix, Figs. S3 and S4). As a final validation, we profiled eight patients (six responders, two nonresponders) with both scRNAseq and flow cytometry (SI Appendix, Fig. S5). An approximate 1:1 correspondence was found between the abundance of immune cell types obtained using each method (SI Appendix, Fig. S6). Immune cell numbers were quantified in two ways: 1) The frequency of cells refers to the percentage of cells in a sample and 2) the abundance refers to the measured number of cells per unit of peripheral blood.

Signaling Activation in Responders’ T Cells upon Initiation of Immunotherapy.

Signaling dynamics upon initiation of immunotherapy were examined through single-cell pathway activity analysis, using single-sample gene set enrichment analysis (ssGSEA) scores (33) of C2-level and hallmark pathway signatures (34, 35). Pathway differences before therapy, during chemotherapy, and during the early immunotherapy phase of the trial were identified using a random effects linear modeling framework (Fig. 4). This approach partitioned the effects of chemotherapy and immunotherapy on pathway activity while accounting for individual variability in expression. The statistical significance of P values was corrected using Holm’s conservative multiple comparison correction procedure.
Fig. 4.
Pathway signaling activation of multiple immune cell types in responders but not nonresponders following initiation of immunotherapy. (A) The number of molecular pathways impacted by chemotherapy and PD-1 immunotherapy and whether PD-1 immunotherapy effects are specific to responders (black bars) or common across patients. The “chemotherapy all patients” panel shows the numbers pathways changing expression between time C1 and C3 in different cell types. The “immunotherapy all patients” panel shows the numbers of pathways showing trends in expression between C3 and C5 which are common to responders and nonresponders. Finally, the “immunotherapy responders” panel shows the numbers of pathways with trends in expression that are unique to responder patients. Pathways with very differing trends in responders and nonresponders are exemplified on the right side. NK, natural killer. (B) IFN and inflammatory signaling of CD4+ and CD8+ T cells is up-regulated in responders more than nonresponders. GSEA pathway categories reflect the most enriched types of pathways for each cell type. Individual GSEA pathways exhibiting differential trends in expression between responders and nonresponders are shown (dashed lines). Overall trends of pathways within each cellular process (solid lines) and variation (shaded regions) are overlain. (C) Heat map of changes in gene expression of responder and nonresponder CD4+ and CD8+ T cells over time. IFN, cell death, NF-κB, MHC I and II, and migration signature genes are displayed as the proportion of maximum level of each gene. Genes not detected in a cell type are shaded gray. (D) Differences in inflammatory signaling, differentiation, and growth-factor production between the monocytes of responders and nonresponders showing overall trends of pathways within each cellular process (solid lines) and variation (shaded regions). Trends of pathways exhibiting differential expression patterns in responders and nonresponders are indicated by dashed lines. (E) Heat map of changes in gene expression of responder and nonresponder monocytes over time. IFN, cell death, NF-κB, TNF-α, growth-factor production, and migration signature genes are displayed as the proportion of maximum level of each gene. Statistical significance of differences between responders and nonresponders was determined for each gene and corrected for multiple comparisons. C1 = cycle 1: baseline; C3 = cycle 3: chemotherapy mFOLOFX6 regimen; C5 = cycle 5: chemotherapy + anti–PD-1 immunotherapy. One cycle = 14 d.
Overall, immune cell gene expression was not greatly altered during chemotherapy treatment (Fig. 4 A, Left). In contrast, after the start of anti–PD-1 treatment there were a subset of pathway changes common to both responder and nonresponder’s monocytes and T cells (Fig. 4 A, Middle). Further, a majority of signaling changes were identified that were specific to responders (Fig. 4 A, Right and SI Appendix, Table S3). For each immune cell type, the most significantly altered GSEA pathways were classified into categories reflecting major biological processes.
Strikingly, IFN signaling pathway activity was significantly up-regulated in CD4+ and CD8+ T cells of responder patients following the initiation of anti–PD-1 treatment (C3 through C5) (CD4+: t = 19.00, P < 0.001; CD8+: t = 16.00, P < 0.001) (Fig. 4B and SI Appendix, Fig. S7). CD8+ T cells of nonresponders showed a lesser up-regulation of IFN signaling after the start of anti–PD-1 (t = 7.61, P < 0.001), while CD4+ T cells showed no such increase. Upon initiation of immunotherapy, a range of IFN-related genes were up-regulated in the CD8+ and CD4+ T cells of just the responders (Fig. 4C and SI Appendix, Fig. S8). Responders’ CD8+ cells showed greater up-regulation of the IFN-γ gene (P < 0.01) and IFN target genes (IRF1/2/7, STAT1/2, and IFN-stimulated genes; SI Appendix, Table S4). In contrast, nonresponders’ CD4+ and CD8+ T cells had greater up-regulation of IFN-repressing genes (e.g., SOCS1 and SOCS2) (P < 0.05), indicating impaired transduction of IFN signaling upon anti–PD-1 treatment (36). Inflammatory response pathways were also up-regulated in T cells of responders (Fig. 4B), including CD8+ T cells of responders prior to the onset of any treatment (t = 5.14, P < 0.001) and after addition of anti–PD-1 (t = 3.8, P < 0.001). Inflammatory genes induced with anti–PD-1 include major histocompatibility complex (MHC class I/II) sorting and processing genes (e.g., CD74, HLA-A/B/C, and PSM) as well as nuclear factor κB (NF-κB) pathway genes (NFKB1, IKBKB, and MYD88) in responders’ CD8+ and CD4+ T cells (Fig. 4C and SI Appendix, Table S4). The NF-κB activation of responders’ T cells may suggest a shift to a prosurvival state. Overall, this shows the activation of these peripheral cells and the increased signal transduction in responders.

Patients Responsive to Therapy Exhibit Changes in Monocyte Signaling during Treatment.

Monocytes also exhibited different phenotypes in responders versus nonresponders but with distinct signaling changes from those of T cells. Before treatment (C1), monocytes from responders had significantly higher activation of three pathways representing related but distinct measures of monocyte developmental states: growth factor production (t = 9.2, P < 0.001), inflammation (t = 6.1, P < 0.001), and differentiation (t = 6.3, P < 0.001) (Fig. 4D). While chemotherapy decreased each of these pathway scores in both responders and nonresponders, patients responsive to anti–PD-1 treatment exhibited a significant reduction in all three pathways after anti–PD-1 treatment (P < 0.001 for each pathway), while nonresponders showed a significant increase (P < 0.001 for each pathway). During immunotherapy, responders and nonresponders’ monocytes showed specific gene dysregulation of growth factor, IFN, tumor necrosis factor (TNF), NF-κB, and MHC genes (Fig. 4E and SI Appendix, Fig. S9). In addition, genes promoting the migration and recruitment of other immune cells types were initially up-regulated in responders’ monocytes (CXCR4, CCR, and CCL family members) (37) (SI Appendix, Table S4). Overall, monocytes showed pretreatment differences in signaling and divergent developmental trajectories in responders versus nonresponders. Activation of monocytes after the start of anti–PD-1 may reflect responses to the up-regulation of IFN and cytokine gene expression observed in responders’ T cells.

During Therapy, T Cells of Responders Differentiate, While Nonresponder CD8 T Cells Lose Cytotoxicity.

The major phenotypic differences within each immune type were identified, using pseudotime reconstruction of scRNAseq profiles (SI Appendix, Fig. S10). By overlaying the cellular phenotype scores onto a UMAP of the expression profile, we validated that the phenotypes reflect the key sources of transcriptional variation within immune cell types (Fig. 3D). The CD4+ T cell phenotypic gradient captured the continuum of differentiation from naïve to effector helper T cells (Fig. 3 D, Left). Similarly, the CD8+ T cell phenotype gradient captured differentiation from a naive to highly cytotoxic state. In both cases, naive, central memory, and effector T cell subtypes aligned clearly along the continuous phenotype gradient and in the expected order.
We next evaluated the distribution of T cell phenotypes in the peripheral blood of responders and nonresponders and examined how they shifted during the course of therapy (Fig. 5 A and B). Before treatment (C1), responders had a higher frequency of undifferentiated (naive) CD4+ T cells, which may have been symptomatic of the tumor-mediated immune suppression (Fig. 5B). In contrast, nonresponders had more differentiated CD4+ T cells, especially CD4+ effector memory (EM) cells (t = −7.5, P < 0.001) (Fig. 5B). This difference remained following the onset of chemotherapy (C3); however, after immunotherapy (C5) the CD4+ T cells of responders showed a significant shift toward increased differentiation (t = 9.9, P < 0.001) and converged with nonresponders (SI Appendix, Fig. S11A). Interestingly, responders had a higher frequency of cytotoxic differentiated CD8+ T cells than nonresponders, both before and during treatment (Fig. 5B and SI Appendix, Fig. S11B) (F = 16.8, P < 0.001). With the addition of anti–PD-1, responders’ CD8+ T cells became even more cytotoxic (t = 3.9, P < 0.001), while nonresponder’s CD8+ T cells shifted to a less cytotoxic state (t = −4.0, P < 0.001).
Fig. 5.
Peripheral blood immune cell phenotypes linked to patients’ immune cell function and immunotherapy responsiveness. Responsiveness to immunotherapy depends on circulating memory T cell differentiation and monocyte IFN activation prior to therapy. (A) Comparison of CD4+ and CD8+ T cell subtype differentiation scores (all subtypes differ with a Tukey test). EMRA, effector memory CD45RA+; CM, central memory. (B) Frequency of CD4+ and CD8+ T cells with different states of differentiation/cytotoxicity in responders and nonresponders at each treatment time point. (C) Frequency of monocytes with different IFN activation states in responders and nonresponders at each time point. (D and E) The ability of patients’ immune cells to attack cancer cells and also the tumor’s sensitivity to chemotherapy was linked to immune cell signaling and differentiation phenotypes. For each patient, the single-cell variability in immune cell phenotypes is presented as an individual violin densities. The black line indicates the relationship between a patient’s average immune cell phenotype and the strength of immune cell attack/chemotherapy sensitivity. Shaded regions show credible intervals for the predicted range of phenotypes of 95% of the immune cells, given the strength of immune cell attack/chemotherapy sensitivity.

Monocytes of Responders Were Activated after the Start of Anti–PD-1 Therapy and the Frequency of Classical Monocytes Was Associated with Response.

Within monocytes, the expression of IFN response genes was the major axis of phenotypic variation (Fig. 3 C and D, Right). Monocytes with high IFN response scores (including dendritic cells) had up-regulation of IFN stimulation genes (e.g., IFIT1/3, PSME2, and ISG15) and higher MHC class II expression (e.g., HLA.DPA1, HLA.DPB1, and HLA.DMA). In contrast, cells with low IFN scores had up-regulation of proliferation (e.g., FOS, JUN, and JUNB), differentiation (e.g., BTG1, RGS2, and DDX17), inflammation (e.g., SELL, S100A12, and CD36), and migration (e.g., VCAN and VIM) genes. After immunotherapy, monocytes with the highest IFN scores became prevalent in responders (t = 15.463, P < 0.001) (Fig. 5C and SI Appendix, Fig. S14D). Responder patients shifted from having the lowest to the highest average level of IFN activation and MHC class II gene expression (SI Appendix, Fig. S12). In contrast, the distribution of IFN response in nonresponder monocytes remained relatively constant across the trial period.

Linking Immune Function and Phenotypes of the Peripheral Blood.

Finally, we linked the patient-specific estimates of immune attack and chemotherapy sensitivity to the single-cell transcriptomic observations of increased immune cell signaling and phenotypic differentiation states in responders (Fig. 5 D and E). Patients whose immune population had a greater ability to attack tumor cells and response to immunotherapy were found to have CD8+ T cells with higher activity of IFN-γ signaling pathways and more differentiated cytotoxic CD8+ T immune cells (Fig. 5D). Finally, patients whose monocytes showed lower activity of IFN-γ response pathways (classical monocyte differentiation score) before treatment had tumor cells that were significantly less sensitive to chemotherapy (Fig. 5E).


Our findings indicate that peripheral blood immune cell phenotypes reflect the strength of tumor–immune interactions before or early in the course of immunotherapy, and these phenotypes are indicative of patient responsiveness. By combining scRNAseq analysis of peripheral immune phenotypes with dynamical models of patient-specific clinical data, we linked peripheral immune cell phenotypes with the strength of patients’ tumor–immune cell interactions. Increased IFN signaling and differentiation of T cells was related to an increased ability of immune cells to attack cancer cells, regulate tumor growth, and drive patient responsiveness to anti–PD-1 therapy. These results provide motivation for studies interrogating the utility of peripheral blood phenotypes as a biomarker of patient responsiveness to therapy.
Although mathematical modeling has provided important insights into cancer–immune cell interactions and cancer immunotherapy, models incorporating patient-specific clinical or phenotypic data had not previously been developed (3845). Previous theoretical models that do not include patient data have described the potential for cancer–immune interactions to act as “predator–prey-like” systems (reviewed in ref. 45). This study is a step forward in that it uses temporal clinical and single-cell immune phenotyping for data driven ecological modeling of patient-specific responses during treatment.
The cancer–immune interaction model predicts that, in general, patients whose tumors have an immunosuppressive phenotype (e.g., expressing high levels of PD-L1) will have a lower immune cell count prior to treatment, as immune activation and proliferation is inhibited. Hence, we expect that patients with a low PBMC abundance should benefit most from anti–PD-1 immune reactivation therapy. In agreement, we observed significantly lower PBMC abundances in responders at the onset of therapy (SI Appendix, Fig. S15A). These patients showed gradually increasing immune counts during therapy, in contrast to declines observed in nonresponders. Model analysis indicated that, at the onset of the trial, the immune cells of responders had a substantially weaker effect of tumor regulation compared to those in nonresponders, primarily due to the low immune cell count (SI Appendix, Fig. S15B). During immunotherapy, the responders’ immune population gradually increased and their tumor regulatory effect increased toward the level of the nonresponders. This leads to the prediction that, unlike chemotherapy, the tumor’s response to immunotherapy will be delayed. This is a general prediction that emerges from predator–prey models. Due to fewer immune cells present and few cancer antigens being presented to initiate further immune response prior to therapy, several rounds of the cancer–immune response cycle are needed for the immune population to rebuild following PD-L1 suppression.
Our model also predicts that chemotherapy acts as a double-edged sword when used as a combination therapy with immunotherapy. It has the positive effect of inducing tumor cell death and promoting immune cell recruitment; however, it also kills immune cell progenitors, reducing the active immune cell abundance. Therefore, too high a chemotherapy dosage may inhibit the effectiveness of immunotherapy, while too low a level may not promote immune reactivation.
The analyses of T cell and monocyte signaling states, before and during therapy, suggest that circulating immune cells rapidly shift phenotypes during the treatment in GI cancer patients. We suggest that this peripheral immune signaling activation is a valuable early marker of patient responsiveness. The IFN surge after initiation of anti–PD-1 therapy, seen only in responders’ T cells and monocytes, indicates that treatment with anti–PD-1 is promoting differentiation and activation of T cells, resulting in antitumor activity, cytokine release, and stimulation of the immune system. In particular, only responders’ CD8+ T cells up-regulate IFN-γ signaling and immune cell activation and antitumor effect (46). Despite PD-1 blockade, nonresponders’ immune cells were not fully activated, indicating that they struggle to detect cancer cells. Possibly, low cancer antigen release reduced activation of antigen-presenting cells and T cells and prevented initiation of an immune response. Additional studies support an interaction of chemotherapy with immunotherapy in some settings (4751). Using our scRNAseq time courses, we also detected that immunotherapy induces a shift to a more differentiated CD4+ T cell state. Long-term chemotherapy may increase the production of PD-1–expressing regulatory CD4+ EM cells, diminishing pembrolizumab availability to tumor-specific CD8+ T cells (SI Appendix, Fig. S16).
Additionally, patients may have been nonresponsive because cancer cells had PD-1–independent resistance mechanisms of immune avoidance. Indeed, we found that nonresponders’ classical monocytes had low MHC II receptor expression, suggesting lower antigen recognition and presentation. They also developed a more immunosuppressed phenotype, with up-regulation of CD86, a ligand of both PD-1 and CTLA-4, and CD28, a costimulatory signal for activation of T cells. Contrastingly, under anti–PD-1 therapy responders’ monocytes showed activation of costimulatory immune function (up-regulated ISG and MHC).
Overall, we find that the abundance, signaling activity, and differentiation state of peripheral immune cells reflect tumor–immune cell interactions and patient response to immunotherapy. The combination of total PBMC abundance and the relative infrequency of differentiated/activated effector T cells likely provides a noninvasive up-front marker of therapeutic responsiveness. Models of tumor–immune cell interactions, which use clinical and phenotype data, allow quantification of the immune system’s effectiveness in regulating tumor growth and demonstrate the potential of using peripheral blood-based models to assess the dynamics of the immune and tumor cell interactions during treatment.

Materials and Methods

Study Design.

Cryopreserved PBMC samples from patients with advanced (stage 3/4) GI cancers were collected from patients in a clinical trial (NCT02268825) and were treated with a modified FOLFOX6 regimen every 2 wk (i.e., one cycle) until disease progression, death, or completion of the study. After 4 wk of mFOLFOX6 (cycle 3), pembrolizumab (200 mg intravenously every 2 wk) was added to mFOLFOX6. Before treatment and then every 2 wk patients’ blood was collected and PBMCs were isolated and cryopreserved. All human biological samples were collected after written informed patient consent and ethics committee approval, following federal and institutional guidelines. The University of Utah Institutional Review Board and the Huntsman Cancer Institute Protocol Review and Data and Safety Monitoring Committee approved and monitored this study.
The primary outcomes of this phase I study were safety and dose-limiting toxicities. Patients were excluded if they had active infection, autoimmune disease, or were on chronic systemic steroids or immunosuppressant’s. Samples from 13 patients (responder n = 7, nonresponder n = 6) were used for scRNAseq analysis at C1, C3, and C5 time points. Samples from eight patients were utilized for both flow cytometry and scRNAseq analysis (responder n = 6, nonresponder n = 2) to validate the consistency of inferences. Single-cell transcriptional profiling provided information for a total of 70,781 cells from 13 patients.
Clinical response was measured by computed tomography scans and assessed according to RECIST 1.1 and immune-related response criteria every 12 wk. Responders were defined as patients with clinical benefit at 24 wk (CR, PR, or SD). Nonresponders included patients with progressive disease (PD, defined as >20% increase in tumor volume or appearance of new metastatic lesions) between 12 and 24 wk after the trial began. Median of previous history of chemotherapy treatment for responders was 101 d and 42 d for nonresponders (SI Appendix, Table S1).

scRNAseq and Annotation.

PBMC samples analyzed using a Chromium 10X Cell Instrument (10X Genomics) (1,200 to 2,000 cells per sample) and sequenced on an Illumina HiSEq. 2500 with 2 × 125 paired-end reads. Raw BCL sequencing files were processed using Cell Ranger Single Cell Software Suite and samples were aligned to hg19 using the STAR aligner (52). Count tables were generated for 70,781 cells and used as input into Seurat v2 (53). No batch effects were found corresponding to time, patient, or cancer (SI Appendix, Fig. S2 BD).
To identify cell types, variable genes (n = 1,000) and nonoverlapping known immune cell marker genes (n = 1,480) were used for principal component analysis (PCA) (5456). The first 25 PCs captured significant variation, based on Seurat’s jackstraw analysis, and were used for graph-based clustering and UMAP visualization (57). Major T cell clusters were identified using CD3D, CD4, and CD8 expression along with 500 T cell-specific variable genes and 273 known T cell markers (56). Differential expression markers for each cluster were generated using MAST (58, 59). Pathway ssGSEA enrichment scores were generated using the R package GSVA 1.30.0 (33). Immune cell annotations were verified using two public datasets (31, 32) (SI Appendix, Figs. S3 and S4) using training and classification to measure similarity of annotation.

Identifying Gene Set Expression Differences between Responders and Nonresponders.

Differences in the gene set expression of immune cell types were examined between responder and nonresponder patients (R). For each immune cell type, we examine the changes in pathway (X) expression over time (T) and with the addition of the anti–PD-1 (P). A random effects model with the following linear predictor (η) and error structure was constructed for each pathway:
Initial differences in gene set expression between immune cells from responders and nonresponders, at the pretreatment time point (C1), were captured by the group-specific intercepts (β0vs.βR). Differential trends in expression over the first five treatment cycles were described by the group-specific slope terms of responders and nonresponders (βTvsβTR). Differential effects of the addition of anti–PD-1 on gene expression, over cycles C3 to C5, were described by the group-specific anti–PD-1 treatment effect terms of responders and nonresponders (βPTvsβPTR).
Background individual variability in gene expression, independent of therapy impacts, was accounted for by allowing the model intercept to vary among patients (ui). Significant differences in 1) initial pathway scores, 2) temporal trend, and 3) anti–PD-1 treatment effects between nonresponders and responders were assessed using likelihood ratio tests. Multiple comparison corrections were made using Holm’s P value correction.

Quantifying Immune Cell Phenotypes.

Major axes of phenotypic variation were identified separately for CD4+/CD8+ T cells and monocytes using affinity-based pseudotime reconstruction of cell states (6062). This allowed the description of continuous spectrums of cellular states, as is produced by differentiation and activation processes (SI Appendix). These phenotypic axes were validated using comparisons to PCA, zinbwave, and UMAP dimension reduction (57, 63). Random effects linear regression was used to test the statistical differences in immune population phenotype distributions between responders and nonresponders, while accounting for patient-specific random effects.

Modeling and Measuring Tumor–Immune Cell Interactions.

Overall measures of tumor burden.

We assessed the strength of tumor–immune cell interactions and the predictability of responsiveness to therapy by fitting a coupled tumor–immune population model to clinical patient data (Dataset S1). For each patient, a time series of tumor burden was first constructed, by combining RECIST 1.1 measurements, from computed tomography scans, with information from tumor burden biomarkers (CA 19-9 and CEA), using a Gaussian process model (64). Gaussian process models probabilistically combine these tumor burden data sources, allowing inference of tumor burden (SI Appendix).

Tumor–immune interaction model.

The dynamics of tumor and immune cell abundance were coupled with the immunotherapy and chemotherapy dosing schedules, using a patient-specific tumor–immune population dynamic model. The ecologically inspired model (Eq. 1) describes the patient-specific changes in tumor (T) and immune cell (I) abundance over time. Over short periods of time, the increase or decrease in tumor and immune cell abundance was measured by the population’s RGR (RGRT for tumor and RGRIfor immune cells). Positive RGR values indicate population growth, while negative values show population decline. The data-driven model decomposed this population growth rate into effects of different concurrent biological processes. Tumor and immune cells interact in two main ways, with tumor cells being attacked by immune cells (α) and also inducing increased immune cell recruitment (λ). Therapeutic dosing impacts the cell populations and the strength of their interactions, with chemotherapy (C) killing both tumor (μT) and immune cells (μI), while PD-1i immunotherapy (P) influences immune proliferation (βr), recruitment (βλ), and cytotoxic tumor-killing activity (βφ). Both tumor and immune cells experience density-dependent population growth (γT and γI), reflecting competition for resources or growth-stimulating molecules. This leads to the equations
We simultaneously fitted this model to all of the patients’ time-course tumor and immune data and accounted for the differing dosages and timings of therapy. To capture interpatient biological differences, patient-specific parameters were assumed to be drawn from a hyperdistribution of parameters, creating a hierarchical model structure. Model parameters were estimated using Bayesian inference in Stan (65).

Linking Immune Phenotypes and Model-Estimated Biological Processes.

Immune cell phenotypes were related to the model estimates of 1) the effectiveness of immune cells at attacking tumor cells and 2) the tumor cell sensitivity to chemotherapy. These biological estimates of immune and chemotherapy function (X) were regressed against the peripheral immune cell phenotypes identified in 1) the GSEA pathway analysis and 2) the pseudotime analysis of the major phenotypic variation within cell types. For each phenotype, the significance of the relationship between single-cell peripheral immune phenotypes (Y) and immune/chemotherapy function (X) was assessed. A patient-specific intercept was added to account for nonindependence of cell phenotypes within a patient. The random effects regression model was simply
The significance of the relationship between peripheral phenotypes and immune/chemotherapy function was assessed using a likelihood ratio test, with the sample size corrected for the nonindependence of data points.

Data and Materials Availability.

Raw scRNAseq data have been deposited in the Gene Expression Omnibus (GEO) database under accession code GSE130157. Tumor and immune abundance clinical time courses are provided in Dataset S1.

Data Availability

Data deposition: Raw single-cell RNA-sequencing data have been deposited in the Gene Expression Omnibus (GEO) database, (accession no. GSE130157).


We thank the anonymous patients used in our study and S. Brady and B. Nelson for helpful feedback. The High-Throughput Genomics Shared Resource was supported by NIH Award P30CA042014. The content is solely the authors’ responsibility and does not necessarily represent the official views of the NIH. J.T.C. was supported by Cancer Prevention Research Institute of Texas Core Facility Support Award RP170668. A.H.B., J.T.C., J.I.G., P.A.C., P.W., X.L., P.J.M., J.A.M., and F.R.A. were supported by National Cancer Institute, NIH Award U54CA209978.

Supporting Information

Appendix (PDF)
Dataset_S01 (CSV)
Dataset_S02 (XLSX)
Dataset_S03 (XLSX)
Dataset_S04 (XLSX)


M. Loos et al., Clinical significance and regulation of the costimulatory molecule B7-H1 in pancreatic cancer. Cancer Lett. 268, 98–109 (2008).
T. Nomi et al., Clinical significance and therapeutic potential of the programmed death-1 ligand/programmed death-1 pathway in human pancreatic cancer. Clin. Cancer Res. 13, 2151–2157 (2007).
Y. Ohigashi et al., Clinical significance of programmed death-1 ligand-1 and programmed death-1 ligand-2 expression in human esophageal cancer. Clin. Cancer Res. 11, 2947–2953 (2005).
E. Oki et al., Protein expression of programmed death 1 ligand 1 and HER2 in gastric carcinoma. Oncology 93, 387–394 (2017).
M. Song et al., PTEN loss increases PD-L1 protein expression and affects the correlation between PD-L1 expression and clinical parameters in colorectal cancer. PLoS One 8, e65821 (2013).
L. Wang et al., Clinical significance of B7-H1 and B7-1 expressions in pancreatic carcinoma. World J. Surg. 34, 1059–1065 (2010).
Y. Ye et al., Interaction of B7-H1 on intrahepatic cholangiocarcinoma cells with PD-1 on tumor-infiltrating T cells as a mechanism of immune evasion. J. Surg. Oncol. 100, 500–504 (2009).
M. Y. Zhang, Y. Y. Yang, X. H. Wang, X. F. Li, [Expression of Bcl-2, PD-L1 and its clinical significance in colorectal cancer]. Sichuan Da Xue Xue Bao Yi Xue Ban 43, 827–829, 859 (2012).
S. M. Ansell et al., PD-1 blockade with nivolumab in relapsed or refractory Hodgkin’s lymphoma. N. Engl. J. Med. 372, 311–319 (2015).
E. B. Garon et al.; KEYNOTE-001 Investigators, Pembrolizumab for the treatment of non-small-cell lung cancer. N. Engl. J. Med. 372, 2018–2028 (2015).
P. T. Nghiem et al., PD-1 blockade with pembrolizumab in advanced merkel-cell carcinoma. N. Engl. J. Med. 374, 2542–2552 (2016).
O. Hamid et al., Safety and tumor responses with lambrolizumab (anti-PD-1) in melanoma. N. Engl. J. Med. 369, 134–144 (2013).
S. L. Topalian et al., Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N. Engl. J. Med. 366, 2443–2454 (2012).
M. Ayers et al., IFN-γ-related mRNA profile predicts clinical response to PD-1 blockade. J. Clin. Invest. 127, 2930–2940 (2017).
T. Powles et al., MPDL3280A (anti-PD-L1) treatment leads to clinical activity in metastatic bladder cancer. Nature 515, 558–562 (2014).
N. J. Llosa et al., The vigorous immune microenvironment of microsatellite instable colon cancer is balanced by multiple counter-inhibitory checkpoints. Cancer Discov. 5, 43–51 (2015).
D. T. Le et al., Mismatch repair deficiency predicts response of solid tumors to PD-1 blockade. Science 357, 409–413 (2017).
J. Grosso et al., Association of tumor PD-L1 expression and immune biomarkers with clinical activity in patients (pts) with advanced solid tumors treated with nivolumab (anti-PD-1; BMS-936558; ONO-4538). JTO 31, 3016 (2013).
D. P. Carbone et al.; CheckMate 026 Investigators, First-line nivolumab in stage IV or recurrent non-small-cell lung cancer. N. Engl. J. Med. 376, 2415–2426 (2017).
N. A. Rizvi et al., Cancer immunology. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer. Science 348, 124–128 (2015).
R. M. Samstein et al., Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat. Genet. 51, 202–206 (2019).
A. C. Huang et al., T-cell invigoration to tumour burden ratio associated with anti-PD-1 response. Nature 545, 60–65 (2017).
C. S. Garris et al., Successful anti-PD-1 cancer immunotherapy requires T-cell-dendritic cell crosstalk involving the cytokines IFN-gamma and IL-12. Immunity 49, 1148–1161.e7 (2018).
T. R. Cottrell, J. M. Taube, PD-L1 and emerging biomarkers in immune checkpoint blockade therapy. Cancer J. 24, 41–46 (2018).
S. B. Coffelt et al., IL-17-producing γδ T cells and neutrophils conspire to promote breast cancer metastasis. Nature 522, 345–348 (2015).
A. R. Abbas, K. Wolslegel, D. Seshasayee, Z. Modrusan, H. F. Clark, Deconvolution of blood microarray data identifies cellular activation patterns in systemic lupus erythematosus. PLoS One 4, e6098 (2009).
J. K. Neuenburg et al., T-cell activation and memory phenotypes in cerebrospinal fluid during HIV infection. J. Acquir. Immune Defic. Syndr. 39, 16–22 (2005).
S. Chen et al., Increased abundance of myeloid-derived suppressor cells and Th17 cells in peripheral blood of newly-diagnosed Parkinson’s disease patients. Neurosci. Lett. 648, 21–25 (2017).
A. R. Whitney et al., Individuality and variation in gene expression patterns in human blood. Proc. Natl. Acad. Sci. U.S.A. 100, 1896–1901 (2003).
C. Krieg et al., High-dimensional single-cell analysis predicts response to anti-PD-1 immunotherapy. Nat. Med. 24, 144–153 (2018).
E. Azizi et al., Single-cell map of diverse immune phenotypes in the breast tumor microenvironment. Cell 174, 1293–1308.e36 (2018).
M. Sade-Feldman et al., Defining T-cell states associated with response to checkpoint immunotherapy in melanoma. Cell 175, 998–1013.e20 (2018).
S. Hänzelmann, R. Castelo, J. Guinney, GSVA: Gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14, 7 (2013).
A. Liberzon et al., The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 1, 417–425 (2015).
A. Liberzon et al., Molecular signatures database (MSigDB) 3.0. Bioinformatics 27, 1739–1740 (2011).
K. Schroder, P. J. Hertzog, T. Ravasi, D. A. Hume, Interferon-gamma: An overview of signals, mechanisms and functions. J. Leukoc. Biol. 75, 163–189 (2004).
C. Shi, E. G. Pamer, Monocyte recruitment during infection and inflammation. Nat. Rev. Immunol. 11, 762–774 (2011).
G. E. Mahlbacher, K. C. Reihmer, H. B. Frieboes, Mathematical modeling of tumor-immune cell interactions. J. Theor. Biol. 469, 47–60 (2019).
J. Ozik et al., High-throughput cancer hypothesis testing with an integrated PhysiCell-EMEWS workflow. BMC Bioinformatics 19 (suppl. 18), 483 (2018).
M. Robertson-Tessi, A. El-Kareh, A. Goriely, A mathematical model of tumor-immune interactions. J. Theor. Biol. 294, 56–73 (2012).
K. J. Mahasa, R. Ouifki, A. Eladdadi, L. Pillis, Mathematical model of tumor-immune surveillance. J. Theor. Biol. 404, 312–330 (2016).
A. Carbo et al., Systems modeling of molecular mechanisms controlling cytokine-driven CD4+ T cell differentiation and phenotype plasticity. PLOS Comput. Biol. 9, e1003027 (2013).
L. G. de Pillis, A. E. Radunskaya, C. L. Wiseman, A validated mathematical model of cell-mediated immune response to tumor growth. Cancer Res. 65, 7950–7958 (2005).
H. C. Wei, J. L. Yu, C. Y. Hsu, Periodically pulsed immunotherapy in a mathematical model of tumor, CD4+ T cells, and antitumor cytokine interactions. Comput. Math. Methods Med. 2017, 2906282 (2017).
R. Eftimie, J. L. Bramson, D. J. Earn, Interactions between the immune system and cancer: A brief review of non-spatial mathematical models. Bull. Math. Biol. 73, 2–32 (2011).
J. L. Reading, S. A. Quezada, Too much of a good thing? Chronic IFN fuels resistance to cancer immunotherapy. Immunity 45, 1181–1183 (2016).
M. Dosset et al., PD-1/PD-L1 pathway: An adaptive immune resistance mechanism to immunogenic chemotherapy in colorectal cancer. OncoImmunology 7, e1433981 (2018).
I. Péguillet et al., High numbers of differentiated effector CD4 T cells are found in patients with cancer and correlate with clinical response after neoadjuvant therapy of breast cancer. Cancer Res. 74, 2204–2216 (2014).
R. Verma et al., Lymphocyte depletion and repopulation after chemotherapy for primary breast cancer. Breast Cancer Res. 18, 10 (2016).
M. Palma et al., T cells in chronic lymphocytic leukemia display dysregulated expression of immune checkpoints and activation markers. Haematologica 102, 562–572 (2017).
P. Yang, J. Ma, X. Yang, W. Li, Peripheral CD4+ naïve/memory ratio is an independent predictor of survival in non-small cell lung cancer. Oncotarget 8, 83650–83659 (2017).
A. Dobin et al., STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).
A. Butler, P. Hoffman, P. Smibert, E. Papalexi, R. Satija, Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36, 411–420 (2018).
I. Tirosh et al., Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science 352, 189–196 (2016).
A. M. Newman et al., Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457 (2015).
M. De Simone et al., Transcriptional landscape of human tissue lymphocytes unveils uniqueness of tumor-infiltrating T regulatory cells. Immunity 45, 1135–1147 (2016).
J. H. L. McInnes, J. Melville, UMAP: Uniform manifold approximation and projection for dimension reduction. arXiv:1802.03426 (9 February 2018).
G. Finak et al., MAST: A flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biol. 16, 278 (2015).
B. Bischl et al., mlr: Machine learning in R. J. Machine Learning Res. 17, 1–5 (2016).
K. Moon et al., Manifold learning-based methods for analyzing single-cell RNA-sequencing data. Curr. Opin. Syst. Biol. 7, 36–46 (2018).
D. van Dijk et al., Recovering gene interactions from single-cell data using data diffusion. Cell 174, 716–729 e27 (2018).
K. Moon et al., Visualizing structure and transitions for biological data exploration. SSRN Electronic J., (2018).
D. Risso, F. Perraudeau, S. Gribkova, S. Dudoit, J. P. Vert, A general and flexible method for signal extraction from single-cell RNA-seq data. Nat. Commun. 9, 284 (2018).
N. Lawrence, Probabilistic non-liner component analysis with Gaussian process latent variable models. JMLR 6, 1783–1816 (2005).
B. Carpenter et al., Stan: A probabilistic programming language. J. Stat. Software 76, 32 (2017).

Information & Authors


Published in

Go to Proceedings of the National Academy of Sciences
Proceedings of the National Academy of Sciences
Vol. 117 | No. 27
July 7, 2020
PubMed: 32571915


Data Availability

Data deposition: Raw single-cell RNA-sequencing data have been deposited in the Gene Expression Omnibus (GEO) database, (accession no. GSE130157).

Submission history

Published online: June 22, 2020
Published in issue: July 7, 2020


  1. cancer
  2. immunotherapy
  3. ecological models
  4. single-cell RNA sequencing
  5. phenotypes


We thank the anonymous patients used in our study and S. Brady and B. Nelson for helpful feedback. The High-Throughput Genomics Shared Resource was supported by NIH Award P30CA042014. The content is solely the authors’ responsibility and does not necessarily represent the official views of the NIH. J.T.C. was supported by Cancer Prevention Research Institute of Texas Core Facility Support Award RP170668. A.H.B., J.T.C., J.I.G., P.A.C., P.W., X.L., P.J.M., J.A.M., and F.R.A. were supported by National Cancer Institute, NIH Award U54CA209978.


This article is a PNAS Direct Submission.



Department of Medical Oncology & Therapeutics Research, City of Hope National Medical Center, Duarte, CA 91010;
Department of Mathematics, University of Utah, Salt Lake City, UT 84112;
Department of Medical Oncology & Therapeutics Research, City of Hope National Medical Center, Duarte, CA 91010;
Lance T. Pflieger1
Department of Medical Oncology & Therapeutics Research, City of Hope National Medical Center, Duarte, CA 91010;
College of Pharmacy, University of Minnesota, Duluth, MN 55812;
Xuan Liu
Department of Integrative Biology and Pharmacology, School of Medicine, School of Biomedical Informatics, The University of Texas Health Science Center at Houston, Houston, TX 77030;
Department of Medical Oncology & Therapeutics Research, City of Hope National Medical Center, Duarte, CA 91010;
Neena A. Leggett
Department of Integrative Biology and Pharmacology, School of Medicine, School of Biomedical Informatics, The University of Texas Health Science Center at Houston, Houston, TX 77030;
Jasmine A. McQuerry
Department of Medical Oncology & Therapeutics Research, City of Hope National Medical Center, Duarte, CA 91010;
Department of Oncological Sciences, School of Medicine, University of Utah, Salt Lake City, UT 84112;
Department of Pharmacology and Toxicology, College of Pharmacy, University of Utah, Salt Lake City, UT 84112;
Department of Pharmacology and Toxicology, College of Pharmacy, University of Utah, Salt Lake City, UT 84112;
Maura Rossetti
Department of Pathology and Laboratory Medicine, David Geffen School of Medicine, Los Angeles, CA 90095;
Department of Pathology and Laboratory Medicine, David Geffen School of Medicine, Los Angeles, CA 90095;
Philip J. Moos
Department of Pharmacology and Toxicology, College of Pharmacy, University of Utah, Salt Lake City, UT 84112;
Frederick R. Adler
Department of Mathematics, University of Utah, Salt Lake City, UT 84112;
Jeffrey T. Chang
Department of Integrative Biology and Pharmacology, School of Medicine, School of Biomedical Informatics, The University of Texas Health Science Center at Houston, Houston, TX 77030;
Sunil Sharma2 [email protected]
Translational Oncology Research & Drug Discovery, Translational Genomics Research Institute, Phoenix, AZ 85004
Andrea H. Bild2 [email protected]
Department of Medical Oncology & Therapeutics Research, City of Hope National Medical Center, Duarte, CA 91010;


To whom correspondence may be addressed. Email: [email protected] or [email protected].
Author contributions: J.I.G., D.S., P.J.M., F.R.A., J.T.C., S.S., and A.H.B. designed research; J.I.G., P.W., L.T.P., P.A.C., N.A.L., J.A.M., G. Shrestha, M.R., G. Sunga, P.J.M., F.R.A., S.S., and A.H.B. performed research; P.W. and P.A.C. contributed new reagents/analytic tools; J.I.G., L.T.P., X.L., N.A.L., J.T.C., and A.H.B. analyzed data; J.I.G., P.W., L.T.P., D.S., X.L., F.R.A., J.T.C., S.S., and A.H.B. wrote the paper; and L.T.P. conducted the bioinformatics pipeline.
J.I.G., P.W., and L.T.P. contributed equally to this work.

Competing Interests

Competing interest statement: S.S. reports clinical trial funding from Merck for this study, which did not support the research in this paper; clinical research funding from Novartis, GSK, Millennium, MedImmune, Johnson & Johnson, Gilead Sciences, Plexxikon, Onyx, Bayer, Blueprint Medicines, XuanZhu, Incyte, Toray Industries, Celgene, Hengrui Therapeutics, OncoMed, Tesaro, AADi, and Syndax outside the submitted work; equity from Salarius Pharmaceuticals, Iterion Therapeutics, Proterus Therapeutics, ConverGene, and Stingray Therapeutics outside the submitted work; and honoraria from Blend Therapeutics, Foundation Medicine, Guardant Health, Novartis, ARIAD, US Oncology, Exelixis, Genesis Biotechnology Group, LSK, Natera, Loxo Oncology, Hengrui Therapeutics, Tarveda Therapeutics, and Dracen Pharmaceuticals outside the submitted work. D.S. reports grants from Novartis, BMS, Bioverativ, and AstraZeneca outside the submitted work and personal fees from Salarius Pharmaceuticals, Iterion Therapeutics, GlycosBio, and BMS outside the submitted work.

Metrics & Citations


Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.

Citation statements



If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Cited by


    View Options

    View options

    PDF format

    Download this article as a PDF file


    Get Access

    Login options

    Check if you have access through your login credentials or your institution to get full access on this article.

    Personal login Institutional Login

    Recommend to a librarian

    Recommend PNAS to a Librarian

    Purchase options

    Purchase this article to access the full text.

    Single Article Purchase

    Circulating immune cell phenotype dynamics reflect the strength of tumor–immune cell interactions in patients during immunotherapy
    Proceedings of the National Academy of Sciences
    • Vol. 117
    • No. 27
    • pp. 15375-16086







    Share article link

    Share on social media