Main<\/h2>\n\nChimeric antigen receptor T cells (CAR-Ts) have shown remarkable activity in the treatment of B cell malignancies1<\/a><\/sup>. With six approved therapies and hundreds in clinical development for other hematological and solid tumors, genetically engineered T cells represent a therapeutic modality changing the drug development landscape2<\/a><\/sup>. However, T cells bring unique challenges to therapeutic development. These so-called \u2018living drugs\u2019 proliferate, differentiate, actively traffic between tissues and engage in two-way communication with the patient immune system. The resultant pharmacology is different from that of small molecules or biologics, as there is little relationship between administered dose and exposure2<\/a><\/sup>.<\/p>\nThe cellular kinetics (pharmacokinetics) of circulating CAR-Ts are characterized by three distinct phases: initial expansion, followed by a rapid contraction and then slow, long-term decay3<\/a><\/sup>. The degree of cell expansion (Cmax) and long-term exposure (area under the curve (AUC)) vary widely among patients (approximtely three orders of magnitude) and are predictive of both efficacy (tumor size reduction) and toxicity4<\/a><\/sup>. However, the product-intrinsic and host-intrinsic factors mediating this pharmacology remain poorly defined. An empirical, non-linear mixed-effects model was developed to quantify the pharmacokinetics of Kymriah (tisagenlecleul, CTL019)5<\/a><\/sup> and provided as part of the biologics license application (BLA)4<\/a><\/sup>. This formulation has proven applicable to multiple other CAR-T therapies in a variety of indications6<\/a><\/sup> and has been adopted by the FDA for benchmarking7<\/a>,8<\/a><\/sup>. Although a useful tool for quantifying clinical data, the empirical equations do not account for the underlying biology and, thus, are of limited value in simulating the effects of alternate CAR-T designs, cell sources or treatment regimens. A mathematical model capable of quantitatively describing clinical data that is also based on sound biological mechanisms would be useful for the development of novel CAR-T products, as systems pharmacology modeling has proven for other therapeutic modalities9<\/a><\/sup>.<\/p>\nMathematical models of T cell\u2013tumor interactions have a long history7<\/a><\/sup> and have been adapted to describe various aspects of CAR-T pharmacology, such as antigen binding8<\/a>,10<\/a><\/sup>, intercellular signaling11<\/a><\/sup> cytokine release12<\/a><\/sup>, tissue distribution13<\/a><\/sup> and competition with host T cells for immune system reconstitution14<\/a><\/sup>. However, none of the above models adequately defines what limits cell expansion nor what underlies the wide variability in exposure and tumor response observed among patients15<\/a><\/sup>.<\/p>\nInsights can be gleaned by examining T cell dynamics in response to viral infection. Upon viral antigen encounter, antigen-specific T cells clonally expand and differentiate into cytotoxic effectors, which clear infected cells. After elimination of the pathogen, effector cells undergo a precipitous contraction phase, and a small percentage survive to form long-term memory T cells capable of self-renewal and recall responses. However, if the infection fails to resolve, chronic antigen stimulation leads to T cell exhaustion, wherein remnant T cells lose the ability to produce cytokines, kill target cells or proliferate in response to antigen16<\/a>,17<\/a><\/sup>. We hypothesize that an analogous process underlies the pharmacology of CAR-Ts.<\/p>\nWe tested this hypothesis using a conceptually simple mathematical model of T cell differentiation control, wherein an antigen-driven toggle switch regulates cell fate transitions among memory, effector and exhausted T cells. We found that the model is capable of quantitatively describing CAR-T pharmacokinetic and tumor dynamic data from multiple clinical trials and deconvolutes biological mechanisms underlying clinical variance. Specifically, we identified cell-intrinsic differences in the proliferation rate of memory cells and cytotoxic potency of effectors as the primary determinants of exposure and response, and we confirmed these mathematical inferences via analysis of bulk and single-cell RNA sequencing (scRNA-seq) data. Population exposure and response predictions were validated against registrational data from Kymriah and Yescarta. Furthermore, we demonstrate that these cell-intrinsic response-mediating differences originate in the CAR-T product using a machine learning workflow that accurately predicts patient outcomes using pre-infusion product transcriptomes. We found that functional gene signatures outperform standard T cell immunophenotyping in predictive accuracy for two CD19-targeted CAR-T products in three indications, and we summarize the relative expression of these signatures across datasets via a CAR-T response scorecard. In summary, the model predicts, de novo, clinical variance in exposure, covariates of response and the biological mechanisms underlying the pharmacology of CAR-Ts.<\/p>\n<\/div>\n<\/div>\n
\nResults<\/h2>\n\nModel structure<\/h3>\n
We consider T cells (and CAR-T products) to comprise three functionally distinct cell populations: T memory cells (T<\/i>M<\/i><\/sub>), capable of long-term self-renewal and immunological memory; T effectors (T<\/i>E<\/i><\/sub>), responsible for target-mediated cell killing; and exhausted T cells (T<\/i>X<\/i><\/sub>), lacking both killing potential and proliferative capacity. An antigen-sensing toggle switch coordinately regulates the decision of memory cells to self-renew versus differentiate, the rate of effector proliferation, exhaustion and the rate of memory cell regeneration from effectors (Methods<\/a>). This represents a conceptually simple yet biologically sound description of T cell function and regulatory control in response to immunological need, as determined by systemic antigen burden (Fig. 1a<\/a>).<\/p>\n\nFig. 1: An antigen toggle switch model of T cell regulation quantitatively describes pharmacokinetics\/pharmacodynamics behavior of CR, PR and NR patient population response to Kymriah in CLL.<\/b><\/figcaption>\n![\"Science](\"http:\/\/media.springernature.com\/lw685\/springer-static\/image\/art%3A10.1038%2Fs41587-023-01687-x\/MediaObjects\/41587_2023_1687_Fig1_HTML.png\")
<\/picture><\/a><\/div>\na<\/b>, Cartoon depiction of the model structure, comprising three populations of T cells\u2014T memory cells (T<\/i>M<\/i><\/sub>), T effector cells (T<\/i>E<\/i>1<\/sub> and T<\/i>E<\/i>2<\/sub>) and exhausted T cells (T<\/i>X<\/i><\/sub>)\u2014and B cell tumors (B<\/i>). Tumor cells express B cell antigen (B<\/i>A<\/i><\/sub>), which stimulates T cell proliferation and differentiation and inhibits the formation of T memory cells. b<\/b>, We fit the model to published pharmacokinetics\/pharmacodynamics profiles separated by response category (CR\/PR\/NR) from Fraietta et al.18<\/a><\/sup> using PSO. Model fits (curves: mean of 12 parameter sets; dark shaded areas: middle 90%) agree with both CAR-T and B cell tumor dynamics over time (dots: mean data; light shaded areas: range of data) for each of the three prototypic populations. c<\/b>, PCA plot of the logarithm of the best-fitting parameters colored by population. PC1 captures 35.3% of the variability, and PC2 captures 21.7% of the variability. d<\/b>, Sorted PC1 coefficients suggest that TK<\/i>50 (highlighted pink bar) and k<\/i>kill<\/sub>, \u03bc<\/i>M<\/i><\/sub> and d<\/i>M<\/i><\/sub> (highlighted blue bars) are the largest sources of variation between CR and NR populations. These parameters correspond to cytotoxic potency, tumor cell lysis rate, memory cell proliferation and death rates, respectively.<\/p>\n<\/div>\nFull size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\nModel parameterization: patients with CLL treated with Kymriah and grouped by response<\/h3>\n
We first sought to determine whether the mathematical description of T cell regulatory control could quantitatively capture characteristic CAR-T pharmacokinetic and tumor dynamic profiles and whether parameter estimates reveal anything about biological underpinnings of clinical variability. Fraietta et al.18<\/a><\/sup> reported mean pharmacokinetic and tumor dynamic profiles of patients with chronic lymphocytic lymphoma (CLL) treated with Kymriah (CTL019, a CD19-targeted CAR-T), grouped by complete responders (CRs), partial responders (PRs) and non-responders (NRs). We digitized the data (mean \u00b1 s.d.) and used particle swarm optimization (PSO) to estimate model parameters characterizing the three population archetypes (Fig. 1b<\/a>). Parameters were estimated 12 times per patient group. Although parameters are non-identifiable (Supplementary Information<\/a>), the clinical data were captured with good accuracy (Supplementary Fig. 5<\/a>).<\/p>\nBiological mechanisms differentiating CR, PR and NR populations<\/h3>\n
To decipher the biological mechanisms underlying the differing patient response profiles, parameter estimates from the three patient populations were first decomposed into principal components (PCs) (Fig. 1c<\/a>). Note that the three populations form relatively distinct clusters in parameter space, wherein the x<\/i> axis depicting PC1 (accounting for 35.3% of the variance) separates virtual patients by response, and the y<\/i> axis depicting PC2 (accounting for 21.7% of the variance) separates CR and NR groups from PRs. Examining the coefficients of PC1 (Fig. 1d<\/a>), the lowest value (associated with NR) is TK<\/i>50 (cytotoxic potency of effectors), and the largest positive contributions (associated with CR) is memory and effector cell turnover (proliferation and death rates; \u03bc<\/i>M<\/i><\/sub>, d<\/i>M<\/i><\/sub> and d<\/i>E<\/i>2<\/sub>). That is, in responding patients, CAR-T effectors lyse target tumor cells much more efficiently, and both memory and effector cells cycle at a higher rate. These findings are consistent with local parameter sensitivity analysis (Supplementary Fig. 6<\/a>).<\/p>\nIt is established that frequency of memory cells in CAR-T infusion products, as assessed by standard T cell immunophenotyping, is predictive of clinical response19<\/a>,20<\/a><\/sup>. This was one of the primary conclusions of Fraietta et. al.18<\/a><\/sup>. However, the PC1 loadings (Fig. 1d<\/a>) suggest that cell-intrinsic differences in memory cell function (\u03bc<\/i>M<\/i><\/sub> and d<\/i>M<\/i><\/sub>) rather than frequency (f<\/i>Tm<\/i><\/sub>) are more important determinants of response. To discern the importance of memory cell frequency versus function, we preformed two experiments. First, we attempted to fit the data under the hypothesis that the only difference between CR\/PR\/NR populations was the composition of the product (frequency of T<\/i>M<\/i><\/sub>, T<\/i>E<\/i><\/sub> and T<\/i>X<\/i><\/sub> cells), whereas the cell-intrinsic kinetic parameters are conserved (Supplementary Fig. 7<\/a>). The model does capture differences in pharmacokinetics and tumor dynamics between the populations, and the inferred CAR-T product composition is consistent with that reported by Fraietta et al.18<\/a><\/sup>. However, the magnitude of differences between the populations cannot be fully explained by this hypothesis. That is, CAR-T cell composition as defined by memory and exhausted cell frequencies alone is insufficient to explain the variance in clinical activity.<\/p>\nTo directly compare the inferred differences in memory cell function among CR\/PR\/NR groups, we simulated a dose-ranging study using purified memory cell populations from CR\/PR\/NR archetypes (Supplementary Fig. 8<\/a>). The CR memory cells produced robust and dose-dependent CAR-T expansion, persistence and tumor reduction, whereas the NR cells showed very little expansion or anti-tumor activity, and the PR memory cells display somewhat intermediate function. In sum, these results imply that, although memory cell frequency in CAR-T infusion products contributes to exposure and response, cell-intrinsic features, such as proliferative capacity, are necessary to account for the variance clinical outcomes. We next sought to identify molecular signatures that underly these cell-intrinsic features and resultant clinical variance.<\/p>\nMolecular and cellular features differentiating CR, PR and NR populations<\/h3>\n
To examine the molecular and cellular features underlying these functional differences, we used bulk RNA-seq data from the same trial18<\/a><\/sup> wherein pre-infusion CAR-T products were sequenced and annotated by response category. Differential expression analysis on the CR versus NR populations revealed biological features (gene signatures) consistent with inferred functional differences (Supplementary Figs. 9<\/a> and 10<\/a>). We confirmed findings from the original report and additionally found that the CR population is enriched in CD4+<\/sup> and CD8+<\/sup> memory cell gene signatures (defined by single-cell sequencing of thymic tissue21<\/a><\/sup>) and display heightened expression of signatures characterizing T cell proliferation, effector cytokine (interferon) signaling and IL2RB, IL7 and JAK\/STAT signaling (defined by curated pathway databases22<\/a>,23<\/a>,24<\/a><\/sup>). CAR-T cells from NR patients show heightened p53 (ref. 25<\/a><\/sup>) and DNA damage26<\/a><\/sup> signaling, pathways that may underly the proliferative deficit.<\/p>\nSingle-sample gene set enrichment analysis (ssGSEA) was subsequently used to examine distribution of the pathway and cell signatures in individual samples. The CR population is significantly enriched in the \u2018non-exhausted T cell\u2019 signature (Fig. 2a<\/a>), consistent with simulations, wherein the fraction of non-exhausted cells at day 60 (peak of anti-tumor effects) is significantly higher in the CR group (Fig. 2b<\/a>), whereas cells from the NR patients rapidly progress to exhaustion (Supplementary Fig. 11<\/a>). The simulations also align with clinical reports that CAR-T products that fail to expand in vivo show heightened expression of exhaustion markers LAG3 and PD1 (ref. 27<\/a><\/sup>).<\/p>\n\nFig. 2: ssGSEA estimates the activity of signaling pathways and enrichment of cell populations in CAR-Ts, separated by response.<\/b><\/figcaption>\n![\"Science](\"http:\/\/media.springernature.com\/lw685\/springer-static\/image\/art%3A10.1038%2Fs41587-023-01687-x\/MediaObjects\/41587_2023_1687_Fig2_HTML.png\")
<\/picture><\/a><\/div>\na<\/b>,c<\/b>\u2013f<\/b>, ssGSEA reveals differences in cell populations and signaling pathways between populations for selected cell signatures and signaling pathways (panel titles). n<\/i>\u2009=\u200931 independent samples\u2014five CR, five PR and 21 NR. b<\/b>, Using the 12 best-fitting parameter sets for each population and model simulations, we calculated the percentage of the T cell population at day 60 that is non-exhausted. The median non-exhausted T cell population at day 60 (over the 12 parameter sets) is near 100% for both CR and PR populations, whereas the median is approximately 50% for the NR population. Differences between populations were assessed using an unequal variances two-sided t<\/i>-test (P<\/i> values shown). Box plots represent median \u00b1 25th percentiles, with whiskers representing min\/max values.<\/p>\n<\/div>\n
Full size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\nWe found that CRs are differentially enriched in both CD8+<\/sup> and CD4+<\/sup> memory T cell signatures (Fig. 2c,d<\/a>), consistent with the necessity of memory cells for mediating sustained responses28<\/a><\/sup>. Note, however, that bulk sequencing data cannot resolve cell population frequencies nor discern between transcriptionally similar versus co-varying cell types (Supplementary Fig. 12<\/a>). That is, CR products may have higher frequencies of CD4+<\/sup> and CD8+<\/sup> memory cells or may contain cells with more \u2018memory-like\u2019 transcriptomes at similar frequencies. The CR population also shows heightened IL2RB and IL7R signaling (Fig. 2e,f<\/a>), indicating that the CR cell products may show heightened sensitivity to the correspondent cytokines. Notably, IL2 and IL7 are common components of CAR-T expansion media29<\/a><\/sup>, and peak serum IL7 concentration is predictive of CD19 CAR-T exposure and progression-free survival30<\/a><\/sup>. Although the results shown in Fig. 2<\/a> are statistically significant, the ssGSEA distributions overlap between response categories. Thus, in addition to the limitations of bulk sequencing data, none of the gene signatures assessed could serve as univariable predictors of patient response.<\/p>\nCell-intrinsic functional differences mediating CAR-T clinical response<\/h3>\n
To deconvolute the role of cell frequency versus function in mediating response, we leveraged two recently published clinical studies containing scRNA-seq data of pre-infusion, autologous CD19 CAR-T products matched with clinical outcomes. Bai et al.31<\/a><\/sup> reported data for 12 patients with acute lymphoblastic leukemia (ALL) treated with a CD19 CAR-T product analogous to Kymriah\u2014five CRs, two NRs and five patients who relapsed (RL). Haradhvala et al.32<\/a><\/sup> reported data for 32 patients with large B cell lymphoma (LBCL) treated with either Kymriah (n<\/i>\u2009=\u200913) or Yescarta (n<\/i>\u2009=\u200919). For the Kymriah-treated group, there were six CRs and seven NRs; for the Yescarta-treated group, there were 11 CRs, one PR and seven NRs.<\/p>\nExamination of uniform manifold approximation and projection (UMAP) projections of the three datasets (Kymriah in ALL, Kymriah in LBCL and Yescarta in LBCL) reveals some separation of response categories in transcriptome space, particularly in ALL (Fig. 3a,d,g<\/a>). To assess whether response separation is attributable to differences in T cell composition, we assigned cell type labels by mapping expression profiles of the individual cells to annotated tumor-infiltrating lymphocyte populations via ProjecTILs33<\/a><\/sup>. Most CD8+<\/sup> cells in all three datasets are classified as T effector memory (Tem) or T exhausted (Tex), but there are no consistent differences in composition by response category (Supplementary Fig. 13a\u2013c<\/a>). For example, the frequency of cells annotated as exhausted is significantly higher in the NR\/RL categories as compared to CR in the ALL data (P<\/i>\u2009<\u20090.05, mean 4.4% versus 8.7%, respectively; Fig. 3b,e,h<\/a>). However, this pattern does not hold for the LBCL data, and the modest effect size is insufficient to account for the vast disparity in clinical outcomes. We used the cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq) antibody tag data provided by Bai et al.34<\/a><\/sup> to assign early memory (Tmem: CD8+<\/sup>CD45RO\u2212<\/sup>CD27+<\/sup>) and exhausted (CD8+<\/sup>PD1+<\/sup>) cell annotations by immunophenotype, reported to be predictive of response in CLL18<\/a><\/sup>. Although exhausted cell annotations by ProjecTILs and immunophenotype were notably concordant (6.7% versus 5.9% of total cells), cell frequencies did not differ by response category in ALL (Supplementary Fig. 13d,e<\/a>).<\/p>\n\nFig. 3: scRNA-seq of pre-infusion CAR-T products reveals cell-intrinsic defects associated with non-durable response.<\/b><\/figcaption>\n![\"Science](\"http:\/\/media.springernature.com\/lw685\/springer-static\/image\/art%3A10.1038%2Fs41587-023-01687-x\/MediaObjects\/41587_2023_1687_Fig3_HTML.png\")
<\/picture><\/a><\/div>\nUMAP projections of three datasets representing Kymriah in ALL (a<\/b>\u2013c<\/b>), Kymriah in LBCL (d<\/b>\u2013f<\/b>) and Yescarta in LBCL (g<\/b>\u2013i<\/b>). a<\/b>,d<\/b>,g<\/b>, UMAP projections annotated by response category. b<\/b>,e<\/b>,h<\/b>, UMAP projections annotated as exhausted using ProjectTILs33<\/a><\/sup>. c<\/b>,f<\/b>,i<\/b>, UMAP projections annotated for high (above mean) or low (below mean) CAR-T cell dysfunction signature from Good et al.35<\/a><\/sup>. j<\/b>, GSEA for select pathways, comparing both exhausted versus non-exhausted and CR versus PR\/RL\/NR categories within cells annotated as T effector memory (Tem via ProjecTILs) or early memory (Tmem; CD8+<\/sup>CD45RO\u2212<\/sup>CD27+<\/sup> via CITE-seq). A positive normalized enrichment score (NES, blue) indicates higher enrichment in CR\/non-exhausted cells. *NR\u2009=\u2009NR\/RL or NR\/PR. P<\/i> values were calculated by Kolmogorov\u2013Smirnov tests implemented in GSEA.<\/p>\n<\/div>\nFull size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\nTo probe cell-intrinsic function, we annotated cells using a \u2018CAR-T dysfunction\u2019 signature, characteristic of functionally exhausted CAR-T cells with reduced proliferative and cytotoxic capacity35<\/a><\/sup>. Visually, the dysfunction signature is dispersed throughout response categories and not restricted to exhausted regions (Fig. 3g,h,i<\/a>). Interrogating cell-intrinsic functional differences at a deeper resolution, we performed differential gene expression analysis on T sub-cell populations (annotated both by transcriptome and immunophenotype), followed by pathway enrichment for select gene signatures (Fig. 3j<\/a>). As a control, we first assessed differences between cells annotated as exhausted versus non-exhausted. Exhausted cells are consistently enriched in the CAR-T dysfunction signature across datasets, whereas the \u2018exhausted T cell\u2019 and \u2018P53 signaling\u2019 signatures appear specific to the ALL-exhausted cells. Conversely, non-exhausted cells show disparate enrichment for the \u2018early memory T cell\u2019 signature as well as cytokine production and inflammatory response signatures, hallmarks of T cell functional potency.<\/p>\nComparing cell populations from the CR versus NR\/PR\/RL categories reveals a consistent pattern across datasets. Focusing either on effector memory or early memory (CD8+<\/sup>CD45RA\u2212<\/sup>CD27
Chimeric antigen receptor T cells (CAR-Ts) have shown remarkable activity in the treatment of B cell malignancies1<\/a><\/sup>. With six approved therapies and hundreds in clinical development for other hematological and solid tumors, genetically engineered T cells represent a therapeutic modality changing the drug development landscape2<\/a><\/sup>. However, T cells bring unique challenges to therapeutic development. These so-called \u2018living drugs\u2019 proliferate, differentiate, actively traffic between tissues and engage in two-way communication with the patient immune system. The resultant pharmacology is different from that of small molecules or biologics, as there is little relationship between administered dose and exposure2<\/a><\/sup>.<\/p>\n The cellular kinetics (pharmacokinetics) of circulating CAR-Ts are characterized by three distinct phases: initial expansion, followed by a rapid contraction and then slow, long-term decay3<\/a><\/sup>. The degree of cell expansion (Cmax) and long-term exposure (area under the curve (AUC)) vary widely among patients (approximtely three orders of magnitude) and are predictive of both efficacy (tumor size reduction) and toxicity4<\/a><\/sup>. However, the product-intrinsic and host-intrinsic factors mediating this pharmacology remain poorly defined. An empirical, non-linear mixed-effects model was developed to quantify the pharmacokinetics of Kymriah (tisagenlecleul, CTL019)5<\/a><\/sup> and provided as part of the biologics license application (BLA)4<\/a><\/sup>. This formulation has proven applicable to multiple other CAR-T therapies in a variety of indications6<\/a><\/sup> and has been adopted by the FDA for benchmarking7<\/a>,8<\/a><\/sup>. Although a useful tool for quantifying clinical data, the empirical equations do not account for the underlying biology and, thus, are of limited value in simulating the effects of alternate CAR-T designs, cell sources or treatment regimens. A mathematical model capable of quantitatively describing clinical data that is also based on sound biological mechanisms would be useful for the development of novel CAR-T products, as systems pharmacology modeling has proven for other therapeutic modalities9<\/a><\/sup>.<\/p>\n Mathematical models of T cell\u2013tumor interactions have a long history7<\/a><\/sup> and have been adapted to describe various aspects of CAR-T pharmacology, such as antigen binding8<\/a>,10<\/a><\/sup>, intercellular signaling11<\/a><\/sup> cytokine release12<\/a><\/sup>, tissue distribution13<\/a><\/sup> and competition with host T cells for immune system reconstitution14<\/a><\/sup>. However, none of the above models adequately defines what limits cell expansion nor what underlies the wide variability in exposure and tumor response observed among patients15<\/a><\/sup>.<\/p>\n Insights can be gleaned by examining T cell dynamics in response to viral infection. Upon viral antigen encounter, antigen-specific T cells clonally expand and differentiate into cytotoxic effectors, which clear infected cells. After elimination of the pathogen, effector cells undergo a precipitous contraction phase, and a small percentage survive to form long-term memory T cells capable of self-renewal and recall responses. However, if the infection fails to resolve, chronic antigen stimulation leads to T cell exhaustion, wherein remnant T cells lose the ability to produce cytokines, kill target cells or proliferate in response to antigen16<\/a>,17<\/a><\/sup>. We hypothesize that an analogous process underlies the pharmacology of CAR-Ts.<\/p>\n We tested this hypothesis using a conceptually simple mathematical model of T cell differentiation control, wherein an antigen-driven toggle switch regulates cell fate transitions among memory, effector and exhausted T cells. We found that the model is capable of quantitatively describing CAR-T pharmacokinetic and tumor dynamic data from multiple clinical trials and deconvolutes biological mechanisms underlying clinical variance. Specifically, we identified cell-intrinsic differences in the proliferation rate of memory cells and cytotoxic potency of effectors as the primary determinants of exposure and response, and we confirmed these mathematical inferences via analysis of bulk and single-cell RNA sequencing (scRNA-seq) data. Population exposure and response predictions were validated against registrational data from Kymriah and Yescarta. Furthermore, we demonstrate that these cell-intrinsic response-mediating differences originate in the CAR-T product using a machine learning workflow that accurately predicts patient outcomes using pre-infusion product transcriptomes. We found that functional gene signatures outperform standard T cell immunophenotyping in predictive accuracy for two CD19-targeted CAR-T products in three indications, and we summarize the relative expression of these signatures across datasets via a CAR-T response scorecard. In summary, the model predicts, de novo, clinical variance in exposure, covariates of response and the biological mechanisms underlying the pharmacology of CAR-Ts.<\/p>\n<\/div>\n<\/div>\n We consider T cells (and CAR-T products) to comprise three functionally distinct cell populations: T memory cells (T<\/i>M<\/i><\/sub>), capable of long-term self-renewal and immunological memory; T effectors (T<\/i>E<\/i><\/sub>), responsible for target-mediated cell killing; and exhausted T cells (T<\/i>X<\/i><\/sub>), lacking both killing potential and proliferative capacity. An antigen-sensing toggle switch coordinately regulates the decision of memory cells to self-renew versus differentiate, the rate of effector proliferation, exhaustion and the rate of memory cell regeneration from effectors (Methods<\/a>). This represents a conceptually simple yet biologically sound description of T cell function and regulatory control in response to immunological need, as determined by systemic antigen burden (Fig. 1a<\/a>).<\/p>\n a<\/b>, Cartoon depiction of the model structure, comprising three populations of T cells\u2014T memory cells (T<\/i>M<\/i><\/sub>), T effector cells (T<\/i>E<\/i>1<\/sub> and T<\/i>E<\/i>2<\/sub>) and exhausted T cells (T<\/i>X<\/i><\/sub>)\u2014and B cell tumors (B<\/i>). Tumor cells express B cell antigen (B<\/i>A<\/i><\/sub>), which stimulates T cell proliferation and differentiation and inhibits the formation of T memory cells. b<\/b>, We fit the model to published pharmacokinetics\/pharmacodynamics profiles separated by response category (CR\/PR\/NR) from Fraietta et al.18<\/a><\/sup> using PSO. Model fits (curves: mean of 12 parameter sets; dark shaded areas: middle 90%) agree with both CAR-T and B cell tumor dynamics over time (dots: mean data; light shaded areas: range of data) for each of the three prototypic populations. c<\/b>, PCA plot of the logarithm of the best-fitting parameters colored by population. PC1 captures 35.3% of the variability, and PC2 captures 21.7% of the variability. d<\/b>, Sorted PC1 coefficients suggest that TK<\/i>50 (highlighted pink bar) and k<\/i>kill<\/sub>, \u03bc<\/i>M<\/i><\/sub> and d<\/i>M<\/i><\/sub> (highlighted blue bars) are the largest sources of variation between CR and NR populations. These parameters correspond to cytotoxic potency, tumor cell lysis rate, memory cell proliferation and death rates, respectively.<\/p>\n<\/div>\n Full size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\n We first sought to determine whether the mathematical description of T cell regulatory control could quantitatively capture characteristic CAR-T pharmacokinetic and tumor dynamic profiles and whether parameter estimates reveal anything about biological underpinnings of clinical variability. Fraietta et al.18<\/a><\/sup> reported mean pharmacokinetic and tumor dynamic profiles of patients with chronic lymphocytic lymphoma (CLL) treated with Kymriah (CTL019, a CD19-targeted CAR-T), grouped by complete responders (CRs), partial responders (PRs) and non-responders (NRs). We digitized the data (mean \u00b1 s.d.) and used particle swarm optimization (PSO) to estimate model parameters characterizing the three population archetypes (Fig. 1b<\/a>). Parameters were estimated 12 times per patient group. Although parameters are non-identifiable (Supplementary Information<\/a>), the clinical data were captured with good accuracy (Supplementary Fig. 5<\/a>).<\/p>\n To decipher the biological mechanisms underlying the differing patient response profiles, parameter estimates from the three patient populations were first decomposed into principal components (PCs) (Fig. 1c<\/a>). Note that the three populations form relatively distinct clusters in parameter space, wherein the x<\/i> axis depicting PC1 (accounting for 35.3% of the variance) separates virtual patients by response, and the y<\/i> axis depicting PC2 (accounting for 21.7% of the variance) separates CR and NR groups from PRs. Examining the coefficients of PC1 (Fig. 1d<\/a>), the lowest value (associated with NR) is TK<\/i>50 (cytotoxic potency of effectors), and the largest positive contributions (associated with CR) is memory and effector cell turnover (proliferation and death rates; \u03bc<\/i>M<\/i><\/sub>, d<\/i>M<\/i><\/sub> and d<\/i>E<\/i>2<\/sub>). That is, in responding patients, CAR-T effectors lyse target tumor cells much more efficiently, and both memory and effector cells cycle at a higher rate. These findings are consistent with local parameter sensitivity analysis (Supplementary Fig. 6<\/a>).<\/p>\n It is established that frequency of memory cells in CAR-T infusion products, as assessed by standard T cell immunophenotyping, is predictive of clinical response19<\/a>,20<\/a><\/sup>. This was one of the primary conclusions of Fraietta et. al.18<\/a><\/sup>. However, the PC1 loadings (Fig. 1d<\/a>) suggest that cell-intrinsic differences in memory cell function (\u03bc<\/i>M<\/i><\/sub> and d<\/i>M<\/i><\/sub>) rather than frequency (f<\/i>Tm<\/i><\/sub>) are more important determinants of response. To discern the importance of memory cell frequency versus function, we preformed two experiments. First, we attempted to fit the data under the hypothesis that the only difference between CR\/PR\/NR populations was the composition of the product (frequency of T<\/i>M<\/i><\/sub>, T<\/i>E<\/i><\/sub> and T<\/i>X<\/i><\/sub> cells), whereas the cell-intrinsic kinetic parameters are conserved (Supplementary Fig. 7<\/a>). The model does capture differences in pharmacokinetics and tumor dynamics between the populations, and the inferred CAR-T product composition is consistent with that reported by Fraietta et al.18<\/a><\/sup>. However, the magnitude of differences between the populations cannot be fully explained by this hypothesis. That is, CAR-T cell composition as defined by memory and exhausted cell frequencies alone is insufficient to explain the variance in clinical activity.<\/p>\n To directly compare the inferred differences in memory cell function among CR\/PR\/NR groups, we simulated a dose-ranging study using purified memory cell populations from CR\/PR\/NR archetypes (Supplementary Fig. 8<\/a>). The CR memory cells produced robust and dose-dependent CAR-T expansion, persistence and tumor reduction, whereas the NR cells showed very little expansion or anti-tumor activity, and the PR memory cells display somewhat intermediate function. In sum, these results imply that, although memory cell frequency in CAR-T infusion products contributes to exposure and response, cell-intrinsic features, such as proliferative capacity, are necessary to account for the variance clinical outcomes. We next sought to identify molecular signatures that underly these cell-intrinsic features and resultant clinical variance.<\/p>\n To examine the molecular and cellular features underlying these functional differences, we used bulk RNA-seq data from the same trial18<\/a><\/sup> wherein pre-infusion CAR-T products were sequenced and annotated by response category. Differential expression analysis on the CR versus NR populations revealed biological features (gene signatures) consistent with inferred functional differences (Supplementary Figs. 9<\/a> and 10<\/a>). We confirmed findings from the original report and additionally found that the CR population is enriched in CD4+<\/sup> and CD8+<\/sup> memory cell gene signatures (defined by single-cell sequencing of thymic tissue21<\/a><\/sup>) and display heightened expression of signatures characterizing T cell proliferation, effector cytokine (interferon) signaling and IL2RB, IL7 and JAK\/STAT signaling (defined by curated pathway databases22<\/a>,23<\/a>,24<\/a><\/sup>). CAR-T cells from NR patients show heightened p53 (ref. 25<\/a><\/sup>) and DNA damage26<\/a><\/sup> signaling, pathways that may underly the proliferative deficit.<\/p>\n Single-sample gene set enrichment analysis (ssGSEA) was subsequently used to examine distribution of the pathway and cell signatures in individual samples. The CR population is significantly enriched in the \u2018non-exhausted T cell\u2019 signature (Fig. 2a<\/a>), consistent with simulations, wherein the fraction of non-exhausted cells at day 60 (peak of anti-tumor effects) is significantly higher in the CR group (Fig. 2b<\/a>), whereas cells from the NR patients rapidly progress to exhaustion (Supplementary Fig. 11<\/a>). The simulations also align with clinical reports that CAR-T products that fail to expand in vivo show heightened expression of exhaustion markers LAG3 and PD1 (ref. 27<\/a><\/sup>).<\/p>\n a<\/b>,c<\/b>\u2013f<\/b>, ssGSEA reveals differences in cell populations and signaling pathways between populations for selected cell signatures and signaling pathways (panel titles). n<\/i>\u2009=\u200931 independent samples\u2014five CR, five PR and 21 NR. b<\/b>, Using the 12 best-fitting parameter sets for each population and model simulations, we calculated the percentage of the T cell population at day 60 that is non-exhausted. The median non-exhausted T cell population at day 60 (over the 12 parameter sets) is near 100% for both CR and PR populations, whereas the median is approximately 50% for the NR population. Differences between populations were assessed using an unequal variances two-sided t<\/i>-test (P<\/i> values shown). Box plots represent median \u00b1 25th percentiles, with whiskers representing min\/max values.<\/p>\n<\/div>\n Full size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\n We found that CRs are differentially enriched in both CD8+<\/sup> and CD4+<\/sup> memory T cell signatures (Fig. 2c,d<\/a>), consistent with the necessity of memory cells for mediating sustained responses28<\/a><\/sup>. Note, however, that bulk sequencing data cannot resolve cell population frequencies nor discern between transcriptionally similar versus co-varying cell types (Supplementary Fig. 12<\/a>). That is, CR products may have higher frequencies of CD4+<\/sup> and CD8+<\/sup> memory cells or may contain cells with more \u2018memory-like\u2019 transcriptomes at similar frequencies. The CR population also shows heightened IL2RB and IL7R signaling (Fig. 2e,f<\/a>), indicating that the CR cell products may show heightened sensitivity to the correspondent cytokines. Notably, IL2 and IL7 are common components of CAR-T expansion media29<\/a><\/sup>, and peak serum IL7 concentration is predictive of CD19 CAR-T exposure and progression-free survival30<\/a><\/sup>. Although the results shown in Fig. 2<\/a> are statistically significant, the ssGSEA distributions overlap between response categories. Thus, in addition to the limitations of bulk sequencing data, none of the gene signatures assessed could serve as univariable predictors of patient response.<\/p>\n To deconvolute the role of cell frequency versus function in mediating response, we leveraged two recently published clinical studies containing scRNA-seq data of pre-infusion, autologous CD19 CAR-T products matched with clinical outcomes. Bai et al.31<\/a><\/sup> reported data for 12 patients with acute lymphoblastic leukemia (ALL) treated with a CD19 CAR-T product analogous to Kymriah\u2014five CRs, two NRs and five patients who relapsed (RL). Haradhvala et al.32<\/a><\/sup> reported data for 32 patients with large B cell lymphoma (LBCL) treated with either Kymriah (n<\/i>\u2009=\u200913) or Yescarta (n<\/i>\u2009=\u200919). For the Kymriah-treated group, there were six CRs and seven NRs; for the Yescarta-treated group, there were 11 CRs, one PR and seven NRs.<\/p>\n Examination of uniform manifold approximation and projection (UMAP) projections of the three datasets (Kymriah in ALL, Kymriah in LBCL and Yescarta in LBCL) reveals some separation of response categories in transcriptome space, particularly in ALL (Fig. 3a,d,g<\/a>). To assess whether response separation is attributable to differences in T cell composition, we assigned cell type labels by mapping expression profiles of the individual cells to annotated tumor-infiltrating lymphocyte populations via ProjecTILs33<\/a><\/sup>. Most CD8+<\/sup> cells in all three datasets are classified as T effector memory (Tem) or T exhausted (Tex), but there are no consistent differences in composition by response category (Supplementary Fig. 13a\u2013c<\/a>). For example, the frequency of cells annotated as exhausted is significantly higher in the NR\/RL categories as compared to CR in the ALL data (P<\/i>\u2009<\u20090.05, mean 4.4% versus 8.7%, respectively; Fig. 3b,e,h<\/a>). However, this pattern does not hold for the LBCL data, and the modest effect size is insufficient to account for the vast disparity in clinical outcomes. We used the cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq) antibody tag data provided by Bai et al.34<\/a><\/sup> to assign early memory (Tmem: CD8+<\/sup>CD45RO\u2212<\/sup>CD27+<\/sup>) and exhausted (CD8+<\/sup>PD1+<\/sup>) cell annotations by immunophenotype, reported to be predictive of response in CLL18<\/a><\/sup>. Although exhausted cell annotations by ProjecTILs and immunophenotype were notably concordant (6.7% versus 5.9% of total cells), cell frequencies did not differ by response category in ALL (Supplementary Fig. 13d,e<\/a>).<\/p>\n UMAP projections of three datasets representing Kymriah in ALL (a<\/b>\u2013c<\/b>), Kymriah in LBCL (d<\/b>\u2013f<\/b>) and Yescarta in LBCL (g<\/b>\u2013i<\/b>). a<\/b>,d<\/b>,g<\/b>, UMAP projections annotated by response category. b<\/b>,e<\/b>,h<\/b>, UMAP projections annotated as exhausted using ProjectTILs33<\/a><\/sup>. c<\/b>,f<\/b>,i<\/b>, UMAP projections annotated for high (above mean) or low (below mean) CAR-T cell dysfunction signature from Good et al.35<\/a><\/sup>. j<\/b>, GSEA for select pathways, comparing both exhausted versus non-exhausted and CR versus PR\/RL\/NR categories within cells annotated as T effector memory (Tem via ProjecTILs) or early memory (Tmem; CD8+<\/sup>CD45RO\u2212<\/sup>CD27+<\/sup> via CITE-seq). A positive normalized enrichment score (NES, blue) indicates higher enrichment in CR\/non-exhausted cells. *NR\u2009=\u2009NR\/RL or NR\/PR. P<\/i> values were calculated by Kolmogorov\u2013Smirnov tests implemented in GSEA.<\/p>\n<\/div>\n Full size image<\/span><\/a><\/p>\n<\/figure>\n<\/div>\n To probe cell-intrinsic function, we annotated cells using a \u2018CAR-T dysfunction\u2019 signature, characteristic of functionally exhausted CAR-T cells with reduced proliferative and cytotoxic capacity35<\/a><\/sup>. Visually, the dysfunction signature is dispersed throughout response categories and not restricted to exhausted regions (Fig. 3g,h,i<\/a>). Interrogating cell-intrinsic functional differences at a deeper resolution, we performed differential gene expression analysis on T sub-cell populations (annotated both by transcriptome and immunophenotype), followed by pathway enrichment for select gene signatures (Fig. 3j<\/a>). As a control, we first assessed differences between cells annotated as exhausted versus non-exhausted. Exhausted cells are consistently enriched in the CAR-T dysfunction signature across datasets, whereas the \u2018exhausted T cell\u2019 and \u2018P53 signaling\u2019 signatures appear specific to the ALL-exhausted cells. Conversely, non-exhausted cells show disparate enrichment for the \u2018early memory T cell\u2019 signature as well as cytokine production and inflammatory response signatures, hallmarks of T cell functional potency.<\/p>\n Comparing cell populations from the CR versus NR\/PR\/RL categories reveals a consistent pattern across datasets. Focusing either on effector memory or early memory (CD8+<\/sup>CD45RA\u2212<\/sup>CD27Results<\/h2>\n
Model structure<\/h3>\n
<\/picture><\/a><\/div>\nModel parameterization: patients with CLL treated with Kymriah and grouped by response<\/h3>\n
Biological mechanisms differentiating CR, PR and NR populations<\/h3>\n
Molecular and cellular features differentiating CR, PR and NR populations<\/h3>\n
<\/picture><\/a><\/div>\nCell-intrinsic functional differences mediating CAR-T clinical response<\/h3>\n
<\/picture><\/a><\/div>\n