DATA AND IMPLICATION BASED COMPARISON OF TWO CHRONIC MYELOID LEUKEMIA MODELS

Chronic myeloid leukemia, a disorder of hematopoietic stem cells, is currently treated using targeted molecular therapy with imatinib. We compare two models that describe the treatment of CML, a multi-scale model (Model 1) and a simple cell competition model (Model 2). Both models describe the competition of leukemic and normal cells, however Model 1 also describes the dynamics of BCR-ABL, the oncogene targeted by imatinib, at the sub-cellular level. Using clinical data, we analyze the differences in estimated parameters between the models and the capacity for each model to predict drug resistance. We found that while both models fit the data well, Model 1 is more biologically relevant. The estimated parameter ranges for Model 2 are unrealistic, whereas the parameter ranges for Model 1 are close to values found in literature. We also found that Model 1 predicts long-term drug resistance from patient data, which is exhibited by both an increase in the proportion of leukemic cells as well as an increase in BCR-ABL/ABL Model 2, however, is not able to predict resistance and accurately model the clinical data. These results suggest that including sub-cellular mechanisms in a mathematical model of CML can increase the accuracy of parameter estimation and may help to predict long-term drug resistance.


Introduction
Chronic myeloid leukemia (CML) is a cancer of the white blood cells. It is a disorder of hematopoietic stem cells characterized by the increased growth of myeloid cells in the bone marrow and the excessive presence of these cells in the blood. CML can be molecularly diagnosed by detecting the presence of the Philadelphia (Ph) chromosome and the fusion oncogene BCR-ABL. This oncogene is the result of translocation of the BCR, or breakpoint cluster, gene located on chromosome 22 and the ABL, or Ableson leukemia virus, gene located on chromosome 9 [3] . The progression of CML consists of three phases. The first phase, called the benign chronic phase, is typically asymptomatic and can last for several years untreated. The accelerated phase then follows and leads to the last phase, called blast crisis, which is characterized by an abnormally high number of stem cells and precursor cells in the blood or bone marrow [1]. CML can advance from the chronic phase to the fatal blast crisis phase in a timespan of 3 to 5 years [3].
For CML patients in which the BCR-ABL oncogene is detected, targeted molecular therapy can be used to inhibit the growth of stem cells. Imatinib, also known as STI-571 and Gleevec [11], is a tyrosine kinase that binds to the ATP binding site of BCR-ABL kinase, stopping cell-growth signals and decreasing cell proliferation [5]. Previously, treatment options included other drugs, such as hydroxyurea or interferon alpha, and allogeneic bone marrow transplants. Imatinib, which is effective in all phases of CML progression, is now a widely used primary treatment option for BCR-ABL positive CML patients [1].
BCR-ABL transcript levels are obtained using quantitative reverse transcription polymerase chain reaction for diagnosis and to determine the effect of treatment [17]. The amount of BCR-ABL transcript is then normalized using some control gene, usually BCR or ABL. Thus the data is presented as BCR-ABL/control gene percentages (BCR-ABL/ABL%) [17,14].
Most CML patients with imatinib treatment exhibit a biphasic profile, which means that the patients exhibit an initial rapid decline followed by a gradual decline in BCR-ABL/ABL%. However, some patients exhibit a monophasic or triphasic profile. The BCR-ABL/ABL% for a monophasic profile gradually decline over time. Patients with a triphasic profile can exhibit a rapid BCR-ABL/ABL% decline followed by a relatively gradual BCR-ABL/ABL% decline, followed by a rapid BCR-ABL/ABL% increase [18]. The increase in the triphasic profile is most likely caused by mutations in the BCR-ABL gene that encode resistance to Imatinib [18,5,6,7], although gene amplification is also a possible cause for resistance [6]. When treatment is stopped, the BCR-ABL/ABL% of some patients rapidly increases to levels at or beyond pre-treatment baseline [13]. There are different hypotheses as to the cause of this increase and mathematical modeling techniques have the potential to be helpful in elucidating the underlying mechanisms of therapy resistance.
Several groups have utilized mathematical modeling to study the effect of targeted treatment and imatinib on CML, including Roeder et al. and Michor et al. [17,13,12]. Roeder  Stein et al. [18] compared different hypotheses of the models described above. Assuming a biphasic decline in BCR-ABL/ABL%, which was demonstrated by both models, there are two slopes: α, which corresponds to the initial rapid decrease and β, which corresponds to the long-term response. One hypothesis, supported by Roeder et al., was the proliferating-quiescent hypothesis, where α is due to the proliferating stem cells and β is due to the quiescent stem cells. Another was the late-early progenitors hypothesis, supported by Michor et al., where α is due to the late progenitor cells and β is due to the early progenitor cells. Stein et al. also considered a third hypothesis, the early stem cell hypothesis, which states that α is due to a decline of early progenitor cells and β is due to a decline of late progenitor cells. Stein et al. rejected the late-early progenitors hypothesis and concluded that β is due to late progenitor depletion. However, the factors contributing to the parameter α are still unknown. These previous CML models described different cell environments and different cell populations, but not sub-cellular dynamics, which may provide insights into the sub-cellular origins of proliferation and resistance. Recently, Portz et al. [16] described a cell quota model that describes a treatment for patients with prostate cancer called intermittent androgen suppression, a hormone therapy. Normal prostate cells as well as most prostate cancer cells depend on androgen signaling for survival and proliferation. Androgen suppression treatment lowers the androgen levels, which prevents the growth of cancer cells. The treatment can be initially successful, however most patients experience a relapse. Portz et al. suggest that during the relapse, androgen-independent cells (AI), which can grow in low-androgen levels, replace androgen-dependent cells (AD). In their final model, the growth rate of both the AD and AI populations are described by Droops cell quota models, which introduce two new variables to represent the cell quotas for androgen. The switching rates between the two populations are modeled using the hill equations. The marker for prostate cancer, prostate-specific antigen (PSA), is assumed to be dependent on androgen levels in the model. This assumption was an important addition in the final model when comparing the model to clinical data. It is not always possible to determine the dependence of cancer cell phenotypes on sub-cellular factors. Portz et al. exemplified how mathematical modeling can yield insights into how sub-cellular dynamics may contribute to malignant cell growth.
Similar to prostate cancer, the proliferation of malignant cells in CML is dependent on the production of sub-cellular factors, namely the BCR-ABL protein. Therefore, it is natural to adapt the cell quota modeling approach of Portz et al. to CML. Additionally, the increases in BCR-ABL% that occur in some CML patients are comparable to the increases in androgen levels seen in prostate cancer patients following cessation of androgen therapy. Portz et al. showed that a cell quota model can accurately capture such increases in sub-cellular molecular factors that contribute to malignant proliferation [16]. Here, we compare two mathematical models, a cell quota model similar to that of Portz et al. [16] and a density dependent model based on a model described in Michor et al. [12,11], for the treatment of CML. Our results show that additional insights into imatinib treatment for CML patients can be gained by accounting for the dynamics of BCR-ABL at the sub-cellular level.

