Combination of Two but Not Three Current Targeted Drugs Can Improve Therapy of Chronic Myeloid Leukemia

Chronic myeloid leukemia (CML) is a cancer of the hematopoietic system and has been treated with the drug Imatinib relatively successfully. Drug resistance, acquired by mutations, is an obstacle to success. Two additional drugs are now considered and could be combined with Imatinib to prevent resistance, Dasatinib and Nilotinib. While most mutations conferring resistance to one drug do not confer resistance to the other drugs, there is one mutation (T315I) that induces resistance against all three drugs. Using computational methods, the combination of two drugs is found to increase the probability of treatment success despite this cross-resistance. Combining more than two drugs, however, does not provide further advantages. We also explore possible combination therapies using drugs currently under development. We conclude that among the targeted drugs currently available for the treamtent of CML, only the two most effective ones should be used in combination for the prevention of drug resistance.

Suppose that cancerous cells divide with rate l s and die with rate D s . These coefficients can depend on the treatment dose; in particular, the death rate of cells is comprised of the "natural" rate of cell death in an untreated tumor and the action of the drug(s), if any, upon the cells. The mutation rates are defined for all the edges of diagram of figure 1, by u i→j , where the subscript for each mutation rate u represents the binary indices of the starting and the end point of the mutation process. here, u out s→j denotes the mutation rates for all arrows originating at s. Note that here we do not consider the processes of quiescence and awakening of cells (see [4], [3]). While important for the timing of treatment, these processes do not contribute to the overall probability of treatment success. Therefore, for the purposes of this study, they are omitted. Following the standard technique for transport-type equations, we have the following system of equations for characteristics: If initially there are M 0 cells of resistance type 0, then we have that is, to find the function Ψ(ξ 0 , . . . ,ξ n ;t), we need to solve the characteristics equations with the initial conditions ξ s (0) =ξ s , and evaluate Ψ(ξ 0 , . . . ,ξ n ;t) = [ξ 0 (t)] M0 .
In order to calculate the probability of treatment success, which is the same as the probability of extinction of the colony, we need to find the function value: lim t→∞ ϕ 0,...,0 (t) = lim t→∞ Ψ(0, . . . , 0; t).
The probability of treatment success. We assume that treatment starts at time t * . The time t * is assumed to correlate with the size N given by N = s x s (t * ), where the values x s are the expected numbers of cells of each type. These are found from the following ODEs for the mean values: In this paper we make the assumption that all the division and death rates in the absence of treatment are equal to each other. This allows us to solve the above equations very easily to obtain the following simple connection between the time and the colony size: In order to calculate the probability of treatment success, we need to solve the equations for characteristics twice. The time axes is divided into two regions: 0 < t < t * is pre-treatment, then t > t * is the treatment stage. First, we solve the equations for the treatment stage, that is, we find the limit lim t→∞ ξ s (t). Then we use these values as initial conditions and solve the system for the pretreatment phase, for the duration t = t * . The resulting value of ξ 0 (t * ) is used to calculate lim The systems for characteristics are solved numerically by using a standard Runge-Kutta ODE solver. The probability of treatment success given that the colony has grown to size N , is obtained as follows: This is because the probability of extinction ((ξ 0 (t * )) M0 ) equals probability of going extinct before reaching size N ((d/l) M0 ) plus the probability to grow to size N and to be treated (P (treat. succ.)).
Symmetry assumptions. In our calculation we will use several symmetry assumptions, which allow for a very concise representation of results, but are not necessary in principle. As mentioned above, we assume that all the division and natural death rates are equal to each other. The other assumption concerns the death rates of cells. Let us suppose that we apply treatment with m drugs.
In principle, the more drugs the cell is susceptible to, the higher is its death rate. However in this paper we will assume that if a cell is susceptible to at least one of the m drugs, it dies with the same rate as a fully susceptible cell. This death rate is equal to the natural death rate of untreated cells, d, plus a drug-related factor, h. On the other hand, if a cell is resistant to all the drugs, then its death rate is equal to that of the untreated cells, d.
Modeling cross-resistant mutations. The existing medical knowledge on the types of mutations for CML drugs is incorporated in the rates of mutations u. Here we illustrate it with the example of the drugs Imatinib, Dasatinib and Nilotinib. We assume that k i mutations give rise to resistance to drug i without affecting resistance properties with respect to the other drugs (the index i takes values 1, 2, 3 for the three drugs). This translates into u i = k i µ, where µ is the basic point mutation rate. In our simulations we took k i to be the same for all three drugs, k i = k. This is another symmetry assumption which simplifies the system but is not necessary in principle. For the lack of exact measurements of this parameter, we experimented with different values of k i in a biologically plausible range. Further we include the cross-resistant mutation T315 by setting u 000→111 = µ. There are no known mutations which confer resistance to only two of the three drugs; therefore, we set the rate of mutation to be zero for the following processes: 000 → 011, 000 → 101, 000 → 110. Note that the T315 mutation contributes to the following processes: 001 → 111, 010 → 111, 100 → 111 (the corresponding mutation rates are set to be equal µ), and 011 → 111, 101 → 111, 110 → 111 (the corresponding mutation rates are set to be equal (k + 1)µ).
We have also investigated the role of two-drug cross-resistant mutations. We assumed that there is a mutation event which gives rise to resistance against drugs 2 and 3; another mutation event -against drugs 1 and 2, and a third one -against drugs 1 and 3. The corresponding mutation rates are denoted as u T , u B and u A , see diagram of figure 2. For simplicity we assume that all onedrug resistant mutations happen at the same rate, u 1 . In general, there are five possible three-drug cross-resistance networks, see figure 3. In the figure, each drug is represented as a node of the network, and a connection corresponds to the existence of a mutation which confers cross-resistance to both drugs. Different types of line (dashed, solid, double) correspond to different mutants. If the lines are the same then triple cross-resistance takes place. In figure 3 we provide examples of the drugs which participate in each cross-resistance network. There, "I" stands for Imatinib, "D" for Dasatinib and "N" for Nilotinib. "K" stands for one of the drugs still in development which is active against cells with a T315I mutation, see main text. By setting some of the cross-resistance rates to zero, we can model four of the five possible cross-resistance networks Back-mutations. No back mutations are included in the above equations. We have performed a separate study where all possible back-mutations were taken into account. This results in additional terms in the equations for characteristics. However, for biologically reasonable values of the mutation rates, the inclusion of back mutation did not make a measurable difference in the probability of treatment success.
A note on the statistical significance of our results. The stochastic model we use in this paper is a very convenient tool to convey comparative studies of different scenarios and treatment strategies. In this context, our method has significant advantage compared to, say, directly running stochastic simulations. If we were running stochastic simulations, and obtained two different results for two treatment strategies, then the question would be if the difference is "real" (or significant), or if it is only a result of stochastic fluctuations. This question does not arise with our method. Instead of running stochastic simulations, we operate directly with probabilities. For example, in figure 1 of the main paper we plot the probability of treatment success as a function of the tumor load. The difference in this quantity for the two closest curves in that figure, say for log 10 N = 13.5, is about 0.15, which is a 15% increase in the probability of treatment success (here, we compare three-drug combinations with two and three pairs of cross-resistant drugs). There is no stochastic noise here to make the difference an artifact of fluctuations, so we can interpret it directly. Therefore, the issue of "significance" in this context is not a mathematical or statistical question anymore, but rather a moral and financial issue, to be decided in accordance with some guidelines.