Keywords

1 Introduction

Leukemia is a cancer of the blood and bone marrow, characterized by large and uncontrolled growth of white blood cells. The most studied type of leukemia, Chronic myelogenous leukemia (CML), involves granular leukocyte precursors, namely the myelocyte line. The trigger of CML is a chromosomal abnormality, called the Philadelphia chromosome, that occurs in all cell lineages in about 90% of cases. The product of this chromosome is the formation of the BCR-ABl fusion protein which is thought to be responsible for the dysfunctional regulation of myelocyte proliferation and other features of CML. The standard treatment of CML in recent years is Imatinib, a molecular targeted drug [15], that has the effect that almost all patients attain hematological remission [13] and 75% attain cytogenetic remission.

Nowadays, it is well known that both innate and adaptive immunity are implicated in the defense mechanisms against cancer and recent progress in cancer immunology suggest that the immune system plays a fundamental role in tumor progression [23]. In CML, the biological literature reveals that T cells may play an important role in stemming the expansion of leukemic cells. This happens because leukemic cells express antigens that are immunogenic and can be recognized by cytotoxic T cells (CD8+ T cells or CTLs). In the paper [22], the authors found that leukemia-specific effectors CTLs were able to eliminate leukemic stem cells (LSCs) in vitro and in vivo in a setting with minimal leukemia load. The role of CD4+ T cells in leukemia is less clear, although in [3] the authors ascertained that some CML patients under imatinib-induced remission develop an anti-leukemia immune response involving both CD4+ and CD8+ T cells. Hence, our goal in this paper is to capture in a mathematical model the underlying dynamics of this disease by considering the evolution of healthy and leukemic cell populations along with one of the most important component of the cellular immune response to CML, namely T cell response.

Even if a variety of mathematical papers have applied a range of modeling approaches to study tumor-immune interactions in general (see, for example the recent review [7]), only a few described the specific leukemia-immune interaction. Leukemia-immune models have been formulated using mostly ordinary differential equations (ODE) [16, 17] or delay differential equations (DDE) [2, 5, 12, 19]. Some models that specifically study the immune response to CML are [2, 12, 16, 18]. However, none of the above papers have considered competition between healthy and leukemic cell populations, which is an important factor in CML dynamics.

2 Assumptions on the Model

In this paper, we use a five-dimensional system of DDEs. The first four equations describing the healthy and leukemic cell populations are based on the Mackey and collaborators models of hematopoiesis [14]. In the present model, we consider two types of hematopoietic cell populations: stem-like cell populations consisting of stem cells and progenitors with self-renew ability and mature cell populations formed by differentiated cells without self-renew ability. We include two types of cell populations, healthy and leukemic, each of them with two subpopulations of cells: stem-like and mature. In this way, we introduce space competition between the normal cell population and the CML one. We underline that in the case of leukemic cell populations mature cells are mostly unable to perform their functions.

The main difference from the Mackey model is considering the competition between healthy and CML cell populations and the fact that three types of division of a stem-like cell are considered: asymmetric division, symmetric self renewal and differentiation. In this paper we use the notation \(\alpha =h,l\) with h for healthy and l for leukemia. Consequently, we assume that a fraction \(\eta _{1\alpha }\), \(\alpha =h,l\), of stem-like cell population is susceptible to asymmetric division: one daughter cell proceeds to differentiate and the other re-enters the stem cell compartment. A fraction \(\eta _{2\alpha }\), \(\alpha =h,l\), is susceptible to differentiate symmetrically with both cells that result following a phase of maturation and the fraction \(1-\eta _{1\alpha }-\eta _{2\alpha }\), \(\alpha =h,l\), is susceptible to self-renewal so both cells that results after mitosis are stem-like cells. This four-dimensional competition model was introduced by Radulescu et al. in [20].

The fifth equation of the system models the anti-leukemia T cell immune response. The T cell population considered in this paper consists only of activated anti-leukemia T cells (involving both CD4+ and CD8+ T cells), which actively interact with CML cells. We do not consider other parts of the immune response, as the population of antigen presenting cells or the levels of citokines production. We assume that after encountering a leukemic cell, a T cell has two possibilities: either it inhibits the leukemic cell and activates a feedback function to stimulate the production of new T cells, or it is inhibited itself by the leukemic cell.

3 Description of the Model

The state variables of the model are the healthy cell populations \(x_{1}\) - stem-like and \(x_{2}\) - mature, the CML cell populations \(x_{3}\) - stem-like and \(x_{4}\) - more differentiated and \( x_{5}\) - the population of anti-leukemia T cells. The delays for healthy and leukemia stem-like cells are \(\tau _{1}\) and \(\tau _{3}\) for the duration of the cell cycle, independent of the type of division, and \(\tau _{2}\) and \(\tau _{4}\) for the time necessary for differentiation into mature leukocytes for healthy and, respectively, leukemia cells. \(\tau \) is the duration of the cell cycle for T-cells and \(\tau _{5}=n\tau \) with n the number of antigen depending divisions. We denote \(X_{\tau _{i} }=X(t-\tau _{i} ),\) where \(X=(x_{1},x_{2},x_{3},x_{4},x_{5})\).

The optimal control model is

$$\begin{aligned} \begin{aligned} \dot{x}_{1}&= f_{1}(x_{1},x_{2},x_{3},x_{4},x_{1\tau _{1}},x_{2\tau _{1}},x_{3\tau _{1}},x_{4\tau _{1}}) \\ \dot{x}_{2}&= f_{2}(x_{2},x_{1\tau _{2}},x_{2\tau _{2}},x_{4\tau _{2}}) \\ \dot{x}_{3}&= f_{3}(x_{1},x_{2},x_{3},x_{4},x_{5},x_{1 \tau _{3}},x_{2\tau _{3}},x_{3\tau _{3}},x_{4\tau _{3}},u_{1},k_{1l},u_{1\tau _{3}}) \\ \dot{x}_{4}&= f_{4}(x_{3},x_{4},x_{5},x_{2\tau _{4}},x_{3\tau _{4}},x_{4\tau _{4}},k_{2l},u_{1\tau _{4}}) \\ \dot{x}_{5}&= f_{5}(x_{4},x_{5},x_{4\tau _{5}},x_{5\tau _{5}},u_{2}) \end{aligned} \end{aligned}$$
(1)

where

$$\begin{aligned} \begin{aligned} f_{1}&= -\gamma _{1h}x_{1}-(\eta _{1h}+\eta _{2h})k_{h}(x_{2}+x_{4})x_{1}-(1-\eta _{1h}-\eta _{2h})\beta _{h}(x_{1}+x_{3})x_{1}+ \\&\quad +2e^{-\gamma _{1h}\tau _{1}}(1-\eta _{1h}-\eta _{2h})\beta _{h}(x_{1\tau _{1}}+x_{3\tau _{1}})x_{1\tau _{1}}+\eta _{1h}e^{-\gamma _{1h}\tau _{1}}k_{h}(x_{2\tau _{1}}+x_{4\tau _{1}})x_{1\tau _{1}} \\ f_{2}&= -\gamma _{2h}x_{2}+A_{h}(2\eta _{2h}+\eta _{1h})k_{h}(x_{2\tau _{2}}+x_{4\tau _{2}})x_{1\tau _{2}} \\ f_{3}&= -(\gamma _{1l}+\mathbf {f}_{1a})x_{3}-\left[ (\eta _{1l}+\eta _{2l})k_{l}(\left( x_{2}+x_{4}\right) \mathbf {f} _{u_{1}})+(1-\eta _{1l}-\eta _{2l})\beta _{l}(\left( x_{1}+x_{3}\right) \mathbf {f}_{u_{1}})\right] x_{3}+ \\&\quad +[2e^{-\left( \gamma _{1l}+\mathbf {\tilde{f}} _{1a}\right) \tau _{3}}(1-\eta _{1l}-\eta _{2l})\beta _{l}(\left( x_{1\tau _{3}}+x_{3\tau _{3}}\right) \mathbf {f}_{\mathbf {u}_{1\tau _{3}}})+ \\&\quad +\eta _{1l}e^{-\left( \gamma _{1l}+\mathbf {\tilde{f}} _{1a}\right) \tau _{3}}k_{l}(\left( x_{2\tau _{3}}+x_{4\tau _{3}}\right) \mathbf {f}_{\mathbf {u}_{1\tau _{3}}})]x_{3\tau _{3}}-b_{1}x_{3}x_{5}l_{1}(x_{3}+x_{4}) \\ f_{4}&= -\left( \gamma _{2l}+\mathbf {f}_{2a}\right) x_{4}+A_{l}(2\eta _{2l}+\eta _{1l})k_{l}(\left( x_{2\tau _{4}}+x_{4\tau _{4}}\right) \mathbf {f}_{\mathbf {u}_{\tau _{4}}})x_{3\tau _{4}}-b_{2}x_{4}x_{5}l_{1}(x_{3}+x_{4}) \\ f_{5}&= a_{1}-a_{2}x_{5}-a_{3}\mathbf {f} _{u_{2}}x_{5}l_{2}(x_{4})+2^{n_{1}}a_{4}x_{5\tau _{5}}l_{2}(x_{4\tau _{5}}) \end{aligned} \end{aligned}$$

