Research article Special Issues

Design of personalized cancer treatments by use of optimal control problems: The case of chronic myeloid leukemia

  • The advances in the mathematical explanation of the dynamics underlying treated cancer has opened the door to the mathematical design of optimal therapies. In parallel, the improvements and cost reductions in experimentation and data analysis techniques have made the formulation of personalized therapies possible. However, the design of cancer therapies making use of optimal control theory has not fully considered this possibility in detail. In this paper we contribute to the existing literature by analyzing the diverse alternatives that optimal therapy models offer to design personalized treatments. Taking as the starting point the Chronic Myeloid Leukemia (CML) optimal therapy model in [25], we design personalized optimal therapy models for patients with: CML; CML with intrinsic and/or induced resistance to the administered drug; CML and suffering high drug toxicity and/or allergy to the administered drug; and CML with presence of adverse factors. Along the paper we show that the clinical and medical applicability -the ultimate objective of this biomathematical research- of our proposed personalized models relies on the joint and proper use of the implemented calibration, simulation, and mathematical approaches and techniques. All the theoretical results generated by our personalized optimal therapy models are corroborated by clinical evidence.

    Citation: Pedro José Gutiérrez-Diez, Jose Russo. Design of personalized cancer treatments by use of optimal control problems: The case of chronic myeloid leukemia[J]. Mathematical Biosciences and Engineering, 2020, 17(5): 4773-4800. doi: 10.3934/mbe.2020261

    Related Papers:

    [1] Rafael Martínez-Fonseca, Cruz Vargas-De-León, Ramón Reyes-Carreto, Flaviano Godínez-Jaimes . Bayesian analysis of the effect of exosomes in a mouse xenograft model of chronic myeloid leukemia. Mathematical Biosciences and Engineering, 2023, 20(11): 19504-19526. doi: 10.3934/mbe.2023864
    [2] Natalia L. Komarova . Mathematical modeling of cyclic treatments of chronic myeloid leukemia. Mathematical Biosciences and Engineering, 2011, 8(2): 289-306. doi: 10.3934/mbe.2011.8.289
    [3] R. A. Everett, Y. Zhao, K. B. Flores, Yang Kuang . Data and implication based comparison of two chronic myeloid leukemia models. Mathematical Biosciences and Engineering, 2013, 10(5&6): 1501-1518. doi: 10.3934/mbe.2013.10.1501
    [4] Elena Fimmel, Yury S. Semenov, Alexander S. Bratus . On optimal and suboptimal treatment strategies for a mathematical model of leukemia. Mathematical Biosciences and Engineering, 2013, 10(1): 151-165. doi: 10.3934/mbe.2013.10.151
    [5] Katrine O. Bangsgaard, Morten Andersen, Vibe Skov, Lasse Kjær, Hans C. Hasselbalch, Johnny T. Ottesen . Dynamics of competing heterogeneous clones in blood cancers explains multiple observations - a mathematical modeling approach. Mathematical Biosciences and Engineering, 2020, 17(6): 7645-7670. doi: 10.3934/mbe.2020389
    [6] Dan Zhu, Qinfang Qian . Optimal switching time control of the hyperbaric oxygen therapy for a chronic wound. Mathematical Biosciences and Engineering, 2019, 16(6): 8290-8308. doi: 10.3934/mbe.2019419
    [7] Haifeng Zhang, Jinzhi Lei . Optimal treatment strategy of cancers with intratumor heterogeneity. Mathematical Biosciences and Engineering, 2022, 19(12): 13337-13373. doi: 10.3934/mbe.2022625
    [8] Ghassen Haddad, Amira Kebir, Nadia Raissi, Amira Bouhali, Slimane Ben Miled . Optimal control model of tumor treatment in the context of cancer stem cell. Mathematical Biosciences and Engineering, 2022, 19(5): 4627-4642. doi: 10.3934/mbe.2022214
    [9] B. M. Adams, H. T. Banks, Hee-Dae Kwon, Hien T. Tran . Dynamic Multidrug Therapies for HIV: Optimal and STI Control Approaches. Mathematical Biosciences and Engineering, 2004, 1(2): 223-241. doi: 10.3934/mbe.2004.1.223
    [10] Taeyong Lee, Adrianne L. Jenner, Peter S. Kim, Jeehyun Lee . Application of control theory in a delayed-infection and immune-evading oncolytic virotherapy. Mathematical Biosciences and Engineering, 2020, 17(3): 2361-2383. doi: 10.3934/mbe.2020126
  • The advances in the mathematical explanation of the dynamics underlying treated cancer has opened the door to the mathematical design of optimal therapies. In parallel, the improvements and cost reductions in experimentation and data analysis techniques have made the formulation of personalized therapies possible. However, the design of cancer therapies making use of optimal control theory has not fully considered this possibility in detail. In this paper we contribute to the existing literature by analyzing the diverse alternatives that optimal therapy models offer to design personalized treatments. Taking as the starting point the Chronic Myeloid Leukemia (CML) optimal therapy model in [25], we design personalized optimal therapy models for patients with: CML; CML with intrinsic and/or induced resistance to the administered drug; CML and suffering high drug toxicity and/or allergy to the administered drug; and CML with presence of adverse factors. Along the paper we show that the clinical and medical applicability -the ultimate objective of this biomathematical research- of our proposed personalized models relies on the joint and proper use of the implemented calibration, simulation, and mathematical approaches and techniques. All the theoretical results generated by our personalized optimal therapy models are corroborated by clinical evidence.


    After the work by Lotka ([1][3]) and Volterra ([4][6]), biomedical phenomena, characterized by the existence of complex dynamic interrelationships between the entities involved, began to be satisfactorily described by systems of differential/difference equations. The behavior of a disease under treatment is one of these biomedical situations susceptible to being explained through the use of a system of equations. There are two types of variables in this system. On the one hand, the state variables, whose trajectories are described by the difference/differential equations, thus being free magnitudes with unknown values until determined by the system. On the other hand, the control variables, whose values are not determined by the system but previously fixed, that influence the evolution of the state variables through the cross effects specified in the system of differential equations. It is clear that there exists a complete parallelism with the description of a treated disease: the control variables are the administered drugs, whose values and concentrations are decided by the physicians; while the state variables are those biomedical magnitudes of interest, whose behavior depends on both the administered drugs and the metabolic interrelationships described by the system of equations. Following this parallel, in the same way physicians manage drug doses to produce the best results in the patient, a control problem is able to deduce the optimal therapy according to a previously defined objective. In Biomathematics, this type of control problem is called an optimal therapy problem.

    When applied to cancer, if there are no mutations of the considered tumor cells, the system of equations describes the evolution of the mutually dependent number of tumor and normal cells in the course of a treatment with a drug. As is well known, cancer cells exhibit multiple genomic variations, both primary and derived from their considerable mutation capability (see for instance [7]), and each type has its therapeutic profile. As a result, in cancer treatment, physicians use drugs that are specific to a type of cancer or tumor -as opposed to radiation or non-selective chemotherapy-, the drug effectiveness being limited to the considered line of cancer cells. In the absence of mutations, it therefore becomes possible to mathematically identify the observed stable relationships between the existing normal cells, cancer cells and drug doses, opening up the formulation of an optimal therapy problem. Given that in this case the behavior of normal and cancer cells is a steady function of the administered drug dose, it becomes possible to govern the number of cancer and normal cells according to an objective by adequately manipulating the drug doses. From the biomedical perspective, this objective is clear: to eliminate the disease to the greatest extent possible while causing the minimum possible drug toxicity. This is the idea underlying a cancer optimal therapy model: to identify the therapy minimizing the damage to health caused by cancer and drug toxicity, damage measured and quantified by an objective function. As commented before, note that each type of cancer cell requires a specific drug and a particular optimal therapy problem, and that optimal therapies for a particular cancer cell do not work for other kinds of cancer cells because of distinct responses to the treatment, or even complete ineffectiveness of the drug. All these questions will be discussed and clarified along this paper.

    It is therefore clear that there are two dimensions to be considered when designing a cancer optimal therapy problem. First, the description of the dynamic behavior of the disease under any feasible treatment through the system of difference/differential equations. Second, the measurement of the malignancy of each specific status of the treated disease, that is, the formulation of the objective function. Regarding the first aspect, biomedical researchers and practitioners have provided information on the nature of the interrelationships between cancer and normal cells, their dynamics, and the value of the parameters involved in the arising system of difference/differential equations. Concerning the second, biologists and physicians have guided mathematicians in the process of congruously measuring the malignancy of treated cancer, that is, in formulating objective functions. The numerous textbooks and articles published during the last 15 years dealing with optimal therapy modeling reflect this collaboration and the subsequent advances and progress in this field (see for instance [8][10] or the textbooks [11][16]; for the specific case of cancer see [17][25] and the references provided by these authors).

    In parallel, as a result of the improvements and cost reductions in experimentation and data analysis techniques, today it is possible to obtain biomedical data specific to each patient and to formulate personalized therapies. Indeed, in clinical practice, personalized treatments are already heuristically applied, and Statistics and Mathematics have begun to explore this new research avenue, mainly related to genomic data. However, to the best of our knowledge, the design of cancer therapies making use of optimal control theory has not in detail fully considered this possibility. In this paper we contribute to the existing literature by analyzing the diverse alternatives that optimal therapy models offer to design personalized treatments. To this end, taking as the starting point the Chronic Myeloid Leukemia (CML) optimal therapy model in [25], we examine in detail the potential of cancer optimal therapy models to formulate personalized treatments, both through the modification of the system of equations describing the dynamics of the treated disease, and through the alteration of the objective function measuring malignancy.

    This paper is organized as follows. First, we describe the proposed optimal therapy problem for imatinib treated CML. After assigning referential values to the model parameters through the calibration of the model for the most general case, we discuss the distinct alternatives to design personalized therapies. More specifically, we consider the medical cases of patients resistant to imatinib, suffering abnormal toxicity or allergy, or presenting adverse prognostic factors. After the calibration and simulation of the personalized treatment models, we analyze their numerical results. In this respect, the comparison of the generated optimal therapies with the clinical evidence is particularly important. Physicians elaborate and decide treatments on the basis of heuristic procedures and continuous medical monitoring, and, therefore, this trial and error method must be envisaged as an approximation of the exact personalized optimal therapy arising from our optimal control problems. Finally, in section?? we comment on the main conclusions of the research.

    As commented in the introductory section, the first step in the mathematical modeling of this disease is the concretization of the particular considered case. In this respect, our proposed mathematical description is specific to a situation in which CML responds to imatinib treatment, not being applicable to other types of CML, which would require both alternative treatments and mathematical models. In particular and on the basis of work by [26] and [27], we therefore exclude the presence of mutations in the BCR-ABL1 gene, responsible for imatinib treatment failure, or, in the less unfavorable case, of a very high resistance to this drug that discourages imatinib treatment. Our benchmark model is that in [25], which describes the behavior of CML susceptible to imatinib treatment making use of a discrete time dynamic model. Although the most frequent biomedical mathematical models are formulated in continuous time, here we consider a discrete time framework. Since the biomedical literature on CML and clinical practice make use of variables expressed in per day values (drug doses, proliferation and decrease rates, carrying capacity, ...), by using a daily discrete dynamic formulation we allow for clearer comparisons and interpretations. In addition, given that the model involves the use of a system of non-linear equations with no algebraic solution, we avoid the discretization inherent to a continuous time framework.

    Our model of non-treated CML reproduces the biological interactions between the different types of normal and cancer cells observed in this disease. These interactions, described by [28], have been mathematically interpreted by [21], [22] and [25], to which we refer the interested reader. There are two different populations: that of hematopoietic stem cells (HSC), and that of differentiated cells (DC). When the disease is present, each of these populations is divided into normal cells and cancer cells. Then, denoting by T the treatment duration measured in days, at time instant t, t=0,1,,T+1, we distinguish between four different populations: normal HSC, denoted by x0(t); cancer HSC, denoted by y0(t); normal DC, represented by x1(t); and cancer DC, denoted by y1(t).

    The evolution over time of CML is described by a system of difference equations, which incorporates the most relevant biomedical facts. First, the populations of all types of cells naturally decrease at fairly constant rates. In this respect, let d0, g0, d and g be, respectively, the per day decrease rates of normal and cancer HSC, and normal and cancer DC. Cancer and normal DC are produced not only by the proliferation of DC but also by HSC, it being necessary to take into account these two mechanisms of increase. Let d2 and g2 be the per day rates at which normal and cancer DC proliferate and originate, respectively, normal and cancer DC; and let r and q denote the rates at which normal and cancer HSC produce normal and cancer DC, in this order. Through the self-renewal process, normal and cancer HSC produce similar cells by division, at rates n and m per day, respectively. As described by the biomedical literature ([28]), this self-renewing activity of HSC relies on a homeostatic process. Accordingly, the division of x0 is represented by a positive decreasing function Φ, depending on the total level of HSC (x0+y0), and given by

    Φ(x0+y0)=1x0+y0K, (2.1)

    where K represents the carrying capacity of bone marrow. Homeostasis for cancer HSC, y0, is governed by a positive decreasing function Ψ, which depends on (x0+αy0), where α(0,1] measures the fall in the homeostatic efficiency due to the disease (see [17] and [23] for an analysis of this fall),

    Ψ(x0+αy0)=1x0+αy0K. (2.2)

    Imatinib doses are described by a time-dependent positive sequence u(t), t=0,1,,T, which is bounded in [0,umax] due to the drug's toxicity. As in [25], drug effects are represented by the function h(u)=¯hln(u+1), satisfying: h(0)=0 (this means that cancer HSC declines at the baseline natural rate g0 without treatment); h(u) is non linear, since imatinib effects are non linear; and dh(u)du>0 and d2h(u)du2<0, given that higher doses imply higher effectiveness, these gains of effectiveness being decreasing as the dose increases. Since there exists empirical evidence suggesting a logarithmic relationship between biological responses and their causing factors (see [29]-[31]), this expression for h(u) appears as appropriate. In addition, and as shown in the next sections, it provides a good empirical fit to the clinical data on therapies. Given that imatinib is a targeted drug that only affects cancer HSC (see [32]), we limit our study to this biomedical situation. The referential behavior under treatment of CML is therefore given by the following system of difference equations:

    x0(t+1)=x0(t)+nΦ(x0(t)+y0(t))x0(t)d0x0(t), (2.3)
    x1(t+1)=x1(t)+rx0(t)dx1(t)+d2x1(t), (2.4)
    y0(t+1)=y0(t)+mΨ(x0(t)+αy0(t))y0(t)g0y0(t)h(u(t))y0(t), (2.5)
    y1(t+1)=y1(t)+qy0(t)gy1(t)+g2y1(t), (2.6)

    for t=0,1,,T, which, jointly with the initial values for normal and cancer HSC and DC x0(0), x1(0), y0(0), y1(0), describes the baseline dynamics of imatinib treated CML, to be conveniently modified according to the characteristics of each patient.

    To properly formulate the objective function measuring the malignancy of treated CML at each date, we take into account the fact that the aim of the treatment is to kill or limit the growth of cancer cells while keeping the drug toxicity as low as possible. Therefore, denoting the objective function to minimize by N, it must positively depend on the number of cancer cells and on the toxicity of the drug. In addition, these functions measuring malignancy adopt different specifications depending on the assumed malignancy trade-off between cancer HSC, cancer DC and the imatinib dose. We refer the interested reader to [17] for a detailed analysis of this question. Given its satisfactory performance in a wide variety of contexts, here we consider the objective function in [22] and [25]:

    N(u,y0,y1)=Tt=0[u2(t)+y20(t)+y21(t)]+[y20(T+1)+y21(T+1)]. (2.7)

    As explained above, the possible values of the imatinib dose are given by

    U={{u(t)}Tt=0|0u(t)umax,t=0,1,,T},

    where umax is the allowed maximum drug dose. The baseline optimal therapy problem is then:

    min{u(t)}Tt=0UN(u,y0,y1)=Tt=0[u2(t)+y20(t)+y21(t)]+[y20(T+1)+y21(T+1)]
    s.t.{x0(t+1)=x0(t)+nΦ(x0(t)+y0(t))x0(t)d0x0(t),x1(t+1)=x1(t)+rx0(t)dx1(t)+d2x1(t),y0(t+1)=y0(t)+mΨ(x0(t)+αy0(t))y0(t)g0y0(t)h(u(t))y0(t),y1(t+1)=y1(t)+qy0(t)gy1(t)+g2y1(t),x0(t)0,x1(t)0,y0(t)0,y1(t)0, (2.8)
    t=0,1,,T+1,x0(0),x1(0),y0(0),y1(0),initially given.

    The dynamics of the treated CML in this problem, described by equations (2.3)-(2.6), is that observed in the most general situation when CML is susceptible to imatinib treatment. In addition, the objective function (2.7) is also a universal measure of malignancy for treated CML. Then, by assigning to the parameters their most frequent or regular values, the reference problem (2.8) allows the optimal therapy to be determined for a representative patient. This optimal therapy, as well as the subsequent evolution of the cancer and normal cells, will be the reference to design, calibrate and compare personalized optimal therapies. This baseline optimal therapy problem is widely discussed in [25], so we remit the reader interested in going deeper on the solution, calibration, and simulation procedures to the said reference work. Here we simply collect the calibrated parameter values for this representative average patient in table 1.

    Table 1.  Calibrated values for the parameters in the baseline model.
    Parameter Description Value Units
    n Normal HSC division rate 0.00357 /day
    m Cancer HSC division rate 0.0037 /day
    d0 Normal HSC mortality rate 0.0005 /day
    g0 Cancer HSC mortality rate 0.0003 /day
    r Normal DC production rate 1011.5 /day
    q Cancer DC production rate 1011.5 /day
    d Normal DC mortality rate 1.25 /day
    d2 Normal DC proliferation rate 0.25 /day
    g Cancer DC mortality rate 1.1 /day
    g2 Cancer DC proliferation rate 0.5 /day
    K Carrying capacity of bone marrow 12791 HSC/day
    α Fall in homeostatic efficiency 0.1 dimensionless
    ¯h Constant in drug effect function h(u) 0.00455 dimensionless

     | Show Table
    DownLoad: CSV

    As predicted, once the baseline optimal therapy problem is solved, we confirm all the clinical salient features empirically observed for the considered patient. The therapy generated by the problem is specific to each patient and condition, and depending on the initial number of cancer cells, the model allows the optimal therapy to be calculated for any phase of CML (see [26] or [33]), namely chronic (less than 10% of cancer cells), accelerated (between 10% and 29%), and blast crisis (more than 30%). Since the most frequent case is CML at chronic phase (around 90% of patients) and we are interested in showing the versatility and applicability of our model, this will be the usual referential situation. In addition, BCR-ABL1 mutations have been not detected in the large majority of CML patients presenting chronic phase (see for instance [34]), so this seems to be the most suitable situation to design personalized optimal therapies with imatinib. Nevertheless, when necessary for illustrative purposes, the accelerated phase will also be considered as the initial situation of therapy.

    As is obvious from the mathematical formulation of our model, which describes the population of the four different types of cells in CML under imatinib treatment and computes the optimal therapy according to these variables, the only possible criterion to define the remission of the disease must be based on the resulting numbers of cancer and normal cells. This is indeed one of the clinical criteria followed to monitor patients and to define CML remission, but it is not the only one. As [26] and [35] describe, jointly with the hematological criterion of cell complete blood count adopted here, physicians also use cytogenetic and molecular tests. Today, the recommended criterion by the European LeukemiaNet 2020 is based on the quantification by PCR analysis of the BCR-ABL1 fusion gene, a molecular marker of CML. However, although this molecular test is the most sensitive assay available at present and leads to very accurate measures of the response to treatment and disease remission, the unknowns regarding BCR-ABL transcription and action mechanisms (see [36][42]) and the lack of the necessary volume of reliable and homogeneous data as a result of the novelty of this method (see [43] and [44]), currently make the formulation of a useful mathematical model explaining BCR-ABL1 levels impossible. Accordingly, we follow the hematological criterion (see [35] and [45]) to define the total remission of CML. This hematological analysis is accomplished by obtaining a complete blood count, and defines the eradication of CML as the non-existence of cancer HSC and DC, and the recovery of the number of normal HSC and DC in the blood, called complete hematological response (CHR). According to the clinical evidence (see for instance [46]), for CML at chronic phase, CHR is the most frequent result with around 95% of the cases. For our purposes, we will consider that CHR is reached when the number of cancer HSC is lower than 1, the number of cancer DC is undetectable and then lower than e-4 of the total number of DC (see [47][51]), and the number of normal HSC and DC are those in a healthy individual, given by 11000 and 3.48e+15 (see [28]). For this representative patient and according to clinical evidence, this happens before the third month of treatment (see [26], [35], [52][55]). Finally, according to the quoted references, the standard recommended treatments consist of daily doses of imatinib ranging from 100 mg to 400 mg, and not exceeding 800 mg.

    In this respect, assuming initial values of cancer HSC and DC y0(0)=10 and y1(0)=5.2e12, representative of CML at chronic phase (see [26] or [33]), we replicate the therapy and the disease evolution observed in clinical practice for the most general case. More specifically, our optimal therapy entails: an average dose of 407 mg/day; a treatment length (with doses higher than 1 mg) of 1331 days, that is about 3 years and 8 months; and CHR attained at day 77. In this respect, we solve the problems identified by [56] for models based on [28], namely the unrealistic values for the parameters necessary to replicate the observed clinical data. These results substantiate the use of this baseline optimal therapy problem (2.8) as the reference to design personalized therapies by modifying and adjusting functions and parameters to the particular situation of each patient. These personalized therapies will be the subject of the next sections. For further analyses of this baseline model, we refer the reader to [25].

    With respect to the typical behavior of treated CML, the clinical evidence shows that patients can develop resistance to imatinib in a low percentage of cases. Clinically, it is considered that resistance appears when patients do not achieve CHR under standard imatinib treatment after 3 months ([26], [52], [53], [55]). More specifically, as the eight-year International Randomized Study of Interferon and STI571 (IRIS) trial concluded (see [57]), resistence to imatinib appears for around 6%-25% of patients (see also [53] and [58]). Biomedical studies also establish that there exist two types of resistance to imatinib, namely intrinsic and induced ([46], [59], [60]). Primitive resistance to imatinib, known as de novo or intrinsic resistance, previously exists in the patient, while primary cytogenetic resistance or induced resistance, the more common resistance to imatinib, is developed by the patient once imatinib has been administered. These two types of resistance manifest through both BCR-ABL1 dependent and independent mechanisms ([26], [27], [42] and [53]). As commented in section 2.1, mutations in the BCR-ABL1 gene cause the lack of substantial and durable response to imatinib and treatment failure, a situation excluded in our model, and then our optimal therapy problem only applies to BCR-ABL1 independent imatinib resistance, susceptible to treatment with imatinib. In this case, clinical evidence shows that this resistance is solved by increasing imatinib doses within the tolerable range ([59]) until the complete removal of the disease. This is a quite frequent situation: According to the IRIS trial, in those patients who did not achieve CHR after 3 months, a high percentage (86%) attained CHR after dose escalation 3 months later.

    The Biomathematical literature on CML, for its part, has introduced cancer resistance to imatinib by assuming the appearance of a subpopulation of cancer cells resistant to the administered drug. This is the case of, among others, [22], [28] and [61]. This is obviously a pertinent approach to analyze resistance to drugs caused by the mutation of cancer cells and the appearance of a subpopulation of cancer cells totally resistant to the administered drug, and is probably a suitable proposal to tackle the effects associated to the presence of mutations in the BCR-ABL1 gene. Nevertheless and as pointed out before, obtaining realistic optimal therapies for this case requires a more perfect numerical knowledge of the biomedical characteristics of such a subpopulation, mainly clarifying how they originated and how they are related with the other types of cell, knowledge that does not exist today. Our objective, conditioned by the existence of reliable data on the biomedical mechanisms underlying our mathematical formulation, is to explain the resistance to imatinib in CML that can be solved by increasing drug doses. This fact suggests that the considered resistance to imatinib is not related to the appearance of a resistant subpopulation of cancer cells, but to an intrinsic or induced decrease in the effects of imatinib. In this respect, let hR(u) denote the new function that captures the drug effects when resistance is present. Taking as the starting point our baseline model, the following two subsections formulate and simulate personalized optimal therapies able to consider the two types of imatinib resistance, obtained by solving the optimal problem under resistance

    min{u(t)}Tt=0UN(u,y0,y1)=Tt=0[u2(t)+y20(t)+y21(t)]+[y20(T+1)+y21(T+1)]
    s.t.{x0(t+1)=x0(t)+nΦ(x0(t)+y0(t))x0(t)d0x0(t),x1(t+1)=x1(t)+rx0(t)dx1(t)+d2x1(t),y0(t+1)=y0(t)+mΨ(x0(t)+αy0(t))y0(t)g0y0(t)hR(u(t))y0(t),y1(t+1)=y1(t)+qy0(t)gy1(t)+g2y1(t),x0(t)0,x1(t)0,y0(t)0,y1(t)0, (3.1)
    t=0,1,,T+1,x0(0),x1(0),y0(0),y1(0),initially given.

    This type of resistance is previous to, and independent of, imatinib administration. According to the literature ([61]), in the absence of BCR-ABL1 mutations, this de novo resistance is due to diverse factors already present in the patient before treatment with imatinib, such as increased expression of P-glycoprotein, decreased expression of the drug uptake transporter human organic cation transporter 1 (hOCT1), sequestration of imatinib by increased serum protein α acid glycoprotein, or alternative signaling pathway activation. In any case, these factors result in a decrease in the drug effectiveness not related to the administered drug doses. Here we propose a simple modification of the referential optimal therapy problem to take this type of resistance into account. More specifically, we modify the baseline function measuring the drug's effectiveness h(u) by introducing a constant that captures the existence of resistance to imatinib and the constant decrease in its effects, according to the expression

    hR(u)=A¯hln(1+u),

    where 0<A<1 measures the patient's idiosyncratic resistance to imatinib. Obviously, after this substitution in our personalized optimal therapy model (3.1), it is necessary to calibrate the new parameter A for the considered patient. In this respect, our approach makes the valuation of the resistance of each individual possible, something that constitutes an important additional advantage of our proposal.

    As explained above, resistance to imatinib appears when the patient does not achieve a complete hematological response under standard imatinib treatment after 3 months. Our quantitative measure of resistance is implicit in this definition. Let us explain how this personalized calibration of A is possible, taking as example an individual identical to our representative patient except for the presence of intrinsic resistance. In a standard treatment of CML at chronic phase, an imatinib dose of 400 is administered daily, and a blood cell count is carried out the third month of treatment ([26], [35], [52]-[55]). Therefore, we know the administered dose and the number of cancer cells at the day of counting, provided by the physicians. Now, with the parameter values in table 1, if we run the system of difference equations in the personalized problem (3.1) for the administered imatinib concentration and different values of A, we can find the value of the resistance parameter A implying the observed numbers of cancer cells at the day of counting. Without any loss of generality, let us assume that the administered dose is 400 mg/day, and that the blood cell count is carried out in the 3rd month of treatment. The physicians find that, at day 90, the number of cancer HSC is between 1 and 2 and the number of cancer DC is around 6.5e10. According to the medical criteria, this patient presents resistance to imatinib, since after 3 months (90 days), CHR is not achieved: there are more than 1 cancer HSC, and the number of cancer DC is greater than e-4 of the total number, given by 3.4e15. In addition, the physicians can compare these numbers with those for the typical patient with no resistance by solving the baseline problem, and they can identify the existence of resistance from an alternative perspective. Indeed, for the patient free of resistance, the number of cancer HSC at day 90 according to the baseline optimal therapy is lower than 1, and the existence of resistance can therefore be concluded. If we assume that the parameter values for this patient are those in table 1, and we run equations (23)-(27) for A=0.87, we obtain a close approximation to these numbers (1.22 and 6.71e10, respectively) and then we can conclude that A=0.87 is a quantitative measure of the imatinib resistance for that patient. Note that, if our individual is not our representative patient, the calibration procedure is exactly the same after substituting values in table 1 for those obtained for this particular patient. In any case, our approach allows this intrinsic resistance to imatinib to be measured for any patient, something that constitutes an important application of the proposed model. Indeed, our model can also be applied to ascertain, for each resistant patient, his/her degree of intrinsic imatinib resistance.

    Once the resistance to imatinib is quantified for the considered patient, the subsequent personalized optimal therapy can be obtained by solving and simulating the optimal control problem (3.1). The personalized treatment, as well as its consequences, can then be compared to those attained for the representative patient (with no resistance), as well as to those observed in clinical practice. As an illustrative example, let us consider the baseline results and a patient with de novo resistance to imatinib given by A=0.87. Since the parameter values and characteristics are those of the average patient, this case, as our baseline model, can also be interpreted as a referential description of intrinsic resistance to imatinib. Indeed, as figures 1-5 show, we confirm all the regular clinical practices when resistance to imatinib is detected for a patient and that are heuristically applied by the physicians.

    Figure 1.  Personalized optimal therapies with (blue) and without (red) de novo resistance.
    Figure 2.  Difference in cancer HSC, optimal therapy with resistance minus optimal therapy without resistance.
    Figure 3.  Difference in cancer DC, optimal therapy with resistance minus optimal therapy without resistance.
    Figure 4.  Difference in normal HSC, optimal therapy with resistance minus optimal therapy without resistance.
    Figure 5.  Difference in normal DC, optimal therapy with resistance minus optimal therapy without resistance.

    More specifically, as a consequence of resistance, an increase in the optimal imatinib dose appears, which is able to suppress the disease as observed in clinical practice. The biomedical interpretation of these modifications in the optimal drug dose is clear: resistance to imatinib implies a higher number of cancer cells for a given drug dose, and then the necessity of higher imatinib doses in order to reduce cancer cells. With respect to the typical patient with CML at chronic phase, for this patient with resistance to imatinib, the maximum dose is administered for 107 more days, and the personalized optimal therapy entails an average daily dose of 469 mg, that is 63 additional mg. of imatinib per day. This personalized treatment is that depicted in figure 1, and would ensure the total disappearance of chronic phase-CML between days 89 and 106, i.e., a maximum of 29 days later than the equivalent patient with no resistance. Indeed, at day 89, the number of cancer HSC is lower than 1, and the number of cancer DC is lower than 3.4e-11 and undetectable at day 106. In addition, the number of normal HSC and DC are practically those in the healthy individual after day 90. It is worth noting the optimal nature of this therapy: when the drug doses are calculated according to the optimal problem (3.1), the number of cancer HSC is lower than a week and a half before than when the equivalent continuous dose of 469 mg/day is administered. Figures 2-5 depict the difference in cancer and normal cells implied by the personalized optimal therapy applied after detecting intrinsic resistance, computed as cells for the optimal therapy under resistance minus cells for the baseline therapy. In any case and as explained before, these differences are irrelevant from the clinical point of view: cancer and normal HSC differ in less than a cell, and the differences between cancer and normal DC are below the detectable value.

    The biomedical literature has also identified an induced resistance to imatinib, i.e., a decrease in the drug's effectiveness caused by the administration of imatinib and the accommodation of the disease to imatinib (see [59]-[63]). Biomedical evidence also shows that, if mutations in the BCR-ABL1 gene are not present, this induced resistance directly depends on the number of days with imatinib therapy, but not on drug doses, the recommended modification in therapy being a sustained increase in the drug doses ([58], [62] and [64]). These empirical data suggest that induced resistance is increasing in time and independent of imatinib doses, and must be modelled in a different way than intrinsic resistance. In this respect, the following modification of the function measuring the drug's effectiveness h(u) is proposed:

    hR(u(t))=h(u(t))Aln(t+1),

    where A>0. According to this formulation, the term Aln(t+1) measures the patient's induced resistance to imatinib. This resistance depends on idiosyncratic factors specific to the patient, captured by A, and also positively on the treatment elapsed time t. Through the logarithm function, we captures a likely feature of induced resistance: resistance increases as time passes, but these increases are progressively lower. Once this modification has been introduced in the optimal problem under resistance (3.1), the procedures to quantitatively measure this induced resistance and to obtain, simulate and compare the personalized therapy and its results, are the same as those explained for the intrinsic resistance.

    As described in the previous section, we count on data concerning the administered doses of imatinib and the number of cancer cells at the day of counting, provided by the physicians. With the biomedical parameter values of the considered patient, and making use of the system of difference equations describing the disease under induced resistance (that is, after substituting hR(u) in system (2.3)-(2.6), the value of A that best fits the observed cell count can be found. In this respect, and in order to gain insights into the differences between the two type of resistance, let us consider the same situation discussed for the intrinsic resistance: the physicians find that, at day 90, the patient presents resistance to imatinib, since the number of cancer HSC is between 1 and 2 and the number of cancer DC is around 6.2e10. Assuming that the parameter values for this patient are those in table 1, and after running the system of equations describing the dynamics of treated CML under induced resistance, we find that, for A=0.345, the number of cancer HSC and DC at day 90 are very similar to those provided by the physicians: 1.12 and 6.07e10. Therefore, if induced resistance is identified for this patient, it can be concluded that A=0.35 is a quantitative measure of his/her induced imatinib resistance. Note again that, if our individual is not our representative patient, the calibration procedure is exactly the same after substituting values in table 1 for those obtained for that particular patient. As for intrinsic resistance, our approach allows induced resistance to imatinib to be measured for any patient. In addition and as we will analyze in the following paragraphs, it has another relevant application, since it enables physycians to identify which type of resistance the patient presents.

    Once the induced resistance to imatinib has been quantified for the considered patient, the subsequent personalized optimal therapy can be obtained. This personalized treatment and its consequences are susceptible to a triple comparison: they can be compared to those attained for the patient without resistance, to those for the patient with intrinsic resistance, and to those observed in clinical practice when resistance is detected. The comparison with the baseline patient is qualitatively similar to that obtained for intrinsic resistance, so the most interesting is that established with the same patient assuming equivalent intrinsic resistance, depicted in figures 6-10.

    Figure 6.  Personalized optimal therapies, induced (red) and intrinsic (blue) resistance.
    Figure 7.  Difference in cancer HSC, induced resistance under optimal therapy for intrinsic resistance.
    Figure 8.  Difference in cancer DC, induced resistance under optimal therapy for intrinsic resistance.
    Figure 9.  Difference in normal HSC, induced resistance under optimal therapy for intrinsic resistance.
    Figure 10.  Difference in normal DC, induced resistance under optimal therapy for intrinsic resistance.

    As figure 6 shows, although for both patients a similar failure of CHR is clinically observed, the optimal therapies are very different depending on whether resistance to imatinib is induced (in blue) or intrinsic (in red). As we know, at day 90, the counted numbers of cancer HSC and DC were the same in both cases. However, the behavior of the disease clearly differs in the two types of resistance. With intrinsic resistance, independent of the previous administration of imatinib, the loss of the drug's effectiveness is constant along the period of treatment. Nevertheless, if resistance is induced, imatinib resistance, as well as the loss in the drug's effectiveness, increase as time passes. Then the number of cancer HSC and DC cells are greater as resistance rises, and higher drug doses are necessary in comparison with intrinsic constant resistance, as depicted in figure 6. Clinically, the observed results are also confirmed. As documented in [58], [62] and [64], when induced resistance that is independent of BCR-ABL1 mechanisms exists, higher doses of imatinib are able to completely remove the disease. Since the number of cancer HSC and DC are greater than for intrinsic resistance, additional increases in drug doses are necessary to eliminate the disease with respect to the first type of resistance. For our specific patient with induced resistance, the optimal therapy implies an average daily dose of 797 mg, which is 327 more mg per day, and after administering these higher imatinib doses, CHR is attained between days 84 (cancer HSC lower than 1) and 108 (cancer DC lower than 3.4e-11). Note that, concerning the number of cancer cells and the objective of achieving CHR, the effectiveness of this optimal therapy under induced resistance is virtually equal to that under intrinsic resistance. The important difference is the doses of administered drug, much higher with induced resistance. This is a logical consequence of the relatively greater number of cancer HSC when induced resistance exists. To see this, it is interesting to replicate the behavior of the disease for one type resistance case when the administered optimal therapy is that corresponding to the other type of resistance, since this replication constitutes a criterion to elucidate the type of resistance that the patient presents.

    As depicted in figures 7-10, if the type of resistance is misidentified by the physicians, the subsequent wrong personalized optimal therapy entails a higher number of some type of cancer cells and/or a lower number of normal cells, as well as an altered behavior of the disease. Here, to gain insights, we have considered the patient described in the former sections. Let us assume that, although this individual actually presents induced resistance, the physicians interpret that there exists intrinsic resistance. The personalized therapy obtained for this intrinsic case and computed in section 3.1 is therefore applied. As a result of this misidentification, in comparison with the evolution of the disease if the right type of resistance (induced) had been identified, a substantial decrease in the number of normal DC is observed, represented in figure 10. In addition, the logical dynamic relationships between cancer and normal cell disappear. In this case, during the first days of the wrong therapy, the observed number of cancer HSC is lower than expected, but then it slightly increases over the predicted value. Concerning the number of cancer DC, it is always under the expected value, so, in principle, the wrong identification of the resistance type could seem not to be excessively important. However, although the number of normal HSC is greater than expected, there is a dramatic decrease in the number of normal DC, persistent over time and taking very high values.

    In this respect, it is worth noting that these theoretical effects associated to the misidentification of the actual type of resistance constitute a guide to correctly determining it. In fact, if the physicians wrongly identify the type of resistance, they will find not only anomalous numbers of cancer and normal HSC and DC cells, especially the latter, but also altered dynamics. These are clues suggesting the need to reconsider the diagnostic and redesign the personalized optimal therapy.

    As explained in the previous sections, when intrinsic or induced resistance to imatinib exists for our considered patient, CHR is attained 29 and 30 days later than for the baseline case with no resistance, respectively at days 106 and 107. Indeed, after day 200, the difference between the numbers of the same type of cells are undetectable, so we restrict the graphic representation to the first 200 days. Moreover, for normal cells, the evolution is so similar than the representation does not provide any useful information. The specific dynamics for cancer cells in these 3 situations are those represented in figures 11-12. As depicted, although along the first days of therapy there exist some relevant differences, they vanish gradually, being undetectable after day 300.

    Figure 11.  Cancer HSC evolution, baseline (yellow), intrinsic resistance (blue), induced resistance (red).
    Figure 12.  Cancer DC evolution, baseline (yellow), intrinsic resistance (blue), induced resistance (red).

    As with every drug, the administration of imatinib presents adverse effects. Some patients may experience lethargy, dermatological allergies, hepato-toxicity, and vascular toxicity, as well as bone pain and cytopenias. In general, as documented by [65]-[69], all these side effects are reversible on diminishing the drug doses, and the lower administered dosages are able to completely remove the disease. The existence of allergy, which is the prejudicial reaction of the patient's immune system to imatinib, as well as of the drug's toxicity, defined as the occurrence of damages to the organism of the patient caused by imatinib, can be incorporated in our baseline model in order to design personalized therapy models. To contemplate these situations, the formal element to be modified is not the system of equations (2.3)-(2.6) describing the dynamics of treated CML, but the expression of the objective function (2.7). In biomedical terms, if for a patient these deleterious effects associated to the administration of imatinib are observed, the relative malignancy of the drug relative to that of cancer cells must increase, and we must simply increase the weight of u(t) in the objective function keeping all the remaining weights (for y0 and y1) constant. From our objective function (2.7) in the baseline problem, the former reasonings lead to the following proposed modification:

    N(u,y0,y1)=Tt=0[Bu2(t)+y20(t)+y21(t)]+[y20(T+1)+y21(T+1)], (4.1)

    where B>1 is the constant allowing the toxicity to be modulated. However, unlike resistance to imatinib, this allergy or toxicity degree cannot objectively be measured, since it does not depend on variables included in the model but on external appreciations made by the physicians. As we will explain in the next paragraphs, the solution is to calibrate the additional weight for u(t) under the restriction of imposing, under the medical criterion, a limit to the increase in the number of cancer cells.

    As for the former personalized therapies and without any loss of generality, let us assume a patient showing the baseline parameters in table 1. If this patient develops allergy to imatinib, or if imatinib implies a higher than usual toxicity for this individual, the malignancy of a cancer cell is now equivalent to a lower drug dose, and it becomes necessary to decide to what extent imatinib decreases are beneficial in global terms. Note that this question, relative to the calibration of the toxicity of imatinib, is of a different nature than the calibration of resistance to imatinib, since the latter is objectively measured from the counted number of cancer cells, a possibility that now does not exist. As commented on before, the solution passes through the subjective medical criterion of imposing a limit to the increases in the number of cancer HSC and DC that appear as a consequence of lower imatinib doses with respect to the baseline referential situation. For illustrative purposes, here we will consider a patient with CML at accelerated phase, with a much higher number of cancer cells, since it allows these increases to be better visualized and interpreted. More specifically, the number of cancer HSC and DC are, respectively, y0(0)=103 and y1(0)=1010. In this case, if the physicians detect allergies or toxicity and decide to decrease drug doses allowing a maximum increase in cancer DC of e5, the associated value for B is B=106. Analogously, if the drug's toxicity is greater and the medical criterion is to additionally decrease the imatinib dose allowing for a higher increase in cancer DC of e8, the subsequent value for B is B=1012. The biomedical and mathematical analyses and concepts concerning this trade-off between cancer malignancy and drug toxicity are exhaustively discussed in [17], to which we refer the interested reader. Here we present a graphic representation of how toxicity is calibrated in 13 to 17, which will be commented in the following paragraphs.

    As observed in clinical practice, the appearance of toxicity and/or allergies after imatinib administration leads to personalized optimal therapies entailing lower drug doses. In this respect, for our CML patient at accelerated phase, the two considered levels of imatinib toxicity, namely high (B=1012) and low (B=106), imply personalized optimal therapies with substantial reductions in drug doses. Under low toxicity, the maximum dose is administered 225 days less than in the baseline situation with no toxicity, and the average per day dose is 124 mg. lower. When toxicity is high, the maximum dose is given 450 days less, and the reduction in the average per day dose is 248 mg. lower. Figure 13 depicts the reference, high and low toxicity optimal therapies, and clearly shows how, under imatinib toxicity, the personalized optimal therapies entail decreases in drug doses.

    Figure 13.  Personalized optimal therapies with high (red) and low (yelow) imatinib toxicity, compared to the baseline case (blue).

    As a consequence of the lower imatinib doses, there appear increases in the number of cancer cells with respect to the non-toxicity case. As explained above, these relative increases allow the physicians to calibrate the toxicity parameter. Logically, the higher the toxicity, the higher the increases in cancer HSC and DC allowed by the physycians, as figures 14 to 17 show. More specifically, when the toxicity is low, the maximum increases in the number of cancer HSC and DC with respect to the baseline case are 2.1e-7 and 1.1e5, respectively. If the toxicity is high, these maximum allowed increases are 2.1e-4 and 1.1e+08, in the same order. In any case, it is worth noting that these lower doses of imatinib are able to control the disease: as happens in clinical practice, when toxicity is detected and lower drug dosages are administered, these doses are sufficient to attain CHR. Indeed, our personalized optimal therapies imply the complete recovery of the patient suffering imatinib toxicity a few days after the CHR for the baseline situation, and the total control of the disease. Indeed, as represented in figures 14-17, even where the maximum increase in cancer cells is attained (at day 473 for high toxicity and day for 699 for low toxicity), the cancer cells are undetectable and do not imply the infringement of CHR criteria.

    Figure 14.  Increase in cancer HSC, optimal therapy with high toxicity minus baseline optimal therapy.
    Figure 15.  Increase in cancer HSC, optimal therapy with low toxicity minus baseline optimal therapy.
    Figure 16.  Increase in cancer DC, optimal therapy with high toxicity minus baseline optimal therapy.
    Figure 17.  Increase in cancer DC, optimal therapy with low toxicity minus baseline optimal therapy.

    The therapeutical success of imatinib is evident. Jointly with its character of targeted drug and its high effectiveness, the reversibility and manageability of its adverse side effects make imatinib the most suitable drug for CML treatment. Nevertheless, and as happens for every drug, there exist some adverse factors that diminish the probability of achieving CHR for some patients. As explained in the introductory section, the main adverse factor is the presence of mutations in the BCR-ABL1 fusion gene, responsible for imatinib treatment failure or for a very high resistance to this drug that discourages imatinib treatment. In these cases, given the lack of response to imatinib, second (nilotinib, dasatinib and bosutinib) and third (ponatinib) generation drugs have demonstrated to be the most suitable treatments. As our model is appropriate to describe CML susceptible to effective imatinib treatment, we exclude this situation, which would require a different formulation and calibration. In this respect, the biomedical evidence shows that the presence of damaged areas in bones, age of 60 years and older, multiple chromosome changes in cancer HSC, and especially enlarged spleen and an increased number of basophils and eosinophils, constitute adverse prognosis factors in CML that do not totally block imatinib effectiveness. We refer the interested reader to the Sokal, Euro and Euto risk scores ([70]) and to the numerous studies on these predictive indexes (see [71]-[74]).

    For our purposes, the presence of these adverse factors for a patient implies the need for a modification of the baseline optimal therapy problem in order to obtain her/his personalized therapy. In principle, these changes are not associated with losses in the drug's effectiveness or with the emergence of toxicity side effects. The most plausible biomedical explanation is that they are the consequence of the appearance of extramedullary hematopoiesis associated to certain characteristics of cancer HSC (see [73] and [75]). This implies that this type of cancer cells are not only produced in bone marrow. Indeed, the presence of an enlarged spleen and higher than usual numbers of basophils and eosinophils, the best predictors of CML remission, suggest that cancer HSC proliferate, differentiate, and mature in the spleen under these circumstances. To contemplate this fact and to design the appropriate personalized optimal therapy, the modification of the baseline problem requires the introduction of an additional rate of increase for HSC, that coming from the altered spleen. More specifically, equation (2.5) becomes

    y0(t+1)=(1+s)y0(t)+mΨ(x0(t)+αy0(t))y0(t)g0y0(t)h(u(t))y0(t), (5.1)

    where s>0 is the proliferation rate in the spleen of cancer HSC. As we know, this new parameter must be calibrated in order to compute the personalized optimal therapy.

    As explained above, the presence of adverse factors translates into the appearance of an additional proliferation rate for cancer HSC. As happened with intrinsic resistance, our quantitative measure of this growth rate is implicit in its definition, and in fact, the calibration procedure is identical. Since we know the administered dose and the number of cancer cells at the day of counting, provided by the physicians, if we run the system of difference equations describing the dynamics of the treated disease under adverse factors

    x0(t+1)=x0(t)+nΦ(x0(t)+y0(t))x0(t)d0x0(t),x1(t+1)=x1(t)+rx0(t)dx1(t)+d2x1(t),y0(t+1)=(1+s)y0(t)+mΨ(x0(t)+αy0(t))y0(t)g0y0(t)h(u(t))y0(t),y1(t+1)=y1(t)+qy0(t)gy1(t)+g2y1(t),

    for the administered imatinib concentration and different values of s, we can calibrate the value of the spleen proliferation rate s by looking for the best replication of the observed numbers of cancer cells at the day of counting.

    Without any loss of generality, let us consider again the patient with CML at chronic phase taken as reference to calibrate the parameters measuring resistance and toxicity. The administered imatinib dose is 400 mg/day, and the blood cell count, carried out in 3rd month of treatment, provides a number of cancer HSC and DC between 1 and 2 and a number of cancer DC around 6.2e10. Since after 90 days CHR is not achieved, two possible scenarios are considered by the physicians: resistance to imatinib, or the presence of adverse factors. Further medical analyses allows the latter to be concluded (for instance after identifying the existence of an enlarged spleen), and biomathematicians find that the above system replicates these numbers when s=0.004 (cancer HSC 1.11 and cancer DC 6.1e10 at day 90). Once this in spleen-proliferation rate of cancer HSC has been calibrated, the following step is to compute the personalized optimal therapy for this patient presenting adverse factors.

    Once the consequences of the presence of adverse factors have been quantified for the considered patient, the subsequent personalized optimal therapy can be obtained by solving and simulating the associated optimal control problem. The personalized treatment, as well as its effects, can then be compared to those attained for the representative patient without adverse factors, as well as to those observed in clinical practice. Once again, let us consider the baseline results and a patient with adverse factors leading to an additional in-spleen growth rate for cancer HSC s=0.004. As with the previous personalized therapies, this case can be interpreted as a referential description of the personalized treatment under the presence of adverse factors. Indeed, we confirm the regular clinical practice when these factors are detected, namely the increase of imatinib doses heuristically applied by the physicians (see [76]).

    With respect to the no adverse factors situation, this increase in imatinib doses entails an administration of the maximum dose of 800mg/day for 185 more days, and an average daily dose higher in 10 mg. Total remission of the disease appears between days 84 (cancer HSC lower than 1) and 100 (cancer DC lower than 3.4e11, that is days between 7 and 23 days later than for the baseline case. Once the personalized optimal therapy for this patient with adverse factors is administered, the dynamics of the disease is very similar to that obtained when there exists intrinsic resistance. Thus, we only represent here the baseline and personalized optimal therapies between days 400 and 1200 (figure 18). Differences in cancer and normal cells are completely analogous to those displayed in figures 2-5.

    Figure 18.  Personalized optimal therapies with (blue) and without (red) adverse factors.

    The advances in the mathematical explanation of the dynamics underlying treated cancer have opened the door to the mathematical design of optimal therapies. Indeed, the application of optimal control theory to fight cancer is an important research avenue in today's Biomathematics. Two are the dimensions to be considered when designing a cancer optimal therapy problem. First, the description through a system of difference/differential equations of the dynamic behavior of the disease under treatment. Second, the measurement of the malignancy of each specific status of the treated cancer, that is, the formulation of the objective function. By solving the subsequent optimal control problem, the therapy providing the drug doses that minimize the damage to health caused by the cancer and the drug toxicity, known as optimal therapy, is identified. Given the particularities of cancer, whose behavior is very specific to each of the numerous associated molecular and genetic changes that are possible in this disease, the formulated optimal therapy problem is only susceptible to application for the specific considered situation. Indeed, from a strict clinical perspective, there exist as many cancer types as cancer patients, each one requiring its own and specific mathematical formulation. Not surprisingly, this has led to an increasing interest in the mathematical analysis and design of personalized cancer optimal therapies. Separately and in parallel, the improvements and cost reductions in experimentation and data analysis techniques have made the formulation of personalized therapies possible by providing reliable values for the biomedical parameters involved. As a result, in clinical practice, personalized treatments are already heuristically applied, and Statistics and Mathematics have begun to explore this new research avenue, mainly related to genomic data. In this paper, we contribute to the existing literature on the mathematical design of personalized therapies by analyzing the diverse alternatives that optimal therapy models offer to design personalized treatments. To this end, taking as the starting point the Chronic Myeloid Leukemia (CML) optimal therapy model in [25], we examine in detail the potential of cancer optimal therapy models to formulate personalized treatments with imatinib.

    As explained above and in the previous sections, this mathematical model only applies to the specific variant of CML susceptible to imatinib treatment, and therefore excludes patients presenting mutations in the BCR-ABL1 fusion gene, responsible for imatinib treatment failure. In particular, through the modification of the system of equations describing the dynamics of the treated disease, we design optimal therapies for patients with CML affected by resistance to imatinib. Two are the proposed modifications, which correspond to the two types of resistance identified by physicians, namely intrinsic and induced. Interestingly, we show that the subsequent personalized optimal therapy problems can be applied not only to compute the best therapy for each patient, but also to determine the specific type of resistance suffered by the patient, and to quantify the degree of resistance of each patient.

    The existence of allergies or abnormal toxicity can also be introduced to design personalized therapies. In this case, given that the patient experiences higher deleterious effects associated to imatinib administration, the change affects the objective function measuring malignancy. In this respect, as with the previous personalized problems, the optimal therapy model that arises is useful for several purposes. First, it provides a guide for physicians to determine the extent to which a decrease in the drug dose must be recommended for each patient. Second, it allows the patient's specific toxicity degree to be quantitatively measured. And third, it makes the computation of the personalized optimal therapy possible.

    Our baseline optimal problem can also contemplate the presence of adverse factors different from BCR-ABL1 mutations. As happened for the other personalized optimal therapies, after properly modifying the system of equations describing the dynamics of the treated disease, the corresponding optimal therapy problem is applied to quantify the effects of these factors for each patient and to compute the personalized optimal therapy.

    Summing up, we have shown the great potential and versatility that optimal control theory has in formulating personalized therapies, a question still not sufficiently explored and with multiple clinical and medical applications of interest. On this point, it is worth noting that the success of optimal therapy models in allowing patient specific therapies and medical analyses is the result of not only purely mathematical issues, but also of the connections and complementarity between all the applied techniques. Indeed, as we have illustrated along this paper, the clinical and medical applicability -the ultimate objective of any biomathematical research- of our proposed personalized models relies on the joint and proper use of calibration, simulation, and mathematical approaches and techniques. This complementarity allows the personalized optimal therapy to be computed for all the relevant cases, a possibility that does not exist in other models (see [56] and [77]).

    Finally, it is necessary to mention that all the theoretical results generated by our personalized optimal therapy models are corroborated by clinical evidence on CML patients susceptible to treatment with imatinib. Physicians elaborate and decide treatments on the basis of heuristic procedures and continuous medical monitoring. This trial and error method and the successive corrections originate the observed clinical treatments, which must therefore be envisaged as an approximation of the exact personalized optimal therapy arising from our optimal control problems. As shown along the paper, this is in fact just what happens: qualitatively, all the features observed for the clinical treatments are present in our personalized optimal therapies, and in addition, are also quantitatively reproduced on average. Concerning the biomedical and clinical interpretation and applicability of our optimal therapy model for imatinib treated CML, an important last comment is necessary. As shown along the paper, the mathematical model satisfactorily describes CML susceptible to imatinib treatment under different situations of the patient. Indeed, this is the main purpose of the research. In line with our reasonings on the specificity of the model to its objective, it must not be applied to explain other situations. One of these conditions is the already mentioned presence of BCR-ABL1 mutations, but there are others of great clinical interest that are also excluded, related to the evolution of the disease after the discontinuation of the imatinib treatment. Obviously, our model is able to describe the behavior of the disease once imatinib treatment is stopped, but only under the assumption of the existence of four homogeneous types of cells, namely cancer and normal HSC and DC. However, as shown by [78] and [79], the safe discontinuation of treatment or the relapse of CML are determined by the presence of a subpopulation of cancer HSC at a quiescent stage, and also by the characteristics of the natural killer cells in the patient's immune system. Accordingly, any mathematical description of the evolution of CML after the discontinuation of imatinib treatment should include the dynamics of these two additional cell populations, specifying not only how they are produced and disappear, but also their relationships with the other relevant cell populations. Without any doubt, this will constitute an important objective of the future research agenda on optimal therapies in CML, jointly with the consideration of the effects and consequences of the presence of BCR-ABL1 mutations. In this respect, to the extent that medical and biological researchers identify, clarify and quantify the mechanisms of BCR-ABL1 actions and their relationships with the relevant cell populations, and, on the other hand, those concerning the quiescent cancer HSC and the killer cells responsible for CML relapse, biomathematicians will progressively incorporate the obtained knowledge to the existing optimal therapy models.

    PJ Gutiérrez gratefully acknowledges financial support from projects MTM2017-85476-C2-1-P (Spanish Ministry of Economy and Competitiveness and European FEDER Funds) and VA148G18 (Castile and Leon Autonomous Government and European FEDER Funds). We thank Alan F. Hynds, B.A. Dip. TEFL, for manuscript preparation. The suggestions and comments of an anonymous reviewer noticeably improved the final outcome and are highly appreciated.

    The authors declare no conflicts of interest in this paper.



    [1] A. J. Lotka, Contribution to the Theory of Periodic Reaction, J. Phys. Chem., 14 (1910), 271-274.
    [2] A. J. Lotka, Analytical Note on Certain Rhythmic Relations in Organic Systems, Proc. Natl. Acad. Sci. U.S.A., 6 (1920), 410-415.
    [3] A. J. Lotka, Elements of Physical Biology, Williams and Wilkins, 1925.
    [4] V. Volterra, Variazioni e fluttuazioni del numero d'individui in specie animali conviventi, Mem. R. Accad. Naz. Lincei, 2 (1926), 31-113.
    [5] V. Volterra, Fluctuations in the abundance of a species considered mathematically, Nature, 118, (1926), 558-560.
    [6] V. Volterra, Leçons sur la théorie mathématique de la lutte pour la vie, Paris: Gauthier-Villars, 1931.
    [7] L. A. Loeb, K. R. Loeb, J. P. Anderson, Multiple mutations and cancer, Proc. Natl. Acad. Sci. U.S.A., 100 (2003), 776-781.
    [8] G. Magombedze, W. Garira, E. Mwenje, C. P. Bhunu, Optimal control for HIV-1 multi-drug therapy, Int. J. Comp. Math., 88 (2011), 314-340.
    [9] N. Tarfulea, A mathematical model for HIV treatment with time-varying antiretroviral therapy, Int. J. Comp. Math., 88, (2011), 3217-3235.
    [10] Y. Koizumi, S. Iwami, Mathematical modeling of multi-drugs therapy: A challenge for determining the optimal combinations of antiviral drugs, Theor. Biol. Med. Mod., 11 (2014), 41.
    [11] H. T. Banks, Modeling and Control in the Biomedical Sciences, Springer, 1975.
    [12] J. D. Murray, Mathematical Biology I: An Introduction, Springer, 2002.
    [13] J. D. Murray, Mathematical Biology II: Spatial Models and Biomedical Applications, Springer, 2003.
    [14] C. W. Clark, Mathematical Bioeconomics. The Optimal Management of Renewable Resources, John Wiley & Sons, Inc., 1990.
    [15] U. Ledzewicz, H. Schättler, A. Friedman, E. Kashdan, Mathematical methods and models in Biomedicine, Springer, 2013.
    [16] M. Eisen, Mathematical Methods and Models in the Biological Sciences, Prentice Hall, 1988.
    [17] P. J. Gutiérrez Diez, I. H. Russo, J. Russo, The Evolution of the Use of Mathematics in Cancer Research, Springer, 2012.
    [18] A. Świerniak, U. Ledzewicz, H. Schättler, Optimal Control for a Class of Compartmental Models in Cancer Chemotherapy, Int. J. Appl. Math. Comp. Sci., 13 (2003), 357-368.
    [19] J. M. Murray, Optimal control for a cancer chemotherapy problem with general growth and cost functions, Math. Biosci. 98 (1990), 273-287.
    [20] J. M. Murray, Some optimal control problems in cancer chemotherapy with a toxicity limit, Math. Biosci., 100 (1990), 49-67.
    [21] S. Nanda, H. Moore, S. Lenhart, Optimal control of treatment in a mathematical model of chronic myelogenous leukemia, Math. Biosci., 210 (2007), 143-156.
    [22] B. E. Aïnseba, C. Benosman, Optimal control for resistance and suboptimal response in CML, Math. Biosci., 227 (2010), 81-93.
    [23] K. Schepers, E. M. Pietras, D. Reynaud, J. Flach, M. Binnewies, T. Garg, et al., Myeloproliferative neoplasia remodels the endosteal bone marrow niche into a self-reinforcing leukemic niche, Cell Stem Cell, 13 (2013), 285-299. doi: 10.1016/j.stem.2013.06.009
    [24] S. Kapoor, V. P. Rallabandi, C. Sakode, R. Padhi, P.K. Roy, A patient-specific therapeutic approach for tumour cell population extinction and drug toxicity reduction using control systems-based dose-profile design, Theor. Biol. Med. Mod., 10 (2013), 68.
    [25] P. J. Gutiérrez-Diez, M. Á.,López-Marcos, J. Martínez-Rodríuez, J. Russo, The effects of time valuation in cancer optimal therapies: A study of chronic myeloid leukemia, Theor. Biol. Med. Mod., 16 (2019), 10.
    [26] M. Baccarani, F. Castagnetti, G. Gugliotta, G. Rosti, A review of the European LeukemiaNet recommendations for the management of CML, Ann. Hematol., 94 (2015), S141-S147.
    [27] M. W. Deininger, J. G. Hodgson, N. P. Shah, J. E. Cortes, D. W. Kim, F. E. Nicolini, et al., Compound mutations in BCR-ABL1 are not major drivers of primary or secondary resistance to ponatinib in CP-CML patients, Blood, 127 (2016), 703-712.
    [28] F. Michor, T. P. Hughes, Y. Iwasa, S. Branford, N. P. Shah, C. L. Sawyers, et al., Dynamics of chronic myeloid leukaemia, Nature, 435 (2005), 1267-1270.
    [29] D. Murray Lyon, Does the reaction to adrenalin obey Weber's Law?, J. Pharmacol. Exp. Ther., 21 (1923), 229-235.
    [30] M. Inoue, K. Kaneko, Weber's law for biological responses in autocatalytic networks of chemical reactions, Phys. Rev. Lett., 107 (2011), 048301.
    [31] H. W. Sinn, Weber's law and the biological evolution of risk preferences: The selective dominance of the logarithmic utility function, Gen. Pap. Risk Ins. Theor., 28 (2002), 87-100.
    [32] R. Capdeville, S. Silverman, Imatinib: A targeted clinical drug development, Semin. Hematol., 40 (2003), 15-20.
    [33] M. Bonifacio, F. Stagno, L. Scaffidi, M. Krampera, F. Di Raimondo, Management of Chronic Myeloid Leukemia in Advanced Phase, Front. Oncol., 9 (2019), 1132.
    [34] S. Soverini, A. Hochhaus, F. E. Nicolini, F. Gruber, T. Lange, G. Saglio, et al., BCR-ABL kinase domain mutation analysis in chronic myeloid leukemia patients treated with tyrosine kinase inhibitors: Recommendations from an expert panel on behalf of European LeukemiaNet. Blood, 118 (2011), 1208-1215.
    [35] I. Galinsky, S. Buchanan, Guide to Interpreting Disease Responses in Chronic Myeloid Leukemia, J. Adv. Prac. Oncol., 3, (2012), 225-236.
    [36] S. Prabhu, D. Saadat, M. Zhang, L. Halbur, J. P. Fruehauf, S. T. Ong, A novel mechanism for BCR-ABL action: BCR-ABL-mediated induction of the eIF4F translation initiation complex and mRNA translation, Oncogene, 26 (2007), 1188-1200.
    [37] L. Hu, L. Pu, D. Yang, C. Zhang, H. Wang, Y. Ding, et al., How to detect the rare BCR-ABL (e14a3) transcript: A case report and literature review, Oncol. Lett., 14 (2017), 5619-5623.
    [38] S. I. Ismail, R. G. Naffa, A. F. Yousef, M. T. Ghanim, Incidence of BCR-ABL fusion transcripts in healthy individuals, Mol. Med. Rep. 9 (2014), 1271-1276.
    [39] T. Hughes, G. Saglio, A. Quintás-Cardama, M. J. Mauro, D. W. Kim, J. H. Lipton, et al., BCR-ABL1 mutation development during first-line treatment with dasatinib or imatinib for chronic myeloid leukemia in chronic phase, Leukemia, 29 (2015), 1832-1838. doi: 10.1038/leu.2015.168
    [40] J. Kaeda, D. O'Shea, R. M. Szydlo, E. Olavarria, F. Dazzi, D. Marin, et al., Serial measurement of BCR-ABL transcripts in the peripheral blood after allogeneic stem cell transplantation for chronic myeloid leukemia: An attempt to define patients who may not require further therapy, Blood, 107 (2006), 4171-4176.
    [41] M. Houshmand, G. Simonetti, P. Circosta, V. Gaidano, A. Cignetti, G. Martinelli, et al., Chronic myeloid leukemia stem cells, Leukemia, 33 (2019), 1543-1556.
    [42] F. Loscocco, G. Visani, S. Galimberti, A. Curti, A. Isidori, BCR-ABL Independent Mechanisms of Resistance in Chronic Myeloid Leukemia, Front. Oncol., 9 (2019), 939.
    [43] H. Kitamura, Y. Tabe, T. Ai, K. Tsuchiya, M. Yuri, S. Misawa, et al., A new highly sensitive real-time quantitative-PCR method for detection of BCR-ABL1 to monitor minimal residual disease in chronic myeloid leukemia after discontinuation of imatinib, PLoS One, 14 (2019), e0207170.
    [44] R. Arora, R. D. Press, Measurement of BCR-ABL1 transcripts on the International Scale in the United States: Current status and best practices, Leuk. Lymphoma, 58 (2017), 8-16.
    [45] Physycians' Desk Reference, PDR Network (ed) 64 edition, LLC at Montvale, 2010.
    [46] A. Quintás-Cardama, J. Cortes, H. Kantarjian, Biology of Chronic and Acute Myeloid Leukemia, in The Molecular Basis of Cancer, Elsevier, (2008), 371-383.
    [47] L. Löf, L. Arngården, U. Olsson-Strömberg, B. Siart, M. Jansson, J. S. Dahlin, et al., Flow Cytometric Measurement of Blood Cells with BCR-ABL1 Fusion Protein in Chronic Myeloid Leukemia, Sci. Rep. 7 (2017), 623.
    [48] D. Raspadori, P. Pacelli, A. Sicuranza, E. Abruzzese, A. Iurlo, D. Cattaneo, et al., Flow Cytometry Assessment of CD26+ Leukemic Stem Cells in Peripheral Blood: A Simple and Rapid New Diagnostic Tool for Chronic Myeloid Leukemia, Cytometry Part B, 96 (2019), 294-299.
    [49] F. E. Craig, K. A. Foon, Flow cytometric immunophenotyping for hematologic neoplasms, Blood, 111 (2008), 3941-3967.
    [50] D. Campana, Minimal residual disease in acute lymphoblastic leukemia, Semin. Hematol., 46 (2009), 100-106.
    [51] D. Campana, Role of minimal residual disease monitoring in adult and pediatric acute lymphoblastic leukemia, Haematol./Oncol. Clin. N.A., 23 (2009), 1083-1098.
    [52] N. Shah,Medical management of CML, Hematology, 1 (2007), 371-375.
    [53] P. K. Bhamidipati, H. Kantarjian, J. Cortes, A. M. Cornelison, E. Jabbour, Management of imatinib-resistant patients with chronic myeloid leukemia, Ther. Adv. Hematol., 4 (2013), 103-117.
    [54] M. S. Marcolino, E. Boersma, N. C. D. Clementino, A. V. Macedo, A. D. Marx-Neto, M. H. C. R. Silva, et al., Imatinib treatment duration is related to decreased estimated glomerular filtration rate in chronic myeloid leukemia patients, Ann. Oncol., 22 (2011), 2073-2079.
    [55] M. Baccarani, J. Cortes, F. Pane, D. Niederwieser, G. Saglio, J. Apperley, et al., Chronic myeloid leukemia: An update of concepts and management recommendations of European LeukemiaNet, J. Clin. Oncol., 27 (2009), 6041-6051.
    [56] R. A. Everett, Y. Zhao, K. B. Flores, Y. Kuang, Data and implication based comparison of two chronic myeloid leukemia models, Math. Biosci. Eng., 10 (2013), 1501-1518.
    [57] C. Fava, G. Saglio, Can we and should we improve on frontline imatinib therapy for chronic myeloid leukemia?, Semin. Hematol., 47 (2010), 319-326.
    [58] P. Valent, Imatinib-resistant chronic myeloid leukemia (CML): Current concepts on pathogenesis and new emerging pharmacologic approaches, Biol. Targets Ther., 1 (2007), 433-448.
    [59] R. Bhatia, Chronic Myeloid Leukemia, in Hematology, seventh edition, Elsevier, (2018), 1055-1070.
    [60] T. R. Randolph, Myeloproliferative neoplasms, in Rodak's Hematology Clinical Applications and Principles, St. Louis, Missouri: Saunders, (2015), 561-590.
    [61] N. Bouizem, M. Helal, B. Aïnseba, A. Lakmeche, Leukemia mathematical model, ITM Web Conf., 4 (2015), 01006.
    [62] A. Hochhaus, S. Kreil, A. S. Corbin, P. La Rosée, M. C. Müller, T. Lahaye, et al., Molecular and chromosomal mechanisms of resistance to imatinib (STI571) therapy, Leukemia, 16 (2002), 2190-2196.
    [63] V. Karavasilis, A. Reid, R. Sinha, J. S. De Bono, Cancer drug resistance, in Cancer Drug Design and Discovery, (S. Neidle ed.), (2008), Elsevier.
    [64] M. Baccarani, G. Saglio, J. Goldman, A. Hochhaus, B. Simonsson, F. Appelbaum, et al., Evolving concepts in the management of chronic myeloid leukemia: recommendations from an expert panel on behalf of the European Leukemia Net, Blood, 108 (2006), 1809-1820.
    [65] R.P. Nelson Jr, K. Cornetta, K. E. Ward, S. Ramanuja, C. Fausel, L.D. Cripe, Desensitization to imatinib in patients with leukemia, Ann. Allergy, Asthma, Immunol., 97 (2006), 216-222.
    [66] M. Albayrak, H. Celebi, A. Albayrak, E. S. Can, V. Aslan, B. Onec, et al., Serious skin reaction associated with imatinib in a patient with chronic myeloid leukemia, Eurasian J. Med., 43 (2011), 192-195.
    [67] V. Chou, S. McClelland, D. Resnick, M. Lee-Wong, Successful Desensitization of an Adult with Type I Hypersensitivity to Imatinib, Internet J. Asthma, Allergy Immunol., 4 (2004), 2.
    [68] E. Faber, M. Divoká, I. Skoumalová, M. Novák, I. Marešová, K. Mičová, et al., A lower dosage of imatinib is sufficient to maintain undetectable disease in patients with chronic myeloid leukemia with long-term low-grade toxicity of the treatment, Leuk. Lymphoma, 57 (2016), 370-375.
    [69] Y. Zhu, S. X. Qian, Clinical efficacy and safety of imatinib in the management of Ph(+) chronic myeloid or acute lymphoblastic leukemia in Chinese patients, OncoTargets Ther., 7 (2014), 395-404.
    [70] European Leukemia Net, 2020. Available from: https://www.leukemia-net.org/content/leukemias/cml/euro_and_sokal_score/index_eng.html.
    [71] L. C. Kuntegowdanahalli, G. B. Kanakasetty, A. H. Thanky, L. Dasappa, L. A. Jacob, S. B. Mallekavu, et al., Prognostic and predictive implications of Sokal, Euro and EUTOS scores in chronic myeloid leukaemia in the imatinib era-experience from a tertiary oncology centre in Southern India, Ecancermedicalscience, 10 (2016), 679.
    [72] S. Chhikara, S. Sazawal, K. Singh, R. Chaubey, H. Pati, S. Tyagi, et al., Comparative analysis of the Sokal, Euro and European Treatment and Outcome Study score in prognostication of Indian chronic myeloid leukemia-chronic phase patients on imatinib, S.A. J. Cancer, 7, (2018), 258-262.
    [73] J. Hasford, M. Baccarani, V. Hoffmann, J. Guilhot, S. Saussele, G. Rosti, et al., Predicting complete cytogenetic response and subsequent progression-free survival in 2060 patients with CML on imatinib treatment: the EUTOS score, Blood, 118 (2011), 686-692. doi: 10.1182/blood-2010-12-319038
    [74] M. Baccarani, F. Castagnetti, G. Gugliotta, G. Rosti, A review of the European Leukemia Net recommendations for the management of CML, Ann. Hematol., 94 (2015), S141-S147.
    [75] M. Schemionek, T. Spieker, L. Kerstiens, C. Elling, M. Essers, A. Trumpp, et al., Leukemic spleen cells are more potent than bone marrow-derived cells in a transgenic mouse model of CML, Leukemia, 26, (2012), 1030-1037.
    [76] L.O. Ngolet, I. Kocko, A. Elira-Dokekias, Pregnancy and Accelerated Phase of Myeloid Chronic Leukemia Treated with Imatinib: A Case Report from a Developing Country, Case Rep. Hematol., (2016), 6104948.
    [77] E. Fimmel, Y. S. Semenov, A. S. Bratus, On optimal and suboptimal treatment strategies for a mathematical model of leukemia, Math. Biosci. Eng., 10 (2013), 151-165.
    [78] M. Bocchia, A. Sicuranza, E. Abruzzese, A. Iurlo, S. Sirianni, A. Gozzini, et al., Residual Peripheral Blood CD26+ Leukemic Stem Cells in Chronic Myeloid Leukemia Patients During TKI Therapy and During Treatment-Free Remission, Front. Oncol., 8 (2018), 194
    [79] G. Caocci, B. Martino, M. Greco, E. Abruzzese, M. M. Trawinska, S. Lai, et al., Killer immunoglobulin-like receptors can predict TKI treatment-free remission in chronic myeloid leukemia patients, Exp. Hematol., 43 (2015), 1015-1018.
  • This article has been cited by:

    1. Samara Sharpe, Hana M. Dobrovolny, Predicting the effectiveness of chemotherapy using stochastic ODE models of tumor growth, 2021, 101, 10075704, 105883, 10.1016/j.cnsns.2021.105883
    2. Hannah G. Anderson, Gregory P. Takacs, Jeffrey K. Harrison, Libin Rong, Tracy L. Stepien, Optimal control of combination immunotherapy for a virtual murine cohort in a glioblastoma-immune dynamics model, 2024, 595, 00225193, 111951, 10.1016/j.jtbi.2024.111951
  • Reader Comments
  • © 2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(4264) PDF downloads(202) Cited by(2)

Figures and Tables

Figures(18)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog