CARTmath-A mathematical model of CAR-T immunotherapy in preclinical models

Immunotherapy has gained great momentum with chimeric antigen receptor T cell 1 (CAR-T) therapy, in which patient’s T lymphocytes are genetically manipulated to recognize 2 tumor-specific antigens increasing tumor elimination efficiency. In the last years, CAR-T cell 3 immunotherapy for hematological malignancies achieved a great response rate on patients and is a 4 very promising therapy for several other malignancies. Each new CAR design requires a preclinical 5 proof-of-concept experiment using immunodeficient mouse models. The absence of a functional 6 immune system in these mice makes them simple and suitable to be mathematically modeled. 7 In this work, we developed a three population mathematical model to describe tumor response 8 to CAR-T cell immunotherapy in immunodeficient mouse models, encompassing interactions 9 between a non-solid tumor and CAR-T cells (effector and long-term memory). We account for 10 several phenomena, such as tumor-induced immunosuppression, memory pool formation, and 11 conversion of memory into effector CAR-T cells in the presence of new tumor cells. Individual 12 donor and tumor specificities were considered as uncertainties in the model parameters. Our 13 model is able to reproduce several CAR-T cell immunotherapy scenarios, with different CAR 14 receptors and tumor targets reported in the literature. We found that therapy effectiveness mostly 15 depends on some specific parameters such as the differentiation of effector to memory CAR-T 16 cells, CAR-T cytotoxic capacity, tumor growth rate, and tumor-induced immunosuppression. In 17 summary, our model can contribute to reduce and optimize the number of in vivo experiments 18 with in silico tests to select specific scenarios that could be tested in experimental research. Such 19 in silico laboratory was made available in a Shiny R-based platform called CARTmath. It is an 20 open-source, easy to run simulator, available at github.com/tmglncc/CARTmath or directly on 21 the webpage cartmath.lncc.br, containing this manuscript results as examples and documentation. 22 The developed model, together with the CARTmath platform, provides potential use for assessing 23 different CAR-T cell immunotherapy protocols and associated efficacy, becoming an accessory 24 towards in silico trials. 25


Introduction
Adoptive cell therapies have been considered a major advance in the fight against 29 several cancers, especially those associated with the hematopoietic system [1]. Chimeric 30 antigen receptor T (CAR-T) cell immunotherapy is an adoptive cellular therapy in which 31 T lymphocytes are taken from the patient's blood, genetically modified to recognize 32 in vitro/in vivo experiments in substitution to laboratory xenograft mouse studies. This 48 requires, at a first step, the development of mathematical models to accurately describe 49 experimental data already present in the literature. 50 Recent in vitro/in vivo experimental studies investigated the relationship between 51 immunotherapy with CAR-T cells and the development of immunological memory 52 [7][8][9][10]. Using an immunodeficient mouse model, [7] showed that CAR-T 123 therapy can 53 eliminate HL and provide long-term immunity against a challenge of the same tumor. 54 Immune checkpoint blockade (ICB) associated with CAR-T cell therapy is also under 55 investigation in mouse models where CAR-T cell therapy fails. With the use of immunod- 56 eficient mouse models, [11] showed that tumor expressing indoleamine 2,3-dioxygenase 57 (IDO) activity, an intracellular enzyme that has an inhibitory effect on T cells, can be 58 better controlled by combining the CAR-T cell therapy with 1-methyl-tryptophan (1-MT), 59 an IDO inhibitor. By the end of 2016, four different ICB drugs were also approved for 60 the treatment of lymphoma, melanoma, among other cancers. Although the success 61 of CAR-T cell therapy against hematologic cancers is promising, the mechanisms as- 62 sociated with failures have been reported and are the subject of recent investigations 63 [12]. Notably, many challenges remain to be addressed to improve response rates such 64 as minimum effective CAR-T cell dose, selection of CAR-T subtypes, adverse effects 65 management, combination of therapies, formation and maintenance of immunological 66 memory, suppressive microenvironment, and patient specificity, to mention a few [13]. In 67 this context, mathematical models may contribute to understanding the factors involved 68 in malignant transformation, invasion, and metastasis, as well as to examine responses 69 to therapies [8,[14][15][16][17][18][19], confronting hypotheses and testing different settings [20][21][22]. 70 Simplified mathematical models can be used to investigate some of those issues 71 and have several advantages, such as reduced simulation time, which allows testing 72 several experiments in a relatively short period, and the gain of interpretability, that 73 is, understanding all terms of the model and their impact on the results. Several math-74 ematical models in the literature use predator-prey dynamics to explore CAR-T cells' 75 kinetics [23,24]. However, most of these models do not consider the complex dynamics 76 of effector CAR-T cell differentiation into memory CAR-T cells and then back to effector 77 cells after antigen recognition. In this work, we focus on the development of a simple 78 mathematical model using a system of three ordinary differential equations (ODEs) 79 to describe CAR-T and tumor cell populations dynamics in immunodeficient mouse VP heterogeneity allowed exploring the factors that impact therapy outcomes. We also 87 used the model to retrieve a different CAR-T cell immunotherapy scenario, using data 88 from [11] for CAR-T 19 immunotherapy on RAJI tumors. We remark that changing    Table 1. Figure 1. Schematic description of the model structure. Effector CAR-T cells proliferate, have a cytotoxic effect on tumor cells, differentiate into memory CAR-T cells, and die naturally or are impaired by tumor-induced immunosuppressive mechanisms. The long-term memory CAR-T cells also die naturally and are readily responsive to the tumor-associated antigen and when they interact with tumor cells, they differentiate back into effector CAR-T cells, producing a rapid immune response against the tumor. Tumor cells grow subject to available resources in the microenvironment and are killed by effector CAR-T cells.

