Effect of Cellular Quiescence on the Success of Targeted CML Therapy

Background Similar to tissue stem cells, primitive tumor cells in chronic myelogenous leukemia have been observed to undergo quiescence; that is, the cells can temporarily stop dividing. Using mathematical models, we investigate the effect of cellular quiescence on the outcome of therapy with targeted small molecule inhibitors. Methods and Results According to the models, the initiation of treatment can result in different patterns of tumor cell decline: a biphasic decline, a one-phase decline, and a reverse biphasic decline. A biphasic decline involves a fast initial phase (which roughly corresponds to the eradication of cycling cells by the drug), followed by a second and slower phase of exponential decline (corresponding to awakening and death of quiescent cells), which helps explain clinical data. We define the time when the switch to the second phase occurs, and identify parameters that determine whether therapy can drive the tumor extinct in a reasonable period of time or not. We further ask how cellular quiescence affects the evolution of drug resistance. We find that it has no effect on the probability that resistant mutants exist before therapy if treatment occurs with a single drug, but that quiescence increases the probability of having resistant mutants if patients are treated with a combination of two or more drugs with different targets. Interestingly, while quiescence prolongs the time until therapy reduces the number of cells to low levels or extinction, the therapy phase is irrelevant for the evolution of drug resistant mutants. If treatment fails as a result of resistance, the mutants will have evolved during the tumor growth phase, before the start of therapy. Thus, prevention of resistance is not promoted by reducing the quiescent cell population during therapy (e.g., by a combination of cell activation and drug-mediated killing). Conclusions The mathematical models provide insights into the effect of quiescence on the basic kinetics of the response to targeted treatment of CML. They identify determinants of success in the absence of drug resistant mutants, and elucidate how quiescence influences the emergence of drug resistant mutants.

The Kolmogorov forward equation for this system is: (1) The initial condition is such that the number of cycling cells is given by I 0 and the number of quiescent cells is given by J 0 at time t = 0, i.e. ϕ I0,J0 (0) = 1, and ϕ ij (0) = 0 for all other values of i and j.

Equations for the averages
Let us denote by x the average number of cycling cells, and y the average number of quiescent cells: The equations for x and y are as follows: x = (l − d − α)x + βy, (3) y = αx − βy, (4) x(0) = I 0 , y(0) = J 0 . (5) This linear system of ODEs can be solved analytically. In particular, the eigenvalues of the corresponding matrix are: such that |λ + | ≤ |λ − |.
The total number of cells, n tot (t) = x(t) + y(t), as a function time, is given by with the constants The function log n tot (t) has two asymptotes: log n tot → log g − + λ − t as t →inf ty, and log n tot → log g + + λ + t as t → +∞, see figure 1. It is possible to show that g + ≥ 0. There are three cases: (a) g − > g + , (b) 0 ≤ g − ≤ g + and (c) g − < 0. If condition (a) is satisfied, then the two asymptotes intersect at a positive t, and we have a biphasic decline (see below), where an initial fast phase of decay is followed by a slow phase. If (b) holds, there is essentially one wave of decline (the first wave is not pronounced). Finally, if (c) holds, we have a "reversed biphasic decline", where the function log n tot decays slower at first and then faster.

Biphasic decline
If the following inequality holds 1 : then the function n tot (t) has a characteristic shape of a "biphasic decline". There are two waves of decay, first the fast one at the rate λ − , and then the slower one, at the rate λ + . Note that if α (d + β − l), the two eigenvalues are −β and −(d − l). The time of switching from the first to the second phase is obtained from the equation g − e λ−t = g + e λ+−t , such that log n tot log n tot log n tot Figure 1: The three types of population decline: (a) Biphasic decline, (b) One phase of decline and (c) Reverse biphasic decline. The thick line represents the function log ntot, the thin solid line is log g + + λ + t, and the dashed line is log g − + λ − t. In (c), the dashed line is absent because g − < 0.
In the case where J 0 I 0 , we have In particular, if α is small, we have the following simplified expression:

The time until extinction
The time it takes to reduce the total population to small numbers can be estimated by solving the equation x(t) + y(t) = 1 for t. Approximately, this is given by g + e λ+t = 1, or A similar expression for the time of treatment is obtained from studying the probability of extinction, ϕ 0,0 (t), equation (1).

