A Multiscale Protein Abundance Structured Population Kinetic Model Systematically Explores the Design Space of Constitutive and Inducible CAR-T cells

  1. Harshana Rajakaruna
  2. Milie Desai
  3. Jayajit Das

Reviewed by

Read the full article See related articles

Listed in

Log in to save this article

Abstract

Engineered chimeric antigen receptor (CAR)-T cells are designed to bind to antigens overexpressed on the surface of tumor cells and induce tumor cell lysis. However, healthy cells can express these antigens at lower abundances and can get lysed by CAR-T cells. A wide variety of CAR-T cells have been designed that increase tumor cell elimination while decreasing destruction of healthy cells. However, given the cost and labor-intensive nature of such efforts, a systematic exploration of potential hypotheses becomes limited. To this end, we develop a framework (PASCAR) by combining multiscale population dynamic models and multi-objective optimization approaches with data obtained from published cytometry and cytotoxicity assays to systematically explore design space of constitutive and tunable CAR-T cells. We demonstrate 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 CAR affinities in an intermediate range of dissociation constants (K D ) in constitutive and tunable CAR-T cells can dramatically decrease healthy cell lysis but sustain a high rate of tumor cell killing. In addition, our modeling provides guidance towards optimal tuning of CAR expressions in synNotch CAR T cells. The proposed framework can be extended for other CAR immune cells.

