Combined Application of CAR-T Cells and Chlorambucil for CLL Treatment: Insights from Nonlinear Dynamical Systems and Model-Based Design for Dose Finding

1. Introduction

When addressing the most difficult challenges faced by the scientific community, cancer stands at the forefront. Among its various manifestations, blood cancers include Chronic Lymphocytic Leukemia (CLL), which is one of the most prevalent subtypes in Western countries, accounting for $20\%$ to $30\%$ of all leukemia cases. In 2019 alone, over 40,000 deaths were attributed to CLL worldwide. The incidence and mortality rates are notably higher in males and in regions such as North America, Europe, and Australia. Additionally, the global burden of CLL may be underestimated due to limited case reporting in underdeveloped countries [1].

As a blood cancer involving non-solid tumors, Chronic Lymphocytic Leukemia (CLL) is characterized by the accumulation of malignant B lymphocytes in the peripheral blood, bone marrow, lymph nodes, and spleen. Regarding risk factors, there is no evidence linking environmental causes; however, genetic and familial predispositions are well-documented. Patients undergoing therapies such as chemotherapy, immunotherapy, or their combination generally exhibit favorable survival rates. Nevertheless, cases of aggressive CLL necessitate referral to specialized oncology centers. Hematologists must adhere to standard guidelines to determine the optimal timing for initiating therapy and selecting treatment options [2]. Advances in CLL treatment have led to the development of specialized and targeted therapies [3], with chimeric antigen receptor cell therapy (CAR-T) emerging as a highly promising option for high-risk patients who develop resistance to novel agents.

CAR-T cells for cancer treatment involves the precise modification of a patient’s immune system to enable it to recognize and attack cancer cells. This therapy requires a comprehensive understanding of the patient’s immune health and includes the in vitro culture of the patient’s own T cells. Additionally, patients must undergo lymphodepletion to reduce cell population competition and enhance the therapy’s effectiveness. Preparing T cells to target specific tumor cells also necessitates specialized equipment [4]. CAR-T therapy has emerged as a chemo-free alternative; however, one of its primary challenges is resistance. As a result, physicians often adopt combined therapy approaches [5]. Among available therapies, chlorambucil remains a standard due to its low toxicity, affordability, and ease of oral administration. However, when used as a single agent, chlorambucil has shown limited efficacy [6].

Further investigations have been crucial in understanding the evolution of CAR-T cells therapy combined with chlorambucil. Despite the challenges, health specialists, in collaboration with mathematical biologists, have already developed strategies to enhance the effectiveness of therapies and prevention efforts [7]. These advancements have been made possible through the application of mathematical tools, as modeling via Ordinary Differential Equations (ODEs), which captures the dynamics of this disease with remarkable detail.

With advanced tools in traditional mechanistic modeling, researchers are gaining deeper insights into tumor growth, innovative therapies, prevention strategies, early detection, and more [8]. Furthermore, mathematical biologists are integrating cutting-edge technologies to tackle the complex behavior of cancer. In this context, digital twins have emerged as a transformative approach. These high-precision models, implemented on high-performance computing platforms, enable in silico experimentation, drastically reducing reliance on expensive and time-consuming in vitro or in vivo experimental setups.

Digital twins not only enhance the precision of clinical trials but also allow for the exploration of alternative pathways that are challenging or even impossible to study in real-world scenarios with remarkable accuracy [9]. Additionally, the successful application of mathematical biology to the treatment and understanding of cancer has given rise to Systems Pharmacology. This interdisciplinary approach provides specific tools to optimize clinical practice. By leveraging data collected from in vitro or in vivo experiments, in silico models can meet the stringent requirements of governmental health agencies, paving the way for highly personalized therapies with impressive outcomes [10]. Within this framework, model-based design emerges as a key methodology, offering a structured approach to therapy development by systematically integrating mathematical models with clinical and experimental data. By iteratively exploring different treatment scenarios, model-based design allows for precise optimization of therapeutic protocols [11,12]. The synergy between model-based design, mechanistic models, and Systems Pharmacology enables the development of tailored strategies that address the unique complexities of each individual patient [13].

This work is outlined as follows. Section 2 presents a detailed description of our mechanistic model formulated to explore both short- and long-term effects of a combined chemoimmunotherapy strategy involving CAR-T cells and chlorambucil for the treatment of CLL. Additionally, parameter values and assumptions are provided, along with corresponding references and the mathematical oncology background. In Section 3, nonlinear system theories are applied to analyze the local and global dynamics of the system. Results allow us to establish sufficient conditions for therapy dosing and to design appropriate treatment protocols. Section 4 focuses on in silico experimentation, where we illustrate and compare the results of dose-escalation and dose de-escalation protocols, as well as scenarios in which therapies may fail to completely eradicate the cancer, and eventually leading to relapse of the disease. Conclusions of this work are presented in Section 5. Finally, an Appendix is included with notations and mathematical foundations related to nonlinear system theories.

2. Materials and Methods

In this section, we formulate our mechanistic model to explore the effects of both single and combined applications of chemotherapy and immunotherapy on CLL cancer cells within the bloodstream compartment. The mathematical foundations of the model are presented along with the biological implications of each term. Additionally, the numerical values of all parameters are derived, and the corresponding constraints are established to illustrate different clinically expected behavior of tumor-immune dynamics.

2.1. The CLL Mathematical Model

The mathematical model proposed in this work aims to explore the response of CLL cells $\left[x\left(t\right)\right]$ to the combined application of CAR-T cells $\left[y\left(t\right)\right]$ and chlorambucil $\left[z\left(t\right)\right]$ as a chemoimmunotherapy treatment strategy. The system is formulated by the following three nonlinear ODEs:

$$\begin{array}{ccc}\hfill \dot{x}& =& {\kappa }_{1}x\left(1-{\displaystyle \frac{x}{{\kappa }_2}}\right)-{\kappa }_{3}xy-{\displaystyle \frac{{\kappa }_{4}xz}{{\kappa }_5+z}},\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \dot{y}& =& \psi -{\psi }_{1}y+{\psi }_{2}xy-{\displaystyle \frac{{\psi }_{3}yz}{{\psi }_4+z}},\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \dot{z}& =& ϵ-{ϵ}_{1}z,\hfill \end{array}$$

(3)

where the time-evolution $\left(t\right)$ is measured in days and interactions among variables are modeled within the bloodstream compartment. Hence, to estimate the total concentration of cells and drugs we will assume a blood volume of 5 liters in the human body (see Chapter 6 in [14]). Descriptions and units of all parameters are given in Table 1, while their values and references will be discussed in the next section.

Now, let us describe the system as follows. Equation (1) models the dynamics of CLL cells $\left[x\left(t\right)\right]$ under the combined treatment of chemoimmunotherapy. The first term describes the development of leukemia cells using a logistic growth law (see Section 4.3 in [15]), which implies a growth rate and a maximum carrying capacity. The second term represents the overall immune response of CAR-T cells against the leukemia cells based on the law of mass action (see Section 2.3 in [16]). The third term accounts for the chemotherapy effect through a pharmacodynamics of maximum response [17]. In the scope of our model, CAR-T cells immunotherapy is considered as the primary treatment strategy, hence, we assume an eradication of cancer cells that follows the log-kill hypothesis [18], i.e., each consecutive application of the treatment should eradicate a fraction of leukemia cells regardless of its total concentration. Chlorambucil drug is used as a secondary treatment, therefore, its cytotoxicity on the leukemia cells is modeled by the Michaelis-Menten kinetics [19]. Equation (2) models the overall behavior of CAR-T cells $\left[y\left(t\right)\right]$. The first term accounts for the external application of the effector cells, while the second term captures their natural death. The third term describes activation and proliferation of CAR-T cells upon interaction with leukemia cells in the bloodstream. The fourth term reflects the cytotoxic effect of the chlorambucil drug on CAR-T cells when both therapies are administered into the system at simultaneous periods of time. Equation (3) models the pharmacokinetics of chlorambucil $\left[z\left(t\right)\right]$ as a first-order process measured in $mg$. The first term represents the external administration of the chemotherapy, while the second term accounts for its natural decay over time [20]. The schematic representation of the system is illustrated in Figure 1.

Concerning the overall dynamics, one may apply the positivity property for nonlinear systems established by De Leenheer & Aeyels in [21] and conclude the following:

Remark 1.  

Positive system. All forward solutions $\varphi \left(t,{\tau }_0\right)$ of the system are nonnegative, and all semi-trajectories $\varphi \left(·,\tau \right)$ of Equations (1)–(3) will be positively forward invariant for nonnegative initial conditions $\left({\tau }_0\ge 0\right)$ inside the following orthant $\tau ={\left[x,y,z\right]}^T$:

$${\mathbf{R}}_{+,0}^3=\left\{x\left(t\right),y\left(t\right),z\left(t\right)\ge 0\right\}.$$

Now, let us consider the cell eradication threshold assumption proposed in [22] which states: “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”. This leads to following remark:

Remark 2.  

Clinically significant biological behavior. All solutions for the CLL cells, $x\left(t\right)$, and CAR-T cells, $y\left(t\right)$, will exhibit biologically meaningful dynamics if their populations persist above the value of 1 cell. Otherwise, one can assume complete eradication. Hence,

$$x\left(t\right)=0\forall x\left(t\right)<1,$$

and

