## Highlights

- •Model of glioma and CAR T-cells interaction that includes multi-cell conjugates.
- •Enhanced model fitness and accuracy in describing the experimental data.
- •Derivation of CAR T-cell expansion rate that yields treatment success or failure.
- •Study of interaction rates to model different levels of antigen receptor density.

## Abstract

Chimeric antigen receptor (CAR) T-cell based immunotherapy has shown its potential in treating blood cancers, and its application to solid tumors is currently being extensively investigated. For glioma brain tumors, various CAR T-cell targets include IL13R$\alpha $2, EGFRvIII, HER2, EphA2, GD2, B7-H3, and chlorotoxin. In this work, we are interested in developing a mathematical model of IL13R$\alpha $2 targeting CAR T-cells for treating glioma. We focus on extending the work of Kuznetsov et al. (1994) by considering binding of multiple CAR T-cells to a single glioma cell, and the dynamics of these multi-cellular conjugates. Our model more accurately describes experimentally observed CAR T-cell killing assay data than the models which do not consider multi-cellular conjugates. Moreover, we derive conditions in the CAR T-cell expansion rate that determines treatment success or failure. Finally, we show that our model captures distinct CAR T-cell killing dynamics from low to high antigen receptor densities in patient-derived brain tumor cells.

## Graphical abstract

## Keywords

## 1. Introduction

Adoptive cell-based immunotherapy has shown to be successful in treating patients with cancer. In particular, Chimeric Antigen Receptor $T$ cell (CAR T-cell) therapy is one of the adoptive immunotherapies that has been successful in clinical and pre-clinical models, and has been FDA-approved since 2017 [

1

, 2

]. In this therapy, a patient or donor’s $T$ cells are collected and genetically engineered to express a receptor specific to an antigen found on cancer cells, thus improving the ability of T-cells to eradicate the target cancer cells. Finally, these CAR T-cells are cultured to large numbers, then introduced back to the patient [3

, 4

]. CAR T-cell therapy has shown its potential in blood cancer, and solid tumors [3

, 5

]. However, the success of CAR T-cell therapy for solid tumors has been challenging due to difficulties in (1) trafficking CAR T-cells into solid tumors, (2) hostile tumor microenvironment that suppresses $T$ cell activity, and (3) tumor antigen heterogeneity [5

, 6

]. Since CAR T-cells mostly exist in bloodstream and lymphatic system, which makes CAR T-cell therapy a great weapon for hematological tumors (e.g. blood tumor cells), it may be hard for CAR T-cells to penetrate tumor tissue through the vascular endothelium. Moreover, various types of cells infiltrate solid tumors to support tumor growth and restrict the efficacy of CAR T-cell therapy [[6]

].One of the first mathematical models describing the interaction between immune cells and cancer cells is from Kuznetsov et al. (1994) [

[7]

], which is a dynamical system model with two populations, tumor cells and cytotoxic $T$ cells, or $T$ lymphocytes. The model can describe the formation of a tumor dormant state and evasion of the immune system. A subsequent model developed by Kirschner and Panetta (1998) [[8]

] considered the cytokine interleukin-2 in addition to the dynamics between tumor cells and immune effector cells, and was able to model short-term tumor oscillations as well as long-term tumor relapses. In order to model persistent oscillations that were observed in immune systems, periodic treatment and time delay was added to the model in [[9]

], and stability analysis of the model was done in [[10]

]. Thereafter, the later built models were improved by adding new types of cells, such as natural killer cells, normal cells, and different kinds of cytokines [11

, 12

]. These models not only capture tumor immune escape, but also explain multiple equilibrium phases of coexisting immune cells and cancer cells. Other than dynamical system models, spatial models that describe the spatial temporal interaction are developed as well. [[13]

] develops a spatial temporal version of [[7]

] using partial differential equations, [[14]

] develops a hybrid cellular automata-partial differential equation model, and [[15]

] develops a hybrid off-lattice agent-based and partial differential equation model.The surge of clinical trials and the success of CAR T-cell therapy also drew a lot of interest in mathematical modeling of CAR T-cell therapy. This includes modeling CD19 CAR T-cell therapy targeting acute lymphoblastic leukemia in [

[16]

] as a dynamical system, which includes healthy B cell populations and circulating lymphocytes. Later study in [[17]

] further showed relationships between CAR T-cell doses and diseases burden by the observed clinical data. In order to study cytokine release syndrome, which is one of the primary side effects of CAR T-cell therapy, another dynamical system of nine cytokines responding to CAR T-cell therapy was proposed in [[18]

] as well. More recent work in [[19]

] attempted to understand the dynamics of CAR T-cell therapy by considering not only tumor cells and CAR T-cells but also normal $T$ cells. Meanwhile, the model introduced in [[20]

] includes long-term memory CAR T-cells, which are produced by memory pool formation of effector CAR T-cells. The corresponding stability analysis of this model was done in [[21]

].Here we develop a mathematical model of glioma cells and CAR T-cells inspired by the experimental data provided in [

[22]

]. The experiments study the interaction between glioma cell lines, derived from glioblastoma patients undergoing tumor resections at City of Hope [23

, 24

], and IL13R$\alpha $2 targeting CAR T-cells. As shown in Fig. 1, cells were co-cultured in vitro and images were taken under a light microscope over a 72 h period. The experimental images illustrate two aspects of the glioma cell killing of CAR T-cells. First, the elimination process of glioma by CAR T-cell is not immediate, and second, multiple CAR T-cells can interact with glioma cells. Such phenomena can be also observed in a close-up video in []. In subsequent experiments, the glioma cells and CAR T-cells are mixed at different ratios (CAR T-cell to glioma cell ratios of 1:5, 1:10, and 1:20), and glioma cells with different antigen receptor density levels (low, medium, and high) were tested. We remark that number of CAR T-cells binding to glioma cells will depend on the antigen receptor density levels since glioma cells with high antigen receptor density levels will have more capacity for CAR T-cells to bind. Throughout the course of experiment, real-time monitoring of glioma population size was performed by using xCELLigence cell analyzer system [[26]

]. This system quantifies the glioma cell population with a dimensionless number referred to as cell-index (CI) (1 CI $\approx 1{0}^{4}$ cells [[22]

]), where the population size is tracked every 15 min. See appendix Fig. C.1 for the full set of data.The rest of the paper is structured as follows. In Section 2, we describe the proposed model, where we assume that one glioma cell may interact with multiple CAR T-cells, and further come up with a system of ODE that describes not only the size of tumor cells and CAR T-cells, but also the size of multiple CAR T-cell binding conjugates. We denote the model with conjugates up to $n$ CAR T-cells binding to glioma cell as the $n$-binding model. The stability analysis of one-binding slow reaction model is presented in Section 3, and we provide parameter conditions that guarantee either CAR T-cell treatment success or failure. In Section 4, we compare the accuracy between the slow reaction and fast reaction one-binding model, and show that the slow reaction model more accurately describes the experimental data. We also simulate the multiple binding slow reaction models from one- to five-binding conjugates, and study the hypotheses in reaction rates that captures the experimental result regarding low to high antigen receptor density levels of glioma cells. Summary of our findings and future work is discussed in Section 5.

## 2. Mathematical model

This section summarizes the mathematical models that we study in this work. The motivation of building these models comes from the experiments that show multiple CAR T-cells interacting with glioma cells [] and non-instantaneous killing process that can take hours. We start by deriving the dynamical system model proposed in Kuznetsov et al. (1994) [

[7]

] that considers one CAR T-cell interacting with one glioma cell. The two versions of the model – slow and fast reaction models – are compared. Then we extend the models to consider the multi-cellular conjugate of multiple CAR T-cells and glioma cell.### 2.1 One CAR T-cell bound to one glioma cell model

Among the mathematical models considering interaction of glioma cells and immune cells, one of the most recognized models is the model developed in Kuznetsov et al. (1994) [

Here, parameter $\rho $ is the maximal growth rate of the glioma cells assuming logistic growth; parameter $K$ is the maximal carrying capacity of the biological environment for the glioma cells; $p$ is the rate that the CAR T-cells accumulate in the region due to the presence of the tumor, in other words, CAR T-cell expansion rate; $g$ is a concentration of glioma cells that halves the maximum rate $p$; $\theta $ represents the death rate of the CAR T-cells. Lysed cancer cells and CAR T-cells will be formed at the end of this interaction, but since their dynamics are simply decaying, we do not include those equations in our model as in [

where

The new parameter $\alpha $ represents the rate of net inflow of CAR T-cells from the conjugate ${I}_{1}$, $\beta $ represents the inflow rate for glioma cells, and $\gamma $ stands for the net decay rate of the conjugate ${I}_{1}$ due to either detachment or lysed reaction. We will refer to Eqs. (1) or (2) as the

[7]

], modeling the one to one binding of cancer and cytotoxic immune cells. The glioma cells are subject to be attacked by cytotoxic effector cells, e.g. cytotoxic $T$ lymphocytes and natural killer cells. Here, we consider CAR T-cells instead of the general effector cells. The killing process of CAR T-cell has mainly 4 stages: binding, recognition, lethal hit, and target cell rounding [5

, 27

]. More explicitly, tumor cells are first scanned by CAR T-cells, then antigens are recognized through proteins or glycans on the surface. After that, tumor cells are killed by CAR T-cells and cell conjugates lose adhesion and commit to cancer cells’ death. In our model, we describe the interaction between the glioma cancer cells $C\left(t\right)$ and CAR T-cells $T\left(t\right)$ by the kinetic scheme illustrated in Fig. 2 (top), where ${I}_{1}\left(t\right)$ is the conjugate of one CAR T-cell and one glioma cell. The kinetic parameter ${k}_{1}^{\left(1\right)}$ describes the rate that the glioma cell binds to the CAR T-cell, and ${k}_{-1}^{\left(1\right)}$ is the rate that the conjugate detaches without damaging the cells. ${k}_{3}^{\left(1\right)}$ is the rate that the CAR T-cell and glioma cell interaction kills the glioma cell, and ${k}_{2}^{\left(1\right)}$ corresponds to the rate of reaction that the CAR T-cells become no longer active, for example, due to cell death or exhaustion [[28]

]. Hence ${k}_{3}^{\left(1\right)}$ stands for the successful killing process, and ${k}_{2}^{\left(1\right)}$ and ${k}_{-1}^{\left(1\right)}$ indicates that some stage in the killing process did not work properly. In addition to the interaction, the growth dynamics of the glioma cells and CAR T-cells are included in the model, which yields in the following system of differential equation. $$\stackrel{\u0307}{T}=pT\frac{C}{g+C}-\theta T-{k}_{1}^{\left(1\right)}CT+{k}_{-1}^{\left(1\right)}{I}_{1}+{k}_{3}^{\left(1\right)}{I}_{1}$$$$\stackrel{\u0307}{C}=\rho C(1-\frac{C}{K})-{k}_{1}^{\left(1\right)}CT+{k}_{-1}^{\left(1\right)}{I}_{1}+{k}_{2}^{\left(1\right)}{I}_{1}$$$$\stackrel{\u0307}{{I}_{1}}={k}_{1}^{\left(1\right)}CT-{k}_{-1}^{\left(1\right)}{I}_{1}-{k}_{3}^{\left(1\right)}{I}_{1}-{k}_{2}^{\left(1\right)}{I}_{1}.$$

Here, parameter $\rho $ is the maximal growth rate of the glioma cells assuming logistic growth; parameter $K$ is the maximal carrying capacity of the biological environment for the glioma cells; $p$ is the rate that the CAR T-cells accumulate in the region due to the presence of the tumor, in other words, CAR T-cell expansion rate; $g$ is a concentration of glioma cells that halves the maximum rate $p$; $\theta $ represents the death rate of the CAR T-cells. Lysed cancer cells and CAR T-cells will be formed at the end of this interaction, but since their dynamics are simply decaying, we do not include those equations in our model as in [

[7]

]. We reparameterize the equation as $$\stackrel{\u0307}{T}=pT\frac{C}{g+C}-\theta T-kCT+\alpha {I}_{1}$$$$\stackrel{\u0307}{C}=\rho C\left(1-\frac{C}{K}\right)-kCT+\beta {I}_{1}$$$$\stackrel{\u0307}{{I}_{1}}=kCT-\gamma {I}_{1},$$

where

$$k={k}_{1}^{\left(1\right)},\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\alpha ={k}_{-1}^{\left(1\right)}+{k}_{3}^{\left(1\right)}$$$$\beta ={k}_{-1}^{\left(1\right)}+{k}_{2}^{\left(1\right)},\phantom{\rule{1em}{0ex}}\gamma ={k}_{-1}^{\left(1\right)}+{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)}.$$