Model 1: A Cell Quota Model
Our goal is to gain a more in-depth understanding of CML and imatinib treatment by developing a model which may produce plausible solutions that reasonably match clinical data. Our model is based on the model framework of Portz, Kuang, and Nagy [16]. The growth rate of the BCR-ABL dependent and independent populations are modeled using the Droop's cell quota model, where Q(t) represents the cell quota for BCR-ABL. The BCR-ABL dependent, BCR-ABL independent, and normal populations are modeled respectively by the following: The BCR-ABL dependent population is equivalent to non-resistant cells and the BCR-ABL independent population is equivalent to the resistant cells. We assume that the proliferation rates, r i (1 − q i Q ), i = 1, 2, of both BCR-ABL dependent and independent populations are BCR-ABL cell quota dependent while the proliferation rate, , for the normal population is density dependent. Notice that r i , i = 1, 2, 3 are the corresponding maximum proliferation rates, and p 3 is the parameter that simulates the crowding effect. q i , i = 1, 2 are minimum BCR-ABL cell quota for BCR-ABL dependent and independent cells. We assume q 1 > q 2 since BCR-ABL independent cells are more likely to proliferate than BCR-ABL dependent cells in low BCR-ABL environment. r i (1 − q i Q ), i = 1, 2 implies that at minimum BCR-ABL cell quota (Q = q i ), corresponding leukaemia cells do not proliferate, while the proliferation rate increases and approaches the maximum as the BCR-ABL cell quota increases. The death rate, d 0 is also assumed to be the same for all stem cells.
The mutation or switching rates between the BCR-ABL dependent and independent populations are given by the hill equations The maximum BCR-ABL dependent to independent mutation rate is given by k 1 and similarly, the maximum BCR-ABL independent to dependent mutation rate is given by k 2 . K 1 and K 2 represent the BCR-ABL dependent to independent, and independent to dependent, mutation halfsaturation level respectively.
We assume the cell quotas for BCR-ABL for both the BCR-ABL dependent and independent cells are the same and are modeled by The maximum cell quota is q m1 and the minimum cell quota is q 1 , with q 1 > q 2 . The parameter v m represents the maximum cell quota uptake rate and the parameter v h represents the uptake rate half-saturation level. According to Abbott and Michor [1], the Ph chromosome contributes to the cell growth and so we assume BCR-ABL is used within the cells for growth at rate µ m . We assume BCR-ABL degrades at a constant rate b.
The dynamics of the intracellular BCR-ABL concentrations for the BCR-ABL dependent and independent populations under imatinib treatment are modeled by The state variable u(t) controls the administration of the BCR-ABL inhibitor treatment, which affects the BCR-ABL concentration within the BCR-ABL dependent and independent cells with the same efficacy. We assume that the rate of change of the intracellular BCR-ABL concentration is much higher than the rate of change of the cell populations. Applying the quasi-steady-state argument to B equation yields .
To simplify the model further, we assume that p(x 1 , x 2 ) is a constant.

