Model-based decision rules reduce the risk of molecular relapse after cessation of tyrosine kinase inhibitor therapy in chronic myeloid leukemia

Chronic myeloid leukemia (CML) is a clonal disorder that is cytogenetically characterized by a translocation of chromosomes 9 and 22, resulting in the formation of the BCR-ABL fusion gene on the level of hematopoietic stem cells (HSCs).^1,2  The oncogenic capacity of this gene, which is located on the shortened chromosome 22 (Philadelphia [Ph] chromosome) is well-established. Its product, the BCR-ABL protein, is a constitutively activated tyrosine kinase, which has been shown to be responsible for the pathogenesis of CML.² 

The current therapy of choice for de novo CML is oral administration of tyrosine kinase inhibitors (TKIs), such as imatinib mesylate (IM),³  or recently approved second-generation TKIs dasatinib and nilotinib.^4,5  Main mechanism of action of IM (and other TKIs) is inhibition of the oncogenic BCR-ABL tyrosine kinase,⁶  resulting in the switching-off of downstream signaling pathways promoting leukemogenesis.³  It is well-established that IM selectively acts on leukemia cells inducing a proliferation inhibitory effect.⁷  Furthermore, it has been reported that the apoptotic rate of actively proliferating leukemia cells is increased under IM therapy.^8,9  Although treatment with IM could still not be demonstrated to be curative, it is suited to achieve a sustained control of the disease in the majority of patients.

Molecular monitoring of tumor load in peripheral blood using quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) revealed that in most patients IM monotherapy induces a biphasic decline of BCR-ABL transcript levels, characterized by an initially steep decline, followed by a second moderate decline.^10,11  A sensible explanation for this behavior is the rapid initial depletion of actively cycling BCR-ABL–positive cells, followed by the slow elimination of residual leukemic stem cells (LSCs)¹¹  because of their low turnover.^12,13  Alternatively, this behavior has been interpreted as the result of distinctive IM effects on different hematopoietic cell stages including IM insensitivity of stem cells.¹⁰ 

Although molecular IM response in the population of CML patients can be approximated by an average biphasic decline, a large heterogeneity in response patterns between individuals is observed.¹¹  In the absence of IM resistance, long-term treatment could be shown to be associated with sustained decline of BCR-ABL transcript levels.¹⁴  These levels potentially fall below the detection threshold of PCR techniques, which allow for the detection of 1 CML cell in 10⁵ nucleated blood cells.^15,16 

Side effects of IM therapy as well as economic considerations pose the question whether IM can be safely discontinued after achieving deep molecular remissions (eg, MR^5.0). Uncertainty arises from the fact that on discontinuation of IM treatment, a heterogeneous picture is obtained: although some patients retain previously achieved molecular responses,^17,18  a molecular recurrence of BCR-ABL transcripts is observed in others.¹⁸  Relapses can also be observed in patients lacking any measurable BCR-ABL transcripts in peripheral blood.¹⁷ 

We sought to predict patient-specific long-term time courses of CML under IM treatment and to support decision-making for potential treatment cessation. We could demonstrate that a statistical description of the disease kinetics is not sufficient to correctly estimate numbers of residual LSCs, which are a critical determinant of relapse after therapy discontinuation. Therefore, we adapted our established mathematical model of HSC and leukemia organization,^11,19,20  which was demonstrated to consistently explain CML genesis and IM treatment on the population level, to predict patient-specific stem cell numbers and long-term treatment outcome. Based on 7-year follow-up data of BCR-ABL transcript dynamics from the German cohort of the IRIS trial,¹⁴  we determined model parameters that quantitatively characterize the inter-individual heterogeneity of molecular treatment response. Given a patient's BCR-ABL transcript kinetics, the adapted model generates predictions for patient-specific long-term response to IM as well as individual times to complete eradication of minimal residual disease (MRD). To test the validity of our approach, we used an independent dataset from the CML IV trial.²¹  Furthermore, we derived a model-based predictor for the individual risk of molecular relapse on treatment cessation, which is complemented by a model-independent approximation that can easily be calculated from clinical data. Using a comparative simulation study, we showed that the proposed predictor results in a superior clinical rule to decide on potential discontinuation of therapy compared with relying on a fixed (eg, 2 years) time in sustained deep molecular remission.

The mathematical model underlying our results has originally been developed for the HSC system.¹⁹  We previously demonstrated that the model is well-suited to explain clonal dynamics in the human situation.^11,20  In brief, the model assumes CML to be caused by a competition process between malignant (Ph⁺) and normal (Ph⁻) hematopoietic cells, with quantitative differences in specific model parameters. IM therapy is assumed to induce an apoptotic effect and inhibition of the proliferative activity of Ph⁺ stem cells. The stochastic, single cell-based model is implemented in the programming language C++. The source code can be obtained from the authors. A schematic representation and a detailed technical description of the model implementation can be found in supplemental Methods (available on the Blood Web site; see the Supplemental Materials link at the top of the online article).

As a training set for model development we used the German cohort of the IRIS trial (n = 69; 400 mg IM daily).¹⁴  Each individual was represented by a time series of BCR-ABL/ABL measurements. Six patients (9%) were excluded from the analysis because of insufficient data (defined as less than 5 measurements). Seven patients (10%) suffered from considerable molecular relapse under continued therapy because of IM resistance, whereas 2 (3%) showed treatment failure for unknown reasons (defined as no decline of BCR-ABL levels within first year). In 3 individuals (4%) a uniphasic decline was observed. The remaining 51 (74%) patients, who are characterized by a biphasic decline of BCR-ABL transcript levels, were used for our modeling analyses.

As independent validation data, we used the 400 mg IM arm of the CML IV trial (n = 280).²¹  Because of shorter overall follow-up compared with IRIS, 95 patients (34%) were excluded for reasons of insufficient data. Another 56 (20%) were associated with molecular relapse under IM, whereas 3 (1%) showed treatment failure. In 14 patients (5%) a uniphasic decline was observed. The remaining 112 (40%) showed a biphasic decline of BCR-ABL levels. However, only 31 patients had a minimum of 14 qRT-PCR measurements, required for reliable slope estimation (see “Statistical analysis of clinical data”). Hence, only this subset was used for validation.

The selection process described was necessary because our current model version is not able to quantitatively account for treatment failure (eg, because of IM resistance) and uniphasic decline kinetics. Still, the presented model analysis applies to more than 90% of patients who respond to IM without observed relapse. In addition, we confirmed (data not shown) that exclusion of the minority of patients showing a uniphasic decline does not bias the median behavior within the cohorts.

For this analysis, all BCR-ABL/ABL ratios were standardized according to the International Scale (IS). Each qRT-PCR measurement was preceded by a plasmid dilution series, that is, the smallest absolute number of BCR-ABL transcripts that could be detected in an individual assay was determined. In case of BCR-ABL negativity by qRT-PCR, nested PCR was performed. For negative results, BCR-ABL/ABL was assumed to be zero. For positive results, BCR-ABL/ABL was approximated using the smallest positive BCR-ABL copy number from the dilution series (representing an upper bound for the unknown true value).

On a logarithmic scale, the previously mentioned biphasic decline kinetics of BCR-ABL levels in response to IM are sufficiently described by 2 piecewise linear relationships.^10,11  Herein, the slope of the first, steep decline is denoted as α and the slope of the second, moderate decline as β. Patient-specific values for α and β were calculated using a segmented linear regression model (library “segmented,“ statistical programming environment R). Specifically, we used the Davies test to check the null hypothesis of no slope difference.²²  If the null hypothesis was rejected at the 5% significance level, the breakpoint between first and second slope was estimated by an iterative procedure. As starting value of the iteration, we chose 250 days after therapy initiation based on previous knowledge. For quantification of correlations among slope values and breakpoint, Pearson correlation coefficient was applied.

