## Abstract

CAR T cells are engineered to bind and destroy tumor cells by targeting overexpressed surface antigens. However, healthy cells expressing lower abundances of these antigens can also be lysed by CAR T cells. Various CAR T cell designs increase tumor cell elimination, whereas reducing damage to healthy cells. However, these efforts are costly and labor-intensive, constraining systematic exploration of potential hypotheses. We develop a protein abundance structured population dynamic model for CAR T cells (PASCAR), a framework that combines multiscale population dynamic models and multi-objective optimization approaches with data from cytometry and cytotoxicity assays to systematically explore the design space of constitutive and tunable CAR T cells. PASCAR can quantitatively describe in vitro and in vivo results for constitutive and inducible CAR T cells and can successfully predict experiments outside the training data. Our exploration of the CAR design space reveals that optimal CAR affinities in the intermediate range of dissociation constants effectively reduce healthy cell lysis, whereas maintaining high tumor cell-killing rates. Furthermore, our modeling offers guidance for optimizing CAR expressions in synthetic notch CAR T cells. PASCAR can be extended to other CAR immune cells.

## Introduction

Chimeric antigen receptor (CAR) T cells have been widely successful in treating hematologic cancers (Lim and June, 2017; Ellis et al, 2021). CAR T cells are engineered to express CAR molecules which bind to antigens overexpressed in tumor cells and stimulate cytotoxic T cell responses. However, healthy cells can express such antigens at low copy numbers (Yasui et al, 1988; Klichinsky et al, 2020) and can become targets of CAR T cell cytotoxicity (Morgan et al, 2010; Hernandez-Lopez et al, 2021; Duan et al, 2022). Such off-target destruction to healthy tissues is a major source of severe immune-related adverse effects in patients undergoing CAR T cell therapy (Parkhurst et al, 2011; Guo et al, 2018) and can become a major issue for solid tumors where viral organs are damaged by CAR T cells. Therefore, generating an optimal CAR T cell response that maximizes tumor cell elimination, whereas minimizing off-target destruction has been a longstanding pursuit of CAR T cell therapies. A wide range of design strategies to engineer CAR T cells have been proposed, including constitutive co-expression multiple CARs (Srivastava & Riddell, 2015; Cho et al, 2018) that sense several target antigens overexpressed by cancer cells or generation of inducible CAR expressions that tune copy numbers or abundances of CARs (Roybal et al, 2016; Hernandez-Lopez et al, 2021) depending on the abundances of target antigens present on target cells. However, most of these design strategies are conceived intuitively, and thus limit systematic and wide explorations of the design space for CAR T cells.

Optimal design strategies in economics or engineering often involve optimization of conflicting objectives where the best trade-offs between objectives are explored by subjecting computational/mathematical models to multiple objective optimization or Pareto optimization (Atashkari et al, 2005; Deb, 2011; Dworczak et al, 2021). However, the application of similar approaches for designing CAR T cells face a major challenge because of the lack of experimentally validated and computationally efficient multiscale mathematical/computational models that can describe CAR T cell responses against target cells. Though a large variety of population dynamic models based on ordinary differential equations (ODEs) (Chaudhury et al, 2020) have been developed to describe kinetics of populations of CAR T cells interacting with cancer cells in vitro or in vivo, most of these models do not explicitly include CAR ligand affinities, CAR abundances, co-receptors, and their cognate ligand molecules or biochemical signaling processes initiated upon engagement of CARs with cognate ligands. Because the design of CAR constructs often involves manipulation of the CAR affinities, CAR abundances or T cell signaling processes, it is difficult to systematically explore the design space of CAR constructs with such existing models.

CAR abundances in single CAR T cells can vary over 1,000-fold in CAR T cell populations (Hernandez-Lopez et al, 2021) which could play an important role in regulating the response of CAR T cell populations against target cells. For example, in vitro cytotoxic assays of T cells with inducible CAR expressions showed that an increase in mean CAR abundances by less than 1.3-fold can produce over threefold increase in the rate of lysis of target cells in vitro. Furthermore, abundances of target antigens (e.g., CD19) on tumor cells in patients can also vary over 100-fold and lead to disparate (responders or nonresponders) outcomes in patients undergoing CAR T cell therapy (Majzner et al, 2020). A handful of recent pharmacodynamic models (Singh et al, 2020; Qi et al, 2022) incorporated CAR ligand interactions by considering mean abundances of CARs and their cognate ligands, however, these models are unable to capture the variations of single-cell CAR abundances or the T cell signaling kinetics. Detailed agent-based models incorporating details of receptor–ligand interactions, T cell signaling kinetics and cell metabolism have been developed to describe CAR T cell responses (Prybutok et al, 2022), however, quantitative validation of these models with experiments is challenging because of the presence of many difficult-to-calibrate model parameters and computationally intensive nature of the simulations.

We develop an experimentally validated protein abundance structured population dynamic model for CAR T cells (PASCAR) that integrates molecular receptor–ligand interactions to single CAR T cell signaling and activation to population kinetics of interacting CAR T cells and target cells. PASCAR integrates CAR–ligand interactions and the ensuing signaling kinetics with a minimal but generalizable mechanistic modeling approach using ODEs where model parameters can be well estimated from routinely carried out experiments such as cytotoxicity assays and flow cytometry. We demonstrate PASCAR can quantitatively describe in vitro results for constitutive and inducible CAR T cells and can successfully predict experiments outside the training data. Constitutive CAR T cells constitutively express CAR molecules, whereas inducible CAR T cells generate CAR expressions depending on their interaction with target antigens on target cells. PASCAR is then combined with a Pareto optimization that includes the trade-off between lysis of tumor and healthy cells to explore the design space of CAR constructs. Our investigations show CAR–ligand affinities with dissociation constants in the micromolar range can dramatically decrease healthy cell lysis but sustain a high rate of tumor cell killing. The proposed framework can be extended to model responses in other CAR immune cells (Klichinsky et al, 2020; Liu et al, 2020; Kararoudi et al, 2022).

### Model development