Basic Analysis of Model 1
In the following we show that solutions of (2.1), (2.2), (2.3), and (2.6) with biologically appropriate initial values stay positive. Specifically, we assume , and all the parameters are positive. These assumptions are natural for our application.
Observe that x 1 and x 2 are dynamical variables dependent on the intracellular BCR-ABL concentration Q regulated by B which is assumed to be independent of x 1 and x 2 . Hence they are essentially assumed to grow or die exponentially. Mathematically this means that boundedness of x 1 and x 2 can not be established without additional assumptions.
It is easy to see that q m1 ≥ Q(t) ≥ q 1 µm µm+b for t > 0 with initial condition q m1 ≥ Q(0) ≥ q 1 . A straightforward application of standard comparison argument will establish the positivity of x 1 , x 2 and x 3 .
We consider now the boundedness of x 3 .
From this, we see that lim We are now in a position to consider the uniform boundedness of x 1 and x 2 .
Since Q(t) ≤ q m1 , we have amounts to say that even at the maximum intracellular BCR-ABL concentration q m1 , the populations x 1 and x 2 grow at a rate less than their death rate d 0 , which trivializes this modeling task. A much more natural mechanism that shall ensure the boundedness of solutions is the density dependent death rate. In more plausible CML models with more desirable long term dynamics, one can add an additional term such as d 1 x 2 1 to (2.1) and d 1 x 2 2 to (2.2). In the following, we assume that r 3 > d 0 ) and let The following proposition provides some basic local stability results for Model 1.
It is easy to see that lim . We now only need to consider the stability of (x 1 , x 2 , x 3 ). Routine local stability analysis leads to that the stability of the equilibria depends on eigenvalues such that It is straightforward to conclude the linear stability for the three different cases.
Proposition 3.2 implies that there is no oscillatory behavior on the BCR-ABL/ABL% (see (6.1)), which suggests that the oscillatory nature of some individual patients' data may be caused by stochastic factors not considered here. Also notice that E 1 corresponds to 0% in BCR-ABL/BCR(%), while BCR-ABL/BCR(%) does not apply to E 0 .

Model 2: A Simple Density Dependent Model
The second model, based on a model by Michor [12,11], describes the change in abundances of normal stem cells x and leukemia stem cells y respectively: ] r x Φ, r y φ represent the density dependent cell dividing rates, and c x , c y are parameters that simulate the crowding effect that is seen in the bone marrow microenvironment. The normal and leukemia stem divide at rates at most r x , r y , respectively, per day. The death rate of both normal and leukemia stem cells is represented by d 0 . This model assumes that cells can reproduce both symmetrically and asymmetrically.

Basic Analysis of Model 2
Our first proposition present the positivity and bondedness results for Model 2.
Proof. We show first that x(t) > 0 and y(t) > 0 when exist for t > 0. If not, there is a first time t 1 > 0 such that x(t 1 ) = 0 or y(t 1 ) = 0. Assume first that x(t 1 ) = 0. Then for t ∈ [0, t 1 ], we see that x (t) ≥ −d 0 x(t) and hence x(t 1 ) ≥ x(0)e −d 0 t 1 > 0, a contradiction. Similar contradiction can be obtained by assuming that y(t 1 ) = 0, proving the positivity of the solutions.
Next to establish the boundedness of solutions. Observe that By a comparison argument, we can conclude that x is bounded by max{ 1 d 0 cx (r x − d 0 ), x(0)} and y is bounded by max{ 1 d 0 cy (r y − d 0 ), y(0)}. The next proposition present stability results for Model 2.
, and no interior equilibrium. There are no periodic solutions of (4.1) and (4.2).
(1) If r x < d 0 and r y < d 0 , then E 0 is the unique equilibrium and E 0 is (globally) stable.
(2) If r x > d 0 and r y < d 0 , we have equilibria E 0 and E 2 , while E 0 is unstable and E 2 is (globally) stable. (3) If r x < d 0 and r y > d 0 , we have equilibria E 0 and E 1 , while E 0 is unstable and E 1 is (globally) stable. (4) If r x > d 0 and r y > d 0 , we have all three equilibria: E 0 , E 1 , and E 2 .
, then E 1 is unstable and E 2 is (globally) stable.
Proof. All the equilibria can be easily calculated. Since there is no interior equilibrium, and the solutions are bounded, then by Poincare-Bendixson Theorem and the fact that a periodic orbit must enclose at least one equilibrium, there are no periodic solutions of (4.1) and (4.2).
Routine local stability analysis leads to that the stability of the equilibria depends on eigenvalues for E 0 , E 1 , and E 2 respectively. Then it is straightforward to conclude the linear stability for the four different cases. Since we have eliminated the existence of periodic solutions, local stability implies global stability for case (1), (2), and (3).
Proposition 5.2 implies that there is no oscillatory behavior on the BCR-ABL/ABL%, again suggests that the oscillatory nature of some individual patients' data may be caused by stochastic factors not considered here. Also notice that E 1 corresponds to 100% in BCR-ABL/BCR(%) and E 2 corresponds to 0% in BCR-ABL/BCR(%), while BCR-ABL/BCR(%) does not apply to E 0 .

Data
We used data from a previous study [14,17] that consists of samples from patients who were recruited in Germany between June 2000 and January 2001 and enrolled in the International Randomized Study of Interferon and STI571 (IRIS study). Müller et al. [14] studied 139 patients, who were recently diagnosed BCR-ABL positive chronic phase CML patients. Out of these patients, 69 were treated with imatinib and 70 were treated with interferon (IFN)/Ara-C. Our analysis only considers the 69 patients who were treated with imatinib. These patients received 400mg orally daily. The blood samples were collected, either by mail or locally, after months 1, 2, and 3, and then were collected at three month intervals. The data consists of BCR-ABL/ABL% from months ranging from 0 to 66 from each patient. The BCR-ABL transcripts were obtained using qualitative reverse transcriptasepolymerase chain reaction and ABL was used as the control gene. For more information, see [14].
To compare the clinical data to Model 1, we used the following to approximate the percents: To compare the clinical data to Model 2, we used the following to approximate the percents: In using these approximations, we assume that a BCR-ABL positive cell also contains a non mutated chromosome 9 and 22 and thus the BCR-ABL/ABL% values cannot be over 100. Therefore we did not analyze any patients with BCR-ABL/ABL values over 100% and only analyzed data from the remaining 51 patients. Since the data ranges from 0 to 66 months and after month 3, the samples were collected every 3 months, ideally each patient should have 25 data points. Out of these 51 patients, 11 patients had fewer than 10 data points. Since we are comparing the two models to the data, we only considered the 40 patients with more than 10 data points. We assume that, for each patient, the initial BCR-ABL/ABL% value consists of 99% BCR-ABL dependent cells and 1% BCR-ABL independent cells.

