PASCAR: a multiscale framework to explore the design space of constitutive and inducible CAR T 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.

Figure 1. Multiscale protein abundance structured population dynamic model for CAR T cells (PASCAR) model.

Single chimeric antigen receptor (CAR) T cells interact with single target cells in the PASCAR model. The strength of CAR T cell signaling depends on the abundance of the CAR–human epidermal growth factor receptor 2 (HER2) complex (C₀) in Model NKP or on the abundance of an active complex (C_N) formed because of N number of chemical modifications in the CAR–HER2 complex in Model KP. The abundances of C₀ (Model NKP) or C_N (Model KP) depend on the abundances of CAR (R) and HER2 (H) in single CAR T cell and the target cell, respectively. The rate of lysis (λ_RH) of target cells because of the cytotoxic response or the proliferation rate (ρ_RH) of CAR T cells depends on the abundances of C₀ (Model NKP) or C_N (Model KP) in single CAR T cells. The lysis and the proliferation rates are used to describe the kinetics of populations of CAR T cells and target cells.

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₀. The abundance (C₀) 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 $C_{0} (R, H, K_{D}) = \frac{1}{2} (R + H + K_{D}) (1 - \sqrt{1 - \frac{4 R H}{{(R + H + K_{D})}^{2}}})$

We assume that the rate of CAR T cell proliferation ρ_RH and the rate of target cell lysis λ_RH are given by $ρ_{R H} = ρ_{c} C_{0} (R, H, K_{D})$ $λ_{R H} = λ_{c} C_{0} (R, H, K_{D})$ where ρ_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₀^(KP), k_p, and k_off is given by, $C_{N} (C_{0}^{(K P)}, k_{o f f}, k_{p}) = C_{0}^{(K P)} {β (\frac{k_{p}}{k_{o f f} + k_{p}})}^{N}$ (1)where $β = 1 + \frac{k_{p}}{k_{o f f}} .$ C₀^(KP) is given by, $C_{0}^{(K P)} (R, H, K_{D}, k_{p}, k_{o n}, β) = \frac{1}{2 β} (R + H + K_{D}) (1 - \sqrt{1 - \frac{4 R H}{{(R + H + K_{D})}^{2}}})$ (2)

The derivations of Equations (1) and (2) are shown in Supplemental Data 1. The proliferation and lysis rates in Model KP are given by $ρ_{R H} = ρ_{c} C_{N} (C_{0}^{(K P)}, k_{o f f}, k_{p}, N)$ $λ_{R H} = λ_{c} C_{N} (C_{0}^{(K P)}, k_{o f f}, k_{p}, N)$

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: $\frac{d U_{H}}{d t} = r U_{H} - \sum_{R^{'} = R_{\min}}^{R_{\max}} λ_{R^{'} H} T_{R^{'}} U_{H}$ (3) $\frac{d T_{R}}{d t} = \sum_{H^{'} = H_{\min}}^{H_{\max}} ρ_{R H^{'}} U_{H^{'}} T_{R} - δ T_{R}$ (4)

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⁻⁵ s⁻¹, k_on = 5.1 × 10³ M⁻¹ s⁻¹) and low (K_D = 210 nM, k_off = 6.8 × 10⁻⁴ s⁻¹, k_on = 3.2 × 10⁴ M⁻¹ s⁻¹) 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⁻⁵ s⁻¹, k_on = 3.76 × 10⁻⁵ (molecules/cell)⁻¹ s⁻¹, and K_D = 28.47 molecules/cell, k_off = 6.8 × 10⁻⁴ s⁻¹, k_on = 2.38 × 10⁻⁵ (molecules/cell)⁻¹ s⁻¹, 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⁻⁵ s⁻¹) to low (K_D = 210 nM, k_off = 6.8 × 10⁻⁴ s⁻¹), 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] + 2k 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⁻⁵ s⁻¹) 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² ≈ 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).

Figure 2. Protein abundance structured population dynamic model for CAR T cells modeling of cytotoxic and proliferation responses of constitutive chimeric antigen receptor (CAR) T cell against target cells.

(A) Shows fits for Model NKP to the percentage lysis of target cells 3 d after 10,000 constitutive CAR T cells were incubated with 20,000 target cells. The data and the fits are shown for five different human epidermal growth factor receptor 2 (HER2) expressions with average HER2 abundances at 10^4.7, 10^5.2, 10^5.7, 10^6.2, 10^6.9 molecules/cell. The constitutive CAR T cells express either high affinity (K_D = 17.6 nM, k_off = 9.0 × 10⁻⁵ s⁻¹) or low affinity (K_D = 210 nM, k_off = 6.8 × 10⁻⁴ s⁻¹) CARs. The values of percentage lysis at HER2 abundances that are in between the values mentioned above were calculated in the model by interpolating the means and variances of the HER2 distributions (see Figs S11A–D and S12A and B for details). (A, B) Variation of the CAR ligand abundance (or C₀) with the copy numbers of HER2 ligands expressed on target cells for the high- and low-affinity CARs in (A). (A, C) Shows fits for Model KP to the percentage lysis data described in (A). (A, D) Variation of the active CAR–ligand complex (or C_N) with the copy numbers of HER2 ligands expressed on target cells for the high- and low-affinity CARs in (A). (E) Comparison of the measured means and the variances of CAR abundances at day 3 post co-incubation with the fits obtained from Model KP. The mean HER2 abundances of the target cells used for modeling the co-culture experiments are shown along the x-axis. The means and variances for CAR abundances obtained from Hernandez-Lopez et al (2021) are shown with the green bar. (F) Comparison of the model fits for the percentage lysis, and means and variances of distributions of CAR abundances at day 3 post co-incubation with their experimental counterparts. The variables are made nondimensional by scaling the variables by their respective SDs calculated using the measured values. The goodness of fit is quantified by the correlation coefficient R² which shows an excellent agreement (R² = 0.97) between the model and the data. (G) Shows fits (solid lines) for Model KP to the percentage lysis data at day 3 for co-culture experiments at E:T = 1:0.35 (20,000 target cells and 7,000 CAR T cells) with constitutive CAR T cells expressing high-affinity CAR (K_D = 1.9 nM, k_on = 1.2 × 10⁴ M⁻¹ s⁻¹, k_off = 2.2 × 10⁻⁵ s⁻¹) at high (WT) and low (+degron) abundances (Hernandez-Lopez et al, 2021) with HER2 expression distributions given for target cells at different mean HER2 levels. The predictions (dashed lines) for percentage lysis at day 3 generated from Model KP for co-culture experiments at E:T = 1:1 (20,000 target cells and 20,000 CAR T cells) are compared with available measurements in Hernandez-Lopez et al (2021) corresponding percentage lysis co-culture are constitutive model fitted (think line) to day 3 CAR expression data, estimating only the initial CAR (high affinity) distribution (mean and variance) for high (WT)- and low (+degoron)-expression CARs, with all other parameters best estimated by the earlier models fitted to high- and low-affinity CAR for constitutive and synthetic notch CAR. (G, H) Shows goodness of the fit for Model KP with the data (percentage lysis, mean, and variances of the CAR abundances at day 3) for the co-culture experiments in (G) for E:T = 1:0.35. The comparison is shown for the nondimensional variables where they are scaled by their SDs calculated using the measured values.

Figure S2. Fits to chimeric antigen receptor (CAR) expression and percentage lysis data in constitutive CAR T cells for Model NKP.

(A) Comparison of the model fits for the percentage lysis, and means and variances of distributions of CAR abundances at day 3 post co-incubation with their experimental counterparts. The variables are made nondimensional by scaling the variables by their respective SDs calculated using the measured values. The goodness of fit is quantified by the correlation co-efficient R² which shows an excellent agreement (R² = 0.96) between the model and the data. (B) Comparison of the means and the variances of CAR abundances data with NKP model estimations.

Table 1.

List of estimated and fixed parameter values.

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 $c_{R H}^{(l y s i s)} = λ_{R H} T_{R}$ . The kinetics of $c_{R \bar{H}}^{(l y s i s)}$ for target cells expressing mean H abundance (or $\bar{H}$ ) shows an increase with time resulting from the increase in 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 $c_{R \bar{H}}^{(l y s i s)}$ (Fig S4C). The rate of proliferation of CAR T cell subpopulations as they interact with target cell subpopulations expressing specific abundances of HER2 antigens 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 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 $c_{R H}^{(p r o l i f)} = ρ_{R H} U_{H}$ . For a target cell subpopulation expressing mean abundance $\bar{H}$ , $c_{R H}^{(p r o l i f)}$ is given by $c_{R \bar{H}}^{(p r o l i f)} = ρ_{R \bar{H}} U_{\bar{H}}$ . The proliferation rates of CAR T cell subpopulations increase with increasing abundances of CAR expressions (Figs 3B, S3, and S4) almost linearly with R because ρ_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.

Figure 3. Constitutive chimeric antigen receptor (CAR) T cell subpopulations expressing different CAR abundances show different cytotoxic and proliferative responses.

(A) Shows lysis rates of target cells expressing mean human epidermal growth factor receptor 2 abundances (=10^6.2 molecules/cell) where the lysis is mediated by subpopulations of constitutive CAR T cells expressing different CAR abundances. The lysis rates increase with time as the numbers of cells in the CAR T cell subpopulations increase because of cell proliferation. The lysis rates for the CAR T cell subpopulations corresponding to the mean CAR expression (denoted as R-mean) and to the mode of the CAR distribution (denoted as R-peak) at day 3 are marked. The dotted line shows the average of lysis rates across the CAR T cell subpopulations. (B) Shows proliferation rates of CAR T cells as they interact with target cells expressing mean human epidermal growth factor receptor 2 abundances (=10^6.2 molecules/cell).

Figure S3. Constitutive chimeric antigen receptor (CAR) T cell subpopulations expressing different CAR abundances show different cytotoxic and proliferative responses.

(A, C) Shows lysis rates of target cells expressing different mean human epidermal growth factor receptor 2 (HER2) abundances (10^4.5 versus 10^6.2 molecules/cell) by subpopulations of CAR T cells expressing different CAR abundances at different times post co-incubation. Results are shown for CARs with high (17.6 nM, top panel) and low (210.0 nM, bottom panel) affinity towards HER2. The dotted line shows the average of lysis rates across the CAR T cell subpopulations. (A, B, C, D) Shows proliferation rates of CAR T cells as they interact with target cells expressing different mean HER2 abundances (10^4.5 versus 10^6.2 molecules/cell) for the same assays in (A, C).

Figure S4. Kinetics of the CAR T cell populations at mean and peak values of CAR abundances.

(A) Shows the increase in the size of the subpopulations {T_R} with time. The “mean” and “peak” indicate the subpopulations corresponding to the mean chimeric antigen receptor (CAR) abundance and the CAR abundance corresponding to the peak of the CAR distribution, respectively, at t = 0. (B) Shows increase in λ_HR with R for H=H̅. (C) Shows the CAR distributions at days t = 0, 1.0, 2, and 3 d. Note, that the peak of the R distribution remains at intermediate values of R during the kinetics. This behavior can be explained as follows. Consider, $T_{R_{1}} (0)$ and $T_{R_{2}} (0)$ as the sizes of the subpopulations at intermediate (close to the peak) and larger values of R, respectively. If the proliferation rates for the subpopulations are given by g₁ and g₂, respectively, at t = 0, then assuming exponential growth with constant g₁ and g₂, the subpopulations’ sizes at a later time t is given by $T_{R_{1}} (t) = T_{R_{1}} (0) \exp (g_{1} t)$ and ${T_{R_{2}} (t) = T}_{R_{2}} (0) \exp (g_{2} t)$ . Now, to have $\frac{T_{R_{1}} (t)}{T_{R_{2}} (t)} > 1$ as we find in (C), the condition $\frac{T_{R_{1}} (0)}{T_{R_{2}} (0)} > e^{(g_{2} - g_{1}) t}$ , should hold. In the kinetics considered here, $g_{1,2} \propto ρ_{R_{1,2} H} U_{H} \propto R_{1,2}$ . At t = 0 we find $\frac{T_{R_{1}} (0)}{T_{R_{2}} (0)} > 10$ , and g₂−g₁ ≈ 0.075/day in our model which at t = 1d yields $e^{(g_{2} - g_{1}) t} \sim 1.07$ . Thus the above condition is satisfied and the peak remains at the intermediate values of R at t>0. In addition, the rate g₁ and g₂ decreases with time (Fig 3B) which also helps in keeping the above condition intact, and the peak remains at the intermediate values of R for the entire kinetics.

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⁻¹) is larger than the ligand unbinding rate k_off (≈10⁻⁴ s⁻¹). The estimated Hill function parameters (n_H ≈ 4, K_H ≈ 2 × 10⁵ molecules/cell) imply that CAR expressions are induced sharply as HER2 abundances increase past 2 × 10⁵ 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 ρ₀ 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⁴ M⁻¹ s⁻¹, k_off = 2.2 × 10⁻⁵ s⁻¹) 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.

Figure 4. Protein abundance structured population dynamic model for CAR T cell (PASCAR) modeling of cytotoxic and proliferation responses of constitutive and synthetic notch (synNotch) chimeric antigen receptor (CAR) T cell against target cells.

(A) Shows fits for Model KP to the percentage lysis of target cells at 3 d after 10,000 synNotch or constitutive CAR T cells were incubated with 20,000 target cells. The data and the fits are shown for five different human epidermal growth factor receptor 2 (HER2) expressions with average HER2 abundances at 10^4.7, 10^5.2, 10^5.7, 10^6.2, 10^6.9 molecules/cell. The comparisons of the distributions of CAR abundances between the model and the measured values at day 3 are shown in Fig S13. The synNotch and constitutive CAR T cells express either high affinity (K_D = 17.6 nM, k_off = 9.0 × 10⁻⁵ s⁻¹) or low affinity (K_D = 210 nM, k_off = 6.8 × 10⁻⁴ s⁻¹) CARs. (B) Comparison of the model fits for the percentage lysis, and means and variances of distributions of CAR abundances at day 3 post co-incubation with their experimental counterparts. The variables are made nondimensional by scaling the variables by their respective SDs calculated using the measured values. The goodness of fit is quantified by the correlation co-efficient R² which shows an excellent agreement (R² = 0.97) between the model and the data. (C) Comparison of the measured means and the variances of CAR abundances at day 3 post co-incubation with the fits obtained from Model KP. The mean HER2 abundances of the target cells used for modeling the co-culture experiments are shown along the x-axis. The means and variances for CAR abundances obtained from Hernandez-Lopez et al (2021) are shown with green bars. (D) PASCAR model predicted percentage lysis (solid and dashed lines) of target cells at 3 d after 6,000, 12,000, and 15,000 synNotch CAR T cells were incubated with target cells. The data for 10,000 synNotch CAR T cells were used for training the model which are also shown as reference. (E) The predictions in (D) are in good agreement (R² = 0.98) with the data obtained from Hernandez-Lopez et al (2021). (F) PASCAR model well-predicted percentage lysis (solid line) of target cells at 3 d for a cytotoxic assay at E:T = 0.3 where synNotch CAR T cells expressing high-affinity CAR (K_D = 1.9 nM, k_on = 1.2 × 10⁴ M⁻¹ s⁻¹, k_off = 2.2 × 10⁻⁵ s⁻¹) were incubated with target cells. (G) Shows lysis rates of target cells expressing mean HER2 abundances (=10^6.2 molecules/cell) where the lysis is mediated by subpopulations of synNotch CAR T cells expressing different CAR abundances. (H) Shows proliferation rates of synNotch CAR T cells as they interact with target cells expressing mean HER2 abundances (=10^6.2 molecules/cell).

Figure S5. Parameter dependence of the cost function in PASCAR model.

(A) Correlations between PASCAR model parameters and their significance at α = 0.05 (P ≤ 0.05) level estimated by 1,000 MCMC simulations with replacements. Most parameters are not significantly correlated with Pearson’s correlation coefficients remaining less than 0.5, highest being 0.49. (B) Variations of the cost function (Equation (5)) as pairs of parameters indicated along the axes are varied. The other values are kept fixed at the best estimate values.