$$y\left(t\right)=0\forall y\left(t\right)<1.$$

These remarks will be helpful in drawing conclusions from the in silico experimentation results. First, it is essential to recognize that with positive initial conditions, representing biologically feasible cell and drug concentrations, all solutions will remain within the nonnegative orthant ${\mathbf{R}}_{+,0}^3$. Additionally, Remark 2 offers a practical threshold to assess cell eradication, emphasizing the interplay between nonlinear systems theory and oncology in the context of a mathematical model describing the combined application of chemotherapy and immunotherapy treatments.

2.2. Parameter Estimation

Altough system (1)–(3) was formulated based on a set of well-established assumptions from the mathematical oncology literature, the values for each parameter were carefully selected or estimated to accurately represent tumor dynamics within the bloodstream of the human body.

First, let us discuss parameters from Equation (1) which models the dynamics of CLL cells $\left[x\left(t\right)\right]$. The leukemia cells growth rate, ${\kappa }_1$, was derived from the experimental validation perform by Guzev et al. in [23], authors determined the growth rate of A20 murine leukemic cells to be $0.07$$h^{-1}$, which, when converted to days, yields a value of $1.680$$day{s}^{-1}$. It should be noted that this rate was estimated from cells isolated in vitro, i.e., without the influence of the immune response or therapeutic interventions. Additionally, in vivo studies conducted by the authors estimated a lower growth rate of $0.01$$h^{-1}$$\left(0.240day{s}^{-1}\right)$ [24]. A20 cells have been widely used to study the tumor microenvironment in different tissues such as the spleen, brain, and eyes [25], as well as their ability to infiltrate other organs like the liver, lungs, spleen, and bone marrow [26]. Hence, the growth rate of CLL cancer cells can be defined within the following interval:

$${\kappa }_1\in \left[0.240,1.680\right]day{s}^{-1}.$$

The maximum leukemia cells carrying capacity, ${\kappa }_2$, is derived as follows. Previous studies have estimated around of $1.8-2\times {10}^{12}$ immune cells in the human body [27,28]. Among these, four main lymphocyte types can be considered: T cells, B cells, NK cells, and plasma cells. According to Sender et al. [27], the total number of B cells is estimated to be $3\times {10}^{11}$ with a $95\%$ confidence interval of $\left[2\times {10}^{11},4\times {10}^{11}\right]$. Therefore, we set the value of ${\kappa }_2$ as follows:

$${\kappa }_2=3\times {10}^{11}cells.$$

Porter et al. published a case report in [29] for a clinical trial involving a patient diagnosed with stage I CLL. In this trial, they designed a lentiviral vector expressing a CAR specific to the B cell antigen CD19. Some key results include the following: CAR-T cells expanded to over $1,000$ times the initial level in vivo; complete remission was achieved and sustained for at least 10 months post-treatment; CAR-T cells persisted at high levels in both the blood and bone marrow for 6 months; a more than 3-log expansion of CAR-T cells was observed by day 21 post-infusion; the doubling time of CAR-T cells in the blood was approximately $1.2$ days, with an elimination half-life $\left(t_{1/2}\right)$ of 31 days. Assuming an exponential decay, the death rate of CAR-T cells, ${\psi }_1$, is estimated as follows:

$${\psi }_1={\displaystyle \frac{ln2}{t_{1/2}}}={\displaystyle \frac{ln2}{31}}=22.360\times {10}^{-3}day{s}^{-1}.$$

Similarly, the decay rate of the chlorambucil drug, ${ϵ}_1$, which has a half-life of $1.5$ hours [30], is calculated as:

$${ϵ}_1={\displaystyle \frac{ln2}{t_{1/2}}}={\displaystyle \frac{ln2}{1.5/24}}=11.090day{s}^{-1}.$$

Regarding the activation rate of CAR-T cells due to encounters with CLL cells, ${\psi }_2$, our mathematical analysis (discussed further in the next section) leads to the following conclusion:

Remark 3.  

Persistence of CAR-T cells.Long-term persistence of CAR-T cells in the system (1)–(3) is ensured if and only if the following condition is satisfied:

$${\psi }_2>{\displaystyle \frac{{\psi }_1}{{\kappa }_2}},$$

(4)

otherwise, in the absence of continuous immunotherapy administration, these cells will eventually be depleted.

Therefore, to fulfill condition (4), the value for the activation rate of CAR-T cells is set to at least ${\psi }_2=1.1\times {\psi }_1/{\kappa }_2$ in order to illustrate CAR-T cells persistence in the in silico experimentation. Now, concerning the cytotoxicity effect of both CAR-T cells$\left({\kappa }_3\right)$ and chlorambucil $\left({\kappa }_4\right)$ on leukemia cells, we have formulated the following two equations, recently published in [31]:

$$\begin{array}{ccc}\hfill {\kappa }_3& =& {\gamma }_1+{\displaystyle \frac{{\psi }_1{\psi }_2}{{\gamma }_2+{\psi }_1{\psi }_2}}+{\displaystyle \frac{{\kappa }_1{\psi }_2}{{\gamma }_3+{\kappa }_1{\psi }_2}},\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill {\kappa }_4& =& {log}_{10}{\left(1+{\epsilon }_1\varkappa \right)}^{1/{\epsilon }_2},\hfill \end{array}$$

(6)

with

$$\varkappa ={\displaystyle \frac{ϵ}{5L}},$$

where $\varkappa$ represents the chlorambucil concentration in $mg/L$, and the values of the parameters for Equations (5)–(5) are provided below in Table 2.

Hence, by substituting the corresponding values for ${\psi }_1$ and ${\psi }_2$, the value of ${\kappa }_3$ may be approximated as follows for both the minimum and maximum growth rate of CLL cells:

$${\kappa }_3=1.256\times {10}^{-7}{\left(cells\times days\right)}^{-1}.$$

The values for the cytotoxicity rate of the chemotherapy on CLL cells, ${\kappa }_4$, will be further discussed in the in silico experimentation section, as different concentrations of the chlorambucil drug $\left(\varkappa \right)$ will be used in the protocol administration design, and results from Guzev et al. [23] indicate that cytotoxicity increases when higher doses of the drug are applied, i.e., greater cell growth inhibition. However, as indicated in off-label dosing guidelines [32,33], total concentrations values should be set within the following range:

$$ϵ\in \left[0.4,0.8\right]{\displaystyle \frac{mg}{day}}\mathrm{per}1kg\mathrm{of}\mathrm{the}\mathrm{patient},$$

for up to 24 cycles. Furthermore, as it has been vastly established in the mathematical oncology literature [19,34], the chemotherapy drug should exhibit a greater capacity for killing cancer cells than effector cells for the treatment to be successful. The latter is to be expected as cancer cells proliferate at a much faster rate than effector cells, hence, the drug must be more effective in killing these malignant cells than healthy or immune cells. Therefore, we set the following constraint for the chlorambucil cytotoxicity rate on CAR-T cells, ${\psi }_3$:

$${\psi }_3<{\kappa }_4,$$

(7)

and by following results from de Pillis et al. [19] we set ${\psi }_3=0.054\times {\kappa }_4$.

Now, concerning the chlorambucil concentration that produces $50\%$ of the maximum response, ${\kappa }_5$, we found an in vivo study conducted by Silber et al. [35], in which the chemosensitivity of lymphocytes from patients with B-cell CLL to chlorambucil was determined. The study involved 65 patients diagnosed with CLL, each with at least $5\times {10}^3/\mu L$ of monoclonal B lymphocytes positive for CD5, CD19, and CD20 markers. The overall results are summarized as follows: The median $I{C}_{50}$ for the total patient group was 61 $\mu mol/L$. This value was used as a threshold to categorize the B lymphocytes as sensitive $\left(I{C}_{50}<61\mu mol/L\right)$ or resistant $\left(I{C}_{50}\ge 61\mu mol/L\right)$. Therefore, the value of ${\kappa }_5$ can be calculated in $mg$ as $I{C}_{50}\times C_{mol}\times 5$L, where the molecular weight of chlorambucil is $C_{mol}=304.212\times {10}^3$$mg/mol$ [36]:

$$\begin{array}{ccc}\hfill {\kappa }_5& <& 92.785mg\mathrm{for}\mathrm{the}\mathrm{sensitive}\mathrm{CLL}\mathrm{cancer}\mathrm{cells},\hfill \\ \hfill {\kappa }_5& \ge & 92.785mg\mathrm{for}\mathrm{the}\mathrm{resistant}\mathrm{CLL}\mathrm{cancer}\mathrm{cells}.\hfill \end{array}$$

For the concentration of chlorambucil that produces $50\%$ of the maximum cytotoxicity on CAR-T cells, ${\psi }_4$, we follow the same line of thougth as in condition (7) and establish the next condition:

$${\psi }_4>{\kappa }_5,$$

(8)

based on the assumption that a higher dose of the chemotherapy drug is required to achieve greater cytotoxicity in CAR-T cells. Similar results were reported by Guzev et al. [24].

3. Results: Nonlinear Dynamical Properties of the System

In this section, we investigate local and global dynamics of Equations (1)–(3) by applying nonlinear system theories such as the Localization of Compact Invariant Sets (LCIS) method to determine lower and upper bounds on all forward solutions of the CLL mechanistic model; Lyapunov’s direct and indirect methods to establish sufficient conditions for eradication of leukemia cells and persistence in both CAR-T cells and cancer cells; and Lipchitz conditions for existence and uniqueness of solutions. Fundamentals and mathematical preliminaries are shown in Appendix A.

