Modeling CAR T-Cell Therapy with Patient Preconditioning

Abstract

The Federal Drug Administration approved the first Chimeric Antigen Receptor T-cell (CAR T-cell) therapies for the treatment of several blood cancers in 2017, and efforts are underway to broaden CAR T technology to address other cancer types. Standard treatment protocols incorporate a preconditioning regimen of lymphodepleting chemotherapy prior to CAR T-cell infusion. However, the connection between preconditioning regimens and patient outcomes is still not fully understood. Optimizing patient preconditioning plans and reducing the CAR T-cell dose necessary for achieving remission could make therapy safer. In this paper, we test treatment regimens consisting of sequential administration of chemotherapy and CAR T-cell therapy on a system of differential equations that models the tumor-immune interaction. We use numerical simulations of treatment plans from within the scope of current medical practice to assess the effect of preconditioning plans on the success of CAR T-cell therapy. Model results affirm clinical observations that preconditioning can be crucial for most patients, not just to reduce side effects, but to even achieve remission at all. We demonstrate that preconditioning plans using the same CAR T-cell dose and the same total concentration of chemotherapy can lead to different patient outcomes due to different delivery schedules. Results from sensitivity analysis of the model parameters suggest that making small improvements in the effectiveness of CAR T-cells in attacking cancer cells will significantly reduce the minimum dose required for successful treatment. Our modeling framework represents a starting point for evaluating the efficacy of patient preconditioning in the context of CAR T-cell therapy.

Appendices

Appendix A: Model Analysis

Here, we present a linear stability analysis for our model post-treatment. Because the treatment plans we consider are generally confined to within 5–20 days total for both chemotherapy and CAR T-cell injection, the dose schedules, $\nu _M(t)$ and $\nu _C(t)$, are zero for t greater than 5–20 days. Furthermore, the chemotherapy plans used for lymphodepletion are on the shorter end of this, spanning only 3–5 days. When $\nu _M(t) = 0$, the concentration of chemotherapy rapidly declines to zero and the only possible chemotherapy concentration at equilibrium is $M(t) = 0$. Thus, we can characterize the relevant equilibria for our model scenarios by analyzing Sys. (2a-c) with $\nu _M(t) = 0$ and $\nu _C(t) = 0$,

$$\begin{aligned} \dfrac{{\text{ d }}{T}}{{\text{ d }}{t}}&= aT(1 - T) - D_E - D_C, \end{aligned}$$

(6a)

$$\begin{aligned} \dfrac{{\text{ d }}{E}}{{\text{ d }}{t}}&= g - j_E\dfrac{D_E^2}{k + D_E^2}\ln \bigg (\frac{E+C}{K}\bigg )E - q_EET- m_EE, \end{aligned}$$

(6b)

$$\begin{aligned} \dfrac{{\text{ d }}{C}}{{\text{ d }}{t}}&= -j_C\dfrac{D_C^2}{k + D_C^2}\ln \bigg (\frac{E+C}{K}\bigg )C- q_CCT - m_CC, \end{aligned}$$

(6c)

where $D_E$ and $D_C$ are still the tumor cell lysis terms,

$$\begin{aligned} D_E = d_E\dfrac{(E/T)^l}{s + (E/T)^l}T ~~~~~\text {and}~~~~~ D_C = d_C\dfrac{(C/T)^l}{s + (C/T)^l}T. \end{aligned}$$

(7)

To begin our analysis, we first nondimensionalize Sys. (6) as follows. Let $x = bT$, $y = aE/g$, $z = aC/g$, $D_y = D_E/a$, $D_z = D_C/a$ and $t^{*} = at$. This results in twelve nondimensional parameters,

$$\begin{aligned} \begin{array}{c c c c c} d_y = d_E/a, &{}{} d_z = d_C/a, &{}{} j_y = j_E/a, &{}{} j_z = j_C/a, &{}{} k^*=b^2k/a^2, \\ [.2cm] K^{*} = aK /g, &{}{} m_y = m_E/a, &{}{} m_z = m_C/a, &{}{} q_y = q_E/(ab), &{}{} q_z = q_C/(ab), \\ [.15cm] \end{array} \end{aligned}$$

and $s^* = s(a/(gb))^l$. The parameter l is nondimensional in the original system, and so remains unchanged. Dropping stars for notational simplicity, the nondimensionalized system without chemotherapy or additional CAR T-cell doses is given by

$$\begin{aligned} \dfrac{{\text{ d }}{x}}{{\text{ d }}{t}}&= x(1 - x) - D_y - D_z, \end{aligned}$$

(8a)

$$\begin{aligned} \dfrac{{\text{ d }}{y}}{{\text{ d }}{t}}&= 1 -y \bigg (j_y\dfrac{D_y^2}{k + D_y^2}\ln \bigg (\frac{y+z}{K}\bigg ) + q_yx + m_y\bigg ), \end{aligned}$$

(8b)

