3.1. 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
) on CLL cells. As one can see below in
Table 2, increasing the drug concentration also increases cytotoxicity. Further, concentration was converted from µM (
mol/L) to mg/L, and cytotoxicity rate from h
to days
to standardize the time units of the mathematical model under study in this work.
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
values from the column ‘Experimental cytotoxicity in days
’ in
Table 2,
where
is the chlorambucil concentration in mg/L, and
represents its approximated cytotoxicity rate in days
. Experimental and fitted data is illustrated in
Figure 1, and the corresponding values are shown in the column ‘Fitted cytotoxicity
in days
’ in
Table 2 with a coefficient of determination
. Further, parameters
and
from Equation (
8) were computed with a
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]:
where
and
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
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) as follows
and by simplifying the latter
the next constant is identified
representing the total number of molecules in 1 mg of chlorambucil. At this step, Equation (
9) may be simplified
where
is the chlorambucil dose in mg. Hence, for 32 mg we have the following result
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
Now, given that the biological half-life
of chlorambucil has been reported as
days [
50,
51], one can determine its deactivation/decay rate (denoted as
) from the first-order pharmacokinetics Equation [
55]
with the following solution
for the initial condition
. Thus, from the latter we have the following
where
represents
of the initial dose at
. Therefore, by isolating
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
’ in
Table 4 were computed with a coefficient of determination
by applying the next Equation
which was formulated by taking into account both the mitosis stimulation rate
and natural death rate of CAR-T cells
, as well as the leukemia growth rate
. Further discussion on these parameters is provided in
Table 1 in
Section 2. Parameters
were computed with a
CI as it is shown in
Table 5.
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.
3.2. Localization of Compact Invariant Sets
Below, we provide the mathematical background that allows us to determine the localizing domain in the non-negative orthant
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:
Each compact invariant set of is contained in the localizing domain:From the latter we have that is a differentiable vector function where is the state vector. is a differentiable function called localizing function. denotes the restriction of on a set with , and is the Lie derivative of . Hence, one can define and . Furthermore, if all compact invariant sets are contained in the set and in the set then they are contained in as well. Nonexistence of compact invariant sets can be considered for a given set if , then the system 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
.
Maximum upper bound for the live leukemia cells population
. In order to determine this bound, the following localizing function is proposed
and its Lie derivative is defined as follows
now, one can write set
as shown below
at this step, negative terms of
can be discarded. Hence, non-negative boundaries for all
limit sets of
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 ) 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
. Let us explore the following localizing function
whose Lie derivative is defined as follows
then set
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
therefore,
now, due to
we may conclude the following result
Thus, from the latter, upper bounds for both the live
and dead
leukemia cells populations can be deducted as follows
Localizing set for the CAR-T cells therapy
. 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
now, one can write set
as indicated below
thus, by condition (
5), i.e.,
, we can complete the square
and rewrite
where
hence, the lower bound is determined by disregarding the non-negative function
and isolating
T, then, we can conclude on the following lower bound for all solutions of
from the latter, it is evident that
, 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
and write it as follows by considering
now, let us complete the square
to rewrite
where
therefore, we have the following upper bound for the localizing function
Now, one can formulate the following lower and upper bounds for all solutions of the live leukemia
and CAR-T
cells populations
Localizing set for the chemotherapy drug concentration
. Let us exploit the localizing function denoted by
where its Lie derivative is defined as follows
now, set
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
. However, if condition (
13) does not hold, then one should consider
as there is no biological sense for negative values for the chemotherapy drug concentration,
. Now, let us take again set
as follows
from which we can determine the next upper bound by disregarding the negative terms as given below
Hence, all solutions
will be bounded by the following domain
where
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
. Now, let us apply the Iterative Theorem to the set
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
as indicated below
which is solved for the CAR-T cells therapy treatment
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 for the leukemia and chemoimmunotherapy system (1)–(4). 3.3. 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
and compute its time derivative
as follows
Then, in order to fulfill Lyapunov’s asymptotic stability conditions
, the derivative is evaluated at the localizing domain, i.e.,
, and the next upper bound is determined
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:
Theorem 2. Leukemia cancer cells eradication. If the CAR-T cells therapy dose meets condition (14), then one can ensure the complete eradication of the leukemia cancer cells population described by the system (1)–(4). Thus, 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.,
. Therefore, Equations (
1)–(
4) become as follows
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:
Corollary 2. Leukemia-free equilibrium point. If condition (14) from Theorem 2 is fulfilled and the chemoimmunotherapy treatments are stopped when the leukemia cells are completely eradicated
, then the tumor-free equilibrium point (19) is globally asymptotically stable. 3.4. 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 nonlinear 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
of chlorambucil is
days, which implies that this drug should be depleted in the patient in a short period of time, i.e.,
. 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.,
∀
with
,
with
Now, as the method requires, the Jacobian matrix
is computed below
which is evaluated at the Equilibrium point (
23) as follows
at this step, eigenvalues are computed from the Jacobian determinant
, results are shown below
and by condition (
24) it is evident that
Re,
. Therefore, conditions for depletion (
6) and persistence (
7) are derived from
as follows. If
Re, then by Theorem
in [
61]: “
The equilibrium is locally asymptotically stable”, this holds when
which biologically implies depletion of the CAR-T cells. Thus, if
Re, then by Theorem
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., Now, following the latter, let us explore leukemia cells persistence under the immunotherapy treatment, i.e.,
. 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)–(
23)
where
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
,
. Thus, local asymptotic stability of (
26) follows from the next condition on
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):
Theorem 4. Leukemia cells persistence. If the CAR-T cells dose meets condition (28), then the immunotherapy treatment will not be able to eradicate the leukemia cells population. Therefore, cancer persits, i.e.,