On the Controllability of a System Modeling Cell Dynamics Related to Leukemia

Abstract

In this paper, two control problems for a symmetric model of cell dynamics related to leukemia are considered. The first one, in connection with classical chemotherapy, is that the evolution of the disease under treatment should follow a prescribed trajectory assuming that the drug works by increasing the cell death rates of both malignant and normal cells. In the case of the second control problem, as for targeted therapies, the drug is assumed to work by decreasing the multiplication rate of leukemic cells only, and the control objective is that the disease state reaches a desired endpoint. The solvability of the two problems as well as their stability are proved by using a general method of analysis. Some numerical simulations are included to illustrate the theoretical results and prove their applicability. The results can possibly be used to design therapeutic scenarios such that an expected clinical evolution can be achieved.

1. Introduction

1.1. Biological Background

Hematological malignancies were the first types of cancer for which chemotherapy was attempted and then refined [1]. The blood stem cell was the first stem cell ever discovered (reviewed in [2]), and the well-established developmental hierarchy of hematopoiesis allowed for its rigorous functional characterization. Leukemias were the first malignancies in which cancer stem cells were documented (reviewed in [3]).

Classical cancer chemotherapeutics have strong cell cycle dependence [4] and, unlike targeted therapies, act preferentially on cancer cells rather than on normal cells, mainly because of the higher division rate of cancer cells. Any normal cell line with a high turnover rate will be also affected to a degree, as is the case with the highly active normal hematopoietic cells in the bone marrow. From the point of view of chemotherapy, the population of leukemia cells (to eliminate) and the population of normal cells (to be spared) are pharmacologically much closer and harder to discriminate between than, say, a colon cancer entoured by normal colonic mucosa or a hepatocarcinoma surrounded by normal hepatic parenchyma [5]. This narrower error margin makes leukemias an ideal proving ground for models of safety in chemotherapy and highlights the need for an accurate control of therapy. The diffuse and homogenous distribution of leukemic infiltrate within the marrow of axial bones [6] simplifies the assumptions of pharmacokinetic models, eliminates confusion factors, and allows focus on the study on the antineoplastic effect itself.

Unlike in solid tumors, leukemic cells do not possess the intracellular and cell surfaces necessary for tumor cohesion. The cellular mass of leukemias is fluidic, lacking, to a degree [7], a fixed and meaningful positional relationship between individual cells. As such, the spatial aspects of tumorigenesis can be reasonably abstracted away, which renders leukemias, as in the case of other types of population dynamics, especially suited for modeling with systems of ordinary differential equations.

In our model, the chemotherapeutics are not represented directly (as a system variable); they are instead modeled as an effect on the growth rates and removal rates of cells. It is unclear how to modify such parameters in response to a drug, as there is a conceptual gap between individual cell models and global tumor models. Taking the case of the common anti-tumor drug Paclitaxel [8], its prime effect on the cell is mitotic arrest: the affected cell will stop dividing and often (but not always) suffer apoptotic death. The way to proceed from this point is dependent on opinion more than fact. One can model its effect as an increase in the tumor cell death rate (the most common way), or as a combination of growth rate decrease and death rate increase. Our first problem, considered in Section 3, concerns a hypothetical drug that acts by increasing solely the death rates of both malignant and normal cells, while for the second problem (Section 4), the drug is assumed to act only on the growth rate of malignant cells.

One can expect that one of the directions of improvement in cancer therapeutics in the near future will consist of treatment algorithms. At the limit, adjusting the drug delivery to a very fine-grained schedule (see Remark 1 (b)) means identifying a continuous change of one or more kinetic parameters so that a clinical evolution can be achieved. If this detailed level of control proves to be beneficial, chemotherapy could be delivered via solutions that are already in clinical practice, such as infusion pumps, either extracorporeal or implanted.

There are two metrics of control over the course of a disease, relating to the disease trajectory and to the disease endpoint, respectively. Traditionally, the endpoint refers to therapeutic efficiency (reduction of the leukemic burden below a set threshold, inducing remission, etc.) and the parcourse to therapeutic safety (such as the total amount of toxic chemotherapy to which the patient is exposed over the entire treatment time course). Still, the delegation of efficiency only to the endpoint criteria is limiting. The recognition of cancer as a dynamic disease makes it obvious that some critical aspects, such as the risk of relapse, can only be taken into account by an analysis of its full trajectory in time [9]. As an immediate example, if the total leukemic cell number is massive for the largest part of the treatment, then the probability of inducing drug resistance is greater, compared to the situation where the leukemic cell number is rapidly reduced from the onset, even if in both cases the treatment manages to bring the leukemic cell number down to the same final size.

Consequently, it may be desirable not only to reach a given endpoint but also to reach it along a specific route in the space of possible disease states, which is complex and only partially observable. We have formulated these two situations separately: one situation is where the dynamics matter as well as the endpoint, and the tumoral mass must decrease under treatment along a prespecified curve (the first control problem, Section 3), and the other situation is where only the end state is of interest (i.e., bringing the tumoral mass to a given low value; Section 4—the second control problem).

In the first case, the drug works by increasing the cell death rate, acts on both normal and leukemic cells, and can model any classical chemotherapy in exclusive use for any acute form of leukemia (myeloid or lymphoblastic). The case cannot apply, however, to the modern treatment of blast crises occurring in chronic myeloid leukemia (acute-like form), in which a targeted drug is used in addition to classical therapy.

In the second case, the drug works by decreasing the multiplication rate, acts on leukemic cells only, and models the effect of targeted therapies used in isolation. This case refers specifically to the chronic or accelerated phases of chronic myeloid leukemia, involving a targeted drug from the tyrosine kinase inhibitor family (TKI) (Imatinib, Dasatinib, Nilotinib, Bosutinib, Ponatinib). These drugs have only marginal toxicity on normal cells from the bone marrow, and a TKI is the single antileukemic treatment administered.

Note that an interesting application of the first control problem lies outside of the field of leukemia treatment. If we replace the “normal” and “malignant” cell populations with two—both normal—cell lineages, with different responses to a growth signal, the problem relates to bone marrow homeostasis. Some other applications of our control problems, or of some of their variants, to population dynamics in general could be also imagined.

1.2. Mathematical Model and Approach

Several mathematical models for cell dynamics in hematology have been proposed in the last decades [10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]. Most of these, as is the case with the pioneering model of blood production from Mackey–Glass [27], use differential equations involving sigmoid or Hill functions, since hematopoiesis is a self-limiting process. The kind of model expressed as

$$\begin{array}{ccc}\hfill x^{\prime }\left(t\right)& =& \frac{ax\left(t\right)}{1+b\left(x\left(t\right)+y\left(t\right)\right)}-cx\left(t\right)\hfill \\ \hfill y^{\prime }\left(t\right)& =& \frac{Ay\left(t\right)}{1+B\left(x\left(t\right)+y\left(t\right)\right)}-Cy\left(t\right)\hfill \end{array}$$

was introduced by Dingli and Michor [16] in order to describe the time competition between the populations $x\left(t\right)$ and $y\left(t\right)$ of normal and leukemic hematopoietic stem cells. This system has only two nonzero equilibria (steady states) that correspond to the normal hematological state and the acute leukemic state, respectively, and thus the model does not support a stable equilibrium in which both cell populations are non-zero, as in the case of the chronic phase of the disease. This can, however, be achieved, as recently shown in [22], if more realistically, instead of a unique sensitivity coefficient b for both normal and leukemic populations, different sensitivity parameters $b_1,$$b_2$ with $b_1>b_2$ are considered. Thus, in [22], the authors arrived at the following model:

$$\begin{array}{ccc}\hfill x^{\prime }\left(t\right)& =& \frac{ax\left(t\right)}{1+b_{1}x\left(t\right)+b_{2}y\left(t\right)}-cx\left(t\right)\hfill \\ \hfill y^{\prime }\left(t\right)& =& \frac{Ay\left(t\right)}{1+B(x\left(t\right)+y\left(t\right))}-Cy\left(t\right)\hfill \end{array}$$

(1)

where $x\left(t\right),y\left(t\right)$ are the normal and malignant cell populations at time $t,$ respectively; $a,A$ are the growth rates; $c,C$ are the cell death rates; and $b_1,b_2,B$ are sensitivity parameters that govern the self-limiting process.

As is the case for many systems in population dynamics, the structure of this model is symmetrical. The fact that the sensitivity coefficients $b_1,b_2$ in the first equation are different $(b_1>b_2)$, unlike those in the second equation, where both are equal to $B,$ is only a simplification of the work, not a structural one, as explained on page 5 of [22]. The structural symmetry of the model is further found in its qualitative analysis as performed in [22]. Indeed, letting

$$d=\frac{1}{b_1}\left(\frac{a}{c}-1\right),D=\frac{1}{B}\left(\frac{A}{C}-1\right),d_1=\frac{1}{b_2}\left(\frac{a}{c}-1\right),$$

where the case $D<d$ corresponds to the normal hematologic state; the case $d<D<d_1$ corresponds to the chronic phase of the disease; and the case $d_1<D$ stands for the acute phase. In the first case, the system has the unique stable equilibrium $\left(d,0\right)$; in the second case, the equilibrium is of the form $\left(x^{*},y^{*}\right)$, where both $x^{*},y^{*}$ are positive; and in the third case, the stable equilibrium is $\left(0,D\right)$. We may say that the chronic case is the expression of symmetry, while the normal and acute cases introduce asymmetry.

The present paper reports on the controllability of system (1). More exactly, two controllability problems are considered. For the first control problem, in Section 3, the control is given by a function $\lambda \left(t\right)$ (which quantifies the effect of the administered drug dose, which is adjustable over time), which is intended to change the death rates $c,C$ of both normal and malignant cells such that a desired clinical evolution is achieved over a given period of time. The model is related to classical chemotherapy. For the second control problem, in Section 4, the control is a constant $\lambda$ (connected to a constant drug dose) intended to decrease the growth rate A of malignant cells. This model could be related to targeted therapies. The stability of the two control problems is also established. Due to this stability property, 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 the parameters of the problem.

The analysis is based on a new method for the controllability of an abstract fixed point equation and the stability of the general control problem, which is first introduced in Section 2. Some numerical simulations accompany the theoretical results from Section 3 in order to demonstrate their practical applicability.

2. Controllability of a Fixed Point Equation

The two control problems inspired by hematology that are considered in this paper are directly related to the general problem of the control systems as defined, for example, by V. Barbu [28] (p. 34): “The general problem of control theory is that of reconstituting a differential system (as a matter of fact, some of its parameters viewed as control variables) from certain properties of the solution”. In this section, we first give a mathematical formulation of the general problem of control theory, a solving principle, and a notion of stability. All of these are stated in terms of nonlinear analysis, more exactly of the fixed point theory.

2.1. A General Controllability Principle

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