3.1. Localizing Domain

The localizing domain containing all compact invariant sets of the system characterizes all biologically meaningful behaviors that may emerge during tumor evolution, such as various concentrations of long-persisting tumor burdens, tumor-free states, periods of cancer cells dormancy followed by heightened effector cells activity, and active immune responses exhibiting abrupt changes in effector and cancer cell concentrations, among other dynamics. Now, by following the LCIS methodology one can establish the next result:

Result I: Clinically significant bounds.  Bounds regarding all compact invariant sets of the CLL system (1)–(3) describing clinically significant behavior when applying either single or combined therapy are given within the next domain:

$$K_{xyz}=K_x\cap K_y\cap K_z,$$

where each set is defined as follows:

$$\begin{array}{ccc}\hfill K_x& =& \left\{0\le x\left(t\right)\le x_{sup}\right\},\hfill \\ \hfill K_y& =& \left\{y_{inf}\le y\left(t\right)\le y_{sup}\right\},\hfill \\ \hfill K_z& =& \left\{z_{inf}\le z\left(t\right)\le z_{sup}\right\},\hfill \end{array}$$

and expressions for the lower and upper bounds of CAR-T cells and CLL cells populations are defined with respect to therapy administration cases as shown below:

Chemotherapy off	Chemotherapy on
$y_{inf}={\displaystyle \frac{\psi }{{\psi }_1}}$	$y_{inf}={\displaystyle \frac{\psi }{{\psi }_1+{\psi }_3}}$
$y_{sup}={\displaystyle \frac{\psi }{{\psi }_1}}+{\displaystyle \frac{{\kappa }_2{\left({\kappa }_1+{\psi }_1\right)}^2}{4{\psi }_1{\kappa }_1}}$	$y_{sup}={\displaystyle \frac{\psi }{\alpha }}+{\displaystyle \frac{{\kappa }_2{\left({\kappa }_1+{\psi }_1\right)}^2}{4\alpha {\kappa }_1}}$ where $\alpha ={\psi }_1+{\displaystyle \frac{{\psi }_3{z}_{inf}}{{\psi }_4+z_{inf}}}$

Immunotherapy	Chemotherapy	Chemoimmunotherapy
$x_{sup}={\kappa }_2\left(1-{\displaystyle \frac{{\kappa }_3\psi }{{\kappa }_1{\psi }_1}}\right)$	$x_{sup}={\kappa }_2\left(1-{\displaystyle \frac{{\kappa }_4{z}_{inf}}{{\kappa }_1\left({\kappa }_5+z_{inf}\right)}}\right)$	$x_{sup}={\kappa }_2\left(1-{\displaystyle \frac{{\kappa }_3\psi }{{\kappa }_1\left({\psi }_1+{\psi }_3\right)}}-{\displaystyle \frac{{\kappa }_4{z}_{inf}}{{\kappa }_1\left({\kappa }_5+z_{inf}\right)}}\right)$

Proof.  

Now, in order to demonstrate Result I , we first compute the bounds concerning all forward solutions $z\left(t\right)$ by integrating Equation (3) with the initial condition $z\left(0\right)=z_0\in {\mathbf{R}}_{+,0}$ as follows:

$$z\left(t\right)={\displaystyle \frac{ϵ}{{ϵ}_1}}-\left({\displaystyle \frac{ϵ}{{ϵ}_1}}-z_0\right)e^{-{ϵ}_{1}t},$$

where it is evident that

$$\underset{t\to \infty }{lim}z\left(t\right)={\displaystyle \frac{ϵ}{{ϵ}_1}},$$

which also represents the equilibrium point for the constant infusion of the drug into the system. Now, bounds for $z\left(t\right)$ may be defined as shown below

$$K_z=\left\{z_{inf}\le z\left(t\right)\le z_{sup}\right\},$$

from which different cases should be considered regarding the initial condition:

Case 1:$z_0=0$.

In this case chemotherapy is not considered at the beginning of the treatment, therefore:

$$z_{inf}=0\mathrm{and}{z}_{sup}={\displaystyle \frac{ϵ}{{ϵ}_1}}.$$

Case 2:$z_0<ϵ/{ϵ}_1$.

In this case, bounds are given as follows:

$$z_{inf}=z_0\mathrm{and}{z}_{sup}={\displaystyle \frac{ϵ}{{ϵ}_1}}.$$

Case 3: $z_0>ϵ/{ϵ}_1$.

In this case, bounds are defined as follows:

$$z_{inf}={\displaystyle \frac{ϵ}{{ϵ}_1}}\mathrm{and}{z}_{sup}=z_0.$$

Case 4: $z_0=ϵ$.

In this case, the initial condition is given by the chemotherapy dose. Thus, when considering parameter values, bounds are given by:

$$z_{inf}={\displaystyle \frac{ϵ}{{ϵ}_1}}\mathrm{and}{z}_{sup}=ϵ.$$

Here, it is important to note that in all numerical simulations, only Case 1 arises, as chemotherapy is not administered as the initial treatment in the designed protocols. These bounds will be necessary to compute lower and upper bounds for both cell’s populations when applying the Iterative Theorem from the LCIS method. First, let us explore the next localizing function:

$$h_1=y,$$

by computing its Lie derivative

$$L_f{h}_1=\psi -{\psi }_{1}y+{\psi }_{2}xy-{\displaystyle \frac{{\psi }_{3}yz}{{\psi }_4+z}},$$

and constructing the set $S\left(h_1\right)=\left\{L_f{h}_1=0\right\}$

$$S\left(h_1\right)=\left\{\psi -\left({\psi }_1+{\psi }_3\right)y+{\psi }_{2}xy+{\displaystyle \frac{{\psi }_3{\psi }_{4}y}{{\psi }_4+z}}=0\right\},$$

which is rewritten as follows

$$S\left(h_1\right)=\left\{y={\displaystyle \frac{\psi }{{\psi }_1+{\psi }_3}}+{\displaystyle \frac{{\psi }_{2}xy}{{\psi }_1+{\psi }_3}}+{\displaystyle \frac{{\psi }_3{\psi }_{4}y}{\left({\psi }_4+z\right)\left({\psi }_1+{\psi }_3\right)}}\right\},$$

from the latter, the lower bound for the CAR-T cells forward solution $\left[y\left(t\right)\right]$ is given by

$$K\left(h_1\right)=\left\{y\left(t\right)\ge y_{inf}\right\},$$

where the next two cases should be differentiated

$$\begin{array}{ccc}\hfill Case1\left[\mathrm{Chemotherapy}\mathrm{off}\right]& :& \mathrm{If}z\left(t\right)=0\Rightarrow y_{inf}={\displaystyle \frac{\psi }{{\psi }_1}},\hfill \\ \hfill Case2\left[\mathrm{Chemotherapy}\mathrm{on}\right]& :& \mathrm{If}z\left(t\right)>0\Rightarrow y_{inf}={\displaystyle \frac{\psi }{{\psi }_1+{\psi }_3}}.\hfill \end{array}$$

Hence, it is evident that

$${\displaystyle \frac{\psi }{{\psi }_1+{\psi }_3}}<{\displaystyle \frac{\psi }{{\psi }_1}},$$

as it is to be expected because the chemotherapy drug is cytotoxic to both cancer cells and CAR-T cells. Now, an upper bound for the CAR-T cells population may be estimated from the next localizing function:

$$h_2=x+y,$$

whose Lie derivative is calculated as

$$L_f{h}_2={\kappa }_{1}x\left(1-{\displaystyle \frac{x}{{\kappa }_2}}\right)-{\kappa }_{3}xy-{\displaystyle \frac{{\kappa }_{4}xz}{{\kappa }_5+z}}+\psi -{\psi }_{1}y+{\psi }_{2}xy-{\displaystyle \frac{{\psi }_{3}yz}{{\psi }_4+z}},$$

then, set $S\left(h_2\right)=\left\{L_f{h}_2=0\right\}$ is given below

$$S\left(h_2\right)=\left\{\psi +{\kappa }_{1}x-{\displaystyle \frac{{\kappa }_1}{{\kappa }_2}}{x}^2-\left({\psi }_1+{\displaystyle \frac{{\psi }_{3}z}{{\psi }_4+z}}\right)y-\left({\kappa }_3-{\psi }_2\right)xy-{\displaystyle \frac{{\kappa }_{4}xz}{{\kappa }_5+z}}=0\right\},$$

from the latter, the following condition must be satisfied, implying that the cytotoxic rate of CAR-T cells against cancer cells exceeds the inactivation rate of cancer cells on CAR-T cells,

$${\kappa }_3>{\psi }_2,$$

(9)

and set $S\left(h_2\right)$ is rewritten as follows

$$S\left(h_2\right)\cap \left\{z=z_{inf}\right\}\subset \left\{{\alpha }_1{h}_2\le \psi +{\displaystyle \frac{{\kappa }_2{\left({\kappa }_1+{\psi }_1\right)}^2}{4{\kappa }_1}}-f_1\left(x,y,z\right)\right\},$$

with

$$\alpha ={\psi }_1+{\displaystyle \frac{{\psi }_3{z}_{inf}}{{\psi }_4+z_{inf}}}>0,$$

