Chemoimmunotherapy Administration Protocol Design for the Treatment of Leukemia through Mathematical Modeling and In Silico Experimentation

Cancer with all its more than 200 variants continues to be a major health problem around the world with nearly 10 million deaths recorded in 2020, and leukemia accounted for more than 300,000 cases according to the Global Cancer Observatory. Although new treatment strategies are currently being developed in several ongoing clinical trials, the high complexity of cancer evolution and its survival mechanisms remain as an open problem that needs to be addressed to further enhanced the application of therapies. In this work, we aim to explore cancer growth, particularly chronic lymphocytic leukemia, under the combined application of CAR-T cells and chlorambucil as a nonlinear dynamical system in the form of first-order Ordinary Differential Equations. Therefore, by means of nonlinear theories, sufficient conditions are established for the eradication of leukemia cells, as well as necessary conditions for the long-term persistence of both CAR-T and cancer cells. Persistence conditions are important in treatment protocol design as these provide a threshold below which the dose will not be enough to produce a cytotoxic effect in the tumour population. In silico experimentations allowed us to design therapy administration protocols to ensure the complete eradication of leukemia cells in the system under study when considering only the infusion of CAR-T cells and for the combined application of chemoimmunotherapy. All results are illustrated through numerical simulations. Further, equations to estimate cytotoxicity of chlorambucil and CAR-T cells to leukemia cancer cells were formulated and thoroughly discussed with a 95% confidence interval for the parameters involved in each formula.


Introduction
Leukemias are non-solid tumors usually known as cancer of the blood. These hematologic malignancies can be classified as acute or chronic depending on the time-evolution of the disease, and as myelogenous or lymphocytic as they arise from the dysfunctional proliferation of developing leukocytes (white blood cells). Hence, the main classifications are as follows: Acute Myelogenous Leukemia (AML), Chronic Myelogenous Leukemia (CML), Acute Lymphocytic Leukemia (ALL), and Chronic Lymphocytic Leukemia (CLL) [1]. However, regardless of the case, the lifespan of leukemia cancer cells is longer than normal cells and they do not have the ability to fight pathogens as effectively as a normal white blood cell could.
Concerning overall statistics, the Global Cancer Observatory (GLOBOCAN) estimated 474,519 new leukemia cases in 2020 with 311,594 deaths in that year. Here, it should be noted that GLOBOCAN groups leukemias from C91 to C95 (according to the International Classification of Diseases) [2,3]. Nonetheless, lymphocytic (C91) and myelogenous (C92) are the most worrying, where ALL is most common in pediatrics, CLL in the elderly, AML in adults, and, although CML remains the least common, it typically affects older adults and rarely occurs in children.
In the particular case of CLL, since it is the topic of interest in this research, is most frequently diagnosed among people aged 65-74, with a median age at diagnosis of 70 years old, it is rarely seen in people under the age of 40, and is extremely rare in children [4,5]. The diagnosis is established by blood counts with at least 5000 monoclonal lymphocytes per mm 3 (or µL) for a minimum of three consecutive months; blood smears to look for abnormalities in appearance, number, and shape in these white blood cells; and immunophenotyping of circulating B-lymphocytes, which identify a clonal B-cell population carrying tumor-specific antigens as well as typical B-cell markers [4,6]. CLL is a disease whose progression responds to a clinically heterogeneous picture, making the existence of a curative treatment difficult as the affected individuals have comorbidities due to their advanced age. The latter also implies the appearance of periodic sequences of both response and relapse in many patients. Five-year overall survival rate is estimated in the range of 23-93%, with an overall survival of 2 to more than 20 years. However, survival statistics are based on large groups of people, they cannot be used to predict exactly what will happen to one individual as treatment and response can vary greatly between two patients [5,7,8].
As of today, CLL remains incurable with conventional therapies, and disease progression is inevitable. Concerning chemotherapy treatment regimens for CLL, this may be conducted by several drugs, most notably chlorambucil in the elderly as it is not as toxic as fludarabine, cyclophosphamide, and rituximab. Although chlorambucil remains the treatment of choice for this disease it has shown limited efficacy [8]. Recent clinical advances are aiming for personalized therapy strategies as the new path to follow in cancer treatment. In patients with certain hematologic malignancies such as ALL and CLL, the use of autologous T cells genetically modified to express chimeric antigen receptors (CARs), such as the CD-19, has led to unprecedented clinical responses opening the door to a new era of personalized cancer therapy. Anti-CD19 CAR-T cells may be manufactured from both CD4 + and CD8 + T cell subsets to treat adults with relapsed or refractory CLL [9][10][11].
Despite the excitement around CAR-T cells for the treatment of hematologic malignancies, this therapy has come under criticism for its cost, which in the case of the most recently CAR T-cell therapy approved by the Food and Drug Administration (FDA) is more than USD 450,000 [12]. Nonetheless, clinical trials around the world have been developed to better understand and optimize the application of CAR-T cells [13][14][15][16][17]. On the latter, mathematical and computational modeling coupled with in silico experimentation and nonlinear dynamical systems theories may be a powerful tool in designing personalized chemoimmunotherapy treatment strategies, computer simulations are intended to replace the laborious efficacy testing in real humans and reduce the likelihood of drug failure [18][19][20].
In recent years, nonlinear system theories such as the Localization of Compact Invariant Sets (LCIS) method, the direct and indirect methods of Lyapunov, and in silico experimentation have been applied to gain some insights in the designing of personalized schedules for the administration of chemotherapy and immunotherapy. However, combined application strategies are still to be explored, which is the main objective of this work. The latter is performed in a mathematical model composed of four first-order ODEs describing the dynamics of alive and dead leukemia cancer cells, CAR-T cells as the immunotherapy treatment, and chlorambucil as the chemotherapy drug. We were able to develop and illustrate diverse protocols of chemoimmunotherapy treatment to control cancer growth and achieve the complete eradication of the leukemia cells in the proposed mathematical model.
The remainder of this paper proceeds as follows. In Section 2, the mathematical model of leukemia and chemoimmunotherapy is formulated, values and units of the parameters are given, and both biological and mathematical assumptions are described. In Section 3, we formulate two equations to estimate cytotoxicity to cancer cells of the chlorambucil chemotherapy drug and the CAR-T cells immunotherapy, as well as how to compute the number of molecules from chlorambucil dose in mg, its concentration in the blood system, and its decay rate by assuming a first-order pharmacokinetics. Further, we present the localizing domain, following with global asymptotic stability conditions that imply the leukemia cells eradication, and persistence of both CAR-T cells and cancer cells. In Section 4, results are illustrated by means of the in silico experimentation, where different cases with treatments and without treatments are explored. In Section 5, both mathematical and numerical simulation results are discussed. Finally, in Section 6, we present the conclusions of our work.