First, we determined the heterogeneity within the patient population of the German IRIS cohort with respect to the biphasic decline of BCR-ABL transcript levels (Figure 1). Because of the association of the second decline with long-term treatment success, it would be desirable to predict this decline from the early response. However, we did not find a correlation between strength or duration of the first decline kinetics (characterized by slope and breakpoint, respectively) and the secondary decline (Figure 2). Hence, prediction of long-term outcome is not possible based on early decline kinetics alone, but requires a sufficient number of measurements of the secondary slope. Interestingly, we found a significant correlation between first decline and breakpoint (supplemental Figure 3 and supplemental Results).

Second, we explored how the estimated kinetic parameters, that is, first (α) and second slope (β) of BCR-ABL transcript levels, translate into the parameter space of the dynamic model. Specifically, we used 2 model parameters (degradation rate r_deg and transition characteristic f_ω, which were previously found to sensitively affect IM therapy effects; see supplemental Methods for details) as patient-specific parameters to account for the heterogeneity in IM response patterns.^11,20  Other model parameters remained unchanged with respect to the values estimated for IM-treated CML (supplemental Table 1). Using an iterative algorithm we calculated model parameters r_deg and f_ω for each patient based on the observed BCR-ABL decline kinetics (see supplemental Results).

Based on patient-specific parameters of the dynamic model we generated predictions of long-term BCR-ABL transcript kinetics for individual patients. Using the example patient from Figure 1A, we demonstrate the application of the procedure: the patient is characterized by slope parameters α = −4.08/year and β = −0.24/year, from which we calculated r_deg = 2.50 and f_ω = 0.055, the parameters of the dynamic model that optimally account for the particular BCR-ABL decline kinetics. Using these parameter estimates, we simulated this particular patient in silico, both for the time of available follow-up data (Figure 3A) and, most importantly, thereafter (Figure 3B). The model predicts that after approximately 15 years the BCR-ABL transcript levels of this particular patient permanently fall below the detection threshold of qRT-PCR (Figure 3B), which we assume to be MR^5.0, that is, reduction of BCR-ABL transcript levels below 0.001%.²³  After approximately 25 years it can be expected that no BCR-ABL–positive cells will be detected in peripheral blood at all, irrespective of the accuracy of the applied measuring method. This is because after this time only very few (mostly quiescent) residual LSCs, which are located in the bone marrow, remain in the system (Figure 3C). Note that the (absolute) stem cell number is not accessible in the clinic, but is a result of the in silico model. For the example patient the model predicts complete eradication of CML cells after 29.4 years of continuing IM therapy, with a standard deviation of 0.78 years.

Table 1 summarizes the model predictions for the considered patients. Even under the ideal circumstances assumed in our model, that is, uninterrupted IM administration and permanent absence of both IM resistance and disease-unrelated complications, the majority of individuals are predicted to require several decades of IM therapy before MRD is completely eradicated. The cumulative rates of complete eradication after 15 and 30 years of treatment are estimated to be 14% and 31%, respectively. In approximately 67% of the patients, residual leukemic cells are predicted throughout the remaining lifetime, assuming life expectancy of 80 years.

As this procedure relies on the full 7-year follow-up data, we examined whether shorter observation times are sufficient for generation of adequate model predictions, that is, whether actual therapeutic predictions can be made after less than 7 years (supplemental Figure 7). We found that in this patient cohort, for obtaining reliable predictions, at least 14 PCR measurements (6-8 per slope) are required, irrespective of follow-up duration. Consequently, it can be speculated that the accuracy of the predictions for fixed observation times can be enhanced if PCR measurements are performed more frequently (see supple-mental Results).

To validate the described method of calculating parameters of the dynamic model from clinically determined BCR-ABL kinetics, we analyzed a second, independent dataset (CML IV trial), which is comparable with the IRIS data with respect to patient characteristics, IM dosage, and response patterns (for selection criteria see “Methods”). We can show that the algorithm derived based on the IRIS data can be applied without changes to the independent dataset. The model results consistently account for the median behavior within the patient population as well as for interpatient heterogeneity of the CML IV cohort (see supplemental Results and supplemental Figures 9 and 10). In addition, we generated model predictions on long-term IM therapy. For the CML IV cohort, the results are summarized in Table 1. The cumulative cure rates after 15 and 30 years of IM treatment (16% and 42%, respectively) are statistically indistinguishable from the predictions generated for the training dataset from the IRIS trial (Table 1).

For the 51 patients from the IRIS trial we simulated the situation of potential IM cessation. In this scenario, each patient was treated with IM in silico until BCR-ABL transcript levels fell below the MR^5.0 detection threshold for 2 consecutive years (assuming trimonthly measurements). For those 42 patients who reached at least a sustained 2-year MR^5.0 within 20 years or less of simulated IM treatment we generated model predictions for the time to molecular relapse after therapy discontinuation (solid line in Figure 4). It is predicted that the majority of molecular relapses is observed within the first 6 months with only few relapses occurring more than 1 year after treatment cessation. Twelve of 42 patients are predicted not to relapse at all. This model prediction quantitatively resembles recently published clinical results (Figure 4).^17,18  Please note that stopping rules and measurement protocols are not exactly comparable: whereas the simulation and the Australian trial¹⁷  only require (at least) 2 years in MR^5.0 before IM discontinuation, the French trial¹⁸  also requires a total IM treatment time of at least 3 years. Furthermore, in our simulation and in the Australian trial molecular relapse is defined as detectable BCR-ABL transcripts at any level in 2 consecutive PCR measurements, whereas the French trial also requires the BCR-ABL to ABL ratio to be at least 0.001%.

Model-predicted relapsing and nonrelapsing subgroups are depicted in Figures 5A and B with respect to IM response kinetics. Qualitatively, a relatively steep second decline is necessary but not sufficient for sustained MR^5.0 after therapy cessation. In contrast and contrary to intuition, a relatively steep first decline is adverse to achieving this goal. This is because of the fact that a steep first decline is associated with rapid depletion of leukemic cells in peripheral blood (where PCR measurements are performed) but not necessarily in the bone marrow. As a consequence, the treatment time required to substantially reduce MRD is potentially underestimated in the clinic. We found a linear relation (log scale) between the ratio of first and second slopes and the number of residual LSCs at the moment of treatment cessation (Figure 5C). In contrast to relapsing patients (circles), nonrelapsing patients (triangles) are characterized by log₁₀(number of LSCs) ≤ 2.3 (horizontal line). This characteristic, which provides an optimal discrimination of the 2 groups, is based on the prediction of the dynamic model. A good approximation of this classification (approximately 5% misclassifications) can be obtained if the criterion α/β ≤ 16 (vertical line) is applied. We call the latter predictor “model-independent” as it does not require model simulations to be calculated. In contrast to the “model-based” predictor that uses the model-predicted number of residual LSCs, it can directly be calculated from clinically measureable variables.

The “model-independent” predictor, which relates to the condition of 2-year MR^5.0 before treatment cessation, can be generalized. Let t describe the duration of stable MR^5.0 under IM treatment until therapy discontinuation (given in years). Then a patient is predicted to remain in MR^5.0 after treatment cessation if the relation α/β ≤ 8 × t holds. That is, the molecularly relapsing subgroup is characterized by a first decline that is more than 8 × t greater than the second decline. This implies that shorter MR^5.0 times require steeper second declines (given fixed first declines) to predict sustained MR^5.0 after therapy discontinuation. The predictor based on the number of residual LSCs can be applied in any case, irrespective of MR^5.0 duration.