(10)

and

$$f_1\left(x,y,z\right)={\displaystyle \frac{{\kappa }_1}{{\kappa }_2}}{\left(x-{\displaystyle \frac{{\kappa }_2\left({\kappa }_1+{\psi }_1\right)}{2{\kappa }_1}}\right)}^2+\left({\kappa }_3-{\psi }_2\right)xy+{\displaystyle \frac{{\kappa }_{4}xz}{{\kappa }_5+z}}.$$

Therefore, by neglecting the negative function in the right-hand side of the inequality, the next result is obtained

$$K\left(h_2\right)=\left\{x\left(t\right)+y\left(t\right)\le {\displaystyle \frac{\psi }{\alpha }}+{\displaystyle \frac{{\kappa }_2{\left({\kappa }_1+{\psi }_1\right)}^2}{4\alpha {\kappa }_1}}\right\},$$

hence, the following upper bound for the CAR-T cells can be extracted

$$K_y\left(h_2\right)=\left\{y\left(t\right)\le y_{sup}={\displaystyle \frac{\psi }{\alpha }}+{\displaystyle \frac{{\kappa }_2{\left({\kappa }_1+{\psi }_1\right)}^2}{4\alpha {\kappa }_1}}\right\},$$

where parameter ${\alpha }_1$ can have one of the following two values

$$\begin{array}{ccc}\hfill Case1\left[\mathrm{Chemotherapy}\mathrm{off}\right]& :& \mathrm{If}z\left(t\right)=0\Rightarrow \alpha ={\psi }_1,\hfill \\ \hfill Case2\left[\mathrm{Chemotherapy}\mathrm{on}\right]& :& \mathrm{If}z\left(t\right)>0\Rightarrow \alpha ={\psi }_1+{\displaystyle \frac{{\psi }_3{z}_{inf}}{{\psi }_4+z_{inf}}}.\hfill \end{array}$$

Now, let us exploit the localizing function:

$$h_3=x,$$

to derive results concerning the ultimate upper bound for the CAR-T cells in the case of single therapy and combined therapy. Thus, the Lie derivative is computed as follows

$$L_f{h}_3={\kappa }_{1}x\left(1-{\displaystyle \frac{x}{{\kappa }_2}}\right)-{\kappa }_{3}xy-{\displaystyle \frac{{\kappa }_{4}xz}{{\kappa }_5+z}},$$

then, set $S\left(h_3\right)=\left\{L_f{h}_3=0\right\}$ is written as given below

$$S\left(h_3\right)=\left\{{\displaystyle \frac{{\kappa }_1}{{\kappa }_2}}x={\kappa }_1-{\kappa }_{3}y-{\displaystyle \frac{{\kappa }_{4}z}{{\kappa }_5+z}}\right\}\cup \left\{x=0\right\}.$$

Thus, the Iterative Theorem can be applied to determine the following cases concerning single or combined treatment strategies.

Case 1: Immunotherapy. 

The Iterative Theorem of the LCIS method is applied as follows:

$$S\left(h_3\right)\cap \left\{y=y_{inf}\right\}\subset \left\{x\le {\kappa }_2-{\displaystyle \frac{{\kappa }_2}{{\kappa }_1}}{\kappa }_3{y}_{inf}\right\},$$

and by substituting the corresponding value of Case 1 $\left[\mathrm{Chemotherapy}\mathrm{off}\right]$ for the CAR-T cells lower bound one can get the next result:

$$K_1\left(h_3\right)=\left\{x\left(t\right)\le x_{sup}={\kappa }_2\left(1-{\displaystyle \frac{{\kappa }_3\psi }{{\kappa }_1{\psi }_1}}\right)\right\}.$$

Case 2: Chemotherapy. 

When only the chlorambucil drug is being considered, the Iterative Theorem is applied as shown below:

$$S\left(h_3\right)\cap \left\{z=z_{inf}\right\}\subset \left\{x\le {\kappa }_2-{\displaystyle \frac{{\kappa }_2}{{\kappa }_1}}{\displaystyle \frac{{\kappa }_4{z}_{inf}}{{\kappa }_5+z_{inf}}}\right\},$$

and the lower bound is defined as follows

$$K_2\left(h_3\right)=\left\{x\left(t\right)\le x_{sup}={\kappa }_2\left(1-{\displaystyle \frac{{\kappa }_4{z}_{inf}}{{\kappa }_1\left({\kappa }_5+z_{inf}\right)}}\right)\right\}.$$

Case 3: Chemoimmunotherapy. 

Now, let us explore the case when the combined therapy is applied, thus, the Iterative Theorem yields the next result:

$$S\left(h_3\right)\cap \left\{y=y_{inf}\right\}\cap \left\{z=z_{inf}\right\}\subset \left\{x\le {\kappa }_2-{\displaystyle \frac{{\kappa }_2}{{\kappa }_1}}{\kappa }_3{y}_{inf}-{\displaystyle \frac{{\kappa }_2}{{\kappa }_1}}{\displaystyle \frac{{\kappa }_4{z}_{inf}}{{\kappa }_5+z_{inf}}}\right\},$$

then, when substituting the corresponding value of Case 2 $\left[\mathrm{Chemotherapy}\mathrm{on}\right]$ for the CAR-T cells lower bound the next result can be establish:

$$K_3\left(h_3\right)=\left\{x\left(t\right)\le x_{sup}={\kappa }_2\left(1-{\displaystyle \frac{{\kappa }_3\psi }{{\kappa }_1\left({\psi }_1+{\psi }_3\right)}}-{\displaystyle \frac{{\kappa }_4{z}_{inf}}{{\kappa }_1\left({\kappa }_5+z_{inf}\right)}}\right)\right\}.$$

Therefore, the proof for Result I is now complete. □

3.2. Eradication Conditions

In this section we establish sufficient conditions by means of Lyapunov’s Direct method to ensure the complete eradication of the CLL cells population described by system (1)–(3). These conditions are formulated on the treatment parameters when considering single and combined application of immunotherapy and chemotherapy. Hence, conditions for CLL cells eradication are given in the next result:

Result II: Minimum concentrations of therapies to ensure a continuous decrease in tumor burden.  If only CAR-T cells are applied, the concentration of this therapy must satisfy the following condition:

$$\psi >{\displaystyle \frac{{\kappa }_1{\psi }_1}{{\kappa }_3}}.$$

On the other hand, if only the chlorambucil drug is applied$\left(z_0=ϵ\right)$and the inequality${\kappa }_4>{\kappa }_1$holds, then the chemotherapy concentration must fulfill the following condition:

$$ϵ>{\displaystyle \frac{{ϵ}_1{\kappa }_1{\kappa }_5}{{\kappa }_4-{\kappa }_1}}.$$

When combined therapy is applied and the concentration of chlorambucil is fixed, the amount of CAR-T cells that needs to be administered into the CLL system (1)–(3) to ensure leukemia cells eradication is given by:

$$\psi >{\displaystyle \frac{\left[{\kappa }_1{\kappa }_5+\left({\kappa }_1-{\kappa }_4\right)z_{inf}\right]\left({\psi }_1+{\psi }_3\right)}{{\kappa }_3\left({\kappa }_5+z_{inf}\right)}}.$$

Proof.  

Let us apply the next Lyapunov’s candidate function

$$V\left(x\right)=x,$$

where it is necessary to highlight that we can explore this function because the lower and upper bounds have already been calculated for both CAR-T cells and the chlorambucil drug in the localizing domain $K_{xyz}$ shown in Result I . Thus, inequalities can be constructed on the time-derivative of $V\left(x\right)$, which is given as follows:

$$\dot{V}\left(x\right)={\kappa }_{1}x\left(1-{\displaystyle \frac{x}{{\kappa }_2}}\right)-{\kappa }_{3}xy-{\displaystyle \frac{{\kappa }_{4}xz}{{\kappa }_5+z}}.$$

Therefore, the following three cases can be formulated based on the application of a single therapy or combined therapy.

Case 1: Immunotherapy. 

When only the CAR-T cells are applied for cancer treatment the time-derivative can be bounded from above as follows:

$$\dot{V}\left(x\right)={\kappa }_{1}x\left(1-{\displaystyle \frac{x}{{\kappa }_2}}\right)-{\kappa }_{3}xy-{\displaystyle \frac{{\kappa }_{4}xz}{{\kappa }_5+z}}<\left({\kappa }_1-{\displaystyle \frac{{\kappa }_1}{{\kappa }_2}}x-{\kappa }_3{y}_{inf}\right)x<0,$$

where $y_{inf}=\psi /{\psi }_1$, and the next constraint can be formulated on the immunotherapy concentration:

$$\psi >{\displaystyle \frac{{\kappa }_1{\psi }_1}{{\kappa }_3}}.$$

(11)

Case 2: Chemotherapy. 

When only the chlorambucil drug is being considered from the beginning of the treatment the time-derivative is bounded from above as indicated below:

$$\dot{V}\left(x\right)={\kappa }_{1}x\left(1-{\displaystyle \frac{x}{{\kappa }_2}}\right)-{\kappa }_{3}xy-{\displaystyle \frac{{\kappa }_{4}xz}{{\kappa }_5+z}}<\left({\kappa }_1-{\displaystyle \frac{{\kappa }_1}{{\kappa }_2}}x-{\displaystyle \frac{{\kappa }_4{z}_{inf}}{{\kappa }_5+z_{inf}}}\right)x<0,$$

