# Mathematical deconvolution of CAR T-cell proliferation and exhaustion from real-time killing assay data

## Abstract

Chimeric antigen receptor (CAR) T-cell therapy has shown promise in the treatment of haematological cancers and is currently being investigated for solid tumours, including high-grade glioma brain tumours. There is a desperate need to quantitatively study the factors that contribute to the efficacy of CAR T-cell therapy in solid tumours. In this work, we use a mathematical model of predator–prey dynamics to explore the kinetics of CAR T-cell killing in glioma: the Chimeric Antigen Receptor T-cell treatment Response in GliOma (CARRGO) model. The model includes rates of cancer cell proliferation, CAR T-cell killing, proliferation, exhaustion, and persistence. We use patient-derived and engineered cancer cell lines with an *in vitro* real-time cell analyser to parametrize the CARRGO model. We observe that CAR T-cell dose correlates inversely with the killing rate and correlates directly with the net rate of proliferation and exhaustion. This suggests that at a lower dose of CAR T-cells, individual T-cells kill more cancer cells but become more exhausted when compared with higher doses. Furthermore, the exhaustion rate was observed to increase significantly with tumour growth rate and was dependent on level of antigen expression. The CARRGO model highlights nonlinear dynamics involved in CAR T-cell therapy and provides novel insights into the kinetics of CAR T-cell killing. The model suggests that CAR T-cell treatment may be tailored to individual tumour characteristics including tumour growth rate and antigen level to maximize therapeutic benefit.

### 1. Statement of significance

We use a mathematical model to deconvolute the nonlinear contributions of chimeric antigen receptor (CAR) T-cell proliferation and exhaustion to predict therapeutic efficacy and dependence on CAR T-cell dose and target antigen levels.

### 2. Introduction

CAR T-cell therapy is a targeted immunotherapy, demonstrating remarkable anti-tumour efficacy, particularly in the treatment of haematologic cancers [1,2]. CAR T-cell therapy is a specific type of immunotherapy where T-cells are genetically modified to recognize a tumour antigen thereby specifically redirecting T-cell cytolytic activity. Inspired by the success of CAR T-cell therapy in liquid tumours, there has been great interest in expanding the use of CAR T-cells for the treatment of solid tumours, such as glioblastoma (GBM), a highly aggressive form of primary brain cancer. Several clinical trials using CAR T-cells to treat GBM have been initiated all over the world [3–6]. At this early stage of clinical development, CAR T-cells offer much promise in solid tumours. However, the diversity of current clinical trials employing varying types of CARs for different solid tumours, target patient populations, preconditioning regimes and cell origins (autologous versus allogeneic) presents a significant challenge in identifying which aspects of a given CAR T-cell treatment protocol are most critical for its effectiveness. An additional critical challenge for CAR T-cell therapy is the potential for transient progression, where the cancer appears to progress before eventually responding to the treatment [7,8].

In order to address these challenges in CAR T-cell therapy for solid tumours, we endeavoured to study the kinetics of CAR T-cell killing with an *in vitro* system and a mathematical model. Mathematical models are useful to describe, quantify and predict multifaceted behaviour of complex systems, such as interactions between cells. A mathematical model is a formalized method to hypothesize systems dynamics, and yield solutions that predict the system's behaviour with a given set of parameters and initial conditions. Mathematical models can be versatile and tested with clinical data which may be obtained *in vivo* from non-invasive imaging [9–11] and the models can be refined when additional information about the system becomes available. Many mathematical models have been developed to understand tumour progression to guide refinement of cancer therapy regimens [12–14].

As CAR T-cell therapy is a newly advanced treatment modality, relatively few studies have used computational modelling to understand and improve this cell-based therapy. Recently, computational models have been developed to investigate cytokine release syndrome for toxicity management [15–17], effect of cytokine release syndrome on CAR T-cell proliferation [18], mechanisms of CAR T-cell activation [19,20], and dosing strategies [21]. However, it remains an open challenge how to use mathematical modelling to study and ultimately predict dynamics of CAR T-cell mediated cancer cell killing with respect to CAR T-cell dose, donor-dependent T-cell differences, cancer cell proliferation, target antigen expression, and how these factors contribute to the overall effectiveness of CAR T-cell therapy.

Based upon our pre-clinical and clinical experience with our well-characterized IL13Rα2-targeted CAR T-cell therapy for recurrent GBM [22,23], we have identified several factors which contribute to the effectiveness of CAR T-cells, namely: rates of proliferation, exhaustion, persistence and target cell killing. To study these various facets of CAR T-cell killing kinetics, we modelled the dynamics between cancer cells and CAR T-cells as a predator–prey system with a mathematical model we call CARRGO: Chimeric Antigen Receptor T-cell treatment Response in GliOma. We use a real-time cell analyser experimental system to estimate parameters of the mathematical model and then apply the model to *in vivo* human data. The long-term aim of this work is to develop a model which could be used to predict and eventually to optimize response to CAR T-cell therapy.

### 3. Methods

The CARRGO mathematical model is a variation on the classic Lotka–Volterra [24,25] predator–prey equations:

*X*represents the density of cancer cells,

*Y*is the density of CAR T-cells,

*ρ*is the net growth rate of cancer cells,

*K*is the cancer cell carrying capacity,

*κ*

_{1}is the killing rate of the CAR T-cells,

*κ*

_{2}is the net rate of proliferation including exhaustion of CAR T-cells when encountered by a cancer cell and

*θ*is the death rate of CAR T-cells. The parameters

*ρ*,

*K*,

*κ*

_{1},

*θ*are constants and assumed to be non-negative except for

*κ*

_{2}which can be either positive or negative (table 1). A positive value of

*κ*

_{2}indicates an increased rate of CAR T-cell proliferation when stimulated by interaction with a cancer cell. A negative value of

*κ*

_{2}indicates exhaustion or limited activation of CAR T-cells resulting from interaction with a cancer cell. Exhaustion and hypoactivation of CAR T-cells are combined into a single value and are not modelled individually.

parameter | description | unit |
---|---|---|

ρ |
cancer cell net growth rate | day^{−1} |

K |
carrying capacity | cell |

κ_{1} |
CAR T-cell killing rate | day^{−1} cell^{−1} |

κ_{2} |
net rate of proliferation and exhaustion of CAR T-cells when stimulated by cancer cells | day^{−1} cell^{−1} |

θ |
CAR T-cell death rate (persistence) | day^{−1} |

We chose to model the net number of cancer cells and simple interactions between cancer cells and CAR T-cells because the output data from the culture system are limited to cell number over time. We therefore are only able to infer dynamics at this scale and dimension (i.e. number of cells and time). Moreover, we performed a system identifiability analysis to demonstrate the parameters of the model can be uniquely determined from the data in this experiment (see electronic supplementary material, methods) [26–29]. Future studies may examine more complex dynamics such as individual cell antigen levels, heterogeneity, resistant and sensitive sub-populations, repeated treatments, etc. with other experimental designs which directly measure these features.

#### 3.1. Model assumptions

The CARRGO model treats cancer cell–CAR T-cell dynamics in this experimental condition as a closed predator–prey system. The model assumes (1) the populations are well mixed, (2) cancer cell growth is limited by space and nutrients (culture media) in the *in vitro* culture system and therefore grow logistically, (3) CAR T-cells kill cancer cells when they interact via the law of mass action, (4) the CAR T-cell killing rate does not explicitly assume a dependence on antigen density, (5) CAR T-cells may be stimulated to proliferate or to undergo loss of effector function—defined as exhaustion—upon contact with a cancer cell [30], and (6) the CAR T-cell death rate is independent of cancer cell density. We chose the logistic growth model for the cancer cell population because the fixed growth rate and carrying capacity parameters were the biological quantities of interest when comparing CAR T-cell killing kinetics across cell lines. Witzel *et al.* compared several sigmoidal growth laws including logistic, Gompertz and Richards, and showed that all these models can be fitted equally well to this form of experimental data [31]. Data supporting our model assumptions are given in electronic supplementary material, figures S1 and S2.

#### 3.2. Dynamical system analysis of the CARRGO model

Closed form solutions cannot be obtained for the relatively simple CARRGO model. To study the possible dynamics of the CARRGO model, we perform classical dynamical system analysis. Detailed mathematical analysis of this model can be found in several textbooks on dynamical systems [25,32]. In the interest of informing the reader, we briefly summarize the main points here. We begin by (1) scaling (non-dimensionalizing) the variables in the system and then (2) identify stationary points and classify their stability and finally (3) interpret the stationary points and system dynamics in terms of the initial numbers of cancer cells and CAR T-cells.

First, we scale the variables in the CARRGO model to obtain a model without physical units in order to study the intrinsic dynamics of the system. We scale time, the cancer cell and CAR T-cell populations as $\tau =t\rho ,\phantom{\rule{1em}{0ex}}y=\hspace{0.17em}({\kappa}_{1}/\rho )Y,\phantom{\rule{1em}{0ex}}x=X/K.$

These variables are substituted into the CARRGO model (equations (3.1) and (3.2)) to obtain the scaled dimensionless system

The steady-state solutions of this system (equations (3.3) and (3.4)) are obtained by setting the time derivatives equal to zero. The values of the dimensionless parameters *A* and *B* determine the dynamics of the system which may be represented as trajectories in a two-dimensional phase-space of cancer cells and CAR T-cells (*x*, *y*). Three stationary points corresponding to steady-state solutions are denoted by ${P}_{i}=(x,y)$: ${P}_{1}=(0,0)$ where both the cancer cells and CAR T-cells are eliminated, ${P}_{2}=(1,0)$ where the cancer cell population reaches the carrying capacity and CAR T-cells are eliminated, and a coexistence of both populations, ${P}_{3}=(A/B,(\hspace{0.17em}1-(A/B)))$. For a given initial condition, three possible dynamics can result from this model depending on the values of *A* and *B* (figure 1).

##### 3.2.1. Case 1: successful CAR T-cell treatment $(A=0,B>0)$

This situation occurs when the death rate of CAR T-cells is negligible $(\theta \approx 0)$ relative to the proliferation rate of cancer cells (*ρ*), and the CAR T-cells are stimulated to proliferate when encountering a cancer cell $({\kappa}_{2}>0)$. In this case, the equilibrium points are when the cancer cells are eliminated with some remaining CAR T-cells (0,*y*) and when the cancer cells reach carrying capacity with no remaining CAR T-cells (*K*,0). Cancer cell elimination (0,*y*) is stable only if the product of the rate of CAR T-cell killing and the number of CAR T-cells is greater than the proliferation rate of the cancer cells $(Y>\rho /{\kappa}_{1})$. If the initial CAR T-cell population satisfies the inequality $Y<(\rho /{\kappa}_{1})-(\rho X/K{\kappa}_{1})$, the CARRGO model predicts a *transient progression* of cancer cells before eventual response (grey region). The point (*K*,0) is an unstable repulsor state (figure 1*a*).

##### 3.2.2. Case 2: CAR T-cell treatment failure $(A=0,B<0)$

This situation occurs when the proliferation rate of CAR T-cells is less than the exhaustion rate of CAR T-cells due to interaction with cancer cells. In this case, the fixed points (0,*y*) and (*K*,0) are the same as in case 1. However, the point (*K*,0) is now a stable attractor state, corresponding to the extinction of CAR T-cells and the cancer cells eventually growing to carrying capacity. Again, cancer cell elimination (0,*y*) is stable only if the CAR T-cell population is larger than the ratio of cancer cell proliferation and CAR T-cell killing rate $(Y>\rho /{\kappa}_{1})$. If the initial CAR T-cell population satisfies the inequality $Y>(\rho /{\kappa}_{1})-(\rho X/K{\kappa}_{1})$, the CARRGO model predicts a *transient regression* of cancer cells before eventual rapid progression (grey region). This case is a failure of CAR T-cell treatment (figure 1*b*).

##### 3.2.3. Case 3: pseudo-failure or pseudo-response $(A>0,B>0)$

In this situation, the third stationary point ${P}_{3}$ corresponding to cancer cell and CAR T-cell coexistence lies in the first quadrant (positive numbers of cancer cells and CAR T-cells) only if $A\le B$. The point ${P}_{3}=A(B-A)/B$ is then a stable sink (figure 1*c*). This case results in an oscillating behaviour of an increase in cancer cells corresponding to tumour progression followed by a decrease in tumour cells corresponding to treatment response. The transient and oscillatory nature of these dynamics may be interpreted as a ‘pseudo’-failure and ‘pseudo’-response to the therapy. We note that cancer progression and treatment occur on finite and sometimes small time scales and therefore oscillatory dynamics may not be observed *in vivo* due to insufficient time to observe these changes.

#### 3.3. Cell lines

Low-passage primary brain tumour (PBT) lines were derived from GBM patients that had undergone tumour resections at City of Hope as previously described [23,33]. Fibrosarcoma line HT1080 was obtained from the American Tissue Culture Collection (ATCC) and maintained according to recommendations. The cell line PBT030 endogenously expresses high level of IL13Rα2. HT1080 and PBT138 do not express IL13Rα2 and were lentivirally engineered to express varied levels based on different promoter strengths to investigate the relationship between killing kinetics and antigen expression level: high (greater than 70%^{+}) driven by the EF1α promoter, medium (between 40%^{+} and 70%^{+}) driven by the PGK promoter, low (less than 20%^{+}) driven by the attenuated PGK100 promoter [34,35]. These cell lines are denoted with H, M, L, respectively (e.g. HT1080-H). These tumour cell lines were selected because they differ in aggressiveness (proliferation rates) and antigen expression levels (endogenous or engineered).

CAR T-cells were derived from healthy donor CD62 L+CD45RO+ central memory T-cell population and lentivirally transduced with second-generation IL13Rα2-targeting CARs: IL13BBζ or IL1328ζ [23,33,36,37]. Transduced product was enriched for CAR and expanded in X-Vivo media with 10% fetal bovine serum until 17 days in culture and cryopreserved. Non-transduced T-cells expanded under the same condition were used as mock control.

#### 3.4. Experimental design

Real-time monitoring of cancer cell growth was performed by using xCELLigence cell analyser system [38]. This system uses electrical impedance to non-invasively quantify adherent cell density with a dimensionless number referred to as cell-index (CI). The cell-index read-out is strongly positively correlated with the number of cells in the well (*r*^{2} > 0.9) and can be used as a measure of cell number [31]. We therefore report CARRGO parameter values in units CI which can be translated into units per cell based on the linear relation. Real-time cytotoxicity assay was performed using xCELLigence system in disposable 96 well E-Plates. Prior to seeding, tumour cells were enzymatically single-celled and seeded at 25 × 10^{3}, 12.5 × 10^{3} or 2 × 10^{3} cells per well depending on the cell line. Cells were either left untreated (triplicates per cell line) or treated with CAR T-cells at effector to target ratios (E : T) of 1 : 5, 1 : 10 and 1 : 20. CAR T-cells were added to the wells about 24 h after cancer cell seeding. Growth curves were recorded over 4 days with temporal resolution of 15 min (figure 2). Each cell line was treated with three IL13Rα2-targeted CAR T-cells: BBζ, 28ζ and mock. At the end of the experiment, flow cytometry was performed to measure the residual CAR T-cells, cancer cells and IL13Rα2 expression level. The details of cancer cell seeding and effector to target ratios used for the experiments are given in electronic supplementary material, table S1. The cancer cell dynamics of all the wells of 96 well E-plate for all cell lines are given in electronic supplementary material, figure S3.

#### 3.5. CARRGO model fitting to experimental data

The first 24 h of the time series describes the process of cell attachment to the bottom of the plate (figure 2). The spatial process of cell adhesion and spreading in the well can be modelled as a reaction–diffusion process, described in electronic supplementary material, figure S4. Since we are interested in cell growth kinetics, we omitted the data from the first 24 h during the attachment process. The final time point is determined when the cells reach confluency which varies by cell line. Observed time of confluency for the PBT cell line was around 120 h from the time of seeding while HT1080 reached confluency within 80 h. The data points from CAR T-cell administration (24 h after seeding) up to 80% of maximum CI value (confluency) were used for model fitting to estimate the parameters. At greater than 80% of the maximum CI, the linear relationship between CI and cell number no longer holds [31].

Cancer cell net growth rate *ρ* and carrying capacity *K* (equation (3.1)) were computed by fitting logistic growth to untreated cancer cell time-series data. The CAR T-cell killing rate *κ*_{1}, growth rate *κ*_{2} and death rate *θ* were computed by fitting the solution of the CARRGO model (equations (3.1) and (3.2)) to treated cancer cell time-series data by minimizing the root mean square error. Linear regression was used to determine the quality of model fitting (*R*^{2}). All optimization computations were performed in Matlab with *fmincon*.

#### 3.6. CARRGO model fitting to human data