The new parameter $\alpha $ represents the rate of net inflow of CAR T-cells from the conjugate ${I}_{1}$, $\beta $ represents the inflow rate for glioma cells, and $\gamma $ stands for the net decay rate of the conjugate ${I}_{1}$ due to either detachment or lysed reaction. We will refer to Eqs. (1) or (2) as the

*slow reaction*model, since following the dynamics of the conjugate ${I}_{1}$ allows the interaction to be in a comparable time scale as the other dynamics, in contrast to the assumption that will be made in the following model.#### 2.1.1 One-to-one binding model with fast reaction

The model proposed in [

Then, Eq. (2) reduces to the Kuznetsov et al. (1994) model as

where we define

We will refer to this model as the

[7]

] further reduces Eq. (1) by a separation of time scale between the dynamics of glioma cells and CAR T-cells compared to the dynamics of the conjugates ${I}_{1}$. Assuming that the dynamics of conjugate ${I}_{1}$ reaches equilibrium quickly, we can take ${\stackrel{\u0307}{I}}_{1}=0$, and we have ${k}_{1}^{\left(1\right)}CT-{k}_{-1}^{\left(1\right)}{I}_{1}-{k}_{3}^{\left(1\right)}{I}_{1}-{k}_{2}^{\left(1\right)}{I}_{1}=0$, that is, $${I}_{1}=\frac{{k}_{1}^{\left(1\right)}}{{k}_{-1}^{\left(1\right)}+{k}_{3}^{\left(1\right)}+{k}_{2}^{\left(1\right)}}CT.$$

Then, Eq. (2) reduces to the Kuznetsov et al. (1994) model as

$$\stackrel{\u0307}{T}=pT\frac{C}{g+C}-\theta T-mCT$$$$\stackrel{\u0307}{C}=\rho C\left(1-\frac{C}{K}\right)-\ell CT,$$

where we define

$$\ell =\frac{{k}_{1}^{\left(1\right)}{k}_{3}^{\left(1\right)}}{{k}_{-1}^{\left(1\right)}+{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)}},\phantom{\rule{1em}{0ex}}m=\frac{{k}_{1}^{\left(1\right)}{k}_{2}^{\left(1\right)}}{{k}_{-1}^{\left(1\right)}+{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)}}.$$

(4)

We will refer to this model as the

*fast reaction*model of one CAR T-cell binding case, as opposed to the slow reaction model in Eq. (2). The fast reaction scenario can happen, for example, when the glioma cell is mixed with a large number of CAR T-cells or the receptor density levels of glioma cells are high, which may lead to a fast elimination of glioma cells.### 2.2 Two CAR T-cells bound to one glioma cell model

The following system, building up on the one-binding models in Eqs. ((2), (3)), will govern the population dynamics of the scenario where up to two CAR T-cells can interact and bind to one glioma cell. As in the diagram of Fig. 2 (middle), conjugate ${I}_{1}$ of one CAR T-cell and one glioma cell can interact with another CAR T-cell ($T$) and compose conjugate ${I}_{2}$ of two CAR T-cells and one glioma cell. The kinetic interaction rates are defined similarly as before, where ${k}_{1}^{\left(2\right)}$ and ${k}_{-1}^{\left(2\right)}$ will describe the binding and detachment rate of cells without damage, and ${k}_{i}^{\left(2\right)}$ (i=2,3,4) will be the rates of the interaction events. The system of equations can be written as follows.

This model follows the dynamics of ${I}_{1}$ and ${I}_{2}$ conjugates, so we will refer to this as a two-to-one or two binding slow reaction model. This model describes the scenario that if one CAR T-cell cannot kill a glioma cell immediately, then another CAR T-cell may come to react with the conjugate ${I}_{1}$ formed by previous reaction.

$$\stackrel{\u0307}{T}=pT\frac{C}{g+C}-\theta T-kCT+\alpha {I}_{1}\phantom{\rule{0.16667em}{0ex}}-{k}_{1}^{\left(2\right)}{I}_{1}T+({k}_{-1}^{\left(2\right)}+2{k}_{4}^{\left(2\right)}+{k}_{2}^{\left(2\right)}){I}_{2}$$$$\stackrel{\u0307}{C}=\rho C\left(1-\frac{C}{K}\right)-kCT+\beta {I}_{1}+({k}_{3}^{\left(2\right)}+{k}_{2}^{\left(2\right)}){I}_{2}$$$$\stackrel{\u0307}{{I}_{1}}=kCT-\gamma {I}_{1}-{k}_{1}^{\left(2\right)}{I}_{1}T+{k}_{-1}^{\left(2\right)}{I}_{2}$$$$\stackrel{\u0307}{{I}_{2}}={k}_{1}^{\left(2\right)}{I}_{1}T-({k}_{-1}^{\left(2\right)}+{k}_{4}^{\left(2\right)}+{k}_{3}^{\left(2\right)}+{k}_{2}^{\left(2\right)}){I}_{2}$$

This model follows the dynamics of ${I}_{1}$ and ${I}_{2}$ conjugates, so we will refer to this as a two-to-one or two binding slow reaction model. This model describes the scenario that if one CAR T-cell cannot kill a glioma cell immediately, then another CAR T-cell may come to react with the conjugate ${I}_{1}$ formed by previous reaction.

#### 2.2.1 Two-to-one binding model with fast reaction

Using a similar approach when we derived the Kuznetsov et al. (1994) model, the two-binding slow reaction system can be reduced to the following two-binding fast reaction system.

Note that we added the assumption that the detach rate is small, that is, ${k}_{-1}^{\left(1\right)}\approx {k}_{-1}^{\left(2\right)}\approx 0$, for simplicity. Derivation of model details are included in Appendix A.

$$\stackrel{\u0307}{T}=pT\frac{C}{g+C}-\theta T-\frac{{k}_{1}^{\left(1\right)}{k}_{2}^{\left(1\right)}}{{k}_{1}^{\left(2\right)}T+{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)}}CT-\phantom{\rule{0.16667em}{0ex}}\frac{{k}_{1}^{\left(1\right)}{k}_{1}^{\left(2\right)}({k}_{2}^{\left(2\right)}+2{k}_{3}^{\left(2\right)})}{({k}_{2}^{\left(2\right)}+{k}_{3}^{\left(2\right)}+{k}_{4}^{\left(2\right)})({k}_{1}^{\left(2\right)}T+{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)})}C{T}^{2}$$$$\stackrel{\u0307}{C}=\rho C(1-\frac{C}{K})-\frac{{k}_{1}^{\left(1\right)}{k}_{3}^{\left(1\right)}}{{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)}+{k}_{1}^{\left(2\right)}T}CT-\phantom{\rule{0.16667em}{0ex}}\frac{{k}_{4}^{\left(2\right)}{k}_{1}^{\left(2\right)}{k}_{1}^{\left(1\right)}}{({k}_{2}^{\left(2\right)}+{k}_{3}^{\left(2\right)}+{k}_{4}^{\left(2\right)})({k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)}+{k}_{1}^{\left(2\right)}T)}C{T}^{2}.$$

Note that we added the assumption that the detach rate is small, that is, ${k}_{-1}^{\left(1\right)}\approx {k}_{-1}^{\left(2\right)}\approx 0$, for simplicity. Derivation of model details are included in Appendix A.

### 2.3 Multiple CAR T-cells bound to one glioma cell model

Building up on the two-binding model, we can further extend the model to $n$ numbers of CAR T-cells binding to one glioma cell. We denote ${I}_{j}$ as the conjugate of $j$ CAR T-cells and one glioma cell. The notation ${k}_{i}^{\left(n\right)}$ denotes the rate of different reactions, where the superscript $\left(n\right)$ denotes the $n$-binding and the subscript $i$ denotes the $i$th reaction in the $n$-binding reaction. The governing equations of the $n$-binding model that includes the conjugates ${I}_{1}$, …, ${I}_{n}$ can be written as follows.

where ${\alpha}_{j}={\sum}_{i=1}^{j}i{k}_{i+2}^{\left(j\right)}+{k}_{-1}^{\left(j\right)}$, ${\beta}_{j}={\sum}_{i=1}^{j}{k}_{i+1}^{\left(j\right)}$, and ${\gamma}_{j}={\sum}_{i=2}^{j+2}{k}_{i}^{\left(j\right)}+{k}_{-1}^{\left(j\right)}$. The parameters and their biological meanings are summarized in Table 1.

$$\stackrel{\u0307}{T}=pT\frac{C}{g+C}-\theta T-{k}_{1}^{\left(1\right)}CT+\sum _{j=1}^{n}{\alpha}_{j}{I}_{j}-\sum _{j=1}^{n-1}{k}_{1}^{(j+1)}{I}_{j}T$$$$\stackrel{\u0307}{C}=\rho C(1-\frac{C}{K})-{k}_{1}^{\left(1\right)}CT+{k}_{-1}^{\left(1\right)}{I}_{1}+\sum _{j=1}^{n}{\beta}_{j}{I}_{j}\phantom{\rule{45.52458pt}{0ex}}$$$$\stackrel{\u0307}{{I}_{1}}={k}_{1}^{\left(1\right)}CT-({\gamma}_{1}+{k}_{1}^{\left(2\right)}T){I}_{1}+{k}_{-1}^{\left(2\right)}{I}_{2}\phantom{\rule{78.24507pt}{0ex}}$$$$\stackrel{\u0307}{{I}_{j}}={k}_{1}^{\left(j\right)}{I}_{j-1}T-({\gamma}_{j}+{k}_{1}^{(j+1)}T){I}_{j}+{k}_{-1}^{(j+1)}{I}_{j+1}\phantom{\rule{44.10185pt}{0ex}}$$$$\stackrel{\u0307}{{I}_{n}}={k}_{1}^{\left(n\right)}{I}_{n-1}T-{\gamma}_{n}{I}_{n}\phantom{\rule{142.26378pt}{0ex}}$$

