# Two numerical methods for solving a nonlinear system of integral equations of mixed Volterra-Fredholm type arising from a control problem related to leukemia

## Abstract

The aim of this paper is to present two algorithms for numerical solving of a fixed final state control problem in connection with the leukemia treatment strategy. In the absence of the controllability condition, our model leads to a nonlinear integral system of Volterra type to whom explicit iterative techniques apply and converge. Once using the controllability condition, the control variable is expressed in terms of the state variables and the integral system changes to a mixed Volterra-Fredholm type one making direct iterative techniques inoperative. However, two paths can be followed. One consists in an iterative procedure where at each step the control variable is calculated using the approximate values of the state variables from the previous step. The other looks for the numerical value of the control variable by using the bisection method. Numerical simulations, error analysis and biological interpretation are given.

## Authors

Lorand Gabriel Parajdi
Department of Mathematics, West Virginia University, Morgantown, WV 26506, USA & Department of Mathematics, Babeş-Bolyai University, Cluj-Napoca, Romania

Flavius Pătrulescu
Department of Mathematics, Technical University of Cluj-Napoca, Cluj-Napoca, Romania;

Faculty of Mathematics and Computer Science and Institute of Advanced Studies in Science and Technology, Babeş-Bolyai University, Cluj-Napoca 400084, Romania
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, Cluj-Napoca 400110, Romania

Ioan Ştefan Haplea
Department of Internal Medicine, Iuliu Haţieganu University of Medicine and Pharmacy, Cluj-Napoca, Romania

## Keywords

Control problem; myeloid leukemia; Volterra-Fredholm integral equation; numerical method.

## Paper coordinates

L.G. Parajdi, F. Pătrulescu, R. Precup, I. Şt. Haplea, Two numerical methods for solving a nonlinear system of integral equations of mixed Voltera-Fredholm type arising from a control problem related to leukemia, Journal of Applied Analysis & Computation, 13 (2023) no. 4, pp. 1797-1812,  http://doi.org/10.11948/20220197

## PDF

##### Journal

Journal of Applied Analysis & Computation

Springer

## Paper (preprint) in HTML form

Two numerical methods for solving a nonlinear system of integral equations of mixed Volterra-Fredholm type arising from a control problem related to leukemia

# Two numerical methods for solving a nonlinear system of integral equations of mixed Volterra-Fredholm type arising from a control problem related to leukemia

Lorand Gabriel Parajdi, Flavius Olimpiu Patrulescu, Radu Precup, Ioan Ştefan Haplea L.G. Parajdi, Department of Mathematics, West Virginia University, Morgantown, WV 26506, USA & Department of Mathematics, Babeş–Bolyai University, 400084 Cluj-Napoca, Romania F.O. Patrulescu, Department of Mathematics, Technical University of Cluj-Napoca, 28 Memorandumului Street, 400114 Cluj-Napoca, Romania R. Precup, Faculty of Mathematics and Computer Science and Institute of Advanced Studies in Science and Technology, Babeş-Bolyai University, 400084 Cluj-Napoca, Romania & Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, P.O. Box 68-1, 400110 Cluj-Napoca, Romania I.Ş. Haplea, Department of Internal Medicine, Iuliu Haţieganu University of Medicine and Pharmacy, 400012 Cluj-Napoca, Romania
###### Abstract.

The aim of this paper is to present two algorithms for numerical solving of a fixed final state control problem in connection with the leukemia treatment strategy. In the absence of the controllability condition, our model leads to a nonlinear integral system of Volterra type to whom explicit iterative techniques apply and converge. Once using the controllability condition, the control variable is expressed in terms of the state variables and the integral system changes to a mixed Volterra-Fredholm type one making direct iterative techniques inoperative. However, two paths can be followed. One consists in an iterative procedure where at each step the control variable is calculated using the approximate values of the state variables from the previous step. The other looks for the numerical value of the control variable by using the bisection method. Numerical simulations, error analysis and biological interpretation are given.

###### Key words and phrases:
Control problem, dynamic system, myeloid leukemia, Volterra-Fredholm integral equation, numerical method
###### 1991 Mathematics Subject Classification:
45G15, 34H05, 65R20, 92C50, 93B05

## 1. Introduction

Human activity in various fields often involves some control over the processes in order to achieve the desired result. Thus, in particular, control intervenes in medicine and is exercised in order to bring the patient to the desired parameters.

The control becomes a mathematical one as soon as a mathematical model is generated for the investigated process by taking into account specific laws such as conservation laws. From a mathematical point of view, the control involves a certain modification of some parameters of the model such that the solution of the problem satisfies a certain requirement.

For such a problem, it is important to demonstrate qualitatively the controllability of the system, i.e., the possibility of achieving the desired goal, but just as important for the implementation of theoretical predictions are the numerical methods designed to provide quantitative results. Therefore, numerical approximation and calculation algorithms are especially important for control problems.

The purpose of this paper is precisely to develop computer-based numerical algorithms for the numerical solution of a control problem related to the treatment of chronic myeloid leukemia. Even if they are presented for a very specific problem, these algorithms can be easily adapted for the numerical treatment of many other classes of control problems.

### 1.1. Biological Background

Chronic myeloid leukemia (CML) is a slowly progressing malignancy of the blood marrow, derived from the granulocyte cell line. Left untreated, it causes moderate symptoms on its own, but can nonetheless progress over years to aggressive, life-threatening forms (accelerated or blast phase) (see E. Jabbour and H. Kantarjian [1]). The molecular culprit (and hallmark) of CML is the oncoprotein BCR-ABL1, generated by an abnormal fusion of two genes, which normally belong to distinct chromosomes. BCR-ABL1 has intrinsic tyrosine kinase activity and switches the molecular machinery of cell signal transduction to an ”always on” state. This increases cell division and proliferation, largely irrespective of stimulation by growth factors and confers a malignant phenotype (see S. Faderl et al. [2]).

Frontline therapy for CML relies on tyrosine kinase inhibitors (TKI). These are targeted anticancer drugs, designed to block the BCR-ABL1 protein, specifically, potent and selective enough to drive the leukemic cell population closer to extinction while sparing the normal granulocyte line. From the first molecule (Imatinib, see S.G. O’Brien et al. [3]), TKIs have diversified to date to second and third generation agents, with increased potency. All TKIs are very effective in treating CML, with rates of major response at 10 years approaching 95% (see A. Hochhaus et al. [4]).

The state of the disease and the response to treatment are monitored primarily by the number of mARN copies of the BCR-ABL1 gene, as measured by qRT-PCR from the patient’s blood, normalized to a control (see T. Hughes et al. [5]). Presently, the life expectancy of patients with CML approaches the average of the general population, especially in developed countries (see M.W. Deininger et al. [6]). In addition to survival, the focus of research in CML is quality of life, the effect of comorbidities, toxicity profile of TKI, and the discontinuation of drugs once major treatment milestones have been met (see A. Hochhaus et al. [7]). These concerns are especially prominent since the median age of CML onset is close to 70 years of age. The side effects of TKIs can be mitigated by dose reduction (for Imatinib, from 400 mg per day to 300 mg per day) or a brief drug holiday.

All these suggest the necessity for an optimized and theoretically grounded dosing schedule in CML that reaches the treatment targets in an efficient yet parsimonious way. We have explored the problem of TKI dosing in CML by a quantitative dynamical model.

### 1.2. Mathematical Model and Approach

Over the years, many mathematicians and biologists have proposed and studied different types of models and problems inspired by life phenomena, such as Y. Cherruault [8], C.J.S. Alves et al. [9], Y. Wang et al. [10] and G. Cedersund et al. [11]. From the numerous works in which there are studied optimization problems and control problems applied to chronic myeloid leukemia (also known as chronic myelogenous leukemia), we mention those of S. Nanda et al. [12], B. Ainseba and C. Benosman [13], Q. He et al. [14], S.B. Mendrazitsky and B. Shklyar [15] and S.B. Mendrazitsky et al. [16].

In a paper of M.C. Mackey and L. Glass [17], it is introduced a mathematical model of blood production using differential equations involving sigmoid or Hill functions, given the fact that hematopoiesis is a self-limiting process. Based on this model of blood production from Mackey-Glass [17], in L.G. Parajdi et al. [18] it is considered the following mathematical model

 (1.1) $\begin{cases}x^{\prime}(t)=\frac{ax(t)}{1+b_{1}x(t)+b_{2}y(t)}-cx(t)\\[5.69054% pt] y^{\prime}(t)=\frac{Ay(t)}{1+B\left(x(t)+y(t)\right)}-Cy(t)\end{cases}$

where $x(t),\ y(t)$ are the normal and leukemic cell populations at time $t$, respectively. Here, the parameters $a,A$ are the growth rates; $c,C$ are the cell death rates (or cell turnover rates); and $b_{1},b_{2},B$ are the sensitivity parameters that govern the self-limiting process. For both cell populations, we assume that the growth rate is greater than the death rate, i.e., $a>c$ and $A>C$. The eventual advantage of leukemic cells of being less sensitive to the microenvironment than normal cells is expressed by

 $b_{1}\geq b_{2}>B.$

The case $\ b_{1}=b_{2}$, was studied by Dingli and Michor [19] (see also the paper of A. Cucuianu and R. Precup [20]), in order to describe the time competition between normal and leukemic hematopoietic stem cells. In this case, there exist only two non-zero steady states $(d,0)$ and $(0,D)$ of the system, where $d$ and $D$ represent the homeostatic amounts of normal and leukemic (or abnormal) cells and are given by

 $d=\frac{1}{b_{1}}\left(\frac{a}{c}-1\right)\ \ \ \text{and}\ \ \ D=\frac{1}{B}% \left(\frac{A}{C}-1\right).$

The case $b_{1}>b_{2}$, was introduced and studied in the paper of L.G. Parajdi et al. [18] and the book of L.G. Parajdi [21]. In this case, besides the two non-zero steady states $(d,0)$, $(0,D)$ there could also exist a third steady state $(x^{\ast},y^{\ast})$, where

 $x^{\ast}=\frac{b_{1}}{b_{1}-b_{2}}d-\frac{b_{2}}{b_{1}-b_{2}}D\ \ \ \text{and}% \ \ \ y^{\ast}=\frac{b_{1}}{b_{1}-b_{2}}(D-d).$

We can observe that under the assumption $b_{1}>b_{2}$, both $x^{\ast}$ and $y^{\ast}$ are positive if and only if $d The case when $D corresponds to the normal hematologic state, the case $d corresponds to the chronic phase of the disease, and the case $(b_{1}/b_{2})d stands for the acute/blast phase. In the first case (normal state), system (1.1) has the unique stable equilibrium $(d,0)$, in the second case (the chronic phase) the stable equilibrium is of the form $(x^{\ast},y^{\ast})$, where both $x^{\ast}$ and $y^{\ast}$ are positive, and in the third case (the acute/blast phase) the stable equilibrium is $(0,D)$.

In this paper, we give two numerical methods for solving a nonlinear system of integral equations of mixed Volterra–Fredholm type, arising from a control problem related to myeloid leukemia. This system of integral equations, presented in Section 2, is related to the control problem where the control is a constant $\lambda$ (connected to a constant drug dose) intended to decrease the growth rate $A$ of malignant cells. The model could be related to targeted therapies. This control problem was introduced and studied in the paper of I.Ş. Haplea et al. [22] and its analysis, reproduced in Section 2, is based on a new method for the controllability of abstract fixed point equations. The stability of the control problem is also reminded. In Section 3, we introduce and describe two algorithms that are used to approximate numerically the solution of the integral system. In this context, by Theorem 3.2, it is proved the continuous dependence of the solution on the control parameter, and by Theorem 3.3, the convergence of our first algorithm. Numerical simulations and error analysis are then performed to illustrate the theoretical results and prove their applicability. Finally, in Section 4, we discuss the biological significance of our theoretical and numerical results.

## 2. A Control Problem for the Normal-Leukemic System

The control problem that is considered in this paper is inspired from hematology, more exactly from the treatment of the chronic phase of chronic myeloid leukemia (CML). It was introduced and studied in I.Ş. Haplea et al. [22] and directly folds on the general control problem for differential systems as stated in V. Barbu [23] (page 34). In this section, we recall the general problem of control theory, the general controllability principle, and a result about the stability of the control problem. All these notions and results are stated in terms of fixed point theory.

### 2.1. A General Controllability Principle

Consider a general control problem that consists in finding a solution $\left(w,\lambda\right)$ to the following system

 (2.1) $\begin{cases}w=N\left(w,\lambda\right)\\[5.69054pt] w\in W,\ \lambda\in\Lambda,\ \left(w,\lambda\right)\in D\end{cases}$

where $w$ is the state variable, $\lambda$ denotes the control variable, $W$ is the domain of the states, $\Lambda$ represents the domain of controls and finally, $D$ is the controllability domain, usually given by means of some condition imposed to $w,$ or to both $w$ and $\lambda.$ Therefore we take into consideration a general formulation of the control problem, only in terms of sets, since $W,\ \Lambda$ and $D\subset W\times\Lambda$ are not necessarily structured sets and $N$ is any mapping from $W\times\Lambda$ to $W$.

We say that the equation $w=N(w,\lambda)$ is controllable in $W\times\Lambda$ with respect to $D$, if problem (2.1) has a solution.

Let $\Sigma$ be the set of all possible solutions $\left(w,\lambda\right)$ of the fixed point equation and $\Sigma_{1}$ be the set of all $w$ that are the first components of some solutions of the fixed point equation, that is

 $\displaystyle\Sigma$ $\displaystyle=$ $\displaystyle\left\{\left(w,\lambda\right)\in W\times\Lambda:\ w=N\left(w,% \lambda\right)\right\},$ $\displaystyle\Sigma_{1}$ $\displaystyle=$ $\displaystyle\left\{w\in W:\exists\,\lambda\in\Lambda\text{ with }(w,\lambda)\in\Sigma\right\}.$

Obviously, the set of all solutions of the control problem (2.1) is given by $\ \Sigma\cap D.$

Consider the set-valued map $\ F:\Sigma_{1}\rightarrow\Lambda$  defined as

 $F(w)=\{\lambda\in\Lambda:\ \text{ }w=N(w,\lambda)\text{ and }(w,\lambda)\in D\}.$

We have the following general principle for solving the control problem (2.1).

###### Proposition 2.1.

If for some extension $\ \widetilde{F}:W\rightarrow\Lambda$ of $F$ from $\Sigma_{1}$ to $W,$ there exists a fixed point $\ w\in W$ of the set-valued map

 $\widetilde{N}\left(w\right):=N\left(w,\widetilde{F}\left(w\right)\right),$

i.e.,

 $w=N\left(w,\lambda\right),\$

for some $\ \lambda\in\widetilde{F}\left(w\right),$ then the couple $\left(w,\lambda\right)$ is a solution of the control problem (2.1).

For the proof of Proposition 2.1, see the paper of I.Ş. Haplea et al. [22].

###### Remark 2.1.

In particular, $F$ and $\widetilde{F}$ can be single-valued maps.

The above general principle of controllability can be applied to control problems related to the normal-leukemic system.

### 2.2. A Control Problem for the Normal-Leukemic System

Consider as a condition of controllability, the decrease to a certain acceptable level of the malignant cell population and as a control variable, the factor of decreasing the proliferation rate of cancer cells. Therefore, we consider the following initial value problem

 (2.2) $\begin{cases}x^{\prime}(t)=\frac{ax(t)}{1+b_{1}x(t)+b_{2}y(t)}-cx(t)\\[5.69054% pt] y^{\prime}(t)=\frac{\lambda Ay(t)}{1+B\left(x(t)+y(t)\right)}-Cy(t)\\[5.69054% pt] x(0)=x_{0},\ y(0)=y_{0}\end{cases}$

where $\lambda>0$ represents the control parameter.

Let us change the variables as follows:

 $u:=ln\ x\ \ \text{and}\ \ v:=ln\ y$

and denote

 $u_{0}=u(0)=ln\ x(0)\ \ \text{and}\ \ v_{0}=v(0)=ln\ y(0).$

Dividing the first equation by $x$ and the second equation by $y$ in the system (2.2), and then integrating, we obtain the integral system equivalent to the initial value problem (2.2), namely

 (2.3) $\begin{cases}u\left(t\right)=u_{0}-ct+\displaystyle\int_{0}^{t}\frac{a}{1+b_{1% }e^{u(s)}+b_{2}e^{v(s)}}\,ds\\[5.69054pt] v\left(t\right)=v_{0}-Ct+\lambda\displaystyle\int_{0}^{t}\frac{A}{1+B\left(e^{% u(s)}+e^{v(s)}\right)}\,ds.\end{cases}$

This is our fixed point equation $w=N(w,\lambda),$ where $w=(u,v),$ $N=(N_{1},N_{2})$ and

 $\displaystyle N_{1}\left(u,v\right)\left(t\right)$ $\displaystyle=$ $\displaystyle u_{0}-ct+\int_{0}^{t}\frac{a}{1+b_{1}e^{u(s)}+b_{2}e^{v(s)}}ds,$ $\displaystyle N_{2}\left(u,v,{\color[rgb]{1,0,0}\lambda}\right)\left(t\right)$ $\displaystyle=$ $\displaystyle v_{0}-Ct+\lambda\int_{0}^{t}\frac{A}{1+B\left(e^{u(s)}+e^{v(s)}% \right)}\,ds.$

The objective condition associated to the integral system (2.3) is

 (2.4) $v(T)=v_{T},$

where $v_{T}$ is the expected level of leukemic cells after a period of time $T$.

The biological interpretation is the following: the system (2.2) extends the basic system (1.1) by adding the effect of the treatment. We have assumed that TKIs, being a targeted drug, acts exclusively on the malignant cell (absolute specificity), and solely by decreasing the proliferation rate of malignant cells (from $A$ to $\lambda A$, where $\lambda$ is a positive subunit quantity, i.e., $0<\lambda<1$) (see F. Michor et al. [24]). A smaller value for $\lambda$ translates to a more efficient treatment (a higher dose or a more potent drug). The parameter $\lambda$ stands for the effect of TKI on the growth rate of malignant cells, rather than the actual drug dosage; $\lambda$ is completely determined by the drug dosage but relates to it in a nontrivial way. As a simplification, we have assumed a dosing scheme that maintains the drug effect $\lambda$ constant, for the whole duration of the treatment; as such, all pharmacokinetic details could be safely excluded from the analysis. The goal of the optimization procedure is to reach the desired (very low) number of malignant cells, after a given treatment time $T$, $v(T)=v_{T}$. This is qualitatively in line with current clinical guidelines, which prescribe that BCR-ABL1 transcript levels should be reduced strongly after 12 months of treatment (major molecular response, MMR) and, in the long run, be brought next to or below detection levels (deep molecular remission, DMR) (see the ESMO Clinical Practice Guidelines: A. Hochhaus et al. [25]).

Compared with the abstract general controllability principle from the previous subsection, here we have

 $W=C\left([0,T],\mathbb{R}^{2}\right),\ \Lambda=\mathbb{R}_{+},\ D=\{(u,v,% \lambda)\in C\left([0,T],\mathbb{R}^{2}\right)\times\mathbb{R}_{+}:v(T)=v_{T}\}$

and the integral system (2.3) stands for the fixed point equation $w=N(w,\lambda)$, with $w=(u,v)$.

We recall the following controllability result obtained using the general principle given by Proposition 2.1 from I.Ş. Haplea et al. [22]. This result shows the solvability of the control problem (2.3)-(2.4).

###### Theorem 2.2.

For each number  $v_{T}$ with  $\max\left\{0,v_{0}-CT\right\} the control problem (2.3)-(2.4) has a solution $\ \left(u,v,\lambda\right)\in C\left(\left[0,T\right],\mathbb{R}^{2}\right)% \times\left(0,+\infty\right).$

###### Proof.

Let $(u,v)\in\Sigma_{1}$, that is $(u,v)\in C\left(\left[0,T\right],\mathbb{R}^{2}\right)$ which solves (2.3) for some $\lambda\in\mathbb{R}_{+}$. From the second equation of (2.3) and using the controllability condition (2.4), we obtain

 (2.5) $\lambda=\frac{v_{T}-v_{0}+CT}{\displaystyle\int_{0}^{T}\frac{A}{1+B(e^{u(s)}+e% ^{v(s)})}\,ds}.$

We consider a function $\ F:\Sigma_{1}\rightarrow\mathbb{R}_{+}$ given by

 $F\left(u,v\right)=\frac{1}{\displaystyle\int_{0}^{T}\frac{A}{1+B(e^{u(s)}+e^{v% (s)})}ds}\left(v_{T}-v_{0}+CT\right).$

Next, we can extend this function $F$ from $\Sigma_{1}$ to $C\left(\left[0,T\right],\mathbb{R}^{2}\right)$ by using the same expression, namely

 $\widetilde{F}\left(u,v\right)=\frac{1}{\displaystyle\int_{0}^{T}\frac{A}{1+B(e% ^{u(s)}+e^{v(s)})}ds}\left(v_{T}-v_{0}+CT\right)$

for $\left(u,v\right)\in C\left(\left[0,T\right],\mathbb{R}^{2}\right)$. We note that $F$ and $\widetilde{F}$ are single-valued maps.

Thus, our fixed point equation $w=\widetilde{N}(w),$ with $\widetilde{N}(w)=N(w,\widetilde{F}(w))$, becomes

 (2.6) $\left\{\begin{array}[]{l}u\left(t\right)=u_{0}-ct+\displaystyle\int_{0}^{t}% \frac{a}{1+b_{1}e^{u(s)}+b_{2}e^{v(s)}}ds\\[8.0pt] v\left(t\right)=v_{0}-Ct+\frac{v_{T}-v_{0}+CT}{\displaystyle\int_{0}^{T}\frac{% A}{1+B(e^{u(s)}+e^{v(s)})}\,ds}\displaystyle\int_{0}^{t}\frac{A}{1+B\left(e^{u% (s)}+e^{v(s)}\right)}ds,\end{array}\right.$

which is a system of mixed Volterra–Fredholm type integral equations, where

 $\displaystyle\widetilde{N}_{1}\left(u,v\right)\left(t\right)$ $\displaystyle=$ $\displaystyle u_{0}-ct+\int_{0}^{t}\frac{a}{1+b_{1}e^{u(s)}+b_{2}e^{v(s)}}ds,$ $\displaystyle\widetilde{N}_{2}\left(u,v\right)\left(t\right)$ $\displaystyle=$ $\displaystyle v_{0}-Ct+\frac{v_{T}-v_{0}+CT}{\displaystyle\int_{0}^{T}\frac{A}% {1+B(e^{u(s)}+e^{v(s)})}\,ds}\int_{0}^{t}\frac{A}{1+B\left(e^{u(s)}+e^{v(s)}% \right)}ds.$

The existence of a fixed point $w$ of $\widetilde{N}$ can be guaranteed by using Schauder’s fixed point theorem (see R. Precup [26], page 33). Indeed, the operator $\widetilde{N}$ is completely continuous, as follows immediately from the Arzelà–Ascoli theorem (see R. Precup [26], pages 2, 3 and 15). On the other hand, we have

 $\displaystyle u_{0}-cT$ $\displaystyle\leq$ $\displaystyle\widetilde{N}_{1}\left(u,v\right)\left(t\right)\leq u_{0}+aT,$ $\displaystyle v_{0}-CT$ $\displaystyle\leq$ $\displaystyle\widetilde{N}_{2}\left(u,v\right)\left(t\right)\leq v_{T}+CT.$

Consequently, there are numbers $R_{1},R_{2}>0$ with

 $\left\|\widetilde{N}_{1}\left(u,v\right)\right\|_{\infty}\leq R_{1},\ \ \ % \left\|\widetilde{N}_{2}\left(u,v\right)\right\|_{\infty}\leq R_{2}$

for all $u,v\in C\left[0,T\right].$ Hence, $\widetilde{N}\left(S\right)\subset S,$ where

 $S:=\{\left(u,v\right)\in C\left(\left[0,T\right],\mathbb{R}^{2}\right):\left\|% u\right\|_{\infty}\leq R_{1},\ \left\|v\right\|_{\infty}\leq R_{2}\}.$

Thus, Schauder’s fixed point theorem applies and guarantees the existence in $S$ of a fixed point $w=\left(u,v\right)\$of $\widetilde{N}.$ Finally, the control $\lambda$ is calculated using formula (2.5) with $u$ and $v$ thus determined.

###### Remark 2.2.

Since (2.3) is a Volterra system with Lipschitz continuous nonlinearities, for any given $\lambda$, it has a unique solution $\left(u,v\right).$ In particular, for the value $\lambda$ corresponding to a solution $\left(u,v,\lambda\right)$ of the control problem, the trajectory $\left(u,v\right)$ is unique. We may interpret this fact in the following way: for a prescribed drug dose, the patient’s evolution is uniquely determined.

### 2.3. Stability of the Control Problem

Let us also mention the following stability result of the control problem given in I.Ş. Haplea et al. [22].

###### Theorem 2.3.

Let $\ \left(u,v,\lambda\right)$ be a solution of the control problem (2.3)-(2.4). For a given $\ \overline{\lambda}\in\mathbb{R}_{+},$ denote by $\ \left(\overline{u},\overline{v}\right)$ the corresponding solution of system (2.3) and $\ \overline{v}_{T}:=\overline{v}\left(T\right).$ Then, for any small $\ \varepsilon>0,$ we have that

 $\left|\overline{\lambda}-\lambda\right|\leq\frac{\varepsilon}{AT}e^{-\left(a+% \lambda A\right)T}=:\delta\ \ \ \text{implies\ \ \ }\left|\overline{v}_{T}-v_{% T}\right|\leq\varepsilon.$

In its proof, there are used Lipschitz estimates and Gronwall’s inequality in order to obtain the estimate

 (2.7) $\left|u\left(t\right)-\overline{u}\left(t\right)\right|+\left|v\left(t\right)-% \overline{v}\left(t\right)\right|\leq AT\delta e^{\left(a+\lambda A\right)t},$

which in view of the expression of $\delta$, immediately implies

 $\ \left|\overline{v}_{T}-v_{T}\right|\leq\varepsilon.$
###### Remark 2.3.

Estimate (2.7) shows us that for a treatment $\overline{\lambda}$ close enough to the prescribed treatment $\lambda,$ the patient’s evolution $\left(\overline{u},\overline{v}\right)$ remains in the vicinity of the prescribed evolution $\left(u,v\right).$

Note that, due to this stability result, the actual administration of only an approximate dose leads, however, to a result close to the expected one. The level of freedom to choose an approximate dose is exactly established in terms of model parameters.

## 3. Numerical Algorithms

In this section, we introduce our numerical algorithms for solving problem (2.3)-(2.4) which can be put under the form of the Volterra-Fredholm integral system (2.6). At any step of each one of the two algorithms, we need to solve numerically a Volterra type integral system and this is done with the Picard-Lindel${\rm\ddot{o}}$f iteration technique (see E. Hairer et al. [27]) combined with the successive approximation method. The description of the algorithms is followed by some implementation details. Finally, some numerical results in the form of tables and figures are given, and their biological interpretation is the subject of the discussion section, Section 4. All numerical simulations were carried out on the Kotys HPC (High Performance Computing) infrastructure of ”Babeş–Bolyai” University of Cluj-Napoca (see D. Bufnea et al. [28]), using MATLAB software.

### 3.1. The First Algorithm

In this subsection, we introduce our first algorithm. Therefore, we prove the continuous dependence result, see Theorem 3.2, using the definition of the Bielecki norm and the definition of the spectral radius plus the lemma of convergence of the power sequence of a matrix, all mentioned below.

###### Definition 3.1.

Let $(C([0,T],\mathbb{R}),||\cdot||_{\theta})$ be a Banach space, where $f\in C([0,T],\mathbb{R})$, we have the Bielecki norm as follows:

 $||f||_{\theta}=\max_{t\in[0,T]}\left(|f(t)|e^{-\tau|t-t_{0}|}\right),\ \text{% where}\ \tau>0\ \text{and}\ t_{0}\in[0,T].$
###### Definition 3.2.

The spectral radius of a square matrix $M$ is

 $\rho(M)=\max\ \{|\lambda|:\lambda\ \text{is an eigenvalues of matrix M}\}.$

For an $n\times n$ matrix $M,$ let $||M||=\max\ \{|M_{ij}|:1\leq i,j\leq n\}$.

###### Lemma 3.1.

Let $M$ be a square matrix, with spectral radius $\rho(M).$ If $\rho(M)<1$, then $\displaystyle{\lim_{k\to\infty}M^{k}=0}$. On the other hand, if $\rho(M)>1$, then $\displaystyle{\lim_{k\to\infty}||M^{k}||=\infty}$. The statement holds for any choice of matrix norm.

###### Proof.

We recall that the matrix $M=CJC^{-1}$ for a matrix $J$ in the Jordan normal form and regular matrix $C$, and we know that $M^{k}=CJ^{k}C^{-1}$. Next, if $\rho(M)=\rho(J)<1$, then $J^{k}$ converges to the zero matrix, and thus $M^{k}$ converges to the zero matrix as well. If $\rho(M)>1$, then $J^{k}$ has a diagonal entry $(J^{k})_{ii}=\lambda^{k}$ for an eigenvalue $\lambda$ such that $|\lambda|>1$, and if $w$ is the $i$-th column of $C$ and $w^{\prime}$ the $i$-th row of $C^{-1}$, then $w^{\prime}M^{k}w=w^{\prime}CJ^{k}C^{-1}w=e^{T}_{i}J^{k}e_{i}=\lambda^{k}$. Therefore, $\displaystyle{\lim_{k\to\infty}|w^{\prime}M^{k}w|=\infty}$, and thus $\displaystyle{\lim_{k\to\infty}||M^{k}||=\infty}$. ∎

###### Theorem 3.2.

The solution $(u,v)$ of system (2.3) depends continuously on parameter $\lambda$.

###### Proof.

We denote by $S_{1}(\lambda)$ and $S_{2}(\lambda)$ the components $u$ and $v$ of the solution that corresponds to $\lambda$. From the integral system, using the Lipschitz continuity of the nonlinearities and the Volterra property of the equations, we deduce the estimates of the form

 $\|S_{1}(\lambda)-S_{1}(\mu)\|_{\theta}\leq\alpha_{11}\|S_{1}(\lambda)-S_{1}(% \mu)\|_{\theta}+\alpha_{12}\|S_{2}(\lambda)-S_{2}(\mu)\|_{\theta}$

and

 $\|S_{2}(\lambda)-S_{2}(\mu)\|_{\theta}\leq\alpha_{21}\|S_{1}(\lambda)-S_{1}(% \mu)\|_{\theta}+\alpha_{22}\|S_{2}(\lambda)-S_{2}(\mu)\|_{\theta}+\beta|% \lambda-\mu|.$

Here $\|\cdot\|_{\theta}$ denotes the Bielecki norm on $C[0,T]$ with respect to a sufficiently large $\theta>0$. Moreover, $\beta>0$ represents some constant and the spectral radius $\rho(M)$ of matrix $M:=[\alpha_{ij}]_{1\leq i,j\leq 2}$ is subunitary, i.e.,

 $\rho(M)<1.$

It follows that

 $\begin{bmatrix}\|S_{1}(\lambda)-S_{1}(\mu)\|_{\theta}\\[5.69054pt] \|S_{2}(\lambda)-S_{2}(\mu)\|_{\theta}\end{bmatrix}\leq(I_{2}-M)^{-1}\begin{% bmatrix}0\\[5.69054pt] \beta|\lambda-\mu|\end{bmatrix}$

or, equivalently,

 $\displaystyle\|S_{1}(\lambda)-S_{1}(\mu)\|_{\theta}\leq\beta\gamma_{12}|% \lambda-\mu|$ $\displaystyle\|S_{2}(\lambda)-S_{2}(\mu)\|_{\theta}\leq\beta\gamma_{22}|% \lambda-\mu|$

where $\gamma_{ij}$ are the entries of the matrix $(I_{2}-M)^{-1}$. Clearly, from these inequalities, we may conclude the Lipschitz continuity on $\lambda$ of the solution.

Next, we assume that the following conditions hold

 $v(T)

From a biological point of view, this is a rational assumption meaning that the target value $v_{T}$ for the leukemic level is chosen below the level that would occur in the absence of any control and is above the level corresponding to a hypothetical zero growth rate of leukemia cells.

Taking into account the above arguments, we define the following iterative algorithm. It is aimed to bring us as close as possible to a value of $\lambda$ that corresponds to a solution of the control problem.

Step. 1.1. We initialize $\underline{\lambda}_{0}:=0$, $\overline{\lambda}_{0}:=1$

Step. 1.2. At any iteration $k\geq 1,$ we define $\lambda_{k}:=\frac{\underline{\lambda}_{k-1}+\overline{\lambda}_{k-1}}{2}$ and solve system (2.3) for $\lambda:=\lambda_{k}$. The numerical solution is obtained

 $(u_{k},v_{k})=(S_{1}(\lambda_{k}),S_{2}(\lambda_{k})),$

where $S_{1}$ and $S_{2}$ are given in the proof of Theorem 3.2.

If $v_{k}(T) then put $\underline{\lambda}_{k}=\lambda$, $\overline{\lambda}_{k}=\overline{\lambda}_{k-1}$,

otherwise, we take $\underline{\lambda}_{k}=\underline{\lambda}_{k-1}$, $\overline{\lambda}_{k}=\lambda$.

Step. 1.3. The algorithm stops if

 $|v_{k}(T)-v_{T}|<\delta,$

where $0<\delta\ll 1$. Otherwise, it continues with Step 1.2.

Concerning the above algorithm, we have the following convergence result.

###### Theorem 3.3.

The iterative algorithm Step 1.1–Step 1.3 is convergent to a solution of the control problem (2.3)-(2.4).

###### Proof.

For $k\geq 1$ we have the solution $(u_{k},v_{k})$ corresponding to $\lambda=\lambda_{k}$. Moreover, we have $v_{k}=S_{2}(\underline{\lambda}_{k})$ if $v_{k}(T) or $v_{k}=S_{2}(\overline{\lambda}_{k})$ if $v_{k}(T)\geq v_{T}$, where $S_{1}$ and $S_{2}$ are the operators defined in the proof of Theorem 3.2. Therefore, we obtain an increasing sequence ($\underline{\lambda}_{k}$) and a decreasing sequence ($\overline{\lambda}_{k}$) with the following properties

 (3.1) $S_{2}(\underline{\lambda}_{k})(T)

and

 (3.2) $\overline{\lambda}_{k}-\underline{\lambda}_{k}=\frac{1}{2^{k}}.$

The two sequences being monotone and bounded are convergent. Moreover, from (3.2) they have the same limit $\lambda^{*}$. Using the continuity of $S_{2}$ with respect to $\lambda$ and (3.1) we deduce that

 (3.3) $S_{2}(\lambda^{*})(T)=v_{T}.$

Finally, we use again the continuity of the operators $S_{1}$ and $S_{2}$ with respect to $\lambda$ to obtain $S_{1}(\lambda^{*})=u^{*}$ and $S_{2}(\lambda^{*})=v^{*}$. These together with (3.3) show that $(u^{*},v^{*},\lambda^{*})$ is a solution of the control problem (2.3)-(2.4).

Further, we provide some computational details concerning the implementation of this algorithm. More exactly, we consider the following $n$-point discretization of the interval $[0,T],$

 $\Delta_{n}=\{0=t_{1}

and we define the following subsets

 $\Delta_{i}=\{t_{1},t_{2},\ldots,t_{i}\}\subseteq\Delta_{n},\quad i=\overline{1% ,n}.$

We mention that at Step 1.2 we must solve, for each $k\geq 1$, a system of nonlinear Volterra integral equations, namely

 $\begin{cases}u_{k}(t)=u_{0}-ct+\displaystyle\int_{0}^{t}\frac{a}{1+b_{1}e^{u_{% k}(s)}+b_{2}e^{v_{k}(s)}}\,ds\\[5.69054pt] v_{k}(t)=v_{0}-Ct+\lambda_{k}\displaystyle\int_{0}^{t}\frac{A}{1+B\left(e^{u_{% k}(s)}+e^{v_{k}(s)}\right)}\,ds.\end{cases}$

To this end, we introduce the notation $Q(f,\Delta_{i})$ for a general quadrature formula that uses the nodes from $\Delta_{i}$, i.e,

 $\displaystyle\int^{t_{i}}_{0}f(s)\,ds\approx\sum^{i}_{j=1}\alpha_{j}f(t_{j}),$

where $f:[0,T]\to{\mathbb{R}_{+}}$ and $\alpha_{j}\in{\mathbb{R}}$ represent the coefficients. Moreover, we denote by $(u^{i}_{k},v^{i}_{k})$ the numerical approximations of $(u_{k},v_{k})$ at the point $t_{i}\in\Delta_{n},\ i=\overline{1,n}\$, where $u_{k}^{1}:=u_{0}$ and $v_{k}^{1}:=v_{0}$. They are the solutions of the nonlinear systems

 $\begin{cases}u^{i}_{k}=u^{1}_{k}-ct_{i}+Q(\frac{a}{1+b_{1}e^{u}+b_{2}e^{v}},% \Delta_{i})\\[5.69054pt] v^{i}_{k}=v^{1}_{k}-Ct_{i}+\lambda_{k}Q(\frac{A}{1+B(e^{u}+e^{v})},\Delta_{i})% \end{cases}\quad i=\overline{2,n}.$

In the case of the trapezoidal rule for equidistant nodes, the previous systems take the form

 (3.4) $\begin{cases}X=Const_{1}(u,v,i)+\frac{h}{2}\frac{a}{1+b_{1}e^{X}+b_{2}e^{Y}}\\% [5.69054pt] Y=Const_{2}(u,v,i)+\frac{h}{2}\frac{\lambda_{k}A}{1+B(e^{X}+e^{Y})}\end{cases}$

where

 $Const_{1}(u,v,i)=u^{1}_{k}-ct_{i}+\frac{h}{2}\frac{a}{1+b_{1}e^{u^{1}_{k}}+b_{% 2}e^{v^{1}_{k}}}+h\displaystyle\sum^{i-1}_{j=2}\frac{a}{1+b_{1}e^{u^{j}_{k}}+b% _{2}e^{v^{j}_{k}}}$
 $Const_{2}(u,v,i)=v^{1}_{k}-Ct_{i}+\frac{h}{2}\frac{\lambda_{k}A}{1+B(e^{u^{1}_% {k}}+e^{v^{1}_{k}})}+h\displaystyle\sum^{i-1}_{j=2}\frac{\lambda_{k}A}{1+B(e^{% u^{j}_{k}}+e^{v^{j}_{k}})}$

and $h$ denotes the constant length between discretization points, i.e., $h=t_{i}-t_{i-1},\,i=\overline{2,n}$. To solve the nonlinear systems (3.4) we can use a specific numerical method, such as the method of successive approximations.

### 3.2. The Second Algorithm

In this subsection, we propose a second algorithm which has the following structure.

Step. 2.1. We consider two start functions $\overline{u},\,\overline{v}:[0,T]\to{\mathbb{R}_{+}}$ satisfying

 $\overline{u}(0)=u_{0},\ ,\overline{v}(0)=v_{0},\ ,\overline{v}(T)=v_{T}.$

Step 2.2. Using (2.5) we compute

 $\overline{\lambda}:=\frac{v_{T}-v_{0}+CT}{\displaystyle\int_{0}^{T}\frac{A}{1+% B(e^{\overline{u}(s)}+e^{\overline{v}(s)})}\,ds}.$

Step 2.3. We solve Volterra system (2.3) for $\lambda=\overline{\lambda}$ and obtain its solution $(u,v).$

Step 2.4. The algorithm stops if the following conditions are satisfied

 $\|u-\overline{u}\|_{\infty}<\delta,\quad\|v-\overline{v}\|_{\infty}<\delta,$

where $0<\delta\ll 1$. Otherwise, there are performed Steps 2.2 and 2.3 for $\overline{u}:=u$,  $\overline{v}:=v$.

### 3.3. Numerical Results

We finish this section with some numerical results which are illustrated in the figures below and whose biological analysis is given in the next section. The numerical results are obtained using a grid $\Delta_{n}$ with equidistant points. We use the trapezoidal method and we take $\delta=10^{-12}$ in the stop condition. Concerning the parameters, the following values are taken (see I.Ş. Haplea et al. [22]):

 $\displaystyle a=5,\,A=7,\,b_{1}=$ $\displaystyle 0.75\times 10^{-5},\,b_{2}=0.38\times 10^{-5},\,B=0.19\times 10^% {-5},\,c=0.05,\,C=0.2$ $\displaystyle u_{0}=\ln(7\times 10^{6}),\quad v_{0}=\ln(10^{7}),\quad T=2000$

while for $v_{T}$ we consider values from the interval $[11.1181,15.6181]$.

Figures (a) and (b) give the graphs of the functions $u$ and $v$ for $v_{T}\in\{11.11,14.01,15.61\}$. Notice the asymptotic behavior of these functions and their inverse monotony with respect to $v_{T}.$ Indeed, at any time $t$, the value $v(t)$ decreases and the value $u(t)$ increases as the target value $v_{T}$ decreases.

Figure (c) shows the variation of the control parameter $\lambda$ with respect to $v_{T}$.

Figures (d)-(f) have the role to give some details about the efficiency of the algorithm Step 1.1-Step 1.3. Thus, the histogram in Figure (d) shows the frequency of the number of iterations in Step 1.2 necessary to solve the nonlinear integral system, in the case when $v_{T}$ takes $46$ equidistant values in the interval $[11.1181,15.6181]$, that is when the difference between two consecutive values of $v_{T}$ is constant and equal with $0.1$. It is observed that around $45-49$ iterations are necessary to reach the desired error, i.e., $|v(T)-v_{T}|<10^{-12}$.

Figures (e) and (f) show the behavior of the error $\ lg|v(T)-v_{T}|$ with respect to the number of iterations, for $v_{T}=11.1181$ and $v_{T}=15.6181$.

In Figures (g) and (h), there are represented the graphs of functions $x=e^{u}$ and $y=e^{v}$. Some biological considerations about them are given in the next section.

Finally, we mention that similar results are obtained using the second algorithm Step 2.1 -Step 2.4. In this case, the number of iterations strongly depends on the start functions taken in Step 2.1. Also, the approximation of $\overline{\lambda}$ in Step 2.2 increases the complexity of the algorithm. However, in contrast to the first algorithm, the error decreases in a monotone manner.

## 4. Discussion

The results of the numerical simulations (Figures (b) and (h)) agree with the data of CML cases treated with TKIs (see M.C. Müller et al. [29] and F. Michor et al. [24]); malignant cell counts decrease monotonically in time, with a sharp decline for the first $150$ to $200$ time units (days of treatment), followed by a slower phase that eventually plateaus out. Normal cells expand in compensation (Figures (a) and (g)). As expected, the more stringent the target (a lower $v_{T}$), the longer it takes to approach and converge to it (Figure (h)). To answer the question: what drug effect is required to suppress the malignant population down to a given level, we have explored numerically the relationship between $\lambda$ and $v_{T}$ (Figure (c)). We have found a monotonic relationship, where $\lambda$ increases faster than linearly with $v_{T}$, and in consequence incrementally larger drug loads are required for every additional $log$ unit of reduction in the malignant cell population. This confirms that hematological cure by dose escalations alone can be prohibitively costly in terms of drug toxicity for unresponsive patients, and alternative management (like switching to an alternative TKI) is best considered in such cases.

Normal and leukemic cell populations are internally structured, with subpopulations of stem, progenitor and variably differentiated cells. We have collated all the cell subtypes within a single normal cell phenotype and another single malignant cell phenotype. Also, myeloid leukemias develop functionally in two compartments, the hematopoietic bone marrow and the peripheral blood. While the bone marrow (the production compartment) is severely space constrained, the more compliant circulatory blood can accommodate massive increases in myeloid cell counts. These increases range from severalfold (physiologically) to almost two orders of magnitude (pathologically) (see J. Gong et al. [30]). In our model, we have compacted the generative central (medullary) compartment and the transit peripheral (blood) compartment in a single compartment. This amounts to assuming constant rates of cell transport between compartments. While the numbers of modeled cells, normal and malignant, need not add up to a constant ($x+y$ can vary in time), there is an implicit pressure in the model toward lower counts through allo- and hetero feedback. This is in qualitative accord with actual biology, as even the bone marrow has some margin for variable cell density (histologically, it can be of ”higher” or ”lower cellularity”).

In current hematology practice, the CML state is assessed by the ratio of BCR-ABL1 transcripts to ABL1 transcripts, rather than by the absolute levels of BCR-ABL1. We have departed from this standard, as we have chosen for the goal of the optimization, the absolute number of leukemic cells $exp(v_{T})$. This makes our method more general, as it is readily extendable to solid malignancies where absolute cell counts can be derived from imaging studies.

The goal of CML treatment with TKIs is to achieve long-term molecular remission. Earlier evaluations at shorter terms ($3$ and $6$ months of treatment) serve only to predict the long-term response (see B. Hanfstein et al. [31]), but otherwise have no clinical utility on their own. If the treatment fails to lower BCR-ABL1$\%$ enough at three months, the patient is classified as a nonresponder and switched to an alternative, likely more efficient TKI. Imatinib (a first-generation molecule) is the agent most often used as a first-line for CML, while the more potent second and third-generation TKIs may be withheld due to the cost and toxicity profile. A direction for future research is to model the treatment with different TKIs in succession, each with its specific $\lambda$ value, and to identify the optimal time points of the switch from one drug to the other, in view of some goal of persistent remission.

In the work I.Ş. Haplea et al. [22] we have also optimized the (modeled) treatment of acute leukemia, using as a control objective the whole timecourse of the disease under treatment (i.e., the leukemic cells should decrease along a predefined curve), which reflects the need to monitor tightly an acute disease. For chronic leukemia, it appears reasonable to use solely the long-term endpoint as a goal. Nonetheless, as the existence of the solution $(u,v,\lambda)$ has been proved for the chronic form, the endpoint will be reached, in every case, along a single uniquely determined curve. It would be of interest to check the form of the so-established curve against the intermediate targets from clinical guidelines. Our model may predict that the final target can be attained even if early targets are not met. Both coincidence or disparity between the model and reality would lead to some insight.

## Acknowledgment:

The authors would like to thank the reviewers for their careful reading of the manuscript and for their valuable comments and suggestions.

L.G. Parajdi benefited from financial support through the project ”The Development of Advanced and Applicative Research Competencies in the Logic of STEAM+ Health”, POCU/993/6 /13/153310, a project co-financed by the European Social Fund through The Romanian Operational Programme Human Capital 2014-2020.

## References

• [1] E. Jabbour and H. Kantarjian, Chronic Myeloid Leukemia: 2020 Update on Diagnosis, Therapy and Monitoring. Am. J. Hematol. 95 (6) (2020), 691–709.
• [2] S. Faderl, M. Talpaz, Z. Estrov, S. O’Brien, R. Kurzrock and H.M. Kantarjian, The Biology of Chronic Myeloid Leukemia. N. Engl. J. Med. 341 (3) (1999), 164–172.
• [3] S.G. O’Brien, F. Guilhot, R.A. Larson, I. Gathmann, M. Baccarani, F. Cervantes, J.J. Cornelissen, et al., Imatinib Compared with Interferon and Low-Dose Cytarabine for Newly Diagnosed Chronic-Phase Chronic Myeloid Leukemia. N. Engl. J. Med. 348 (11) (2003), 994–1004.
• [4] A. Hochhaus, R.A. Larson, F. Guilhot, J.P. Radich, S. Branford, T.P. Hughes, M. Baccarani, et al. Long-Term Outcomes of Imatinib Treatment for Chronic Myeloid Leukemia. N. Engl. J. Med. 376 (10) (2017), 917–927.
• [5] T. Hughes, M. Deininger, A. Hochhaus, S. Branford, J. Radich, J. Kaeda, M. Baccarani, et al., Monitoring CML Patients Responding to Treatment with Tyrosine Kinase Inhibitors: Review and Recommendations for Harmonizing Current Methodology for Detecting BCR-ABL Transcripts and Kinase Domain Mutations and for Expressing Results. Blood 108 (1) (2006), 28–37.
• [6] M.W. Deininger, N.P. Shah, J.K. Altman, E. Berman, R. Bhatia, B. Bhatnagar, D.J. DeAngelo, et al., Chronic Myeloid Leukemia, Version 2.2021, NCCN Clinical Practice Guidelines in Oncology. J. Natl. Compr. Canc. Netw. 18 (10) (2020), 1385–1415.
• [7] A. Hochhaus, M. Baccarani, R.T. Silver, C. Schiffer, J.F. Apperley, F. Cervantes, R.E. Clark, et al., European LeukemiaNet 2020 Recommendations for Treating Chronic Myeloid Leukemia. Leukemia 34 (4) (2020), 966–984.
• [8] Y. Cherruault, Global optimization in biology and medicine. Mathl. Comput. Modelling 20 (6) (1994), 119–132.
• [9] C.J.S. Alves, P.M. Pardalos and L.N. Vicente, Optimization in Medicine. Springer, New York, 2008.
• [10] Y. Wang, X-S. Zhang and L. Chen, Optimization meets systems biology. BMC Syst. Biol. 4 (suppl 2) (2010), 1–4.
• [11] G. Cedersund, O. Samuelsson, G. Ball, J. Tegnér and D. Gomez-Cabrero, Optimization in biology parameter estimation and the associated optimization problem. In: Uncertainty in Biology, A Computational Modeling Approach. Springer, Chem, (2016), 177–197.
• [12] S. Nanda, H. Moore and S. Lenhart, Optimal control of treatment in a mathematical model of chronic myelogenous leukemia. Math. Biosci. 210 (1) (2007), 143–156.
• [13] B. Ainseba and C. Benosman, Optimal control for resistance and suboptimal response in CML. Math. Biosci. 227 (2) (2010), 81–93.
• [14] Q. He, J. Zhu, D. Dingli, J. Foo and K.Z. Leder, Optimized treatment schedules for chronic myeloid leukemia. PLoS Comput. Biol. 12 (10) (2016), e1005129.
• [15] S.B. Mendrazitsky and B. Shklyar, Optimization of combined leukemia therapy by finite-dimensional optimal control modeling. Springer, J. Optim. Theory. Appl. 175 (1) (2017), 218–235.
• [16] S.B. Mendrazitsky, N. Kronik and V. Vainstein, Optimization of interferon-alpha and imatinib combination therapy for chronic myeloid leukemia: a modeling approach. Adv. Theory. Simul. 2 (1) (2019), 1800081–8.
• [17] M.C. Mackey and L. Glass, Oscillation and chaos in physiological control systems. Science 197 (4300) (1977), 287–289.
• [18] L.G. Parajdi, R. Precup, E.A. Bonci and C. Tomuleasa, A mathematical model of the transition from normal hematopoiesis to the chronic and accelerated-acute stages in myeloid leukemia. Mathematics 8 (3) (2020), 376.
• [19] D. Dingli and F. Michor, Successful therapy must eradicate cancer stem cells. Stem Cells 24 (12) (2006), 2603–2610.
• [20] A. Cucuianu and R. Precup, A hypothetical-mathematical model of acute myeloid leukemia pathogenesis. Comput. Math. Methods Med. 11 (1) (2010), 49–65.
• [21] L.G. Parajdi, Analysis of Some Mathematical Models of Cell Dynamics in Hematology. Casa Cărţii de Ştiinţă, Cluj-Napoca, 2021.
• [22] I.Ş. Haplea, L.G. Parajdi and R. Precup, On the controllability of a system modeling cell dynamics related to leukemia. Symmetry, 13 (10) (2021), 1867.
• [23] V. Barbu, Mathematical Methods in Optimization of Differential Systems. Springer Science+Business Media: Dordrecht, 1994.
• [24] F. Michor, T.P. Hughes, Y. Iwasa, S. Branford, N.P. Shah, C.L. Sawyers and M.A. Nowak, Dynamics of Chronic Myeloid Leukaemia. Nature 435 (7046) (2005), 1267–1270.
• [25] A. Hochhaus, S. Saussele, G. Rosti, F.-X. Mahon, J.J.W.M. Janssen, H. Hjorth-Hansen, J. Richter, C. Buske, et al., Chronic myeloid leukaemia: ESMO Clinical Practice Guidelines for diagnosis, treatment and follow-up. Ann. Oncol. 28 (supp$\_$4) (2017), iv41-iv51.
• [26] R. Precup, Methods in Nonlinear Integral Equations. Kluwer Academic, Dordrecht, 2002.
• [27] E. Hairer, S.P. Nørsett and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems. 2nd revised edn., Springer Series in Computational Mathematics 8, Springer-Verlag, Berlin, Heidelberg, 1993.
• [28] D. Bufnea, V. Niculescu, G. Silaghi, A. Sterca, Babeş-Bolyai University’s High Performance Computing Center. Stud. Univ. Babeş-Bolyai Inf., 61 (2) (2016), 54–69.
• [29] M.C. Müller, N. Gattermann, T. Lahaye, M.W.N. Deininger, A. Berndt, S. Fruehauf, A. Neubauer, et al. Dynamics of BCR-ABL mRNA Expression in First-Line Therapy of Chronic Myelogenous Leukemia Patients with Imatinib or Interferon Alpha/Ara-C. Leukemia 17 (12) (2003), 2392–2400.
• [30] J. Gong, B. Wu, T. Guo, S. Zhou, B. He and X. Peng, Hyperleukocytosis: A Report of Five Cases and Review of the Literature. Oncol. Lett. 8 (4) (2014), 1825-1827.
• [31] B. Hanfstein, M.C. Müller, R. Hehlmann, P. Erben, M. Lauseker, A. Fabarius, S. Schnittger, et al. Early Molecular and Cytogenetic Response Is Predictive for Long-Term Progression-Free and Overall Survival in Chronic Myeloid Leukemia (CML). Leukemia 26 (9) (2012), 2096–2102.

[1]  B. Ainseba and C. Benosman, Optimal control for resistance and suboptimal response in CML, Math. Biosci., 2010, 227(2), 81–93. doi: 10.1016/j.mbs.2010.06.005
[2]  C. J. S. Alves, P. M. Pardalos and L. N. Vicente, Optimization in Medicine, Springer Science & Business Media, New York, 2008.
[3]  V. Barbu, Mathematical Methods in Optimization of Differential Systems, Springer Science & Business Media, Dordrecht, 1994. Google Scholar
[4]  D. Bufnea, V. Niculescu, G. Silaghi and A. Sterca, Babeş-Bolyai University’s high performance computing center, Stud. Univ. Babeş-Bolyai Inf., 2016, 61(2), 54–69.
[5]  G. Cedersund, O. Samuelsson, G. Ball, J. Tegnér and D. Gomez-Cabrero, Optimization in biology parameter estimation and the associated optimization problem. In: Uncertainty in Biology, A Computational Modeling Approach, Springer, Cham, 2016, 177–197.
[6]  Y. Cherruault, Global optimization in biology and medicine, Mathl. Comput. Modelling, 1994, 20(6), 119–132. doi: 10.1016/0895-7177(94)90027-2
[7]   A. Cucuianu and R. Precup, A hypothetical-mathematical model of acute myeloid leukaemia pathogenesis, Comput. Math. Methods Med., 2010, 11(1), 49–65. doi: 10.1080/17486700902973751
[8]  M. W. Deininger, N. P. Shah, et al., Chronic myeloid leukemia, Version 2.2021, NCCN Clinical practice guidelines in oncology, J. Natl. Compr. Canc. Netw., 2020, 18(10), 1385–1415. doi: 10.6004/jnccn.2020.0047
[9]  D. Dingli and F. Michor, Successful therapy must eradicate cancer stem cells, Stem Cells, 2006, 24(12), 2603–2610. doi: 10.1634/stemcells.2006-0136
[10]  S. Faderl, M. Talpaz, Z. Estrov, S. O’Brien, R. Kurzrock and H. M. Kantarjian, The biology of chronic myeloid leukemia, N. Engl. J. Med., 1999, 341(3), 164–172. doi: 10.1056/NEJM199907153410306
[11]  J. Gong, B. Wu, T. Guo, S. Zhou, B. He and X. Peng, Hyperleukocytosis: A report of five cases and review of the literature, Oncol. Lett., 2014, 8(4), 1825–1827. doi: 10.3892/ol.2014.2326
[12]  E. Hairer, S. P. Nørsett and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd edn., Springer-Verlag, Berlin, Heidelberg, 1993.
[13]  B. Hanfstein, M. C. Müller, et al., Early molecular and cytogenetic response is predictive for long-term progression-free and overall survival in chronic myeloid leukemia (CML), Leukemia, 2012, 26(9), 2096–2102. doi: 10.1038/leu.2012.85
[14] I. Ş. Haplea, L. G. Parajdi and R. Precup, On the controllability of a system modeling cell dynamics related to leukemia, Symmetry, 2021, 13(10), 1867. doi: 10.3390/sym13101867
[15] Q. He, J. Zhu, D. Dingli, J. Foo and K. Z. Leder, Optimized treatment schedules for chronic myeloid leukemia, PLoS Comput. Biol., 2016, 12(10), e1005129. doi: 10.1371/journal.pcbi.1005129
[16] A. Hochhaus, R. A. Larson, et al., Long-term outcomes of Imatinib treatment for chronic myeloid leukemia, N. Engl. J. Med., 2017, 376(10), 917–927. doi: 10.1056/NEJMoa16093
[17]  A. Hochhaus, S. Saussele, et al., Chronic myeloid leukaemia: ESMO Clinical Practice Guidelines for diagnosis, treatment and follow-up, Ann. Oncol., 2017, 28, iv41-iv51.
[18]  A. Hochhaus, M. Baccarani, et al., European LeukemiaNet 2020 recommendations for treating chronic myeloid leukemia, Leukemia, 2020, 34(4), 966–984. doi: 10.1038/s41375-020-0776-2
[19]  T. Hughes, M. Deininger, et al., Monitoring CML patients responding to treatment with Tyrosine Kinase Inhibitors: review and recommendations for harmonizing current methodology for detecting BCR-ABL transcripts and kinase domain mutations and for expressing results, Blood, 2006, 108(1), 28–37. doi: 10.1182/blood-2006-01-0092
[20]  E. Jabbour and H. Kantarjian, Chronic myeloid leukemia: 2020 update on diagnosis, therapy and monitoring, Am. J. Hematol., 2020, 95(6), 691–709. doi: 10.1002/ajh.25792
[21]  M. C. Mackey and L. Glass, Oscillation and chaos in physiological control systems, Science, 1977, 197(4300), 287–289. doi: 10.1126/science.267326
[22]  S. B. Mendrazitsky and B. Shklyar, Optimization of combined leukemia therapy by finite-dimensional optimal control modeling, J. Optim. Theory. Appl., 2017, 175(1), 218–235. doi: 10.1007/s10957-017-1161-9
[23]  S. B. Mendrazitsky, N. Kronik and V. Vainstein, Optimization of interferon-alpha and imatinib combination therapy for chronic myeloid leukemia: a modeling approach, Adv. Theory. Simul., 2019, 2(1), 1800081–8. doi: 10.1002/adts.201800081
[24]  F. Michor, T. P. Hughes, Y. Iwasa, S. Branford, N. P. Shah, C. L. Sawyers and M. A. Nowak, Dynamics of chronic myeloid leukaemia, Nature, 2005, 435(7046), 1267–1270. doi: 10.1038/nature03669
[25] M. C. Müller, N. Gattermann, et al., Dynamics of BCR-ABL mRNA expression in first-line therapy of chronic myelogenous leukemia patients with Imatinib or Interferon Alpha/Ara-C, Leukemia, 2003, 17(12), 2392–2400. doi: 10.1038/sj.leu.2403157
[26] S. Nanda, H. Moore and S. Lenhart, Optimal control of treatment in a mathematical model of chronic myelogenous leukemia, Math. Biosci., 2007, 210(1), 143–156. doi: 10.1016/j.mbs.2007.05.003
[27] S. G. O’Brien, F. Guilhot, et al., Imatinib compared with interferon and low-dose cytarabine for newly diagnosed chronic-phase chronic myeloid leukemia, N. Engl. J. Med., 2003, 348(11), 994–1004. doi: 10.1056/NEJMoa022457
[28] L. G. Parajdi, R. Precup, E. A. Bonci and C. Tomuleasa, A mathematical model of the transition from normal hematopoiesis to the chronic and accelerated-acute stages in myeloid leukemia, Mathematics, 2020, 8(3), 376. doi: 10.3390/math8030376
[29] L. G. Parajdi, Analysis of Some Mathematical Models of Cell Dynamics in Hematology, Casa Cărţii de Ştiinţă, Cluj-Napoca, 2021.
[30] R. Precup, Methods in Nonlinear Integral Equations, Kluwer Academic, Dordrecht, 2002.
[31] Y. Wang, X. Zhang and L. Chen, Optimization meets systems biology, BMC Syst. Biol., 2010, 4(2), 1–4.