where $z_{inf}=ϵ/{ϵ}_1$ as $z_0=ϵ$, and the next constraint is formulated on the chemotherapy concentration

$$ϵ>{\displaystyle \frac{{ϵ}_1{\kappa }_1{\kappa }_5}{{\kappa }_4-{\kappa }_1}},$$

(12)

however, it is evident from the denominator that the next condition relating the rates of CLL cells growth and cytotoxicity of the drug should be fulfilled:

$${\kappa }_4-{\kappa }_1>0,$$

(13)

or else, condition (12) will not have a biologically feasible solution.

Case 3: Chemoimmunotherapy. 

When the combined therapy is applied the time-derivative is bounded as follows:

$$\dot{V}\left(x\right)={\kappa }_{1}x\left(1-{\displaystyle \frac{x}{{\kappa }_2}}\right)-{\kappa }_{3}xy-{\displaystyle \frac{{\kappa }_{4}xz}{{\kappa }_5+z}}<\left({\kappa }_1-{\displaystyle \frac{{\kappa }_1}{{\kappa }_2}}x-{\kappa }_3{y}_{inf}-{\displaystyle \frac{{\kappa }_4{z}_{inf}}{{\kappa }_5+z_{inf}}}\right)x<0,$$

where $y_{inf}=\psi /\left({\psi }_1+{\psi }_3\right)$. In order to properly solve this condition, we set the CAR-T cells as the control therapy and come to the following result:

$$\psi >{\displaystyle \frac{\left[{\kappa }_1{\kappa }_5+\left({\kappa }_1-{\kappa }_4\right)z_{inf}\right]\left({\psi }_1+{\psi }_3\right)}{{\kappa }_3\left({\kappa }_5+z_{inf}\right)}},$$

(14)

in this case it is important to consider two paths. If condition (13) is satisfied, chemotherapy alone could eradicate the leukemia cells population described by the CLL mechanistic model (1)–(3). However, both therapies may be combined in varying proportions to design a successful protocol that ultimately achieves cancer eradication. Conversely, if condition (13) is not fulfilled, then a biologically feasible solution can be computed to estimate the concentration of CAR-T cells required to ensure the complete eradication of cancer cells in the corresponding scenarios.

Therefore, the proof for Result II is now complete. □

3.3. Persistence Conditions

In order to derive results regarding the so-called persistence conditions for both CLL cells and CAR-T cells, we compute the core equilibrium points of system (1)–(3) by excluding the therapy parameters and setting the equations to zero as follows:

$$\begin{array}{ccc}\hfill {\kappa }_{1}x-{\displaystyle \frac{{\kappa }_1{x}^2}{{\kappa }_2}}-{\kappa }_{3}xy-{\displaystyle \frac{{\kappa }_{4}xz}{{\kappa }_5+z}}& =& 0,\hfill \\ \hfill -{\psi }_{1}y+{\psi }_{2}xy-{\displaystyle \frac{{\psi }_{3}yz}{{\psi }_4+z}}& =& 0,\hfill \\ \hfill -{ϵ}_{1}z& =& 0.\hfill \end{array}$$

It should be noted that treatments are not considered because, in real-life clinical scenarios, anti-cancer therapies are administered in short cycles over limited periods of time. The latter, within the scope of nonlinear system theory, are referred to as perturbations and their primary characteristic is their ability to alter both the short- and long-term dynamics of a system. These changes may be minor, disrupting only short-term behavior, or major, affecting long-term dynamics and transitioning the system through completely different states. Hence, by performing the corresponding arithmetical operations, the following three sets of equilibrium points are identified:

$$\begin{array}{ccc}\hfill \left(x_0^{*},y_0^{*},z_0^{*}\right)& =& \left(0,0,0\right),\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill \left(x_1^{*},y_1^{*},z_0^{*}\right)& =& \left({\displaystyle \frac{{\psi }_1}{{\psi }_2}},{\displaystyle \frac{{\kappa }_1\left({\kappa }_2{\psi }_2-{\psi }_1\right)}{{\kappa }_2{\kappa }_3{\psi }_2}},0\right),\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill \left(x_2^{*},y_0^{*},z_0^{*}\right)& =& \left({\kappa }_2,0,0\right).\hfill \end{array}$$

(17)

Now, let us discuss the biological interpretation of each equilibrium point. Equilibrium (15) corresponds to the tumor-free state, representing the desired outcome after a successful treatment strategy that eradicates the CLL cells population. Concerning equilibrium (16), we can see that it is biologically feasible if the persistence condition (4) from Remark 3 is satisfied. Furthermore, positive values are required to describe real-life cells concentrations. This equilibrium represents a state where CAR-T cells persist in the long-term, enabling them to mount a consistent immune response against cancer cells, thereby reducing the tumor burden in the system to a value directly related to the ratio ${\psi }_1/{\psi }_2$. Equilibrium (17) represents a scenario in which leukemia cells reach their maximum carrying capacity. Within the framework of system (1)–(3) this implies the death of the subject, as the tumor burden becomes unsustainable for survival. Now, let us present the following result regarding the long-term concentration of cells populations:

Result III: Persistence and depletion of CAR-T cells.  If condition (4) from Remark III is satisfied, i.e.,

$${\psi }_2>{\displaystyle \frac{{\psi }_1}{{\kappa }_2}},$$

CAR-T cells will persist in the long-term after a single application of the therapy. Under this condition, cell populations will reach the following steady-state concentrations where both coexist:

$$\underset{t\to \infty }{lim}x\left(t\right)={\displaystyle \frac{{\psi }_1}{{\psi }_2}},\underset{t\to \infty }{lim}y\left(t\right)={\displaystyle \frac{{\kappa }_1\left({\kappa }_2{\psi }_2-{\psi }_1\right)}{{\kappa }_2{\kappa }_3{\psi }_2}}.$$

However, if condition (4) is not satisfied, CAR-T cells will eventually be depleted, and CLL cells will grow to their maximum carrying capacity, i.e.,

$$\underset{t\to \infty }{lim}x\left(t\right)={\kappa }_2,\underset{t\to \infty }{lim}y\left(t\right)=0.$$

Proof.  

Local dynamics of the CLL mechanistic model can be explored by applying the Indirect Lyapunov’s method which requires the Jacobian matrix of the system which is computed as follows:

$$J=\left[\begin{array}{ccc}{\kappa }_1-{\displaystyle \frac{2{\kappa }_{1}x}{{\kappa }_2}}-{\kappa }_{3}y-{\displaystyle \frac{{\kappa }_{4}z}{{\kappa }_5+z}}& -{\kappa }_{3}x& -{\displaystyle \frac{{\kappa }_4{\kappa }_{5}x}{{\left({\kappa }_5+z\right)}^2}}\\ {\psi }_{2}y& -{\psi }_1+{\psi }_{2}x-{\displaystyle \frac{{\psi }_{3}z}{{\psi }_4+z}}& -{\displaystyle \frac{{\psi }_3{\psi }_{4}y}{{\left({\psi }_4+z\right)}^2}}\\ 0& 0& -{ϵ}_1\end{array}\right].$$

(18)

From the latter, one can determine the corresponding eigenvalues to each equilibrium point as indicated in the cases below.

Case 1: Equilibrium$\left(x_0^{*},y_0^{*},z_0^{*}\right).$ 

The Jacobian matrix has the following form:

$${\left.J\right|}_{\left(x_0^{*},y_0^{*},z_0^{*}\right)}=\left[\begin{array}{ccc}{\kappa }_1& 0& 0\\ 0& -{\psi }_1& 0\\ 0& 0& -{ϵ}_1\end{array}\right].$$

Hence, it is evident that this equilibrium is unstable, as one of the eigenvalues is positive,

$$\begin{array}{ccc}{\lambda }_1={\kappa }_1,& {\lambda }_2=-{\psi }_1,& {\lambda }_3=-{ϵ}_1.\end{array}$$

This outcome is to be expected, as once cancer cells evade the immune response, either by suppressing the immune system or by remaining undetected by effector cells, they will grow to their maximum carrying capacity or to the extent that the subject can endure. This scenario also indicates that treatment strategies were either unsuccessful or not administered.

Case 2: Equilibrium$\left(x_1^{*},y_1^{*},z_0^{*}\right).$ 

The Jacobian matrix evaluated at this equilibrium point is given by:

$${\left.J\right|}_{\left(x_1^{*},y_1^{*},z_0^{*}\right)}=\left[\begin{array}{ccc}-{\kappa }_1{\displaystyle \frac{{\psi }_1}{{\kappa }_2{\psi }_2}}& -{\kappa }_3{\displaystyle \frac{{\psi }_1}{{\psi }_2}}& -{\displaystyle \frac{{\kappa }_4}{{\kappa }_5}}{\displaystyle \frac{{\psi }_1}{{\psi }_2}}\\ -{\displaystyle \frac{{\kappa }_1}{{\kappa }_2{\kappa }_3}}\left({\psi }_1-{\kappa }_2{\psi }_2\right)& 0& {\kappa }_1{\psi }_3{\displaystyle \frac{{\psi }_1-{\kappa }_2{\psi }_2}{{\kappa }_2{\kappa }_3{\psi }_2{\psi }_4}}\\ 0& 0& -{ϵ}_1\end{array}\right],$$