Figure S6. Protein abundance structured population dynamic model for CAR T cells modeling of cytotoxic and proliferation responses of synthetic notch (synNotch) chimeric antigen receptor (CAR) T cell against target cells.

(A, C) Shows lysis rates of target cells expressing different mean human epidermal growth factor receptor 2 (HER2) abundances (10^4.5 versus 10^6.2 molecules/cell) by subpopulations of synNotch CAR T cells expressing different CAR abundances at different times post co-incubation. Results are shown for CARs with high (17.6 nM, top panel) and low (210.0 nM, bottom panel) affinity towards HER2. The dotted line shows the average of lysis rates across the synNotch CAR T cell subpopulations. (A, B, C, D) Shows proliferation rates of synNotch CAR T cells as they interact with target cells expressing different mean HER2 abundances (10^4.5 versus 10^6.2 molecules/cell) for the same assays in (A, C).

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.

Figure 5. Pareto fronts revealing optimal responses of constitutive and synthetic notch (synNotch) chimeric antigen receptor (CAR) T cells against tumor and healthy cells.

(A) Pareto fronts for constitutive CAR T cells in the plane spanned by % lysis of healthy cells and 1/(% lysis of tumor cells). The Pareto fonts are calculated for different dissociation constants K_D = k_off/k_on where for each K_D value, k_off and k_on = k_off/K_D are varied to obtain the corresponding front. The Pareto fronts are calculated 5 d after the CAR T cells were incubated with a 1:1 mixture of tumor and healthy cells (Fig S10) in silico (details in the main text) with tumor cell number being 20,000. (B) Variations of % lysis of healthy cells with K_D for specific % lysis of tumor cells by constitutive CAR T cells along the Pareto fronts. The % lysis of healthy cells decrease with increasing K_D until a certain value (the end points shown on the graph); increasing the K_D further decreases % lysis of tumor cells as well (not shown on the graph). (C) Pareto fronts for synNotch CAR T cells for different dissociation constants K_D = k_off/k_on where for each K_D value, k_off and k_on = k_off/K_D are varied to obtain the corresponding front. The other parameters K_H and n_H are fixed throughout at 2.42 × 10⁵ molecules/cell and 3.97, respectively. (A) The in silico cytotoxic assay is set up the same way as in (A). (D) Pareto fronts for synNotch CAR T cells for different fixed k_off and K_D values at (9.0 × 10⁻⁵ s⁻¹, 17.6 nM), (5.0 × 10⁻⁴ s⁻¹, 3.7 × 10⁴ nM), (5.0 × 10⁻⁴ s⁻¹, 1.2 × 10⁵ nM), (7.0 × 10⁻⁴ s⁻¹, 1.0 × 10⁶ nM) and (1.5 × 10⁻³ s⁻¹, 1.1 × 10⁶ nM), where k_on = k_off/K_D. K_H and N are varied within (1 × 10³, 5 × 10⁷) molecules/cell and (1, 8) on each Pareto front. The symbol size is proportional to N and the shades filling the symbols are proportional to K_H. (A) The in silico cytotoxic assays are set up the same way as in (A).

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).

Figure S7. Variation of the steady state abundances of C0 and CN with KD.

Shows the variation of C₀ and C_N with K_D where the abundances are computed at mean chimeric antigen receptor = 3,038 and mean human epidermal growth factor receptor 2 abundances for healthy (=10^4.5) and tumor (=10^6.2) cells. The k_off and k_p are set to 9 × 10⁻⁵ s⁻¹ and 0.0072. C₀ and C_N do not change appreciably when chimeric antigen receptor T cells interact with healthy or tumor cells for smaller K_D values (<100 nM).

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⁻⁵ s⁻¹) 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⁻⁴ s⁻¹) 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).

Figure S8. Pareto fronts revealing optimal responses of constitutive and synthetic notch (synNotch) chimeric antigen receptor (CAR) T cells against tumor and healthy cells.

(A, B, C) Pareto fronts for constitutive CAR T cells in the plane spanned by % lysis of healthy cells and 1/(% lysis of tumor cells). The Pareto fonts are calculated for different dissociation constants K_D = k_off/k_on, where for each K_D value, k_off and k_on = k_off/K_D are varied to obtain the corresponding front. The Pareto fronts are calculated 5 d after the CAR T cells were incubated with a 1:4 or 1:1 or 4:1 mixture of healthy and tumor cells in silico (details in the main text) with tumor cell number being 20,000. (D, E, F) Pareto fronts for synNotch CAR T cells for different dissociation constants K_D = k_off/k_on, where for each K_D value, k_off and k_on = k_off/K_D are varied to obtain the corresponding front. The other parameters K_H and n_H are fixed throughout at 2.42 × 10⁵ molecules/cell and 3.97, respectively. (A, B, C) The in silico cytotoxic assay is set up the same way as in (A, B, C). (G, H, I) Pareto fronts for synNotch CAR T cells for different fixed k_off and K_D values at (9.0 × 10⁻⁵ s⁻¹, 17.6 nM), (5.0 × 10⁻⁴ s⁻¹, 3.7 × 10⁴ nM), (5.0 × 10⁻⁴ s⁻¹, 1.2 × 10⁵ nM), (7.0 × 10⁻⁴ s⁻¹, 1.0 × 10⁶ nM), and (1.5 × 10⁻³ s⁻¹, 1.1 × 10⁶ nM), where k_on = k_off/K_D. K_H and N are varied within (1 × 10³, 5 × 10⁷) molecules/cell and (1, 8) on each Pareto front with healthy to tumor cell ratios varying from 1:4, 1:1, 4:1 with tumor cell number being 20,000. The symbol size is proportional to N and the shades filling the symbols are proportional to K_H. (A) The in silico cytotoxic assays are set up the same way as in (A).

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⁻⁵ s⁻¹) and low-affinity (K_D = 210 nM, k_off = 6.8 × 10⁻⁴ s⁻¹) 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.

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)⁻¹ s⁻¹ (or (μM)⁻¹ s⁻¹), s⁻¹, and nM (or μM), respectively. We convert k_on and K_D rates into units of (# of molecules)⁻¹ s⁻¹ 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)³ using 1 nM = 600 × 10⁻³ (# of molecules)/(μm)³ = 0.6 (# of molecules)/(μm)³. 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)² μm². 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)⁻¹ s⁻¹ are obtained by the following relations: K_D [# of molecules] = K_D [(# of molecules)/(μm)³]/(Ad)[(μm³)] and k_on [(# of molecules)⁻¹ s⁻¹] = k_on [(# of molecules)⁻¹s⁻¹/(μm)³]/(Ad) [(μm³)].

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. $C o s t f u n c t i o n = \sum_{a l l e x p e r i m e n t a l c o n d i t i o n s (e . g ., H E R 2 a b u n d a n c e s)} \frac{{({(% l y s i s)}_{e x p t} - {(% l y s i s)}_{m o d e l}))}^{2}}{η_{(% l y s i s)}^{2}} + \frac{{({(μ_{R})}_{e x p t} - {(μ_{R})}_{m o d e l}))}^{2}}{η_{(R)}^{2}} + \frac{{({(σ_{R}^{2})}_{e x p t} - {(σ_{R}^{2})}_{m o d e l}))}^{2}}{η_{(R^{2})}^{2}}$ (5)where μ_R and σ_R denote the mean and the SD of the CAR abundances at day 3 post incubation. $η_{(% l y s i s)}^{2}$ denotes the variance in the % lysis in the experiments which were calculated from Hernandez-Lopez et al (2021) (Fig 2A in that reference). $η_{(R)}^{2}$ and $η_{(R^{2})}^{2}$ denote the variance and the fourth cumulant in CAR expressions, respectively, which were calculated from Hernandez-Lopez et al (2021). The confidence intervals were estimated by the 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 $μ_{R} (0) = μ_{s} \frac{{(μ_{H})}^{n_{H}}}{{(μ_{H})}^{n_{H}} + {(K_{H})}^{n_{H}}}$ . n_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.

Figure S9. Interpolation of the means and variances of human epidermal growth factor receptor 2 (HER2) abundances expressed by target cells.

The points show the measured values and the dashed line shows the interpolation. We used the interpolated means and variances to estimate HER2 abundances following log-normal distributions for HER2 expressions not measured in Hernandez-Lopez et al (2021). The estimated HER2 distributions are used in the protein abundance structured population dynamic model for CAR T cell model to evaluate percentage lysis at shown in Figs 2 and 4 in the main text.

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⁻⁵ s⁻¹). 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 = $μ_{s} \frac{{μ_{H}}^{n_{H}}}{{μ_{H}}^{n_{H}} + {K_{H}}^{n_{H}}}$ where μ_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₁ = % lysis of healthy cells at day 5, and f₂ = 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₁ and f₂ as functions of k_off within ranges (1.0 × 10⁻⁵, 1.0 × 10⁻²) s⁻¹ such that k_on = k_off/K_D for different values of K_D fixed at (17.6, 1.76 × 10⁵, 1.76 × 10⁶, 1.76 × 10⁷, 1.76 × 10⁸) nM with all the other parameters fixed at the estimated values given in Table 1 for constitutive and synNotch CAR T cells.

Figure S10. Simulated initial distribution of human epidermal growth factor receptor 2 densities per target cell in a population of mixture of healthy and tumor cells.

Shows normalized distribution of HER to in a mixture of healthy (10^4.5 molecules/cell on average) and tumor cells (10^6.2 molecules/cell on average) at t = 0. The distribution is a superposition of two log-normal distributions with means at 10^4.5 molecules/cell and 10^6.2 molecules/cell on average and SD of 0.3.

Figure S11. Digitized chimeric antigen receptor (CAR) expression data and discretization of the probability distribution function.

(A) Shows digitized distributions of CAR expressions in synthetic notch CAR T cells obtained from Hernandez-Lopez et al (2021). The units for the abundances are given in terms of AU the fluorescence intensity. (B) Shows the histograms for the CAR distributions used in our modeling. The bin sizes are shown in the figure panels. The AU of fluorescence intensity is converted to the units of molecules/cell based on the calibration of the mean fluorescence intensity to the mean number of CAR abundances performed by Hernandez-Lopez et al (2021). (C) Similar to (A) for constitutive CAR T cells. (D) Similar to (B) for constitutive CAR T cells.

Figure S12. Digitized human epidermal growth factor receptor 2 (HER2) expression data and discretization of the probability distribution function.

(A) Shows digitized distributions of HER2 expressions in the target cells obtained from Hernandez-Lopez et al (2021). The units for the abundances are given in terms of AU the fluorescence intensity. (B) Shows the histograms for the HER2 distributions used in our modeling. The bin sizes are shown in the figure panels. The AU of fluorescence intensity is converted to the units of molecules/cell based on the calibration of the mean fluorescence intensity to the mean number of HER2 abundances performed by Hernandez-Lopez et al (2021).

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₁ and f₂ varied with n_H and K_H within ranges (1, 8) and (1 × 10³, 5 × 10⁷), respectively, with k_off and K_D values fixed at (9.0 × 10⁻⁵ s⁻¹, 17.6 nM), (5.0 × 10⁻⁴ s⁻¹, 3.7 × 10⁴ nM), (5.0 × 10⁻⁴ s⁻¹, 1.2 × 10⁵ nM), (7.0 × 10⁻⁴ s⁻¹, 1.0 × 10⁶ nM), and (1.5 × 10⁻³ s⁻¹, 1.1 × 10⁶ 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.

Legal Disclaimer:

EIN Presswire provides this news content "as is" without warranty of any kind. We do not accept any responsibility or liability for the accuracy, content, images, videos, licenses, completeness, legality, or reliability of the information contained in this article. If you have any complaints or copyright issues related to this article, kindly contact the author above.

You just read:

PASCAR: a multiscale framework to explore the design space of constitutive and inducible CAR T cells

Distribution channels:

EIN Presswire's priority is author transparency. We do our best to weed out false and misleading content. The content above is the sole responsibility of the author who makes it available. If you have any complaints, kindly contact the author above.