# A Control Delay Differential Equations Model of Evolution of Normal and Leukemic Cell Populations Under Treatment

I. Rodica Rădulescu, Doina Cândea, Andrei Halanay

► To cite this version:

I. Rodica Rădulescu, Doina Cândea, Andrei Halanay. A Control Delay Differential Equations Model of Evolution of Normal and Leukemic Cell Populations Under Treatment. Christian Pötzsche; Clemens Heuberger; Barbara Kaltenbacher; Franz Rendl. System Modeling and Optimization: 26th IFIP TC 7 Conference, CSMO 2013, Klagenfurt, Austria, September 9–13, 2013, Revised Selected Papers, AICT-443, Springer Berlin Heidelberg, pp.257-266, 2014, IFIP Advances in Information and Communication Technology, 978-3-662-45503-6. 10.1007/978-3-662-45504-3_25. hal-01286684

HAL Id: hal-01286684

https://inria.hal.science/hal-01286684

Submitted on 11 Mar 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A control delay differential equations model of evolution of normal and leukemic cell populations under treatment*

I. Rodica Rădulescu**, Doina Cândea, and Andrei Halanay

POLITEHNICA University of Bucharest
Department of Mathematics and Informatics
Spiaul Independentei 313 RO-060042 Bucharest, Romania

**Abstract.** The dynamics and evolution of leukemia is determined by the interactions between normal and leukemic cells populations at every phase of the development of hematopoietic cells. For both types of cell populations, two subpopulations are considered, namely the stem-like cell population (i.e. with unlimited self-renew ability) and a more mature, differentiated one, possessing only the capability to undergo limited reproduction. Treatment effects are included in the model as functions of time and a cost functional is considered. The optimal control is obtained using a discretization scheme. Numerical results are discussed in relation to the medical interpretation.

**Keywords:** leukemia, asymmetric division, competition, optimal control, treatment

**AMS Classification:** 34K35, 37N25, 92C50, 93C23.

## 1 Introduction

For the description of biological processes implied in hematopoiesis, a mathematical model that includes time delays will be used. It is based on the mass action principle, in the spirit of [1, 2, 4, 12, 16] and [20]. Other authors ([14, 15]) used the more simple model from [18] for the dynamics of hematopoietic stem cells (HSC).

Chronic Myelogenous Leukemia (CML), also known as Chronic Granulocytic Leukemia, is a cancer of white blood cells. It is a clonal marrow stem cells disorder in which the main characteristic is the proliferation of granulocytes (neutrophils, eosinophils and basophils) and of their precursors in the bone marrow and the accumulation of these cells in the blood. It is a type of myeloproliferative disease associated to a chromosomal translocation called the Philadelphia chromosome (see also [6, 19]) presenting the oncogene BCR-ABL that encodes

<footnote>* This work was supported by CNCS-ROMANIA Grant ID-PCE-2011-3-0198.</footnote>

<footnote>** Corresponding author. E-mail: nicola.rodica@yahoo.com</footnote>
a tyrosine kinase protein. Tyrosine kinases are enzymes that play an important
role in tumor development by supporting cell growth through phosphorylation
of signaling proteins [11]. Understanding the molecular mechanism of CML per-
mitted the development of specific tyrosine kinase inhibitors (TKIs) as imatinib
(Gleevec), dasatinib or nilotinib. The standard first line therapy is nowadays
imatinib, which acts through competitive inhibition at the ATP-binding site of
the BCR-ABL enzyme, leading to the inhibition of tyrosine phosphorylation of
proteins involved in BCR-ABL signal transduction [7]. The molecular effect of
imatinib is mainly the inhibition of cell proliferation of BCR-ABL-positive cells
but, there are experimental evidences that, in the mature cell lines, the inhibition
of cell proliferation is followed by apoptosis [11]. Although imatinib has a very
good successful rate, there are many experimental evidences attesting it does
not affect quiescent stem cells deep in the bone marrow, and the consequence is
the disease reapers after the treatment is stopped.