with the next eigenvalues:

$$\begin{array}{ccc}\hfill {\lambda }_1& =& -{\displaystyle \frac{1}{2{\kappa }_2{\psi }_2}}\left({\kappa }_1{\psi }_1-\sqrt{{\kappa }_1^2{\psi }_1^2-4{\kappa }_1{\kappa }_2{\psi }_1{\psi }_2\left({\kappa }_2{\psi }_2-{\psi }_1\right)}\right),\hfill \\ \hfill {\lambda }_2& =& -{\displaystyle \frac{1}{2{\kappa }_2{\psi }_2}}\left({\kappa }_1{\psi }_1+\sqrt{{\kappa }_1^2{\psi }_1^2-4{\kappa }_1{\kappa }_2{\psi }_1{\psi }_2\left({\kappa }_2{\psi }_2-{\psi }_1\right)}\right),\hfill \\ \hfill {\lambda }_3& =& -{ϵ}_1.\hfill \end{array}$$

It is evident that the radicands in ${\lambda }_1$ and ${\lambda }_2$ will be positive and real if the persistence condition (4) from Remark III holds, i.e., ${\psi }_2>{\psi }_1/{\kappa }_2$,which allows us to conclude that equilibrium (16) is locally asymptotically stable. Conversely, if condition (4) is not satisfied, equilibrium (16) becomes unstable and biologically infeasible, as $y_1^{*}<0$. In this case, if therapies are unsuccessful, cancer cells will grow to their maximum carrying capacity ${\kappa }_2$, making (17) the only locally asymptotically stable equilibrium point of the system as it is shown in the next case.

Case 3: Equilibrium$\left(x_2^{*},y_0^{*},z_0^{*}\right).$ 

The Jacobian matrix is upper triangular, as shown below:

$${\left.J\right|}_{\left(x_2^{*},y_0^{*},z_0^{*}\right)}=\left[\begin{array}{ccc}-{\kappa }_1& -{\kappa }_2{\kappa }_3& -{\displaystyle \frac{{\kappa }_2{\kappa }_4}{{\kappa }_5}}\\ 0& -\left({\psi }_1-{\kappa }_2{\psi }_2\right)& 0\\ 0& 0& -{ϵ}_1\end{array}\right].$$

Therefore, the eigenvalues correspond to the diagonal elements of the matrix:

$$\begin{array}{ccc}{\lambda }_1=-{\kappa }_1,& {\lambda }_2=-\left({\psi }_1-{\kappa }_2{\psi }_2\right),& {\lambda }_3=-{ϵ}_1.\end{array}$$

From these results, two scenarios are identified. First, equilibrium point (17) is locally asymptotically stable if

$${\psi }_2<{\displaystyle \frac{{\psi }_1}{{\kappa }_2}}.$$

This indicates that the persistence condition (4) from Remark III is not satisfied. In this case, long-term persistence of CAR-T cells cannot be expected, as equilibrium (16) does not exist within the positive and biologically feasible domain of the system. Nonetheless, if the persistence condition (4) is satisfied, equilibrium (17) becomes unstable, and the only locally asymptotically stable equilibrium point is (16). This implies that both cell populations described by the CLL mechanistic model (1)–(3) coexist. However, the in silico experimentation will reveal whether this sustained immune response by the CAR-T cells is sufficient to reduce the tumor burden to a level that is manageable for the subject’s health. Results from these three cases allows us to conclude that if the chemoimmunotherapy treatment strategy is not successful, leukemia cells will always persist.

Therefore, the proof for Result III is now complete. □

3.4. Existence and Uniqueness

A system of first-order nonlinear ODEs may exhibit multiple solutions for a given set of initial conditions. This ambiguity poses significant challenges when modeling specific scenarios in tumor-immune dynamics, as these dynamics are intrinsically linked to particular conditions and parameter values. However, based on the results presented in this section, we establish the following conclusions regarding the existence and uniqueness of solutions for the CLL system (1)–(3) proposed in this work:

Result IV: Ensuring a reliablein silicoexperimentation for tumor-immune dynamics. Let vector fields$f_i\left(t,\rho \right)$, where$\rho ={\left[x,y,z\right]}^T$, be piecewise continuous intand locally Lipchitz in$\rho$for all$t\ge t_0$and all$\rho$in$K_{xyz}\subset {\mathbf{R}}_{+,0}^3$. Since every solution of the system (1)–(3), with initial conditions$\rho \left(t_0\right)={\rho }_0\in K_{xyz}$, remains entirely within$K_{xyz}$, there exist a unique solution defined for all$t\ge t_0$ .

Proof.  

First, note that vector fields $f_i\left(t,\rho \right)$ of system (1)–(3) are given by

$$\begin{array}{ccc}\hfill f_1\left(t,\rho \right)& =& {\kappa }_{1}x\left(1-{\displaystyle \frac{x}{{\kappa }_2}}\right)-{\kappa }_{3}xy-{\displaystyle \frac{{\kappa }_{4}xz}{{\kappa }_5+z}},\hfill \\ \hfill f_2\left(t,\rho \right)& =& \psi -{\psi }_{1}y+{\psi }_{2}xy-{\displaystyle \frac{{\psi }_{3}yz}{{\psi }_4+z}},\hfill \\ \hfill f_3\left(t,\rho \right)& =& ϵ-{ϵ}_{1}z,\hfill \end{array}$$

where all elements are continuous in $\rho ={\left[x,y,z\right]}^T$ by Remark I , i.e., ${\mathbf{R}}_{+,0}^3=\left\{x\left(t\right),y\left(t\right),z\left(t\right)\ge 0\right\}$, and piecewise continuous in t. This is crucial because therapies are considered perturbations, applied as impulse changes in short cycles over limited periods of time. Furthermore, one can see that all components of the Jacobian Matrix $\left[\partial f/\partial \rho \right]\left(t,x\right)$ (18) are also piecewise continuous in t and bounded within the domain $K_{xyz}\subset {\mathbf{R}}_{+,0}^3$, as established by Results I and II . Consequently, $f_i\left(t,\rho \right)$ is locally Lipchitz in $\rho$ in the domain $\mathsf{\Omega }=$$\left[t,t_1\right]\times K_{xyz}$, where $\left[t_0,t_1\right]\subset \left[t_0,\infty \right)$, as there exists a Lipschitz constant $\ell \ge 0$ such that

$$∥{\displaystyle \frac{\partial f}{\partial \rho }}\left(t,x\right)∥\le \ell$$

on $\mathsf{\Omega }$, and $f_i\left(t,\rho \right)$ satisfies the Lipchitz condition:

$$∥f\left(t,{\rho }_1\right)-f\left(t,{\rho }_2\right)∥\le \ell ∥{\rho }_1-{\rho }_2∥.$$

Therefore, the proof for Result IV is now complete. □

4. Discussion: In Silico Experimentation

In this section, we explore the global dynamics of our mechanistic model (1)–(3) as well as its capabilities to investigate different scenarios of tumor-immune dynamics through in silico experimentation. First, Figure 2 illustrates outcomes related to condition (4) from Result III , and the eradication condition (11) from Result II through the following three cases:

Case 1: Depletion of CAR-T cells. 

Depletion of CAR-T cells in the short-term happens when condition (4) is not satisfied. As expected, the concentration of CLL cells $\left[x\left(t\right)\right]$ grows to their maximum carrying capacity, while CAR-T cells $\left[y\left(t\right)\right]$ eventually fall below the threshold of clinically significant biological behavior, defined as fewer than 1 cell. Hence, solutions go to the equilibrium point (17), i.e, $\left(x_2^{*},y_0^{*},z_0^{*}\right)=\left({\kappa }_2,0,0\right)$. For the sake of numerical simulations, dynamics are simulated following a single application of the therapy under the next set of initial conditions: $x\left(0\right)=1\times {10}^{10}$ CLL cells, $y\left(0\right)=$$1\times {10}^7$ CAR-T cells (unique dose), and $z\left(0\right)=0$ mg of chlorambucil. Parameter values are as follows: ${\kappa }_1=1.680$, ${\kappa }_2=3\times {10}^{11}$, ${\kappa }_3=1.256\times {10}^{-7}$, ${\psi }_1=2.236\times {10}^{-2}$, ${\psi }_2=5.217\times {10}^{-14}$. Other parameters are set as zero, as chemotherapy is not considered in this scenario. Results for this in silico experimentation are shown in the top panels of Figure 2.

Case 2: Persistence of CAR-T cells. 

Long-term persistence of CAR-T cells $\left[y\left(t\right)\right]$ is observed when condition (4) is fulfilled. In this case, cells concentrations stabilize at a steady-state defined by the equilibrium point (16) $\left(x_1^{*},y_1^{*},z_0^{*}\right)=\left({\displaystyle \frac{{\psi }_1}{{\psi }_2}},{\displaystyle \frac{{\kappa }_1\left({\kappa }_2{\psi }_2-{\psi }_1\right)}{{\kappa }_2{\kappa }_3{\psi }_2}},0\right)$. However, this outcome does not represent a viable treatment strategy for real-world clinical scenarios, because the tumor burden $\left[x\left(t\right)\right]$ remains near the maximum carrying capacity. For the sake of numerical simulations, dynamics are simulated following a single application of the therapy under the next set of initial conditions: $x\left(0\right)=1\times {10}^{10}$ CLL cells, $y\left(0\right)=$$1\times {10}^7$ CAR-T cells (unique dose), and $z\left(0\right)=0$ mg of chlorambucil. Parameter values are as follows: ${\kappa }_1=1.680$, ${\kappa }_2=3\times {10}^{11}$, ${\kappa }_3=1.256\times {10}^{-7}$, ${\psi }_1=2.236\times {10}^{-2}$, ${\psi }_2=8.199\times {10}^{-14}$. Other parameters are set as zero, as chemotherapy is not considered in this scenario. Results for this in silico experimentation are shown in the middle panels of Figure 2.

