Over the past decade, several targeted therapies (e. g. imatinib, dasatinib, nilotinib) have been developed to treat Chronic Myeloid Leukemia (CML). Despite an initial response to therapy, drug resistance remains a problem for some CML patients. Recent studies have shown that resistance mutations that preexist treatment can be detected in a substantial number of patients, and that this may be associated with eventual treatment failure. One proposed method to extend treatment efficacy is to use a combination of multiple targeted therapies. However, the design of such combination therapies (timing, sequence, etc.) remains an open challenge. In this work we mathematically model the dynamics of CML response to combination therapy and analyze the impact of combination treatment schedules on treatment efficacy in patients with preexisting resistance. We then propose an optimization problem to find the best schedule of multiple therapies based on the evolution of CML according to our ordinary differential equation model. This resulting optimization problem is nontrivial due to the presence of ordinary different equation constraints and integer variables. Our model also incorporates drug toxicity constraints by tracking the dynamics of patient neutrophil counts in response to therapy. We determine optimal combination strategies that maximize time until treatment failure on hypothetical patients, using parameters estimated from clinical data in the literature.
Chronic Myeloid Leukemia (CML) is an acquired hematopoietic stem cell disorder leading to the over-proliferation of myeloid cells and an increase in cellular output from the bone marrow that is often associated with splenomegaly. The most common driving mutation in CML is a translocation between chromosomes 9 and 22 that produces a fusion gene known as BCR-ABL. The BCR-ABL protein promotes proliferation and inhibits apoptosis of myeloid progenitor cells and thereby drives expansion of this cell population. By targeting the BCR-ABL oncoprotein, imatinib (brand name Gleevec) is able to induce a complete cytogenetic remission in the majority of chronic phase CML patients. A minority of patients, however, either fail to respond or eventually develop resistance to treatment with imatinib [1]. It is thought that a primary driver of this resistance to imatinib is point mutations within the BCR-ABL gene. A recent study utilizing sensitive detection methods demonstrated that a small subset of these mutations may exist before the initiation of therapy in a significant fraction of patients, and that this status is correlated with eventual treatment failure [2]. Second generation agents such as dasatinib and nilotinib have been developed and each has shown efficacy against various common mutant forms of BCR-ABL. This leads to the observation that the various mutant forms of BCR-ABL result in CML that have unique dynamics under therapy, and that combinations of these inhibitors may be necessary to effectively control a rapidly evolving CML population. Patients with CML often die due to transformation of the disease into an acute form of leukemia known as blast crisis. It has been shown that blast crisis is due to the accumulation of additional mutations in CML progenitor cells [3]. The goal of this work is to leverage the differential responses of CML mutant strains to design novel sequential combination treatment schedules using dasatinib, imatinib and nilotinib that optimally control leukemic burden and delay treatment failure due to preexisting resistance. We develop and parametrize a mathematical model for the evolution of both wild-type (WT) CML and mutated (resistant) CML cells in the presence of each therapy. Then we formulate the problem as a discrete optimization problem in which a sequence of monthly treatment decisions is optimized to identify the temporal sequence of imatinib, dasatinib and nilotinib administration that minimizes the total CML cell population over a long time horizon. There has been a significant amount of work done in the past to mathematically model CML. For example, in [4] the authors developed a system of ordinary differential equations (ODEs) that model both the normal progression from stem cell to mature blood cells and abnormal progression of CML. A hierarchical system of differential equations was used to model the response of CML cells to imatinib therapy in [5]; this model fit the biphasic nature of decline in BCR-ABL positive cells during imatinib treatment. In [6] the authors investigated the number of different resistant strains present in a newly diagnosed chronic phase CML patient. An optimal control approach was utilized to optimize imatinib scheduling in [7]. Particularly relevant to our work is [8,9] where the authors investigated simultaneous continuous administration of dasatinib, nilotinib and imatinib; in particular, they explored the minimal number of drugs necessary to prevent drug resistance. In the current work, we focus on understanding the optimal administration schedule of multiple therapies to prevent resistance, and studying the impact of toxicity constraints on optimal scheduling. Since several of the available tyrosine kinase inhibitors (TKI) share similar toxicities (in particular neutropenia, see e. g. , [10–12]) combining them together can lead to elevated risk of adverse events. Thus we consider sequential combination therapies in which only one TKI may be administered at a time. Moreover, it has been shown that the risk of treatment failure and blast crisis are highest within the first 2 years from diagnosis [1]. Therefore it is possible that optimized, sequential single agent therapy may be sufficient to minimize the risk of treatment failure. Allowing only one treatment at a time leads to a complex, time-dependent discrete optimization problem. Another line of research closely related to the current work is the use of optimal control techniques in the design of optimal temporally continuous drug concentration profiles (see, e. g. , review articles [13,14] and the textbook [15]). In this field the tools of optimal control such as the Pontryagin principle and the Euler-Lagrange equations are used to find drug concentration profiles that result in minimal tumor cell populations under toxicity constraints. Particularly relevant to the current work is [16] where the authors searched for optimal anti-HIV treatment strategies. They dealt with the similar problem of treating heterogeneous populations with multiple drugs. One major drawback of these works is the fact that it is nearly impossible to achieve a specific optimal continuous drug-concentration profile in patients, since drug concentration over time is a combined result of a treatment schedule (e. g. sequence of discrete oral administrations) and pharmacokinetic processes in the body including metabolism, elimination, etc. Thus the clinical utility of an optimal continuous drug concentration profile is limited. In contrast to these previous works, here we model the optimization problem as a more clinically realistic sequence of monthly treatment decisions. Imposition of this fixed discrete set of decision times leads to a challenging optimization problem. Such dynamical systems are referred to as ‘switched nonlinear systems’ in the control community [17], and our problem additionally imposes fixed switching times. In this work we will leverage the system structure and tools from mixed-integer linear optimization [18] to solve this problem numerically, resulting in optimal therapy schedules that are easy to implement in practice.
We consider an ODE model of the differentiation hierarchy of hematopoietic cells, adapted from [5,19,20]. Stem cells (SC) on top of the hierarchy give rise to progenitor cells (PC), which produce differentiated cells (DC), which in turn produce terminally differentiated cells (TC). This differentiation hierarchy applies to both normal and leukemic cells [21]. We consider in our model leukemic WT cells as well as preexisting BCR-ABL mutant cell types. We use type 1, type 2, and type i (3 ≤ i ≤ n) cells to denote normal, leukemic WT, and (n − 2) leukemic mutant cells; layer 1,2, 3,4 cells to denote SC, PC, DC, and TC; and drug 0,1, 2,3 to denote a drug holiday, nilotinib, dasatinib, and imatinib, respectively. Let xl, i (t) denote the abundance of type i cell at layer l and time t, and x (t) = (xl, i (t) ) be the vector of all cell abundance at time t. If drug j ∈ {0,1, 2,3} is taken from month m to month m + 1, then the cell dynamics are modeled by the following set of ODEs. x ˙ (t) = f j (x (t) ), t ∈ [ m Δ t, (m + 1) Δ t ], x (m Δ t) = x m, for some function fj, where Δt = 30 days and xm is the cell abundance at the beginning of month m. The concrete form of function fj under drug j is described as follows. SC level x ˙ 1, i = (b 1, i j ϕ i - d 1, i j) x 1, i, i = 1, …, n PC level x ˙ 2, i = b 2, i j x 1, i - d 2, i j x 2, i, i = 1, …, n DC level x ˙ 3, i = b 3, i j x 2, i - d 3, i j x 3, i, i = 1, …, n TC level x ˙ 4, i = b 4, i j x 3, i - d 4, i j x 4, i, i = 1, …, n. (1) See Fig 1 for an illustration of the differentiation hierarchy of hematopoietic cells, including neutrophils as part of the TC. Here we describe the function of each parameter of this model. For a detailed discussion of how these parameters were estimated from biological data, please see Section A of S1 Text. Type i stem cells divide at rate b 1, i j per day under drug j. The production rates of type i progenitors, differentiated cells and terminally differentiated cells under drug j are b l, i j per day for l = 2,3, 4, respectively. The type i cell at layer l dies at rate d l, i j per day under drug j, for each i, l and j. The competition among normal and leukemic stem cells is modeled by the density dependence functions ϕi (t), where ϕ i (t) = 1 / (1 + p i ∑ i = 1 n x 1, i (t) ) for each i; these functions ensure that the normal and leukemic stem cell abundances remain the same once the system reaches a steady state. The parameter p1 (resp. p2) is computed from the equilibrium abundance of normal (resp. leukemic WT) stem cells assuming only normal (resp. leukemic WT) cells are present, and we set pi = p2 for each i ≥ 3. In particular, p 1 = (b 1,1 0 / d 1,1 0 - 1) / K 1 and p 2 = (b 1,2 0 / d 1,2 0 - 1) / K 2, with K1 (resp. K2) being the equilibrium abundance of normal (resp. leukemic WT) stem cells assuming only normal (resp. leukemic WT) cells are present.
Note the total leukemic cell abundance at day t is given by ∑l≥1∑i≥2 xl, i (t). The OTP can be formulated as the following mixed-integer optimization problem with ODE constraints. min ∑ l ≥ 1 ∑ i ≥ 2 x l, i (M Δ t) (2) s. t. x ˙ (t) = ∑ j = 0 3 z m, j f j (x (t) ), t ∈ [ m Δ t, (m + 1) Δ t ], m = 0,1, …, M - 1, (3) ∑ j = 0 3 z m, j = 1, m = 0,1, …, M - 1, (4) z m, j ∈ { 0,1 }, j = 0,1, 2,3, m = 0,1, …, M - 1, (5) x (0) = x 0. (6) To summarize the previous display, in eq (2) we state that our objective is to minimize the leukemic cell population at the end of the treatment horizon. In eq (3) we stipulate that the cell dynamics are governed by the system of differential equations given by eq (1). Together eqs (4) and (5) stipulate that during each time period we administer either one drug or no drug. The OTP problem is a mixed-integer nonlinear optimization problem, in which some constraints are specified by the solution to a nonlinear system of ODEs eq (3). This optimization problem is beyond the ability of state-of-the-art optimization software. However, if we assume the TKI therapies do not affect the stem cell compartment, then it is possible to handle the ODE constraints numerically. This is because the non-linearities in the ODE model are only present in the stem cell compartment, and the remaining compartments are modeled by linear differential equations. Thus we are able to build a refined linear approximation to the ODE constraints (see Section C of S1 Text), and recast the problem as a mixed-integer linear optimization problem, which can be solved efficiently for the size of our problem (see Section B of S1 Text). It should be noted that our model does not explicitly consider the phenomena of TKI resistance acquired during therapy. Our model suggests a method for designing optimized treatment plans, beyond monotherapies, at the beginning of treatment. In addition, this optimization procedure can be re-run and modified during the course of treatment, with updated inputs from each patient’s response to the treatment—including the presence of acquired mutations. Below we summarize our notation for the ease of the reader.
We first utilize the model to demonstrate the dynamics of CML populations with preexisting BCR-ABL mutations under monotherapy with the standard therapies imatinib, dasatinib and nilotinib. Recall that the standard dosing regimens are 300mg twice daily for nilotinib, 100mg once daily for dasatinib, and 400mg once daily for imatinib [22]. The birth rate parameters for each leukemic mutant type under different TKIs (b l, i j for j = 1,2, 3, l = 1,2, 3,4, and i ≥ 3) in the model are estimated using in vitro IC50 values reported in [23] for each drug. The initial cell populations at the start of therapy are derived by running the model starting from clonal expansion of a single leukemic cell in a healthy hematopoietic system at equilibrium [19] until CML detection (when the total leukemic burden reaches approximately 1012 cells [24]). At this point the total cell burden is 2–3 times the normal cell burden in a healthy individual and thus the total leukemic cells make up approximately 77% of the total cell population; this is consistent with clinical reports [25]. Details on deriving the initial cell abundances at diagnosis are provided in Section A of S1 Text. In the first example we consider a patient harboring a low level of the BCR-ABL mutant F317L before the initiation of TKI therapy. According to the in vitro IC50 value reported in [23], F317 is resistant to dasatinib, and moderately resistant to nilotinib and imatinib. The initial population conditions are given in Table 1 with the leukemic WT and F317L cells taking up 95% and 5% of the leukemic cells, respectively. We plot in Fig 2 the cell dynamics over 120 months for four treatment plans: (1) nilotinib monotherapy (2) dasatinib monotherapy, (3) imatinib monotherapy, (4) no therapy—control. We observe that as predicted, the disease burden responds well to imatinib and nilotinib; the percentage of cancerous cells after a 24 month treatment drops to 0. 19% with nilotinib and 0. 26% with imatinib, respectively. However, the F317L mutant population is fairly resistant to dasatinib; we observe that the percentage of cancerous cells after 24 months is 58. 1% with dasatinib and 95. 4% with no treatment. Over the 120 month period dasatinib treatment provides only modest improvement over the ‘no drug’ option in controlling the F317L population; however, dasatinib remains quite effective in controlling the WT leukemic population. It is interesting to note that overall, nilotinib is the most effective in controlling both the WT and F317L leukemic populations. However, nilotinib also negatively impacts the healthy cell population more severely than imatinib, which is slightly less effective in controlling the leukemic populations. This suggests that some trade-offs between these drugs exist, and these trade-offs may be exploited in designing combination therapies. In the next example we consider a patient with BCR-ABL mutant type M351T preexisting therapy. In contrast to the previous example, this commonly occurring mutant has been found to be partially sensitive in varying degrees to all three therapies. The initial conditions are given in Table 2. Once again we have assumed that WT and M351T cells take up 95% and 5% of total leukemic cells, respectively. In Fig 3 the cell dynamics over 120 months for the four standard treatment plans are plotted: (1) nilotinib monotherapy (2) dasatinib monotherapy, (3) imatinib monotherapy, (4) no therapy—control. Since the M351T mutant is responsive to each drug in contrast to the previous example, the percentage of cancerous cells after a 24 month treatment drops to 0. 18% with nilotinib, 0. 18% with dasatinib, and 0. 25% with imatinib, respectively. Without treatment, the percentage of cancerous cells after 24 months is 95. 4%. Here, we observe that although nilotinib is more effective than dasatinib in controlling the total mutant M351T burden, the effect is reversed in the progenitor population. Higher levels of stem and progenitor populations will lead to faster rebound during treatment breaks, suggesting another trade-off to consider in the combination setting. We next solve the discrete optimization problem to identify sequential combination therapies utilizing imatinib, dasatinib and nilotinib to optimally treat CML patients with preexisting BCR-ABL mutations. We consider schedules in which a monthly treatment decision is made between one of four choices: imatinib, dasatinib, nilotinib, and drug holiday. During months in which one of the three drugs is administered, the dosing regimen is fixed at 300mg twice daily for nilotinib, 100mg once daily for dasatinib, and 400mg once daily for imatinib. In the following we optimize over feasible treatment decision sequences that result in a minimal leukemic cell burden after 3 years. Each treatment plan is completely characterized by a temporal sequence of drugs over a long time horizon. The extent to which TKIs affect leukemic stem cells is currently unknown. Several studies have demonstrated that these cells may in fact be resistant to TKIs, see e. g. [27] and references therein. Based on these results we have assumed throughout this work that TKIs do not affect the leukemic stem cell population. Given that there is uncertainty regarding this issue, in this subsection we investigate how our optimal schedules will perform if we assume TKIs impact the leukemic stem cells. In Tables 9 and 10, we look at the leukemic cell burden after a 36-month treatment with our optimal therapy and the best monotherapy for preexisting M351T and F317L, respectively, assuming that each TKI reduces the production rate of leukemic stem cells by a certain level. In both tables we observe that the optimized schedule outperforms the monotherapy even if leukemic stem cells are susceptible to TKI. These result imply that the optimal treatment schedules we derive here outperform a traditional monotherapy regardless of the impact of TKIs on leukemic stem cells.
In this work we have considered the problem of finding optimal treatment schedules for the administration of a variety of TKIs for treating chronic phase CML. We modeled the evolution of wild-type and mutant leukemic cell populations with a system of ordinary differential equations. We then formulated an optimization problem to find the sequence of TKIs that lead to a minimal cancerous cell population at the end of a fixed time horizon of 36 months. The 36-month therapeutic horizon is clinically meaningful since it appears that the risk of therapeutic failure and disease progression to blast crisis is highest within the first two years from diagnosis [1]. At first glance the optimization problem studied in this work (OTP) is quite challenging. It is a mixed-integer nonlinear optimization problem, where the nonlinear constraints are specified by the solution to a nonlinear system of differential equations. However, one factor mitigating the complexity of the problem is the assumption that the TKIs do not effect the stem cell compartment. This has the effect of making the evolution of the stem cell compartment independent of the TKI schedule chosen. In addition, the remaining layers in the cellular hierarchy are modeled by linear differential equations. We can thus numerically solve the differential equation governing the stem cell layer, and treat this function as an inhomogeneous forcing term in the linear differential equation governing the progenitor cells. This allows us to approximate the nonlinear constraints specified by the differential equations by linear constraints with high accuracy. Then the OTP problem can be approximated by a mixed-integer linear optimization problem, which we are able to solve with state-of-the-art optimization software CPLEX [28] within one hour. Importance of minimizing progenitor cell population. We first aimed to minimize leukemic cell burden at 36 months after initiation of therapy, starting with an initial leukemic population of wild-type CML cells and either M351T (sensitive to all three therapies) or F317L (resistant to dasatinib) mutant leukemic cells. For both starting mutant populations, we observed that the optimal schedule involves initiating therapy with dasatinib and later switching to nilotinib, although the timing of the switch differed. To further understand this result, we noted that within this parameter regime, dasatinib is the most effective of the three TKIs at controlling leukemic progenitor cells, while nilotinib is the most effective at controlling the differentiated cells, which comprise most of the total leukemic burden. Thus, we note that controlling the leukemic progenitor cell population is important in long-term treatment outcome. This is further supported by the observation that blast crisis emerges due the acquisition of additional mutations in CML progenitor cells (not stem cells) [3]. Our approach suggests that using combination TKI therapies may be a viable method of controlling this population. Our modeling suggests that it is best to reduce the progenitor cells early and then reduce the differentiated cells towards the end of the treatment planning horizon. An early reduction in progenitor cells pays off in later stages of the treatment planning horizon, since a small progenitor cell population results in a lower growth rate for differentiated cells which leads to a greater response to subsequent TKI therapy. Effects of toxicity constraints. We also imposed a toxicity constraint on therapy optimization procedure by mandating that patient ANC levels stay above a given threshold that reduces the risk of infections. We observed that incorporating this toxicity constraint does impact the structure of the optimal schedules significantly, resulting in mandated treatment breaks as well as switching some months to imatinib therapy, which has a lower toxicity effect. We also noted that the choice of treatment breaks occurring also in one-month intervals may result in dangerous rebound of leukemic burden to levels close to pre-treatment, suggesting that shorter breaks to combat toxicity may be recommended. Although the model we have used for describing the dynamics of the ANC levels is simple, our findings demonstrate that incorporating a mechanistically modeled toxicity constraint into optimization of therapy scheduling is both feasible and important in determining optimal scheduling. Multiple preexisting mutant types. While some previous studies have suggested that the majority of CML patients are diagnosed with 0 or 1 preexisting BCR-ABL mutations, some patients do harbor multiple mutants at the start of therapy [2,6]. Thus we also investigated the impact of having 2 mutant types present (M351T and F317, or E255K and F317L) at the start of therapy, on optimal schedules. We observed the number and specific combination of preexisting mutants present can significantly impact the optimization results. This points to the importance of determining which BCR-ABL mutations preexist in patients at diagnosis, before treatment planning is done. Throughout this work we have observed that the structure of the optimal therapy depends heavily on model parameters, e. g. , cellular growth rates and ANC decay rates. It is likely that each individual patients will have unique model parameters, and therefore a unique best schedule. An exciting application of this work would be the development of personalized optimal therapeutic schedules. Determination of (i) the mutant types (if any) present in a patient’s leukemic cell population, (ii) growth kinetics of their leukemic cell populations, and (iii) patient ANC level responses under various TKIs, would enable our optimization framework to build treatment schedules in a patient-specific setting. At the start of treatment for acute lymphoblastic leukemia (ALL), Quantitative RT-PCR or similar techniques are sometimes used to perform mutational analysis to identify the preexisting mutant types. Indeed, studies have demonstrated that BCR-ABL mutants are present at the time of diagnosis in many ALL patients, and as sequencing technologies improve, smaller and smaller subclones with resistant phenotypes will likely be discovered [29]. In addition, the ANC level can be tested weekly by a routine complete blood count and in principle, the growth kinetics of each leukemic cell type could be analyzed in the laboratory using standard techniques such as flow cytometry or quantitative imaging. However, very few laboratories currently have the ability to analyze the growth kinetics of each type of leukemic cell population in a reproducible way. We believe that the testing procedure needs to be standardized before it becomes helpful for the treatment of CML. If such a standardized growth kinetic analysis can be realized, during treatment a patient’s response to various TKIs could be monitored so that the impact of TKI on growth kinetics of leukemic cell population and ANC levels could be quantified. Then the optimization procedure could be re-run on the fly with updated patient parameters, providing dynamic feedback into each patient’s optimal therapy schedule.
Targeted therapy using imatinib, nilotinib or dasatinib has become standard treatment for chronicle myeloid leukemia. A minority of patients, however, fail to respond to treatment or relapse due to drug resistance. One primary driving factor of drug resistance are point mutations within the driving oncogene. Laboratory studies have shown that different leukemic mutants respond differently to different drugs, so a promising way to improve treatment efficacy is to combine multiple targeted therapies. We build a mathematical model to predict the dynamics of different leukemic mutants with imatinib, nilotinib and dasatinib, and employ optimization techniques to find the best treatment schedule of combining the three drugs sequentially. Our study shows that the optimally designed combination therapy is more effective at controlling the leukemic cell burden than any monotherapy under a wide range of scenarios. The structure of the optimal schedule depends heavily on the mutant types present, growth kinetics of leukemic cells and drug toxicity parameters. Our methodology is an important step towards the design of personalized optimal therapeutic schedules for chronicle myeloid leukemia.
medicine and health sciences
pathology and laboratory medicine
cancer treatment
cell differentiation
toxicology
oncology
optimization
toxicity
developmental biology
mathematics
pharmaceutics
stem cells
animal cells
stem cell therapy
differential equations
clinical genetics
cell biology
biology and life sciences
cellular types
physical sciences
drug therapy
Optimized Treatment Schedules for Chronic Myeloid Leukemia