Fitting the model with the data
We have fitted out model with the experimental data of [4] and [6]. The data shows the relative number of CML cells in patient's blood. More precisely, [4] plots the ratio of BCR-ABL transcripts to BCR transcripts, while [6] shows the ratio of BCR-ABL transcripts to ABL transcripts in peripheral blood of patients treated with Imatinib. Formula (7) was used to fit the data. Although the formula was derived for the absolute number of cells, it is capable of measuring the ratios, as long as the "base" measurement remains constant. Thus it is assumed that in the experiments, the number of non-leukemic cells remains roughly constant throughout the experiment.
We used a standard least-square algorithm to approximate the logarithm of the data with a piecewise-linear function representing the two stages of decline. To simplify the fitting procedure, we used the assumption that α (d + β − l). In this case the first slope yields the value of d − l, and the second slope yields the parameter β. Next, we had to extract the value of the two remaining parameters, α and J 0 (the value of I 0 is obtained from I 0 = N 0 − J 0 , where N 0 is the measurement at time zero). As α decreases, the error of the fit first decreases and then flattens out and stays essentially constant (and small) for small values of α. The corresponding values of J 0 grow as α decays, but also saturate at a certain low level. This suggests that the shape of the curve is not sensitive to the exact value of α, as long as it is sufficiently small.
These observations make sense because smaller values of α (a smaller production rate of quiescent cells) can be somewhat compensated by higher initial numbers of quiescent cells. Moreover, even though in formula (7), the values of α and J 0 are independent parameters, in reality they are strongly linked. As the colony grows, the amount of quiescent cells at size N , J 0 , is defined by the rate of quiescence, α, according to (20) below. Thus in order to correctly determine the parameter α, we would need to use the growth-phase of the cell colonies.
To conclude, we note that the least-square procedure yields estimates for the decay rate of cells, d − l, the rate of awakening, β, and the ranges for α and J 0 . These calculations demonstrate that the model is capable of capturing the behavior of CML cells treated by Imatinib, and the parameters of the biphasic decay can be determined from the quiescence rates. However, the analysis presented here is only an illustration of the method. The fitting procedure should be repeated for future data sets which show the mean rather than the median, and also include information on the absolute numbers of non-CML cells.
2 Modeling resistance to drugs in a population with quiescence Here we include the possibility of mutations in the birth-death-quiescence process. This is an extension of the framework of [1,2], where generation of resistance by cells without quiescence was studied.

Stochastic model
A cell can be either resistant or susceptible to each of the m drugs. In total, there are 2 m − 1 resistant types. They can be enumerated by a binary index s, which is a string of m components with numerical values from the set {0, 1}. We denote by i s the number of cycling cells of resistance type s, and j s the number of quiescent cells of resistance type s. Each type s is characterized by its division rate, l s , and death rate, d s . Also, all the types are naturally placed on a binary directed network where each edge connects two types which only differ by their resistance properties with respect to one drug, see [1]. A single mutation separates each pair of connected types. The Kolmogorov forward equation for this system can be written (not shown). The following probability generating function can be defined: . . ξ in n η jn n (the expression above can be treated as a transformation from discrete variables i 0 , j 0 , . . . i n , j n to continuous variables ξ 0 , η 0 , . . . , ξ n , η n ). This function satisfies the PDE (see e.g. [5], [7], Chapter 3; [2]): For each type s, the quantities u s→j are the mutation rates corresponding to the transformations of type s into each of the types j, with which s is directly connected. The types j are assumed to be downstream from type s. The linear transport-type PDE can be solved by the standard method of characteristics: If initially we have M 0 cycling cells of resistance type 0, then we have Ψ(ξ 0 , η 0 , . . . , ξ n , η n ; 0) = ξ M0 0 .
All the resistant types can be separated into classes such that in each class k, all the types are resistant to exactly k drugs (and susceptible to m − k drugs).
Let us suppose that within each class, the birth and death rates are equal, and also that all mutation rates are equal to each other. Then the equations for all types within the same class are identical, and we have for each class k: The coefficients in equations (11-12) depend on whether treatment takes place or not. We assume that before the start of treatment, the growth and death rates of all types are the same: After the treatment starts, we assume that the fully-susceptible and all the partially resistant types are killed by the drugs with the intensity H, and the fully-resistant type is not affected by treatment. This translates into the following values for the coefficients:

Probability of resistance generation during treatment
Suppose that at the start of treatment we have I 0 cells of class 0, . . ., I m cells of class m. Then the probability of treatment success is given by where ξ ∞ k = lim t→∞ ξ k (t). From studying the fixed points of system (11-12), we can see that the limiting values ξ ∞ k do not depend on the quiescence parameters α and β. This means that the probability of treatment success (for a given colony) is not affected by the process of quiescence.