subject to minimization of the cost functional

$$\begin{aligned} \min J(u), \end{aligned}$$
(2)

where

$$\begin{aligned} J(u)=g(x(T))+L(t,u(t),x_{5}(t)) \end{aligned}$$

with \(g(x(T))=A_{1}x_{3}(T)+A_{2}x_{4}(T)+E_{1}x_{3}(T)/x_{1}(T)+E_{2}x_{4}(T)/x_{2}(T)\) – being the weighted sum of the final tumor population and the ratio the ratio of leukemia cells and the healthy ones and

\(L(t,u(t),x_{5}(t))=\int \limits _{0}^{T}\left[ B_{1}u_{1}(t)+B_{2}u_{2}(t)+C_{1}k_{1}(t)+C_{2}k_{2}(t)-Dx_{5}(t)\right] dt\) – the cumulative drug toxicity and T cell amount.

The history of the state variables is given by \(X(\theta )=\varphi (\theta ), \theta \in [-\tau _{\max },0], \tau _{\max }=\max (\tau _{1},\tau _{2},\tau _{3},\tau _{4},\tau _{5})\).

The healthy and leukemic blood cell populations are seen in competition for resources and this is reflected in the fact that both feedback laws for self-renewal and differentiation depend on the sum of healthy and leukemia cells. Consequently, the rate of self-renewal is \(\beta _{\alpha }(x_{1}+x_{3})=\beta _{0\alpha }\dfrac{\theta _{1\alpha }^{m_{\alpha }}}{\theta _{1\alpha }^{m_{\alpha }}+(x_{1}+x_{3})^{m_{\alpha }}}\), with \(\beta _{0\alpha }\) the maximal rate of self-renewal and \(\theta _{1\alpha }\) half of the maximal value and the rate of differentiation is \(k_{\alpha }(x_{2}+x_{4})=k_{0\alpha }\dfrac{\theta _{2\alpha }^{n_{\alpha }}}{\theta _{2\alpha }^{n_{\alpha }}+(x_{2}+x_{4})^{n_{\alpha }}}\), with \(k_{0\alpha }\) the maximal rate of differentiation and \(\theta _{2\alpha }\) is half of the maximal value. The rest of the parameters for healthy and CML cell populations are: \(\gamma _{1\alpha }\) - the natural apoptosis, \(A_{\alpha }\) - an amplification factor and \( m_{\alpha },n_{\alpha }\) parameters that control the sensitivity of \(\beta _{\alpha }\) respectively \(k_{\alpha }\) to changes in the size of stem-like and respectively mature populations. Table 1 contains a complete description of parameters of the model.

To model the influence of T cells on CML cells, we consider the feedback function \(l_{1}(y)=\dfrac{1}{b_{5}+y}.\) Consequently, the last terms of the third and fourth equations represent the inhibition of CML cells by anti-leukemia T cells. We assumed that the inhibition of CML cell population by T cells increases with the number of leukemic cells up to a certain level and then reaches a maximal value of inhibition. A further increase in CML population will not modify this value.

As concerns the fifth equation, the first term \(a_{1}\) is the natural supply of naive T cells, while the second term \(-a_{2}x_{5}\) indicates that T cells exit the population through death at the rate \(a_{2}\). Leukemia cells suppress anti-leukemia immune response. The precise mechanism is unknown, but it is assumed that the level of down regulation depends on the current leukemia population, so we consider that the immune system is regulated by the feedback function \(l_{2}(y)=\dfrac{y}{b_{5}+y^{2}}.\) This function ensures that T cells are stimulated by CML cells only if leukemia cell population has values in a certain range, called “the optimal load zone” (see [12]). We take the rate of antigen stimulation as a feedback function depending on the level of the mature leukemic population, \( l_{2}(x_{4})\) and the third and the fourth terms, \( -a_{3}u_{2}x_{5}l_{2}(x_{4})\) and \(2^{n_{1}}a_{4}x_{5\tau _{6}}l_{2}(x_{4\tau _{5}})\) gives the rate at which naive T cells leave and re-enter the effector state after finishing the minimal developmental program of \(n_{1}\) cell divisions (due to antigen stimulation). The time delay \(\tau _{5}=n_{1}\tau \) is the duration of this program. These terms represent the loss and respectively the production of T cells due to the competition with leukemic cells.

The treatment targeting the BCR-ABL gene is supposed to affect the apoptosis and the proliferation rates of leukemia cells [11]. In view of this fact, we consider the treatment functions \(f_{u_{1}}=\dfrac{1}{1-u_{1}}\), \(f_{1a}=\left( \gamma _{1h}-\gamma _{1l}\right) k_{1l}\) and \(f_{2a}=c\gamma _{2h}k_{2l}\), with \(u_{1},k_{1l},k_{2l}:[0,T]\rightarrow [0,1]\), where \(u_{1}(t),k_{1l}(t),k_{2l}(t)\) are the treatment effects.

The action of treatment on the proliferation rate will be considered through \(f_{u_{1}}\) in the function of self-renew \(\beta _{l}\) and in the function of differentiation or asymmetric division \(k_{l}\). Note that, in this way, both \(\beta _{l}\) and \(k_{l}\) became decreasing functions of \(u_{1}\). For more details, see [21]. The treatment acts on the apoptosis through the function \(f_{1a}\) on the stem cells and through the function \(f_{2a}\) on the mature ones. Also, from the law of the mass, we have \(\tilde{f}_{1a}=\int \limits _{t-\tau _{1}}^{t}k_{1l}(s)ds\). Moreover, it seems that, in vivo, Imatinib is the trigger of complex mechanisms, some of them able to promote T cell expansion (see [22]). Imatinib’s effect on T cell population is introduced in the form of a treatment function \(f_{u_{2}}=1-u_{2}\), with the stimulatory effect \( u_{2}:[0,T]\rightarrow [0,1]\). If no drug is given, then \(f_{u_{2}}=1\) and a maximal effect takes place for \(u_{2}(t)=1\), when T cell population is no longer inhibited by CML cell population.

The existence of an optimal control follows since one can transform the given problem into an optimal control problem for a system of ODEs whose solutions will be bounded together with their derivatives on compact intervals (see [1]).

4 Discretization of the Optimal Problem

In this section, we apply the numerical procedure from Gollmann and Maurer [10], in order to solve the delay optimal control problem (1) + (2) (see also [8, 9]). For that matter, we write the cost functional in the Mayer form

$$\begin{aligned} J(u,x)=h(x(T)),\,\, x=(x_{1},x_{2},x_{3},x_{4},x_{5})\in R^{5}. \end{aligned}$$

We introduce the additional state variable z through the equation

$$\begin{aligned} \dot{z} (t)=B_{1}u_{1}(t)+B_{2}u_{2}(t)+C_{1}k_{1}(t)+C_{2}k_{2}(t)-Dx_{5}(t),\,\, z(0)=0. \end{aligned}$$

Then, the cost functional (2) is rewritten as

$$\begin{aligned} J(u,x,z)=g(x(T))+z(T). \end{aligned}$$

In the following, let \(\tau >0\) be such that \(\tau _{1}=j_{1}\tau \), \(\tau _{2}=j_{2}\tau , \tau _{3} =j_{3}\tau \), \(\tau _{4}=j_{4}\tau ,\tau _{5}=j_{5}\tau , j_{i}\in N^{*}\), \(i=\overline{1,5}\), \(T=N\tau \) and use the Euler integration method with a uniform step size \(\tau >0\). Of course, \(\tau \) can be refined in order to obtain an appropriate smaller step-size. Using the grid points \(t_{i}=i\tau \), \(i=\overline{0,N}\) and the approximations \(x_{1}(t_{i})\simeq x_{1i}\in R\), \(x_{2}(t_{i})\simeq x_{2i}\in R\), \(x_{3}(t_{i})\simeq x_{3i}\in R\), \(x_{4}(t_{i})\simeq x_{4i}\in R\), \(x_{5}(t_{i})\simeq x_{5i}\in R\), \(u_{1}(t_{i})\simeq u_{1i}\), \(u_{2}(t_{i})\simeq u_{2i}\) and \(k_{1}(t_{i})\simeq k_{1i}\), \(k_{2}(t_{i})\simeq k_{2i}\) the treatment function \(f_{1a}\) becomes \( \sum \limits _{j=1}^{j_{3}}k_{1l_{i-j}}\tau \) and the delay control problem (1) + (2) is transformed into the nonlinear programming problem (NLP)

$$\begin{aligned} Minimize\,J=g(x_{N})+z_{N} \end{aligned}$$
(3)

subject to