We developed a protein abundance structured model (PASCAR) for describing interacting populations of CAR CD8^{+} T cells and target cells. In the model, individual CAR T cells interact with single-target cells where the interaction between a CAR CD8^{+} T cell and a target cell is initiated by binding of the CAR with its cognate ligand such as HER-2. Once the CAR–ligand complex is formed, it goes through a series of chemical modifications such as phosphorylation of tyrosine residues in the CAR-associated CD3ζ adaptors because of signaling reactions (Harris et al, 2018; Rohrs et al, 2020). These modifications eventually lead to the release of lytic granules by the CAR T cell which induces disintegration of the target cell membrane and eventual death of the target cell (Davenport et al, 2018). The signaling reactions can also induce proliferation of the CAR CD8^{+} T cells (Sahoo et al, 2020; Faude et al, 2021). There are a multitude of interconnected signaling processes involving many proteins, lipids, and ions that link the formation of the CAR–ligand complex to the lysis of target cells or to CAR T cell proliferation. In the model, we make simplifying assumptions to model these signaling reactions to relate the rate of target cell lysis and CAR CD8^{+} T cell proliferation to the abundances of CAR–ligand signaling complexes. Because signaling reactions occur at faster time scales (∼minutes) (Huse et al, 2007) compared with the time scales (∼hours) of target cell lysis (Weigelin et al, 2021) or CAR T cell proliferation (Obst, 2015), we assume that the rates of target cell lysis and T cell proliferation are influenced by the steady state values of the abundances of CAR–ligand complexes in our model. To simplify the notation, we denote CAR and its cognate ligand by R and H, respectively, and denote the single-cell copy numbers or abundances of these proteins by the italicized versions of the symbols. In PASCAR, a CAR T cell with *R* number of CARs interacts with a target cell with expressing *H* number of cognate ligands. R binds with H at a rate *k*_{on} to create the complex, R–H, which unbinds at a rate *k*_{off} (Fig 1). The formation of the complex R–H induces a series of signaling reactions in the CAR T cell that eventually leads to the activation of the CAR T cell.

We propose two models, Model NKP and Model KP, to investigate different mechanisms of CAR T cell activation because of CAR binding to cognate ligands. In Model NKP, CAR molecules in CAR T cells interact with cognate ligands on target cells to form a R–H complex and activation of a CAR T cell is assumed to be proportional to the steady state abundance of the R–H complex (Fig 1). In Model KP, we use an approximate model for CAR T cell signal transduction based on McKeithan’s kinetic proofreading (KP) model (McKeithan, 1995). In this model, the R–H complex formed in the CAR T cell undergoes *N* number of modifications representing chemical modifications by downstream signaling reactions to create an active complex C_{N} (Fig 1). Production of C_{N} leads to cytotoxic response and proliferation of CAR T cells. The above series of reactions approximate signaling reactions in CAR T cells initiated by the formation of the CAR–ligand complex that produces activation (e.g., phosphorylation) of key signaling proteins (Salter et al, 2018, 2021) such as Zap70 or PLCγ or SLP-76 crucial for CAR T cell activation. The chemical modifications of the R–H complex are described by first-order reactions (Fig 1) with rate *k*_{p}. At any chemically modified state of the complex, unbinding the ligand leads to immediate reversion of all the modifications in the complex and the complex reverts to the native unbound ligand and receptor state. This step is known as the KP step which creates a sharp increase in the steady state number of C_{N} as *k*_{off} decreases. The KP model gives rise to a waiting time (τ_{w}) that the receptor–ligand complex should last to generate productive signaling—receptor–ligand complexes with lifetimes (∼1/k_{off}) larger than this waiting time (i.e., 1/k_{off} > τ_{w} ∝1/k_{p}) generates active complexes C_{N} at greater abundances compared with short-lived (1/k_{off} ≪ τ_{w}) receptor complexes. Below, we describe the kinetics of populations of interacting CAR T cells and target cells for Model NKP and Model KP using ODEs.

#### Model NKP