In this paper, we study an optimal control delay differential equation model
of four cell populations, namely two healthy and two leukemic. For these classes
of cells, we consider a population of mature cells which lost their self-renew
ability and a population of stem-like cells involving a larger category consisting of
proliferating stem and progenitor cells with self-renew capacity. The emphasis in
this optimal control model is on establishing treatment strategies, considering the
competition of healthy vs. CML cell populations and three types of division that
a stem-like cell can exhibit: self-renew, asymmetric division and differentiation
([4, 17, 20, 21, 24]).

Of course, besides a correct mathematical model for the time evolution of
the studied cell populations, it is very important to model the treatment effect
as accurate as one can. Obviously, different drugs have different effects: some
affect not only leukemic populations but also healthy ones (the cytotoxic ones),
some kill the cells while others only delay or stop the division process (the
cytostatic ones). We take into account here the standard treatment protocol
with imatinib which, as we specified, acts by inhibiting the BCR-ABL signal
transduction. In this way, it restrains the proliferative advantage of the CML
cells and healthy cells regain their advantage. Moreover, it is straightforward
to say that imatinib restores most of the abnormal functions of the CML cells,
and the most important function affected by the drug is the division process of
CML cell population. However, it is uncertain to what extent it affects all three
kinds of division and therefore, we consider here the hypothesis that imatinib
influences self-renew, asymmetric division and differentiation equally.

## 2 Description of the model

In the present paper, is assumed that the hematopoietic stem cells that are
considered are in the proliferative phase or spend a short time into the resting
phase. These cells are called, following [17], Short-Term Hematopoietic Stem
Cells (ST-HSC). In what follows $x_1$ denotes the density of short-term stem-like
healthy cells, $x_2$ the density of mature healthy cells, $x_3$ the density of short-term stem-like leukemic cells, $x_4$ the density of mature leukemic cells.

The time necessary for a ST-HSC to complete a cycle of self-renewal, asymmetric division or differentiation is $\tau_{1l}$ for leukemic cells and $\tau_{1h}$ for the healthy ones, while the time necessary for the maturation of leukocytes is denoted by $\tau_{2l}$ in the case of leukemic cells and $\tau_{2h}$ for the healthy ones.

As we mentioned in the introduction, experimental evidences attest that the imatinib therapy affects primarily the proliferation rate and secondary, the apoptosis rate. In view of this fact, we consider the treatment functions $f_u = \frac{1}{1-u}$ and $f_{1a} = (\gamma_{1h} - \gamma_{1l}) u_1$, $f_{2a} = c\gamma_{2h}u_2$, with $u, u_1, u_2: [0, T] \to [0, 1]$, where $u(t), u_1(t), u_2(t)$ are the treatment effects. The action of treatment on the proliferation rate will be considered through $f_u$ in the function of self-renewal $\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$. If no drug is given (i.e. $u(t) = 0$) then $\beta_l((x_1 + y_1)f_u) = \beta_l(x_1 + y_1)$, $k_l((x_1 + y_1)f_u) = k_l(x_2 + y_2)$ and also, a maximal effect happens for $u(t) = 1$ when the process of division essentially stops ($\beta_l \equiv k_l \equiv 0$).

Treatment will be considered to act only on the leukemic stem cells compartment. The treatment acts on the apoptosis of mature CML cells through $f_{2a}$ and on the apoptosis of stem-like CML cells through $f_{1a}$, restoring this rate to a value closed to the mortality rate of healthy cells. From the law of the mass, we have $\tilde{f}_{1a} = \int_{t-\tau_{1l}}^{t} u_1(s) ds$. The optimal control model is

$$ 
\begin{aligned}
\dot{x}_1 &= f_1(x_1, x_2, y_1, y_2, x_{1\tau_{1h}}, x_{2\tau_{1h}}, y_{1\tau_{1h}}, y_{2\tau_{1h}}) \\
\dot{x}_2 &= f_2(x_2, x_{1\tau_{2h}}, x_{2\tau_{2h}}, y_{2\tau_{2h}}) \\
\dot{y}_1 &= f_3(t, x_1, x_2, y_1, y_2, x_{1\tau_{1l}}, x_{2\tau_{1l}}, y_{1\tau_{1l}}, y_{2\tau_{1l}}, u_1, u, u_{\tau_{1l}}) \\
\dot{y}_2 &= f_4(y_2, x_{2\tau_{2l}}, y_{1\tau_{2l}}, y_{2\tau_{2l}}, u_2, u_{\tau_{2l}})
\end{aligned}
\quad (1) $$

where

$$ 
\begin{aligned}
f_1 &= -\gamma_{1h}x_1 - (\eta_{1h} + \eta_{2h})k_h(x_2 + y_2)x_1 - (1 - \eta_{1h} - \eta_{2h})\beta_h(x_1 + y_1)x_1 + \\
&\quad +2e^{-b_{1h}\tau_{1h}}(1 - \eta_{1h} - \eta_{2h})\beta_h(x_{1\tau_{1h}} + y_{1\tau_{1h}})x_{1\tau_{1h}} + \\
&\quad +\eta_{1h}e^{-b_{1h}\tau_{1h}}k_h(x_{2\tau_{1h}} + y_{2\tau_{1h}})x_{1\tau_{1h}} \\
f_2 &= -\gamma_{2h}x_2 + A_h k_h(x_{2\tau_{2h}} + y_{2\tau_{2h}})x_{1\tau_{2h}} \\
f_3 &= -(\gamma_{1l} + \mathbf{f}_{1a})y_1 - [(\eta_{1l} + \eta_{2l})k_l((x_2 + y_2)\mathbf{f}_u) + (1 - \eta_{1l} - \eta_{2l})\beta_l((x_1 + y_1)\mathbf{f}_u)]y_1 + \\
&\quad +[2e^{-b_{1l}\tau_{1l}} - \tilde{\mathbf{f}}_{1a}(1 - \eta_{1l} - \eta_{2l})\beta_l((x_{1\tau_{1l}} + y_{1\tau_{1l}})\mathbf{f}_{u_{\tau_{1l}}}) + \\
&\quad +\eta_{1l}e^{-b_{1l}\tau_{1l}} - \tilde{\mathbf{f}}_{1a}k_l((x_{2\tau_{1l}} + y_{2\tau_{1l}})\mathbf{f}_{u_{\tau_{1l}}})]y_{1\tau_{1l}} \\
f_4 &= -(\gamma_{2l} + \mathbf{f}_{2a})y_2 + A_l k_l((x_{2\tau_{2l}} + y_{2\tau_{2l}})\mathbf{f}_{u_{\tau_{2l}}})y_{1\tau_{2l}}
\end{aligned}
$$

subject to minimization of the cost functional

$$ \min J(u), \quad (2) $$
where

$$J(u) = ay_1(T) + by_2(T) + \int_{0}^{T} [u^2(t) + u_1^2(t) + u_2^2(t)] dt$$

with $g(y(T)) = ay_1(T) + by_2(T)$ the weighted sum of the final tumor population and $L(u(t)) = \int_{0}^{T} [u^2(t) + u_1^2(t) + u_2^2(t)] dt$, the cumulative drug toxicity.

In this paper we denote $X_{\tau} = X(t - \tau)$, where $X = (x_1, x_2, y_1, y_2)$. The history of the state variables is given by

$$X(\theta) = \varphi(\theta), \theta \in [-\tau_{\max}, 0], \tau_{\max} = \max(\tau_{1h}, \tau_{2h}, \tau_{1l}, \tau_{2l}).$$

It is not difficult to see that if the initial conditions have all components positive functions, the solutions of the system will have positive components on all the interval of existence. Indeed, if $x_1(\theta) > 0$ for any $\theta \in [-\tau, 0]$ and there exists $T > 0$ such that $x_1(T) = 0$ the derivative $\dot{x}_1(T)$ will be positive and this leads to a contradiction. The same argument works for the other components, too.

We assume that: a percentage $\eta_{1\alpha}$, $\alpha = h, l$ of stem-like cells population is supposed to undergo asymmetric division; a percentage $\eta_{2\alpha}$, $\alpha = h, l$ of the population differentiate symmetrically and the percentage $(1-\eta_{1\alpha}-\eta_{2\alpha})$, $\alpha = h, l$ of the population is supposed to self-renew (see also [13]).