7.
A comparison of the two models 7.1. Simulations. To compare the two models we ran simulations with MATLAB using the clinical data of the 40 patients from the earlier study [14]. We used the MATLAB built-in function fminsearch to find the optimum parameters for each model for each patient. We calculated the error using the following equation: where N represents the total number of data points, y i represents the actual value, andŷ i represents the estimated value from the models.
Parameter Meaning r 1 Maximum proliferation rate of BCR-ABL D population r 2 Maximum proliferation rate of BCR-ABL I population r 3 Maximum proliferation rate of Normal population p 3 Parameter that simulates the crowding effect n Hill coefficient q 1 Minimum BCR-ABL D cell quota q 2 Minimum BCR-ABL I cell quota k 1 Maximum BCR-ABL D to BCR-ABL I mutation rate k 2 Maximum BCR-ABL I to BCR-ABL D mutation rate K 1 BCR-ABL D to BCR-ABL I mutation half-saturation level K 2 BCR-ABL I to BCR-ABL D mutation half-saturation level q m1 Maximum BCR-ABL D cell quota v m Maximum cell quota uptake rate v h Uptake rate half-saturation level b Cell quota degradation rate µ m Rate at which BCR-ABL is used within the cell for growth p Average intracellular production rate of BCR-ABL d 1 Average intracellular degradation rate of BCR-ABL u Rate of BCR-ABL inhibition by imatinib treatment After comparing the errors for each patient from each of the models, 26 out of 40 patients had a smaller error associated with Model 1 compared to Model 2. Table 1 contains statistical information about the errors of the two models. The median error for Model 1 was 0.9939 whereas the median error for Model 2 was 1.080. The average error for Model 1 was 2.030 whereas the average error for Model 2 was 2.181. When comparing the difference between the errors for each model for each patient, only 4 patients out of the 40 patients had a difference in error that was greater than 1. Figure 1 contains the simulations for three patients whose Model 1 error was smaller than the Model 2 error. Figure 2 contains the simulations for three patients whose Model 2 error was smaller than the Model 1 error. Figure 3 contains the simulations for three patients where the two models were similar in terms of error.  Tables 3 and 5 contain statistical information about the parameters for the 40 patients for Model 1 and Model 2 respectively. We used the fminsearch function in Matlab to find optimal model parameters with respect to the error defined in (7.1). The initial guesses for the parameters of the models were initially fit by hand to provide good qualitative agreement with the clinical CML data. Note that both models use the stem cell death rate, d 0 = 0.003/day [13]. We can see that, although Model 2 has fewer parameters, the range of the values is extremely large and biologically unrealistic. The maximum value for r 3 in Model 1 is about 0.015 per day, whereas the maximum dividing rate for normal stem cells in Model  Parameter that simulates the crowding effect Table 4.

Model 2 Parameter Meanings
In Model 1, the growth rate of the nonresistant leukemic stem cells has a maximum value of about 0.021 per day, which is relatively close to the value of 0.008 per day used by Michor [13]. However, the maximum value of r y in Model 2 is about 4.392 × 10 7 per day, which is biologically unrealistic. For Model 1, the median and average values for r 1 are 0.0036 and 0.0054 respectively, while the median and average values for r y in Model 2 are 0.011 and 1.738 × 10 6 respectively. Although Model 2 seems to be a much simpler model and fit the data similarly to Model 1, we can see that the parameter ranges are biologically unrealistic, suggesting that Model 1 to be biologically more plausible.  Table 5. Model 2 Parameter Statistics most likely due to resistance to Imatinib. Although both models were able to show resistance for at least one patient, the resistance described by Model 1 seems more biologically relevant. The BCR-ABL dependent population in Model 1 represents the non-resistant cells while the BCR-ABL independent population represents the resistant cells. Model 1 suggests that a relapse occurs when the BCR-ABL dependent population is replaced by the BCR-ABL independent population. Abbott and Michor describe a slightly more complex model [1], which also describes resistance, however the model was not available nor completely described in the paper, so we were not able to compare Model 1 to their model with resistance. Models that describe resistance are important biologically since resistance is a common problem in cancer treatments. A model by Foo et. al predicted that, for every 100 patients treated with only imatinib, 89 will eventually develop resistance [4]. Models can be used to determine when the patient will stop responding to treatment based on their previous data. For chronic phase patients who start imatinib treatment early, only 12% develop resistance within the first two years of treatment [13]. This implies that resistance is probably not an immediate occurrence.
We searched for signatures of resistance using the parameters estimated from patient data. Although clinical data only contains values up to about 5.5 years, we ran the simulations for a longer time span to compare the ability of Model 1 and Model 2 to predict long-term resistance. The simulation for patient 1 in Model 2 showed an increase in BCR-ABL/ABL% around day 1000 (Figure 4), suggesting that the CML cell population started to outgrow the normal cell population. This increase in the leukemia cell population was then verified by the simulation using Model 2 for the proportion of the cell populations ( Figure 4). However, patient 1 had the largest error out of all of the other patients for both models. The error for patient 1 from Model 1 was 18.14 and the error from Model 2 was 17.06. The next largest error out of all the patients for Model 1 was 6.858 and for Model 2 was 7.533, which are relatively small errors. (Errors for patient 20) Thus, neither model accurately describes the patient 1 data, so it seems irrelevant that Model 2 describes resistance for patient 1. Figure 4 also contains graphs where resistance was predicted by Model 1 by an increase in BCR-ABL/ABL%, an increase in the proportion of CML cells, and a decrease in the proportion of normal cells. The errors for these patients were relatively smaller than the errors for patient 1.