where ${\alpha}_{j}={\sum}_{i=1}^{j}i{k}_{i+2}^{\left(j\right)}+{k}_{-1}^{\left(j\right)}$, ${\beta}_{j}={\sum}_{i=1}^{j}{k}_{i+1}^{\left(j\right)}$, and ${\gamma}_{j}={\sum}_{i=2}^{j+2}{k}_{i}^{\left(j\right)}+{k}_{-1}^{\left(j\right)}$. The parameters and their biological meanings are summarized in Table 1.

In the following sections, we will test various hypotheses on the reaction rates to reproduce the phenomena of saturation of CAR T-cell therapy efficacy when the antigen receptor density of cancer increases.

Table 1Model parameters and their biological interpretation.

Parameter | Biological meaning |
---|---|

$\rho $ | Proliferation rate of glioma cells |

$K$ | Carrying capacity of glioma cells |

$p$ | Rate of CAR T-cell expansion induced by glioma |

$\theta $ | Death rate of CAR T-cells |

$g$ | Steepness coefficient of CAR T-cell expansion |

${k}_{1}^{\left(n\right)}$ | Binding rate of CAR T-cell and conjugate ${I}_{n-1}$ |

${k}_{-1}^{\left(n\right)}$ | Detaching rate of CAR T-cell and conjugate ${I}_{n}$ |

${k}_{n+2}^{\left(n\right)}$ | Death rate of glioma cells from conjugate ${I}_{n}$ |

${k}_{i}^{\left(n\right)}{|}_{i=2}^{n+1}$ | Death rate of CAR T-cells from conjugate ${I}_{n}$ |

## 3. Stability analysis

In order to better understand the dynamics of CAR T-cell therapy with the slow and fast reaction models, we present stability analysis of the one-to-one binding model.

### 3.1 One-to-one binding slow reaction model

In this section, we study the stability of one-binding slow reaction model (2). There are four steady states $(T,C,{I}_{1})$ in the slow reaction model

where

and

Note that we order the points as ${C}_{1}<{C}_{2}$, so that ${C}_{1}$ is the state with a smaller cancer size. Among these four states, we are interested in $(0,K,0)$ and $({T}_{1},{C}_{1},\frac{k}{\gamma}{C}_{1}{T}_{1})$, which represent CAR T-cell treatment failure and success, respectively. The steady state $(0,0,0)$ is always an unstable saddle point (See Appendix B), thus it is not of our interest. We aim to find the parameter conditions that yield these two states, especially focusing on the CAR T-cell expansion rate induced by the presence of cancer, $p$.

$$(0,0,0),\phantom{\rule{1em}{0ex}}(0,K,0),\phantom{\rule{1em}{0ex}}({T}_{1},{C}_{1},\frac{k}{\gamma}{C}_{1}{T}_{1}),\phantom{\rule{1em}{0ex}}({T}_{2},{C}_{2},\frac{k}{\gamma}{C}_{2}{T}_{2}),$$

where

$${T}_{1}=\frac{\rho (1-\frac{{C}_{1}}{K})}{k-\frac{\beta k}{\gamma}},\phantom{\rule{1em}{0ex}}{T}_{2}=\frac{\rho (1-\frac{{C}_{2}}{K})}{k-\frac{\beta k}{\gamma}},$$

and

$${C}_{1}=\frac{(p-\theta -kg+\frac{\alpha kg}{\gamma})}{2(k-\frac{\alpha k}{\gamma})}\phantom{\rule{0.16667em}{0ex}}-\frac{\sqrt{{(p-\theta -kg+\frac{\alpha kg}{\gamma})}^{2}-4(\frac{\alpha k}{\gamma}-k)(-\theta g)}}{2(k-\frac{\alpha k}{\gamma})},$$$${C}_{2}=\frac{(p-\theta -kg+\frac{\alpha kg}{\gamma})}{2(k-\frac{\alpha k}{\gamma})}\phantom{\rule{0.16667em}{0ex}}+\frac{\sqrt{{(p-\theta -kg+\frac{\alpha kg}{\gamma})}^{2}-4(\frac{\alpha k}{\gamma}-k)(-\theta g)}}{2(k-\frac{\alpha k}{\gamma})}.$$

Note that we order the points as ${C}_{1}<{C}_{2}$, so that ${C}_{1}$ is the state with a smaller cancer size. Among these four states, we are interested in $(0,K,0)$ and $({T}_{1},{C}_{1},\frac{k}{\gamma}{C}_{1}{T}_{1})$, which represent CAR T-cell treatment failure and success, respectively. The steady state $(0,0,0)$ is always an unstable saddle point (See Appendix B), thus it is not of our interest. We aim to find the parameter conditions that yield these two states, especially focusing on the CAR T-cell expansion rate induced by the presence of cancer, $p$.

#### 3.1.1 CAR T-cell treatment failure

One of the steady states is $(0,K,0)$ that represents the tumor reaching its maximal capacity, or in other words, the CAR T-cell treatment failure. This steady state becomes stable if

or using the kinetic reaction parameters,

This condition implies that if the engineered CAR T-cell expansion rate $p$ is less than $\left(g+K\right)\left(k(1-\frac{\alpha}{\gamma})+\frac{\theta}{K}\right)$, the tumor reaches its maximum capacity; meanwhile, if the engineered CAR T-cells can expand enough so that $p\ge \left(g+K\right)\left(k(1-\frac{\alpha}{\gamma})+\frac{\theta}{K}\right)$, then the CAR T-cell therapy will prevent the cancer from growing to its maximum capacity. Therefore, this condition provides a minimal level of CAR T-cell expansion rate which can prevent cancer from growing to its maximal size.

$$p<\left(g+K\right)\left(k(1-\frac{\alpha}{\gamma})+\frac{\theta}{K}\right).$$

(8)

or using the kinetic reaction parameters,

$$p<\left(g+K\right)\left({k}_{1}^{\left(1\right)}\left(\frac{{k}_{2}^{\left(1\right)}}{{k}_{-1}^{\left(1\right)}+{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)}}\right)+\frac{\theta}{K}\right).$$

(9)

This condition implies that if the engineered CAR T-cell expansion rate $p$ is less than $\left(g+K\right)\left(k(1-\frac{\alpha}{\gamma})+\frac{\theta}{K}\right)$, the tumor reaches its maximum capacity; meanwhile, if the engineered CAR T-cells can expand enough so that $p\ge \left(g+K\right)\left(k(1-\frac{\alpha}{\gamma})+\frac{\theta}{K}\right)$, then the CAR T-cell therapy will prevent the cancer from growing to its maximum capacity. Therefore, this condition provides a minimal level of CAR T-cell expansion rate which can prevent cancer from growing to its maximal size.

Let us further examine the condition in Eqs. ((8), (9)). Observe that as the reaction rate ${k}_{3}^{\left(1\right)}$ increases, or equivalently, $\alpha /\gamma $ becomes closer to 1, the right-hand side of the inequality decreases. Thus, if the killing of cancer becomes more effective, the minimum expansion rate of the CAR T-cell required to prevent treatment failure becomes smaller. It means that a CAR T-cell with a high cancer killing rate can prevent the tumor from growing to its maximum size despite a relatively low expansion rate. On the other hand, as ${k}_{2}^{\left(1\right)}$ increases, or $\alpha /\gamma $ gets closer to 0, the right-hand-side of the inequality increases. Thus if the inactivation rate of the CAR T-cells from interaction increases, CAR T-cells will need a higher expansion rate to successfully treat the tumor. As for the parameter $K$, which stands for the capacity of cancer cells, the right-hand-side is an increasing function with respect to $K$ on the range of $K>\sqrt{g\theta /(k-k\alpha /\gamma )}$. This means that a tumor with a larger capacity will need a larger value of expansion rate $p$ to prevent CAR T-cell treatment failure.

#### 3.1.2 Potential CAR T-cell treatment success

Another steady state of our interest is the CAR T-cell treatment success state, $({T}_{1},{C}_{1},\frac{k}{\gamma}{C}_{1}{T}_{1})$. Then, the tumor is reduced to size ${C}_{1}$ compared to reaching its maximal capacity $K$. For this equilibrium point to exist, the following condition needs to be satisfied,

or

This condition provides a minimum level of CAR T-cell expansion rate that makes it possible to avoid treatment failure that results in tumor growing to its maximal capacity, although the treatment success will depend on other parameters and initial CAR T-cell dosage. When ${k}_{2}^{\left(1\right)}$ is large, the term ${k}_{2}^{\left(1\right)}/({k}_{-1}^{\left(1\right)}+{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)})$ will be closer to 1, and our lower bound for $p$ will approach ${(\sqrt{\theta}+\sqrt{g{k}_{1}^{\left(1\right)}})}^{2}$. This implies that with a less effective cancer killing rate, the CAR T-cell requires a higher expansion rate to achieve potential treatment success. Also if the reaction rate ${k}_{3}^{\left(1\right)}$ has a large value, our lower bound for $p$ will decrease to $\theta $. It follows that the interaction requires a lower expansion rate to achieve a potential treatment success with a more effective cancer killing rate. Detailed analysis can be found in Appendix B, and further discussion on the stable nonzero CAR T-cell population is in Section 5.

$${\left(\sqrt{\theta}+\sqrt{gk\left(1-\alpha /\gamma \right)}\right)}^{2}\le p.$$

(10)

or

$${\left(\sqrt{\theta}+\sqrt{g{k}_{1}^{\left(1\right)}\frac{{k}_{2}^{\left(1\right)}}{{k}_{-1}^{\left(1\right)}+{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)}}}\right)}^{2}\le p$$

(11)

This condition provides a minimum level of CAR T-cell expansion rate that makes it possible to avoid treatment failure that results in tumor growing to its maximal capacity, although the treatment success will depend on other parameters and initial CAR T-cell dosage. When ${k}_{2}^{\left(1\right)}$ is large, the term ${k}_{2}^{\left(1\right)}/({k}_{-1}^{\left(1\right)}+{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)})$ will be closer to 1, and our lower bound for $p$ will approach ${(\sqrt{\theta}+\sqrt{g{k}_{1}^{\left(1\right)}})}^{2}$. This implies that with a less effective cancer killing rate, the CAR T-cell requires a higher expansion rate to achieve potential treatment success. Also if the reaction rate ${k}_{3}^{\left(1\right)}$ has a large value, our lower bound for $p$ will decrease to $\theta $. It follows that the interaction requires a lower expansion rate to achieve a potential treatment success with a more effective cancer killing rate. Detailed analysis can be found in Appendix B, and further discussion on the stable nonzero CAR T-cell population is in Section 5.