Furthermore, it is assumed that homeostatic mechanisms maintain the hematopoietic stem cell population at a constant level. In this respect, the rate of self-renewal is given by a Hill function

$$\beta_{\alpha}(X) = \beta_{0\alpha} \frac{\theta_{1}^{m}}{\theta_{1}^{m} + X^{m}}, \alpha = h, l$$

and the rate of differentiation, through symmetric or asymmetric division is supposed to be dictated, through a feedback law, by

$$k_{\alpha}(X) = k_{0\alpha} \frac{\theta_{2}^{n}}{\theta_{2}^{n} + X^{n}}, \alpha = h, l.$$

Because in this paper we consider competition between healthy and CML cell populations, both this rates will depend on the sum of stem-like respectively mature populations (similar approaches on competition were modeled in [23, 22]).

For $\alpha = h, l$, the other parameter are defined as follows: $b_{1\alpha}$ accounts for the death rate of stem cells and a positive $K_{\alpha}$ for the loss rate due to differentiation into other cell lines - the resulting loss rate is denoted as $\gamma_{1\alpha} = K_{\alpha} + b_{1\alpha}$; $\beta_{0\alpha}$ and $k_{0\alpha}$ represent the maximal rate of self-renewal, respectively of asymmetric division or differentiation into leukocyte line; $\theta_i, i = 1, 2$, is the value for which $\beta_{\alpha}$, respectively $k_{\alpha}$ attains half of their maximum value; $\gamma_{2\alpha}$ is the mortality of mature cells; $A_{\alpha}$ is an amplification factor of mature cells due to differentiation; $m$ is the parameter controlling the sensitivity of the mitotic re-entry rate $\beta_{\alpha}$ to changes in the size of $G_0$ and $n$ is the parameter controlling the sensitivity of the
asymmetric division or differentiation rate $k_\alpha$ to changes in the size of mature population.

**Existence of an optimal control** The existence of an optimal control results from transforming the problem into an optimal control problem for a system of ordinary differential equations (see next section) whose solutions will be bounded together with their derivatives on compact intervals (see [5]).

### 3 Discretization of the optimal control problem