In this model, the rate of CAR T cell proliferation or the rate at which a CAR T cell lyses an interacting target cell is proportional to the steady state abundance of the R–H complex or C_{0}. The abundance (*C*_{0}) of the R–H complex formed at the steady state as a function of *R*, *H*, and the dissociation constant *K*_{D} = *k*_{off}/*k*_{on} is given by

We assume that the rate of CAR T cell proliferation ρ_{RH} and the rate of target cell lysis λ_{RH} are given by_{c} and λ_{c} are the proportionality constants.

#### Model KP

The rates of CAR T cell proliferation and target cell lysis are proportional to the steady state abundance of the end complex C_{N} which as a function of the steady state abundance of R-H or C_{0}^{(KP)}, k_{p}, and k_{off} is given by,*C*_{0}^{(KP)} is given by,

The derivations of Equations (1) and (2) are shown in Supplemental Data 1. The proliferation and lysis rates in Model KP are given by

### Supplemental Data 1.

Derivation of Equations (1) and (2).[LSA-2023-02171_Supplemental_Data_1.pdf]

Now, we set up kinetic equations for a population of target cells interacting with a population of CAR T cells. The target cells can represent healthy or tumor cells. We consider a population of CAR T cells where individual CAR T cells express *R* copy number of CARs within a range [*R*_{min}, *R*_{max}]. The CAR T cells interact with a population of target cells where any single-target cell expresses *H* copy number of cognate ligands within a range [*H*_{min}, *H*_{max}]. The target cells replicate with a rate *r* and the target cell population decreases as they are lysed by interacting CAR T cells. The CAR T cells proliferate because of their interaction with the target cells, and the CAR T cells die with a constant rate δ. The population kinetics can be described in terms of the number (*U*_{H}) of target cells each carrying *H* copy number of cognate ligands and the number (*T*_{R}) of CAR T cells each carrying *R* copy number of CARs and as follows:

Given the initial condition {T_{R} (0), U_{H}(0)}, the above system of ODEs describe the composition or the structure of the populations of CAR T cells and target cells at any later time.

##### Parameter estimation

The affinity parameters (*k*_{on}, *k*_{off}) are usually measured in in vitro using experiments such as surface plasmon resonance (Ghorashian et al, 2019) or titrations using flow cytometry (Sharma et al, 2020). The single-cell abundances of CARs and their cognate ligands are measured in quantitative flow cytometry experiments (Hernandez-Lopez et al, 2021) which can be used to estimate {T_{R} (0), U_{H}(0)}. The distributions of R are modeled as log-normal distributions where the mean (μ_{R} (0)) and the SD (σ_{R} (0)) for the distributions at t = 0 are estimated. The signaling parameter k_{p}, the proportionality constants, λ_{c} and ρ_{c}, can be estimated from cytotoxicity assays where CAR T cells and target cells are co-cultured for different CAR T cell and target cell ratios (or E:T ratios) and the fraction of lysed target cells is measured after few days (e.g., 3 d). We used data available from Hernandez-Lopez et al (2021) for constitutive and synthetic notch (synNotch)-CAR T cells to estimate the parameters in our model. Further details are provided in the Materials and Methods section.

##### Pareto optimization

There is a trade-off between maximizing lysis of cancer cells and minimizing lysis of healthy cells by CAR T cells. For example, CAR T cells that lyse target cells expressing a small number of cognate ligands can lyse cancer cells efficiently but can also produce large off-target killing of healthy cells. We set up a multi-objective optimization problem to systematically explore the space of optimal parameter values in our PASCAR model and calculate the Pareto front, which represents the set of optimal parameter values where any objective cannot be optimized further without choosing less than optimal values for the other objectives. We set up a two-objective optimization problem where the percentage of lysed healthy cells and the inverse of the percentage of the lysed cancer cells at a fixed time (e.g., 5 d) post co-incubation of the CAR T cells and the target cells are minimized simultaneously. We consider a mixture of healthy and cancer cells where healthy and cancer cells express on average 10^{4.5} and 10^{6.5} human epidermal growth factor receptor 2 (HER2) molecules/cell (Hernandez-Lopez et al, 2021). The CAR T cells are introduced at t = 0 and interact with the target cells following the PASCAR model, and after a fixed time interval τ, we evaluate the total number of lysed healthy and cancer cells. The parameters in the PASCAR model for constitutive and synNotch-CAR T cells are varied to evaluate the Pareto fronts. Further details of the calculation are provided in the Materials and Methods section.

## Results

### Model with kinetic proofreading captures lysis of target cells by constitutive CAR T cells

We evaluated the capability of Model NKP and Model KP to describe the response of constitutive CAR T cells expressing anti-HER2 single-chain antibody (scFv) against target cells expressing HER-2 in vitro. We fitted both the models to the percentage lysis data obtained from cytotoxicity assays of target cells co-cultured with constitutive CAR T cells and to abundances of CARs in the T cells in the co-culture assayed after 3 d (Hernandez-Lopez et al, 2021). We converted the binding rate k_{on} and the dissociation constant K_{D} from the units of nM to molecules/cell. The high (K_{D} = 17.6 nM, k_{off} = 9.0 × 10^{−5} s^{−1}, k_{on} = 5.1 × 10^{3} M^{−1} s^{−1}) and low (K_{D} = 210 nM, k_{off} = 6.8 × 10^{−4} s^{−1}, k_{on} = 3.2 × 10^{4} M^{−1} s^{−1}) affinitiy CARs used in the experiments by Hernandez-Lopez et al (2021) are converted to the units of K_{D} = 2.39 molecules/cell, k_{off} = 9.0 × 10^{−5} s^{−1}, k_{on} = 3.76 × 10^{−5} (molecules/cell)^{−1} s^{−1}, and K_{D} = 28.47 molecules/cell, k_{off} = 6.8 × 10^{−4} s^{−1}, k_{on} = 2.38 × 10^{−5} (molecules/cell)^{−1} s^{−1}, respectively. The details of the conversion are shown in the Materials and Methods section. We also evaluated the effect of membrane diffusion on estimation of k_{on}, k_{off}, and K_{D} in Equation (1) to find that the values are changed by a small amount (≤20%) from their well-mixed counterpart (Supplemental Data 2) which induces negligible changes in the values of C_{N} (Fig S1). Therefore, we used the rate constants without considering diffusion for simplicity. We found that Model NKP is unable to describe the increase in percentage lysis (Fig 2A) as the affinity of the CAR increases from high (K_{D} = 17.6 nM, k_{off} = 9.0 × 10^{−5} s^{−1}) to low (K_{D} = 210 nM, k_{off} = 6.8 × 10^{−4} s^{−1}), though the model fits the means and the variances of CAR abundances in CAR T cells at t = 3 d reasonably well (Figs S2A and B). This is because the abundances of CAR–HER2 complexes formed in Model NKP for the high- and low-affinity CARs are roughly similar (Fig 2B) which produces almost equal rates of lysis and proliferation for the CAR T cells bearing the high- and low-affinity CARs. Next, we fitted Model KP to the same percentage lysis and CAR expression data which successfully captured the increase in the lysis by the CAR T cells expressing high-affinity CARs (Fig 2C). In the model, the increase in percentage lysis with increasing HER2 density closely follows the increase in C_{N} (Fig 2D) with the HER2 density. The model also reasonably fitted the means and variances of the CARs expressed at t = 3 d (Figs 2E and F). In addition, we computed the Akaike Information Criterion (AIC) for the fits for the KP and the NKP models using AIC = n [*ln* (SSR/n) + *ln* (2π) + 1] + 2*k* where n (=36) is the number of data points, *k* is the number of model parameters, and SSR is the sum of squared residuals. We used the cost function in Equation (5) to calculate the SSR and *k* = 4 and 6 for the NKP and the KP models, respectively. The AIC values for the NKP (AIC = 135.66) and the KP (AIC = 112.05) models show that the KP model is favored substantially (ΔAIC ≫ 2 [Burnham et al, 2011]) over the NKP model. The estimated parameters (Table 1) show that inclusion of kinetic proofreading in the signaling kinetics with and active complex formed at N ≈ 7 steps can separate the CAR T cell responses between high and low-affinity CARs for the same HER2 concentrations (Fig 2D). We further used Model KP to describe the percentage lysis of target cells when constitutive CAR T cells expressed high affinity (K_{D} = 17.6 nM, k_{off} = 9.0 × 10^{−5} s^{−1}) CARs at higher and lower abundances than that considered above. We fitted the percentage lysis for the higher and lower CAR abundances for the cytotoxicity assay performed at the 1:0.35 E:T ratio to estimate the distributions of CAR abundances at t = 0 and kept all other parameters fixed at the values estimated (Table 1) in previously from low and high-affinity CARs. Next, we predicted the percentage lysis for cytotoxic assays of at E:T ratio 1:1 for target cells expressing different HER2 abundances (Fig 2G). The predictions agreed with data reasonably (R^{2} ≈ 0.9) (Fig 2H); however, the model systematically underpredicted the percentage lysis by ∼20% at higher HER2 abundances. This could point to parameters pertaining to signaling kinetics affected by CAR abundances due differences in receptor clustering giving rise to immunological synapse formation at different CAR abundances (see the Discussion section).

### Supplemental Data 2.

Effect of diffusion on the receptor-ligands binding/unbinding rates.[LSA-2023-02171_Supplemental_Data_2.pdf]

Our PASCAR model allows us to analyze how CAR T cell subpopulations expressing different CAR abundances respond to target cells. We used Model KP with best fit parameters (Table 1) to carry out the analysis. The ability of a CAR T cell subpopulation expressing a specific CAR abundance to lyse target cells will depend on (i) the affinity (k_{on}, k_{off}) of CARs for HER2 and the strength of the ensuing signaling in the CAR T cell, and (ii) the number of member CAR T cells in the subpopulation. Therefore, a CAR T cell subpopulation expressing higher CAR abundances but containing a small number of member cells could lyse target cells at a lower rate than a subpopulation expressing lower CAR abundances with a larger number of member cells. We quantified the rate of lysis of target cells expressing abundances of HER2 antigens of magnitude *H* by a CAR T cell subpopulation expressing *R* number of CARs with *T*_{R} number of member cells by *H* abundance (or *T*_{R} because of proliferation (Figs 3A and S3A–D). The CAR T cell subpopulations with intermediate range of CAR abundances ∼(500–2,500 molecules/cell at t = 0) induce a larger rate of lysis of target cells compared with the subpopulations with higher or lower CAR expressions (Fig 3A). This is because of the following reason. The size of the subpopulation T_{R} increases as R increases from small to (∼10 molecules/cell) to intermediate (∼500 molecules/cell) and then starts decreasing as R increases further to larger values (∼2,500 molecules/cell) (Fig S4A). In contrast, λ_{RH} (∝ R, as K_{D} ≪ R or K_{D} ≪ H) increases monotonically with R (Fig S4B). However, the decrease in the subpopulation size (>10-fold) outweighs the increase (∼fivefold) in λ_{RH} as R increases from intermediate to larger values resulting in the decrease of _{on}, k_{off}) of CARs for HER2 and the strength of the ensuing signaling in the CAR T cell, and (ii) the number of target cells expressing a particular abundance (H) of the HER2 antigens. The proliferation rate of a CAR T cell subpopulation expressing abundances of magnitude *R* and containing T_{R} number of member cells as they interact with a target cell population of size *U*_{H} expressing abundance of magnitude *H* of HER antigens is quantified by _{RH} ∝ R when CAR and HER2 interact with high affinity, that is, K_{D} ≪ R and K_{D} ≪ H. The initial CAR distribution estimated by our method follows a lognormal distribution, that is, the T cell subpopulation size for intermediate R is larger than that for larger R. However, the increase in the proliferation rate at larger CAR abundances is unable to increase the size of the subpopulation such that the population peaks at these large values of CAR abundances. In addition, the proliferation rate of the CAR T cells decreases with time as the number of target cells decreases over time because of the lysis of the target cells by the CAR T cells (Figs 3B, S3, and S4). Therefore, the peak of the population still remains at intermediate values of R values at later times.

### Model with kinetic proofreading describes lysis of target cells by synNotch-CAR T cells

We applied Model KP to describe the cytotoxic response and proliferation of synNotch-CAR T cells when the CAR T cells were co-cultured with target cells in vitro. The mean CAR abundance of the synNotch-CAR T cells in Model KP was assumed to be proportional to a Hill function of the mean HER2 abundance of the target cells (details in the Materials and Methods section). We reasoned that CAR T cell signaling and activation in constitutive and synNotch-CAR T cells should involve the same physiochemical processes; therefore, we fitted the percentage lysis data (Hernandez-Lopez et al, 2021) and the CAR expression data (Hernandez-Lopez et al, 2021) for constitutive and synNotch-CAR T cells simultaneously with Model KP. Model KP fits the percentage lysis and the means and variances of CAR and HER2 abundances reasonably well (Fig 4A–C). The estimated parameters (Table 1) show a signaling cascade of N ≈ 7 best fit the data. This result is consistent with recent experiments (Britain et al, 2022) where clustering of LAT in T cell signaling is achieved in N = 7.8 ± 1.1 kinetic proofreading reaction steps. The estimated value of the phosphorylation rate, k_{p} (≈0.007 s^{−1}) is larger than the ligand unbinding rate k_{off} (≈10^{−4} s^{−1}). The estimated Hill function parameters (n_{H} ≈ 4, K_{H} ≈ 2 × 10^{5} molecules/cell) imply that CAR expressions are induced sharply as HER2 abundances increase past 2 × 10^{5} molecules/cell giving rise to the almost binary cytotoxic response (off/on) against healthy cells (∼10^{4.5} HER2 molecules/cell) or tumor cells (>10^{6.5} HER2 molecules/cell). Next, we tested the ability of Model KP to predict synNotch-CAR T cell response for experiments (Hernandez-Lopez et al, 2021) not included in training the model. We used Model KP to predict percentage lysis of target cells by synNotch-CAR T cells at day 3 post-incubation where the synNotch-CAR T cells were co-incubated with different concentrations of target cells not included in model training. The model predictions captured the dependency of the cytotoxic response on the initial effector target ratio observed in experiments (Hernandez-Lopez et al, 2021) reasonably well (Fig 4D and E). We evaluated the correlations between the fitted parameters (Fig S5A) which showed low correlation (|r| < 0.5) between the model parameters; parameter ρ_{0} and N showed the largest dependency (r = −0.49). We also predicted the variation of the percentage lysis for a cytotoxic assay performed at E:T ratio of 1:0.3 for target cells displaying different HER2 abundances interacting syn-Notch CAR T cells expressing the highest affinity CAR (K_{D} = 1.9 nM, k_{on} = 1.2 × 10^{4} M^{−1} s^{−1}, k_{off} = 2.2 × 10^{−5} s^{−1}) which showed excellent agreement (Fig 4F). We also computed the variations in the cost function (Equation (5)) with all pairs of parameters (Fig S5B). The analysis of the synNotch-CAR T cell subpopulations to induce cytotoxicity showed that similar to constitutive CAR T cells, synNotch-CAR T cells with an intermediate range of CAR abundances generated the larger response compared with the subsets with low and high CAR abundances (Figs 4G and S6A–D). The proliferation rates of synNotch-CAR T cell subpopulations increase with increasing CAR expressions (Figs 4H and S6A–D) similar to that of constitutive CAR T cells.

### Optimal design of constitutive and inducible CAR T cells

We performed Pareto optimization in the space of parameters that can be potentially manipulated in experiments. For constitutive CAR T cells, we carried out the analysis for CAR ligand-binding affinity parameters, k_{on} and k_{off}, and, for synNotch-CAR T cells, in addition to the CAR affinity parameters, the threshold (K_{H}) and sharpness (n_{H}) of CAR expression, were also considered. To evaluate the optimal parameters, we set up in silico cytotoxicity assays, where a mixture of 20,000 healthy and tumor target cells with 10,000 constitutive or synNotch-CAR T cells were co-incubated in vitro (Fig 5A). The CAR T cells and the target cells interacted following Model KP set at the best fit parameter values (Table 1) and the percentages of healthy and tumor cells lysed after 5 d of culture were computed. When synNotch-CAR T cells were present, we assumed the CAR expressions in the synNotch-CAR T cells are determined by the HER2 expressions of the tumor cells. The multi-objective optimization was performed in MATLAB where the competing objectives of percentage lysis of healthy cells and 1/(% lysis of tumor cells) were minimized simultaneously. The calculations of Pareto fronts showed that for constitutive CAR T cells, decreasing k_{off} increases lysis of both healthy and tumor cells on a Pareto front at a fixed dissociation rate K_{D} = k_{off}/k_{on}. Decreasing k_{off} increases the lifetime of the CAR–HER2 complexes allowing for the complexes to last longer than the waiting time required for the signaling reactions to generate the active signaling complex required to activate the CAR T cell. Thus, decreasing k_{off} leads to increase in the destruction of tumor and healthy target cells.

However, Pareto fronts at different K_{D} values revealed that increasing K_{D} until a certain limit can increase the lysis of tumor cells, whereas decreasing healthy cell lysis (Fig 5A). For example, for a fixed % lysis of tumor cells (e.g., 66%), increasing K_{D} can decrease % lysis of healthy cells from 60% to 20% (Fig 5A and B). However, increasing K_{D} further starts decreasing lysis of tumor cells as well. This behavior can be explained by the dependency of the abundances of CAR–ligand complexes with K_{D}. Because the average abundances (∼5,000 molecules/cell) of CARs are ∼6 and ∼600 times smaller than that of HER2 on tumor and healthy cells, respectively, most of the CARs in the CAR T cells form complexes with HER2 antigens displayed by the healthy and the tumor cells for high-affinity CARs (K_{D} ≪ 5,000 molecules/cell). As, K_{D} is increased (K_{D} > 5,000 molecules/cell), larger numbers of CAR-HER2 complexes are formed when CAR T cells interact with tumor cells compared with the healthy cells because of the availability of 100 times more HER2 antigens on tumor cells. However, as K_{D} is increased further (K_{D} ≫ 5,000 molecules/cell), the copy numbers of CAR–HER2 complexes decrease substantially for CAR T cells interacting with tumor and healthy cells as the majority of the CARs are unable to form complexes because of the weak affinity of the binding (Fig S7). This results in inefficient lysis of healthy and tumor cells by CAR T cells. The above increase in discrimination between healthy and tumor cells with increasing K_{D} upto a certain limit is qualitatively consistent with experiments (Caruso et al, 2015) with CAR T cells interacting with target cells expressing low (∼30,899 molecules/cell [Caruso et al, 2015]) and high (∼628,265 molecules/cell [Caruso et al, 2015]) EGFR abundances where the CAR was generated from high-affinity cetuximab (K_{D} = 1.9 nM [Talavera et al, 2009]) or low-affinity nimotuzumab (K_{D} = 21 nM [Talavera et al, 2009]). The experiments showed that the CAR T cells with nimotuzumab produced a lower killing of target cells with low EGFR abundances compared with the CAR T cells with cetuximab (Caruso et al, 2015). Both CAR T cells produced similar killing of target cells with high EGFR abundances (Caruso et al, 2015).

Next, we carried out our analysis for synNotch-CAR T cells. Similar to constitutive CAR T cells, the Pareto fronts for synNotch-CAR T cells for fixed n_{H} and K_{H} values showed increased tumor and decreased healthy cell lysis when K_{D} increased for a fixed unbinding rate k_{off} until a particular limit (Fig 5C). For a fixed HER2 affinity (fixed k_{on} k_{off}), increasing n_{H} and decreasing K_{H} increased lysis of tumor and healthy cells (Fig 5D). This is because, increasing n_{H} and decreasing K_{H} produce higher CAR expressions in the synNotch-CAR T cells when they interact with tumor cells, and those CAR T cells induce greater lysis of tumor cells. However, the same CAR T cells become efficiently activated by healthy cells for higher affinity CARs (K_{D} ≪ 5,000 molecules/cell, k_{off} ∼ 10^{−5} s^{−1}) and thus induce increased lysis of healthy cells. Thus, an optimal design of these inducible CAR T cells could be generation of CAR expressions with moderate affinities (K_{D} ∼ 1 μM, k_{off} ∼ 10^{−4} s^{−1}) with higher values of n_{H} and lower values of K_{H}.

We also carried out the above Pareto front analysis for other ratios of healthy and tumor cells (1:4 and 4:1) for constitutive and syn-Notch CAR T cells which showed similar behavior (Figs S8A–I).

## Discussion

We developed a multiscale mechanistic model PASCAR that integrates processes across the scales of molecules of CARs and antigens to the populations of CAR T- and target cells. This multiscale approach is particularly relevant for modeling response of heterogeneous populations of CAR T cells expressing a wide range of CAR abundances against target cells displaying antigen abundances that can vary over 1,000-fold across target cells. We pursued an approximate or coarse-grained modeling approach where many microscopic details of CAR T cell signaling and activation were incorporated implicitly in model parameters. This reduced approach allowed us to quantitatively explore the roles of CAR and HER2 (CAR antigen) abundances, CAR affinities, and the ensuing cell signaling in regulating CAR T cell responses against healthy and tumor cells. The ability of this approach to describe the percentage lysis data and the proliferation of T cell subpopulations with different CAR abundances in cytotoxic assays and to predict results outside model training shows the success of our approach in modeling antitumor responses of constitutive and inducible CAR T cells. We were also able to estimate model parameters reasonably well using flow cytometry measurements of CAR and HER2 expressions and percentage lysis data from cytotoxic assays. Therefore, the PASCAR model combined with data from standard immunoassays can be used to make quantitative predictions for various experimental conditions investigating CAR T cell responses in vitro.

Application of PASCAR to describe cytotoxic responses of constitutive CAR T cells showed the relevance of kinetic proofreading in CAR T cell signaling in discriminating response between high- and low-affinity CARs. Model NKP which did not contain any kinetic proofreading was unable to separate the cytotoxic responses mediated by high- (K_{D} = 17.6 nM, k_{off} = 9.0 × 10^{−5} s^{−1}) and low-affinity (K_{D} = 210 nM, k_{off} = 6.8 × 10^{−4} s^{−1}) CARs. Including kinetic proofreading steps at early time signaling events in time scales of 1/k_{p} (∼2.4 min) in Model KP was able to separate the cytotoxic responses. There are several differences in signaling events induced by TCR and CAR stimulation, for example, LAT is weakly phosphorylated in CAR T cells (Salter et al, 2021), whereas LAT is robustly phosphorylated and forms condensates in TCR signaling within 40 s (McAffee et al, 2022). Therefore, signaling proteins that could correspond to the active complex in Model KP could represent different signaling proteins in TCR and CAR signaling, for example, experiments with light-gated immunoreceptors show the active complex corresponding the N ≈ 7.8 in TCR signaling could represent LAT (Britain et al, 2022); however, the active complex in Model KP for N = 7 is likely to represent a different signaling protein. Formation of immunological synapse where CARs along with the cognate ligands and associated siganling molecules spatially cluster at the interface of CAR T cell and the target cell plays an important role in CAR T cell activation. The degree of spatial clustering of these proteins can influence signal transduction and be affected by the CAR abundances and affinities. This could be a reason why our model underpredicted percentage lysis for the constitutive CAR T cells (Fig 2G) at low and high CAR abundances where PASCAR used the same signaling rates estimated for intermediate CAR abundances. A potential extension of PASCAR could be to include such details of signaling and spatial clustering that can be developed by building on several existing CAR T cell (Rohrs et al, 2018) and signaling models in lymphocytes (Čemerski et al, 2008; Grewal & Das, 2022).

We carried out a multi-objective optimization that minimizes destruction of healthy cells, whereas maximizing the elimination of tumor cells when a select set of model parameters that can be manipulated in experiments were allowed to vary. The calculation of Pareto fronts in our multi-objective optimization showed that intermediate values of CAR affinities led to an increase in tumor cell killing, whereas decreasing healthy cell killing. This finding is supported in previous experiments which showed reduction of CAR affinity reduced healthy cell killing but increased tumor cell killing (Caruso et al, 2015). For inducible CAR T cells, the increase in the sharpness of CAR up-regulation in conjunction with intermediate range of CAR affinities can produce a desirable amount of tumor cell and healthy cell killing.

An exciting extension of PASCAR would be to describe CAR T cell response in vivo. CAR T cells undergo a maturation process over a longer duration (∼weeks) than that considered here to give rise to exhausted (Good et al, 2021) and memory phenotypes (Ghorashian et al, 2019; Wilson et al, 2021 *Preprint*). Cytokines and chemokines in the tumor microenvironment contribute to this maturation process (De Boer et al, 2001; Bell & Gottschalk, 2021). These processes have been modeled for T cells (Carbo et al, 2013) and CAR T cells (Liu et al, 2021) using mechanistic or data-driven models which could be potentially incorporated in the PASCAR model.

### Limitations

The PASCAR model approximates signaling using a minimal model which might not be able to describe engineered CAR adaptors that manipulates the number of ITAMs (Majzner et al, 2020), or engineered T cells where specific signaling proteins such as RasGTPase-activating protein (Carnevale et al, 2022) are knocked out. However, some of these changes could be potentially included in our minimal signaling network by implementing signaling models that have described similar effects in other contexts (Das et al, 2009). The kinetics of the induction and the decay (Roybal et al, 2016) of synNotch CARs are not considered in the current model which could be relevant to describe responses of synNotch CAR T cells as they transit from microenvironments rich in tumor cells to healthy cells in vivo. An extension of the current model to include different compartments with explicit kinetics of synNotch CAR abundances could potentially address this issue.

PASCAR assumes that the CAR abundances in single CD8^{+} T cells do not mix as they divide, that is, daughter cells have the same CAR abundances as the mother cell. However, proteins in human cells (e.g., H1299 non-small cell lung carcinoma cell line) have been observed to mix because of cell division (Sigal et al, 2006) in time scales longer than two cell generations. It is unclear if the CARs follow the similar pattern as the CD8^{+} T cells proliferate. The doubling time scale for the faster proliferating CD8^{+} T cells in our model is ∼1.7 d and the mean doubling time of the CD8^{+} T cell population is ∼3 d; therefore, there will a negligible amount of mixing in the system because of cell proliferation if a similar mixing time scale as in Sigal et al (2006) occurs for the CAR CD8^{+} T cells.

PASCAR also assumes that the target and the T cells are well mixed. In vitro cytotoxicity experiments are carried out in culture wells and for the experiments in Hernandez-Lopez we estimated 99.8% of the T cells were partnered with target cells (details in Supplemental Data 3). However, depending on the number of target and T cells, the number of target cells in the immediate vicinity of a T cell is likely to be varied (details in Supplemental Data 3), therefore, a weighted sum for the target and T cells in Equations (3), (4), and (5) would be more appropriate. We plan to include that in a future study.

### Supplemental Data 3.

Estimation of the probability of contact for chimeric antigen receptor T cells and target cells in cytotoxic assays in vitro.[LSA-2023-02171_Supplemental_Data_3.pdf]

## Materials and Methods

### Solution of the ODEs

We set up a rectangular lattice for the abundances *R* and *H* with lattice constants Δ_{R} and Δ_{H}, respectively. *T*_{R} and *U*_{H} in the ODEs in Equations (3) and (4) denote the numbers of CAR T- and tumor cells with *R* and *H* abundances between *R* to *R* + Δ_{R} and *H* to *H* + Δ_{H}, respectively. The ranges of R, Δ_{R}, H, and Δ_{H}, are chosen based on the ranges ([*R*_{min}, *R*_{max}] and [*H*_{min}, *H*_{max}]) and the distributions of R and H measured in flow cytometry experiments in Hernandez-Lopez et al (2021). Given the parameters, k_{on}, k_{off}, k_{p}, λ_{c}, ρ_{c}, and the initial distributions of *R* and *H*, the nonlinear system of ODEs is solved numerically in MATLAB using the *ode45* function with Runge-Kutta 4 numerical method. The units of *R* and *H* in the ODEs are in (# of molecules)/cell.

### Unit conversion of kinetic rates

The rates k_{on}, k_{off}, and K_{D} are usually provided in the literature in units of (nM)^{−1} s^{−1} (or (μM)^{−1} s^{−1}), s^{−1}, and nM (or μM), respectively. We convert k_{on} and K_{D} rates into units of (# of molecules)^{−1} s^{−1} and (#of molecules) to use in the ODEs where the units for CAR and HER2 abundances are given by (# of molecules)/cell. The unit conversion is carried out as follows. The nM unit is changed to (# of molecules)/(μm)^{3} using 1 nM = 600 × 10^{−3} (# of molecules)/(μm)^{3} = 0.6 (# of molecules)/(μm)^{3}. We assume CAR and HER2 molecules form complexes when these molecules are separated by a distance (*d*) of 2 nm (Helm et al, 1991) or smaller, and the CAR and HER2 molecules interact in the immunological synapse formed at the interface of the interacting CAR T cell and the target cell. The area (*A*) of the synapse region is taken as 1/2 times the area of a T cell (Teague et al, 1993) = ½ × 4π (7/2)^{2} μm^{2}. The number of CAR and HER2 molecules present in the area *A* is assumed to be half of the total numbers of these molecules present in individual cells. The parameters K_{D} and k_{on} in the units of (#of molecules) and (# of molecules)^{−1} s^{−1} are obtained by the following relations: K_{D} [# of molecules] = K_{D} [(# of molecules)/(μm)^{3}]/(*Ad*)[(μm^{3})] and k_{on} [(# of molecules)^{−1} s^{−1}] = k_{on} [(# of molecules)^{−1}s^{−1}/(μm)^{3}]/(*Ad*) [(μm^{3})].

### Parameter estimation

We estimated model parameters k_{p}, λ_{c}, ρ_{c}, and the mean (μ_{R} (0)) and the variance (σ_{R} (0)) of CAR abundances at t = 0 for constitutive and synNotch CAR T cells by fitting percentage lysis data and means and variances of CAR T cells obtained at day 3 post-incubation with target cells in cytotoxic assays. For synNotch-CAR T cells, additional parameters describing up-regulation of CARs because of binding of synNotch receptors with HER2 ligands on target cells were estimated. We minimized the cost function (Wu et al, 2022 *Preprint*) below using *levenberg-marquardt* algorithm in MATLAB._{R} and σ_{R} denote the mean and the SD of the CAR abundances at day 3 post incubation. *nlparci* function in Matlab using the residuals and covariance matrices given by the *nlinfit* in Matlab.

We provide specific details of parameter estimation for constitutive CAR and synNotch-CAR T cells below. (1) *Constitutive CAR T cells*: The data are obtained from Hernandez-Lopez et al (2021) where human CD8^{+} T cells were engineered to express CARs constitutively that bind HER2 with high (K_{D} = 17.6 nM) and low (K_{D} = 210 nM) affinities. In their experiments, human leukemia K562 cell lines were engineered to express five different average concentrations (10^{4.2}, 10^{5.2}, 10^{5.7}, 10^{6.2}, 10^{6.7} molecules/cell) of HER2 molecules which were used as target cells in cytotoxic assays. We fitted (*nlinfit* function in MATLAB) flow cytometry data (Fig S9, and Fig 1C in Hernandez-Lopez et al [2021]) with log-normal distributions to estimate means and variances of HER2 in target cells used in our ODE models. Similarly, the distributions of CAR abundances at day 3 post co-incubation were obtained by fitting the flow cytometry data (Fig S1C in Hernandez-Lopez et al [2021]) for constitutive CARs with log-normal distributions. We assumed the same CAR distributions for high- and low-affinity CARs. Means and variances for the CAR abundances at day 3 in the co-culture experiments were calculated from the estimated log-normal distributions. (2) *synNotch-CAR T cells*: Hernandez-Lopez et al (2021) developed tunable CAR T cells by engineering a synNotch receptor in CD8^{+} T cells. The synNotch receptor binds with HER2 on target cells with low affinity (210 nM) and acts as a high-density antigen filter for inducing CAR expressions in the T cells. The generation of CARs by the synNotch circuit was modeled implicitly. We assume that the changes in the CAR expression in the syn-Notch CAR T cells occur at a faster time scale than that of target cell lysis and T cell proliferation. Thus, the mean CAR abundance (*μ*_{R} (0)) in synNotch-CAR T cells at the start of a co-culture experiment is assumed to be a Hill function of the mean HER2 abundance (μ_{H}) of the target cells given by _{H} is the Hill coefficient and for n_{H} ≥ 2, μ_{R} (0) changes in a switch like all or none fashion as μ_{H} increases beyond the threshold K_{H}. The parameters μ_{s}, n_{H}, and K_{H} were estimated by fitting the percentage lysis (Fig 2A in Hernandez-Lopez et al [2021]) and the CAR expression data (Fig S1C in Hernandez-Lopez et al [2021]) obtained at 3 d after co-incubating target cells and CAR T cells in cytotoxic assays. The estimated distributions of HER2 in target cells described for constitutive CAR T cells were used for modeling experiments with synNotch-CAR T cells as well.

### Pareto optimization

A target cell population of 20,000 cells composed of a 1:1 mixture of healthy and tumor cells was taken at t = 0. The healthy and tumor cells expressed HER2 molecules following a linear superposition of two lognormal distributions where mean values of the HER2 molecules were set to 10^{4.5} molecules/cell and 10^{6.9} molecules/cell for the healthy and tumor cells, respectively. The SDs of each lognormal distribution were set to 0.3 to avoid any substantial overlap between the distributions (Fig S10). In our in silico cytotoxicity assays, we considered the above 20,000 target cells were mixed with 10,000 constitutive or synNotch CAR T cells at t = 0 expressing a high-affinity CAR (K_{D} = 17.6 nM, k_{off} = 9 × 10^{−5} s^{−1}). The distribution of CARs for the constitutive CAR T cells was constructed following lognormal distributions with the parameters estimated in Table 1. For the synNotch CAR T cells, we assumed all the 10,000 CAR T cells expressed CARs following a lognormal distribution with a mean value = *μ*_{H} = 10^{6.9} molecules/cell. The values of n_{H} and K_{H} were set to the values estimated in Table 1 or varied for the Pareto front calculations. We used *gamultiobj* routine in Matlab to compute the Pareto front by optimizing two conflicting objective functions: f_{1} = % lysis of healthy cells at day 5, and f_{2} = 1/(% lysis of the tumor cells) at day 5. The Pareto fronts were calculated for fixed K_{D} values where k_{off} were varied. These calculations (Fig 5C) for the synNotch CAR T cells fixed the values of n_{H} and K_{H} at the values shown in Table 1. For these calculations, we varied f_{1} and f_{2} as functions of k_{off} within ranges (1.0 × 10^{−5}, 1.0 × 10^{−2}) s^{−1} such that k_{on} = k_{off}/K_{D} for different values of K_{D} fixed at (17.6, 1.76 × 10^{5}, 1.76 × 10^{6}, 1.76 × 10^{7}, 1.76 × 10^{8}) nM with all the other parameters fixed at the estimated values given in Table 1 for constitutive and synNotch CAR T cells.

Another set of Pareto fronts were calculated for synNotch CAR T cells (Fig 5D) where n_{H} and K_{H} values were varied and all the other parameters including K_{D} and k_{off} were fixed. In this case, f_{1} and f_{2} varied with n_{H} and K_{H} within ranges (1, 8) and (1 × 10^{3}, 5 × 10^{7}), respectively, with k_{off} and K_{D} values fixed at (9.0 × 10^{−5} s^{−1}, 17.6 nM), (5.0 × 10^{−4} s^{−1}, 3.7 × 10^{4} nM), (5.0 × 10^{−4} s^{−1}, 1.2 × 10^{5} nM), (7.0 × 10^{−4} s^{−1}, 1.0 × 10^{6} nM), and (1.5 × 10^{−3} s^{−1}, 1.1 × 10^{6} nM).

We computed the Pearson’s correlation coefficients between the parameter values as follows. We estimated the model parameters for the SynNotch CAR T cells in 10,000 MCMC simulations and sampled the % lysis data with replacement. We then computed the mean % lysis and the SD (Fig S5). The MATLAB function *corrcoef* was used to compute the Pearson’s correlation coefficients and the *P*-values.

## Data Availability

The MATLAB codes and the data used are available at the GitHub link https://github.com/Harshana4532/CART_Project_Matlab_03Jan2023.

## Acknowledgements

This work is supported by the Research Institute at the Nationwide Children’s Hospital. We thank members of Das laboratory for discussions and feedback. We also thank Dr. Rogelio Hernandez-Lopez and Dr. Wendell Lim for providing clarification for the data shown in Hernandez-Lopez et al (2021).

### Author Contributions

H Rajakaruna: data curation, software, formal analysis, investigation, visualization, methodology, and writing—original draft, review, and editing.

M Desai: data curation, software, investigation, methodology, and writing—original draft.

J Das: conceptualization, resources, formal analysis, supervision, funding acquisition, investigation, methodology, project administration, and writing—original draft, review, and editing.

### Conflict of Interest Statement

The authors declare that they have no conflict of interest.

- Received May 18, 2023.
- Revision received July 8, 2023.
- Accepted July 12, 2023.

- © 2023 Rajakaruna et al.

This article is available under a Creative Commons License (Attribution 4.0 International, as described at https://creativecommons.org/licenses/by/4.0/).