Discussion
The two models compared in this paper both describe the treatment of chronic myeloid leukemia but do so in different ways. Model 2 describes the competition of leukemic and normal stem cells. In Model 1, normal stem cells are in competition with leukemic cells and the growth of leukemic cells depends on the concentration of BCR-ABL. Model 1 also incorporates more biological detail than Model 2 by describing the subcellular dynamics of BCR-ABL and allowing for phenotypic switching between BCR-ABL dependent and independent leukemic cell populations. We compared these two models using clinical data and simulations in order to gain insights into how adding these biological details can more accurately describe and predict resistance to imatinib treatment for CML patients. We found that although Model 2 is a much simpler model, it still describes the data well for some patients. However, the parameter ranges for Model 2 are extremely large and biologically unrealistic. In contrast, Model 1 fit the clinical data better for more patients (26/40) and the estimated parameter ranges were more realistic. This result suggests that the additional biological mechanisms described in Model 1 are relevant, since they can increase the accuracy in data fitting in a majority of patient data.
A mathematical model that predicts treatment resistance in cancer can be a valuable tool to increase the effectiveness of treatment strategies. This is especially relevant for CML where 62% of accelerated phase patients treated with imatinib develop resistance within 2 years of treatment [2]. Using simulations with parameter estimates from patient data, Model 1 was able to predict resistance to imatinib in terms of both an increase in BCR-ABL/ABL% and in the proportion of leukemic cells in the total stem cell population. Model 2 predicted resistance in a single patient, however, the error associated with Model 2 for this patient was relatively large (2.265 times greater than any other patient). Since Model 2 is unable to accurately fit the clinical data for this patient, it is unlikely that the consequent prediction of resistance is accurate. Thus, Model 1 was able to show resistance in some patients in a biologically meaningful way whereas Model 2 was not able to show resistance and accurately model the clinical data.
The results we have discussed for Model 1, although promising, are mainly computational and in need of further exploration. A thorough mathematical analysis of Model 1 can provide additional insights into how the subcellular regulation of BCR-ABL levels dictates the long-term transition of CML cells to an imatinib resistant phenotype, i.e. BCR-ABL independent. In future computational work we can further evaluate the accuracy of Model 1 to predict resistance by using patient data that exhibits long-term (i.e. > 2 years) resistance to imatinib. For example, the simulations for Model 1 where resistance occurs suggest that a more optimal patient data set for A novel mechanism encoded in Model 1 is the BCR-ABL dependent switching between BCR-ABL dependent and independent populations. We speculate that these transitions could have an epigenetic basis. This is in contrast to previous CML models that have only considered transitions due to genetic mutations [1] or switching between a proliferative and nonproliferative state [17,10]. Indeed, recent studies have elucidated important epigenetic changes that may cause resistance to the imatinib drug. For example, imatinib therapy could cause drug resistance by affecting epigenetic alterations in cells that down-regulates tumor suppressor genes [15]. Such alterations could lead to a reduced dependence of leukemic cells on BCR-ABL to express a malignant phenotype, i.e. a BCR-ABL independent cell population. Another study showed that aberrant changes in DNA methylation could be an epigenetic marker associated with imatinib resistance [8]. Our computational work here highlights the importance of further experimental work to ascertain the rate at which epigenetic transitions occur in CML, how this rate is related to imatinib dosage, and how it affects imatinib resistance.
Kareva, Berezovskaya, and Castillo-Chavez [9] analyze the balance between immature and mature myeloid cells and how this balance effects tumors. They claim that if there is a small enough population of cancer cells, then there is a small region of initial conditions where the immune system alone will cure the cancer and the patient will not need treatment. Future work could look into incorporating the mature and immature myeloid cell populations into the normal cell population in Model 1 as well as combining intermittent imatinib therapy with the immune system's defense. The work can also be expanded in the future by using more biologically relevant function forms of p(x) and also comparing Model 1 to the model by Roeder et al. [17].

Acknowledgments
We thank Dr. Ingo Roeder for providing the data and the referee for many helpful suggestions. The work of Rebecca A. Everett, Yuqin Zhao and Yang Kuang is supported in part by NSF DMS-0920744.