A patient with recurrent glioma received CAR T-cells engineered for IL13Rα2 and showed complete tumour regression, which was published as a brief report by Brown *et al.* in 2016 [23]. We retrospectively collected the magnetic resonance imaging (MRI) data of this patient and calculated tumour volumes to be used to fit the CARRGO model. Three lesions were selected using the lesion labelling reported in Brown *et al.* [23]: lesions T6, T7 which responded to IL13Rα2 targeted CAR T-cells and lesion T9 which was a lesion that appeared later on which did not respond to the therapy. Tumour volumes for each lesion were estimated by manual segmentation of contrast-enhancing lesions from T1-weighted post-contrast MRIs. The number of cancer cells (CC) was estimated by calculating CC = [tumour volume (µm^{3})]/[GBM cell size (µm^{3})] with the average cell diameter assumed to be 20 µm [39]. The volume of a spherical cell is then given by $V=(4\pi /3){(10\hspace{0.17em}\mu \mathrm{m})}^{3}$. These relationships were used to estimate the total number of cancer cells in a tumour volume. The tumour growth rate (*ρ*) for lesions T6 and T7 was computed from two subsequent imaging time points following the first appearance of the lesion on MRI. Lesion T9 underwent surgical resection prior to CAR T-cell administration. The tumour growth rate for lesion T9 was computed from the pre-surgical MRI data. CAR T-cells were administered first directly into the tumour tissue and subsequently into the cerebrospinal fluid via the intraventricular injection. Because the CAR T-cells migrated to several tumour foci in the patient, we assumed a small fraction (5–10%) of the infused dose reached each individual lesion at each infusion. The CARRGO model was fitted to the time-series MRI-derived tumour volume data by minimizing the root mean square error for MRIs before and during CAR T-cell treatment to compute the rates of CAR T-cell killing *κ*_{1}, exhaustion *κ*_{2} and death *θ*. For lesion T9, the CARRGO model was fitted to the MRI data following CAR T-cell administration which occurred after the partial surgical resection.

### 4. Results

#### 4.1. Model/data fitting to *in vitro* data

A high goodness of fit of the CARRGO model to the xCELLigence data was observed across all cell lines (${R}^{2}=0.93\pm 0.1$, figure 3; electronic supplementary material, figure S3). To investigate the sensitivity of our model fitting to sampling frequency, we down-sampled the data by taking time intervals of 2, 5 and 10 h. No significant variation was observed in the model parameters *κ*_{1}, *κ*_{2} and *θ* to the down-sampled data (repeated measure ANOVA *p* > 0.1) (electronic supplementary material, figure S5 and S6). We consistently observed very small values of the CAR T-cell death rate $(\theta <{10}^{-3})$. Uniqueness of the parameters was tested by choosing 100 different combinations of values of the parameters across several orders of magnitude for the model fitting optimization procedure. We found that if the optimization converged, it converged to unique values of the parameters, which is a direct consequence of the identifiability analysis of the model and minimum number of points required to resolve the model (see electronic supplementary material, S1, figure S7 and movie S1).

#### 4.2. Validation of xCELLigence dynamics with flow cytometry

Because the xCELLigence system is an indirect measure of cell number, we validated the previously reported linear relationship [38] between the cell index read-out from the machine and the number of cells measured with flow cytometry (*R*^{2} > 0.9; electronic supplementary material, S2 and figure S8). Because cell index measures the change in electrical impedance caused by both cancer cells and CAR T-cells adhering to the well, CAR T-cell dynamics are not directly measured by the system, so we compared the cancer cell (CC) to T-cell ratio (TC) from flow cytometry to that predicted from the CARRGO model. The model predicted ratio CC/TC at end time point shows a similar trend to that measured with flow cytometry, indicating the CARRGO model-predicted CAR T-cell dynamics derived from the xCELLigence data are consistent with flow cytometry measurements. This trend was observed in PBT030 and PBT138 for BBζ and 28ζ CAR T-cells and for all doses (electronic supplementary material, S2 and figure S9).

#### 4.3. CAR T-cell dose-dependent dynamics

We examined the effect of varying the effector to target ratio, i.e. CAR T-cell dose, for all cell lines. The CAR T-cell death rate parameter was found to be very small ($\theta <{10}^{-3}$ day^{−1}) for all cancer cell lines and all CAR T-cells and doses. The killing rate parameter *κ*_{1} shows a negative correlation while *κ*_{2} shows positive correlation with respect to the CAR T-cell dose. The parameter *κ*_{2} was negative only for tumour line HT1080-H which indicates the exhaustion rate being much stronger than the proliferation rate of the CAR T-cells. A positive correlation of *κ*_{1} with CAR T-cell dose indicates that higher dose of CAR T-cell results in a lower killing rate, since each individual CAR T-cell has fewer number of cancer cells to encounter. The range of values for κ_{1} and κ_{2} varied with cancer cell line as the growth rate of each cell line is different from each other; however, the overall trends were preserved across the cell lines. Plots of *θ, κ*_{1} and *κ*_{2} for PBT030, PBT138 (seeding 12.5 × 10^{3}) and HT1080-H treated with BBζ CAR T-cells are shown in figure 4. Other cell lines and parameters for 28ζ CAR T-cells are given in electronic supplementary material, S2 and figure S10.

#### 4.4. Relating *κ*_{1}, *κ*_{2} with tumour growth rate and antigen expression

Tumour growth rate *ρ* varies significantly (*p* < 0.01) among different cell lines and with antigen expression level (see electronic supplementary material S2 and figure S11a). To investigate the relationship between tumour growth rate and CAR T-cell killing *κ*_{1} and exhaustion *κ*_{2}, we evaluated cell lines with antigen levels greater than 80% and treated with BBζ CAR T-cells at an effector to target ratio of 1 : 5. No significant correlation was found between the cancer cell proliferation rate *ρ* and killing rate (*κ*_{1}) (electronic supplementary material, figure S11b). However, the exhaustion rate *κ*_{2} is significantly correlated with tumour growth rate (electronic supplementary material, figure S11c) with Pearson correlation coefficient $r=-0.9$, *p* < 0.001. Similar results were observed for the cells treated with 28ζ IL13Rα2-CARs. Figure 5 shows the density of IL13Rα2 level on cancer cell surface and its relation to CAR T-cell killing for cell line HT1080-H and PBT138-H. We observed that *κ*_{1} shows a decreasing trend from medium to high antigen level (figure 5*c*) suggesting that high levels of antigen expression may not result in faster rates of CAR T-cell killing. The rate constant *κ*_{2} increases from low to medium antigen expression and plateaus with high levels (figure 5*d*). This suggests limited activation of CAR T-cell at lower antigen expression and exhaustion rate from medium to high antigen may not change significantly and may be the result of over-activation of the CAR T-cells.