114
The dynamics of the effector (activated) CAR-T cells is described by: The first right-hand side term of (1) specifies that effector CAR-T cells undergo expansion 115 due to proliferation at a rate of φ. This population is reduced at a rate ρ, which includes 116 the natural death of the effector CAR-T cells and also its differentiation into long-term 117 memory CAR-T cells [25,26] according to the linear progression model described in 118 [8,27]. The term θTC M describes the activation of memory CAR-T cells into the effector 119 state, due to the contact with tumor cells. Indeed, it is well-known that memory CAR-T 120 cells may provide long-lasting protection to the specific tumor/antigen [28,29]. At any 121 future time in which memory CAR-T cells come into contact with the same tumor cells, 122 they can rapidly be converted into effector CAR-T cells, readily activated to prevent 123 tumor progression when enough memory cells are present in the system. It is also known 124 that memory CAR-T cells have a lower activation threshold, which eases the secondary 125 response to a future tumor recurrence [30]. Finally, the term αTC T models the combined immune checkpoint [31]. Many other immune checkpoint molecules have already been 131 described such as IDO, LAG3, and VISTA with high potential to be used as target therapy 132 [32,33]. IDO is an intracellular enzyme that has an inhibitory activity on T cells and is 133 overexpressed in several human cancers [34,35]. In this work, we consider that inhibitory 134 signals prevail, resulting in positive values for α and we will specifically consider the 135 effects of IDO inhibition later on. Moreover, we assume that a given dose of effector 136 CAR-T cells is introduced into the system as an adoptive therapy.

137
The dynamics of the immunological memory CAR-T cells, a key of the adaptive 138 immune system [8,36], is modeled by: Equation (2) assumes that memory CAR-T cells are formed exclusively from the differ-  The response of tumor cells to the CAR-T immunotherapy is modeled by: In the absence of immunosurveillance, we assume a density-dependent growth of 147 cancer cells due to the limitation of available resources in the tumor microenvironment, 148 characterizing the existence of intraspecific tumor cell competition. Tumor growth is 149 described using a logistic growth with maximum growth rate r and carrying capacity 150 1/b [37,38]. Finally, we assume that effector CAR-T cells kill tumor cells upon contact, at 151 a constant per capita rate γ; this anti-tumor cytotoxicity mechanism is modeled by the 152 term γC T T [39][40][41].

153
All parameters assume positive values. Further, based on reasonable biological 154 assumptions, we impose two additional conditions on the model parameters as follows.

155
First, we note that parameter ρ may be written as ρ = η + , where η is the natural 156 mortality rate of effector CAR-T cells and is the rate of memory formation. We will detected on peripheral blood analyses after tumor elimination [7,11], memory CAR-T 164 cells can survive for years [42], providing long-term protection against the target antigen 165 presented by the tumor. This biological behavior is obtained by imposing the restriction 166 φ < ρ, which ensures that effector CAR-T cells decay to zero in absence of tumor cells.
167 Table 1 summarizes the two restrictions imposed on the values of the parameters.
168 Table 1: Summary of the model parameters and the two restrictions imposed among φ, ρ, and . .

