Modeling interaction of Glioma cells and CAR T-cells considering multiple CAR T-cells bindings

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) [], 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 [, ]. 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 []. 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.

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 []. We reparameterize the equation as

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

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.

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.

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.

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

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 [, ]. 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)$,

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.

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 []. 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) [, ] 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.

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 [] 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.

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,L=2$, and $M=2,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,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 [] 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.