#### 4.5. CARRGO model applied to *in vivo* human data

To translate the *in vitro* dynamics of the model to clinical data [23], we fit the CARRGO model to MRI-derived tumour volumes during CAR T-cell treatment (figure 6). The CARRGO model is able to fit the tumour growth dynamics quite accurately for lesions T6 and T7 with the same set of parameters ${\kappa}_{1}=6\times {10}^{-9}$(day^{−1} cell^{−1}), ${\kappa}_{2}=0.3\times {10}^{-10}$ (day^{−1} cell^{−1}), $\theta =0.1\times {10}^{-5}$(day^{−1}) and lesion T9 with ${\kappa}_{1}=9\times {10}^{-8}$(day^{−1} cell^{−1}), ${\kappa}_{2}=-2\times {10}^{-13}$ (day^{−1} cell^{−1}), $\theta =5\times {10}^{-5}$ (day^{−1}). In the case of lesion T9, although the CARRGO model is consistent with the overall tumour dynamics, it does not fit the later time points following CAR T-cell treatment well. This is because lesion T9 received radiation treatment between day 200 and 300, which is not included in the CARRGO model. We note the negative correlation between the tumour growth rate ($\rho =0.06,0.07$ and 0.2 day^{−1} for T6, T7 and T9, respectively) with the CAR T-cell exhaustion rate *κ*_{2} in the patient data, which is consistent with that observed in the experimental data (electronic supplementary material, figure S11c). We remark that the parameters *κ*_{1} and *κ*_{2} are of the order of $O({10}^{-13})$, which appear to be very small; however, these parameters are scaled by the carrying capacity *K* in units of cells, which is of order $O({10}^{9})$. Therefore, these parameter values are comparable with the *in vitro* data when scaled relative to the carrying capacity (figure 4).

### 5. Discussion

The CARRGO model is a simple representation of cancer cell–T-cell interactions. We developed the CARRGO model with the aim of understanding CAR T-cell efficacy in terms of rates of killing, proliferation, exhaustion, and persistence with a real-time cell analyser in a simple, controlled *in vitro* system. The CARRGO model fitted remarkably well to the highly temporally resolved experimental data and as well as to data derived from a patient treated with IL13Rα2 BBζ CAR T-cells. Although the predator–prey mathematical model formalism has been widely used in a number of biological settings, the novelty of this model is in the application to a novel form of cancer therapy with a high temporal resolution cell monitoring experimental design which provides nearly continuous data on killing kinetics.

With the CARRGO model we show that the rate of cancer cell killing by CAR T-cells is inversely related to the CAR T-cell dose. With a fixed number of cancer cells as an initial condition, as the number of CAR T-cells increases (dose), any individual T-cell will encounter fewer number of cancer cells to kill, indicating that increasing dose does not result in a maximal rate of killing on a per T-cell basis. For example, the PBT138 cell line shows complete killing within 80 and 100 h for effector to target ratios 1 : 5 and 1 : 10, respectively. This result suggests that a lower dose of CAR T-cells may change the time to complete cancer cell killing but shows the same overall cancer cell killing effectiveness. Moreover, we observed that *κ*_{2} positively correlated with CAR T-cell dose. Because the parameter *κ*_{2} is a net measure of CAR T-cell proliferation and exhaustion or lack of activation and the T-cell proliferation is not dose-dependent [18], the trend observed in *κ*_{2} with dose is dominated by the exhaustion rate: the higher the dose, the lower the exhaustion rate, resulting in an increased value of *κ*_{2}. The death rate of CAR T-cells was very small when compared with the cancer cell proliferation rate for all conditions. This is likely due to the short time scale of the experiment and because the T-cells were stimulated to proliferate by the presence of cancer cells.

We observed that the cancer cell growth showed no relation with CAR T-cell killing rate and an inverse relationship with *κ*_{2}. This may explain variations in patient-specific responses even for the same CAR T-cell dose. For a fixed CAR T-dose, *κ*_{2} is the principal determinant of treatment failure or success as shown in phase-plane analysis (figure 1), which is also observed in patient data (figure 6). This result, driven by the CARRGO model analysis suggests that the balance between proliferation and exhaustion of CAR T-cells may contribute more than the rate of CAR T-cell mediated cancer cell killing in determining treatment success or failure. Moreover, the CARRGO model predicts transient progression of cancer cells even in the case of successful CAR T-cell therapy. This prediction may be consistent with the clinical phenomenon of pseudo-progression, in which the cancer is seen to progress during therapy before eventually responding [7,8]. Identifying characteristics of the patient and the CAR T-cells which may result in pseudo-progression could have a profound effect on interpretation of these dynamics observed in the clinic.

Interestingly, we found *κ*_{1} decreases and *κ*_{2} plateaued from medium antigen level to higher level of antigen expression. One of the possible explanations of this behaviour could be the antigen density is more heterogeneous in the higher antigen level cell population when compared with medium and low antigen levels (figure 5*a*). More heterogeneity in the density of antigen expression intensity in the cancer cells within the initial population may cause clustering of CAR T-cells resulting in their exhaustion [20,40]. Another confounding factor can be the dependence of the detected antigen signal intensity on both the number of antigen-positive tumour cells and their individual antigen expression intensity. Elucidating these factors individually can better tune the model parameters and the prediction of the tumour response dynamics. However, more studies are required to examine effect of cell and population level antigen density on CAR T-cell killing kinetics.

