Optimized Treatment Schedules for Chronic Myeloid Leukemia

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.

• The birth rates b j 1,i . The estimates b j 1,1 = 8.00 × 10 −3 and b j 1,i = 1.00 × 10 −2 for any cell type i ≥ 2 and drug j. The value 1.00 × 10 −2 is used in [1] for the birth rate of leukemic stem cells without drug. We further assume that this value remains the same under any therapy, which is different from [1].
• The death rates d j 1,i . The estimate d j 1,i = 5.00 × 10 −4 for any i and j, from [1].

Progenitor cell kinetics
• The death rates d j 2,i . The estimates d 1 2,i = 2.80 × 10 −3 and d 2 2,i = 5.30 × 10 −3 for any i, from [1]. The death rate of leukemic progenitor cells under high-dose imatinib (800 mg/day) is 3.50 × 10 −3 in [1]. We consider imatinib with regular dose (400 mg/day) in this paper, so we set d 3 2,i = 3.5 × 10 −3 ÷ 2 = 1.75 × 10 −3 for any i. Note the death rates are the same across all cell types with the same therapy, but vary with different therapies. In addition, we set the death rate of normal progenitor cells d 0 2,i = min{d 1 2,i , d 2 2,i , d 3 2,i } = 1.75 × 10 −3 for any i.
• The differentiation rates b j 2,i .
-For mutants, the differentiation rates are listed in Table 1. Since there are little in vivo data available in the literature related to leukemic mutant birth rate, our estimation is based on in vitro data for these mutants, in particular the IC50 values. We use a piecewise linear interpolation to estimate the differentiation rates, based on the relative IC50 values of mutants under nilotinib, dasatinib, and imatinib reported in [2]. For sensitive or moderately resistant mutants (the relative IC50 value is less than or equal to 4), the differentiation rate of mutant i is estimated using the linear interpolation For resistant mutants (the relative IC50 value is between 4.01 and 10.0), the differentiation rate is estimated with the following linear interpolation: Thus if the relative IC50 value for a resistant mutant is 4.01, then its differentiation rate is 90% of the differentiation rate of the WT cell without any drug (0.700 per day); if the relative IC50 value for a resistant mutant is 10.0, then its differentiation rate is equal to the birth rate of the WT progenitor cell without drug. For highly resistant mutants (the relative IC50 is larger than 10.0), we set its differentiation rate to the differentiation rate of the WT progenitor cells without drug.

Differentiated cell kinetics
• The death rate d j 3,i . The estimates d 1 3,i = 0.0442 and d 2 3,i = 0.0394 for any i, from [1]. The death rate of leukemic differentiation cells under high-dose imatinib (800 mg/day) is 0.0550 in [1]. We consider imatinib with regular dose (400 mg/day) in this paper, so we set d 3 3,i = 0.0550/2 = 0.0275 for any i. In addition, -For normal cells, b j 3,1 = 5.50 for any j. The estimate is from [1].
-For the mutant, if it is sensitive or moderately resistant to drug j (the relative IC50 value is less than or equal to 4), then b j

Terminally differentiated cell kinetics
Using the estimates from [1], we set the differentiation rates b j 4,i = 100 and death rates d j 4,i = 1.00 for any i and j.

Initial cell populations at diagnosis
The normal marrow output in an adult is approximately 3.50 × 10 11 cells per day [3]. To achieve this equilibrium condition, we set K 1 = 8.75 × 10 4 in differential equations from equation (1) of the main text with parameters described in Sections "Stem cell kinetics" to "Terminally differentiated cell kinetics" and in the absence of leukemic cells.
To obtain an estimate of K 2 , we assume that diagnosis of CML occurs once the leukemic cell burden reaches a threshold of 10 12 cells [4], and that the differential equations in equation (1) of the main text have parameters described in Sections "Stem cell kinetics" to "Terminally differentiated cell kinetics" and start with K 1 = 8.75 × 10 4 normal stem cells, one wild-type leukemic stem cell, and no other cells. We set K 2 = 3 × 10 6 so that the patient is diagnosed with CML around 78 months (6.5 years) after the first leukemic stem cell arises. At diagnosis, the normal stem cell, progenitor cell, differentiated cell, and terminally differentiated cell populations are 7.34 × 10 4 , 1.61 × 10 7 , 3.24 × 10 9 , and 3.24 × 10 11 , respectively; the leukemic stem cell, progenitor cell, differentiated cell, and terminally differentiated cell populations are 2.95 × 10 5 , 4.07 × 10 7 , 1.08 × 10 10 , and 1.08 × 10 12 respectively. These are used as the initial cell populations for a patient diagnosed with CML.