$$\left\{\begin{array}{c}w=N\left(w,\lambda \right)\\ w\in W,\lambda \in \mathrm{\Lambda },\left(w,\lambda \right)\in D.\end{array}\right.$$

(2)

where w is the state variable, $\lambda$ is the control variable, W is the domain of the states, $\mathrm{\Lambda }$ is the domain of controls, and D is the controllability domain, usually given by means of some condition imposed to $w,$ or to both w and $\lambda .$ Notice the very general formulation of the control problem, in terms of sets, where $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.$

In this context, we say that the equation $w=N\left(w,\lambda \right)$ is controllable in $W\times \mathrm{\Lambda }$ with respect to $D,$ providing that problem (2) 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,

$$\begin{array}{ccc}\hfill \Sigma & =& \left\{\left(w,\lambda \right)\in W\times \mathrm{\Lambda }:w=N\left(w,\lambda \right)\right\},\hfill \\ \hfill {\Sigma }_1& =& \left\{w\in W:\mathrm{there}\mathrm{is}\lambda \in \mathrm{\Lambda }\mathrm{with}w=N\left(w,\lambda \right)\right\}.\hfill \end{array}$$

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

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

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

Obviously, $F\left(w\right)\ne \varnothing$ for every $w\in {\Sigma }_1.$

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

Proposition 1.

If for some extension $\tilde{F}:W\to \mathrm{\Lambda }$ of F from ${\Sigma }_1$ to $W,$ there exists a fixed point $w\in W$ of the set-valued map

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

i.e.,

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

(3)

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

Proof. 

Since $w\in W$ and $\lambda \in \mathrm{\Lambda },$ one has $w\in {\Sigma }_1.$ Hence, $\tilde{F}\left(w\right)=F\left(w\right).$ Then, $\lambda \in F\left(w\right)$ and from the definition of $F,$ it follows that $\left(w,\lambda \right)\in D.$ Therefore, $\left(w,\lambda \right)$ solves (2). □

Note that F and $\tilde{F}$ can be single-valued maps in particular.

2.2. Stability of the General Control Problem

Consider the control problem (2). Assume that $\left(\mathrm{\Lambda },d_{\mathrm{\Lambda }}\right)$ is a metric space and consider a family $\mathcal{D}=\left\{D_{\epsilon }\right\}$ of subsets of $W\times \mathrm{\Lambda },$ indexed by $\epsilon ,$ where $0\le \epsilon <{\epsilon }_0,$ with the following properties:

$$D_{\epsilon }\subset D_{{\epsilon }^{\prime }}\mathrm{for}0\le \epsilon \le {\epsilon }^{\prime }<{\epsilon }_0,D_0=D,\bigcap _{0<\epsilon \le {\epsilon }_0}{D}_{\epsilon }=D.$$

We say that a solution $\left(w,\lambda \right)$ of the control problem (2) is $\mathcal{D}$-stable if for every $\epsilon \in (0,{\epsilon }_0)$ there exists $\delta =\delta \left(\epsilon ,w,\lambda \right)>0$ such that

$$\left(\overline{w},\overline{\lambda }\right)\in \Sigma ,d_{\mathrm{\Lambda }}\left(\overline{\lambda },\lambda \right)\le \delta ⟹\left(\overline{w},\overline{\lambda }\right)\in D_{\epsilon }.$$

We say that the control problem is $\mathcal{D}$-stable if all its solutions are $\mathcal{D}$-stable.

The following proposition contains a sufficient condition for the $\mathcal{D}$-stability of the control problem (2). It relies on a double continuous dependence of $\lambda$ and of $w.$

Proposition 2.

Assume that $\left(W,d_W\right)$ and $\left(\mathrm{\Lambda },d_{\mathrm{\Lambda }}\right)$ are metric spaces and that the following conditions are satisfied:

(a) 

For each solution $\left(w,\lambda \right)\in \Sigma$ and $\overline{\epsilon }>0,$ there exists $\delta =\delta \left(\overline{\epsilon },w,\lambda \right)>0$ such that

$$\left(\overline{w},\overline{\lambda }\right)\in \Sigma ,d_{\mathrm{\Lambda }}\left(\overline{\lambda },\lambda \right)\le \delta ⟹d_W\left(\overline{w},w\right)\le \overline{\epsilon };$$

(b) 

For each solution $\left(w,\lambda \right)\in \Sigma \cap D$ and $\epsilon \in (0,{\epsilon }_0),$ there exists $\overline{\epsilon }=\overline{\epsilon }\left(\epsilon ,w,\lambda \right)>0$ such that

$$\left(\overline{w},\overline{\lambda }\right)\in \Sigma ,d_W\left(\overline{w},w\right)\le \overline{\epsilon }⟹\left(\overline{w},\overline{\lambda }\right)\in D_{\epsilon }.$$

Then, the control problem (2) is $\mathcal{D}$-stable.

Proof. 

Let $\left(w,\lambda \right)\in \Sigma \cap D$ be any solution of the control problem and let $\epsilon \in (0,{\epsilon }_0).$ Let $\overline{\epsilon }$ and $\delta$ be the numbers guaranteed by conditions (b) and (a), respectively. If $\left(\overline{w},\overline{\lambda }\right)\in \Sigma$ and $d_{\mathrm{\Lambda }}\left(\overline{\lambda },\lambda \right)\le \delta ,$ then from (a), we have $d_W\left(\overline{w},w\right)\le \overline{\epsilon }$ and next from (b), $\left(\overline{w},\overline{\lambda }\right)\in D_{\epsilon }.$ Hence, the solution $\left(w,\lambda \right)$ is $\mathcal{D}$-stable. The arbitrariness of the solution makes the control problem to be $\mathcal{D}$-stable. □

Next, the above general principles of controllability and stability are applied to the normal–leukemic system.

3. First Control Problem for the Normal–Leukemic System

We now apply Proposition 1 to the two-dimensional system of ordinary differential equations modeling normal–leukemic cell dynamics in leukemia, namely

$$\left\{\begin{array}{c}{x}^{\prime }\left(t\right)=\frac{ax\left(t\right)}{1+b_{1}x\left(t\right)+b_{2}y\left(t\right)}-\lambda \left(t\right)cx\left(t\right)\hfill \\ y^{\prime }\left(t\right)=\frac{Ay\left(t\right)}{1+B\left(x\left(t\right)+y\left(t\right)\right)}-\lambda \left(t\right)Cy\left(t\right)\hfill \\ x\left(0\right)=x_0,y\left(0\right)=y_0,\hfill \end{array}\right.$$

(4)

where $a>c,A>C,b_1>b_2>B.$ Here, $\lambda$ is the control function, and it is assumed that the patient starts the treatment with the initial cell populations $x\left(0\right)=x_0$ and $y\left(0\right)=y_0.$

3.1. Solving of the Control Problem

Let us change the variables as follows:

$$u:=lnx,v:=lny$$

and denote $u_0=u\left(0\right)=lnx\left(0\right),$$v_0=v\left(0\right)=lny\left(0\right).$ First, dividing the two equations by $x$ and $y,$ respectively, and then integrating, we obtain the integral system equivalent to the initial value problem (4), namely

$$\left\{\begin{array}{c}u\left(t\right)=u_0+{\int }_0^t\left(\frac{a}{1+b_1{e}^u+b_2{e}^v}-\lambda \left(s\right)c\right)ds\\ v\left(t\right)=v_0+{\int }_0^t\left(\frac{A}{1+B\left(e^u+e^v\right)}-\lambda \left(s\right)C\right)ds.\end{array}\right.$$

(5)

This is our equation $w=N\left(w,\lambda \right),$ where $w=\left(u,v\right),$$N=\left(N_1,N_2\right),$

$$\begin{array}{ccc}\hfill N_1\left(u,v,\lambda \right)\left(t\right)& =& u_0+{\int }_0^t\left(\frac{a}{1+b_1{e}^u+b_2{e}^v}-\lambda \left(s\right)c\right)ds,\hfill \\ \hfill N_2\left(u,v,\lambda \right)\left(t\right)& =& v_0+{\int }_0^t\left(\frac{A}{1+B\left(e^u+e^v\right)}-\lambda \left(s\right)C\right)ds.\hfill \end{array}$$

(6)

Here, $W=C\left(\left[0,T\right],{\mathbb{R}}^2\right),$ $\mathrm{\Lambda }=L^{\infty }\left(0,T\right),$ and the controllability domain D is introduced by means of an objective condition of the type

$$\frac{v\left(t\right)}{u\left(t\right)}=k\left(t\right)\frac{v_0}{u_0}=:r\left(t\right),t\in \left[0,T\right],$$

(7)

for a given function

$$0\le k\left(t\right)\le 1(0\le r\left(t\right)\le \frac{v_0}{u_0})$$

describing the expected response to treatment with $r\left(0\right)=v_0/u_0.$ Assume that $r\in W^{1,\infty }\left(0,T\right).$ For example, we may take the linear function

$$r\left(t\right)=k_0\frac{v_0}{u_0}\frac{t}{T}+\frac{v_0}{u_0}\left(1-\frac{t}{T}\right),t\in \left[0,T\right],$$

for which $r\left(0\right)=v_0/u_0$ and $r\left(T\right)=k_0{v}_0/u_0$ for some number $k_0<1.$ This corresponds to a linear response to the treatment, from the initial patient state $v_0/u_0,$ to the final expected state $k_0{v}_0/u_0.$ More generally, we may take $r\left(t\right)$ of the form

$$r\left(t\right)=\left\{\begin{array}{cc}\frac{v_0}{u_0},\hfill & 0\le t\le T_0\hfill \\ \frac{v_0}{u_0}{\left[1-\left(1-k_0^m\right){\left(\frac{t-T_0}{T-T_0}\right)}^n\right]}^{\frac{1}{m}},\hfill & T_0\le t\le T,\hfill \end{array}\right.$$

(8)

where $0\le T_0<T,$$n\ge 1$ and $m>0$ are fixed numbers, and again $r\left(0\right)=v_0/u_0$ and $r\left(T\right)=k_0{v}_0/u_0.$ Furthermore, one can take a piecewise linear function that gives values at some intermediate times as required by the medical protocol and verified by laboratory tests. Obviously, other objective conditions can be considered instead (see Section 3.4).

The interpretation is as follows: the patient begins the treatment with the initial ratio $ln{y}_0/ln{x}_0$, and during the treatment, it is expected that the ratio $lny\left(t\right)/lnx\left(t\right)$ is $k\left(t\right)$-times smaller than the initial ratio, and at the end T of the treatment time interval, it reaches a prescribed safe level. Of course, the control variable $\lambda$ by which the cell death rates $c$ and $C$ are increased in order to guarantee the desired patient evolution depends on the drug dose.

The choice of $v\left(t\right)/u\left(t\right)$ and $y\left(t\right)/x\left(t\right)$ (Section 3.4), as opposed to other possible tumor measures, is based on reasons grounded in medical practice. In hematology, the treatment outcome is evaluated as a ratio between malignant cells and total (blood) cells. This applies to the older cytology (the percentage of malignant cells out of total cells, as seen under the microscope) and to more recent molecular methods (the number of copies of the malignancy gene BCR-ABL vs. number of copies of normal ABL gene). The use of these ratios originates in the technical impossibility of having an absolute count of leukemic cells (or a proxy thereof) in a living organism. This is unlike the situation in general oncology, where the volume of solid tumors can be assessed by MRI, CT, etc., and where the number of malignant cells alone would have made sense as an objective function. The obtained mathematical predictions can be easily mapped to biomedical data such as BCR-ABL ratios and the percentage of blasts in the bone marrow or peripheral blood.

Theorem 1.

If $r\in W^{1,\infty }\left(0,T\right),$ $r\left(0\right)=v_0/u_0$ and $r\left(t\right)\ne C/c$ for all $t\in \left[0,T\right],$ then the control problem (5)–(7) has a unique solution $\left(u,v,\lambda \right)\in C\left(\left[0,T\right],{\mathbb{R}}^2\right)\times L^{\infty }\left(0,T\right).$

Proof. 

Let $\left(u,v\right)\in {\Sigma }_1,$ that is $\left(u,v\right)\in C\left(\left[0,T\right],{\mathbb{R}}^2\right)$ solves (5) for some $\lambda \in L^{\infty }\left(0,T\right).$ Using (5), replacing in (7) and solving in ${\int }_0^t\lambda \left(s\right)ds,$ we obtain

$$\begin{array}{ccc}& & {\int }_0^t\lambda \left(s\right)ds\hfill \\ & =& \frac{r\left(t\right)u_0-v_0}{r\left(t\right)c-C}+\frac{1}{r\left(t\right)c-C}\left(r\left(t\right){\int }_0^t\frac{a}{1+b_1{e}^u+b_2{e}^v}ds-{\int }_0^t\frac{A}{1+B\left(e^u+e^v\right)}ds\right)\hfill \\ & =& :H\left(u,v\right).\hfill \end{array}$$

(9)

Since $r\in W^{1,\infty }\left(0,T\right),$ one has $H:{\Sigma }_1\to W^{1,\infty }\left(0,T\right).$ Then $F:{\Sigma }_1\to L^{\infty }\left(0,T\right)$ is a single-valued map given by

$$F\left(u,v\right)=\lambda \left(t\right)=\frac{d}{dt}H\left(u,v\right).$$

Clearly, we may extend $H$ from ${\Sigma }_1$ to the whole space $C\left(\left[0,T\right],{\mathbb{R}}^2\right)$ by using the above expression of $H.$ Let $\tilde{H}:$ $C\left(\left[0,T\right],{\mathbb{R}}^2\right)\to W^{1,\infty }\left(0,T\right)$ be this extension. Then we may consider $\tilde{F}:$$C\left(\left[0,T\right],{\mathbb{R}}^2\right)\to L^{\infty }\left(0,T\right),$

$$\tilde{F}\left(u,v\right)=\frac{d}{dt}\tilde{H}\left(u,v\right).$$

Now, our fixed point equation $w=\tilde{N}\left(w\right)$ becomes

$$\left\{\begin{array}{c}u=N_1\left(u,v,\tilde{F}\left(u,v\right)\right)\hfill \\ v=N_2\left(u,v,\tilde{F}\left(u,v\right)\right).\hfill \end{array}\right.$$

(10)

Since $r\left(0\right)=v_0/u_0,$ one has $\tilde{H}\left(u,v\right)\left(0\right)=0$ and so

$${\int }_0^t\tilde{F}\left(u,v\right)\left(s\right)ds=\tilde{H}\left(u,v\right)\left(t\right).$$

Then using (6) we obtain

$$\begin{array}{ccc}\hfill {\tilde{N}}_1\left(u,v\right)\left(t\right)& =& u_0-\tilde{H}\left(u,v\right)c+{\int }_0^t\frac{a}{1+b_1{e}^u+b_2{e}^v}ds,\hfill \\ \hfill {\tilde{N}}_2\left(u,v\right)\left(t\right)& =& v_0-\tilde{H}\left(u,v\right)C+{\int }_0^t\frac{A}{1+B\left(e^u+e^v\right)}ds.\hfill \end{array}$$

Replacing $\tilde{H}\left(u,v\right)$ by its expression (9) gives

$$\begin{array}{ccc}\hfill {\tilde{N}}_1\left(u,v\right)& =& u_0-c\frac{r\left(t\right)u_0-v_0}{r\left(t\right)c-C}+\left(1-\frac{r\left(t\right)c}{r\left(t\right)c-C}\right){\int }_0^t\frac{a}{1+b_1{e}^u+b_2{e}^v}ds\hfill \\ & & +\frac{c}{r\left(t\right)c-C}{\int }_0^t\frac{A}{1+B\left(e^u+e^v\right)}ds,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill {\tilde{N}}_2\left(u,v\right)& =& v_0-C\frac{r\left(t\right)u_0-v_0}{r\left(t\right)c-C}-\frac{r\left(t\right)C}{r\left(t\right)c-C}{\int }_0^t\frac{a}{1+b_1{e}^u+b_2{e}^v}ds\hfill \\ & & +\left(1+\frac{C}{r\left(t\right)c-C}\right){\int }_0^t\frac{A}{1+B\left(e^u+e^v\right)}ds.\hfill \end{array}$$

Since the functions $\frac{a}{1+b_1{e}^u+b_2{e}^v}$ and $\frac{A}{1+B\left(e^u+e^v\right)}$ are Lipschitz continuous in both $u$ and $v$ in ${\mathbb{R}}^2,$ and the operator $N$ is of Volterra type, using Bielecki’s technique of equivalent norms and Banach’s contraction principle, one arrives to the conclusion that system (10) has a unique solution $\left(u,v\right).$ In addition, the solution is the limit of the successive approximation sequence starting with any initial point. Thus, system (4) is controllable. Finally, using the solution $\left(u,v\right),$ we first compute $H\left(u,v\right)$ according to (9) and then the control function $\lambda \left(t\right)$ by

$$\begin{array}{ccc}\hfill \lambda \left(t\right)& =& \frac{d}{dt}H\left(u,v\right)\hfill \\ & =& \frac{r^{\prime }\left(t\right)u_0}{r\left(t\right)c-C}-\frac{\left(r\left(t\right)u_0-v_0\right)r^{\prime }\left(t\right)c}{{(r\left(t\right)c-C)}^2}+\frac{r^{\prime }\left(t\right)}{r\left(t\right)c-C}{\int }_0^t\frac{a}{1+b_1{e}^u+b_2{e}^v}ds\hfill \\ & +& \frac{r\left(t\right)a}{(r\left(t\right)c-C)(1+b_1{e}^u+b_2{e}^v)}-\frac{r\left(t\right)r^{\prime }\left(t\right)c}{{\left(r\left(t\right)c-C\right)}^2}{\int }_0^t\frac{a}{1+b_1{e}^u+b_2{e}^v}ds\hfill \\ & -& \frac{A}{(r\left(t\right)c-C)(1+B\left(e^u+e^v\right))}+\frac{r^{\prime }\left(t\right)c}{{(r\left(t\right)c-C)}^2}{\int }_0^t\frac{A}{1+B\left(e^u+e^v\right)}ds.\hfill \end{array}$$

□

Remark 1.

Let $\left(u,v,\lambda \right)$ be the solution of the control problem.

(a) 

(control regularity) In case $r\in C^1\left[0,T\right],$ one has $\lambda \in C\left[0,T\right].$

(b) 

(continuous dependence) If $u_k\to u$ and $v_k\to v$ as $k\to \infty ,$ uniformly on $\left[0,T\right]$ and ${\lambda }_k\left(t\right):=dH\left(u_k,v_k\right)/dt,$ then ${\lambda }_k\to \lambda$ as $k\to \infty$ in $L^{\infty }\left(0,T\right).$

3.2. Stability of the Control Problem

We now discuss the stability of the control problem (5)–(7) by using Proposition 2.

Let $\left(u,v,\lambda \right)$ be the solution of the control problem with the objective condition $v\left(t\right)/u\left(t\right)=r\left(t\right).$ The stability issue that we consider is the following: how close to the exact drug prescription given by function $\lambda \left(t\right)$ should be the administrated treatment, let it be $\overline{\lambda }\left(t\right),$ in order that the final result $\overline{v}\left(T\right)/\overline{u}\left(T\right)=:\overline{r}\left(T\right)$ differs from the desired one $r\left(T\right)$ by at most $\epsilon ?$ More exactly, for a given $\epsilon >0,$ find $\delta >0$ such that

$${∥\overline{\lambda }-\lambda ∥}_{\infty }\le \delta \mathrm{implies}\left|\overline{r}\left(T\right)-r\left(T\right)\right|\le \epsilon .$$

The even more complete response is given by the following theorem.

Theorem 2.

Let $\left(u,v,\lambda \right)$ be the solution of the control problem (5)–(7). For a given $\overline{\lambda }\in L^{\infty }\left(0,T\right),$ we denote by $\left(\overline{u},\overline{v}\right)$ the corresponding solution of system (5) and by $\overline{r}\left(t\right)$ the ratio $\overline{v}\left(t\right)/\overline{u}\left(t\right).$ Then, for any small $\epsilon >0,$ denoting $\gamma =\left(c+C\right)T{e}^{\left(a+A\right)T}$, one has

$${∥\overline{\lambda }-\lambda ∥}_{\infty }\le \frac{\epsilon u\left(T\right)}{\left(1+r\left(T\right)-\epsilon \right)\gamma }=:\delta \mathit{implies}{∥\overline{r}-r∥}_{\infty }\le \epsilon .$$

(11)

Proof. 

Using the Lipschitz-type estimate

$$\left|\frac{a}{1+b_1{e}^u+b_2{e}^v}-\frac{a}{1+b_1{e}^{\overline{u}}+b_2{e}^{\overline{v}}}\right|\le a\left(\left|u-\overline{u}\right|+\left|v-\overline{v}\right|\right)$$

and the analogue inequality for $A/\left(1+B\left(e^u+e^v\right)\right)$ gives

$$\begin{array}{ccc}\hfill \left|u\left(t\right)-\overline{u}\left(t\right)\right|& \le & a{\int }_0^t\left(\left|u\left(s\right)-\overline{u}\left(s\right)\right|+\left|v\left(t\right)-\overline{v}\left(t\right)\right|\right)ds+\delta cT,\hfill \\ \hfill \left|v\left(t\right)-\overline{v}\left(t\right)\right|& \le & A{\int }_0^t\left(\left|u\left(s\right)-\overline{u}\left(s\right)\right|+\left|v\left(t\right)-\overline{v}\left(t\right)\right|\right)ds+\delta CT.\hfill \end{array}$$

Then, taking the sum

$$\begin{array}{ccc}\hfill \left|u\left(t\right)-\overline{u}\left(t\right)\right|+\left|v\left(t\right)-\overline{v}\left(t\right)\right|& \le & \left(a+A\right){\int }_0^t\left(\left|u\left(s\right)-\overline{u}\left(s\right)\right|+\left|v\left(t\right)-\overline{v}\left(t\right)\right|\right)ds\hfill \\ & & +\delta \left(c+C\right)T\hfill \end{array}$$

and Gronwall’s inequality yields

$$\left|u\left(t\right)-\overline{u}\left(t\right)\right|+\left|v\left(t\right)-\overline{v}\left(t\right)\right|\le \delta \left(c+C\right)T{e}^{\left(a+A\right)T}=\delta \gamma =:\overline{\epsilon }.$$

This implies

$$\begin{array}{ccc}\hfill u\left(t\right)-\overline{\epsilon }& \le & \overline{u}\left(t\right)\le u\left(t\right)+\overline{\epsilon },\hfill \\ \hfill v\left(t\right)-\overline{\epsilon }& \le & \overline{v}\left(t\right)\le v\left(t\right)+\overline{\epsilon },\hfill \end{array}$$

where

$$\frac{v\left(t\right)-\overline{\epsilon }}{u\left(t\right)+\overline{\epsilon }}-r\left(t\right)\le \overline{r}\left(t\right)-r\left(t\right)\le \frac{v\left(t\right)+\overline{\epsilon }}{u\left(t\right)-\overline{\epsilon }}-r\left(t\right),$$

or equivalently

$$-\frac{\overline{\epsilon }\left(1+r\left(t\right)\right)}{u\left(t\right)+\overline{\epsilon }}\le \overline{r}\left(t\right)-r\left(t\right)\le \frac{\overline{\epsilon }\left(1+r\left(t\right)\right)}{u\left(t\right)-\overline{\epsilon }}.$$

Then,

$$\left|\overline{r}\left(t\right)-r\left(t\right)\right|\le \frac{\overline{\epsilon }\left(1+r\left(t\right)\right)}{u\left(t\right)+\overline{\epsilon }}\le \epsilon \left(t\in \left[0,T\right]\right),$$

(12)

where $\epsilon \ge \overline{\epsilon }{max}_{t\in \left[0,T\right]}\frac{1+r\left(t\right)}{u\left(t\right)+\overline{\epsilon }}.$

Thus, the conditions of Proposition 2 are fulfilled for

$$\begin{array}{ccc}\hfill D_{\epsilon }& :& =\{\left(w,\lambda \right)\in C\left(\left[0,T\right],{\mathbb{R}}^2\right)\times W^{1,\infty }\left(0,T\right):\hfill \\ \hfill w& =& \left(u,v\right),\left|v\left(t\right)-r\left(t\right)u\left(t\right)\right|\le \epsilon \left|u\left(t\right)\right|\mathrm{for}\mathrm{all}t\in \left[0.T\right]\},\hfill \end{array}$$

$$\overline{\epsilon }=\epsilon \underset{t\in \left[0,T\right]}{min}\frac{u\left(t\right)}{1-\epsilon +r\left(t\right)},\delta =\frac{\overline{\epsilon }}{\gamma },$$

which shows that the control problem (2) is $\mathcal{D}$-stable.

Finally, conclusion (11) is immediate from (12). □

Remark 2.

The medical interpretation of the result in Theorem 2 is as follows: if the administrated drug dose is close enough to the prescribed dose given by the solution of the control problem, then the result of the treatment applied to the patient does not differ much from the intended result.

3.3. Numerical Simulations

Below, we illustrate the theoretical results obtained for the first controllability problem by numerical simulations using the Maple software and an modified package for numerical and graphical solutions of Volterra integral equations. The initial Maple package for numerical solutions of Volterra integral equations was published by L.C. Backer and M. Wheeler [29]. For more information about numerical methods for Volterra integral equations, see the books by T.A. Burton [30] and P. Linz [31].

In this section, we compute some numerical simulations considering the system of Volterra integral equations (10); more exactly,

$$\left\{\begin{array}{c}u\left(t\right)=u_0-c\frac{r\left(t\right)u_0-v_0}{r\left(t\right)c-C}+{\int }_0^t\left\{\left(1-\frac{r\left(t\right)c}{r\left(t\right)c-C}\right)\frac{a}{1+b_1{e}^{u\left(s\right)}+b_2{e}^{v\left(s\right)}}+\frac{c}{r\left(t\right)c-C}\frac{A}{1+B\left(e^{u\left(s\right)}+e^{v\left(s\right)}\right)}\right\}ds\hfill \\ v\left(t\right)=v_0-C\frac{r\left(t\right)u_0-v_0}{r\left(t\right)c-C}+{\int }_0^t\left\{-\frac{r\left(t\right)C}{r\left(t\right)c-C}\frac{a}{1+b_1{e}^{u\left(s\right)}+b_2{e}^{v\left(s\right)}}+\left(1+\frac{C}{r\left(t\right)c-C}\right)\frac{A}{1+B\left(e^{u\left(s\right)}+e^{v\left(s\right)}\right)}\right\}ds.\hfill \end{array}\right.$$

For these numerical simulations, we may take the expression of the function $r\left(t\right)$ given by (8).

We restrict our simulations to two cases: Case 1 (the chronic phase of the disease), where the parameters of the nonlinear Volterra integral system have values such that the condition of chronicity $(d<D<d_1)$ holds $a=5,$ $b_1=0.75\times {10}^{-5},$ $b_2=0.38\times {10}^{-5},$ $c=0.05,$ $A=7,$ $B=0.19\times {10}^{-5},$ $C=0.2$; and Case 2 (the accelerated–acute phase of the disease), where the parameters of the nonlinear Volterra integral system have values such that the condition of the accelerated–acute phase $(d_1<D)$ holds $a=5,$ $b_1=0.75\times {10}^{-5},$ $b_2=0.38\times {10}^{-5},$ $c=0.05,$ $A=7,$ $B=0.19\times {10}^{-5},$ $C=0.1$. The parameters of the function $r\left(t\right)$ are $k_0=0.1,$ $T=20,$ $n=2,$ $m=1/2,$ $T_0=0$ and the initial conditions for this two simulations are $u_0=15.76142071,$ $v_0=16.11809565$. The above parameter values were first estimated in the paper of Dingli and Michor [16], and we have analyzed this set of parameters in our previous work Parajdi et al. [22]. The step size for the numerical integration is $h=0.1$.

Note that, although no longer relevant to clinical practice, in order to illustrate how the overall appearance of $\lambda \left(t\right)$ changes with different model parameters, we have made computations for a parameter set of chronic myeloid leukemia (CML) as well. This describes a purely hypothetical case in which a chronic disease was treated with conventional chemotherapy (for the chronic phase of CML see Figure 1 and Figure 2, and

Table 1. Numerical simulations of the system of two nonlinear Volterra integral equations, given by (10), obtained in Case 1 (the chronic phase of the disease).

t	u(t)	v(t)	v(t)/u(t)	$\mathit{\lambda }\left(\mathit{t}\right)$
0.0	15.76142	16.11810	1.02263	1.03667861
1.0	15.74559	16.04691	1.01914	1.82681929
2.0	15.69218	15.82858	1.00869	2.72096808
3.0	15.60090	15.46682	0.99141	3.73208846
4.0	15.47141	14.96789	0.96745	4.88700566
5.0	15.30324	14.34054	0.93709	6.21220861
6.0	15.09671	13.59667	0.90064	7.72168000
7.0	14.85444	12.75238	0.85849	9.42156274
8.0	14.58279	11.82826	0.81111	11.31917081
9.0	14.29297	10.84895	0.75904	13.40181913
10.0	14.00196	9.84183	0.70289	15.58459339
11.0	13.73138	8.83391	0.64334	17.68482475
12.0	13.50188	7.84645	0.58114	19.48688614
13.0	13.32522	6.89067	0.51712	20.85563950
14.0	13.20103	5.96906	0.45217	21.77786674
15.0	13.12077	5.08115	0.38726	22.31403995
16.0	13.07384	4.22854	0.32343	22.53852452
17.0	13.05129	3.41688	0.26180	22.51022240
18.0	13.04694	2.65569	0.20355	22.26679888
19.0	13.05703	1.95760	0.14993	21.82786481
20.0	13.07965	1.33756	0.10226	21.19948845

, and in addition for the accelerated-acute phase of CML see Figure 3 and Figure 4, and

Table 2. Numerical simulations of the system of two nonlinear Volterra integral equations, given by (10), obtained in Case 2 (the accelerated–acute phase of the disease).

t	u(t)	v(t)	v(t)/u(t)	$\mathit{\lambda }\left(\mathit{t}\right)$
0.0	15.76142	16.11810	1.02263	3.15804
1.0	15.59321	15.89161	1.01914	6.05741
2.0	15.27793	15.41073	1.00869	9.93469
3.0	14.76795	14.64103	0.99141	16.19792
4.0	13.93492	13.48141	0.96745	29.76677
5.0	12.57768	11.78644	0.93709	57.07073
6.0	11.49280	10.35085	0.90064	66.27579
7.0	11.14817	9.57058	0.85849	66.82106
8.0	11.05090	8.96350	0.81111	67.63740
9.0	11.01412	8.36017	0.75904	68.47741
10.0	10.99174	7.72598	0.70289	69.17105
11.0	10.97355	7.05969	0.64334	69.69562
12.0	10.95824	6.36825	0.58114	70.06188
13.0	10.94651	5.66061	0.51712	70.28064
14.0	10.93937	4.94642	0.45217	70.35631
15.0	10.93770	4.23573	0.38726	70.28740
16.0	10.94224	3.53910	0.32343	70.06809
17.0	10.95364	2.86770	0.26180	69.68924
18.0	10.97248	2.23344	0.20355	69.13897
19.0	10.99933	1.64909	0.14993	68.40266
20.0	11.03477	1.12845	0.10226	67.46265

).

3.4. A Different Objective Condition

As already noted, some other objective conditions can be considered for the first control problem. One of these, with immediate practical importance, consists of following a predefined path for the ratio $y\left(t\right)/x\left(t\right)$ instead of $lny\left(t\right)/lnx\left(t\right)$ from (7). Hence, given a function $r\left(t\right)$ with $k_0\frac{y_0}{x_0}\le r\left(t\right)\le \frac{y_0}{x_0}\left(t\in \left[0,T\right]\right),$ $r\left(0\right)=y_0/x_0$ and $r\left(T\right)=k_0{y}_0/x_0,$ we look for $\lambda \left(t\right)$ such that problem (4) has a solution $\left(x,y\right)$ such that

$$\frac{y\left(t\right)}{x\left(t\right)}=r\left(t\right),t\in \left[0,T\right].$$

With the notations from Section 3.1, from (5), one has

$$\begin{array}{ccc}\hfill v\left(t\right)-u\left(t\right)& =& v_0-u_0+{\int }_0^t\left(\frac{A}{1+B\left(e^u+e^v\right)}-\frac{a}{1+b_1{e}^u+b_2{e}^v}\right)ds\hfill \\ & & +\left(c-C\right){\int }_0^t\lambda \left(s\right)ds.\hfill \end{array}$$

Assuming that $C\ne c$ and taking into account that $v\left(t\right)-u\left(t\right)=ln\left(y\left(t\right)/x\left(t\right)\right),$ we find

$${\int }_0^t\lambda \left(s\right)ds=\frac{1}{c-C}\left(lnr\left(t\right)-lnr\left(0\right)+{\int }_0^t\left(\frac{a}{1+b_1{e}^u+b_2{e}^v}-\frac{A}{1+B\left(e^u+e^v\right)}\right)ds\right).$$

Hence,

$$\lambda \left(t\right)=\frac{1}{c-C}\left(\frac{r^{\prime }\left(t\right)}{r\left(t\right)}+\frac{a}{1+b_{1}x\left(t\right)+b_{2}y\left(t\right)}-\frac{A}{1+B\left(x\left(t\right)+y\left(t\right)\right)}\right).$$

This gives the expression of $\tilde{F}\left(x,y\right)$, and adding this into the integral equivalent version of (4) leads to a fixed point equation in $\left(x,y\right)$ which, as in Section 3.1, has a unique solution.

The stability of the solution and numerical simulations for this problem can be determined similarly.

4. Second Control Problem for the Normal–Leukemic System

Our second problem for the normal–leukemic cell dynamic system is to reduce the leukemic cell population to an acceptable level; this time, by controlling the proliferation rate of leukemic cells. Thus, we consider the system

$$\left\{\begin{array}{c}{x}^{\prime }=\frac{ax}{1+b_{1}x+b_{2}y}-cx\hfill \\ y^{\prime }=\frac{\lambda Ay}{1+B\left(x+y\right)}-Cy\hfill \\ x\left(0\right)=x_0,y\left(0\right)=y_0\hfill \end{array}\right.$$

(13)

where number $\lambda >0$ is the control parameter.

4.1. Solving of the Control Problem

Making the same change of variables as above yields the following system of integral equations:

$$\left\{\begin{array}{c}u\left(t\right)=u_0-ct+{\int }_0^t\frac{a}{1+b_1{e}^u+b_2{e}^v}ds\hfill \\ v\left(t\right)=v_0-Ct+\lambda {\int }_0^t\frac{A}{1+B\left(e^u+e^v\right)}ds\hfill \end{array}\right.$$

(14)

to which we associate the objective condition

$$v\left(T\right)=v_T,$$

(15)

where $v_T$ is the expected level of leukemic cells after a period of time $T.$

Compared with the abstract problem from Section 2, here

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

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

Applying the general principle given by Proposition 1, we obtain the following controllability result.

Theorem 3.

For each number $v_T$ with $max\left\{0,v_0-CT\right\}<v_T<v_0,$ the control problem (14) and (15) 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 $\left(u,v\right)\in {\Sigma }_1,$ that is $\left(u,v\right)\in C\left(\left[0,T\right],{\mathbb{R}}^2\right)$ solves (14) for some $\lambda \in {\mathbb{R}}_{+}.$ From the second equation of (14), we then have

$$\lambda =\frac{1}{{\int }_0^T\frac{A}{1+B(e^u+e^v)}ds}\left(v_T-v_0+CT\right),$$

(16)

hence, $F:{\Sigma }_1\to {\mathbb{R}}_{+}$ is given by

$$F\left(u,v\right)=\frac{1}{{\int }_0^T\frac{A}{1+B(e^u+e^v)}ds}\left(v_T-v_0+CT\right).$$

Clearly, we may extend F from ${\Sigma }_1$ to $C\left(\left[0,T\right],{\mathbb{R}}^2\right)$ by using the same expression, namely

$$\tilde{F}\left(u,v\right)=\frac{1}{{\int }_0^T\frac{A}{1+B(e^u+e^v)}ds}\left(v_T-v_0+CT\right)\mathrm{for}\left(u,v\right)\in C\left(\left[0,T\right],{\mathbb{R}}^2\right).$$

Now, our fixed point equation $w=\tilde{N}\left(w\right),$ with $\tilde{N}\left(w\right)=N\left(w,\tilde{F}\left(w\right)\right),$ becomes

$$\left\{\begin{array}{c}u=u_0-ct+{\int }_0^t\frac{a}{1+b_1{e}^u+b_2{e}^v}ds\hfill \\ \\ v=v_0-Ct+\frac{v_T-v_0+CT}{{\int }_0^T\frac{A}{1+B(e^u+e^v)}ds}{\int }_0^t\frac{A}{1+B\left(e^u+e^v\right)}ds.\hfill \end{array}\right.$$

Hence,

$$\begin{array}{ccc}\hfill {\tilde{N}}_1\left(u,v\right)\left(t\right)& =& u_0-ct+{\int }_0^t\frac{a}{1+b_1{e}^u+b_2{e}^v}ds,\hfill \\ & & \\ \hfill {\tilde{N}}_2\left(u,v\right)\left(t\right)& =& v_0-Ct+\frac{v_T-v_0+CT}{{\int }_0^T\frac{A}{1+B(e^u+e^v)}ds}{\int }_0^t\frac{A}{1+B\left(e^u+e^v\right)}ds.\hfill \end{array}$$

The existence of a fixed point w of $\tilde{N}$ can be guaranteed by using Schauder’s fixed point theorem. Indeed, the operator $\tilde{N}$ is completely continuous, as follows immediately from the Arzelà–Ascoli theorem. On the other hand, one has

$$\begin{array}{ccc}\hfill u_0-cT& \le & {\tilde{N}}_1\left(u,v\right)\left(t\right)\le u_0+aT,\hfill \\ \hfill v_0-CT& \le & {\tilde{N}}_2\left(u,v\right)\left(t\right)\le v_T+CT.\hfill \end{array}$$

Consequently, there are $R_1,R_2>0$ with

$${∥{\tilde{N}}_1\left(u,v\right)∥}_{\infty }\le R_1,{∥{\tilde{N}}_2\left(u,v\right)∥}_{\infty }\le R_2$$

for all $u,v\in C\left[0,T\right].$ Hence, $\tilde{N}\left(S\right)\subset S,$ where $S:=\{\left(u,v\right)\in C\left(\left[0,T\right],{\mathbb{R}}^2\right):$${∥u∥}_{\infty }\le R_1$, ${∥v∥}_{\infty }\le 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 $\tilde{N}.$ Finally, the control $\lambda$ is calculated using Formula (16), with u and v thus determined. □

Remark 3.

Since (13) is a Volterra system with Lipschitz nonlinearities, for any given $\lambda ,$ system (13) has a unique solution $\left(u,v\right).$ In particular, for the value λ 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 note in the following way: for a prescribed drug dose, the patient’s evolution is uniquely determined.

4.2. Stability of the Control Problem

Theorem 4.

Let $\left(u,v,\lambda \right)$ be a solution of the control problem (14) and (15). For a given $\overline{\lambda }\in {\mathbb{R}}_{+},$ denote by $\left(\overline{u},\overline{v}\right)$ the corresponding solution of system (14) and ${\overline{v}}_T:=\overline{v}\left(T\right).$ Then, for any small $\epsilon >0,$ one has

$$\left|\overline{\lambda }-\lambda \right|\le \frac{\epsilon }{AT}{e}^{-\left(a+\lambda A\right)T}=:\delta \mathit{implies}\left|{\overline{v}}_T-v_T\right|\le \epsilon .$$

(17)

Proof. 

As before, using the Lipschitz estimates, one has

$$\left|u\left(t\right)-\overline{u}\left(t\right)\right|\le a{\int }_0^t\left(\left|u\left(s\right)-\overline{u}\left(s\right)\right|+\left|v\left(t\right)-\overline{v}\left(t\right)\right|\right)ds,$$

$$\begin{array}{ccc}& & \left|v\left(t\right)-\overline{v}\left(t\right)\right|\hfill \\ & \le & \left|\lambda -\overline{\lambda }\right|{\int }_0^t\frac{A}{1+B\left(e^{\overline{u}}+e^{\overline{v}}\right)}ds+\lambda {\int }_0^t\left|\frac{A}{1+B\left(e^{\overline{u}}+e^{\overline{v}}\right)}-\frac{A}{1+B\left(e^u+e^v\right)}\right|ds\hfill \\ & \le & AT\delta +\lambda A{\int }_0^t\left(\left|u-\overline{u}\right|+\left|v-\overline{v}\right|\right)ds.\hfill \end{array}$$

Then, taking the sum

$$\begin{array}{ccc}\hfill \left|u\left(t\right)-\overline{u}\left(t\right)\right|+\left|v\left(t\right)-\overline{v}\left(t\right)\right|& \le & \left(a+\lambda A\right){\int }_0^t\left(\left|u\left(s\right)-\overline{u}\left(s\right)\right|+\left|v\left(t\right)-\overline{v}\left(t\right)\right|\right)ds\hfill \\ & & +AT\delta \hfill \end{array}$$

and Gronwall’s inequality yields

$$\left|u\left(t\right)-\overline{u}\left(t\right)\right|+\left|v\left(t\right)-\overline{v}\left(t\right)\right|\le AT\delta e^{\left(a+\lambda A\right)t}.$$

(18)

In view of the expression of $\delta ,$ this implies

$$\left|{\overline{v}}_T-v_T\right|\le \epsilon .$$

□

Remark 4.

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

5. Discussion

Model applicability: In the developmental hierarchy of bone marrow, whether normal or leukemic, there are four cell types (stem cells, progenitors, differentiated, and terminally differentiated). Our model aggregates over the properties of all cell types in the hierarchy to model a unique cell type for the normal phenotype and another unique cell type for the leukemic phenotype. Our abstract cells have regulatory negative feedback on the growth found at the stem cell level $1/(1+b_{1}x+b_{2}y)$ and $1/\left(1+Bx+By\right),$ respectively), the high growth rates of progenitor cells (large a and A), and the differentiated and terminally differentiated bulk mass effect, which eventually accounts for the disease dynamics.

Although our model was developed for chronic myelogenous leukemia, it applies to any leukemia type, acute or chronic, as long as it can be assumed that, over the considered timespan, all basic parameters of cell kinetics are time invariant. It does not apply to the natural course of chronic leukemia as a whole, from subclinical to blast crises, where the growth rate A increases by orders of magnitude with disease progression.

Moreover, our model can be easily extended to other physiological conditions. In our stability analysis, we have defined the satisfactory success criterion as $\left|\overline{r}-r\right|<\epsilon ,$ allowing for small deviances both above and beyond the set point. This relates our model to general models of tissue development and homeostasis, where two normal cell populations (instead of normal and leukemic) compete for space but are both desirable and are regulated by a single molecular signal (instead of chemotherapy).

First control problem: We divided up the treatment protocol into two phases. The first is a burn-in phase, from time $t=0$ to $T_0,$ in which the drug infusion only freezes the cellular status quo, keeping the malignant/normal cells ratio constant in time (at pretreatment value). In practice, this phase should be very short, as it exposes the patient to additional drug toxicity. It could be useful in clinical trial settings to assess the side effects of chemotherapy separately from the effects of the disease progression or to avoid tumor lysis syndrome at the onset of an induction cycle.

In the next phase, from $t=T_0$ to final $T,$ the tumor mass is reduced along a predefined curve, with its final point dependent on the subunitary parameter $k_0,$ and its steepness on the exponent $m.$

Theorem 2 explores the stability conditions of the control problem, namely establishing the maximum allowable deviation of the drug delivery rate from an ideal delivery schedule, such as if the actual endpoint of the treatment remains within the vicinity of the desired endpoint (only a small fraction of malignant cells remaining, $k_0{v}_0/u_0$). Crucially, we found the admissible error for the infusion function $\lambda \left(t\right)$ to decrease exponentially with both the growth rates a (of the normal cells) and A (of leukemic cells). This suggests that it may be much harder to achieve an ideal control of acute leukemia, compared to other tumor types, due to the higher division rates of leukemia and its normal cellular background. The admissible error only depends inversely on a linear combination of the cell decay rates c and $C,$ indicating that the rates of spontaneous cell death only marginally influence the control compared with the growth rates.

We found the same inverse exponential dependence on the total time treatment, $T.$ Acute leukemia cannot be equated with a “time-compressed” version of chronic leukemia, which behaves and is treated quite differently. Still, our results may suggest generally that the control of a slowly growing tumor, over a longer time, is equally hard to obtain as the control of an “acute” tumor over a shorter time span, insofar as the assumptions of our model hold.

In practice, the main sources of oscillations of the drug delivery function $\lambda \left(t\right)$ are the pharmacokinetic excursions of systemic drug levels, since chemotherapy is administered discontinuously. Our results establish a theoretical limit for how large these excursions can be, while still allowing for a cure (and implicitly a limit on the interdose intervals).

Second control problem: We obtained a similar exponential relation for the stability of the second control problem (Theorem 4), where the drug dose lambda is constant in time. The caveat here is that factor is $e^{-\left(a+\lambda A\right)T}$ instead of $e^{-\left(a+A\right)T}.$ This implies that, with faster leukemic cell production rates A and a subsequent need for higher drug levels $\lambda ,$ to achieve the same end result, the control is exponentially harder to obtain.