$$\begin{aligned} \left\{ \begin{array}{l} x_{1i}-x_{1i+1}+\tau f_{1}(x_{1i},x_{2i},x_{3i},x_{4i},x_{1i-j_{1}},x_{2i-j_{1}},x_{3i-j_{1}},x_{4i-k_{1}})=0 \\ x_{2i}-x_{2i+1}+\tau f_{2}(x_{2i},x_{1i-j_{2}},x_{2i-j_{2}},x_{4i-j_{2}})=0 \\ x_{3i}-x_{3i+1}+\tau f_{3}(x_{1i},x_{2i},x_{3i},x_{4i},x_{5i},x_{1i-j_{3}},x_{2i-j_{3}},x_{3i-j_{3}},x_{4i-j_{3}},u_{1i},\\ k_{1i},\sum \limits _{j=1}^{j_{3}}k_{1_{i-j}}\tau ,u_{1i-j_{3}})=0 \\ x_{4i}-x_{4i+1}+\tau f_{4}(x_{4i},x_{5i},x_{2i-j_{4}},x_{3i-j_{4}},x_{4i-j_{4}},k_{2li},u_{1i-j_{4}})=0 \\ x_{5i}-x_{5i+1}+\tau f_{5}(x_{4i},x_{5i},x_{4i-j_{5}},x_{5i-j_{5}},u_{2i})=0 \\ z_{i}-z_{i+1}+\tau (B_{1}u_{1i}+B_{2}u_{2i}+C_{1}k_{1i}+C_{2}k_{2i}-Dx_{5i})=0 \end{array} \right. \end{aligned}$$
(4)
$$\begin{aligned} -u_{1i}\le & {} 0,u_{1i}-1\le 0,-u_{2i}\le 0,u_{2i}-1\le 0, \\ -k_{1i}\le & {} 0,k_{1i}-1\le 0,-k_{2i}\le 0,k_{2i}-1\le 0, i=\overline{0,N-1}.. \nonumber \end{aligned}$$
(5)

Herein, the initial value profiles \(\varphi _{1},\varphi _{2},\varphi _{3},\varphi _{4}\) and \(\varphi _{5}\) give the values

$$\begin{aligned}\begin{gathered} x_{1_{-i}}:=\varphi _{1}(-i\tau ),\ i=\overline{0,l_{1},}x_{2_{-i}}:=\varphi _{2}(-i\tau ),\ i=\overline{0,l_{2}},x_{3_{-j}}:=\varphi _{3}(-i\tau ),\ i= \overline{0,l_{3}} \\ x_{4_{-j}}:=\varphi _{4}(-i\tau ),\ i=\overline{0,l_{4}},x_{5_{-j}}:=\varphi _{5}(-i\tau ),\ i=\overline{0,l_{5}}. \end{gathered}\end{aligned}$$

The variable to be optimized is represented by the vector \(w=(u_{1_{0}},u_{2_{0}},k_{1_{0}},k_{2_{0}}, x_{1_{1}}, .., , x_{5_{1}},z_{1}, .., u_{1_{N-1}},u_{2_{N-1}},k_{1_{N-1}}, k_{2_{N-1}}, x_{1_{N}},..,x_{5_{N}},z_{N}) \in R^{10N}\).

5 Numerical Results

In the following figures, we plotted the trajectories of the healthy, respectively CML cell populations for the competition system, showing a comparison between the dynamics of a system without treatment and the dynamics of a system subject to optimal treatment. In the following simulations two aspects are combined, resulting four distinct manifestations of the disease:

  • starting treatment in two different stages of the disease: a less severe stage when the population of leukemic cells and the healthy cells coexist (the leukemia cell population is still small) (S1) and a stage where healthy cell population disappeared and the number of leukemic cells is already very high (S2);

  • two configurations of parameters describing two different forms of the disease for patients, configuration 1 and configuration 2 (see Table 1). The configuration 2 corresponds to a more serious disease.

Considering multiple effects of treatment, on the apoptosis of leukemic stem cells, leukemic mature cells, proliferation rate and immune system, simulations show the impact of the disease on various optimal control solutions for four hypothetical patients (see Figs. 1, 2, 3 and 4 for comparison between the dynamics of a system without treatment and with optimal reatment).

Fig. 1.
figure 1

Simulations start from S1 for configuration 1 of parameters. Dashed line represents the dynamics of a system without treatment and continuous line represents the dynamics of a system with optimal treatment.

Full size image
Fig. 2.
figure 2

Simulations start from S2 for configuration 1 of parameters. Dashed line represents the dynamics of a system without treatment and continuous line represents the dynamics of a system with optimal treatment.

Full size image
Fig. 3.
figure 3

Simulations start from S1 for configuration 2 of parameters. Dashed line represents the dynamics of a system without treatment and continuous line represents the dynamics of a system with optimal treatment.

Full size image
Fig. 4.
figure 4

Simulations start from S2 for configuration 2 of parameters. Dashed line represents the dynamics of a system without treatment and continuous line represents the dynamics of a system with optimal treatment.

Full size image
Fig. 5.
figure 5

Controls for configuration 1 of parameters, simulations start from S1. The cost function was improved from 4000 to 2700.

Full size image
Fig. 6.
figure 6

Controls for configuration 1 of parameters, simulations start from S2. The cost function was improved from 5100 to 4100.

Full size image
Fig. 7.
figure 7

Controls for configuration 2 of parameters, simulations start from S1. The cost function was improved from 3400 to1600.

Full size image
Fig. 8.
figure 8

Controls for configuration 2 of parameters, simulations start from S2. The cost function was improved from 2700 to 2100.

Full size image
Table 1. Description of parameters
Full size table

In the Figs. 5, 6, 7 and 8 the controls \(k_{1l}\), \(k_{2l}\), \(u_1\), \(u_2\) represent the influence of drug on the apoptosis of leukemic stem cells, leukemic mature cells, proliferation rate and immune system. The value of cost functional was improved in all situations (see figures). To solve the problem of optimal control the Matlab solver for NLP problems fmincon was used, selecting the ‘interior-point’ solver.

6 Conclusions

In this paper, an optimal control model for CML with the influence of the immune system and treatment was investigated. Based on clinical evidences and assumptions, the effects of Imatinib, the current first line treatment in CML, was considered. These effects include the decrease of leukemic proliferation and differentiation, the increase of leukemic apoptosis and some influences on the anti-leukemia immune response.

From Figs. 1, 2, 3 and 4 one can see one can see the decline of leukemic cells (i.e. molecular remission) and the increase of the number of healthy cells after approximately three months of treatment. Depending on the level of cell populations at diagnosis (S1 or S2) and on the leukemia severity (i.e. configuration 1 or configuration 2) the evolution of healthy cell population to a normal amount is more or less rapid.

The plots of optimal controls (Figs. 5, 6, 7 and 8) exhibit an optimal control effect different for four hypothetical patients. One can observe that the drug influence is slightly different for various manifestations of the disease. Consequently, for an optimal effect of treatment, the prescribed dose should be adapted considering the parameter’s disease of a certain patient and the leukemic burden at diagnosis. Although the identification of most of the values parameters of the disease is a daunting task, there are some which can be computed by means of current methods and they might provide an important indication concerning dose adjustment and therapy management.