2.3
The probability of resistance generation before treatment.
The probability of not having resistance in a colony of size N can be estimated as follows. We assume that the size N of the colony and the time t are connected by the deterministic relation, x(t N ) + y(t N ) = N , where, as before, x(t) is the expected total number of cycling cells (of all types), and y(t) is the expected total number of quiescent cells (the validity of this assumption is discussed in [3]). The functions x(t) and y(t) satisfy the system (3-4) with the initial condition x(0) = M 0 , y(0) = 0.
The probability of not having resistance in a colony of size N is then given by This follows from the fact that the probability of having no resistant mutants ([ξ 0 (t N )] M0 ) equals the probability of going extinct before reaching size N ((d/l) M0 ), plus the probability to grow to size N and to have no resistant mutants (P rob(no resistance at size N )).
3 The average number of mutants at size N Another way to study the production of resistant mutants is to look at the expected behavior (the first moments) of various types of cells. We can derive equations for the moments from the Kolmogorov forward equations.

Equations for the moments
To gain intuitive understanding of the dependence of mutant production upon quiescence, we consider a simplified system for m drugs with d = 0: Here, the rates u k are probabilities of mutation from class k to class k + 1.
Because of the combinatorial nature of our problem, we have where u is a basic rate of mutations.

The expected number of mutants for β = 0
Let us assume that β = 0. Then linear system (13-17) can be solved analytically (the simplifying assumption that β = 0 makes the equations for cycling cells independent of the equations for quiescent cells). In particular, for any value of k, we have the exact solution Therefore, the time it takes to reach size N is given by To calculate the fraction of cycling cells of class k at size N , we substitute the above expression for t N into the exact solution for x m (t) + y m (t), and expand in a double Taylor series in terms of small parameters u and 1/N . We obtain: Similarly, the average number of mutants resistant to m drugs is given by Again, here we only present the first term in the u-expansion and the largest term in the 1/ log N expansion of the exact expressions. 2 Exact values of quantities x k (t N ) + y k (t N ) and x k x k +y k are plotted in figure 2(a,b) as functions of α.
Let us first consider the case m = 1. The average number of (cycling and quiescent) 1-hit mutants in this case is an increasing function of α, see equation (21). On the other hand, the average number of cycling 1-hit mutants is obtained from equations (21) and (20) and is given by N u log N , that is, it is independent of α. To understand this result, we note that cycling mutants are produced by cycling wild-type cells and they grow according to the same law as the cells producing them. When α increases, the mutant clones grow more slowly because of quiescence, but at the same time they have more time to grow (it takes more "events" before the number of cells reaches size N ). In other words, the changes in the mutant growth are completely compensated by the change in the time of growth. That is why the quantity x 1 (t N ) is a function of the colony size only, see figure 2(c).
The dependence of the total number of mutants, x m (t N ) + y m (t N ), on α becomes stronger as the number of drugs increases, see figure 2. In particular, both the total number of 2-hit mutants and the number of cycling 2-hit mutants depend on α, formula (21).

The probability of having resistance to m drugs, β = 0
In our studies of resistance generation, two quantities are of particular importance: the expected number of resistant mutants, and the probability of having resistance.
Of course, the latter quantity (the probability of having resistance) cannot be obtained directly when studying system (13-17). However, we can study a related quantity is W m−1 ≡ tN 0 lux k (t) dt, the total number of fully resistant mutants produced by the colony. We have In particular, for m = 1, both W 0 and the probability of generating resistance to one drug is independent of α. However, the dependence on α does take place for larger values of m. For example, the probability of resistance to 2 drugs is an increasing function of α. The case of m = 2 drugs can be explained intuitively. Two-hit mutants are produced by cycling one-hit mutants. At each value of N , the number of cycling wild-type cells is given by that is, it decreases with α, and the number of cycling one-hit mutants is which is α-independent (the last two expressions are obtained by a direct solution of system (13-17) and using equation (19)). We can see that the relative proportion of one-hit mutants for a given colony size is a growing function of α. Therefore, for increasing α, every time a cell divides, it becomes more likely that the dividing cell already has a mutation. Therefore, the production of two-hit mutants is a growing function of α, see the expression for W 1 , equation (22).

The case β > 0
Once we include a nonzero value of β in system (13-17), the equation m k=0 (x k + y k ) = N cannot be analytically resolved for t N . In other words, there is no equivalent of formula (19) with β > 0. In order to see how the results depend on the value of β, we have performed numerical solutions of system (13-17). These are presented in figure 3. We can see that the total number of mutants of class k decays with β ( figure 3(a)), and the fraction of cycling cells increases with β ( figure 3(b)). The number of cycling mutants, x k , is a decaying function of β for k > 1 ( figure 3(c)).
These results show that the dependence of the number of resistant mutants on the rate of cell awakening, β, is the opposite of that on the rate of quiescence, α. In general, the number of resistant mutants increases with the amount of quiescence in the system. The probability of resistance to k drugs (which is related to the function x k−1 , formula (22)), also increases with the amount of quiescence, except for the case k = 1, where it is independent of the amount of quiescence. This illustrates the results of stochastic calculations of Section 2 (see figure 2 of the main paper) at an intuitive level.