Parameter Meaning Unit
φ C T proliferation rate day −1 ρ C T reduction rate, encompassing the C T natural death and their differentiation into Inverse of the tumor carrying capacity with the total number of cells as shown in [43]. The cytotoxic activity of CAR-T on tumor 186 cells was retrieved from a standard in vitro 4-hour chromium-51 release assay [44]. For 187 the RAJI tumor, inference of the CAR-T cell inhibition due to interaction with tumor 188 (α) was performed based on data from [11]. All used data were extracted using the

Mathematical analysis of model dynamics 194
We perform a mathematical analysis of the model long-term dynamics, finding the 195 steady states and characterizing their stability. In order to simplify the calculations, we where X, Y, Z, and τ are dimensionless variables. Note that Z represents a fraction of 198 the T tumor cell population with respect to the carrying capacity. The dimensionless 199 system is given by where p = ρ−φ r , q = θ br , s = α br , u = r , and w = µ r ; note that these parameters are 201 positive, due to the conditions imposed on the original parameters (see Table 1). System 202 (4-6) has the following steady states. The trivial equilibrium point corresponding to 203 tumor elimination is Another equilibrium point, representing the tumor escape, given by Finally, there are also two nontrivial equilibria corresponding to the coexistence between 206 tumor cells, effector, and long-term memory CAR-T cells, given by where ϑ = u − p, with u − p positive due to the first condition imposed on the original 208 parameters, and X 2 and X 3 are the roots of the second-degree equation Assessing the positiveness and stability of the steady states, we found two thresholds (bifurcation points), given by These thresholds determine the following regions in the parameter space where the 211 model presents different dynamic behavior (see Supplementary Material for details): non-negative equilibria are P 0 (which is a saddle point) and P 1 (which is locally 214 asymptotically stable); which are P 0 (saddle point), P 1 (saddle point), and P 3 (locally asymptotically 217 stable); equilibria, which are P 0 (saddle point), P 1 (locally asymptotically stable), P 2 220 (saddle point), and P 3 (locally asymptotically stable).

221
The division of the ϑ × s plane into regions R 1 , R 2 , and R 3 is shown in Figure 2a.

222
(a) ϑ × s plane (b) Phase portrait of the model Figure 2. The CAR-T therapy ODE model presents different dynamical behaviors in each of the three regions R 1 , R 2 , and R 3 , indicated in the ϑ × s plane (a). In the HDLM-2 + CAR-T 123 scenario, the parameter values correspond to region R 3 , and the phase portrait in this case, together with typical model trajectories, is shown in (b). The equilibrium points are indicated by red dots. The yellowish surface represents the separatrix between the basins of attraction of P 1 (escape) and P 3 (stable coexistence). The saddle points are indicated by P 0 and P 2 .
In order to achieve the patient's cure, the system trajectory must be either in the 223 basin of attraction of the tumor elimination equilibrium P 0 or in the basin of attraction 224 a stable coexistence equilibrium where only a harmlessly small amount of tumor cells 225 is present, as described by equilibrium P 3 . Since the point P 0 is always unstable, only 226 the last option is possible, which can be reached in Regions R 2 and R 3 . While in region 227 R 2 the tumor escape equilibrium P 1 is unstable (and all trajectories eventually converge 228 to equilibrium P 3 ), in region R 3 we have bistability between P 3 and P 1 ; in this case, the 229 model outcome (tumor control or escape) depends on the initial conditions. Setting we rely on building a VP that reflects the variability observed in the available data using 243 the strategy similar to that described in [46]. Of note, virtual clinical trials are becoming 244 increasingly popular to represent the heterogeneity of patient cohorts in pharmacology 245 models [47,48]. Here we use the resulting VP to investigate how population heterogene-246 ity impacts overall treatment responses and to identify the most influential parameters 247 for each of them.