$$\begin{aligned} \dfrac{{\text{ d }}{z}}{{\text{ d }}{t}}&= -z \bigg (j_z\dfrac{D_z^2}{k + D_z^2}\ln \bigg (\frac{y+z}{K}\bigg ) + q_zx + m_z\bigg ), \end{aligned}$$

(8c)

where $D_x$ and $D_y$ are now dimensionless ratio terms,

$$\begin{aligned} D_y = d_y\dfrac{(y/x)^l}{s + (y/x)^l}x ~~~~~\text {and}~~~~~ D_z = d_z\dfrac{(z/x)^l}{s + (z/x)^l}x. \end{aligned}$$

(9)

Next, we find the equilibria of the system, which occur at the intersections of the tumor, effector, and CAR T-cell nullclines, where $\dot{x} = 0$, $\dot{y} = 0$, and $\dot{z} = 0$, respectively. From Eq. (8a), we see that the tumor cell nullclines are $x=0$ and

$$\begin{aligned} x = 1- d_y\dfrac{(y/x)^l}{s + (y/x)^l}- d_z\dfrac{(z/x)^l}{s + (z/x)^l}. \end{aligned}$$

(10)

For the effector cell nullcline, we set the right-hand side of Eq. (6b) equal to zero and find that $\dot{y} = 0$ when

$$\begin{aligned} y = \left( m_y + q_yx + j_y\frac{D_y^2}{k+D_y^2}\ln \bigg (\frac{y+z}{K}\bigg )\right) ^{-1}. \end{aligned}$$

(11)

From Eq. (8c), we can see that CAR T-cell population change will be zero when $z = 0$, or when the following implicit equation is satisfied

$$\begin{aligned} 0 = m_z +q_zx + j_z\frac{D_z^2}{k+D_z^2}\ln \bigg (\frac{y+z}{K}\bigg ). \end{aligned}$$

(12)

We will refer to the intersection between the tumor nullcline $x=0$, the effector cell nullcline, Eq. (11), and the CAR T-cell nullcline $z = 0$ as the zero-tumor or tumor-free equilibrium. We refer to equilibria that occur at the intersection between the tumor nullcline Eq. (10), the effector cell nullcline, Eq. (11), and the CAR T-cell nullcline $z = 0$ as nonzero-tumor boundary equilibria. Equilibria that occur at the intersection between the tumor nullcline Eq. (10), the effector cell nullcline, Eq. (11), and CAR T-cell nullcline Eq. (12) will be referred to as interior equilibria.

In order to determine the stability of the steady states, we apply linear stability analysis. The linearized system at point (x, y, z) is summarized by the Jacobian, H, of Sys. (8) evaluated at (x, y, z). Let

$$\begin{aligned} J_y = -j_y\frac{D_y^2}{k + D_y^2}\ln \bigg (\frac{y+z}{K}\bigg )y ~~~\text {and}~~~ J_z = -j_z\frac{D_z^2}{k + D_z^2}\ln \bigg (\frac{y+z}{K}\bigg )z \end{aligned}$$

(13)

Then, the Jacobian of Eq. (8) evaluated at that point (x, y, z) is,

$$\begin{aligned} H|_{(x,y,z)} = \begin{pmatrix} 1- \left( 2x +\frac{\partial D_y}{\partial x} + \frac{\partial D_z}{\partial x}\right) &{} -\frac{\partial D_y}{\partial y} &{} -\frac{\partial D_z}{\partial z} \\ &{}&{} \\ \frac{\partial J_y}{\partial x}-q_yy &{} \frac{\partial J_y}{\partial y}-q_yx - m_y&{} \frac{\partial J_y}{\partial z} \\ &{}&{} \\ \frac{\partial J_z}{\partial x} - q_zz &{} \frac{\partial J_z}{\partial y}&{} \frac{\partial J_z}{\partial z} - q_zx - m_z \end{pmatrix}. \end{aligned}$$

(14)

Let

$$\begin{aligned} H|_{(x,y,z)} = \begin{pmatrix} h_{11}&{} h_{12}&{} h_{13} \\ h_{21}&{} h_{22} &{} h_{23} \\ h_{31} &{} h_{32} &{} h_{33}\\ \end{pmatrix}, \end{aligned}$$

so that substituting the partial derivatives into Eq. (14) gives

$$\begin{aligned} h_{11}&= 1 - \left( 2x + \frac{D_y}{x}\left( 1-\frac{ls}{s + (y/x)^l}\right) +\frac{D_z}{x}\left( 1 - \frac{ls}{s + (z/x)^l}\right) \right) , \end{aligned}$$

(15a)

$$\begin{aligned} h_{12}&= -D_y\bigg (\frac{ls}{s + (y/x)^l}\bigg ), \end{aligned}$$

(15b)

$$\begin{aligned} h_{13}&= -D_z\bigg (\frac{ls}{s + (z/x)^l}\bigg )\end{aligned}$$

(15c)

$$\begin{aligned} h_{21}&= \frac{1}{x}\bigg (\frac{2kJ_y}{k+ D_y^2}\bigg (\frac{ls}{s + (y/x)^l}-1\bigg ) - q_yxy\bigg ), \end{aligned}$$

(15d)

$$\begin{aligned} h_{22}&= \frac{J_y}{y}\left( 1 + \frac{2k}{k+ D_y^2}\frac{ls}{s + (y/x)^l}\right) -m_y -q_yx - \frac{y}{K(y+z)}\frac{j_yD_y^2}{k+D_y^2} \end{aligned}$$

(15e)

$$\begin{aligned} h_{23}&= \frac{-j_yy}{y+z}\frac{D_y^2}{k+D_y^2}\end{aligned}$$

(15f)

$$\begin{aligned} h_{31}&= \frac{1}{x}\bigg (\frac{2kJ_z}{k+ D_z^2}\bigg (\frac{ls}{s + (z/x)^l}-1\bigg ) - q_zxz\bigg ) \end{aligned}$$

(15g)

$$\begin{aligned} h_{32}&=\frac{-j_zz}{y+z}\frac{D_z^2}{k+D_z^2}\end{aligned}$$

(15h)

$$\begin{aligned} h_{33}&= \frac{J_z}{z}\left( 1 + \frac{2k}{k+ D_z^2}\frac{ls}{s + (z/x)^l}\right) -m_z -q_zx - \frac{z}{K(y+z)}\frac{j_zD_z^2}{k+D_z^2} \end{aligned}$$

(15i)

The simplest steady state to characterize is the tumor-free equilibrium. Evaluating the effector cell nullcline at $x=0$ and $z = 0$ yields $y = 1/m_{y}$. Hence, the tumor-free equilibrium occurs at $(x^*_0, y^*_0,z^{*}_0) = (0,1/m_y,0)$. Linearizing about this point, the Jacobian is

$$\begin{aligned} H|_{\left( 0,\frac{1}{m_y},0\right) } = \begin{pmatrix} 1-d_y -d_z &{} 0 &{} 0\\ \frac{-q}{m_y} &{} -m_y &{} 0\\ 0&{}0 &{} -m_z \end{pmatrix}, \end{aligned}$$

(16)

which has eigenvalues $\lambda _1= 1-d_z - d_y$, $\lambda _2 = -m_y$, and $\lambda _3 = -m_z$. The parameter values $d_y$, $d_z$, $m_y$, and $m_z$ will always be positive real numbers. It follows that $\lambda _1$, $\lambda _2$ and $\lambda _3$ are always real and negative when $d_y+d_z>1$. For the biologically relevant cases, we considered, the parameters $d_y > 1$ and $d_z> 1$, so this condition is satisfied. Hence, the equilibrium point will be a stable node, which aligns with the results reported by de Pillis et al. for their 2006 model. The stability of the tumor-free equilibrium reflects the idea that, once activated, the immune system can contain a tumor if it is small enough.

Now we turn our attention to the nonzero-tumor boundary equilibria. For these equilibria, $z=0$ still, but now

$$\begin{aligned} x = 1-D_y - D_z = 1 - d_y\frac{(y/x)^l}{s + (y/x)^l}. \end{aligned}$$

(17)

In this case, we can solve Eq. 17 for y as a function of x along the tumor nullcline,

$$\begin{aligned} L_1(x) = \bigg (\frac{s(1-x)x^l}{d_y - (1-x)}\bigg )^{1/l}. \end{aligned}$$

(18)

Substituting Eq.( 17) and $z = 0$ into the effector cell nullcline yields an implicit curve in the plane $z = 0$,

$$\begin{aligned} L_2(x,y): ~~y = \left( m_y + q_yx+j_y\frac{(1-x)^2x^2}{k+(1-x)^2x^2}\ln \left( \frac{y}{K}\right) \right) ^{-1}. \end{aligned}$$

(19)

The nonzero tumor boundary equilibrium points are the intersections of Eqs. (18) and (19) for each set of parameters. See Fig. 1a for a plot of $L_1(x)$ (tumor cell nullcline) and $L_2(x,y)$ (effector cell nullcline) for the Patient 3 parameters with equilibria marked. For reasonable parameter ranges, reported in Table 1, we observed that the model has two positive nonzero tumor equilibrium points in the plane $z = 0$ because the parameter s is large enough that $L_1(x)$ exceeds $L_2(x,y)$ for most of the interval $x \in (0,1)$, but $L_1(0)=0 < L_2(0,1/m)$ and $L_1(1)=0 < L_2(1,y^*)$.

The first nonzero tumor boundary equilibrium occurs at a large number of tumor cells, $x_1^* \approx 1$, a small number of effector cells, $y_1^* \approx L_1(1) = 1/(q+m)$, and $z_1^{*} = 0$. Due to the location near the tumor cell carrying capacity, we call this the high-tumor equilibrium. The second nonzero tumor boundary equilibrium occurs at a small number of tumor cells, a moderate number of effector cells and $z = 0$. Although it is not possible to find exact, closed-form expressions for the nonzero tumor boundary equilibria because it requires finding the roots of a quintic polynomial, we can use the Routh–Hurwitz criterion to find conditions under which the high-tumor equilibrium is stable  (Wiggins 2003, p. 14). Let $\tau (M)$ denote the trace of a matrix M and $\Delta (M)$ denote the determinant of a matrix M. For a three-dimensional system, the Routh–Hurwitz criterion says that all eigenvalues of matrix H will have negative real part if the following criteria are satisfied,

• 1. $\tau (H) < 0$,

• 2. $\Delta (H) < 0$,

• 3. $\tau (H^2)-\tau (H)^2 < 2\Delta (H)/\tau (H).$

First, we use the nullclines to simplify the Jacobian. Substituting $z = 0$ and $D_y = (1-x)x$ from Eq. (17) into Sys. (16), the Jacobian has entries

$$\begin{aligned} h_{11}&= -x + \frac{ls(1-x)}{s + (y/x)^l}, \end{aligned}$$

(20a)

$$\begin{aligned} h_{12}&= -x\bigg (\frac{ls(1-x)}{s + (y/x)^l}\bigg ), \end{aligned}$$

(20b)

$$\begin{aligned} h_{13}&=0\end{aligned}$$

(20c)

$$\begin{aligned} h_{21}&= \frac{1}{x}\bigg (\frac{2kJ_y}{k+ (1-x)^2x^2}\bigg (\frac{ls}{s + (y/x)^l}-1\bigg ) - q_yxy\bigg ), \end{aligned}$$

(20d)

$$\begin{aligned} h_{22}&= \frac{J_y}{y}\left[ \left( \frac{2k}{k+ (1-x)^2x^2}\right) \left( \frac{ls}{s + (y/x)^l}\right) + 1\right] \nonumber \\&~~~~ -q_yx - m_y - \frac{j_y}{K}\left( \frac{(1-x)^2x^2}{k+ (1-x)^2x^2}\right) \end{aligned}$$

(20e)

$$\begin{aligned} h_{23}&= -j_y\frac{(1-x)^2x^2}{k+(1-x)^2x^2}\end{aligned}$$

(20f)

$$\begin{aligned} h_{31}&=0 \end{aligned}$$

(20g)

$$\begin{aligned} h_{32}&= 0\end{aligned}$$

(20h)

$$\begin{aligned} h_{33}&= -m_z -q_zx \end{aligned}$$

(20i)

Before applying the Routh–Hurwitz criterion, we will need to compute $\tau (H)$, $\Delta (H)$, and $\tau (H^2)$. The trace of the Jacobian is

$$\begin{aligned} \tau (H) = h_{11} + h_{22} + h_{33}. \end{aligned}$$

(21)

This will be negative, satisfying criterion 1, when the three diagonal entries are negative. Next, consider the determinant of the Jacobian

$$\begin{aligned} \Delta (H)&=h_{11}h_{22}h_{33} - h_{12}h_{21}h_{33}, \end{aligned}$$

(22)

due to zero entries in the Jacobian. Furthermore, with $x \approx 1$, $h_{12}$ will be very nearly 0 and so we can consider

$$\begin{aligned} \Delta (H)&= h_{11}h_{22}h_{33}. \end{aligned}$$

(23)

Hence, the determinant will also be negative, satisfying criterion 2, when the three diagonal entries share the same sign. Squaring H and summing the diagonal entries yields

$$\begin{aligned} \tau (H^2)&= h_{11}^2 + h_{22}^2 + h_{33}^2 + 2h_{12}h_{21}\end{aligned}$$

(24)

$$\begin{aligned}&= h_{11}^2 + h_{22}^2 + h_{33}^2, \end{aligned}$$

(25)

again because $h_{12}$ goes to zero for $x \approx 1$. We know that $\Delta (H)/\tau (H) > 0$ when both are negative, so showing $\tau (H^2) < \tau (H)^2$ will be sufficient to satisfy the third criterion if the first two are met. Squaring $\tau (H)$ yields

$$\begin{aligned} \tau (H)^2&= h_{11}^2 + h_{22}^2 + h_{33}^2 + 2h_{11}h_{22} +2 h_{22}h_{33} + 2h_{11}h_{33}. \end{aligned}$$

(26)

All the terms in the sum will be positive when $h_{11} < 0$, $h_{22} < 0$, and $h_{33} < 0$ in which case it is clear that $\tau (H^2) < \tau (H)^2$. Thus, all three of the Routh–Hurwitz criterion will hold when the three diagonal entries of H are negative.

We start with the first diagonal entry. From Eq. (20a), we can see that when $x \approx 1$ it follows that that $h_{11} \approx -1$.

Next, we consider the second diagonal entry, $h_{22}$, defined in Eq. (20e). Substituting in $J_y/y = 1/y - m_y - q_yx$ yields,

$$\begin{aligned} \begin{aligned} h_{22}&= \frac{1}{y}\left[ \left( \frac{2k}{k+ (1-x)^2x^2}\right) \left( \frac{ls}{s + (y/x)^l}\right) + 1\right] \\&~~~~ - (q_yx +m_y)\left[ \left( \frac{2k}{k+ (1-x)^2x^2}\right) \left( \frac{ls}{s + (y/x)^l}\right) + 2\right] \\&~~~~~~~~~~~~~ - \frac{j_y}{K}\left( \frac{(1-x)^2x^2}{k+ (1-x)^2x^2}\right) \end{aligned} \end{aligned}$$

(27)

Then, noting that the final term will be negative, but close to zero and that the first fractional term will approach 2 when $x \approx 1$, $h_{22}$ will be less than 0 when

$$\begin{aligned} \frac{1}{y} < (q_y + m_y)\frac{2\left( \frac{ls}{s + (y/x)^l}\right) + 2}{2 \left( \frac{ls}{s + (y/x)^l}\right) + 1}. \end{aligned}$$

(28)

Hence, when $x \approx 1$, the entry $h_{22}$ will be negative when

$$\begin{aligned} y > \frac{1}{q_y + m_y}. \end{aligned}$$

(29)

Finally, from Eq. (20i), we can see that for any positive value of x, including $x \approx 1$, $h_{33} < 0$ because both parameter values, $m_z$ and $q_z$, are real and positive.

Thus, when $z = 0$, $x \approx 1$, and $y >1/(q_y+m_y)$, the diagonal entries of H, $h_{11}$, $h_{22}$ and $h_{33}$, will all be negative. Hence, the Routh–Hurwitz Criterion will be satisfied and all eigenvalues of the system will have negative real parts. These conditions are satisfied by the high-tumor equilibrium for the patient parameters considered here, so it follows by the Routh–Hurwitz criterion that the high-tumor steady state is stable.

Characterizing the interior equilibria analytically is intractable. However, we located equilibria at the intersection of nullsurfaces, and characterized their stability numerically. For the parameter sets considered here, there are up to four interior equilibria with $x \ge 0$, $y \ge 0$, and $z \ge 0$. To explore these equilibria, we constructed a bifurcation diagram by varying the maximum CAR T-cell recruitment rate, $j_C$, from 0 to 3 while holding all other parameters at their Patient 3 value (Fig. 1c). With strictly positive parameter values, an unstable node and a saddle point are present at a low (but nonzero) population level for all three cell types. The stable manifold of this saddle point separates the basin of attraction for healthy outcomes from the basin of attraction for unhealthy outcomes. Increasing $j_C$ causes this saddle point to move towards larger tumor and CAR T-cell values as the nonzero CAR T-cell nullcline shifts upward with respect to tumor and CAR T-cells. For the patient 3 values, a saddle node bifurcation occurs when $j_C = 0.5105$ and a stable node appears simultaneously with another saddle point near the tumor cell carrying capacity. The stable manifold of this saddle point separates the basins of attraction for the high-tumor equilibrium from the basin of attraction for this new, stable coexistence equilibrium. Figure 1b shows the null surfaces with the unperturbed patient 3 parameter set, i.e., with $j_C = 0.6$, so all four interior equilibria are visible. Continuing to increase $j_C$ past the point $j_C = 2.56$, the nonzero CAR T-cell nullcline shifts far enough upward with respect to tumor and CAR T-cells that the stable coexistence equilibrium and the saddle node at high CAR T-cell values collide and disappear in another saddle-node bifurcation. This leaves two remaining unstable interior equilibria, and the stable manifold of the remaining saddle node now separates the basins of attraction for the low tumor and high-tumor equilibria. Varying other parameters associated with CAR T-cells can also shift the nonzero CAR T-cell nullsurface through this sequence of bifurcations.

The positive quadrant is an invariant set of system 6. In order for a trajectory to leave the first quadrant, it would have to cross an axis, which requires $\dot{x} < 0$ somewhere along the $x=0$ axis, $\dot{y} < 0$ somewhere along the $y=0$ axis or $\dot{z} < 0$ somewhere along the $z = 0$ axis. However, when $x=0$ it is clear from Eq. (6a) that $\dot{x} =0$. From Eq. (6b), we can also see that when $y=0$, it follows that $\dot{y} = 1$, which is positive. And when $z = 0$, $\dot{z} =0$ also.

Appendix B: Parameter Selection

We selected parameters based on previous models in the literature. The metastatic melanoma parameter sets came directly from the parameter values for Patients 9 and 10 in the governing equations for tumor cells and CD8$^+$ T-cells reported by de Pillis et al. (2005). However to counteract the effect of eliminating the other immune cell types (NK cells and circulating lymphocytes) that stimulate CD8$^+$T cells, the base recruitment rate of effector cells from Kuznetsov et al. (1994) was taken to be the base recruitment rate of effector cells for Patients 3 and 4.

We set the parameter values for our diffuse large B-cell lymphoma patient, Patient 1, to the midpoint of the ranges reported by Rösch et al. (2016) when there was an analogous term. We could use most parameters directly from Rösch et al. (2016). However, the parameters involved in the novel tumor-cell lysis term introduced by de Pillis et al. did not correspond directly to parameters in Rösch’s model. We expect endogenous T-cells in lymphoma patients to be less effective at killing tumor cells than CAR T-cells, so we set $d_E = 2.02$, or roughly 90% as effective as CAR T-cells. We set the remaining two parameters, l and s, to the average of the two patients considered in the model from de Pillis et al. (2005). The Rösch et al. model assumed exponential tumor growth, rather than logistic so we set $b = 5 \times 10^{-13}$ because this choice allows tumor cells to exhibit essentially exponential growth towards the number of lymphocytes in the body in the absence of an immune response.

The chronic lymphocytic leukemia patient parameter values were set to the midpoint of the ranges reported by Nanda et al. (2013) when there was an analogous term, adjusting for differences in units. The Nanda et al. model also assumed exponential tumor growth, rather than logistic so we set $b = 5 \times 10^{-13}$. Under the same reasoning as the lymphoma parameters, we set $d_E = 2.02$, and set the final parameters involved in the novel interaction term (l and s) to the average of the two patients considered in the model from de Pillis et al. (2005).

CAR T-cell parameters that differ from their corresponding parameter in the endogenous effector cells $(d_C, m_C, j_C,$ and $q_C)$ were determined from the values reported by Kimmel et al. (2019). For all patients, the CAR T-cell death rate and the maximum CAR T-cell kill rate were set to the median value for their deterministic model, converting units if necessary. Note that the maximum tumor kill rate is lower for CAR T-cells than for endogenous effector cells in the metastatic melanoma patients. This may seem counter-intuitive; however, the parameter values from de Pillis et al. (2005) were fit to clinical trials of successful tumor infiltrating lymphocyte (TIL) treatment. TILs are currently more effective than CAR T-cells against metastatic melanoma (Itzhaki et al. 2020). The CAR T-cell exhaustion rate and the effector cell carrying capacity for blood cancers (Patients 1 and 2), came from converting the corresponding values in the model from Kimmel et. al (2019) to the our units. For the metastatic melanoma parameter sets (Patient 3 and 4), we use a higher inactivation rate so that the order of magnitude of the ratio of the tumor carrying capacity to deactivation is preserved. Numerous studies have shown the microenvironment around solid tumors to be more hostile to immune cells than in the case of hematologic cancers (Villanueva and Herlyn 2008). For metastatic melanoma patients we also was scaled down the carrying capacity of effector cells by $10\%$ to account for the lower access of T-cells to solid tumor cells (Villanueva and Herlyn 2008; Jorge et al. 2020). In their model, Kimmel et al. (2019) use a functional form which places our maximum CAR T-cell recruitment rate, $j_C$, in the interval [0.11, 1.31]. We chose to set the CAR T-cell recruitment rate to be a factor larger than the endogenous cell recruitment rate: $j_C = 22j_E$. The choice of 22 as a scale factor places the resulting values within the range reported by Kimmel et al. and produces simulation results spanning an interesting range of biological outcomes.

The kill parameters relating to chemotherapy were taken from  de Pillis et al. (2006), since the patients in the metastatic melanoma study were also treated with fludarabine. The fractional cell kills used by de Pillis et al. were $K_T = 0.9$ and $K_E=0.6$ which assumes that chemotherapy is one log-kill and that it kills a larger fraction of tumor cells than host cells. However, to reflect the fact that CAR T-cell treatment is commonly used for patients whose disease is less responsive to chemotherapy, we reduced $K_T$ to 0.7. We set $K_C = K_E = 0.6$, therein assuming that CAR T-cells will be impacted by chemotherapy at the same rate as endogenous T-cells. The decay rate $\gamma$ can be calculated from the half-life of a substance, $\tau$, by $\gamma = \ln {(2)}/ \tau$. The choice of strength for chemotherapy regimens was chosen by working backwards to see what values resulted in reasonable concentrations of fludarabine building up in the system, where “reasonable” concentrations matched those reported by  Ju et al. in a pharmacokinetic study of the drug (2014). For 25 mg of fludarabine administered over one half hour each day for 5 days, Ju et al. reported a peak concentration of $C_{\mathrm{max}} = 1222 (668-1732) \text {ng/mL} = 3.34~ \mu$M. We found that a dose strength of $S = 77~ \mu$Mday$^{-1}$ achieved this $C_{\mathrm{max}}$ in our simulations, so this was used as medium strength. The high strength was set to be $S=125~ \mu$Mday$^{-1}$ in order to achieve the same area under the concentration curve with only 3 days of injections. If we chose to instead match the dose administered (area under $\nu _M(t))$ between the two plans, the patient outcomes were qualitatively similar.