Mathematical Model: Leukemia and Chemoimmunotherapy
The leukemia and chemoimmunotherapy mathematical model is formulated by considering the dynamics between CLL cancer cells [A(t)], dead leukemia cells [A d (t)], chemotherapy as chlorambucil [C(t)], and immunotherapy in the form of CAR-T cells [T(t)] in the blood circulatory system by the following four first-order ODEs: where the time unit is considered as days; A(t), A d (t) and T(t) are given in cells; and C(t) is measured in molecules. The concentration of cells and molecules in the system is computed when assuming 5 L of blood in a human adult. Further, all solutions with non-negative initial conditions will be located in the non-negative orthant (see Section II.A in [40]): Interactions among live and dead leukemia cells as well as the chemotherapy units in molecules were introduced by Guzev et al. in their work concerning the experimental validation and in silico experimentation to describe cytotoxicity of three different chemotherapy compounds in leukemia cells [25]. Pérez et al. [27] and León et al. [28] formulated and explored the dynamics of CAR-T cells in the treatment of T-cell and B-cell leukemias, respectively. The latter was further discussed by Valle et al. in [29]. Now, let us discuss interactions in each ODE as follows. The growth of CLL cancer cells [A(t)] in the blood circulatory system is defined by the logistic law in the first term of Equation (1), whereas eradication of these cells is considered by the law of mass action when interacting with dead leukemia cells in the second term; cytotoxicity to leukemia cells from the chemotherapy drug and the CAR-T cells is modeled by the Michaelis-Menten kinetics and the law of mass action in the third and fourth terms, respectively. Further, the maximum tumor carrying capacity is estimated by considering that lymphocytes rep-resent 1.6% of the total amount of cells in a human adult with around 4.5 × 10 13 total cells [41,42]. The evolution of leukemia dead cells [A d (t)] is described by Equation (2), the first term represents dissolution of these cells in the system; remaining terms of this equation indicate accumulation by the eradication of living leukemia cells, the second term describes apoptosis or necrosis as a result of living cells competing for nutrients in the blood circulatory system, while the third and fourth terms are the lysed leukemia cells due to the log-kill effect of the chemotherapy and the CAR-T cells immunotherapy, respectively. Pharmacokinetics of the chemotherapy drug [C(t)] is given in Equation (3), where the first term is the decay rate of the drug, this is computed and further explained in Section 3 from its reported biological half-life, term two was introduced by Guzev et al. when considering the total number of molecules attacking each cancer cell; therefore, we incorporated the third term by assuming that the chemotherapy will inevitably have a cytotoxicity effect on the CAR-T cells when combining these two treatments. Constant or periodic administrations of the chemotherapy can be considered by the fourth term. The dynamics of CAR-T cells [T(t)] are formulated in Equation (4), the first term represents stimulation to mitosis, i.e., activation and proliferation, due to encounters with the target antigen in leukemia cancer cells and other CAR-T cells; the second term indicates natural death and/or exhaustion; fratricide or self-inhibition is formulated by the third term; the fourth term considers the cytotoxicity of the chemotherapy in the case that both treatments are present in the system at the same time; and further infusions of the therapy are considered with the fifth term. Description and values of parameters of the leukemia and chemoimmunotherapy mathematical model (1)-(4) are shown below in Table 1.
Now, by the cells eradication threshold assumption [20,29,43] that establishes the following: "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", one can formulate the next constraints for all solutions of the leukemia and chemoimmunotherapy system (1)-(4): Live CLL cancer cells : Molecules of the chemotherapy drug : Concerning the overall dynamics of CAR-T cells, it should be noted that the parameters regarding their killing efficacy rate (α T ) and the mitosis/activation rate (ρ T ) have the following constraint providing a mathematical restriction on their proliferation once infused into the system, see [27][28][29]. This condition is directly related to the ultimate bounds of the localizing domain and the leukemia eradication by the immunotherapy treatment. Further, there is a threshold for the numerical value of ρ T that implies either depletion or persistence of the CAR-T cells, and as it is shown below, this value is given by a condition in terms of the death rates of both leukemia (µ A ) and CAR-T cells (µ T ), and the dissolution rate of dead leukemia cells (µ d ) in the blood circulatory system mathematical background on these conditions will be provided and discussed in Section 3.
Regarding the depletion/persistence phenomenon, both circumstances have been reported in in vivo clinical studies. Lee et al. [44] analyzed peripheral blood by flow cytometry and qPCR in 21 patients and reported that CAR-T cells were no longer detected (fewer than 10 4 absolute cells) by day 68 after the last infusion. Porter et al. [7] detected, also by means of flow cytometry and qPCR, persistence of CAR-T cells in the range of 14 to 49 months in the four patients who achieved CR (complete response). However, we have identified that recent studies are leaning to the side of CAR-T cells long-term persistence as new generations of these cells are being developed [11,[45][46][47].

Estimating Cytotoxicity of Chemotherapy and CAR-T Cells to Cancer Cells
In this section, nonlinear functions are formulated to estimate the cytotoxicity to cancer cells from different doses of chlorambucil chemotherapy drug and the environment surrounding CAR-T cells immunotherapy. First, let us take the results shown in Table A.1 in Guzev et al. [25]; particularly, chlorambucil concentration (presented in µM), and its corresponding cytotoxicity rate (denoted as µ AC ) on CLL cells. As one can see below in Table 2, increasing the drug concentration also increases cytotoxicity. Further, concentration was converted from µM (10 −6 mol/L) to mg/L, and cytotoxicity rate from h −1 to days −1 to standardize the time units of the mathematical model under study in this work. Table 2. Chlorambucil concentration and its computed and approximated cytotoxicity rates (µ AC ) on leukemia cancer cells. Data from columns 'Concentration in µM' and 'Experimental cytotoxicity in h −1 ' were extracted from Appendix A in Guzev et al. [25]. It is evident that chlorambucil cytotoxicity rates are not increasing linearly. Thus, by following the log-kill hypothesis, which states that when cancer volume is increasing by a constant fraction of itself every fixed unit of time and it is in the presence of effective anticancer drugs, then it will shrink by a constant fraction of itself [48]. In other words, a therapy dose eradicates a constant proportion of a tumor cell population rather than a constant number of cells. Further, the total amount of eradicated cells will be given by a logarithmic function in base 10. Hence, inspired by the "log" part of the log-kill hypothesis we formulate the next nonlinear function to fit the µ AC values from the column 'Experimental cytotoxicity in days −1 ' in Table 2,

Concentration
where κ is the chlorambucil concentration in mg/L, and µ AC (κ) represents its approximated cytotoxicity rate in days −1 . Experimental and fitted data is illustrated in Figure 1, and the corresponding values are shown in the column 'Fitted cytotoxicity [µ AC (κ)] in days −1 ' in Table 2 with a coefficient of determination R 2 = 0.999. Further, parameters ε 1 and ε 2 from Equation (8) were computed with a 95% confidence interval (CI) as it is shown in Table 3.  The mathematical model constructed by Guzev et al. [25] describes the dynamics of the chemotherapy drug in molecules per unit of time. Therefore, let us determine the total number of molecules from its usual dose given in mg as follows [49][50][51]: Chlorambucil average molecular weight mg mol , where Avogadro's number = 6.02214076 × 10 23 molecules mol , and Chlorambucil average molecular weight = 304.212 × 10 3 mg mol .
Now, let us explore results shown above by considering the usual protocol of chlorambucil dosing in previously untreated adult patients with CLL. According to the Access Pharmacy Educational Resource [52], off-label dosingis given as indicated below 0.4 mg kg × day 1 every 2 weeks to a maximum of 0.8 mg kg × day with 24 cycles.
Thus, by assuming a weight of 80 kg (although, in a real-life scenario, the particular weight of each patient can be used) for a CLL adult patient with an average age of ∼70 years [5, 53,54], then, therefore, the total number of molecules can be computed by Equation (9) representing the total number of molecules in 1 mg of chlorambucil. At this step, Equation (9) may be simplified where C mg is the chlorambucil dose in mg. Hence, for 32 mg we have the following result N.Molecules = 32 × ς = 6. 33467792 × 10 19 molecules.
Regarding Equation (8), it is important to remember that chlorambucil cytotoxicity was investigated by Guzev et al. in micromolar drug concentrations, see column 'Concentration in mg/L' in Table 2. Hence, given that an adult has about 5 L of blood [42], in order to properly apply this equation one needs to compute the concentration in the blood circulatory system in mg per liter of blood by dividing as indicated below therefore, one can estimate the cytotoxicity of this concentration of chlorambucil to leukemia cells as follows µ AC (κ) = log 10 1 + 4.799 × 10 −2 (6.4) Now, given that the biological half-life (t 1/2 ) of chlorambucil has been reported as 0.0625 days [50,51], one can determine its deactivation/decay rate (denoted as µ C ) from the first-order pharmacokinetics Equation [55] x = −µ C x, with the following solution for the initial condition x(0) = C mg . Thus, from the latter we have the following where C mg /2 represents 50% of the initial dose at t 1/2 . Therefore, by isolating µ C we obtain the next deactivation/decay rate for chlorambucil Concerning cytotoxicity of CAR-T cells, Kiesgen et al. [56] made a review study on CAR-T cell-mediated cytotoxicity and concluded that the antitumor efficacy of this 'living drug' is influenced by their activation, proliferation and inhibition, as well as their exhaustion. The first three are directly related to their consecutive encounters with cancer and other CAR-T cells, whereas T-cell exhaustion has been reported in many chronic diseases such as human cancer [57]. Furthermore, we were able to estimate cytotoxicity rate values from the Lee et al. [44] phase 1 dose-escalation trial on CAR-T cells for the treatment of leukemia on 21 patients by also considering cancer growth rates. Our results from Table 2 in [29] are summarized below.
Fitted cytotoxicity values of CAR-T cells from column 'Fitted cytotoxicity Table 4 were computed with a coefficient of determination R 2 = 0.995 by applying the next Equation which was formulated by taking into account both the mitosis stimulation rate (ρ T ) and natural death rate of CAR-T cells (µ T ), as well as the leukemia growth rate (ρ A ). Further discussion on these parameters is provided in Table 1 in Section 2. Parameters γ i , i = 1, 2, 3; were computed with a 95% CI as it is shown in Table 5. Table 4. CAR-T cells cytotoxicity rates (denoted as α T and given in days −1 ) estimated by Valle et al. [29] from the Lee et al. [44] phase 1 dose-escalation trial on 21 leukemia patients.
shows the approximated values with Equation (12).

CAR-T Cells Mitosis CAR-T Cells Natural Leukemia Cells Experimental Fitted Cytotoxicity Stimulation Rate
However, accurately estimating cytotoxicity of CAR-T cells remains as an open problem. Additional clinical data is needed to either validate and/or reformulate Equation (12). Following this statement, we provide in the Supplementary Material the necessary code to perform this task in Anaconda 3 with the Jupyter Notebook 6.4.5 (using Python 3.9.7), as all parameters from Equations (8) and (12) were fitted by applying the curve_fit function from scipy.optimize [58] with the SciPy Version 1.7.1.

Localization of Compact Invariant Sets
Below, we provide the mathematical background that allows us to determine the localizing domain in the non-negative orthant R 4 +,0 where all compact invariant sets of the leukemia and chemoimmunotherapy system (1)-(4) are located. These compact invariant sets can include equilibrium points, periodic orbits, limit cycles and chaotic attractors, among others as illustrated in [59] at Section 3. The so-called General Theorem concerning the LCIS method was formalized by Krishchenko and Starkov in [60] (see Section 2) and it states the following: Furthermore, if all compact invariant sets are contained in the set K(h i ) and in the set K h j then they are contained in K(h i ) ∩K h j as well. Nonexistence of compact invariant sets can be considered for a given set Λ ⊂ R n if Λ ∩ K(h) = ∅, then the systemẋ = f (x) has no compact invariant sets located in Λ.
Now, let us explore five localizing functions in order to formulate lower and upper bounds for a localizing domain containing all compact invariant sets of the system (1)-(4) in R 4 +,0 . Maximum upper bound for the live leukemia cells population [A(t)]. In order to determine this bound, the following localizing function is proposed and its Lie derivative is defined as follows at this step, negative terms of S(h 1 ) can be discarded. Hence, non-negative boundaries for all ω−limit sets of A(t) are given in the next domain However, it should be noted that the latter was determined when all the interactions with the other cells and treatments are neglected, thus, the maximum carrying capacity (which is also an equilibrium of A(t)) is located at the upper boundary. Below we will explore another localizing function that considers the effect of the dead leukemia cells in the overall cancer cells population.
Localizing set for the live and the dead leukemia cells populations [A(t) and A d (t)]. Let us explore the following localizing function whose Lie derivative is defined as follows then set S(h 2 ) = L f h 2 = 0 can be formulated and simplified by basic arithmetic operations as follows when considering that and by completing the square in the right side of the equation we obtain the next result now, due to h 2 = A + A d we may conclude the following result Thus, from the latter, upper bounds for both the live [A(t)] and dead [A d (t)] leukemia cells populations can be deducted as follows Localizing set for the CAR-T cells therapy [T(t)]. Two localizing functions are analyzed to determine these bounds. First, let us take the next to formulate the lower bound whose Lie derivative is defined as follows thus, by condition (5), i.e., α T > ρ T , we can complete the square hence, the lower bound is determined by disregarding the non-negative function f 1 (A, C, T) and isolating T, then, we can conclude on the following lower bound for all solutions of T(t) from the latter, it is evident that T inf | ϕ T =0 = 0, which is expected in the case where the CAR-T cells therapy is not constantly or periodically applied into the system. Concerning the upper bound, we explore the following localizing function and by computing its Lie derivative one can define set S(h 4 ) = L f h 4 = 0 and write it as follows by considering now, let us complete the square therefore, we have the following upper bound for the localizing function h 4 Now, one can formulate the following lower and upper bounds for all solutions of the live leukemia [A(t)] and CAR-T [T(t)] cells populations Localizing set for the chemotherapy drug concentration [C(t)]. Let us exploit the localizing function denoted by h 5 = C where its Lie derivative is defined as follows now, set S(h 5 ) = L f h 5 = 0 can be written as indicated below hence, by applying the Iterative Theorem the next lower bound can be determined from the latter, it is evident that the next condition should be fulfilled for C inf > 0. However, if condition (13) does not hold, then one should consider C inf = 0 as there is no biological sense for negative values for the chemotherapy drug concentration, C(t). Now, let us take again set S(h 5 ) as follows from which we can determine the next upper bound by disregarding the negative terms as given below Hence, all solutions C(t) will be bounded by the following domain where ϕ C ≥ 0 as this term represents the application of the chemotherapy drug into the patient.
Nonexistence conditions and supreme upper bound for the live leukemia cells population [A(t)]. Now, let us apply the Iterative Theorem to the set S(h 1 ) as follows from the latter and by assuming (5) is fulfilled, the next ultimate upper bound can be established Furthermore, one can formulate the following nonexistence condition from A sup as indicated below A sup < 0, which is solved for the CAR-T cells therapy treatment (ϕ T ) Therefore, results shown in this section allow us to conclude the following two statements: Theorem 1. Localizing domain. If condition (5) holds, then all compact invariant sets of the leukemia and chemoimmunotherapy system (1)-(4) are located within the following domain: with the next lower and upper bounds Corollary 1. Nonexistence. If condition (14) fulfills, then there are no compact invariant sets outside the plane A = 0 for the leukemia and chemoimmunotherapy system (1)-(4).

Global Asymptotic Stability: Leukemia Cells Eradication
When considering that Equations (1)-(4) describe cancer evolution as a nonlinear dynamical system, one can apply Lyapunov's Direct method (see Chapter 4.1 by Khalil in [61] and Chapter 2 by Hahn in [62]) to investigate the global asymptotic stability of the tumor-free equilibrium point and establish sufficient conditions to ensure the leukemia cells eradication. The latter may be translated into the real-world as immunotherapy doses that could potentially control cancer cells growth. First, let us take the following Lyapunov candidate function (V) V = A, and compute its time derivative V as followṡ Then, in order to fulfill Lyapunov's asymptotic stability conditions [V(0) = 0 anḋ V < 0 ∀ A > 0], the derivative is evaluated at the localizing domain, i.e.,V K Γ , and the next upper bound is determinedV hence, by assuming (5) holds, we formulate the following condition and solve it for the CAR-T cell therapy parameter as shown below Therefore, the nonexistence condition (14) also becomes a sufficient condition for asymptotic stability and the following statement is concluded: (14), then one can ensure the complete eradication of the leukemia cancer cells population described by the system (1)-(4). Thus,

Theorem 2. Leukemia cancer cells eradication. If the CAR-T cells therapy dose meets condition
Now, one can investigate the leukemia and chemoimmunotherapy system when the live leukemia cancer cells are completely eradicated and both therapies are stopped, i.e., [A(t), ϕ C , ϕ T = 0]. Therefore, Equations (1)-(4) become as followṡ where the only biologically meaningful equilibrium point is the leukemia-free state given by hence, the next statement is a direct result from Theorem 2:  (19) is globally asymptotically stable.

Persistence: Leukemia and CAR-T Cells
As stated before, long-term persistence of CAR-T cells in cancer patients has been reported in several papers. Thus, this particular phenomenon can be studied as a local asymptotic stability problem by means of Lyapunov's Indirect method (see Chapter 4.3 by Khalil in [61]) when considering interactions between leukemia and CAR-T cells as a non-linear system. First, let us simplify the leukemia and chemoimmunotherapy system (1)-(4) as followsȦ The latter represents the original system in the short-term as the biological half-life (t 1/2 ) of chlorambucil is 0.0625 days, which implies that this drug should be depleted in the patient in a short period of time, i.e., C(t) = 0. Hence, one can identify the following equilibrium point in order to investigate persistence of CAR-T cells in the mathematical model after the last infusion, i.e., ϕ Now, as the method requires, the Jacobian matrix [∂ f (x)/∂x] is computed below which is evaluated at the Equilibrium point (23) as follows at this step, eigenvalues are computed from the Jacobian determinant det J| (A * 1 ,A * d1 ,T * 1 ) − λI = 0 , results are shown below and by condition (24) it is evident that Reλ i < 0, i = 1, 2. Therefore, conditions for depletion (6) and persistence (7) are derived from λ 3 as follows. If Reλ 3 < 0, then by Theorem 4.7 in [61]: "The equilibrium is locally asymptotically stable", this holds when which biologically implies depletion of the CAR-T cells. Thus, if Reλ 3 > 0, then by Theorem 4.7 in [61]: "The equilibrium is unstable", which holds when and according to Liu and Freedman (see Section 3.1.3 in [63]) this could be biologically interpreted as a "necessary condition for the cell population to grow". Therefore, the following statement is concluded: Theorem 3. CAR-T cells persistence. If conditions (7) and (24) hold, then the CAR-T cells immune response to the leukemia cancer cells persists after at least one infusion into the system, i.e., lim t→∞ T(t) > 0 when A(t) > 0, for T(0) > 0 and/or ϕ T (τ) > 0 with τ ∈ [t 1 , t 2 ]. Now, following the latter, let us explore leukemia cells persistence under the immunotherapy treatment, i.e., ϕ T > 0. It is important to note that this implies the CAR-T cells dose will not be enough to control cancer growth. Hence, in this case one should investigate the local stability of the next equilibrium point from the simplified system (20)- (22) where α T > ρ T by condition (5). Thus, the equilibrium (26) is evaluated at the Jacobian matrix (25) as indicated below as the resulting matrix is lower triangular, all eigenvalues are given by each element of the main diagonal. Therefore, and it is evident that λ j < 0, j = 5, 6. Thus, local asymptotic stability of (26) follows from the next condition on λ 4 which is solved for the immunotherapy treatment parameter In the biological sense, if condition (27) holds, then the immunotherapy could be able to eradicate a sufficiently small initial tumor population, i.e., the equilibrium point (26) is locally asymptotically stable. Hence, given the following condition the next statement can be concluded regarding the persistence of the leukemia cells population in the system (1)-(4): (28), then the immunotherapy treatment will not be able to eradicate the leukemia cells population. Therefore, cancer persits, i.e., lim t→∞ A(t) > 0 ⇔ ϕ T < ϕ prstnc .

Results: In Silico Experimentation
In this section, we will explore by means of numerical simulations the overall dynamics of the leukemia and chemoimmunotherapy system (1)-(4). It is important to note that our mathematical model aims to describe the evolution of leukemia cancer cells in the blood circulatory system in a hypothetical adult patient when considering the application of immunotherapy in the form of CAR-T cells either alone or combined with the chemotherapy drug chlorambucil. Hence, we formulate four scenarios for the in silico experimentation, i.e., no treatments, CAR-T cells depletion and persistence, and both immunotherapy and chemoimmunotherapy protocols. Now, it should be noted that the in silico experimentation was performed in Matlab 2022a in a desktop computer with a Ryzen 7 5800X CPU, 64 GB of RAM DDR4 3200, and a 2 TB M.2 Samsung 980 PRO SSD. The system of ODEs (1)-(4) was solved by means of Euler's method, that is with a step size ∆ t = 1 × 10 −7 to further reduce the intrinsic error in the system solutions.

No Treatments
First, let us consider the 'no treatments case', i.e., C(t), T(t) = 0. Thus, the mathematical model becomes as followsȦ where the leukemia cells growth rate (ρ A ) is set to 0.1680. This represents 10% of the rate value of Guzev et al. as this was estimated in vitro in ideal conditions. Further, values of this order of magnitude have been reported in several works concerning mathematical models formulated with in vivo studies of both solid and non-solid tumors [24,[26][27][28][34][35][36].
Results are illustrated in Figure 2 with the following initial conditions A(0) = 10 10 cells, which are used through all the in silico experimentation, i.e., Figures 2-8. The latter are set to consider a high initial non-solid tumor population and the corresponding 5% of dead leukemia cells proposed by Guzev et al. [25]. A(t) [cells] [cells] Leukemia cells K(h 2 )   [cells] T (t) T 1 Figure 3. CAR-T cells depletion case. Numerical simulations allow us to illustrate that this phenomenon arises when the mitosis/activation rate (ρ T ) of CAR-T cells is below a threshold directly related to the death rates of both leukemia (µ A ) and CAR-T cells (µ T ), and the dissolution rate of dead leukemia cells (µ d ) as given by condition (6   ρ A = 0.1680, ρ T = 2.9409 × 10 −13 , µ T = 1/14, and α T = 1.2560 × 10 −7 by Equation (12). For the values of µ AC and µ TC see columns 'Cytotoxicity µ AC (κ)' and 'Cytotoxicity µ TC ' of Table 6.

CAR-T Cells Depletion and Persistence
Now, regarding depletion and persistence of CAR-T cells, Figures 3 and 4 illustrate, respectively, conditions (6) and (7), i.e., 'CAR-T cells depletion case' and 'CAR-T cells persistence case'. For these two scenarios Equations are as followṡ where ρ T = 2.5714 × 10 −13 for depletion, and ρ T = 2.9409 × 10 −13 for persistence. These values are directly related to the dissolution and dead rates of both leukemia and CAR-T cells as indicated below CAR-T cells depletion case : Further, by following the Lee et al. [44] escalation trial the immunotherapy dose was set to 2.4 × 10 8 CAR-T cells as this is the last infusion applied when considering an 80 kg patient, hence is the initial condition for the in silico experimentation performed in Figures 3 and 4.

CAR-T Cells Treatment Protocols
In the following two cases we explore the leukemia cells eradication when applying the immunotherapy treatment in two different schedules, each one with its corresponding dose. These scenarios are identified as the 'weekly CAR-T cells application case' and the 'fortnight CAR-T cells application case' as illustrated in Figures 5 and 6, respectively. Numerical simulations are performed by solving the next set of Equationṡ The in silico experimentation allowed us to design two immunotherapy administration protocols with the following characteristics: Mathematically, the first dose is considered with the initial condition T(0), whereas the last three consecutive doses were performed with the treatment parameter ϕ T in the form of a delayed pulse train with asymmetrical waves. The experimentation indicates that delaying doses implies an increase in the total number of CAR-T cells that needs to be infused into the patient to achieve cancer eradication, and it should be noted that each immunotherapy dose fulfills the leukemia cancer cells eradication condition (14) by several orders of magnitude as when considering the next set of parameter values: ρ A = 0.1680, ρ T = 2.9409 × 10 −13 , µ T = 1/14, α T = 1.2560 × 10 −7 by Equation (12), and µ TC = 0 as the chemotherapy treatment is not applied in these two cases. Regarding the cytotoxic effect of the therapy, one can see in the upper panel of Figures 5 and 6 that each application produces a 2 to 3 log-kill of leukemia cells. Nonetheless, numerical simulations also illustrate that cancer cells could begin to grow again if the treatment is stopped. These two cases allow us to conclude that when immunotherapy is the only treatment applied, then reducing the period between applications yields a better result concerning the number of CAR-T cells needed to control and eradicate the leukemia cells population described by the mathematical model under study in this work. Further, the lower panel in Figures 5 and 6 shows that CAR-T cells population is eventually depleted once leukemia cancer cells have been eradicated after the last application of the therapy. The latter is to be expected as the system becomes as followṡ when A(t), A d (t) = 0, and this equation has only one biologically feasible equilibrium point given by T * = 0, which is globally asymptotically stable. Hence, long-term persistence of CAR-T cells could be related to the survival of a small tumor population.

Chemoimmunotherapy Treatment Protocols
Now, regarding the combined chemoimmunotherapy treatment strategy, we formulated the next two scenarios: 'constant chemoimmunotherapy case' and 'increasing chemoimmunotherapy case'. As it was stated in Section 3, we considered an 80 kg CLL adult patient with an average age of ∼70 years. In silico experimentations of these two cases were performed with the complete mathematical model (1)-(4), i.e., In order to properly combine the two therapies, we continue with the fortnight CAR-T cells application protocol, i.e., one application every two weeks, and our aim was to eliminate the fourth dose illustrated in Figure 6. First, constant dose applications of the chlorambucil chemotherapy drug were explored. Numerical simulations allowed us to conclude that two chemotherapy administrations after each immunotherapy infusion with a dose of 0.7 mg/kg were necessary to achieve the leukemia cells eradication, as it is shown in Figure 7.
Concerning the increasing dose of chemotherapy, by means of the in silico experimentation we were able to formulate the protocol illustrated in the lower panel of Figure 8, which is as follows Dose 1 at day 6 with 0.5 mg/kg → C mg = 40 mg, Dose 2 at days 10 and 20 with 0.6 mg/kg → C mg = 48 mg, Dose 3 at days 24 and 34 with 0.7 mg/kg → C mg = 56 mg, Dose 4 at day 38 with 0.8 mg/kg → C mg = 64 mg.
From the latter, one is able to compute the number of molecules of chlorambucil (11), the concentration in the circulatory system in mg per liter of blood (κ), and both cytotoxicity to cancer cells (8) and CAR-T cells as indicated below in Table 6. Now, the leukemia cells eradication condition (14) from Theorem 2 changes as the chemotherapy cytotoxicity increases.
342, 971 CAR-T cells, 345, 131 CAR-T cells, 347, 066 CAR-T cells,  [20,29,43]. The cytotoxicity of the chemotherapy drug to the CAR-T cells was estimated by considering the results reported by de Pillis et al. in [36] as a proportion of the cytotoxicity to cancer cells, i.e., µ TC = 0.054µ AC . Hence, as the chlorambucil dose is increased, then the values of µ AC and µ TC increase as well, this is shown in columns 'Cytotoxicity µ AC (κ)' and 'Cytotoxicity µ TC ' of Table 6. These increments over time in the chemotherapy doses and parameter values were incorporated in the in silico experimentation illustrated in Figure 8.

Discussion
First, concerning the analytical results, our methodology is as follows. The LCIS method was applied to determine ultimate bounds to all solutions for the leukemia and chemoimmunotherapy system (1)-(4) as given in the localizing domain K Γ (see (15)). From the latter, nonexistence conditions for the tumor population can be derived. Then, following these results one can establish the sufficient condition (14) on the immunotherapy treatment to ensure the complete eradication of the leukemia cancer cells by means of Lyapunov's direct method. This condition was only formulated for the CAR-T cells as there are already well established off-label dosing protocols for chlorambucil administration with which the immunotherapy treatment was intended to be combined. Nonetheless, conditions on both therapies could be explored through this analytical procedure. In addition to the global asymptotic stability conditions of the tumor-free equilibrium point (19), local stability conditions were calculated with Lyapunov's indirect method to investigate the long-term persistence of both leukemia and CAR-T cells, as given by (28) and (7), respectively.
All conditions derived in this research are given in terms of the system parameters. These are summarized below Localizing domain and stability (5) : α T > ρ T , Depletion of CAR-T cells (6) : Persistence of CAR-T cells (7) : Eradication of leukemia cells (14) : Persistence of leukemia cells (28) : from the latter, the first condition (5) is directly related to the boundedness, local and global asymptotic stability of the system, and it implies that the killing efficacy rate (α T ) of CAR-T cells should be larger than its activation rate (ρ T ). Concerning depletion or long-term persistence of CAR-T cells we found that this phenomenon is directly proportional to the death rates of leukemia (µ A ) and CAR-T cells (µ T ), and inversely proportional to the dissolution rate in the blood circulatory system of dead leukemia cells (µ d ). At the same time, persistence of CAR-T cells could be associated with the existence of a small population of undetectable cancer cells. Now, the last two constraints, (14) and (28), provide sufficient and necessary conditions to ensure the complete eradication of the leukemia cells by the immunotherapy treatment, and a threshold on which the CAR-T cells dose will not be able to control any initial size of the non-solid tumor population.
Regarding the in silico experimentation performed in this work, the leukemia and chemoimmunotherapy system (1)-(4) was explored under different scenarios. First, numerical simulations illustrate that in the absence of therapies, leukemia cells massively accumulate in the blood circulatory system reaching values of 50,000 cells/µL as shown in Figure 2, which is expected to be lethal in CLL patients. Depletion and persistence of CAR-T cells after the last infusion into the system was found to be strictly related to their activation rate (ρ T ). Both scenarios are illustrated in Figures 3 and 4, where recent studies demonstrate that long-term persistence of CAR-T cells is to be expected as new generations continue to be developed and explored in in vivo clinical trials. In Figures 5-8 we design four administration protocols of immunotherapy and chemoimmunotherapy that completely eradicate the CLL cancer population. Here, it is important to remember that we assume that a cells populations described by an ODEs system are eliminated once the corresponding solution goes below the threshold of 1 cell. The latter is essential to properly apply nonlinear systems theory for modeling the dynamics between cells populations, since any value less than one will not represent any biologically meaningful real-life scenario.
When only the immunotherapy treatment is considered, the in silico experimentation illustrates that reducing the period between applications yields a better overall outcome by improving the in vivo toxicity profile of this so-called living drug. Figure 5 shows that by applying the therapy weekly, fewer cells are needed to eradicate an initial leukemia population of 10 10 cells, whereas by increasing the period between applications, a higher dose of immunotherapy should be infused. Furthermore, there is an instant in time when the leukemia population starts to grow again as there are not enough CAR-T cells to control the tumor population, which is illustrated in Figure 6. Given these results, we decided to explore the chemoimmunotherapy scheme for the 'fortnight CAR-T cells case' aiming to stop the regrowth of cancer cells by including the chlorambucil drug between the immunotherapy administrations. Additionally, we are taking advantage of the period between doses to avoid, at least to some extent, the cytotoxicity of the chemotherapy on CAR-T cells. Hence, the in silico experimentation was performed in order to design the combined administration protocols illustrated in Figures 7 and 8. As one can see, we were able to discard the fourth CAR-T cells application by incorporating two chlorambucil intakes after each dose. Both the constant and increasing dosing scenarios demonstrated similar results concerning the time at which complete eradication of leukemia cancer cells was achieved and the total dose of chlorambucil administered to the hypothetical CLL patient.
Another contribution of this work is the formulation of two equations to estimate the cytotoxicity to cancer cells of the chemotherapy drug chlorambucil (8) and the CAR-T cells immunotherapy (12), to the best of our knowledge, equations of this form have not been proposed before. Equation (8) was constructed by fitting a base 10 logarithmic function to the data from the in vitro study of Guzev et al. where cytotoxicity directly depends on the concentration in mg/L of the drug; whereas Equation (12) provides a minimum value for the cytotoxicity of CAR-T cells that could be enhanced by considering the rates of activation and exhaustion of these cells, as well as the tumor growth rate. However, it is important to note that this equation was fitted to the data of a previous work and further research is needed to either validate or reformulate this equation. Additionally, parameters from Equations (8) and (12) were computed with a 95% CI, and they fit the corresponding cytotoxicity data with coefficients of determination R 2 equal to 0.999 and 0.995, respectively.
Nonetheless, it is important to discuss that even if mathematical models are formulated by considering several aspects of the biological or physiological phenomenon under study, which is also known as mechanistic modeling, they still could be considered as an ideal representation of such phenomenon. Hence, time series data from in vivo clinical studies where the evolution of each cell population could be accurately measured or estimated is needed to validate these models or to better fit the values of the proposed parameters in the system. Following this path, we provided in the Supplementary Material the necessary code to fit real-life data to nonlinear equations in Python.

Conclusions
With the continuous advancement of computing power in recent decades, coupled with the increasingly affordable prices, the paradigm of exploring complex biological phenomena such as cancer evolution through mathematical modeling and in silico experimentation has provided interesting and promising results in this field. Particularly, we applied nonlinear system theories and combine them with numerical simulations to explore several scenarios of CLL progression when considering two anticancer therapies: CAR-T cells and chlorambucil. Our methodology allowed us to design treatment protocols for the administration of immunotherapy and chemoimmunotherapy that completely eradicate the leukemia cancer cells population described by our proposed mathematical model.
We expect this research to be useful in the designing of administration protocols for cancer treatment as the in silico experimentation illustrates that mathematical modeling and nonlinear system theories could be applied to obtain insights on tumor evolution and the cytotoxicity of combined anticancer therapies.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/pharmaceutics14071396/s1. Author Contributions: Conceptualization, P.A.V. and R.G.; methodology, P.A.V. and L.N.C.; software, P.A.V. and R.G.; validation, Y.S., L.N.C. and C.P.; formal analysis, P.A.V. and C.P.; investigation, Y.S. and C.P.; resources, P.A.V., Y.S. and L.N.C.; writing-original draft preparation, P.A.V. and R.G.; writing-review and editing, Y.S., L.N.C. and C.P.; visualization, P.A.V. and R.G.; project administration, P.A.V. All authors have read and agreed to the published version of the manuscript. Acknowledgments: Raul Garrido thanks the National Council of Science and Technology of Mexico (Consejo Nacional de Ciencia y Tecnología, CONACyT) for the financial support during his postgraduate studies under his scholarship agreement number 801678. Finally, a special thanks to Eli, a long-time colleague and Engineer, who selflessly taught us the basics of Python.

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

Abbreviations
The following abbreviations are used in this manuscript: