CAR-T Cell Therapy for the Treatment of ALL: Eradication Conditions and In Silico Experimentation

: In this paper, we explore the application of Chimeric Antigen Receptor (CAR) T cell therapy for the treatment of Acute Lymphocytic Leukaemia (ALL) by means of in silico experimentation, mathematical modelling through ﬁrst-order Ordinary Differential Equations and nonlinear systems theory. By combining the latter with systems biology on cancer evolution we were able to establish a sufﬁcient condition on the therapy dose to ensure complete response. The latter is illustrated across multiple numerical simulations when comparing three mathematically formulated administration protocols with one of a phase 1 dose-escalation trial on CAR-T cells for the treatment of ALL on children and young adults. Therefore, both our analytical and in silico results are consistent with real-life scenarios. Finally, our research indicates that tumour cells growth rate and the killing efﬁcacy of the therapy are key factors in the designing of personalised strategies for cancer treatment.


Introduction
Acute leukaemias are classified into myeloid and lymphocytic, depending on the result of the immunophenotype characterisation. Concerning acute lymphocytic leukaemias (ALLs), the latest World Health Organization classification replaced the classic cytological classification into B cell ALL and T cell ALL. Although, both are a malignant transformation and proliferation of lymphoid progenitor cells (B or T lymphoblasts) affecting mainly the bone marrow and peripheral blood. ALL is thought to be the result of a damage to DNA that causes lymphoid cells to undergo a rapidly and uncontrolled growth replacing elements from the bone marrow and other lymphoid organs, ultimately spreading throughout the body [1,2]. Furthermore, leukaemia cells do not have the ability to fight pathogens as well as a normal white blood cell would and possess an increased maximum lifespan, which is one of the main characteristics of a cancer cell.
Although ALL affects mostly children with 80% of all registered cases, the prognosis in paediatric patients is excellent with an overall 5-year survival rate of 69.9% and 90-95% of cases achieving complete response (CR) after chemotherapy. Nonetheless, it represents a devastating disease when it occurs in older adults, where the highest percentages of deaths occur after 65 years of age [2][3][4][5].
Concerning the treatment of ALL [4,[6][7][8], incremental advances in leukaemia therapies have led to marked improvements in survival, abolishing the prognostic influence of clinical and biologic variables that were previously related to a poor outcome. Currently available treatments for ALL include, mainly, chemotherapy and stem cell transplantation. In the particular case of chemotherapy, it consists of induction, consolidation and longterm maintenance. The goal of induction therapy is to achieve CR and restore normal haematopoiesis. Consolidation therapy aims to eradicate the submicroscopic residual disease that remains after CR. Maintenance therapy is the final stage of treatment, this phase has been demonstrated to lower the risk of relapse once CR has been established. However, 15% of children do not respond well to this approach [9], and untreatable relapse remains as a leading cause of cancer deaths in this age group. In adult patients, relapse is more frequent and often unsalvageable. Elderly patients are often unable to tolerate such regimens and carry a particularly poor prognosis. Relapsing patients usually have had already the maximum tolerable dose of chemotherapy. Thus, there is a pressing need for novel strategies in the overall treatment and relapse in ALL.
In recent years, chimeric antigen receptor (CAR) T cells therapy has revolutionised the treatment of resistant haematological malignancies and quickly became a new standard of treatment for relapsed/refractory disease [10]. This therapy has been investigated in both solid and non-solid tumours with promising results due to many studies showing high rates of remission [11][12][13][14][15]. Furthermore, CAR-T cell technology has the potential to be applied for the treatment of numerous diseases such as hemophilia, type 1 diabetes, multiple sclerosis, influenza A, HIV, SARS-CoV-2 and cardiac fibrosis, among others [16]. At a cellular level, CAR-T cells are generated by the T cells from either a patient or a donor's blood. After these cells are expanded and genetically modified they are reinfused into the cancer patients to specifically target and destroy tumour cells [17]. CARs are engineered receptors that can graft an arbitrary specificity onto an immune effector cell [17,18]. Unfortunately, CAR-T cells are associated with serious toxicity that includes cerebral oedema, cytokine release syndrome, on target/off tumour recognition, neurological toxicity, cytopenia and anaphylaxis [1,[19][20][21]. Thus, there are still some challenges that need to be resolved concerning maximum tolerated dose, schedules of application, toxicity management and long-term effects of this treatment.
On the last subject, mathematical modelling can provide a powerful tool to further investigate both short-and long-term effects of CAR-T cells therapy administration. Additionally, mathematical models could be used to solve both the dosing problem and the scheduling for the intervals of application. However, for these models to be applicable in real-life scenarios they need to be formulated from experimental data and/or clinical trials. In the particular case of leukaemia and its treatment, Chulián et al. performed an extensive review concerning several mathematical models of first-order Ordinary Differential Equations (EDOs) in the literature [22]. Among these models, we found two of particular interest where CAR-T cells therapy is applied for the treatment of B cell ALL [23] and T cell ALL [24]. From these models, we formulated two first-order ODEs to explore the periodical applications of CAR-T cells and the corresponding effect on the time-evolution of ALL. In order to validate this approach, we compared our results with those reported by Lee et al. [11] in their clinical trial on CAR-T cells therapy for the treatment of ALL in children and young adults.
In the scope of our system, we were able to solve the dosing problem on the CAR-T cells therapy administration. This was achieved by applying nonlinear systems theories such as the Localization of Compact Invariant Sets (LCIS) [25,26] and Lyapunov's Direct [27,28] methods to establish a sufficient condition on the therapy dose, i.e., the minimum dose required to ensure both CR and complete eradication in the mathematical model under study in this work. Particularly, we designed three personalised therapy administration protocols for the application of CAR-T cells on a hypothetical patient with a given set of characteristics describing its overall health. The latter is illustrated by means of in silico experimentation [29][30][31], which is the process by which multiple numerical simulations are performed to compare different scenarios regarding initial tumour burdens, parameter values, total doses and intervals of therapy applications.
The remainder of this paper proceeds as follows. In Section 2, the necessary background on nonlinear systems theories to determine our mathematical results are presented. In Section 3, the ALL and CAR-T cells mathematical model is thoroughly described. In Section 4, bounds of the localizing domain are computed and sufficient conditions are established on the CAR-T cells therapy dose to ensure nonexistence and eradication of the ALL cancer cells population. In Section 5, three therapy administration protocols are compared and discussed with one of a clinical trial, in silico experimentations illustrate that our mathematical results are consistent with real-life scenarios of cancer treatment. General conclusions of this work are given in Section 6, and we elaborate an Appendix A section to further explore the dynamical properties of the ALL and CAR-T cells system.

Materials and Methods
In this section, we provide the necessary background on nonlinear systems that allows us to derive sufficient conditions to ensure cancer cells eradication in this work. Further, when nonlinear systems theory is combined with systems biology, one can assume that there is a final critical value below which a population of cells cannot longer be considered as biologically meaningful. Therefore, results from previous works have enable us to formulate the next assumption on cell-cell interaction mathematical models. Assumption 1. Cells Eradication Threshold. See Section 4.2 in [31]. If a solution describing the growth of a cell population goes below the value of 1 cell, then it is possible to assume the complete eradication of such population.
Assumption 1 is extremely relevant on the in silico experimentation phase as it is applied to determine the instant on which the cancer cells population is completely eradicated when performing numerical simulations. Therefore, the estimated time provides a threshold to stop the therapy administration protocol.

Localization of Compact Invariant Sets Method
Krishchenko and Starkov proposed the LCIS method in [25,26] to study the short-and long-time dynamics of nonlinear systems of first-order ODEs by computing the so-called localizing domain, which is a bounded region in the state space R n where all compact invariant sets of a system are located.
The method is formulated as follows. Let us take an autonomous nonlinear ODEs system of the formẋ = f (x), where f (x) is a C ∞ -differentiable vector function and x ∈ R n is the state vector. Let h(x) : R n → R be a C ∞ -differentiable function, h| S denotes the restriction of h(x) on a set S ⊂ R n . The function h(x) used in this statement is called localizing and it is assumed that h(x) is not the first integral of f (x). S(h) denotes the set x ∈ R n | L f h(x) = 0 , where L f h(x) represents the Lie derivative (temporal derivative) of f (x) and is given by L f h(x) = (∂h/∂x) f (x). From the latter, one can define h inf = inf{h(x) | x ∈ S(h)} and h sup = sup{h(x) | x ∈ S(h)}. Hence, the General Theorem concerning the localization of all compact invariant sets of a dynamical system establishes the following. Theorem 1. General Theorem. See Section 2 in [26]. Each compact invariant set Γ ofẋ = f (x) is contained in the localizing domain: If the location of all compact invariant sets is inside the domain Λ ⊂ R n , then the set K(h) ∩ Λ may be formulated. It is evident that if all compact invariant sets are located in the sets K(h i ) and K h j , with K(h i ), K h j ⊂ R n , then they are located in the set K(h i ) ∩ K h j as well. Therefore, a refinement of Theorem 1 is realised with help of the Iterative Theorem stated as follows.
Theorem 2. Iterative Theorem. See Section 2 in [26]. Let h m (x), m = 0, 1, 2, . . . be a sequence of C ∞ -differentiable functions. Sets contain any compact invariant set of the systemẋ = f (x) and Nonetheless, if one considers the location of all compact invariant sets inside the domain Λ ⊂ R n and Λ ∩ K(h) = ∅, then it is possible to formulate the following Nonexistence Proposition. Proposition 1. Nonexistence Proposition. See Section 2 in [26]. If Λ ∩ K(h) = ∅, then the systemẋ = f (x) has no compact invariant sets located in Λ.
Equilibrium points, periodic, homoclinic and heteroclinic orbits, limit cycles and chaotic attractors are examples of compact invariant sets. Localizing functions are selected by a heuristic process, this means that one may need to analyse several functions in order to find a proper set that will allow fulfilling Theorem 1.

Lyapunov's Direct Method
Aleksandr Lyapunov [27,28] concluded that a certain type of functions can be analysed to determine stability of an equilibrium point which is denoted as x * ∈ R n and satisfies f (x * ) = 0. Stability in the Lyapunov sense has a central role in the study of autonomous nonlinear ODEs systems of the formẋ = f (x), with x ∈ R n . Lyapunov's stability theorems provide sufficient conditions for stability and asymptotic stability of the equilibrium, both local and global.
In order to apply Lyapunov's direct method, it is necessary to formulate the so-called Lyapunov candidate function, which is usually denoted as V(x) : R n → R, a continuously differentiable function whose temporal derivative is given byV(x) = (∂V/∂x) f (x). This function must be positive definite, i.e., V(0) = 0 and V(x) > 0 for x = 0, while a negative definite function is also V(0) = 0, but V(x) < 0 for x = 0. Further, function V(x) is said to be radially unbounded if V(x) → ∞ as x → ∞.
Now, by considering the information shown above, let us present the following theorem: Theorem 3. Global Asymptotic Stability. See Chapter 4 in [27] and Chapter 2 in [28]. The equilibrium point x * is globally asymptotically stable if there exists a function V(x) positive definite, radially unbounded and decrescent such that its temporal derivativeV(x) is negative definite.
A function V(x) satisfying the properties of Theorem 3 is called Lyapunov function. Nonetheless, there is no explicit method to find a candidate. Therefore, they should be formulated by trial and error. It is also important to mention that the fact that a Lyapunov candidate function fails to satisfy the requirements of Theorem 3 does not imply that the equilibrium must be necessarily unstable, it only means that stability in the sense of Lyapunov cannot be established.

Positiveness of Solutions
A dynamical system is said to be positive [P] if and only if the non-negative orthant R n +,0 is forward invariant, that is, for any non-negative initial conditions all trajectories remain either inside or at the boundaries of R n +,0 for all future times [32,33]. Given an autonomous nonlinear ODEs system of the formẋ = f (x), the following Lemma provides a sufficient and necessary condition to establish its positivity. [32]. The autonomous nonlinear systeṁ x = f (x) is positive if and only if the next holds

Lemma 1. Positiveness of solutions. See Section I I.A in
The latter implies that when evaluating the vector function f (x) at the boundary ∂ of R n +,0 , i.e., f (x)| x=0 , the result must be a non-negative function [ f (x) ≥ 0]. Therefore, given non-negative initial conditions, all solutions will have non-negative real values for all t ≥ 0.

The ALL and CAR-T Cells Mathematical Model
Pérez et al. [24] and León et al. [23] formulated two mathematical models concerning CAR-T cell therapy for the treatment of T-ALL and B-ALL, respectively. These works explored the relation of these two cells populations interactions by means of nonlinear system theory in the form of first-order ODEs. The mathematical properties of each model were thoroughly analysed, including equilibrium points, local stability, the positiveness of solutions and the existence of periodic orbits. Further, parameter values and initial therapy dose incidence on the final dynamics of each system are illustrated and discussed.
The novelty of these systems relies mainly on the modelling of the CAR-T cells evolution. Authors included the following complex phenomena for this population: serial killing of leukaemia cancer cells due to direct encounters; fratricide, mutual killing between the cells preventing generation, expansion, and persistence; clonal expansion due to mitosis stimulation after recognising the target antigen, located in both leukaemia and other CAR-T cells; and a finite lifespan of between two weeks to one month. Additionally, note that the following two assumptions were made: CAR-T cells do not die after killing target cells and they do not expand in vivo. Expansion may be considered in the in vitro phase as cytokines are added externally forcing them to divide. Nonetheless, this process must be modelled separately under different conditions and assumptions.
The ALL cancer cells growth is modelled by means of the exponential law, which is not biologically applicable in the long-term. Therefore, in our approach, we apply a sigmoidal law as this provides an initial exponential growth phase with an upper bound given by the so-called maximum tumour carrying capacity [34]. This implies that a new parameter should be estimated. Concerning the CAR-T cells therapy, only one dose is applied as the initial condition. The latter leads to the conclusion that the initial dose of the therapy does not affect the final outcome on both CAR-T cells and leukaemia populations [23,24]. Therefore, we introduce a protocol administration parameter to consider further applications of the therapy in order to achieve complete eradication of the ALL population on the system. The ALL and CAR-T cells mathematical model is given by the next two first-order ODEs: The dynamic of the CAR-T cells is formulated by Equation (1), and it is described as follows. The first term represents the stimulation to mitosis with a rate ρ C due to encounters with the target antigen in both ALL cancer cells and other CAR-T cells. The natural death rate of these cells is given on the second term by parameter τ C . The α parameter rate represents both fratricide, i.e., third term of Equation (1), and killing efficacy of leukaemia cells from the therapy, i.e., the second term in Equation (2). The value of this parameter is going to be estimated in the Discussion Section from a clinical trial and compared with others available in the literature. Further, of note is the fact that providing a mathematical restriction on the overall proliferation of this cell population. Below, in Table 1, one can see that (3) always fulfils. The therapy administration protocol is given in the fourth term of Equation (1) by parameter ϕ C . Sufficient conditions on the minimum dose will be derived by means of the LCIS and Lyapunov's Direct method. Now, as discussed above, the first term of Equation (2) represents the logistic growth of the ALL cancer cells population with a growth rate ρ L and the inverse of the maximum tumour carrying capacity given by τ L . Units, values and their ranges are shown in Table 1. Further, specific values for each parameter will be thoroughly discussed and justified at Section 5 by taking into account a clinical trial in order to perform the in silico experimentation by means of multiple numerical simulations. The latter is necessary to illustrate our mathematical results and compare them with real-life scenarios.
The range on the maximum tumour carrying capacity τ −1 L is estimated as follows.
The human body contains approximately (3.7 ± 0.8) × 10 13 cells [35,36], where approximately 1.6% represents the population of lymphocytes, see Figure 6.1 in [35]. Thus, by direct proportion, one can compute (5. 92 ± 1. 28) × 10 11 total lymphocytes, and we assume that if the population of leukaemia cells reaches a value within this range, then it will inevitably imply the death of the patient as all healthy lymphocytes would have been replaced by malignant cells. Now, note that by conditions of Lemma 1, the ALL and CAR-T cells system (1) and (2) is positive for any non-negative initial conditions. Therefore, any semi-trajectory is going to be positively forward invariant in the non-negative quadrant R 2 +,0 , and all dynamics are located in the following domain: Furthermore, by Assumption 1, we have the following, and implying that if any solution goes below the value of 1 cell, then it is possible to consider the complete eradication of that population. Now, let us determine the tumour-free equilibrium point of Equations (1) and (2) by equating them to zero and solving for the state variable C(t) when L(t) = 0. The result is shown below, and, as condition (3) always holds, it is evident that the tumour-free equilibrium In Appendix A, all remaining equilibria of the ALL and CAR-T cells system are calculated, i.e., another three equilibriums; conditions for persistence of the ALL cancer population are derived, i.e., necessary conditions to ensure that L(t) > 0; and existence and uniqueness of solutions for any biologically feasible initial conditions are established, i.e., C(0), L(0) ≥ 0.

Localizing Domain and Nonexistence Conditions
In this section, the lower and upper bounds for both CAR-T cells and leukaemia cells populations are formulated by means of three localizing functions. These bounds are given by inequalities in terms of the system parameters. Further, a nonexistence condition of compact invariant sets outside the plane L = 0 (tumour-free equilibrium point) is derived from the upper bound of the ALL cancer cells [L(t)]. This condition is established on the immunotherapy treatment dose [ϕ C ].
Now, let us explore the first localizing function: and compute its Lie derivative as indicated below, from the latter, set S(h 1 ) = L f h 1 = 0 may be written as follows by completing the square now, it possible to estimate a lower bound by rewriting S(h 1 ) thus, results are formulated in the next set it is evident that if the CAR-T cell therapy dose is zero, then the lower bound C inf = 0. Therefore, as long as therapy is being administered to the patient, a concentration of CAR-T cells can be expected in the system. Now, the following localizing function is applied in order to find the upper bound: the Lie derivative is computed as follows: and the next constraint is established, note that the latter always holds, parameter ρ C is a proportion of α (see condition (3) and Table 1). Thus, set S(h 2 ) = L f h 2 = 0 can be written as indicated below, now, by substituting C = h 2 − L into the set, completing the square and performing basic arithmetic operations, we get the following, and the localizing domain of h 2 = C + L is given as follows: therefore, one can estimate the next upper bound of the CAR-T cells Now, by applying the next localizing function, one can determine an ultimate upper bound for the ALL cancer cells population: whose Lie derivative is given by thus, set S(h 3 ) = L f h 3 = 0 is given as follows: at this step, one needs to apply the Iterative Theorem as shown below, from the latter, the next upper bound is determined Results shown above allow us to formulate the next statement: where CAR-T cells : Nonexistence conditions may be computed as a secondary result of Theorem 4 by formulating the following constraint: which is rewritten as follows: and solved for the CAR-T cells therapy parameter [ϕ C ] thus we are able to establish the next statement.

Corollary 1.
Nonexistence. If condition (6) holds, then there are no compact invariant sets for the ALL and CAR-T cells therapy system outside the plane L = 0.

Eradication Conditions
Note that the fulfilment of the nonexistence condition (6) does not imply the leukaemia cancer cells eradication described by the system (1)- (2). Therefore, we apply Lyapunov's Direct Method in order to derive sufficient conditions that can allow us to establish global asymptotic stability of the tumour-free equilibrium point and ensure the complete ALL cancer cells eradication. Thus, let us analyse the following Lyapunov candidate function: and compute its temporal derivative as follows: thus, one can formulate an upper bound by evaluating the function at the localizing domain, i.e.,V Γ CL , and get the next resultV thus, if the next condition holds thenV is negative definite, i.e.,V(0) = 0 andV < 0 ∀ L > 0. Thus, by solving (7) for the CAR-T cells therapy parameter [ϕ C ] we get the following, Therefore, in this particular case, nonexistence and global asymptotic stability conditions are the same for this mathematical model describing cancer evolution. Results shown in this section allow us to formulate the next statement. (1) and (2) is completely eradicated for any initial tumour size [L(0)]. Hence,

Theorem 5. ALL Cancer Cells Eradication. If the CAR-T cells therapy dose [ϕ C ] fulfils condition (6), then the ALL cancer cells population [L(t)] described by the system
and the tumour-free equilibrium point (4) [(C * 0 , L * 0 )] is globally asymptotically stable.
If condition from Theorem 5 holds, then L(t) = 0 and Equation (1) is rewritten as follows: which has a unique biologically meaningful (non-negative) equilibrium given by C * 0 that is globally asymptotically stable, i.e., lim t→∞ C(t) = C * 0 . Additionally, once the ALL cancer cells population is eradicated then the CAR-T cells therapy can be stopped [ϕ C = 0] and C * 0 = 0. This implies that the CAR-T cells population will eventually be depleted in the patient, i.e., lim t→∞ C(t) = 0. Results determined in this section are illustrated and discussed below by means of in silico experimentations.

Discussion and In Silico Experimentation
In this section, our mathematical results are explored by means of the so-called in silico experimentation [30,31] in the form of several numerical simulations under different assumptions regarding parameter values, initial tumour burdens, and therapeutic doses.
The ALL cancer cells eradication condition (6) formulated on the immunotherapy treatment [ϕ C ] (see Corollary 1 and Theorem 5) may be fulfilled for diverse scenarios as it is written in terms of the following three parameters: ρ L , the leukaemia cancer cells growth rate; α, the killing efficacy rate of CAR-T cells; and τ C , its natural death rate. From previous works [30,31,[37][38][39][40], we have found that both the cancer cells growth rate and the killing efficacy of the treatment are consistently involved in tumour clearance conditions. Therefore, scientifically formulating a personalised treatment strategy for each patient requires accurately fitting the values of these two parameters. However, this may be difficult to achieve: clinical trials do not always provide all the required information to estimate these parameters as it is not the main purpose of these studies.
One of particular interest is that of Lee et al. [11], where they performed a phase 1 dose-escalation trial on CAR-T cells for the treatment of B-ALL on 21 patients concluding that this therapy is feasible, safe and mediates potent anti-leukaemic activity in children and young adults. Further, the authors provided important information concerning dose, intervals of the therapy application, thresholds for cancer remission and CAR-T cells detectability in the peripheral blood, overall toxicity and a period for CR. Our main interest relies on the CR period and thresholds for remission and clinical detectability. These are going to be applied with the complete eradication threshold (see Assumption 1) to estimate a set of values for the α parameter, as well as to discuss the feasibility of the system (1)- (2) to reproduce the ultimate ALL cancer dynamics under CAR-T cell therapy.
First, let us summarised Lee et al. methods and results [11]. Using a standard protocol guideline of 3 + 3 to establish the maximum tolerated dose [41], therapy was infused on days 0 and 7. Expansion of CAR-T cells occurs during the first 2 weeks, followed by a rapid CAR-T cell contraction. Dose 1 was 1 × 10 6 per kg and dose 2 was 3 × 10 6 per kg, both of CAR-transduced T cells. Protocol-prescribed doses were successfully produced for 19 of 21 patients for a 90% feasibility rate. Of the 14 responding patients, 12 had undetectable circulating B cells after treatment between days 14 and 28 or shortly thereafter. Information concerning the weight of each patient was not registered in their results.
Therefore, numerical simulations should illustrate that the B-ALL cancer cells population is at least below the remission threshold by day 28 after the initial application of the described two-dose therapy protocol. Nonetheless, we aim for all solutions to be below the eradication value of 1 cancer cell in our in silico experimentation process. First, let us consider the following concerning weight, initial conditions and other parameters shown in Table 1 to estimate the necessary values of α to achieve CR:

1.
Patient weight is increased by 10 kg in each iteration going from 10 to 100 kg.

2.
Six initial tumour burdens of 10 5 − 10 10 are set simultaneously for each weight.

3.
The maximum tumour carrying capacity τ −1 L is set to 4. 64 × 10 11 cells. As shown in Table 1, this is the lower value for this parameter as it is taking into account that B-ALL occurs mostly in children [3] with a median age of diagnosis of 17 years old [4] [23,24].
Numerical simulations were performed by applying Euler's method [42] at Section 1.7] to solve system (1) and (2) with a step size ∆ t = 1 × 10 −5 to further reduce the intrinsic error in Euler's solutions. Table 2 summarises our results for the set parameter values.  The killing efficacy of T cells has been the subject of interest in several works regarding the modelling of cancer evolution and immune response, and the significant results are shown in Table 3.  [51] It is evident that results from Table 2 are close to those presented by Kirschner et al. [43,44] and de Pillis et al. [46,47] when applying an adoptive immunotherapy treatment in the form of tumour infiltrating lymphocytes and Kronik et al. [45] in their clinical trial testing an allogeneic prostate cancer whole-cell vaccine. Therefore, our estimated values of α are consistent with those find in the literature. Now, to further continue with our in silico experimentation in the ALL and CAR-T cells system (1) and (2), a hypothetical 60 kg patient is selected as this is the average weight in both female and male children of 17 years old [52]. The particular characteristics for this patient are set as follows: a CAR-T cells lifetime τ −1 C of 14 days [23,53], and a rapidly growing tumour, i.e., ρ L = 1/20. Therefore, the killing efficacy of the therapy [α] should be 2.154 × 10 −7 to achieve CR by day 28, see Table 2. Further, three initial [L(0)] B-ALL tumour burdens given by 10 7 , 10 8 and 10 9 cells are investigated as most of the patients in the clinical trial of Lee et al. are between these values, see Figure 1 at [11].
First, the four CAR-T cells therapy administration protocols are illustrated in Figure 1. The one defined in the clinical trial performed by Lee et al. [11] is compared with three protocols of our own. Then,   [11] therapy administration protocol, i.e., an initial dose of 6 × 10 7 CAR-T cells at day 0, and a second one of 1.8 × 10 8 at day 7. Figure 2 illustrates results when applying the therapy administration protocol designed by Lee et al. [11]. The protocol is given as follows: an initial dose of 6 × 10 7 at day 0 and a second dose of 1.8 × 10 8 at day 7, which implies that a total of 2.4 × 10 8 CAR-T cells are infused by the end of the treatment. One can see that as the initial leukaemia cancer cells burden increases, then the population takes a few more days to go below the threshold for tumour eradication. The time-evolution of the CAR-T cells remains almost the same through the three iterations with a small peak (lower left panel) due to the high initial tumour burden that stimulates clonal expansion. The latter is observed again in  . The fortnight administration protocol. Two applications were applied, the first one at day 0 and the second at day 14. In silico experimentations allow us to determine that in this approach a dose of 111, 337, 711 CAR-T cells was needed to achieve CR by day 28. Figures 3-5 illustrate results of the three CAR-T cells therapy administration protocols formulated in the form of impulse trains with an amplitude (dose) given by a proportion of the ALL cancer cells eradication condition (6), whose value is computed as follows: where ρ L = 1/20, α = 2.154 × 10 −7 , and τ C = 1/14.
In the first iteration shown in Figure 3, two applications of the therapy are applied at days 0 and 14, i.e., two weeks apart. The necessary dose for each application was estimated as 111, 337, 711 CAR-T cells for the three initial tumour burdens. The latter implies that a total of 222, 675, 422 cells were infused. This value is very close to that obtained by Lee et al. in their clinical trial. In this case, CAR-T cells are below the detectability threshold by day 70. The second iteration considered weekly applications at days 0, 7, 14 and 21 as it is illustrated in Figure 4. In this case, numerical simulations have shown that a dose of 21, 140, 072 CAR-T cells was needed for each application. Therefore, the total amount of cells infused was significantly lower with a final value of 84, 560, 288 CAR-T cells, and the population is below the detectability threshold by day 77.  Figure 5. The daily administration protocol. In this case we aim for CR by day 28 by means of daily therapy applications from day 0 to day 21 as it was the last day of application in the weekly protocol illustrated in Figure 4. The latter allows for one week of rest between the last dose and the peripheral blood test indicated in Lee et al. [11] at day 28.
For the third and final iteration, the therapy was applied every day, and as it is seen in Figure 5, the daily dose was determined to be 2, 226, 755 CAR-T cells, this means that the final amount by the end of day 21 was a total of 48, 988, 610 cells. In this iteration, CAR-T cells are below the detectability threshold by day 77. However, to the best of our knowledge, we should note that this may not be a suitable administration protocol for the treatment in a real-life scenario. Nonetheless, it is evident that shortening the period of the applications ultimately decreases both the dose and the final amount of therapy that needs to be infused into the patient in order to achieve CR.
Concerning the dynamics of CAR-T cells, one can see in  that solutions of C(t) go below the threshold of clinical detectability of 10 4 cells between days 56 and 77 even though the B-ALL cancer cells population is eliminated by day 28. The latter could be explained by the self-stimulation term given by ρ C C × C in Equation (1), indicating that the given value should be lower than the one considered in this and previous works. Nonetheless, further information is required in order to accurately fit parameter ρ C .
Although the results illustrated and discussed in this section are only for the established scenario, note that all sixty values computed for α in Table 2 were estimated in this form. The latter implies that CR by day 28 is achieved in all other cases, and it is straight forward to illustrate that the B-ALL cancer cells population [L(t)] goes below the threshold of tumour eradication for any initial burden in the range of 10 5 -10 10 cells, for any weight going from 10 to 100 kg, and for every combination on the parameter values.

Conclusions
In silico experimentation provides a powerful tool that could potentially be applied in the design of personalised strategies for cancer therapy administration protocols. Numerical simulations may help us to observe beyond the thresholds of clinical detectability in current cancer imaging technologies. Furthermore, this approach allows researchers to explore multiple scenarios at the same time when taking into account different possibilities on the health of the patient. The latter due to the solutions illustrating the efficacy of the cancer therapy in both short-and long-term, information that can be used to decide the best treatment strategy for each case.
In this work, we were able to formulate three CAR-T cell therapy administration protocols (see Figure 1) and compare them with the overall results reported by Lee et al. when exploring a phase 1 dose escalation trial for the treatment of B-ALL. As illustrated in Figures 2-5, CR was achieved in all strategies where different doses and intervals of application were tested. The amount of each dose was computed as a proportion of the ALL cancer cells eradication condition (6), this was established in Corollary 1 and Theorem 5 by means of the LCIS and Lyapunov's direct methods. Therefore, we can conclude that the combined application of nonlinear systems theory and systems biology with in silico experimentation can provide useful information in the designing of cancer therapy administration protocols.
Regarding the dose of the therapy, Figures 3-5 gradually illustrate that more applications of the therapy decrease the total amount of CAR-T cells that needs to be infused into the patient for CR. Nonetheless, the feasibility of these strategies remains as an open question that needs to be further investigated and validated by clinical trials.
Finally, our results indicate that factors such as tumour growth, immune response and the efficacy of the therapy in eliminating malignant cells, among many others, can be translated into specific parameter values that ultimately yield different tumour dynamics. Therefore, by accurately estimating these data, mathematical models that can better describe real-life scenarios of cancer evolution of either solid or non-solid tumours would eventually be formulated and validated at the same time. In the specific case of leukaemia, blood cell counts can be determined by means of peripheral blood flow cytometry which is performed by Lee et al. in their clinical trial. Nonetheless, measurements are acquired for a very limited number of days, and they are not identified for each patient. However, by incorporating daily blood cell counts to these types of clinical trials to properly mea-sured each population, these data can be modelled through genetic algorithms focusing on objective functions such as those from growth laws. Furthermore, parameter values can be estimated with 95% confidence interval through nonlinear curve-fitting [54].

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: remembering that (C * 0 , L * 0 ) was introduced at the end of Section 3 with the tumour-free equilibrium point (4), and α > ρ C as condition (3) always holds. Therefore, the following assertions can be made.
Now, given that f i (t, x), i = 1, 2; and x = [C, L] T , then the Jacobian matrix (see [42] at Section 7.4) [∂ f /∂x](t, x) of (1) and (2) is computed as indicated below and it is evident that both f i (t, x) and [∂ f /∂x](t, x) are continuous and exist on the domain Ω = [t 0 , t 1 ] × Γ CL with [t 0 , t 1 ] ∈ [t 0 , ∞) and Γ CL ⊂ R 2 +,0 . Regarding existence and uniqueness, the latter implies that f i (t, x) is locally Lipschitz in x on Ω, see Lemma 3.2 by Khalil in [27] at Section 3.1. Further, it is possible to ensure that each element of (A4) is bounded by Theorem 4, and we conclude the following. Theorem A1. Existence and uniqueness. There is a Lipschitz constant ≥ 0 such that [∂ f /∂x](t, x) ≤ on Ω. Then, f i (t, x) satisfies the Lipschitz condition and there exists some δ > 0 such that the ALL and CAR-T cells system (1) and (2), given aṡ x = f i (t, x) with x(t 0 ) = x 0 , has a unique solution over [t 0 , t 0 + δ].
In order to derive necessary conditions for ALL cancer cells persistence in the solutions of system (1) and (2), the local stability of the tumour-free equilibrium point (4) should be investigated. Hence, let us evaluate this point at the Jacobian matrix (A4) as follows: and as (A5) is an upper triangular matrix, eigenvalues are given by each element of the main diagonal, see Section 6.4 in [42]. Thus, by substituting C * 0 and performing the corresponding arithmetic operations, results are indicated below: Theorem 4.7 by Khalil in [27] at Section 4.3 allows us to conclude that (4) is locally asymptotically stable if Reλ i < 0, i = 1, 2. Therefore, one can solve λ 2 < 0 for the immunotherapy parameter [ϕ C ] as follows: and a necessary condition for tumor persistence (see Liu and Freedman in [55] at Section 3.1.3) is given by and we formulate the next statement.
Theorem A2. Persistence. If condition (A6) on the CAR-T cells therapy dose holds, then the ALL cancer cell population described by the system (1) and (2) persists.
The latter implies that all solutions L(t) > 0 for all t ≥ 0 as the CAR-T cells therapy dose will not be sufficient to completely eradicate the leukaemia cancer cells population. When substituting the parameter values from the scenario defined in Section 5 one gets the following value ϕ inf = 25, 866 CAR-T cells, which is lower than the therapy dose that was established by means of condition (6) to ensure the complete eradication of the tumour population, i.e., ϕ C > ϕ CART = 28, 187 CAR-T cells, see Corollary 1 and Theorem 5.