Note that while the treatment parameters are well-founded in pharmacokinetic data, the patient parameters have been drawn from previous theoretical studies for different cancer types. Some have been fit to patient data, but not all. In order to strengthen the impact of any suggestions made by this model, it remains to calibrate parameters model to data from previous CAR T cell therapy patients and validate the model by evaluating its predictive performance on untrained data.

Appendix C: Numerical Simulations and Sensitivity Analysis

We implemented numerical simulations with parameters for Patient 1–4 listed in Table 1 using a combination of fourth order Runge–Kutta integration schemes in MATLAB. During chemotherapy, we enforced a fixed step size of ${\text {d}}\!{t} = 0.0208$ corresponding to approximately one half hour or 1/48 of a day using ode4(). The remainder of each simulation was completed with ode45(), which dynamically chooses the time-step.

To discover which components of the model contribute most significantly to the effectiveness of CAR T-cell injections, we performed a numerical sensitivity analysis assessing how the model parameters impact the shape of the separatrix between the basin of attraction for the high-tumor equilibrium and the basin of attraction for the tumor-free equilibrium.

As a baseline for Patients 1–4, we calculated the precise threshold of CAR T-cells above which a patient initial condition will move towards the healthy, tumor-free equilibrium when the system starts at the low-TB initial condition. We call this number of effector cells the CAR T-cell success threshold, because any treatment plan which injects a dose of CAR T-cells exceeding this quantity will be successful (in theory) against tumor burdens up to and including the low TB, (Fig. 7a). For the reasonable patient initial conditions considered in our simulations, post-chemotherapy tumor burdens generally fell below the low TB value at the time of CAR T-cell injection. After calculating the threshold with unperturbed parameters, each parameter was increased by $1\%$ while holding all other parameter values constant. The relative change in the CAR T-cell success threshold was computed by subtracting the unperturbed CAR T-cell success threshold from the perturbed CAR T-cell success threshold and dividing the difference by the unperturbed threshold. The same procedure was followed for a $1\%$ decrease in each parameter. The results of sensitivity analysis for Patient 1 are shown in Fig. 7b and the results for Patients 2–4 are shown in Fig. 8.

For all patients, varying the tumor growth rate, a, had a significant impact on the CAR T-cell success threshold, indicating that cancers characterized by a higher growth rate will require higher CAR T-cell doses if administered at the same patient initial condition. In contrast with Patients 1–2, the tumor carrying capacity, 1/b, also had a non-negligible effect on the success threshold for Patients 3–4. This seems to be a consequence of the fact that the low TB condition is only one order of magnitude below the tumor carrying capacity for Patient’s 3 and 4, so the behavior of the tumor population can enter the saturating regime of logistic growth, shaped by the carrying capacity.

Fig. 8

Results of sensitivity analysis for a Patient 2, b Patient 3, and c Patient 4. The dose of CAR T-cells required to achieve a healthy patient outcome from the low tumor initial condition is most sensitive to the tumor growth rate, a, and the exponent of the tumor cell lysis rate, l (a parameter associated with the cooperativity of CAR T-cells). Other parameters dictating the behavior of CAR T-cells also have a non-negligible impact on the CAR T-cell success threshold. The inverse of the tumor cell carrying capacity, b, has little effect on the threshold for Patients 1–2, but a noticeable effect for Patients 3–4. It is worth noting that the low tumor initial condition is relatively close to the tumors carrying capacity for Patients 3–4 compared to Patients 1–2

Fig. 9

Impact of varying parameters characterizing the lysis of tumor cells by CAR T-cells on patient trajectories initiated at $(T_0 = 5\times 10^8, E_0 = 3\times 10^4, C_0 = 10^6)$ for Patient 1 and $(T_0 = 10^8, E_0 = 3\times 10^4, C_0 = 10^6)$ for Patient 3. The maximum tumor cell lysis rate, $d_C$, was varied between 0 and 10 for a Patient 1 and b Patient 3. At low values of $d_C$, trajectories approach an unhealthy patient outcome. As $d_C$ increases, CAR T-cells are able to mount a larger response against the tumor, and eventually trajectories transition to approaching the healthy outcome. Continuing to increase $d_C$, CAR T-cells eradicate the tumor more quickly and with a smaller peak. The exponent in the tumor cell lysis rate, l, was varied from 0 to 3 for c Patient 1 and d Patient 3. For low values of l, CAR T-cells effectively kill tumor cells independently so trajectories rapidly reach the zero-tumor burden with low spikes in CAR T-cells. As l increases, higher cooperation between CAR T-cells is required and so larger peaks in CAR T-cells are seen. Eventually, the level of cooperation required is too high, and trajectories now approach an unhealthy outcome, though CAR T-cells still increase initially. Once $l > 2$, the CAR T-cells no longer spike after injection, as trajectories quickly approach the high-tumor state. The half-saturation value in the tumor cell lysis rate, s, was varied from 0 to 1 for e Patient 1 and f Patient 3. Increasing s increases the effector-target ratio required to reach the saturation regime of tumor cell lysis. For low values of s the tumor is eradicated quickly with few CAR T-cells. As s increases, more CAR T-cells are required to eliminate the tumor observed in larger CAR T-cell peaks, until eventually this dose of CAR T-cells can no longer mount an effective response to this tumor burden. Further increasing s decreases the peak of the failed response