### 3.2 One-to-one binding fast reaction model

The stability analysis of the one-to-one binding fast reaction model Eq. (3) is comparable to the slow reaction model. The analysis can be found in many literature including [

where

The steady state that the tumor growing to its maximum capacity, $(T,C)=(0,K)$, becomes stable if

This condition is comparable to Eq. (8), which provides the condition that the CAR T-cell expansion rate prevents the cancer from growing to its maximum. The other condition related to the equilibrium point $({T}_{1},{C}_{1})$, the CAR T-cell therapy success case, is

comparable to Eq. (10). The CAR T-cell needs to expand at least this level to have chances for successful CAR T-cell treatment. We remark that the other equilibrium points $(0,0)$ and $({T}_{2},{C}_{2})$ are always saddles, so that they are not of our interest.

7

, 29

]. Here, we briefly summarize it to compare it with the slow reaction model. Identical to the slow reaction model, there are four possible steady states $(T,C)$, $$(0,0),\phantom{\rule{1em}{0ex}}(0,K),\phantom{\rule{1em}{0ex}}({T}_{1},{C}_{1}),\phantom{\rule{1em}{0ex}}({T}_{2},{C}_{2}),$$

where

$${T}_{1,2}=\frac{\rho \left(1-\frac{{C}_{1,2}}{K}\right)}{\ell},$$$${C}_{1,2}=\frac{(p-\theta -mg)\pm \sqrt{{(p-\theta -mg)}^{2}-4mg\theta}}{2m}.$$

The steady state that the tumor growing to its maximum capacity, $(T,C)=(0,K)$, becomes stable if

$$p<\left(mK+\theta \right)\left(\frac{g}{K}+1\right).$$

(12)

This condition is comparable to Eq. (8), which provides the condition that the CAR T-cell expansion rate prevents the cancer from growing to its maximum. The other condition related to the equilibrium point $({T}_{1},{C}_{1})$, the CAR T-cell therapy success case, is

$${\left(\sqrt{\theta}+\sqrt{mg}\right)}^{2}\le p,$$

(13)

comparable to Eq. (10). The CAR T-cell needs to expand at least this level to have chances for successful CAR T-cell treatment. We remark that the other equilibrium points $(0,0)$ and $({T}_{2},{C}_{2})$ are always saddles, so that they are not of our interest.

## 4. Numerical study of glioma cells and CAR T-cells interaction

Our experimental data consist of multiple doses of CAR T-cells and multiple types of glioma cells based on the antigen receptor densities. Specifically, we look into the experiment data using IL13BB CAR T-cells [

[22]

]. The data includes glioma cells measured in a unitless cell index (CI) with low, medium, and high antigen receptor densities, where the latter is likely to have a better response to CAR T-cell treatment, although not strictly better. For each density level, the experiments were initialized with different mixture ratios between CAR T-cells and glioma cells as 1:5, 1:10, and 1:20 (See Fig. C.1). The number of initial CAR T-cells are determined to be proportional to the initial glioma cell number. For instance, for a 1:5 mixture ratio, the glioma cells are mixed with $0.2{C}_{0}$ (CI) $=0.2{C}_{0}\cdot 1{0}^{4}$ (cells) CAR T-cells, where ${C}_{0}$ is the initial glioma size in CI. Using this rich dataset, we numerically study the proposed multiple CAR T-cell binding model to find a better model that describes the experimental data.To numerically solve the models, we use the built-in ODE solver (ode23) in MATLAB, which is based on an explicit Runge–Kutta pair of Bogacki and Shampine [

30

, 31

]. The calibration of our model is done by minimizing the sum of squares norm of the error between the data and our model fit. To improve the accuracy and parameter identifiability of the calibration, we calibrate the parameters sequentially, first using the tumor growth data without treatment, then using the data with treatment. Using the no treatment data, we first calibrate the tumor growth parameter $\rho $ with fixed $K$ in $\stackrel{\u0307}{C}=\rho C(1-C/K)$. We obtain three sets of values of $\rho $ and $K$ given three different types of tumor considering antigen receptor densities — low, medium, and high. Afterward, we calibrate the remaining parameters for each density level while parameters $\rho $, $K$, and $\theta $ are fixed. See Tables in Appendix C for the full set of parameters. We comment that the calibrated values of CAR T-cell expansion rate $p$ change more than a magnitude when mixed with low to high antigen density glioma, where varying levels of T-cell expansion depending on the level of presented antigens has been observed [[32]

]. In the remaining section, we show that the slow reaction model Eq. (2) describes the data more accurately compared to the fast reaction model Eq. (3). Then, we compare different assumptions in the reaction rates to investigate the saturation of CAR T-cell treatment efficacy regarding the antigen density levels.### 4.1 Comparison between fast reaction and slow reaction models

The fast reaction model, for example, Eq. (3) does not describe the dynamics of the conjugates ${I}_{j}$, assuming that they reach the equilibrium state immediately. However, in reality, the reaction is not always fast enough as it takes time for the CAR T-cells to detect and infiltrate the glioma cells, which is also depicted in the experimental data [

[22]

]. Therefore, we expect the slow reaction model (2) to describe the data more accurately. Our hypothesis is confirmed in Fig. 5, where we compare the accuracy of the slow reaction model Eq. (2) and the fast reaction model Eq. (3). The error is computed as the sum of squares differences between the data and the calibrated model fit. In these plots, the errors from the three sets of data with the same receptor density (low and high) and CAR T-cell ratio (1:5, 1:10, 1:20) are marked as black crosses, and the bar shows the average of the errors. We observe that the slow reaction model results in significantly smaller errors for all cases. Let us comment on each dataset more closely. Fig. 3 compares the accuracy of slow binding and fast reaction models calibrated to the

*low*receptor density glioma data. We observe that the slow reaction model does have a better fitting than the fast reaction model. The top row shows the case of initial ratio 1:5, where we can see that the data points are closer to the calibrated model curve especially near $t=0.5$ and $t=1.5$ using the slow reaction model. Similarly in the second row that corresponds to the 1:10 ratio, the fitting of the slow reaction model is significantly better at capturing the saturating tail at the later time points after $t=3$. The last row shows the result of initial ratio 1:20, which is the case with the smallest CAR T-cell dosage. Again, the slow reaction model better describes the non-monotonic data. In all three cases of low antigen receptor density glioma experiments, we observe that the slow reaction model is more accurate, especially when the concavity is nonzero. The reason behind this result is that the reaction speed is slower when the antigen receptor density is low, which makes the slow reaction model more appropriate than the fast reaction model.On the other hand, we hypothesize that the fast reaction model could be more appropriate to fit the experiments with

*high*antigen receptor density glioma. Since the high receptor density glioma cells have more receptors for the CAR T-cells to bind, the reactions can be more efficient. We observe that from the data plotted in Fig. 4 that the tumor size decay faster compared to the low receptor density case. This is the case especially when the glioma cells are mixed with a large number of CAR T-cells, that is, the 1:5 ratio case among our experiments. The top row of Fig. 4 shows the 1:5 ratio mixture, where the CAR T-cells eliminate the glioma cells most rapidly. The model fits using the slow and fast reaction models both look accurate, and we confirm that this is the case that the fast reaction model can accurately capture the data. No big difference can be seen between the slow and fast reaction model fits, although the error is smaller using the slow reaction model (see Fig. 5(b)). Nonetheless, in case of lower dosages of CAR T-cells, for instance, 1:20 mixture, the slow reaction model improves the accuracy by more than 75%, being able to capture the initial bump of the data. Finally, Akaike information criterion (AIC) [33

, 34

] that quantifies the quality of model fitness is computed in Table 2. The AIC values of slow reaction model is smaller than the fast reaction model across all datasets thus confirms the better model fitness of the slow reaction model.We conclude that the slow reaction model is more appropriate to be considered in general, since the model fit is more accurate compared to the fast reaction model, and moreover, it agrees with the experimental observation that showed the reaction between the glioma cells and the CAR T-cells not being immediate. This happens especially when the glioma cells have low antigen receptor density and CAR T-cell numbers are small. In addition to the model fit comparison, we argue that the slow reaction model can either match the fast reaction model or deviate from it based on the parameter values. We remark that the slow reaction model has four interaction parameters that defines two parameters of fast reaction model, $\ell $ and $m$, presented in Eq. (4). Therefore there are multiple sets of parameters of slow reaction model that reduces to an identical fast reaction model. Fig. 6 shows such an example. The bottom figures are computed from two distinct parameter sets of the slow reaction model that yields identical $\ell $ and $m$ values, $\ell =2.785$ and $m=0.0223$, for the fast reaction model. Large values of ${k}_{3}^{\left(1\right)}$ indicate fast reaction, for example, when ${k}_{3}^{\left(1\right)}=100$, the ${I}_{1}$ conjugate dynamics are trivial and the slow reaction model agrees the fast reaction model. However, when ${k}_{3}^{\left(1\right)}=0.01$, we observe the ${I}_{1}$ conjugates number increasing and the slow reaction model show distinct result from the fast reaction model. While the fast reaction model shows the number of glioma cells declining, the slow reaction model with a different choice of ${k}_{3}^{\left(1\right)}$ changes the glioma cell dynamics from decreasing to increasing. This example shows the richer dynamics that the slow reaction model contains. However, the slow reaction model has a drawback that more uncertainty is present in the estimated parameters due to its larger number of parameters, thus interpretation of the fitted parameters needs to be done carefully, and it is our future work to study parameter identifiability in the proposed model when additional data of CAR T-cells and conjugates ${I}_{i}$ are available.

Table 2Comparison of model fitness using Akaike information criterion (AIC)

33

, 34

between fast reaction model Eq. 3 (Fast) and the slow reaction model Eq. 2 (Slow). Mean of AIC values among all datasets and in the parenthesis, minimal and maximal values are shown. The slow reaction model show better model fitness with lower values of AIC across all datasets.Low | Medium | High | |
---|---|---|---|

Fast | −300.0 | −346.9 | −298.1 |

(−337,-285) | (−378,-304) | (−385,-242) | |

Slow | −495.8 | −455.8 | −453.0 |

(−540,-438) | (−590,-340) | (−588,-284) |

### 4.2 Simulation of multiple-binding slow reaction models