Case 3: Eradication of CLL cells. 

Within the scope of system (1)–(3), eradication of CLL cells is ensured when the concentration of CAR-T cells meets condition (11). Hence, the immunotherapy control parameter $\psi$ is set with a constant value of $1.1\times {\kappa }_1{\psi }_1/{\kappa }_3$ which equals to $328,980$ cells. Dynamics are simulated following a constant application of the therapy under the next set of initial conditions: $x\left(0\right)=1\times {10}^{10}$ CLL cells, $y\left(0\right)=328,980$ CAR-T cells (constant dose), and $z\left(0\right)=0$ mg of chlorambucil. Parameter values are as follows: ${\kappa }_1=1.680$, ${\kappa }_2=3\times {10}^{11}$, ${\kappa }_3=1.256\times {10}^{-7}$, ${\psi }_1=2.236\times {10}^{-2}$, ${\psi }_2=8.199\times {10}^{-14}$. Other parameters are set as zero, as chemotherapy is not considered in this scenario. Results for this in silico experimentation are shown in the bottom panels of Figure 2.

All chemoimmunotherapy protocols explored in this work are designed under the next assumptions: CAR-T cells will persist in the long-term after infusion into the subject as it was described in Case 2 of Figure 2, i.e, ${\psi }_2=8.199\times {10}^{-14}$${\left(cells\times days\right)}^{-1}$; CLL cells display the highest growth rate, i.e, ${\kappa }_1=1.680$ $day{s}^{-1}$ and are resistant to the chemotherapy drug, i.e., ${\kappa }_5=92.785$$mg$.

Hence, let us explore the dose-escalation scheme often employed in phase I clinical trials for cancer patients treated with T-cell therapies [37,38,39,40,41]. In such protocols, dose levels are typically defined as an absolute number of cells, ranging from $2\times {10}^7$ to $8\times {10}^7$ cells, or as dose normalized to the patient’s body weight, which may range from $5\times {10}^4$$cells/kg$ up to a maximum dose of $1\times {10}^6$$cells/kg$. It is important to note that dosing regimens can vary between studies, with the optimal dose typically determined through careful escalation in clinical trials to balance efficacy and safety. However, this falls outside the scope of the present research, as the focus is on in silico experimentation rather than in vivo studies. Hence, after extensive iterations of numerical simulations, the scheme of Table 3 was designed using an algorithm developed in MATLAB.

This protocol aims to eliminate CLL cells concentrations under varying initial tumor burden scenarios, defined as $x\left(0\right)={10}^i$ leukemia cells, for $i=1,2,...,11$. For the purposes of the in silicoexperimentation, the maximum initial tumor burden considered represents approximately $33\%$ of the system’s estimated carrying capacity. Results are illustrated in Figure 3, the top panel corresponds to the CLL cells $\left[x\left(t\right)\right]$, demonstrating the combined effects of the CAR-T cells $\left[y\left(t\right)\right]$ and chlorambucil $\left[z\left(t\right)\right]$ doses over time. It is evident that, under the same protocol scheme, a lower initial concentration of leukemia cells is eradicated in a shorter period of time. Panels in the middle row depict the behavior of CAR-T cells and the chlorambucil drug in the bloodstream compartment, whereas panels in the bottom row illustrate the dose-escalation protocol. It should be noted that each dose of CAR-T cells satisfies conditions outlined in Result II , particularly condition (11) as the therapies are administered sequentially without overlapping. Chlorambucil doses were calculated based on the median weight of an adult CLL patient, estimated to be 81 $kg$ [42,43], and drug cytotoxicity on both CLL and CAR-T cells changes for each dose concentration as indicated in Table 4.

Now, following nonlinear dynamical systems theory, we can expect that an initial higher dose of CAR-T cells will have a more significant and long-term effect on the overall behavior of the system. Using our algorithm, a dose de-escalation protocol was designed with the characteristics shown in Table 5.

Results are illustrated in Figure 4, where the arrangement of solutions, i.e., cells and drug concentrations dynamics as well as protocols, is consistent with that in Figure 3. Nonetheless, there are some key differences that should be discussed. First, the dose-escalation protocol required a total of $26,916,516$ CAR-T cells to completely eradicate the highest concentration of CLL cells, i.e., $x\left(0\right)={10}^{11}$ cells. In contrast, the dose de-escalation scheme achieved leukemia eradication with a total of $20,097,666$ CAR-T cells. Numerical simulations indicate that in the dose-escalation protocol, CLL cell growth still occurs after the first dose of immunotherapy, as the initial lower concentration of CAR-T cells is insufficient to mount a successful immune response against cancer, even at lower tumor burdens. However, the dose de-escalation protocol demonstrates a consistent reduction in CLL cell concentrations starting from the first dose. Furthermore, the highest dose in the protocol falls within the previously stated range, with ${10}^7$ CAR-T cells administered in the first dose. However, assessing the maximum tolerated dose lies beyond the scope of this in silico experimentation. Additionally, it should be noted that the effects of the chemotherapy drug on CLL cell concentrations are not perceptible in Figure 3 and Figure 4 because numerical simulations are performed under two key assumptions. First, the leukemia cells are assumed to exhibit drug resistance, and second, condition (13) is not satisfied, as the growth rate of leukemia cells is deliberately set to its highest value.

While protocol failure is undesirable in cancer treatment, a significant proportion of phase I clinical trials encounter such outcomes [41,44,45]. Therefore, we illustrate in Figure 5 and Figure 6 tumor escape in both dose-escalation and dose de-escalation schemes, respectively, emphasizing that any protocol may fail if doses are not accurately calculated, and cancer cell concentrations are not consistently monitored. Parameter values are the same as the successful protocols, i.e., ${\kappa }_1=1.680$, ${\kappa }_2=3\times {10}^{11}$, ${\kappa }_3=1.256\times {10}^{-7}$, ${\kappa }_5=92.785$, ${\psi }_1=2.236\times {10}^{-2}$, ${\psi }_2=8.199\times {10}^{-14}$, ${\psi }_4=231.962$, ${ϵ}_1=11.090$. It is important to note that the CAR-T doses were determined heuristically, aiming for a solution $x\left(t\right)$ that approaches the threshold of 1 cancer cell, as established in Remark 2 for clinically significant biological behavior.

In silico experiments presented in this section were performed on a computer with the following specifications: Ryzen 9 7950X CPU, 192 GB DDR5 CL40 RAM, and a 1 TB Samsung 990 Pro Gen 4 NVMe M.2 drive. The algorithm used for these experiments, designed in MATLAB, is available as supplementary material in a GitHub repository, and certified under INDAUTOR registration number 03-2023-052912465900-01.

5. Conclusions

In this work, we proposed a mechanistic nonlinear model to study the dynamics of CLL under combined chemoimmunotherapy. The model consists of a system of three first-order ODEs that describe the effects of CAR-T cells and chlorambucil on the leukemia cells population over both short- and long-term timescales. Using a model-based design approach and extensive in silico experimentation, we determined sufficient dosing strategies for each therapeutic application. The local and global dynamics of model (1)–(3) were thoroughly analyzed using various nonlinear systems theories. Consequently, we derived conditions to ensure the eradication of cancer cells, with numerical simulations illustrating that the population falls below the threshold for clinically significant biological behavior, defined as one cell.

Furthermore, our algorithm enables the exploration of different scenarios. As illustrated, we compared two chemoimmunotherapy application protocols: dose-escalation and dose de-escalation. While dose-escalation protocols are commonly used in Phase I trials, numerical simulations suggest that the dose de-escalation protocol may require a lower total dose to achieve effective tumor-immune dynamics, as a higher initial dose tends to have a greater impact. However, in silico experiments also show that, even with complete remission (peripheral blood lymphocytes count $<4\times {10}^9$$cells/L$) or partial remission (a decrease in the number of blood lymphocytes to $50\%$ or less from the value before therapy) [46], relapse can occur if the dosing is insufficient or if therapy is administered without ensuring the eradication of minimal residual disease.

Nonetheless, for any model-based design approach to be successful in dose-finding for anti-cancer therapies, one open problem remains: all parameters of a mathematical model, such as the one proposed in this work, must be carefully estimated for each patient. Therefore, clinical data on cancer cell growth rates, CAR-T cell activation rates, and the cytotoxicity rates of therapies are crucial values for protocol dosing design from the perspective of nonlinear dynamical systems and in silico experimentation. Finally, our algorithm could be useful for estimating key parameters through nonlinear regression functions in the MATLAB environment when tumor-immune dynamics data, in the form of cell counts over time, are available for a given patient.

Appendix A.1. Positiveness of Solutions