Article activity feed

  1. Review Commons

    Note: This rebuttal was posted by the corresponding author to Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Reply to the reviewers

    We thank both the reviewers for the positive comments, insightful suggestions, and constructive criticisms. We have added new calculations, analyzed new data, added new figures in the main text and the supplementary material, revised the main text and the supplementary text to address the reviewer’s comments. We believe the modifications have improved the quality of the manuscript. Below we append our point-by-point response to the reviewers’ comments.

    *Reviewer #1 (Evidence, reproducibility and clarity (Required)): ** Summary: The paper by Rajakaruna, Desai, and Das develops a multiscale computational model (PASCAR) to study design of features of CAR-T cells. Their ODE-based model captures cell-cell interactions using equilibrium receptor-ligand binding relations and a simplified representation of intracellular signaling. The cellular-scale model feeds into a population model including populations of CAR-T cells, healthy cells, and target cells. The model accounts for CAR-T cells with varying numbers of CARs; it accounts for target/healthy cells with varying numbers of ligands. The authors use published cytometry and cytotoxicity data to parameterize the model, which they show captures experimental trends well when they use a kinetic proofreading model to phenomenologically represent intracellular signaling. They then use an optimization approach to characterize features of the CAR-T cell design space that maximize killing of target cells while minimizing the destruction of healthy cells. Their conclusions are consistent with physical intuition and provide quantitative and generalizable predictions about biophysical parameters (including equilibrium binding constants, kinetic rates, etc.).

    Major comments:

    1. In general, the paper is logically laid out and easy to follow.*

    Thank you for the positive comment.

    *However, some of the underlying mechanisms were not clear to me: *

    (A) - In Figure 3A, it is not clear why the rate of lysis is greater for populations with an intermediate range of CAR abundances. I understand the authors' statement that "this is because the smaller number of CAR-T cells present in [higher or lower CAR expression groups]..." However, it is not clear why intermediate-level CAR-T cell populations are largest, noting that populations with higher receptor levels have larger proliferation rates. Is this simply reflecting the initial size of the populations, or is another mechanism at play? I think the paper would benefit from discussion/quantification of this.*

    We thank the reviewer for the insightful comment.

    We have edited the main text (pages 8-9) and added a supplementary figure (Fig. S4) to address the above comment. The proliferation rate (ρRHUH) for a subpopulation of size TR increases almost linearly with R since ρRH ∝ R when CAR and HER2 interact with high affinity, i.e., KD≪ R and KD≪H. The distribution of CAR abundances at t=0 estimated by our method follows a lognormal distribution, i.e., the CAR T cell subpopulation size at an intermediate value of R is larger than that at larger values of R. However, given the estimated values for ρRHUH (RHUH decreases with time as the number of target cells are lysed by the CAR T cells. Therefore, the peak of the population remains at intermediate values of R values at later times, and both the initial distribution and the relatively smaller values of the estimated proliferation rate are responsible for the behavior.

    *(B) - In Figure 2, C_N plateaus at HER2 density between 10^3 and 10^4. However, % lysis does not plateau until a density closer to 10^5. What is the underlying mechanism? I think the paper would benefit from its exploration.

    Thank you for the comment.

    We have modified the Figure 2 to show the percentage lysis at intermediate values of the HER2 density. The percentage lysis plateaus at similar HER2 concentrations as CN which is expected as the lysis and proliferation rates are proportional to CN. We have modified the text to explain the behavior.

    • It is likely that CAR-T cells would encounter various ratios of target and healthy cells depending on patient, microenvironment, in vitro experimental design, etc. My understanding is that the optimization was done with equal numbers of healthy and target cells. It would be useful to explore how sensitive the optimization is to the ratio of target and healthy cells. This could provide useful design guidance to experimentalists.

    Thank you for the comment.

    We have carried out our Pareto front calculations at different healthy and tumor cell ratios namely, 1:4 and 4:1. The Pareto fronts show similar behavior as that shown for the 1:1 ratio in Figure 4 with small vertical and horizontal shifts in the curves compared to the 1:1. We have shown these additional results in the main text (page 11) and the Figure S8 in the supplementary material.

    • Two assumptions of the population model initially surprised me. However, given the strong performance of the model, they seem to be well justified. Could the authors address the following with brief discussion?

    (A) - The model essentially assumes a well-mixed population of cells, treating lysis and proliferation as second-order "reactions." However, individual cells likely encounter relatively few cells on the time scale of simulations. *

    This is an excellent point.

    We have added a supplementary text (Supplementary Text 3) and the text below in the Discussion section (page 14) to address the above comment.

    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 et al.6 Our estimates show that 99.8% of the T cells were partnered with target cells (details in Supplementary Text 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 Supplementary Text 3), therefore, a weighted sum for the target and T cells in Eq. (3)-(4) would be more appropriate. We plan to include that in a future study.

    (B) - The model doesn't include "mixing" of CAR-T cell populations upon proliferation (i.e., a cell with R receptors can't divide and result in cell with a different number of receptors). Is this justified by the biology?

    We thank the reviewer for raising this excellent point.

    We have added the text below in the Discussion section (page 14) to address the above comment.

    PASCAR assumes that the CAR abundances in single CD8+ T cells do not mix as they divide, i.e., 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 due to 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 days and the mean doubling time of the CD8+ T cell population is ~3 days, therefore, there will a negligible amount of mixing in the system due to cell proliferation if a similar mixing time scale as in Sigal et al. (Sigal et al., 2006) occurs for the CAR CD8+ T cells.

    *Minor comments

    1. A couple of figures appear to be mis-referenced in the paper: Figure 2D (end of paragraph, page 7) should presumably be Figure 2E, and Figure S3 should be S2 (end of page 7).*

    Done.

    • I don't know what "conv." stands for in the figures (I couldn't find the abbreviation in the captions or main text). *

    We have changed all the “conv” abbreviations to “const.” indicating constitutive CAR T cells.

    • Last pargraph of page 9: "However, decreasing K_D further starts decreasing lysis..." Given the previous sentences, this should presumably say "...increasing K_D further..."*

    Done.

    • I was initially confused by the second sentence of the last paragraph of Results. Based on the previous paragraph, I expected the sentence to read "when K_D increased" (not decreased) because of the statement "similar to constitutive CAR-T cells..." However, I think the underlying argument/conclusion is the same.*

    Done.

    • It is not completely clear to me what is being plotted in Figures 2D, 4B, and S1. Where do each of the points come from? Why do there appear to be sets of 4 points when there are more experimental data points in 2C?*

    Thank you for the comment.

    We have edited the figures and the captions to address the comment.

    • What were the values of $\Delta_R$ and $\Delta_H$? (Or, alternatively, how many bins were chosen?) *

    We have shown the ΔH and ΔR values in Figures S10-11 in the supplementary material.

    7. It would be useful to define $\mu$ and $\sigma$ in the main text before they are introduced in Table 1.

    Done.

    *Reviewer #1 (Significance (Required)):

    The paper addresses an important and timely problem, given the strong interest in designing CAR-T cells for cancer therapy. This paper adds to existing computational approaches (clearly summarized in the intro) by introducing a multiscale framework that includes cellular-level properties - like CAR binding affinities, etc. - with a population model that is important for capturing collective behavior of many cells. It is parameterized by experimental data and provides a framework for optimizing cells to maximize target cell killing while minimizing off-target killing of healthy cells. This will be of interest to computational groups in the field, can be extended to incorporate additional biology and/or data, and has the potential to provide useful guidance to experimental design of CAR-T cells. It could be highly impactful if combined with experiments in the future.*

    We appreciate the positive comment.

    *Reviewer #2 (Evidence, reproducibility and clarity (Required)): ** Summary :

    In this study, Rajakaruna et al. propose a mathematical modeling framework called PASCAR to simulate the response of T cells engineered to express chimeric antigen receptors (CAR) against the oncogenic marker HER2 on target cells. The authors propose that the degree of T cell activation is given by a simple CAR occupancy model in the presence or absence of additional kinetic proofreading steps. A population-level ODE model is then used to describe proliferation and death of both target cells and T cells on longer time-scales as a function of the degree of T cell activation.

    The model results are compared with recent experiments by Fernandez-Lopez et al. 2021 (ref. 6). showing that T cells respond to HER2 abundance on target cells in an ultrasensitive manner when using a circuit where a low-affinity synthetic Notch receptor (synNotch) for HER2 controls the expression of a high-affinity CAR for HER2. With this synNotch-CAR circuit, T cells are able to kill tumor cells with high HER2 abundance while sparing healthy cells with 100 fold less HER2.

    The authors first fit model parameters in the situation where T cells constitutively express CARs of either low (Kd=210nM) or high (Kd=17.6nM) affinity and suggest that 7 kinetic proofreading steps are needed to account for the experimental data. For the situation where T cells are endowed with the synNotch-CAR circuit, the authors implicitly account for ultrasensitivity by assuming that CAR abundance is proportionnal to a Hill function of HER2 abundance on target cells. Using this asumption, they fit model parameters to experimental data. The authors use the model to predict the response (%lysis in target cells) for different initial numbers of T cells not used in the fitting procedure and conclude that the PASCAR model can generalize well to unseen situations.

    The authors finally perform pareto optimization to investigate how optimization of lysis of tumor cell with high HER2 abundance and sparing of healthy cells with low HER2 can be performed simultaneously. They conclude by suggesting parameters values of the synNotch-CAR circuit optimizing the discrimination between healthy and tumor cells.

    Major comments :

    1(A) - In Figure 2A, the error bars are not properly drawn : the whiskers are not horizontal, not aligned and seem to have been placed manually. Some points and error bars are found outside plot limits. This is intriguing and suggests that this figure was edited manually. Error bars do not seem to agree with the original experimental data.

    Thank you for pointing this out. We have modified Figure 2A. The extracted experimental data from Hernandez-Lopez et al.6 and the corresponding error bars are also provided in in an excel data file ”Error_bars_%lysis_&Predictions” available at the GitHub link.

    *(B) The lines corresponding to the model results should be drawn using many more points(as compared to experimental results) to better illustrate the model behavior over the range of HER2 abundance values represented. Here, the representation of the model result is misleading (it suggests a linear increase of the %lysis between the 10^3 and 10^5 HER2 molecules/cell). *

    Done. Reviewer 1 also pointed to this in Major comment #1B. We have modified Figures 2A-D,G and Figures 4A,D,F to address this. Please refer to our response to Reviewer #1 for details.

    (C)The number of target and CAR T cells should be modified to correspond to those used in the original experiment (100000 target cells and 15000 CAR T cells for the model versus 20000 target cells and 10000 CAR T cells in the experiments).

    The numbers (20000 of target cells and 10000 CAR T cells) we used in our model are taken from Hernandez-Lopez (see figure caption of Fig. 2 (page 2) in Hernandez-Lopez et al.). In addition, we reached out to Dr. Hernandez-Lopez and Dr. Wendell Lim to confirm the numbers of target cells and CAR T cells used in the experiments.

    2 - An error is found in the expression of CN in the KP model. Indeed, Kd_tilde should be equal to Kd. Given that the authors do not provide the details of their calculation, it is not clear where the error comes from. A potential source of error could be the absence of the dissociation reaction of the CN complex in Fig.1 (see minor comments). The authors should correct the error and re-run the model with the correct expression for CN.

    Thank you for raising this point.

    The expressions shown in Eqs. 1 and 2 in our original manuscript is correct, however, it can be further simplified to get an expression where Kd_tilde is equal to K_d. We have included a detailed derivation of the equations in the supplementary material (Supplementary Text 1).

    3(A) - The conversion scheme used to convert the association rate k_on (and also of the dissociation Kd) from units 1/(nM . s) to units 1/(molecules per cell . s) is not appropriate. Indeed, association rates depend on the diffusion of molecules which greatly differs between soluble (SPR measurements) and membrane-bound molecules. These different diffusion properties should be taken into account for the conversion.

    We thank the reviewer for raising this point.

    We have modified the main text (page 7) and added additional calculations (Supplementary Text 2) and a figure (Figure S1) in the supplementary material to address the comment. Briefly, we evaluated the role of diffusion in modifying the binding (k_on) and unbinding (k_off) rates following the approach developed by Eigen, Bell, Keizer, and others and found that for the rates used in Hernandez-Lopez et al, the changes in the rates are less than 20%, which does not lead to appreciable changes in in CN (Fig. S1). Therefore, including the effect of diffusion will affect the percentage lysis negligibly in this case.

    *(B) Moreover, the conversion formula was found to give a Kd of 3.2 molecules/cell for the high affinity CARs and of 39 molecules/cell for the low affinity CARs. This information is important to interpret the results and should appear in the main text. As such, these Kd estimates are far below the typical number of CARs on T cells (10^3) and of HER2 ligands (from 10^3 to 10^6) on target cells. In this situation, the number of HER2:CAR complexes is entirely determined by the limiting component (HER2 or CAR) and does not depend on Kd (as shown in Fig.2B). Estimates of Kd in the 10^3 molecules/cell range would produce T cell responses that depend on Kd and potentially provide a good agreement between the NKP model and the experimental results shown in Fig.2A. The authors should therefore re-evaluate the performances of the NKP model. *

    We agree with the reviewer regarding the above comment.

    We have moved the conversion of K_D to the unit of molecules/cell. As we point out in our response to the previous comment, including the effect of diffusion in the reactions does not change the K_d appreciably and therefore, the NKP model would be unable to fit the results (% Lysis) shown in Fig 2A.

    4(A) - KP model : The authors make a confusion about the performance of the KP model when stating : "The estimated value of the phosphorylation rate, kp (≈ 0.007 s -1 ) is larger than the ligand unbinding rate koff (≈ 10^-4 s^-1 ) indicating that the kinetic proofreading scheme is active in separating CAR-T cell responses across high and low affinity CARs". KP ligand discrimination is efficient when kp is much smaller than koff which is not the case here. To compensate for this inefficient discrimination, a large number N of proofreading steps is needed.

    Thank you for the comment.

    We have edited the sentence to address the comment (bottom of page 9 in the main text).

    (B) The rate kp=7e-3 s^-1 sets a time-scale 1/kp = 2.5 min. The typical time scale for activation would then be N/kp = 17.5min which is much larger than the time-scale at which the KP mechanism operates during normal TCR signaling following pMHC-TCR engagement. Indeed, during TCR signaling, KP operates within a few tens of seconds as recently illustrated in a study by Mc Affee et al. (https://doi.org/10.1038/s41467-022-35093-9*) showing that LAT condensates appear typically 20 to 30 seconds following TCR engagement. *

    Thank you for the comment and the reference.

    We have edited the discussion to incorporate the above comment (last paragraph in page 12 in the main text).

    5 - The authors should modify Fig.2D to produce a separate graph for each variable ("%lysis", "mean" and "var"). Indeed, differences are potentially hidden by plotting variables with different scales on the same graph. For the comparison of the % lysis between data and model, it would be helpful to have colors and shapes as in Fig.2A and Fig.2C. The authors should also represent the experimental and theoretical distributions of CAR molecules / T cell for the different experimental conditions (different abundance of HER2/target cell).

    Done.

    We have included Figures 2E and 4C in the main text to show comparisons between the means and variances.

    6 - The "Cost" function in the methods section seems to be defined for only one experimental condition (corresponding to a given HER2 abundance). Please indicate that it is the sum over all experimental conditions that is minimized. Author show in Fig.2A the variable (%lysis) used in the cost function. Importantly, the other variables (mean and variance of CAR abundance at day 3) should also be plotted to help the reader appreciate the agreement between model and data.

    Done.

    We have edited the equation in page 15 in the main text to reflect the sum in the cost function.

    We have included figures (Figs. 2F, 2H, 4B) in the main text and the supplementary material (Fig. S2A ) in the supplementary material quantifying comparisons between the data and our model.

    7 - The author should quantify the goodness of the fit and analyze the sensitivity of the model results (the "Cost" function) with respect to parameter values. An investigation of the sensitivity of the model output as of function of all pairs of parameters would certainly highlight the fact that the estimates of parameters kp and N are correlated.

    Done.

    We have calculated correlations between the estimated parameters to address this. The results are shown in the main text (pages 10 and 17) and in Fig. S5 in the supplementary material.

    8 - To model the ultrasensitive synNotch-CAR circuit, the authors assume that CAR expression is a Hill function of HER2 abundance on target cells. The parameters of this Hill function are estimated by fitting both the % lysis and CAR expression at day 3. The authors should evaluate the agreement between model results and experimental values by plotting CAR expression at day 3 (using CAR expression data from Supplementary Figure 1C in ref.6).

    Done.

    We have included the comparisons in Figure S12 in the supplementary material and mentioned those in the figure legend for Fig. 4 in the main text.

    9 - In Figure 4C, for 12000 and 15000 T cells, experimental results show a non-ultrasensitive response (also Supp.Fig.S4E in ref 6) which does not agree well with the model results. Hence, the authors claim that the agreement is good does not seem to be supported. The model should also be tested using experimental data where synNotch and CAR receptors with different affinities are used (in particular Supp.Fig.S4D in ref.6 where a Hill coefficent of 0.6 is estimated for the med-low circuit). It will also be of interest to show that the model can reproduce results from Supp.Fig.S5A in ref.6 where cells with high and low CAR expression are used.

    Thank you for the comment.

    We have calculated the goodness of the fit R2 for our predictions of Fig. S4E in ref. 6 and our R2 = 0.97 suggest excellent agreement with the data. However, we agree that the data for E:T= 0.75 or E:T=0.6 in Fig. S4E in ref. 6 can be interpreted as a gradual increase, however, additional experiments at HER2 abundances 104 molecules/cell probably would be able to help distinguish a ultrasensitive response vs a gradual response at these E:T ratios.

    We have tested model predictions for synNotch CAR T cells for higher affinity CAR (Fig. S4D, top panel, in ref. 6) which shows excellent agreement (Fig. 4F in our revised manuscript). We were unable to generate any prediction for the med-low synNotch circuit (Fig. S4D, bottom panel, in ref. 6) as the CAR expressions for the med-low synNotch circuit are not available in the manuscript. However, we have tested new model predictions (Figs. 2G-H) for constitutive CAR T cells (Fig. S5A in ref. 6) for high and low CAR expressions at different E:T ratios. The agreement between the data and the predictions is reasonable (R2 = 0.90) . Therefore, we believe we have confronted our model with responses generated by several types of CAR T cells and the excellent to reasonable agreements of the PASCAR model with the data provide confidence in the utility of our framework to investigate CAR T cell responses in vitro.

    Minor comments :

    *Figure 1 : left : The schematic of the CAR receptor is inaccurate. CD3z domains should be shown within the intra-cytoplasmic part of the receptor, and not as separate proteins. In addition to the CD3z domain, the 41BB domain should also be represented (see ref 6 by Fernadez-Lopez at al.). The arrow correspond to the dissociation of the complex C_N is missing. Without this reaction in the KP model, the steady state solution leads to C_i = 0 for i in [0, N-1]. All receptor-ligand complexes eventually reach and stay in the C_N state. *

    We thank the reviewer for catching this.

    We have made changes in the schematic figure to include the changes suggested above.

    middle: In the schematic, add labels "lysis of target cell" and "proliferation" next to the corresponding parameters. The arrows should point from the label ("cancer cell" and "CAR-T cell") to the cell not to the other way around.

    Done.

    Right : this schematic is not useful and should either be removed or completed to provide useful information to the reader.

    Done.

    *Figure 2 : A - A reference to the article with the original data should be added in the legend. In line with the original article, "conv." should be replaced by "const." as an abbreviation for "constitutive". B and E - Please use a log scale for the y-axis. Have the authors represented an average of C0 and CN across the population at day 3? If so, that information should appear in the figure legend. *

    Done.

    Figure 3. To help the reader interpret these graphs, the authors should also show T_R(t) and U_H(t) on separate graphs. It would also be useful to show C_N as a function of R for H=10^6.2 both for high and low affinity CAR.

    Done.

    Figure 4. A - the number of target cells and CAR T cells do not correspond to those used in the original experiment (see major comment 1).

    Done. Please check our response to comment #1(B) for Reviewer #2.

    *B - Same comments as for Fig.2D. The legend appears to be incorrect ("% lysis" should be blue open circle and orange open squares). *

    Done.

    *C - Use more points to plot the model results. *

    These figures have all the points plotted. However, some points are overlapping with some other.

    Reviewer #2 (Significance (Required)):

    The modeling framework is minimal but appropriate to describe CAR T cell activation and subsequent proliferation and lysis of target cells. The advantage of the simplicity of the model formulation is to allow a direct interpretation of the impact of the different parameters on the model output.

    Thank you for the positive remark.

    *However, the authors seldom discuss how the behavior of the model is controlled by the parameters and their values. Accordingly, the analysis of the theoretical results needs to be further developed. The authors should also quantify the goodness of fit and analyze the sensitivity of model results to parameter values. This will allow to evaluate how the experimental data constrain model parameters and to compare the performances of the different models. *

    We have provided further explanation for some of the results obtain using our model (please refer to responses to comments # 1A, 3A, and 2B) for Reviewer 1), added a correlation analysis for the estimated parameters (Fig. S5 in the supplementary material), and included goodness of the fit (R2 values) for the fits and the model predictions. We believe these new results address the reviewers’ comments.

    The authors should also assess the biological significance of their results by confronting their parameter estimates to related parameter estimates used in other models of T cell activation. The model should further be tested in other situations with different receptor affinities and cell numbers.

    Thank you for the comment.

    We have now tested new predictions for a highest affinity synN-CAR (Figure 4F) and for increasing and decreasing CAR abundances at different E:T ratios (Figure 2G-H). There are few models of CAR T cells (e.g., Rohrs et al. iScience, 2020) which use more detailed signaling models. Therefore, it will be difficult to compare the estimated values of parameters in our model and with those models. The time scales of signaling can also depend on the specifics of the CAR construct. Therefore, though it will be useful to compare different models and their parameter values developed so far, we believe this is outside the scope of this manuscript. We have included some potential directions for extending our current framework in the Discussion section (pages 12-14). If the editor and the reviewer strongly feel the need for including this in the current manuscript, we can do that.

    This research should be of interest to readers specialized in the field of mathematical modeling of biological systems and in CAR-T cell immunotherapies.

    Thank you for the positive comment.

    *Fields of expertise of the reviewer : biophysics, mathematical modeling of biological systems, molecular mechanisms of T cell signaling and activation. *

  2. Review Commons

    Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Referee #2

    Evidence, reproducibility and clarity

    Summary:

    In this study, Rajakaruna et al. propose a mathematical modeling framework called PASCAR to simulate the response of T cells engineered to express chimeric antigen receptors (CAR) against the oncogenic marker HER2 on target cells. The authors propose that the degree of T cell activation is given by a simple CAR occupancy model in the presence or absence of additional kinetic proofreading steps. A population-level ODE model is then used to describe proliferation and death of both target cells and T cells on longer time-scales as a function of the degree of T cell activation.

    The model results are compared with recent experiments by Fernandez-Lopez et al. 2021 (ref. 6). showing that T cells respond to HER2 abundance on target cells in an ultrasensitive manner when using a circuit where a low-affinity synthetic Notch receptor (synNotch) for HER2 controls the expression of a high-affinity CAR for HER2. With this synNotch-CAR circuit, T cells are able to kill tumor cells with high HER2 abundance while sparing healthy cells with 100 fold less HER2.

    The authors first fit model parameters in the situation where T cells constitutively express CARs of either low (Kd=210nM) or high (Kd=17.6nM) affinity and suggest that 7 kinetic proofreading steps are needed to account for the experimental data. For the situation where T cells are endowed with the synNotch-CAR circuit, the authors implicitly account for ultrasensitivity by assuming that CAR abundance is proportionnal to a Hill function of HER2 abundance on target cells. Using this asumption, they fit model parameters to experimental data. The authors use the model to predict the response (%lysis in target cells) for different initial numbers of T cells not used in the fitting procedure and conclude that the PASCAR model can generalize well to unseen situations.

    The authors finally perform pareto optimization to investigate how optimization of lysis of tumor cell with high HER2 abundance and sparing of healthy cells with low HER2 can be performed simultaneously. They conclude by suggesting parameters values of the synNotch-CAR circuit optimizing the discrimination between healthy and tumor cells.

    Major comments:

    1. In Figure 2A, the error bars are not properly drawn : the whiskers are not horizontal, not aligned and seem to have been placed manually. Some points and error bars are found outside plot limits. This is intriguing and suggests that this figure was edited manually. Error bars do not seem to agree with the original experimental data. The lines corresponding to the model results should be drawn using many more points(as compared to experimental results) to better illustrate the model behavior over the range of HER2 abundance values represented. Here, the representation of the model result is misleading (it suggests a linear increase of the %lysis between the 10^3 and 10^5 HER2 molecules/cell). The number of target and CAR T cells should be modified to correspond to those used in the original experiment (100000 target cells and 15000 CAR T cells for the model versus 20000 target cells and 10000 CAR T cells in the experiments).
    2. An error is found in the expression of CN in the KP model. Indeed, Kd_tilde should be equal to Kd. Given that the authors do not provide the details of their calculation, it is not clear where the error comes from. A potential source of error could be the absence of the dissociation reaction of the CN complex in Fig.1 (see minor comments). The authors should correct the error and re-run the model with the correct expression for CN.
    3. The conversion scheme used to convert the association rate k_on (and also of the dissociation Kd) from units 1/(nM . s) to units 1/(molecules per cell . s) is not appropriate. Indeed, association rates depend on the diffusion of molecules which greatly differs between soluble (SPR measurements) and membrane-bound molecules. These different diffusion properties should be taken into account for the conversion.

    Moreover, the conversion formula was found to give a Kd of 3.2 molecules/cell for the high affinity CARs and of 39 molecules/cell for the low affinity CARs. This information is important to interpret the results and should appear in the main text. As such, these Kd estimates are far below the typical number of CARs on T cells (10^3) and of HER2 ligands (from 10^3 to 10^6) on target cells. In this situation, the number of HER2:CAR complexes is entirely determined by the limiting component (HER2 or CAR) and does not depend on Kd (as shown in Fig.2B). Estimates of Kd in the 10^3 molecules/cell range would produce T cell responses that depend on Kd and potentially provide a good agreement between the NKP model and the experimental results shown in Fig.2A. The authors should therefore re-evaluate the performances of the NKP model.

    1. KP model : The authors make a confusion about the performance of the KP model when stating : "The estimated value of the phosphorylation rate, kp (≈ 0.007 s -1 ) is larger than the ligand unbinding rate koff (≈ 10^-4 s^-1 ) indicating that the kinetic proofreading scheme is active in separating CAR-T cell responses across high and low affinity CARs". KP ligand discrimination is efficient when kp is much smaller than koff which is not the case here. To compensate for this inefficient discrimination, a large number N of proofreading steps is needed.

    The rate kp=7e-3 s^-1 sets a time-scale 1/kp = 2.5 min. The typical time scale for activation would then be N/kp = 17.5min which is much larger than the time-scale at which the KP mechanism operates during normal TCR signaling following pMHC-TCR engagement. Indeed, during TCR signaling, KP operates within a few tens of seconds as recently illustrated in a study by Mc Affee et al. (https://doi.org/10.1038/s41467-022-35093-9) showing that LAT condensates appear typically 20 to 30 seconds following TCR engagement.

    1. The authors should modify Fig.2D to produce a separate graph for each variable ("%lysis", "mean" and "var"). Indeed, differences are potentially hidden by plotting variables with different scales on the same graph. For the comparison of the % lysis between data and model, it would be helpful to have colors and shapes as in Fig.2A and Fig.2C. The authors should also represent the experimental and theoretical distributions of CAR molecules / T cell for the different experimental conditions (different abundance of HER2/target cell).
    2. The "Cost" function in the methods section seems to be defined for only one experimental condition (corresponding to a given HER2 abundance). Please indicate that it is the sum over all experimental conditions that is minimized. Author show in Fig.2A the variable (%lysis) used in the cost function. Importantly, the other variables (mean and variance of CAR abundance at day 3) should also be plotted to help the reader appreciate the agreement between model and data.
    3. The author should quantify the goodness of the fit and analyze the sensitivity of the model results (the "Cost" function) with respect to parameter values. An investigation of the sensitivity of the model output as of function of all pairs of parameters would certainly highlight the fact that the estimates of parameters kp and N are correlated.
    4. To model the ultrasensitive synNotch-CAR circuit, the authors assume that CAR expression is a Hill function of HER2 abundance on target cells. The parameters of this Hill function are estimated by fitting both the % lysis and CAR expression at day 3. The authors should evaluate the agreement between model results and experimental values by plotting CAR expression at day 3 (using CAR expression data from Supplementary Figure 1C in ref.6).
    5. In Figure 4C, for 12000 and 15000 T cells, experimental results show a non-ultrasensitive response (also Supp.Fig.S4E in ref 6) which does not agree well with the model results. Hence, the authors claim that the agreement is good does not seem to be supported. The model should also be tested using experimental data where synNotch and CAR receptors with different affinities are used (in particular Supp.Fig.S4D in ref.6 where a Hill coefficent of 0.6 is estimated for the med-low circuit). It will also be of interest to show that the model can reproduce results from Supp.Fig.S5A in ref.6 where cells with high and low CAR expression are used.

    Minor comments:

    Figure 1 : left : The schematic of the CAR receptor is inaccurate. CD3z domains should be shown within the intra-cytoplasmic part of the receptor, and not as separate proteins. In addition to the CD3z domain, the 41BB domain should also be represented (see ref 6 by Fernadez-Lopez at al.). The arrow correspond to the dissociation of the complex C_N is missing. Without this reaction in the KP model, the steady state solution leads to C_i = 0 for i in [0, N-1]. All receptor-ligand complexes eventually reach and stay in the C_N state.

    middle: In the schematic, add labels "lysis of target cell" and "proliferation" next to the corresponding parameters. The arrows should point from the label ("cancer cell" and "CAR-T cell") to the cell not to the other way around.

    Right : this schematic is not useful and should either be removed or completed to provide useful information to the reader.

    Figure 2 : A - A reference to the article with the original data should be added in the legend. In line with the original article, "conv." should be replaced by "const." as an abbreviation for "constitutive". B and E - Please use a log scale for the y-axis. Have the authors represented an average of C0 and CN across the population at day 3? If so, that information should appear in the figure legend.

    Figure 3. To help the reader interpret these graphs, the authors should also show T_R(t) and U_H(t) on separate graphs. It would also be useful to show C_N as a function of R for H=10^6.2 both for high and low affinity CAR.

    Figure 4. A - the number of target cells and CAR T cells do not correspond to those used in the original experiment (see major comment 1). B - Same comments as for Fig.2D. The legend appears to be incorrect ("% lysis" should be blue open circle and orange open squares). C - Use more points to plot the model results.

    Significance

    The modeling framework is minimal but appropriate to describe CAR T cell activation and subsequent proliferation and lysis of target cells. The advantage of the simplicity of the model formulation is to allow a direct interpretation of the impact of the different parameters on the model output. However, the authors seldom discuss how the behavior of the model is controlled by the parameters and their values. Accordingly, the analysis of the theoretical results needs to be further developed. The authors should also quantify the goodness of fit and analyze the sensitivity of model results to parameter values. This will allow to evaluate how the experimental data constrain model parameters and to compare the performances of the different models. The authors should also assess the biological significance of their results by confronting their parameter estimates to related parameter estimates used in other models of T cell activation. The model should further be tested in other situations with different receptor affinities and cell numbers.

    This research should be of interest to readers specialized in the field of mathematical modeling of biological systems and in CAR-T cell immunotherapies.

    Fields of expertise of the reviewer: biophysics, mathematical modeling of biological systems, molecular mechanisms of T cell signaling and activation.

  3. Review Commons

    Note: This preprint has been reviewed by subject experts for Review Commons. Content has not been altered except for formatting.

    Learn more at Review Commons

    Referee #1

    Evidence, reproducibility and clarity

    Summary: The paper by Rajakaruna, Desai, and Das develops a multiscale computational model (PASCAR) to study design of features of CAR-T cells. Their ODE-based model captures cell-cell interactions using equilibrium receptor-ligand binding relations and a simplified representation of intracellular signaling. The cellular-scale model feeds into a population model including populations of CAR-T cells, healthy cells, and target cells. The model accounts for CAR-T cells with varying numbers of CARs; it accounts for target/healthy cells with varying numbers of ligands. The authors use published cytometry and cytotoxicity data to parameterize the model, which they show captures experimental trends well when they use a kinetic proofreading model to phenomenologically represent intracellular signaling. They then use an optimization approach to characterize features of the CAR-T cell design space that maximize killing of target cells while minimizing the destruction of healthy cells. Their conclusions are consistent with physical intuition and provide quantitative and generalizable predictions about biophysical parameters (including equilibrium binding constants, kinetic rates, etc.).

    Major comments:

    1. In general, the paper is logically laid out and easy to follow. However, some of the underlying mechanisms were not clear to me:
    • In Figure 3A, it is not clear why the rate of lysis is greater for populations with an intermediate range of CAR abundances. I understand the authors' statement that "this is because the smaller number of CAR-T cells present in [higher or lower CAR expression groups]..." However, it is not clear why intermediate-level CAR-T cell populations are largest, noting that populations with higher receptor levels have larger proliferation rates. Is this simply reflecting the initial size of the populations, or is another mechanism at play? I think the paper would benefit from discussion/quantification of this.
    • In Figure 2, C_N plateaus at HER2 density between 10^3 and 10^4. However, % lysis does not plateau until a density closer to 10^5. What is the underlying mechanism? I think the paper would benefit from its exploration.
    1. It is likely that CAR-T cells would encounter various ratios of target and healthy cells depending on patient, microenvironment, in vitro experimental design, etc. My understanding is that the optimization was done with equal numbers of healthy and target cells. It would be useful to explore how sensitive the optimization is to the ratio of target and healthy cells. This could provide useful design guidance to experimentalists.
    2. Two assumptions of the population model initially surprised me. However, given the strong performance of the model, they seem to be well justified. Could the authors address the following with brief discussion?
    • The model essentially assumes a well-mixed population of cells, treating lysis and proliferation as second-order "reactions." However, individual cells likely encounter relatively few cells on the time scale of simulations.
    • The model doesn't include "mixing" of CAR-T cell populations upon proliferation (i.e., a cell with R receptors can't divide and result in cell with a different number of receptors). Is this justified by the biology?

    Minor comments

    1. A couple of figures appear to be mis-referenced in the paper: Figure 2D (end of paragraph, page 7) should presumably be Figure 2E, and Figure S3 should be S2 (end of page 7).
    2. I don't know what "conv." stands for in the figures (I couldn't find the abbreviation in the captions or main text).
    3. Last pargraph of page 9: "However, decreasing K_D further starts decreasing lysis..." Given the previous sentences, this should presumably say "...increasing K_D further..."
    4. I was initially confused by the second sentence of the last paragraph of Results. Based on the previous paragraph, I expected the sentence to read "when K_D increased" (not decreased) because of the statement "similar to constitutive CAR-T cells..." However, I think the underlying argument/conclusion is the same.
    5. It is not completely clear to me what is being plotted in Figures 2D, 4B, and S1. Where do each of the points come from? Why do there appear to be sets of 4 points when there are more experimental data points in 2C?
    6. What were the values of $\Delta_R$ and $\Delta_H$? (Or, alternatively, how many bins were chosen?)
    7. It would be useful to define $\mu$ and $\sigma$ in the main text before they are introduced in Table 1.

    Significance

    The paper addresses an important and timely problem, given the strong interest in designing CAR-T cells for cancer therapy. This paper adds to existing computational approaches (clearly summarized in the intro) by introducing a multiscale framework that includes cellular-level properties - like CAR binding affinities, etc. - with a population model that is important for capturing collective behavior of many cells. It is parameterized by experimental data and provides a framework for optimizing cells to maximize target cell killing while minimizing off-target killing of healthy cells. This will be of interest to computational groups in the field, can be extended to incorporate additional biology and/or data, and has the potential to provide useful guidance to experimental design of CAR-T cells. It could be highly impactful if combined with experiments in the future.

  4. Version published to 10.1101/2023.01.09.523295v1 on bioRxiv