Having compared the slow and fast reaction models in the previous section, we now study the multiple CAR T-cells binding model Eq. (7). In particular, we focus on the experimental results of glioma cells with different levels of antigen receptor density: low, medium, and high. When mixed with the same number of CAR T-cells, glioma cells with low antigen receptor density respond less than the glioma cells with medium antigen receptor density. In other words, the decay rate of glioma cell is positively correlated with the antigen receptor density in general. However, an interesting observation was made in [

[22]

] that the effectiveness CAR T-cell therapy saturated after a certain receptor density level. As shown in Fig. 7, the high receptor density cells did not respond better than the medium receptor density glioma cells, or responded rather worse. We aim to study these phenomena using the multiple CAR T-cells binding models Eq. (7). Assuming that the antigen receptor density levels of glioma cells determines the number of CAR T-cells that can bind to a single glioma cell, we associate the one CAR T-cell binding model ($n=1$) to the low density receptor glioma, and multiple CAR T-cells binding model up to the five binding ($n=5$) to higher density receptor levels of glioma. One problem we face for the multiple CAR T-cells binding models is that the number of interaction parameters ${k}_{i}^{\left(j\right)}$ increases as the number of binding increases. Since we do not have data for the dynamics of any conjugates ${I}_{j}$ for $j=1,\dots ,n$, it is difficult to estimate parameters from the data. Therefore, using the parameter set estimated for the one-binding model, we test the following two different hypotheses describing the relationships between the reaction rates for the multiple binding models with more than one binding.- •Hypothesis 1: Assume that the reaction rates are
*uniform*across conjugates ${I}_{j}$.$\u2022$ The CAR T-cell attaches to the glioma cell or the conjugate ${I}_{j}$ with an identical rate, i.e.,$${k}_{1}^{\left(1\right)}={k}_{1}^{\left(j\right)},\phantom{\rule{1em}{0ex}}\text{forall}j=2,3,\dots ,n.$$(14)$\u2022$ The reactions that the glioma cell dies from the conjugate ${I}_{j}$ have the same rates across number of bindings, i.e.,$${k}_{1+2}^{\left(1\right)}={k}_{j+2}^{\left(j\right)},\phantom{\rule{1em}{0ex}}\text{forall}j=2,3,\dots ,n.$$(15)$\u2022$ The reactions that glioma cells’ survive also have the same reaction rates with equal chances of CAR T-cells dying i.e,$${k}_{j}^{\left(n\right)}=\frac{\left(\genfrac{}{}{0ex}{}{n}{j-1}\right)}{2{n}^{2}-6n+7}{k}_{2}^{\left(1\right)},\phantom{\rule{0.16667em}{0ex}}\phantom{\rule{0.16667em}{0ex}}\text{for}j=2,3,\dots ,n+1,$$(16)

so that ${k}_{2}^{\left(1\right)}={\sum}_{i=2}^{n+1}{k}_{i}^{\left(n\right)}$.Note that $\left(\genfrac{}{}{0ex}{}{n}{j-1}\right)$ denotes the number of cases when $j-1$ CAR T-cells die after the detachment of ${I}_{n}$, while $2{n}^{2}-6n+7$ denotes the number of all cases when at least one CAR T-cell dies (in other words, glioma cell survives) from the detachment of ${I}_{n}$. - •Hypothesis 2: Instead of assuming that the reaction rates are uniform across the conjugates ${I}_{j}$, in this hypothesis, we assume that the reaction rates are
*non-uniform*, either increasing or decreasing, and saturate after a certain number of CAR T-cells bind.$\u2022$ The reaction rates of a new CAR T-cell attaching to the conjugate ${I}_{j}$ decrease geometrically, i.e.,$${k}_{1}^{\left(j\right)}=\frac{{k}_{1}^{\left(1\right)}}{{M}^{j-1}},\phantom{\rule{1em}{0ex}}\text{forall}j=2,3,\dots ,n.$$(17)

for some positive integer $M\ge 1$.$\u2022$ The glioma cells’ death rate increases geometrically as the number of bindings increases, but saturates after three CAR T-cells binding, i.e.,$$\phantom{\rule{0.16667em}{0ex}}{k}_{4}^{\left(2\right)}={k}_{3}^{\left(1\right)}L,\phantom{\rule{0.16667em}{0ex}}\phantom{\rule{0.16667em}{0ex}}\phantom{\rule{0.16667em}{0ex}}{k}_{5}^{\left(3\right)}={k}_{3}^{\left(1\right)}{L}^{2},\phantom{\rule{0.16667em}{0ex}}{k}_{j+2}^{\left(j\right)}={k}_{5}^{\left(3\right)},\phantom{\rule{1em}{0ex}}\text{forall}j=4,5,\dots ,n.\phantom{\rule{42.67912pt}{0ex}}$$(18)

for some positive integer $L\ge 1$.We observe that the parameter $L$ directly leads to efficiency of cancer killing rate. With a bigger $L$, glioma cells’ death rate gets bigger for the first 3 bindings, which is more likely to bring a successful treatment. To estimate the parameter L directly from experiments, it can be determined by comparing the glioma cells’ death rates between $n+1$-binding interaction and $n$-binding interaction.The reactions that glioma cells’ survive are distributed as in hypothesis 1, Eq. (16).

We compare the two hypotheses to find which assumption yields the experimental data we have considering different levels of antigen receptor densities in glioma cells. In particular, we compare the final tumor size after 3 days of CAR T-cell treatment, i.e. $C\left(t\right)$ at $t=3$. In the top row of Fig. 8, we show the results testing hypothesis 1. The results do not match what we observe in the experiments, as the final tumor size increases from one-binding to five-binding models. We presume that this is due to the relative reaction rate of cancer death decreasing as the number of reaction increases. Hence, by simply considering uniform reaction rates for one to multiple CAR T-cell conjugates as in hypothesis 1, we cannot recover the dynamics of glioma cells responding better in higher antigen receptor density levels.

Nevertheless, results from hypothesis 2 show distinctive outcomes from hypothesis 1. The simulation results taking different values of $M$ and $L$ are shown in the bottom row of Fig. 8, Fig. 9. In particular, Fig. 9 shows the result of $M=1,\phantom{\rule{0.16667em}{0ex}}L=2$, and $M=2,\phantom{\rule{0.16667em}{0ex}}L=2$. It turns out that we no longer have the increase of tumor size as $n$ increases, instead, we have a decrease for both low and high densities. We can further observe that the reaction saturates after the numbers of binding become more than three, as the final tumor size remains very small eventually. However, if we consider $M=2,\phantom{\rule{0.16667em}{0ex}}L=1$, the tumor size increases as $n$ increases, and the results again deviate from the experimental data as shown in the bottom row of Fig. 8. Thus, we conclude that hypothesis 2 on the reaction rates with $L>1$ agrees with the experimental data shown in Fig. 7 (c,d) regarding the decay and saturation of tumor size with respect to increasing antigen density receptor levels. It makes sense that $L$ is the critical parameter for this result, since it is the parameter that makes the death rate of the glioma cells increase as $n$ increases. In addition, parameter sensitivity analysis [

[35]

] is included in Appendix D that shows larger Pearson correlation coefficient of $L$ compared to $M$ to the final tumor volume despite various levels of magnitude.## 5. Summary and future work

Here we developed an ODE model extending [

[7]

]. While the original model considers one conjugate ${I}_{1}$ of one glioma cell and one CAR T-cell, and assumes that the dynamics of ${I}_{1}$ conjugate is in equilibrium, our model considers multiple conjugates ${I}_{j}$, $j=1,\dots ,n$ with more than one CAR T-cells binding to the glioma cell, and follow their dynamics. We denote the original model as one binding fast reaction model, and our models as $n$ binding slow reaction models.First, we study the stability of the one-binding slow reaction ODE system, and compare it with the fast reaction model. We obtain similar equilibrium states, but in terms of the parameter of the slow reaction model. Also, we derive the conditions, especially regarding the range of the CAR T-cell expansion rate parameter $p$, (i) the minimum level of CAR T-cell expansion rate that provides the possibility that the CAR T-cell treatment can be successful, and (ii) the level of CAR T-cell expansion rate that guarantees the escape from reaching the maximum tumor size.

Using our model, we then compare the fast and slow reaction model numerically, and show that the slow reaction model describes the experimental data more accurately, especially when the glioma antigen receptor density is low and/or mixed with relatively small numbers of CAR T-cells. In addition, we use the one-binding to five-binding slow reaction models to simulate low to high receptor density glioma cells, and study the assumption in the reaction rates that yield desired outcome, the decay and saturation of tumor size with respect to the antigen density levels. We come up with two different hypotheses to describe their connections, where the first hypothesis considers homogeneous rates across different numbers of CAR T-cell binding conjugates. More precisely, we assume identical rates for reactions involving glioma cells dying from the conjugates ${I}_{j}$. We show that the second hypothesis, where we consider non-homogeneous reaction rates, in particular, increasing tumor killing rates as the number of bindings increases, reflects the reduced but saturated tumor size.

One of our future works is to calibrate the reaction rate parameters from experimental data, and verify whether our assumptions in hypothesis 2 on the reaction rates of multiple binding model is valid. Currently, our model is calibrated to the time series data of glioma cells, and additional CAR T-cell dynamics data and experiments can be used to infer the reaction rates from data. We also hope to validate the effectiveness of the slow reaction model to in vivo data where we expect the interaction between the glioma cells and CAR T-cells to be slower than our current in vitro data. Different CAR T-cell designs lead to different affinities and cytotoxicities [

36

, 37

, 38

], thus we hope to estimate and compare the kinetic reaction rates for different immunological synapse designs of CAR T-cells. In addition, parameter identifiability study and strategies to find universal parameters that can be fixed across patients are our future work as well.A limitation of our model is that it is fitted to the in vitro data that does not include the human immune system or tumor microenvironment. We propose to study and model the interaction of glioma cells and CAR T-cells with other immune cells in the future. Modeling the interaction with endogenously produced immune cells including non CAR T-cells, natural killer cells, and antigen presenting cells that include macrophages, dendritic cells, and B cells, along with cytokines may help understand how and why CAR T-cell conjugates are created as well as the role of other immune cells in the immune system in this process. This is also related to further investigating the persisting nonzero CAR T-cell population in case of treatment success predicted by our model in the stability analysis. Such a stable CAR T-cell population has been shown to be achieved in the clinics treating blood cancer [

39

, 40

]. In particular, [[41]

] studies the persistence of CD4+ CAR T-cell population over 10 years after therapy treating leukemia, and finds that CD8+ and CD4+ CAR T-cells are critical for elimination of leukemic cells at early stage and in long-term, respectively [[42]

]. This indicates the importance of considering different populations of T-cells to better achieve durable remission, and this is also our future work.## CRediT authorship contribution statement

**Runpeng Li:**Analysis, Investigation, Visualization, Software, Writing – original draft.

**Prativa Sahoo:**Data curation, Interpretation, Writing – review & editing.

**Dongrui Wang:**Data curation.

**Qixuan Wang:**Methodology, Writing – review & editing.

**Christine E. Brown:**Data curation, Writing – review & editing.

**Russell C. Rockne:**Data curation, Interpretation, Writing – review & editing.

**Heyrim Cho:**Conceptualization, Methodology, Analysis, Visualization, Writing – review & editing.

## Data availability

I have shared the link to my code at the Attached File step.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

Research reported in this publication was supported by the California Institute for Regenerative Medicine (CLIN2-10248, CB), the National Cancer Institute of the National Institutes of Health under grant numbers R01CA236500 (CB), P30CA033572 (S. Rosen). 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. We want to thank Renate Starr and Jim O’Hearn for help with the manuscript and figures. This manuscript is dedicated to the memory of Xin (Cindy) Yang and her contributions on this study.

## Appendix A. Computation of two-binding model with fast reaction

Our assumption in this case is that

which means that the two conjugates decompose immediately right before the other CAR T-cells come to react with the conjugates themselves. By rearranging Eq. (19), we obtain

and by rearranging Eq. (20), we obtain

where we further assume that ${k}_{-1}^{\left(1\right)}\approx 0$. Now we look at $\stackrel{\u0307}{T}$

by substituting (21) into our equation and simplify it, we obtain

For $\stackrel{\u0307}{C}$, by using ((21), (22)), we have

which can be rewritten as

Therefore, ((23), (24)) give our two-binding model with fast reaction.

$$\stackrel{\u0307}{{I}_{1}}={k}_{1}^{\left(1\right)}CT-({k}_{-1}^{\left(1\right)}+{k}_{3}^{\left(1\right)}+{k}_{2}^{\left(1\right)}){I}_{1}-{k}_{1}^{\left(2\right)}{I}_{1}T+\phantom{\rule{0.16667em}{0ex}}{k}_{-1}^{\left(2\right)}{I}_{2}=0$$$$\stackrel{\u0307}{{I}_{2}}={k}_{1}^{\left(2\right)}{I}_{1}T-({k}_{-1}^{\left(2\right)}+{k}_{4}^{\left(2\right)}+{k}_{3}^{\left(2\right)}+{k}_{2}^{\left(2\right)}){I}_{2}=0$$

which means that the two conjugates decompose immediately right before the other CAR T-cells come to react with the conjugates themselves. By rearranging Eq. (19), we obtain

$${I}_{1}=\frac{{k}_{1}^{\left(1\right)}CT+{k}_{-1}^{\left(2\right)}{I}_{2}}{{k}_{-1}^{\left(1\right)}+{k}_{3}^{\left(1\right)}+{k}_{2}^{\left(1\right)}+{k}_{1}^{\left(2\right)}T}=\frac{{k}_{1}^{\left(1\right)}CT}{{k}_{3}^{\left(1\right)}+{k}_{2}^{\left(1\right)}+{k}_{1}^{\left(2\right)}T}$$

(21)

and by rearranging Eq. (20), we obtain

$${I}_{2}=\frac{{k}_{1}^{\left(2\right)}{I}_{1}T}{{k}_{-1}^{\left(2\right)}+{k}_{4}^{\left(2\right)}+{k}_{3}^{\left(2\right)}+{k}_{2}^{\left(2\right)}}=\frac{{k}_{1}^{\left(2\right)}{I}_{1}T}{{k}_{4}^{\left(2\right)}+{k}_{3}^{\left(2\right)}+{k}_{2}^{\left(2\right)}}$$

(22)

where we further assume that ${k}_{-1}^{\left(1\right)}\approx 0$. Now we look at $\stackrel{\u0307}{T}$

$$\stackrel{\u0307}{T}=pT\frac{C}{g+C}-\theta T-{k}_{1}^{\left(1\right)}CT+{k}_{3}^{\left(1\right)}{I}_{1}-\phantom{\rule{0.16667em}{0ex}}{k}_{1}^{\left(2\right)}{I}_{1}T+(2{k}_{4}^{\left(2\right)}+{k}_{3}^{\left(2\right)}){I}_{2}=pT\frac{C}{g+C}-\theta T-{k}_{1}^{\left(1\right)}CT+{k}_{3}^{\left(1\right)}{I}_{1}+\phantom{\rule{0.16667em}{0ex}}\frac{{k}_{4}^{\left(2\right)}{k}_{1}^{\left(2\right)}-{k}_{1}^{\left(2\right)}{k}_{3}^{\left(2\right)}}{{k}_{2}^{\left(2\right)}+{k}_{3}^{\left(2\right)}+{k}_{4}^{\left(2\right)}}T{I}_{1}$$

by substituting (21) into our equation and simplify it, we obtain

$$\stackrel{\u0307}{T}=pT\frac{C}{g+C}-\theta T-\frac{{k}_{1}^{\left(1\right)}{k}_{2}^{\left(1\right)}}{{k}_{1}^{\left(2\right)}T+{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)}}CT-\phantom{\rule{0.16667em}{0ex}}\frac{{k}_{1}^{\left(1\right)}{k}_{1}^{\left(2\right)}({k}_{2}^{\left(2\right)}+2{k}_{3}^{\left(2\right)})}{({k}_{2}^{\left(2\right)}+{k}_{3}^{\left(2\right)}+{k}_{4}^{\left(2\right)})({k}_{1}^{\left(2\right)}T+{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)})}C{T}^{2}$$

(23)

For $\stackrel{\u0307}{C}$, by using ((21), (22)), we have

$$\stackrel{\u0307}{C}=\rho (1-\frac{C}{K})-{k}_{1}^{\left(1\right)}CT+{k}_{2}^{\left(1\right)}{I}_{1}+({k}_{3}^{\left(2\right)}+{k}_{2}^{\left(2\right)}){I}_{2}=\rho C(1-\frac{C}{K})-{k}_{1}^{\left(1\right)}CT+{k}_{2}^{\left(1\right)}{I}_{1}+\phantom{\rule{0.16667em}{0ex}}({k}_{3}^{\left(2\right)}+{k}_{2}^{\left(2\right)})\frac{{k}_{1}^{\left(2\right)}{I}_{1}T}{{k}_{4}^{\left(2\right)}+{k}_{3}^{\left(2\right)}+{k}_{2}^{\left(3\right)}}=\rho C(1-\frac{C}{K})-{k}_{1}^{\left(1\right)}CT+\phantom{\rule{0.16667em}{0ex}}({k}_{2}^{\left(1\right)}+\frac{{k}_{1}^{\left(2\right)}({k}_{2}^{\left(2\right)}+{k}_{3}^{\left(2\right)})}{{k}_{4}^{\left(2\right)}+{k}_{3}^{\left(2\right)}+{k}_{2}^{\left(2\right)}}T){I}_{1}=\rho C(1-\frac{C}{K})-{k}_{1}^{\left(1\right)}CT+\phantom{\rule{0.16667em}{0ex}}({k}_{2}^{\left(1\right)}+\frac{{k}_{1}^{\left(2\right)}({k}_{2}^{\left(2\right)}+{k}_{3}^{\left(2\right)})}{{k}_{4}^{\left(2\right)}+{k}_{3}^{\left(2\right)}+{k}_{2}^{\left(2\right)}}T)\frac{{k}_{1}^{\left(2\right)}CT}{{k}_{3}^{\left(1\right)}+{k}_{2}^{\left(1\right)}+{k}_{1}^{\left(2\right)}T}$$

which can be rewritten as

$$\stackrel{\u0307}{C}=\rho C(1-\frac{C}{K})-\frac{{k}_{1}^{\left(1\right)}{k}_{3}^{\left(1\right)}}{{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)}+{k}_{1}^{\left(2\right)}T}CT-\phantom{\rule{0.16667em}{0ex}}\frac{{k}_{4}^{\left(2\right)}{k}_{1}^{\left(2\right)}{k}_{1}^{\left(1\right)}}{({k}_{2}^{\left(2\right)}+{k}_{3}^{\left(2\right)}+{k}_{4}^{\left(2\right)})({k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)}+{k}_{1}^{\left(2\right)}T)}C{T}^{2}$$

(24)

Therefore, ((23), (24)) give our two-binding model with fast reaction.

## Appendix B. Stability analysis of one-binding model with slow reaction

The equilibrium points can be computed by:

Observe that because of $kCT-\gamma {I}_{1}=0$, this ODE system have the equilibrium points stated in Section 3.1. More precisely, using the relation $\frac{k}{\gamma}CT={I}_{1}$, we can rewrite the three dimensional ODE system into a two dimensional ODE system:

The first equation gives us

and

If we let $m\u2254k-\frac{\alpha k}{\gamma}$ and $\ell \u2254k-\frac{\beta k}{\gamma}$, these agree with $({T}_{1},{C}_{1})$ and $({T}_{2},{C}_{2})$ obtained with the fast reaction model. Therefore, again we have 4 equilibrium points

We then notice the Jacobian matrix of the ODE system Eq. (2) is

Let us examine each equilibrium points and compare them with the fast reaction model.

$$pT\frac{C}{g+C}-\theta T-kCT+\alpha {I}_{1}=0$$$$\rho C(1-\frac{C}{K})-kCT+\beta {I}_{1}=0$$$$kCT-\gamma {I}_{1}=0.$$

Observe that because of $kCT-\gamma {I}_{1}=0$, this ODE system have the equilibrium points stated in Section 3.1. More precisely, using the relation $\frac{k}{\gamma}CT={I}_{1}$, we can rewrite the three dimensional ODE system into a two dimensional ODE system:

$$pT\frac{C}{g+C}-\theta T-kCT+\frac{\alpha k}{\gamma}CT=0\rho C(1-\frac{C}{K})-kCT+\frac{\beta k}{\gamma}CT=0$$

The first equation gives us

$${T}_{1,2}=\frac{\rho (1-\frac{{C}_{1,2}}{K})}{k-\frac{\beta k}{\gamma}},$$

and

$${C}_{1,2}=\frac{(p-\theta -kg+\frac{\alpha kg}{\gamma})}{2(k-\frac{\alpha k}{\gamma})}\phantom{\rule{0.16667em}{0ex}}\pm \frac{\sqrt{{(p-\theta -kg+\frac{\alpha kg}{\gamma})}^{2}-4(\frac{\alpha k}{\gamma}-k)(-\theta g)}}{2(k-\frac{\alpha k}{\gamma})}.$$

If we let $m\u2254k-\frac{\alpha k}{\gamma}$ and $\ell \u2254k-\frac{\beta k}{\gamma}$, these agree with $({T}_{1},{C}_{1})$ and $({T}_{2},{C}_{2})$ obtained with the fast reaction model. Therefore, again we have 4 equilibrium points

$$(0,0,0),\phantom{\rule{1em}{0ex}}(0,K,0),\phantom{\rule{1em}{0ex}}({T}_{1},{C}_{1},\frac{k}{\gamma}{C}_{1}{T}_{1}),\phantom{\rule{1em}{0ex}}({T}_{2},{C}_{2},\frac{k}{\gamma}{C}_{2}{T}_{2})$$

We then notice the Jacobian matrix of the ODE system Eq. (2) is