248
To build the VP, we first assume that each model parameter is a random variable 249 following a uniform distribution with a wide plausible range. We take random samples 250 from the parametric space, each one representing a plausible virtual mouse. The set of 251 accepted parameters must satisfy the restrictions φ < ρ and − ρ + φ > 0, as indicated 252 in Table 1. Moreover, each physiologically plausible set of parameters is included in 253 the VP only if it leads to a predefined characteristic or behavior similar to that of the 254 target distribution. Specifically, here we build a VP that matches the overall survival 255 of non-treated NSG mice injected with HDLM-2 cells reported in [7]. It means that we    We first simulate the scenario presented in [7], which consists of CAR-T 123 therapy 308 against HDLM-2 cells. Ruella et al. [7] reported that 2 × 10 6 cells of Hodgkin lymphoma presented in [7]. Our simulation also provides the dynamics of memory CAR-T cells.
315 Figure 3a shows that, as the population of C T cells decreases, phenotypic differentiation 316 occurs giving rise to memory CAR-T cells C M . Our simulation shows that effector 317 CAR-T cell populations remain undetectable until t = 250 days, which agrees on results 318 presented in [7]. Moreover, our model indicates the presence of long-term memory 319 CAR-T cells, which slightly decline in time due to a small mortality rate of µ. Model 320 parameters values used in this simulation are displayed in Table 2. remains until the end of simulation on day 500. As explained in [7], tumor rejection 331 occurs due to the re-activation of previously undetectable memory CAR-T cells.

332
Next, we investigate the model behavior in a different scenario with a fast growth 333 tumor cell. The corresponding experiment is described in [11], which uses RAJI tumor 334 and immunotherapy with CAR-T 19 cells. RAJI tumors are much more aggressive than 335 HDLM-2 tumors and express the CD19 antigen. Ninomiya et al. [11] reported that There is an expansion of effector CAR-T cells, which can reduce tumor growth rate but did not eliminate the tumor. Effector CAR-T cells are extinct at the end of the simulation. There is no memory formation. (c) CAR-T 19 immunotherapy on RAJI-IDO + cells. On day 7, 1 × 10 7 CAR-T 19 cells were introduced and were rapidly eliminated; (d) CAR-T 19 immunotherapy with IDO inhibitor (1-MT) shows a restoration of CAR-T cell dynamics, demonstrating the impact of IDO. The parameter α was estimated for these two cases and was responsible to capture the effect of IDO inhibition due to 1-MT. Its value decreased for the RAJI-IDO + + CAR-T 19 + 1-MT case, being small enough to promote a higher expansion of the effector CAR-T cells, and ultimately leading to a more effective control on the tumor growth. However, both therapies were not able either to eliminate the tumor or build memory cells. Dots and standard deviation correspond to experimental data from [11].

Insights on immune checkpoint inhibitors 348
Our model includes the term αTC T in equation (1) Table 2 for the RAJI-control. According 357 to Figure 2C of [11], it should be noted that RAJI and RAJI-IDO + tumor sizes on the  not grow or shrink significantly on day 500; the tumor is reduced to a very small (but 385 not zero) value, which characterizes a state of a residual disease, as depicted in Figure   386 4b. In this immunotherapy outcome, both C T and C M cells are non-zero, and therefore 387 there is the coexistence of the three cell populations. This is a typical configuration of 388 tumor equilibrium, one of three "Es" of immunoediting [53]. Finally, further reducing    [54]. In all cases (d)-(f), the tumor is eliminated in a few days, followed by a decrease of the effector CAR-T cells. Fractionated infusions lead to the formation of memory CAR-T cells, although the quantity depends on the rest time between doses. The next experiment explores the alternative possibility of a fractionated treatment 405 using CAR-T cells, which is a strategy tested in the clinic aiming to reduce toxicity effects 406 [54]. We selected the same scenario described in Figure 3a Figure   414 3a, immunological memory is formed, and the peak of memory cells is similar to that 415 of a single total dose infusion, although a certain delay is observed due to fractionated 416 dose. Such delay ultimately yields a greater formation of immunological memory on 417 day 200. Specifically, the number of memory CAR-T cells at that time is around 7% 418 and 15% larger for 7 and 14 days rest time between doses, respectively. Alternatively, a 419 simulation is performed for fractionated immunotherapy described in [54]. In that work, 420 patients with relapsed or refractory CD19 + ALL were treated with three fractionated T is the initial tumor burden and C T is the CAR-T 123 cell dose injected on day 42. The usual ranges of T and C T were considered, with the number of tumor cells starting at the detectable limit established in [7] and with the maximum CAR-T cell dose corresponding to the highest value used in [7]. Higher doses (C T ≥ 1.3 × 10 6 cells) are able to eliminate any tumor burden smaller than 10 7 cells. It is worth noting that the CR region decreases with the increase of T.
infusions over 3 consecutive days with increasing doses (10%, 30%, and 60%). It was 422 shown that such treatment protocol does not compromise effectiveness while reducing 423 toxicity effects [54]. Figure 4f shows the in silico predictions using this protocol. Like in a   Table 2 for HDLM-2 + CAR-T 123. This range was 434 crucial to obtain our target VP with mean and median survival around 137 and 128 days, 435 as observed in [7]. The survival curves for both control data from [7] and the VP are 436 depicted in Figure 6a, displaying statistically similar mean and median survival times. 437 We then submitted our VP to CAR-T 123 treatment with different doses, varying from 438 1.5 × 10 6 to 1.0 × 10 5 CAR-T cells. While 100% overall survival was reached with mice 439 treated with 1.5 × 10 6 CAR-T 123 cells in [7], the VP reached 95% overall survival in 300 440 days, corresponding to 4754 VM (see Figure 6b). Such a 5% reduction can be explained   [7] (green) and VP of mice engrafted with HDLM-2 tumor (red). (b) CAR-T dose of 1.5 × 10 6 cells led to 95% overall survival in almost one year, 5% lower than those observed in [7] owed to individual parameter uncertainties. Overall survival decreased significantly with dose reduction. Specifically, the survival rate reached 7% when the dose decreased to 1.0 × 10 5 CAR-T cells. The number of VM that survived for 300 days for each dosing strategy is also indicated.