There are some important limitations to consider with this model and experimental system. Perhaps the most obvious is that the *in vitro* system is not a model for the human immune system or tumour microenvironment. It does not include cytokines, stromal cells or additional immune cells such as myeloid cells which contribute to CAR T-cell activity *in vivo*. Another limitation is the assumption that the populations are well mixed. In practice, this assumption may depend on the route of CAR T-cell administration, with intracavitary and intraventricular injections potentially resulting in spatially heterogeneous densities of CAR T-cells or cancer cells, although methods to assess the distribution of CAR T-cells *in vivo* remain an open challenge [41,42]. To address this limitation, the well-mixed assumption may be relaxed and CAR T-cell killing dynamics interrogated with spatial or agent-based models [43]. Another limitation is with regard to the experimental system: the change in electrical impedance measured by the cell index does not differentiate cell detachment from cell killing. This is only a minor consideration as the cell lines used are very adherent to the plate and were not observed to detach. Finally, the experimental system does not directly measure dynamics of CAR T-cells. However, our model is initialized with known numbers of cancer cells and CAR T-cells and the model-predicted cancer cell to CAR T-cell ratio at the experimental endpoint was validated with flow cytometry, giving confidence to our model predictions and parameter estimates. To address this limitation, the CARRGO model CAR T-cell dynamics can be validated by labelling the CAR T-cells and directly measuring their dynamics with live cell imaging-based methods [44]. Despite these limitations, the CARRGO model succeeded in revealing nonlinear dynamics, quantifying kinetics of killing, and generating hypotheses which may be tested in other *in vitro* systems, and other computational or *in vivo* models.

An interesting application of this model would be to adapt the system to evaluate differences between donor T-cells—either autologous or allogeneic. The intrinsic fitness of the T-cell used for CAR T-cell manufacturing is known to be a critical differentiator between responding and non-responding patients [45], and this model may be able to predict T-cell products with high proliferative potential versus products more prone to exhaustion. While our studies used an allogeneic co-culture platform for the predator–prey modelling of effector activity, we expect these findings to be translatable to autologous CAR T-cell therapies, as CAR T-cells recognize targets in an MHC-independent manner. Future studies could specifically evaluate allogeneic and autologous donor-dependent differences to assess the fit of the mathematical model, and guide potential adjustments to the model. The ability to differentiate key quality attributes of the therapeutic product could potentially provide powerful predictions of CAR T-cell product efficacy and facilitate patient treatment management.

In summary, CAR T-cells have shown promise in hematologic malignancies and are being actively investigated in solid tumours. We aimed to use mathematical modelling to investigate factors which contribute to the kinetics of CAR T-cell mediated cancer cell killing in a simple isolated *in vitro* system. We were able to fit the CARRGO model to *in vitro* and *in vivo* human data with remarkable accuracy. We demonstrated that we can consistently and reproducibly estimate rate constants in the CARRGO model and investigate their dependence on CAR T-cell dose and antigen expression levels. The CARRGO model may be combined with other mathematical models which estimate cancer cell growth and proliferation rates non-invasively with MRI data [9,11,46] to produce a fine-tuned and benchmarked suite of mathematical models, which may aide in optimization of dosing and scheduling of CAR T-cells for greater individualized and personalized therapy.

### Data accessibility

All raw experimental data are included in the electronic supplementary material.

### Authors' contributions

P.S., X.Y., C.B., R.R. conceived of the study and performed the experiments, mathematical modelling and analysis. D.A., D.M., V.A., D.F., H.C., V.M., D.W., M.B., M.G., S.B. provided development, analysis and interpretation of mathematical model and data. All authors contributed to writing and approving the final manuscript.

### Competing interest

Patents associated with CAR design, T-cell manufacturing and delivery have been licensed by Mustang Bio., Inc, for which C.E.B. receives licensing and consulting payments. The other authors declare no potential conflicts of interest.

### Funding

The research reported in this publication was supported by the California Institute for Regenerative Medicine (CLIN2-10248) and the National Cancer Institute of the National Institutes of Health under grant numbers R01CA236500, P30CA033572. The content is solely the responsibility of the authors and does not necessarily represent the official views of the California Institute for Regenerative Medicine or National Institutes of Health.

### References