Solving the inequality for t, the generalized predictor results in a simple decision rule given as t ≥ α/(8 × β). That means the time of sustained MR^5.0, which is predicted to ensure continuing MR^5.0 after treatment cessation, can be estimated for each patient individually. Applying this rule to the in silico IRIS cohort, that is, performing computer simulations for each patient and discontinuing IM after the calculated patient-specific time in MR^5.0, a significantly lower relapse rate (P < .001, log-rank test) with approximately 3 times less relapses at 2 years after IM cessation can be observed (solid line in Figure 6A) compared with the fixed time of 2 years in MR^5.0 (dashed line). Median treatment time in MR^5.0 before IM cessation resulting from this model prediction is estimated to be 2.8 (range: 0.3-20.5) years (Figure 6B). We also considered a scenario where IM is discontinued after a fixed treatment time of 2.8 years in MR^5.0, which equals the median treatment time in the individualized stopping regimen. This result, which is depicted by the dotted line in Figure 6A, still represents an inferior outcome with respect to the proposed patient-specific stopping rule.

The procedure described here allows for a better appraisal of the individual prognosis of CML patients with direct impact on clinical decision-making. In the clinics, the tumor burden cannot be followed beyond the time point when BCR-ABL transcript levels have fallen below the detection threshold of qRT-PCR. In the majority of patients, this threshold is reached only a few years after initiation of IM therapy. However, LSCs, which are believed to drive the disease,²⁴  reside in the bone marrow and can only be detected in peripheral blood samples if they contribute to blood production in sufficient amounts. There is increasing evidence that primitive HSCs, and probably also LSCs, can “hide” from detection in peripheral blood because of their sustained quiescence.^12,25,26  As the proposed model acts on the level of stem cells, particularly accounting for the duality of stem cell quiescence and active proliferation, it offers a valuable tool to quantitatively complement clinically accessible measures.

With respect to the statistical analysis of the data, taken from the German IRIS cohort,¹⁴  we were able to confirm and extend previous reports of a typical biphasic decline of BCR-ABL transcript levels under IM,^10,11  where the first decline probably characterizes initial response to IM and the second decline potentially represents elimination of residual LSCs.^27,28  Furthermore, we observed a strong correlation between steepness of the first decline and breakpoint separating both declines. This notwithstanding, we found that steepness of the first decline as well as the breakpoint are not sufficient to predict the second slope and thus the time point of MRD eradication. That is, to evaluate the possibility of IM cessation, steepness of the second decline needs to be assessed with sufficient accuracy. Herein, the quality of the estimate sensitively depends on “late” data points, especially in case of few available measurements. Note, however, that follow-up time itself is not the most important parameter in this respect. More important is the number of measurements (6-8 within each of the 2 slope phases). As the BCR-ABL transcript levels of the IRIS patients have, in median, been assessed 7 times within the first year of IM therapy but only twice within each year thereafter, the critical number of measurements is reached only after 4 to 5 years. Performing, for example, 2 monthly measurements throughout, 2 to 3 years of follow-up might be sufficient for achieving the same quality of predictions.

The proposed method cannot be applied in case of IM resistance.²⁹  Furthermore, there are few cases of biphasically declining kinetics where the application of the procedure fails. This applies to “bad responders,“ that is, patients characterized by a comparatively flat first decline, and consequently by an unusually late breakpoint (after more than 3 years of IM therapy). Still, our model is a valuable tool to estimate the time to complete eradication of the leukemic clone in the majority of CML patients. Model results are consistent with the clinical experience that the majority of patients need to be treated with IM for a very long time period, before (if at all) complete leukemia eradication can be achieved. However, it should be kept in mind that complete eradication might not be necessary before IM therapy can be safely discontinued. Small BCR-ABL levels can be detected in healthy volunteers³⁰  leading to speculations that there exists a critical threshold below which the leukemic clone is small enough to be intrinsically controlled. It is one of our major results that the proposed model is able to identify those patients for whom treatment discontinuation would be a safe option. Along these lines, the model allows to predict the individual outcome in case of treatment stop. We show that patients with no predicted relapse are associated with a relatively steep second decline resulting in low numbers (≤ 200), or complete eradication, of residual LSCs at the moment of discontinuation. Predicted “late” relapses (more than 1 year after IM cessation) are also characterized by relatively small numbers of LSCs. However, in these cases the leukemic clone is still large enough to slowly expand and repopulate the marrow.