470
CAR-T cell immunotherapies are spreading across hematological cancers and are 471 already products of big pharma companies [55]. On the road, there are new CAR designs, 472 including new antigen targets [6], different CAR affinity [56], and expansion protocols 473 [57]. Mathematical models can be used as accessory tools for new developments [18,19]. 474 Here, we built a three population mathematical model to describe tumor response to 475 CAR-T cell immunotherapy in immunodeficient mouse models (NSG and SCID/beige) 476 based on two published articles from literature [7,11]. Our model was able to represent HDLM-2 target. However, for none of the evaluated RAJI scenarios the formation of a 496 memory pool was observed, due to the rapid growth dynamics of this tumor. 497 We performed in silico studies to highlight how the model could be used as an 498 adjuvant platform to contribute to a better understanding of the underlying processes 499 and for experimental research. Investigating, the application of different dosing pro-500 tocols, we showed that fractionated dose appears to be as effective as a single dose, 501 and the rest periods between infusions might favor long-term immunological memory.

502
These results corroborate previous clinical trials using fractionated CAR-T cell dose 503 with similar effectiveness to single-dose and persistence of CAR-T cells on the blood 504 20 months after therapy [55]. We also found the CAR-T cell dose determination for a 505 given tumor burden is a critical factor for the success of the immunotherapy. A previous 506 model already considered CAR-T cell proliferation in response to antigen burden [26], 507 but memory CAR-T cell was not considered, neither the effect of tumor inhibition of 508 CAR-T cells. A recent paper considered naïve, effector, and long-term memory T cells in 509 a refractory large B cell lymphoma model [10]. We did not include naïve CAR-T cells, 510 because they pass through an in vitro activation protocol, and only activated effector CAR-T cells are present in the treatment [58]. Another interesting mathematical model 512 was made upon tisagenlecleucel-treated patient data [25]. This model was adapted 513 from a previous empirical model of an immune response to bacterial/viral infections.

514
They captured CAR-T cell expansion, contraction, and persistence like our model does, 515 including memory CAR-T cell population. Their model was calibrated on patients' data, 516 and different from ours, no difference in dose-response was detected. They attributed 517 this result to CAR-T cell proliferation capacity in vivo. We partially agree, but there is a  [23]. No spatial distribution was considered in our model, as we are 556 dealing with hematological cancer, but this is required in CAR-T therapy for solid tumors. antigens against glioblastoma [61]. CARs that incorporate multiple target antigens are 561 also the subject of recent research to overcome the mechanism of resistance to CAR-T 562 therapy [13]. Although not completely understood, the incidence of this phenomenon 563 has been linked to antigen escape or lineage switch [62,63] which can be modeled 564 as stochastic events. A recent mathematical model [10] has already pointed out the importance of considering stochastic events to deal with tumor elimination in response to 566 CAR-T cell therapy. It was proposed a hybrid technique that combines deterministic and 567 stochastic events, the latter included only when tumor cells are under a given threshold 568 which ultimately impacts tumor extinction. This strategy reduces the computational 569 burden associated with the higher cost of stochastic models. However, as stochastic 570 events can not be neglected in many situations, further researches are still needed 571 towards accurate and computationally efficient methodologies.