ANC kinetics
• We require that the patient's ANC cannot fall below L anc = 1000/mm 3 . The normal range of the ANC is 1500 to 8000/mm 3 . We assume that the normal level of ANC is U anc = 3000/mm 3 , and the patient's initial ANC is 3000/mm 3 .
• We initially set the monthly decrease rates of ANC to be d anc,1 = 350/mm 3 under nilotinib, d anc,2 = 300/mm 3 under dasatinib, and d anc,3 = 250/mm 3 under imatinib. For imatinib, in [5] it was observed that the median time until the first neutropenia event was 74 days. Given the large range of possible ANC values (1500-8000) at the start of the trial this gives a range of 200/mm 3 to 2800/mm 3 for feasible ANC decrease rates under imatinib. For nilotinib, the reference [6] observes that nilotinib has a higher risk of adverse hematological events than imatinib. Furthermore [6] observes a median time of 24 days until myelosuppresion events. Given that these events may have included lower-grade hematological events we assumed that the median time was an underestimate for the time until neutropenia and thus chose d anc,1 = 350/mm 3 . We were unable to find any published data on the median time until neutropenia during dasatinib therapy. However, it has been observed that dasatinib induces grade 3 and 4 neutropenias at a higher rate than imatinib [7]. Based on this we initially set dastinib to have a toxicity in between imatinib and nilotinib, but then also considered the scenario where dasatinib has a higher toxicity than nilotinib, i.e., d anc,1 = 300/mm 3 and d anc,2 = 350/mm 3 . It usually takes less than a month for the ANC of a patient with grade 3-4 neutropenia to go back to normal level, so the monthly increase rate of ANC level is between 500 and 7000/mm [5]. In our model, we assume that the ANC of a patient increases by b anc = 2000/mm 3 during a drug holiday, before it reaches the normal level 3000/mm 3 .

B: Method to solve the optimization model
We describe the method to solve the optimization model introduced in the Treatment optimization problem section (see equations (2)-(6) of main text) with toxicity constraints given by equation (7) of the main text and y m ≥ L anc for m ∈ M. Our strategy is to build a mixed-integer linear optimization model [8] that approximates the optimization model given by equations (2)-(6) of the main text, and then solve the approximation model to optimality numerically by off-the-shelf optimization software CPLEX [9]. The mixed-integer linear optimization model is built through two steps: (1) we first approximate the ODE constraints (equation (3) of main text) by bilinear constraints; (2) we then transform the bilinear constraints and nonlinear toxicity constraints (equation (7) of the main text) into equivalent linear constraints, by adding auxiliary decision variables. We first describe how to approximate the ODE constraints (equation (3) of main text) by bilinear constraints. Suppose patients take drug j in month m. Since the cell dynamics are modeled by the following set of ODEṡ the cell abundances in month m + 1, x m+1 , are completely determined by the initial cell abundance x m and function f j . Without loss of generality, we assume this relationship is described by x m+1 with some unknown nonlinear function g j l,i : R L×n → R, for each month m, layer l, and cell type i. Recall that L is total number of cell layers (L = 4), and n is the total number of cell types. Then the ODE constraints (equation (3) We will approximate the nonlinear function g j l,i with an affine function g m,j l,i : R L×n → R, for each m, j, l, and i. In particular, the function where a j,l,i is an (Ln)-dimensional vector and does not depend on m. Details of hoŵ g m,j l,i is constructed are provided in Section "linear approximation to the solutions of the ODEs" of the Appendix. Let a j,l,i = [a j,l,i 1,1 , . . . , a j,l,i k,s , . . . , a j,l,i L,n ]. Then constraint (2) can be approximated by the bilinear constraint for each type i cell at layer l in month m. We now describe how to transform bilinear constraints (4) and piecewise linear constraints (equation (7) of the main text) into linear constraints. These are standard techniques in mixed-integer linear optimization [8]. We introduce auxiliary continuous variables v m,j k,s , and set v m,j k,s = z m,j x m k,s . Then bilinear constraints (4) are transformed into the equivalent linear constraints below.
x m+1 , where U k,s is an upper bounds of cell abundance x m k,s for each m. The value of U k,s can be obtained by taking the maximum value of layer k type s cell abundances over the whole planning horizon under all three monotherapies and no treatment. The piecewise linear constraints can be transformed into equivalent linear constraints below, by introducing auxiliary continuous variable u m and binary variable q m for each m.
Overall, the optimization model (equations (2)-(6) of main text) with toxicity constraints (equation (7) of main text and y m ≥ L anc for m ∈ M) is approximated by the following mixed-integer linear optimization model.
x(0) = x 0 , y 0 is given. (6) C: Linear approximation to the solutions of the ODEs In this section, we describe how to construct the affine functionĝ m,j l,i in (3) in Section "Method to solve the optimization model" of the Appendix. If we assume that the drugs do not affect stem cells, we can compute the abundance of stem cells over the planning horizon numerically in advance, regardless of the treatment schedules. Thus we assume x 1,1 (t), . . . , x 1,n (t) are given as data, for any t. We can first eliminate all the variables x m 1,i and constraints containing x m 1,i , for each m and i, in the optimization problem (6). The dynamics of wild-type leukemic cells and each mutant type have no impact on each other. We can decouple the ODEs into a series of linear ODEs as follows, each describing the dynamics for type i cell from layer 2 to layer 4.
Write the above equations in the matrix form, we havė v i (t) = W j i v i (t) + w j i (t), for t ∈ [m∆t, (m + 1)∆t], where v i (t) = [x 2,i (t), x 3,i (t), x 4,i (t)] , w j i (t) = [b j 2,i x 1,i (t), 0, 0] , and W j i is the lower triangular matrix in (7).