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.

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

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

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

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

We can observe that under the assumption $b_1>b_2$, both $x^{\ast }$ and $y^{\ast }$ are positive if and only if The case when corresponds to the normal hematologic state, the case corresponds to the chronic phase of the disease, and the case 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.

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

where $w$ is the state variable, $\lambda$ denotes the control variable, $W$ is the domain of the states, $\mathrm{\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,\mathrm{\Lambda }$ and $D\subset W\times \mathrm{\Lambda }$ are not necessarily structured sets and $N$ is any mapping from $W\times \mathrm{\Lambda }$ to $W$.

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

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

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

Consider the set-valued map $F:{\mathrm{\Sigma }}_1\to \mathrm{\Lambda }$  defined as

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

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

The above general principle of controllability can be applied to control problems related to 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

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

Let us change the variables as follows:

and denote

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

This is our fixed point equation $w=N(w,\lambda ),$ where $w=(u,v),$ $N=(N_1,N_2)$ and

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

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., ) (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

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).

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

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

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

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.

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.

Next, we assume that the following conditions hold

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.

If then put ${\underset{¯}{\lambda }}_k=\lambda$, ${\overline{\lambda }}_k={\overline{\lambda }}_{k-1}$,

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

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

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],$

and we define the following subsets

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

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

where $f:[0,T]\to {ℝ}_{+}$ and ${\alpha }_j\in ℝ$ represent the coefficients. Moreover, we denote by $(u_k^i,v_k^i)$ the numerical approximations of $(u_k,v_k)$ at the point $t_i\in {\mathrm{\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

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

where

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.

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

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 ${\mathrm{\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]):

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., .

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.