In addition to the individual model predictions of long-term BCR-ABL levels in peripheral blood as well as of residual LSC numbers in the bone marrow, we propose a simple “model-independent” predictor that allows for discrimination between relapsing and nonrelapsing patients after treatment cessation based on their individual treatment response kinetics. To investigate the quality of the predictor we performed simulation studies comparing an individualized treatment discontinuation strategy versus fixed time periods in sustained deep molecular remission before treatment cessation and found the former to be superior with respect to observed molecular relapse rates after discontinuation. It should be noted that this gain in the nonrelapsing fraction comes at the cost of increasing the median time to IM cessation. However, increasing the time of IM treatment in deep molecular remission uniformly for all patients does not result in comparable benefit because applying this strategy, according to our model analyses, some patients are treated unnecessarily, whereas others require longer TKI therapy to avoid relapsing on treatment discontinuation.

Our model prediction, regarding the strategy that applies a fixed 2-year period in MR^5.0 before IM discontinuation, closely resembles the outcome from 2 independent clinical trials^17,18  with respect to time to molecular relapse. Although stopping rules and measurement protocols are not identical, it is interesting to note that there are only minor differences between the individual Kaplan-Meier estimates, which we interpret as strong evidence for the validity of our approach. Comparing the 3 scenarios, the highest percentage of events can be observed in the in silico cohort, which is probably because of the fact that the “real” patient cohorts are naturally selected more strongly for good responders, in the sense of patients who reach a sustained deep remission comparatively early. In contrast, the simulation also includes patients reaching the 2-year MR^5.0 goal relatively late (up to 20 years after initiation of IM therapy), which might increase the relapse probability after therapy discontinuation slightly.

The validity of our model predictions sensitively depends on the underlying assumptions, such as the mechanisms of IM action. However, the predicative power of our model with respect to relapses in discontinuation trials (Figure 4) as well as in previous settings^11,20  ensures its applicability. Other mathematical CML models did not assume a direct or indirect IM effect on LSCs.¹⁰  However, these assumptions have been modified over time as they could not account for long-term declines of BCR-ABL levels in IM treated patients. More recent models^31,32  by those authors can be expected to yield qualitatively similar results to our model as they also assume clonal competition between HSCs and LSCs and a direct effect of IM on LSCs.

As our results were obtained by computer simulations, we recommend testing our strategy in a prospective randomized clinical trial that compares the described discontinuation strategies. The described approach is not limited to IM but can also be applied in an analogous way to first-line second generation TKIs, for example, dasatinib or nilotinib.^4,5  The mechanisms of action of these drugs are similar to IM,³³  and molecular response is also characterized by a biphasic behavior.²⁷  Currently, decline parameters cannot yet be determined because of insufficient follow-up, but we plan to adapt our approach in due course.

This work was supported by the Bundesministerium für Bildung und Forschung (BMBF) grant on Medical Systems Biology “HaematoSys” (BMBF-FKZ 0315452) and grant RO3500/1-2 by the German Research Foundation (DFG). Novartis supported the PCR analysis of the IRIS samples. The CML IV study was supported by Deutsche Krebshilfe, Novartis, BMBF, Deutsche José-Carreras-Leukämiestiftung, Roche, and Essex Pharma.

Contribution: M.H., I.R., and M.L. conceived the mathematical model; M.H. conducted the simulations and the model analyses; I.R. and M.L. provided supervision; M.M., A.H., and R.H. provided the clinical data; I.G., I.R., and M.L. provided advice and input for modeling and data analysis; M.H. drafted the paper; I.G., M.M., R.H., A.H., M.L., and I.R. critically reviewed the paper; and all authors read and approved the final version of the paper.

Conflict-of-interest disclosure: A.H. and M.M. receive research support from Novartis, Bristol-Myers Squibb, Pfizer, and Ariad. The remaining authors declare no competing financial interests.

Correspondence: Matthias Horn, Institute for Medical Informatics, Statistics and Epidemiology, University of Leipzig, Haertelstrasse 16-18, 04107 Leipzig, Germany; e-mail: matthias.horn@imise.uni-leipzig.de.