- 1.
Brentjens RJ 2003 Eradication of systemic B-cell tumors by genetically targeted human T lymphocytes co-stimulated by CD80 and interleukin-15.**Nat. Med.**, 279-286. (doi:10.1038/nm827) Crossref, PubMed, Web of Science, Google Scholar**9** - 2.
Milone MC 2009 Chimeric receptors containing CD137 signal transduction domains mediate enhanced survival of T cells and increased antileukemic efficacy in vivo.**Mol. Ther.**, 1453-1464. (doi:10.1038/mt.2009.83) Crossref, PubMed, Web of Science, Google Scholar**17** - 3.
Filley AC, Henriquez M, Dey M . 2018 CART immunotherapy: development, success, and translation to malignant gliomas and other solid tumors.**Front. Oncol.**, 453. (doi:10.3389/fonc.2018.00453) Crossref, PubMed, Web of Science, Google Scholar**8** - 4.
Li J, Li W, Huang K, Zhang Y, Kupfer G, Zhao Q . 2018 Chimeric antigen receptor T cell (CAR-T) immunotherapy for solid tumors: lessons learned and strategies for moving forward.**J. Hematol. Oncol.**, 22. (doi:10.1186/s13045-018-0568-6) Crossref, PubMed, Web of Science, Google Scholar**11** - 5.
Newick K, Moon E, Albelda SM . 2016 Chimeric antigen receptor T-cell therapy for solid tumors.**Mol. Ther. Oncolytics**, 16006. (doi:10.1038/MTO.2016.6) Crossref, PubMed, Web of Science, Google Scholar**3** - 6.
Bagley SJ, Desai AS, Linette GP, June CH, O'Rourke DM . 2018 CAR T-cell therapy for glioblastoma: recent clinical advances and future challenges.**Neuro. Oncol.**, 1429-1438. (doi:10.1093/neuonc/noy032) Crossref, PubMed, Web of Science, Google Scholar**20** - 7.
Ranjan S 2018 Clinical decision making in the era of immunotherapy for high grade-glioma: report of four cases.**BMC Cancer**, 239. (doi:10.1186/s12885-018-4131-1) Crossref, PubMed, Web of Science, Google Scholar**18** - 8.
Thust SC, van den Bent MJ, Smits M . 2018 Pseudoprogression of brain tumors.**J. Magn. Reson. Imaging**, 571. (doi:10.1002/jmri.26171) Crossref, Web of Science, Google Scholar**48** - 9.
Rockne R 2010 Predicting efficacy of radiotherapy in individual glioblastoma patients in vivo: a mathematical modeling approach.**Phys. Med. Biol.**, 3271-3285. (doi:10.1088/0031-9155/55/12/001) Crossref, PubMed, Web of Science, Google Scholar**55** - 10.
Gaw N 2019 Integration of machine learning and mechanistic models accurately predicts variation in cell density of glioblastoma using multiparametric MRI.**Sci. Rep.**, 10063. (doi:10.1038/s41598-019-46296-4) Crossref, PubMed, Web of Science, Google Scholar**9** - 11.
Rayfield CA 2018 Distinct phenotypic clusters of glioblastoma growth and response kinetics predict survival.**JCO Clin. Cancer Inform.**, 1-14. (doi:10.1200/CCI.17.00080) Crossref, Web of Science, Google Scholar**1** - 12.
Eftimie R, Bramson JL, Earn DJD . 2011 Interactions between the immune system and cancer: a brief review of non-spatial mathematical models.**Bull. Math. Biol.**, 2-32. (doi:10.1007/s11538-010-9526-3) Crossref, PubMed, Web of Science, Google Scholar**73** - 13.
de Pillis LG, Radunskaya AE, Wiseman CL . 2005 A validated mathematical model of cell-mediated immune response to tumor growth.**Cancer Res.**, 7950-7958. (doi:10.1158/0008-5472.CAN-05-0564) Crossref, PubMed, Web of Science, Google Scholar**65** - 14.
Frascoli F, Kim PS, Hughes BD, Landman KA . 2014 A dynamical model of tumour immunotherapy.**Math. Biosci.**, 50-62. (doi:10.1016/j.mbs.2014.04.003) Crossref, PubMed, Web of Science, Google Scholar**253** - 15.
Hardiansyah D, Ng CM . 2019 Quantitative systems pharmacology model of chimeric antigen receptor T-cell therapy.**Clin. Transl. Sci.**, cts.12636. (doi:10.1111/cts.12636) Crossref, Web of Science, Google Scholar**12** - 16.
Stein AM 2017 CTL019 model-based cellular kinetic analysis of chimeric antigen receptor (CAR) T cells to characterize the impact of tocilizumab on expansion and to identify correlates of cytokine release syndrome severity.**Blood**, 2561. (doi:10.1182/blood-2017-04-779405) Web of Science, Google Scholar**130** - 17.
Hopkins B, Tucker M, Pan Y, Fang N, Huang ZJ . 2018 A model-based investigation of cytokine storm for T-cell therapy.**IFAC-PapersOnLine**, 76-79. (doi:10.1016/J.IFACOL.2018.09.039) Crossref, Google Scholar**51** - 18.
Stein AM 2019 Tisagenlecleucel model-based cellular kinetic analysis of chimeric antigen receptor–T cells.**CPT Pharmacometrics Syst. Pharmacol.**, 285-295. (doi:10.1002/psp4.12388) Crossref, PubMed, Web of Science, Google Scholar**8** - 19.
Rohrs JA, Wang P, Finley SD . 2019 Understanding the dynamics of T-cell activation in health and disease through the lens of computational modeling.**JCO Clin. Cancer Inform.**, 1-8. (doi:10.1200/CCI.18.00057) Crossref, PubMed, Web of Science, Google Scholar**3** - 20.
Harris DT . 2018 Comparison of T cell activities mediated by human TCRs and CARs that use the same recognition domains.**J. Immunol.**, 1088-1100. (doi:10.4049/jimmunol.1700236.Comparison) Crossref, PubMed, Web of Science, Google Scholar**200** - 21.
Kimmel GJ, Locke FL, Altrock PM . 2019 Evolutionary dynamics of CAR T cell therapy. bioRxiv 717074. (doi:10.1101/717074) Google Scholar - 22.
Brown CE 2018 Optimization of IL13Rα2-targeted chimeric antigen receptor T cells for improved anti-tumor efficacy against glioblastoma.**Mol. Ther.**, 31-44. (doi:10.1016/j.ymthe.2017.10.002) Crossref, PubMed, Web of Science, Google Scholar**26** - 23.
Brown CE 2016 Regression of glioblastoma after chimeric antigen receptor T-cell therapy.**N. Engl. J. Med.**, 2561-2569. (doi:10.1056/NEJMoa1610497) Crossref, PubMed, Web of Science, Google Scholar**375** - 24.
Murray JD . 2002**Mathematical biology I: an introduction**, 3rd edn. New York, NY: Springer. Crossref, Google Scholar - 25.
Edelstein-Keshet L . 2005**Mathematical models in biology**. Philadelphia, PA: Society for Industrial and Applied Mathematics. Crossref, Google Scholar - 26.
Villaverde AF . 2019 Observability and structural identifiability of nonlinear biological systems.**Complexity**, 8497093. (doi:10.1155/2019/8497093) Crossref, Web of Science, Google Scholar**2019** - 27.
Ligon TS, Fröhlich F, Chiş OT, Banga JR, Balsa-Canto E, Hasenauer J . 2018 GenSSI 2.0: multi-experiment structural identifiability analysis of SBML models.**Bioinformatics**, 1421-1423. (doi:10.1093/bioinformatics/btx735) Crossref, PubMed, Web of Science, Google Scholar**34** - 28.
Sontag ED . 2017 Dynamic compensation, parameter identifiability, and equivariances.**PLoS Comput. Biol.**, e1005447. (doi:10.1371/journal.pcbi.1005447) Crossref, PubMed, Web of Science, Google Scholar**13** - 29.
Villaverde AF, Barreiro A, Papachristodoulou A . 2016 Structural identifiability of dynamic systems biology models.**PLoS Comput. Biol.**, e1005153. (doi:10.1371/journal.pcbi.1005153) Crossref, PubMed, Web of Science, Google Scholar**12** - 30.
Mayer A, Zhang Y, Perelson AS, Wingreen NS . 2019 Regulation of T cell expansion by antigen presentation dynamics.**Proc. Natl Acad. Sci. USA**, 5914-5919. (doi:10.1073/pnas.1812800116) Crossref, PubMed, Web of Science, Google Scholar**116** - 31.
Witzel F, Fritsche-Guenther R, Lehmann N, Sieber A, Blüthgen N . 2015 Analysis of impedance-based cellular growth assays.**Bioinformatics**, 2705-2712. (doi:10.1093/bioinformatics/btv216) Crossref, PubMed, Web of Science, Google Scholar**31** - 32.
Murray JD . 1993 Models for interacting populations. In**Mathematical biology I. An introduction**. Interdisciplinary Applied Mathematics, vol. 17, pp. 79–118. New York, NY: Springer. (doi:10.1007/978-0-387-22437-4_3) Google Scholar - 33.
Brown CE, Starr R, Aguilar B, Shami AF, Martinez C, D'Apuzzo M, Barish ME, Forman SJ, Jensen MC . 2012 Stem-like tumor-initiating cells isolated from IL13Rα2 expressing gliomas are targeted and killed by IL13-zetakine-redirected T cells.**Clin. Cancer Res.**, 2199-2209. (doi:10.1158/1078-0432.CCR-11-1669) Crossref, PubMed, Web of Science, Google Scholar**18** - 34.
Priceman SJ 2018 Co-stimulatory signaling determines tumor antigen sensitivity and persistence of CAR T cells targeting PSCA+ metastatic prostate cancer.**Oncoimmunology**, e1380764. (doi:10.1080/2162402X.2017.1380764) Crossref, PubMed, Web of Science, Google Scholar**7** - 35.
Frigault MJ 2015 Identification of chimeric antigen receptors that mediate constitutive or inducible proliferation of T cells.**Cancer Immunol. Res.**, 356-367. (doi:10.1158/2326-6066.CIR-14-0186) Crossref, PubMed, Web of Science, Google Scholar**3** - 36.
Debinski W, Gibo DM, Hulet SW, Connor JR, Gillespie GY . 1999 Receptor for interleukin 13 is a marker and therapeutic target for human high-grade gliomas.**Clin. Cancer Res.**, 985-990. PubMed, Web of Science, Google Scholar**5** - 37.
Kahlon KS, Brown C, Cooper LJN, Raubitschek A, Forman SJ, Jensen MC . 2004 Specific recognition and killing of glioblastoma multiforme by interleukin 13-zetakine redirected cytolytic T cells.**Cancer Res.**, 9160-9166. (doi:10.1158/0008-5472.CAN-04-0454) Crossref, PubMed, Web of Science, Google Scholar**64** - 38.
Roshan MM, Young A, Reinheimer K, Rayat J, Dai L-J, Warnock GL . 2015 Dynamic assessment of cell viability, proliferation and migration using real time cell analyzer system (RTCA).**Cytotechnology**, 379-386. (doi:10.1007/s10616-014-9692-5) Crossref, PubMed, Web of Science, Google Scholar**67** - 39.
Oraiopoulou M-E, Tzamali E, Tzedakis G, Vakis A, Papamatheakis J, Sakkalis V . 2017*In vitro*/*in silico*study on the role of doubling time heterogeneity among primary glioblastoma cell lines.**Biomed Res. Int.**, 8569328. (doi:10.1155/2017/8569328) Crossref, PubMed, Web of Science, Google Scholar**2017** - 40.
Utzschneider DT, Alfei F, Roelli P, Barras D, Chennupati V, Darbre S, Delorenzi M, Pinschewer DD, Zehn D . 2016 High antigen levels induce an exhausted phenotype in a chronic infection without impairing T cell expansion and survival.**J. Exp. Med.**, 1819-1834. (doi:10.1084/jem.20150598) Crossref, PubMed, Web of Science, Google Scholar**213** - 41.
Weist MR 2018 Positron emission tomography of adoptively transferred chimeric antigen receptor T cells with zirconium-89 oxine.**J. Nucl. Med.**, jnumed.117.206714. (doi:10.2967/jnumed.117.206714) Google Scholar**2018** - 42.
Keu KV 2017 Reporter gene imaging of targeted T cell immunotherapy in recurrent glioma.**Sci. Transl. Med.**, eaag2196. (doi:10.1126/scitranslmed.aag2196) Crossref, PubMed, Web of Science, Google Scholar**9** - 43.
Ghaffarizadeh A, Heiland R, Friedman SH, Mumenthaler SM, Macklin P . 2018 PhysiCell: an open source physics-based cell simulator for 3-D multicellular systems.**PLoS Comput. Biol.**, e1005991. (doi:10.1371/journal.pcbi.1005991) Crossref, PubMed, Web of Science, Google Scholar**14** - 44.
Gelles JD, Chipuk JE . 2016 Robust high-throughput kinetic analysis of apoptosis with real-time high-content live-cell imaging.**Cell Death Dis.**, e2493. (doi:10.1038/cddis.2016.332) Crossref, PubMed, Web of Science, Google Scholar**7** - 45.
Fraietta JA 2018 Determinants of response and resistance to CD19 chimeric antigen receptor (CAR) T cell therapy of chronic lymphocytic leukemia.**Nat. Med.**, 563-571. (doi:10.1038/s41591-018-0010-1) Crossref, PubMed, Web of Science, Google Scholar**24** - 46.
Wang CH 2009 Prognostic significance of growth kinetics in newly diagnosed glioblastomas revealed by combining serial imaging with a novel biomathematical model.**Cancer Res.**, 9133-9140. (doi:10.1158/0008-5472.CAN-08-3863) Crossref, PubMed, Web of Science, Google Scholar**69**