$$J(T,C,{I}_{1})=\left(\begin{array}{ccc}\hfill p\frac{C}{g+C}-\theta -kC\hfill & \hfill \frac{pgT}{{(g+C)}^{2}}-kT\hfill & \hfill \alpha \hfill \\ \hfill -kC\hfill & \hfill \rho -2\rho \frac{C}{K}-kT\hfill & \hfill \beta \hfill \\ \hfill kC\hfill & \hfill kT\hfill & \hfill -\gamma \hfill \end{array}\right)$$

Let us examine each equilibrium points and compare them with the fast reaction model.

1. The equilibrium point $(T,C,{I}_{1})=(0,0,0)$ is the tumor-free and CAR T-cellfree case. The Jacobian matrix becomes

The corresponding characteristic polynomial will be

Because this is an upper triangular matrix, the eigenvalues can be observed easily, i.e. ${\lambda}_{1}=-\theta ,{\lambda}_{2}=\rho ,{\lambda}_{3}=-\gamma $. We conclude that $(0,0,0)$ is an unstable saddle point.

$$J(0,0,0)=\left(\begin{array}{ccc}\hfill -\theta \hfill & \hfill 0\hfill & \hfill \alpha \hfill \\ \hfill 0\hfill & \hfill \rho \hfill & \hfill \beta \hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -\gamma \hfill \end{array}\right)$$

The corresponding characteristic polynomial will be

$$det\left(\begin{array}{ccc}\hfill -\theta -\lambda \hfill & \hfill 0\hfill & \hfill \alpha \hfill \\ \hfill 0\hfill & \hfill \rho -\lambda \hfill & \hfill \beta \hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -\gamma -\lambda \hfill \end{array}\right)$$

Because this is an upper triangular matrix, the eigenvalues can be observed easily, i.e. ${\lambda}_{1}=-\theta ,{\lambda}_{2}=\rho ,{\lambda}_{3}=-\gamma $. We conclude that $(0,0,0)$ is an unstable saddle point.

2. The equilibrium point $(T,C,{I}_{1})=(0,K,0)$ is the tumor reaching the maximal capacity with no CAR T-cell surviving. The Jacobian matrix becomes

For simplicity, let $g\left(C\right)\u2254p\frac{C}{g+C}-\theta -kC$, so $g\left(K\right)=\frac{pK}{g+K}-\theta -Kk$. The corresponding characteristic polynomial will be

By setting $P\left(\lambda \right)=0$, we obtain the first eigenvalue ${\lambda}_{1}=-\rho $, and the following quadratic equation,

This gives two additional solutions

Note that the expression inside the square root is always positive since

Therefore, these two eigenvalues are always real. For this equilibrium point to be stable, we need the eigenvalues to be negative, i.e.

or equivalently,

which comes from $\gamma g\left(K\right)<-Kk\alpha $. If Eq. (25) is satisfied, the equilibrium point $(0,K,0)$ is stable. Therefore, this condition provides the minimal level of CAR T-cell expansion rate to avoid treatment failure.

$$J(0,K,0)=\left(\begin{array}{ccc}\hfill \frac{pK}{g+K}-\theta -Kk\hfill & \hfill 0\hfill & \hfill \alpha \hfill \\ \hfill -Kk\hfill & \hfill -a\hfill & \hfill \beta \hfill \\ \hfill Kk\hfill & \hfill 0\hfill & \hfill -\gamma \hfill \end{array}\right)$$

For simplicity, let $g\left(C\right)\u2254p\frac{C}{g+C}-\theta -kC$, so $g\left(K\right)=\frac{pK}{g+K}-\theta -Kk$. The corresponding characteristic polynomial will be

$$P\left(\lambda \right)=det\left(\begin{array}{ccc}\hfill g\left(K\right)-\lambda \hfill & \hfill 0\hfill & \hfill \alpha \hfill \\ \hfill -Kk\hfill & \hfill -a-\lambda \hfill & \hfill \beta \hfill \\ \hfill Kk\hfill & \hfill 0\hfill & \hfill -\gamma -\lambda \hfill \end{array}\right)=-Kk\alpha (-\rho -\lambda )+(-\gamma -\lambda )(gK-\lambda )(-\rho -\lambda )$$

By setting $P\left(\lambda \right)=0$, we obtain the first eigenvalue ${\lambda}_{1}=-\rho $, and the following quadratic equation,

$$-Kk\alpha +(-\gamma -\lambda )(g\left(K\right)-\lambda )=0.$$

This gives two additional solutions

$$\lambda =\frac{g\left(K\right)-\gamma \pm \sqrt{{(-g\left(K\right)+\gamma )}^{2}+4(Kk\alpha +\gamma g\left(K\right))}}{2}.$$

Note that the expression inside the square root is always positive since

$${(-g\left(K\right)+\gamma )}^{2}+4(Kk\alpha +\gamma g\left(K\right))={(\gamma +g\left(K\right))}^{2}+4Kk\alpha >0.$$

Therefore, these two eigenvalues are always real. For this equilibrium point to be stable, we need the eigenvalues to be negative, i.e.

$$g\left(K\right)-\gamma +\sqrt{{(-g\left(K\right)+\gamma )}^{2}+4(Kk\alpha +\gamma g\left(K\right))}<0,$$

or equivalently,

$$p<(\frac{g}{K}+1)(-\frac{k\alpha}{\gamma /K}+\theta +Kk)$$

(25)

which comes from $\gamma g\left(K\right)<-Kk\alpha $. If Eq. (25) is satisfied, the equilibrium point $(0,K,0)$ is stable. Therefore, this condition provides the minimal level of CAR T-cell expansion rate to avoid treatment failure.

3. The equilibrium points $({T}_{1,2},{C}_{1,2},\frac{k}{\gamma}{C}_{1,2}{T}_{1,2})$ include $({T}_{1},{C}_{1})$ that can be denoted as the CAR T-cell therapy success case by ordering the points as $0<{C}_{1}<{C}_{2}$ and ${T}_{1}>{T}_{2}>0$. For these equilibrium points to exist, we can obtain the two conditions, $p-\theta -kg+\frac{\alpha kg}{\gamma}\ge 0$, and ${(p-\theta -kg+\frac{\alpha kg}{\gamma})}^{2}+4(\frac{\alpha k}{\gamma}-k)\left(\theta g\right)\ge 0$, that reduces to

This condition provides a minimum level of CAR T-cell expansion rate that makes treatment success possible. To examine the stability, let us define $h\left(C\right)=-\theta +p\frac{C}{g+C}-kC+\frac{\alpha k}{\gamma}C$. The characteristic polynomial becomes

where the coefficients are computed as following

The stability of the equilibrium points will depend on the roots of this polynomial. Although it is difficult to analyze the condition from this polynomial, we know that for the CAR T-cell therapy success case $({T}_{1},{C}_{1})$, $D>0$ $\beta ={k}_{-1}^{\left(1\right)}+{k}_{2}^{\left(1\right)}$ $<\gamma ={k}_{-1}^{\left(1\right)}+{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)}$ and ${h}^{\prime}\left({C}_{1}\right)<0$, which provide us that $({T}_{1},{C}_{1})$ is either stable or a saddle.

$${\left(\sqrt{\theta}+\sqrt{g\left(-\alpha k/\gamma +k\right)}\right)}^{2}\le p.$$

(26)

This condition provides a minimum level of CAR T-cell expansion rate that makes treatment success possible. To examine the stability, let us define $h\left(C\right)=-\theta +p\frac{C}{g+C}-kC+\frac{\alpha k}{\gamma}C$. The characteristic polynomial becomes

$$P\left(\lambda \right)=\text{det}\left(\begin{array}{ccc}\hfill -\frac{\alpha k}{\gamma}{C}_{e}-\lambda \hfill & \hfill {h}^{\prime}\left({C}_{e}\right){T}_{e}-\frac{\alpha k}{\gamma}{T}_{e}\hfill & \hfill \alpha \hfill \\ \hfill -k{C}_{e}\hfill & \hfill \rho -2\frac{\rho}{K}{C}_{e}-k{T}_{e}-\lambda \hfill & \hfill \beta \hfill \\ \hfill k{C}_{e}\hfill & \hfill k{T}_{e}\hfill & \hfill -\gamma -\lambda \hfill \end{array}\right)=-{\lambda}^{3}+A{\lambda}^{2}+B\lambda +D$$

where the coefficients are computed as following

$$A=\rho -2\rho \frac{{C}_{e}}{K}-\frac{\alpha k}{\gamma}{C}_{e}-k{T}_{e}-\gamma \phantom{\rule{103.28363pt}{0ex}}$$$$B=(\rho -2\rho \frac{{C}_{e}}{K})\left(\gamma +\frac{\alpha k{C}_{e}}{\gamma}\right)+(\beta -\gamma )k{T}_{e}-k{h}^{\prime}\left({C}_{e}\right){T}_{e}{C}_{e}$$$$D=k{C}_{e}{h}^{\prime}\left({C}_{e}\right){T}_{e}(\beta -\gamma ).\phantom{\rule{139.9875pt}{0ex}}$$

The stability of the equilibrium points will depend on the roots of this polynomial. Although it is difficult to analyze the condition from this polynomial, we know that for the CAR T-cell therapy success case $({T}_{1},{C}_{1})$, $D>0$ $\beta ={k}_{-1}^{\left(1\right)}+{k}_{2}^{\left(1\right)}$ $<\gamma ={k}_{-1}^{\left(1\right)}+{k}_{2}^{\left(1\right)}+{k}_{3}^{\left(1\right)}$ and ${h}^{\prime}\left({C}_{1}\right)<0$, which provide us that $({T}_{1},{C}_{1})$ is either stable or a saddle.

Table C.1Parameters for the slow binding model Eq. 2 with low IL13R$\alpha $2 antigen density. The left column lists all the parameters that are considered in the calibration. The top row specifies different CAR T-cell ratios.

1:5 | 1:10 | 1:20 | |
---|---|---|---|

$\rho $ (day^{−1}) | 0.3953 | 0.3953 | 0.3953 |

$K$ (CI) | 6 | 6 | 6 |

$p$ (day^{−1}) | 3.2430 | 8.8245 | 7.0226 |

$g$ (CI) | 16.8454 | 17.0687 | 10.7698 |

${k}_{1}^{\left(1\right)}$ (day${}^{-1}\cdot $ CI^{−1}) | 6.1317 | 14.2485 | 9.9151 |

${k}_{2}^{\left(1\right)}$ (day${}^{-1}\cdot $ CI^{−1}) | 1.1658 | 0.8095 | 1.1463 |

${k}_{3}^{\left(1\right)}$ (day${}^{-1}\cdot $ CI^{−1}) | 9.1820 | 8.4086 | 12.8200 |

${k}_{-1}^{\left(1\right)}$ (day${}^{-1}\cdot $ CI^{−1}) | 0.1017 | 0.0658 | 0.7729 |

$\theta $ (day^{−1}) | 0.0412 | 0.0412 | 0.0412 |

Table C.2Parameters for the fast binding model Eq. 3 with low IL13R$\alpha $2 antigen density. The left column lists all the parameters that are considered in the calibration. The top row specifies different CAR T-cell ratios.

1:5 | 1:10 | 1:20 | |
---|---|---|---|

$\rho $ (day^{−1}) | 0.3953 | 0.3953 | 0.3953 |

$K$ (CI) | 6 | 6 | 6 |

$p$ (day^{−1}) | 0.4351 | 5.3280 | 7.9233 |

$g$ (CI) | 5.2676 | 11.5747 | 11.4983 |

$m$ (day${}^{-1}\cdot $ CI^{−1}) | 0.1179 | 0.4640 | 0.6260 |

$\ell $ (day${}^{-1}\cdot $ CI^{−1}) | 2.3818 | 3.0213 | 3.4057 |

$\theta $ (day^{−1}) | 0.0412 | 0.0412 | 0.0412 |

Table C.3Parameters for the slow binding model Eq. 2 with high IL13R$\alpha $2 antigen density. The left column lists all the parameters that are considered in the calibration. The top row specifies different CAR T-cell ratios.

1:5 | 1:10 | 1:20 | |
---|---|---|---|

$\rho $ (day^{−1}) | 0.2187 | 0.2187 | 0.2187 |

$K$ (CI) | 6 | 6 | 6 |

$p$ (day^{−1}) | 18.5305 | 16.5100 | 16.6778 |

$g$ (CI) | 8.0816 | 1.6240 | 2.3731 |

${k}_{1}^{\left(1\right)}$ (day${}^{-1}\cdot $ CI^{−1}) | 3.7075 | 3.7439 | 1.6125 |

${k}_{2}^{\left(1\right)}$ (day${}^{-1}\cdot $ CI^{−1}) | 4.3446 | 0.2965 | 2.0446 |

${k}_{3}^{\left(1\right)}$ (day${}^{-1}\cdot $ CI^{−1}) | 13.0198 | 1.8054 | 1.8582 |

${k}_{-1}^{\left(1\right)}$ (day${}^{-1}\cdot $ CI^{−1}) | 0.3386 | 0.0495 | 0.7008 |

$\theta $ (day^{−1}) | 0.0412 | 0.0412 | 0.0412 |

Table C.4Parameters for the fast binding model Eq. 3 with high IL13R$\alpha $2 antigen density. The left column lists all the parameters that are considered in the calibration. The top row specifies different CAR T-cell ratios.

1:5 | 1:10 | 1:20 | |
---|---|---|---|

$\rho $ (day^{−1}) | 0.2187 | 0.2187 | 0.2187 |

$K$ (CI) | 6 | 6 | 6 |

$p$ (day^{−1}) | 10.5097 | 16.7512 | 15.2942 |

$g$ (CI) | 5.5172 | 11.8066 | 17.9938 |

$m$ (day${}^{-1}\cdot $ CI^{−1}) | 0.5691 | 0.3766 | 0.0566 |

$\ell $ (day${}^{-1}\cdot $ CI^{−1}) | 1.8054 | 1.2394 | 1.4700 |

$\theta $ (day^{−1}) | 0.0412 | 0.0412 | 0.0412 |

## Appendix C. Data and parameters

The glioma growth data used to calibrate the models are shown in Fig. C.1.

The calibrated parameter values of one-to-one binding slow and fast reaction models to plot Fig. 3, Fig. 4 are provided in Table C.1, Table C.2, Table C.3, Table C.4. The tumor growth rate $\rho $ is calibrated with the tumor data without CAR T-cell treatment, while $K$ is fixed as $6$ CI $\approx 6\cdot 1{0}^{4}$ cells. The other parameters are calibrated with tumor data with CAR T-cell treatment while $\theta $ is fixed [

[7]

] to reduce uncertainty in parameter estimations.## Appendix D. Parameter sensitivity analysis

In Fig. D.1, we list partial correlation coefficients for all parameters in the five-binding model. In Table D.1, we include the partial correlation coefficients for parameters $M$ and $L$ to final tumor volume.

Table D.1Partial correlation coefficients of parameters $M$ and $L$ in multiple binding models 7 and hypothesis 2 in section 4.2. Presented cases of densities and CAR T-cell ratios are from Figs. 8–9. Both $M$ and $L$ are taken to be 2. The mean of partial correlation coefficients with the range of minimums and maximums across 2 bindings to 5 bindings in parentheses are shown.

$M$ | $L$ | |
---|---|---|

Low density, 1:10 | 0.04286 | −0.3575 |

(0.03768, 0.04676) | (−0.3683, −0.3488) | |

High density, 1:5 | 0.01673 | −0.05055 |

(0.01322, 0.02230) | (−0.05823, −0.04516) |

## References

- The prospects and promise of chimeric antigen receptor immunotherapy in multiple Myeloma.
*Br J Haematol.*2016; 173: 350-364 - CAR T cell therapy: A new era for cancer treatment (review).
*Oncol Rep.*2019; 42: 2183-2195 - The cellular immunotherapy revolution: Arming the immune system for precision therapy.
*Trends Immunol.*2019; 40: 292-309 - Molecular pathways: Breaking the epithelial cancer barrier for chimeric antigen receptor and T-cell receptor gene therapy.
*Clin Cancer Res.*2016; 22: 1559-1564 - CAR T cell immunotherapy for human cancer.
*Science.*2018; 359: 1361-1365 - CAR T cells in solid Tumors: Challenges and opportunities.
*Stem Cell Res Therapy.*2021; 12: 1-16 - Nonlinear dynamics of immunogen1c Tumors: Parameter estimation and global bifurcation analysis.
*Bull Math Biol.*1994; 56: 295-321 - Modeling immunotherapy of the Tumor-immune interaction.
*J Math Biol.*1998; 37: 235-252 - Behavior of tumors under nonstationary therapy.
*Physica D.*2003; 178: 242-253 - Metamodeling tumor-immune system interaction, tumor evasion and immunotherapy.
*Math Comput Modelling.*2008; 47: 614-637 - The dynamics of an optimally controlled tumor model: A case study.
*Math Comput Modelling.*2003; 37: 1221-1244 - A mathematical model for chronic myelogenous leukemia (CML) and T cell interaction.
*J Theoret Biol.*2004; 227: 513-523 - Mathematical modelling of spatio-temporal phenomena in tumour immunology.in: Tutorials in mathematical biosciences III. Springer, Berlin, 2006: 131-183
- A cellular automata model of Tumor–immune system interactions.
*J Theoret Biol.*2006; 239: 334-350 - PhysiCell: An open source physics-based cell simulator for 3-D multicellular systems.
*PLoS Comput Biol.*2018; 14e1005991 - Mathematical model of chimeric anti-gene receptor (CAR) T cell therapy with presence of cytokine.
*Numer Algebra Control Optim.*2018; 8: 63-80 - Pharmacology model of chimeric antigen receptor T-cell therapy.
*Clin Transl Sci.*2019; 12: 343-349 - A model-based investigation of Cytokine storm for T-cell therapy.
*IFAC-PapersOnLine.*2018; 51: 76-79 - The roles of T cell competition and stochastic extinction events in chimeric antigen receptor T cell therapy.
*Proc R Soc B.*2021; 28820210229 - CAR-T cell goes on a mathematical model.
*J Cell Immunol.*2020; 2 - CARTmath—A mathematical model of CAR-T immunotherapy in preclinical studies of hematological cancers.
*Cancers.*2021; 13: 2941 - Mathematical deconvolution of CAR T-cell proliferation and exhaustion from real-time killing assay data.
*J R Soc Interface.*2020; 17 - Stem-like Tumor-initiating cells isolated from IL13R$\alpha $2 expressing gliomas are targeted and killed by IL13-zetakine–redirected T cells.
*Clin Cancer Res.*2012; 18: 2199-2209 - Regression of glioblastoma after chimeric antigen receptor T-cell therapy.
*N Engl J Med.*2016; 375: 2561-2569 - CAR-T cells destroying glioblastoma cells.2017 (URL https://www.youtube.com/watch?v=8D6_xWYlLy4)
- Dynamic assessment of cell viability, proliferation and migration using real time cell analyzer system (RTCA).
*Cytotechnology.*2015; 67: 379-386 - CAR-T cells inflict sequential killing of multiple tumor target cells.
*Cancer Immunol Res.*2015; 3: 483-494 - C-jun overexpressing CAR-T cells are exhaustion-resistant and mediate enhanced antitumor activity.2019: 1-33 (BioRxiv)
- Study of dose-dependent combination immunotherapy using engineered T cells and IL-2 in cervical cancer.
*J Theoret Biol.*2020; 505 - A 3(2) pair of Runge - Kutta formulas.
*Appl Math Lett.*1989; 2: 321-325 - The MATLAB ODE suite.
*SIAM J Sci Comput.*1997; 18: 1-22 - Regulation of T cell expansion by antigen presentation dynamics.
*Proc Natl Acad Sci.*2019; 116: 5914-5919 - A new look at the statistical model identification.
*IEEE Trans Automat Control.*1974; 19: 716-723 - Classical mathematical models for description and prediction of experimental tumor growth.
*PLoS Comput Biol.*2014; 10e1003800 - A methodology for performing global uncertainty and sensitivity analysis in systems biology.
*J Theoret Biol.*2008; 254: 178-196 - CAR T cell trogocytosis and cooperative killing regulate Tumour antigen escape.
*Nature.*2019; 568: 112-116 - The role of immunological synapse in predicting the efficacy of chimeric antigen receptor (CAR) immunotherapy.
*Cell Commun Signaling.*2020; 18: 1-20 - Immunological synapse predicts effectiveness of chimeric antigen receptor cells.
*Mol Therapy.*2018; 26: 963-975 - Chimeric antigen receptor T Cells for sustained remissions in Leukemia.
*N Engl J Med.*2014; 371: 1507-1517 - CD19 CAR–T cells of defined CD4+:CD8+ composition in adult B cell all patients.
*J Clin Investig.*2016; 126: 2123-2138 - Decade-long leukaemia remissions with persistence of CD4+ CAR T cells.
*Nature.*2022; 602: 503-509 - CD4+ chimeric antigen receptor T cells in for the long journey.
*Immunol Cell Biol.*2022; 100: 304-307

## Article info

### Publication history

Published online: January 15, 2023

Accepted:
January 9,
2023

Received in revised form:
January 2,
2023

Received:
January 7,
2022

### Identification

### Copyright

© 2023 The Authors. Published by Elsevier B.V.

### User license

Creative Commons Attribution – NonCommercial – NoDerivs (CC BY-NC-ND 4.0) | How you can reuse

Elsevier's open access license policy

Creative Commons Attribution – NonCommercial – NoDerivs (CC BY-NC-ND 4.0)

## Permitted

### For non-commercial purposes:

- Read, print & download
- Redistribute or republish the final article
- Text & data mine
- Translate the article (private use only, not for distribution)
- Reuse portions or extracts from the article in other works

## Not Permitted

- Sell or re-use for commercial purposes
- Distribute translations or adaptations of the article

Elsevier's open access license policy