Consider a nonlinear dynamical system of the form (A2), which is called positive, P, if the nonnegative orthant ${\mathbf{R}}_{+,0}^n$ is forward invariant. This implies that for any initial condition $x_0\in {\mathbf{R}}_{+,0}^n$ at $t=0$, the forward solution $x\left(t,x_0\right)\in {\mathbf{R}}_{+,0}^n$ for all $t\ge 0$. The boundary of the nonnegative orthant ${\mathbf{R}}_{+,0}^n$ is denoted as $\partial {\mathbf{R}}_{+,0}^n$ and consists of all points x where at least one component of x is zero, i.e., $x_i=0$ for some $i,$ and $x_j\ge 0$ for all $j\ne i$. Hence, the following Lemma provides a sufficient and necessary condition to establish positivity of system (A2), as detailed in Section $II.A$ in [21]:

Lemma A1.  

Positiveness of solutions. The nonlinear system $\dot{x}=f\left(x\right)$ is positive if and only if the following condition holds:

$$P\forall x\in \partial {\mathbf{R}}_{+,0}^n|x=0\Rightarrow f\left(x\right)\ge 0.$$

This condition ensures that the vector field $f\left(x\right)$ does not point out of ${\mathbf{R}}_{+,0}^n$ at its boundary $\partial {\mathbf{R}}_{+,0}^n$. Therefore, given nonnegative initial conditions, all solutions $x\left(t\right)$ will remain within ${\mathbf{R}}_{+,0}^n$ for all $t\ge 0$.

Appendix A.2. Localization of Compact Invariant Sets Method

The Localization of Compact Invariant Sets (LCIS) method [47,48] was formulated to study the long-term dynamics of nonlinear systems of first-order ODEs. Examples of compact invariant sets are the following: equilibrium points; periodic, homoclinic, and heteroclinic orbits; limit cycles, and chaotic attractors. The General Theorem is written as follows for a system of the form (A1), as detailed in Section 2 in [49]:

Theorem A1.  

General theorem. Each compact invariant set Γ of $\dot{x}=f\left(t,x\right)$ is contained in the localizing domain:

$$K\left(h\right)=\left\{x\in {\mathbf{R}}^n|h_{inf}\le h\left(x\right)\le h_{\sup }\right\}.$$

In order to apply the LCIS method, first one should propose the so-called localizing function $h\left(x\right)$, $x\in {\mathbf{R}}^n$, which is defined on the space $\mathbf{R}\times {\mathbf{R}}^{\mathbf{n}}=\left\{\left(t,x\right)\right\}$, and it is assumed that $h\left(x\right)$ is not the first integral of (A1). Further, ${h|}_S$ denotes the restriction of $h\left(x\right)$ on a set $S\subset {\mathbf{R}}^n$. Thus, $S\left(h\right)$ denotes the set $\left\{\left(t,x\right)\in \mathbf{R}\times {\mathbf{R}}^{\mathbf{n}}\mid {\left.L_{f}h\right|}_{\left(t,x\right)}=0\right\}$, where $L_{f}h$ represents the Lie derivative (temporal derivative) of $f\left(t,x\right)$ and it is given by ${\left.L_{f}h\right|}_{\left(t,x\right)}={\sum }_{i=1}^n\left[\partial h/\partial x_i\right]f_i\left(t,x\right)$. From the latter, one can define $h_{inf}=inf\left\{h\left(x\right)\mid x\in S\left(h\right)\right\}$ and $h_{sup}=sup\left\{h\left(x\right)\mid x\in S\left(h\right)\right\}$. Now, a refinement of Theorem A1 is realized with help of the Iterative theorem stated as follows:

Theorem A2.  

Iterative theorem. Let $h_m\left(x\right),m=0,1,2,\cdots$ be a sequence of $C^{\infty }-$ differentiable functions. Sets

$$K_0=K\left(h_0\right),K_m=K_{m-1}\cap K_{m-1,m},m>0,$$

with

$$\begin{array}{c}{K}_{m-1,m}=\left\{x:h_{m,inf}\le h_m\left(x\right)\le h_{m,sup}\right\},\hfill \\ h_{m,sup}=\underset{S\left(h_m\right)\cap K_{m-1}}{sup}{h}_m\left(x\right),\hfill \\ h_{m,inf}=\underset{S\left(h_m\right)\cap K_{m-1}}{inf}{h}_m\left(x\right),\hfill \end{array}$$

contain any compact invariant set of the system $\dot{x}=f\left(t,x\right)$ and

$$K_0\supseteq K_1\supseteq \cdots \supseteq K_m\supseteq \cdots .$$

The distinctiveness of the LCIS method lies in its strictly analytical nature, enabling the solution of the problem without solving the system of ODEs.

Appendix A.3. Stability in the Sense of Lyapunov

Lyapunov demonstrated that certain types of functions, called Lyapunov candidate functions, can be used to investigate the stability of the equilibrium point $x^{*}=\mathbf{0}$ of a nonlinear system of the form (A2). A Lyapunov candidate function $V\left(x\right):{\mathbf{R}}^n\to \mathbf{R}$ must satisfy the following properties:

• Positive definiteness: $V\left(0\right)=0$ and $V\left(x\right)>0$ for all $x\ne 0$.
• Radial unboundedness: $V\left(x\right)\to \infty$ as $∥x∥\to \infty$.

The temporal derivative of $V\left(x\right)$ along the trajectories of the system is given by $\dot{V}\left(x\right)=\left[\partial V/\partial x\right]f\left(x\right)$, where $\dot{V}\left(x\right)$ is a negative definite function if $\dot{V}\left(0\right)=0$ and $\dot{V}\left(x\right)<0$ for all $x\ne 0$. Thus, the Theorem of global asymptotic stability in the sense of Lyapunov is formulated as below, as detailed in Section $4.1$ in [50] and Chapter 2 in [51]:

Theorem A3.  

Lyapunov’s direct method.The equilibrium point $x^{*}=\mathbf{0}$ of $\dot{x}=f\left(x\right)$ is globally asymptotically stable if there exists a positive definite and radially unbounded function $V\left(x\right)$ such that $\dot{V}\left(x\right)$ is negative definite.

A function $V\left(x\right)$ that satisfies the properties of this theorem is called a Lyapunov function. However, there is no explicit method to find a candidate, and they should be formulated by trial and error. It is also important to note that the failure of a Lyapunov candidate function to satisfy the conditions of Theorem A3 does not necessarily imply that the equilibrium is unstable. It only indicates that stability in the sense of Lyapunov cannot be established with the chosen candidate.

In a small neighborhood of the origin, a nonlinear system of the form (A2) can be approximated by its linearization about this point:

$$\dot{x}=Ax,\mathrm{where}A={\left.{\displaystyle \frac{\partial f}{\partial x}}\right|}_{x=0}.$$

The following theorem provides conditions under which the stability of the origin for the nonlinear system (A2) can be inferred from the stability of the linearized system, as detailed in Section 4.3 in [50]:

Theorem A4.  

Lyapunov’s indirect method. Let the equilibrium point $x^{*}=\mathbf{0}$ of $\dot{x}=f\left(x\right)$ where $f:D\to {\mathbf{R}}^n$ is continuously differentiable and D is a neighborhood of the origin. Then:

• The origin is asymptotically stable if all eigenvalues of A satisfy $Re{\lambda }_i<0$.
• The origin is unstable if one or more eigenvalues of Asatisfy $Re{\lambda }_i>0$.

Appendix A.4. Lipschitz Condition

In order to predict the future state of a system described by a nonlinear dynamical system of the form (A1) from its present state at $t_0$, the initial value problem:

$$\dot{x}=f\left(t,x\right),x\left(t_0\right)=x_0,$$

must have a unique solution. The existence and uniqueness of such a solution can be guaranteed by imposing a constraint on the right-hand side function $f\left(t,x\right)$, known as the Lipschitz condition. This condition requires $f\left(t,x\right)$ to satisfy the inequality:

$$∥f\left(t,x\right)-f\left(t,y\right)∥\le \ell ∥x-y∥,$$

for all $\left(t,x\right)$ and $\left(t,y\right)$ in some neighborhood of $\left(t_0,x_0\right)$, where ℓ is the Lipschitz constant. The following results, detailed in Section 3.1 of [50], provide further insight into the Lipschitz condition:

Lemma A2.  

Locally Lipschitz. If $f\left(t,x\right)$ and $\left[\partial f/\partial x\right]\left(t,x\right)$are continuous on $\left[a,b\right]\times D$, with $\left[a,b\right]\subset \left[t_0,\infty \right)$ for some domain $D\subset {\mathbf{R}}^n$, then f is locally Lipschitz in x on $\left[a,b\right]\times D$.

Theorem A5.  

Existence and uniqueness. Let $f\left(t,x\right)$ be piecewise continuous in t and locally Lipschitz in x for all $t\ge t_0$ and x in a domain $D\subset {\mathbf{R}}^n$. Let K be a compact subset of D, $x_0\in K$, and suppose it is known that every solution of $\dot{x}=f\left(t,x\right),$$x\left(t_0\right)=x_0$, lies entirely in K. Then, there exists a unique solution defined for al $t\ge t_0$.

The key challenge in applying Theorem A5 lies in verifying that every solution remains within a compact domain K without explicitly solving the system of ODEs.