The parameters characterizing the activity of CAR T-cells also have a nonnegligible impact on the success threshold. For all patients, changing the exponent in the tumor cell lysis rate, l, has the largest impact on the CAR T success threshold, with a $1\%$ change in l leading to a 5–7$\%$ change in the threshold. The parameter l encodes how the lysis rate depends on the ratio of CAR T-cells/tumor cells. Thinking of this ratio as a chemical concentration we can think of varying l as allowing the expression D to range between non-cooperative (completely independent) binding between CAR T-cells and surface antigens on the tumor cells when $l=1$, and cooperative binding when $l>1$. The fact that a $1\%$ change in l can result in a $7\%$ shift in the threshold indicates that small changes in CAR T-cell cooperativity will affect clinical outcomes. For Patients 1–2, and 4, the next most inlfuential parameter is the maximum recruitment rate of CAR T-cells, $j_C$. Larger values of $j_C$ decrease the necessary injected dose of CAR T-cells because those cells that are injected are able to expand further. For Patient 3, perturbations to the maximum tumor kill rate by CAR T-cells, $d_C$, have the most outsized impact after l. This parameter defines the upper limit of the ratio-dependent cell-lysis term, suggesting that small improvements to the effectiveness of CAR T-cell activity could significantly reduce the dose required by patients.

Sensitivity analysis based on the high tumor burden indicated that the parameters associated with tumor cell lysis by CAR T-cells, $d_C$, l, and s, play an important role in the success of CAR T-cell treatment so we performed further numerical investigations to show how varying these parameters impacts patient trajectories. Figure 9 shows patient trajectories resulting from running simulations with the Patient 1 and Patient 3 parameter sets from an initial condition near the separatrix between the basins of attraction for the different equilibria. Either $d_C$, l, or s was varied through a relevant range while holding all other parameter values constant at the value listed in Table 1 and the trajectories were color coded according to the value of the parameter in question.

First the maximum tumor cell lysis rate, $d_C$, was varied between 0 and 10. When tumor cell lysis is zero, the CAR T-cell recruitment is also affected, and CAR T-cells cannot proliferate so the trajectory proceeds in the $C = 0$ plane to the high-tumor equilibrium. As $d_C$ increases a larger spike in CAR T-cells is seen. Though trajectories still reach the high-tumor equilibrium for Patient 1, they now approach the unhealthy coexistence equilibrium for Patient 3. Eventually the maximum kill rate is high enough that a large spike in CAR T-cells is sufficient to move the patient to the healthy zero-tumor outcome. Continuing to increase $d_C$ the system approaches the healthy zero-tumor equilibrium faster and with a shorter and shorter spike in CAR T-cells.

Next, the exponent in the tumor cell lysis rate, l, was varied from 0 to 3. When $l=0$, with even one CAR T-cell present tumor cells are killed at a rate independent of further increases in CAR T-cells. Consequently, the tumor is eradicated very quickly with very little expansion of CAR T-cells beyond the initial dose. As l increases more cooperation is required between CAR T-cells to attack tumor cells, so we see larger and larger spikes in CAR T-cells needed to contain the tumor and reach a healthy outcome. Eventually, the level of cooperativity required is too great for this CAR T-cell dose to handle this tumor burden. At this point, even though CAR T-cells expand post injection, trajectories approach an unhealthy outcome. Once l is greater than 2, the injection of CAR T-cells has virtually no effect on the disease, and CAR T-cells and tumor cells monotonically approach the high-tumor equilibrium.

The half-saturation value in the tumor cell lysis rate, s, was varied from 0 to 1. Increasing s increases the effector target ratio required to reach the saturation regime of tumor cell lysis. When $s=0$, the lysis of tumor cells by CAR T-cells transitions to the saturation regime immediately and so the tumor is eradicated rapidly by a small number of CAR T-cells. As s increases, the CAR T-cell to tumor cell ratio needed to reach the saturation kill rate increases, and the eradication time and fold-increase in CAR T-cells increases along with it. At a critical value around $s = 0.5$ for Patient 1 and $s = 0.25$ for Patient 3, trajectories now approach an unhealthy outcome and the peak CAR T-cell population decreases as s increases. The effector to target ratio required to reach the saturation regime is now too high for this dose to handle this initial condition.