In this section, we apply the numerical procedure from Gollmann et al. [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

$$J(u, y) = h(y(T)), \quad y = (y_1, y_2) \in \mathbb{R}^2.$$

In our case, the reduction of the more general cost functional (2) to Mayer form, proceeds by the introduction of the additional state variable $z$ through the delayed equation

$$\dot{z}(t) = L(u(t)) \quad \text{so} \quad \dot{z}(t) = u^2(t) + u_1^2(t) + u_2^2(t), \quad z(0) = 0.$$

Then, the cost functional (2) is rewritten as

$$J(u, y, z) = g(y(T)) + z(T).$$

In the following, let $\tau > 0$ such that $\tau_{1h} = k_1\tau$, $\tau_{2h} = k_2\tau$, $\tau_{1l} = k_3\tau$, $\tau_{2l} = k_4\tau$, $k_i \in \mathbb{N}^*$, $i = 1, 4$, $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 = 0, \overline{N}$ and the approximations $x_1(t_i) \simeq x_{1i} \in R$, $x_2(t_i) \simeq x_{2i} \in R$, $y_1(t_i) \simeq y_{1i} \in R$, $y_2(t_i) \simeq y_{2i} \in R$, $u(t_i) \simeq u_i$, $u_1(t_i) \simeq u_{1i}$ and $u_2(t_i) \simeq u_{2i}$, the treatment function $f_{1a}$ becomes $\sum_{j=1}^{k_3} u_{1_{i-j}}\tau$ and the delay control problem (1)+(2) is transformed into the nonlinear programming problem (NLP)

$$\text{Minimize } J = g(x_N, y_N) + z_N \quad (3)$$

subject to

$$\left\{
\begin{array}{l}
x_{1i} - x_{1i+1} + \tau f_1(x_{1i}, x_{2i}, y_{1i}, y_{2i}, x_{1i-k_1}, x_{2i-k_1}, y_{1i-k_1}, y_{2i-k_1}) = 0 \\
x_{2i} - x_{2i+1} + \tau f_2(x_{2i}, x_{1i-k_2}, x_{2i-k_2}, y_{2i-k_2}) = 0 \\
y_{1i} - y_{1i+1} + \tau f_3(x_{1i}, x_{2i}, y_{1i}, y_{2i}, x_{1i-k_3}, x_{2i-k_3}, y_{1i-k_3}, y_{2i-k_3}, u_i, u_{1i}, u_{1i-1}, \dots, u_{1i-k_3}, u_{i-k_3}) = 0 \\
y_{2i} - y_{2i+1} + \tau f_4(y_{2i}, x_{2i-k_4}, y_{1i-k_4}, y_{2i-k_4}, u_{2i}, u_{i-k_4}) = 0 \\
z_i - z_{i+1} + \tau (u_i^2 + u_{1i}^2 + u_{2i}^2) = 0
\end{array}
\right.
\quad (4)$$

$$-u_i \le 0, u_i - 1 \le 0, \quad (5)$$

$$-u_{1i} \le 0, u_{1i} - 1 \le 0,$$

$$-u_{2i} \le 0, u_{2i} - 1 \le 0.$$
$$i = \overline{0, N - 1}.$$

Herein, the initial value profile $\varphi_1, \varphi_2, \varphi_3$ and $\varphi_4$ gives the values

$$x_{1-i} := \varphi_1(-i\tau), \quad i = \overline{0, k_1}$$

$$x_{2-i} := \varphi_2(-i\tau), \quad i = \overline{0, k_2}$$

$$y_{1-j} := \varphi_3(-j\tau), \quad j = \overline{0, k_3}$$

$$y_{2-j} := \varphi_4(-j\tau), \quad j = \overline{0, k_4}$$

The variable to be optimized is represented by the vector

$$w = (u_0, u_{1_0}, u_{2_0}, x_{1_1}, x_{2_1}, y_{1_1}, y_{2_1}, z_1, \dots, u_{N-1}, u_{1_{N-1}}, u_{2_{N-1}}, x_{1_N}, x_{2_N}, y_{1_N}, y_{2_N}, z_N) \in \mathbb{R}^{8N}.$$

The numerical procedure described above is applied in the next section and the graphs for states and control for different sets of parameters are obtained.

## 4 Numerical results and simulations

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 control of treatment. To solve the problem of optimal control the Matlab solver for NLP problems fmincon was used, selecting the 'interior-point' solver.

In all figures, for the healthy cell populations, we choose the same set of parameters value: $\eta_{1h} = 0.7, \eta_{2h} = 0.1, \tau_{1h} = 2, \tau_{2h} = 4, \gamma_{1h} = 0.1, \gamma_{2h} = 2.4, A_h = 922, \beta_{0h} = 1.77, k_{0h} = 0.1, \theta_{1h} = 0.5 \cdot 10^6, \theta_{2h} = 0.36 \cdot 10^8$. For leukemic cell populations, we consider alteration of the value of the following parameters:

- smaller percent of asymmetric division ($\eta_{1l}$);

- bigger percent of self-renewal ($1 - \eta_{1l} - \eta_{2l}$);

- lower rate of apoptosis of leukemic stem cells ($\gamma_{1l}$);

- lower rate of apoptosis of leukemic mature cells($\gamma_{2l}$);

- enhanced differentiation ($A_l$).

In the following example (Figures 1, 2) all these features were modified.

If we maintain the configuration of parameters from the previous example but consider that the percentage of self renewal of leukemic cells is the same as the percentage of self renewal of healthy cells, we see another manifestation of the disease, for which the treatment dose is lower than in the previous case (Figures 3, 4).

## Discussion

The plots of optimal controls (Figures 2 and 4) exhibit an optimal control effect almost constant at 0.5 respectively at 0.25 until the 90th day for all controls. One can intuitively expect that the treatment effect is proportional with the
States: T=100

1. Healthy cells parameters: $\eta_{1h}=0.70, \eta_{2h}=0.10, \tau_{1h}=2.00, \tau_{2h}=4.00, \gamma_{1h}=0.10, \gamma_{1h}=2.40,$
$k_{0h}=0.10, A_h=922$

2. Leukemic cells parameters: $\eta_{1l}=0.10, \eta_{2l}=0.50, \tau_{1l}=2.00, \tau_{2l}=4.00, \gamma_{1l}=0.03, \gamma_{2l}=0.80,$
$k_{0l}=0.10, A_l=1843$

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Comparison between the dynamics of a system without treatment and with optimal control of treatment. The results correspond to the best local minimal solution obtained for different initial guesses (the value of cost functional was improved from 1597 to 901)</figcaption></figure>

Control: T=100

1. Healthy cells parameters: $\eta_{1h}=0.70, \eta_{2h}=0.10, \tau_{1h}=2.00, \tau_{2h}=4.00, \gamma_{1h}=0.10, \gamma_{2h}=2.40,$
$k_{0h}=0.10, A_h=922$

2. Leukemic cells parameters: $\eta_{1l}=0.10, \eta_{2l}=0.50, \tau_{1l}=2.00, \tau_{2l}=4.00, \gamma_{1l}=0.03, \gamma_{2l}=0.80,$
$k_{0l}=0.10, A_l=1843$

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Optimal control of treatment. The controls u, u1, u2 represent the influence of drug on the proliferation rate and apoptosis. One can observe that the drug influence is almost constant and the evolution of u, u1 and u2 are similar if parameters a and b of cost functional are 1.</figcaption></figure>
States: T=100

1. Healthy cells parameters: $\eta_{1h}=0.70, \eta_{2h}=0.10, \tau_{1h}=2.00, \tau_{2h}=4.00, \gamma_{1h}=0.10, \gamma_{1h}=2.40,$
$k_{0h}=0.10, A_h=922$

2. Leukemic cells parameters: $\eta_{1l}=0.40, \eta_{2l}=0.40, \tau_{1l}=2.00, \tau_{2l}=4.00, \gamma_{1l}=0.03, \gamma_{2l}=0.80,$
$k_{0l}=0.10, A_l=1843$

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Comparison between the dynamics of a system without treatment and with optimal control of treatment. The results correspond to the best local minimal solution obtained for different initial guesses (the value of cost functional was improved from 641 to 260)</figcaption></figure>

Control: T=100

1. Healthy cells parameters: $\eta_{1h}=0.70, \eta_{2h}=0.10, \tau_{1h}=2.00, \tau_{2h}=4.00, \gamma_{1h}=0.10, \gamma_{2h}=2.40,$
$k_{0h}=0.10, A_h=922$

2. Leukemic cells parameters: $\eta_{1l}=0.40, \eta_{2l}=0.40, \tau_{1l}=2.00, \tau_{2l}=4.00, \gamma_{1l}=0.03, \gamma_{2l}=0.80,$
$k_{0l}=0.10, A_l=1843$

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Optimal control of treatment. The controls u, u1, u2 represent the influence of drug on the proliferation rate and apoptosis. One can observe that the drug influence is almost constant and the evolution of u, u1 and u2 are similar.</figcaption></figure>
cell concentration of imatinib; however, in order to study the influence of the
prescribed drug concentration on the population's time-evolution, imatinib phar-
macokinetics (PK) and pharmacodynamics (PD) need to be taken into account
(see also [20]). Nevertheless, for an optimal effect of treatment, the prescribed
dose should be adapted in view of the disease parameters of a certain patient.
In that respect, a competition model of healthy vs. CML cell population that
takes into account the influence of PK and PD of imatinib, is subject of further
research.

References

1. Adimy M., Crauste F., Halanay A., Neamțu M., Opriș D.: Stability of Limit Cycles in a Pluripotent Stem Cell Dynamics Model. Chaos. Solitons& Fractals 27(4), 1091-1107 (2006)

2. Adimy M., Crauste F., Ruan S.: A mathematical study of the hematopoiesis process with application to chronic myelogenous leukemia. SIAM J. Appl. Math. 65(4), 1328-1352 (2005)

3. Ainseba B., Benosman C.: Optimal control for resistance and suboptimal response in CML. Mathematical Biosciences 227, 81-93 (2010)

4. Beckman J., Scheitza S., Wernet P., Fischer J., Giebel B.: Asymmetric cell division within the human hematopoietic stem and progenitor cell compartment: identification of asymmetrically segregating proteins. Blood 12, 109, 5494-5501 (2007)

5. Benosman C.: Control of the dynamics of Chronic Myeloid Leukemia by Imatinib. Ph.D Thesis (2010)

6. Colijn C., Mackey M.C.: A mathematical model of hematopoiesis - I. Periodic chronic myelogenous leukemia. J. Theor. Biol. 237, 117-132 (2005)

7. Deiniger MWN., Goldman JM, Lydon N., Melo JV.: The tyrosine kinase inhibitor CGP57148B selectively inhibits the growth of BCRABL positive cells. Blood 90, 3691-3698 (1997)

8. Gollmann L., Kern D., Maurer H.: Optimal control problems with control and state delays and applications to growth processes, IIASA Symposium on "Applications of Dynamic Systems to Economic Growth with Environment", Luxemburg, Nov. 7-8, (2008)

9. Gollmann L., Kern D., Maurer H.: Optimal control problems with delays in state and control variables subject to mixed control state constraints. Optim. Control Appl. Meth. 30, 341-365 (2009)

10. Gollmann L., Maurer H.: Theory and applications of optimal control problems with multiple time-delays. Journal of Industrial and Management Optimization 10,2 , 413-441 (2014)

11. Gottschalk S., Anderson N., Hainz C. et al.: Imatinib (STI571)-Mediated Changes in Glucose Metabolism in Human Leukemia BCR-ABL-Positive Cells. Clin Cancer Res 10, 6661-6668 (2004)

12. Halanay A.: Periodic Solutions in Mathematical Models for the Treatment of Chronic Myelogenous Leukemia. Math. Model. Nat. Phenom 7, No.1, 235-244 (2012)

13. Halanay A., Candea D., Radulescu I.R.: Stability analysis of equilibria in a delay differential equations model of CML including asymmetric division and treatment. Math. Comput. Simul., Elsevier, to appear (2014)

14. Kim P., Lee P., Levy D.: Dynamics and Potential Impact of the Immune Response to Chronic Myelogenous Leukemia. PLoS Comput.Biol. 4(6):e 1000095 (2008)
15. Kim P., Lee P., Levy D.: A theory of immunodominance and adaptive regulation. Bull. Math. Biol. 73(7), 1645-65 (2011)

16. Mackey M.C., Ou C., Pujo-Menjouet L., Wu J.: Periodic oscillations of blood cell population in chronic myelogenous leukemia. SIAM J. Math. Anal. 38, 166-187 (2006)

17. Marciniak-Czochra A., Stiehl T., Wagner W.: Modeling of replicative senescence in hematopoietic development. Aging 1(8), 723-732 (2009)

18. Michor F., Hughes T., Iwasa Y., Branford S., Shah N.P., Sawyers C., Novak M.: Dynamics of chronic myeloid leukemia. Nature (2005), 435, 1267-1270

19. Moore H., Li N.K.: A mathematical model for chronic myelogenous leukemia (CML) and T-cell interaction. J. Theor. Biol. 227, 513-523 (2004)

20. Radulescu I.R., Candea D., Halanay A.: Stability and bifurcation in a model for the dynamics of stem-like cells in leukemia under treatment. American Institute of Physics Proceedings, 1493, 758-763 (2012)

21. Reya T.: Regulation of Hematopoietic Stem Cell Self-Renewal. Recent Progress in Hormone Research, 58, 283-295 (2003)

22. Stiehl T., Marciniak-Czochra A.: Mathematical Modeling of Leukemogenesis and Cancer Stem Cell Dynamics. Math. Model. Nat. Phenom. 7, 1, 166-202 (2012)

23. Tang M., Foo J., Gonen M., Mahon F-X., Michor F.: Selection pressure exerted by imatinib therapy leads to disparate outcomes of imatinib discontinuation trials, Haematologica 97(10), 1553-1561 (2012)

24. Tomasetti C., Levi D.: Role of symmetric and asymmetric division of stem cells in developing drug resistance, PNAS 17 (39), 16766-16771 (2010)