Quantitative prediction of long-term molecular response in TKI-treated CML – Lessons from an imatinib versus dasatinib comparison

Longitudinal monitoring of BCR-ABL transcript levels in peripheral blood of CML patients treated with tyrosine kinase inhibitors (TKI) revealed a typical biphasic response. Although second generation TKIs like dasatinib proved more efficient in achieving molecular remission compared to first generation TKI imatinib, it is unclear how individual responses differ between the drugs and whether mechanisms of drug action can be deduced from the dynamic data. We use time courses from the DASISION trial to address statistical differences in the dynamic response between first line imatinib vs. dasatinib treatment cohorts and we analyze differences between the cohorts by fitting an established mathematical model of functional CML treatment to individual time courses. On average, dasatinib-treated patients show a steeper initial response, while the long-term response only marginally differed between the treatments. Supplementing each patient time course with a corresponding confidence region, we illustrate the consequences of the uncertainty estimate for the underlying mechanisms of CML remission. Our model suggests that the observed BCR-ABL dynamics may result from different, underlying stem cell dynamics. These results illustrate that the perception and description of CML treatment response as a dynamic process on the level of individual patients is a prerequisite for reliable patient-specific response predictions and treatment optimizations.

Bi-exponential parameter estimation. We applied a non-linear least-square approach to fit a biexponential function LRATIO(t) to the BCR-ABL transcript measurements up to time t for an individual patient i, thereby obtaining a tuple of estimated statistical parameters " # $ = ( , , , ) along with an estimated residual variance " # . . The same approach can also be applied to simulated (model-generated) BCR-ABL time courses with model parameters = 0f 2 345 , f 6 345 , r 89: , LRATIO ABAC , r CDEBF H, thereby obtaining " I $ . The approach accounts for left-censored observations, i.e. (BCR-ABL-negative) measurements that are below the quantification limit (QL) of the RT-PCR. For each RT-PCR measurement we obtain an individual estimate of its QL, based on the available ABL reads (QL = 3/ABL * 100) (Cross et al., 2015). Different methods for handling censored observations are discussed in Beal, 2001. In order to obtain conservative estimates for the time courses, we implemented a related adaptive replacement approach: we replaced censored observations with their QL value when the fitted value was above QL but considered a censored observation to match the fit when the fit itself was below QL as well (i.e. model predictions below the QL are not penalized). As up to 80% of observations per patient are below QL we generally weighted below QL observations with factor ½ in order to get robust patient fits. For bi-exponential fits of simulated BCR-ABL/ABL time courses, we weighted observations above 1 % and below 0.001 % (MR5) with factor ½ to reduce the influence of stochastic effects and to minimize their impact on the final fits.
Confidence intervals. We used parametric bootstrap to obtain confidence intervals for the mean LRATIO value for a patient i at any given time point. For a given data set, we started from the predicted LRATIO at the observed time points and repeatedly added a random measurement error sampled from a normal distribution with mean 0 and with variance " # . , the estimated residual variance from the bi-exponential model fit. We also applied the observed lower quantification limits at the affected time points. We repeatedly fit the model to the simulated response data and extracted the range of 95% of the predicted mean LRATIO values, thereby defining the confidence interval.
Group Comparison of imatinib versus dasatinib. For the assessment of treatment differences between imatinib and dasatinib therapy with respect to the LRATIO levels in the peripheral blood, we applied a non-linear mixed effect (NLME) model using maximum likelihood (ML) estimation with the stochastic approximation EM algorithm as implemented in Monolix. Model selection was applied as follows: We started with a bi-exponential model that incorporated solely patient-level random effects for the parameters α, B and β. We did not include a patientlevel random effect for the first intercept A as this parameter showed negligible patient to patient variability. Further model selection, i.e. the decision which treatment effects to include in the final NLME-model, was guided by Akaike's information criterion (AIC). The finally selected model contains random effects for all parameters but initial slope parameter A, while fixed treatment group effects for initial slope α and second slope β. We applied Wald tests to assess the statistical significance of fixed-effect group effects. Residual plots and the distribution of the predicted random effects of the final NLME-model did not reveal any evident violations of the assumptions of NLME-models. The model fit was robust with respect to the choice of initial start values for the parameters. Our retained model is summarized in Suppl. Furthermore, as it is known that measurements of high tumor load (i.e. BCR-ABL ratio > 10%) are more uncertain due to the usage of the ABL1-reference gene, we assessed the impact of these high tumor load measurements on our results and performed the identical analysis when BCR-ABL values > 10 % are either given half weight (Suppl. Time to event data. To analyze time-to-event data, such as time to progression or time to major molecular response (MMR, i.e. BCR-ABL below 0.1%), we used cumulative incidence estimates with death as a competing risk event. We defined high and low risk subgroups based on cut-off values of halfing time, initial slope or long-term slope , respectively. The cutoff values were defined by maximizing the average of sensitivity and specificity (Youden-index) on the data. We used the Gray-test to test for significant differences of the cumulative incidence functions between the high and low risk group (R-package cmprsk).

Simulation algorithm & model description
The described model of HSC organization and leukemia is implemented as a single-cell based, stochastic process. We used a discrete time step approach, with Δt = 1 hour, where every single cell is updated according to defined rules including stochastic decisions. The status and state of a model cell, irrespective of its type, is given by its affinity , its membership to a signaling context m ∈ {A, Ω} and its position in the cell cycle c. The cell cycle for an activated cell in the signaling context Ω is modeled as a sequence of G1 phase, S phase, and G2/M phase, in which the duration of latter phases (τ Z and τ [./4 ) is fixed. To realize an update step, the actual total number of cells with > ]AB in signaling context A (N _ (t)) and Ω (N`(t)) is determined. Based on these numbers and according to its actual signaling context, the status of each model stem cell is updated as follows: ( ) The cell changes to signaling context Ω or stays in A with probabilities ω and 1 − ω, respectively. If it stays in A, its cell specific affinity is increased by the factor (regereneration coefficient) up to ]Ee = 1. I the cell changes to signaling context Ω, its position in the cell cycle will be set to the beginning of S-Phase.
( ) The cell changes to signaling context A or stays in Ω with probabilities α and 1 − α, respectively. Herein, a change to signaling context A is only possible if the cell is in G1-Phase and > ]AB , ending up in G0-Phase (cell inactivity). If a cell stays in signaling context Ω, the affinity is decreased by the division of the differentiation coefficient and the cell cycle counter is incremented. In case the cell cycle is completed, i.e. c(t) = τ [ j + τ Z + τ [./4 ), c(t) is set to 0 and a new identical cell is generated (cell division). In case < ]AB cells can not reenter into the signaling context A and continue development in Ω (differentiation/maturation). iv All cells with > ]AB are denoted as stem cells, whereas cells with a ≤ a ]AB corresponds to differentiating cells. The differentiating cells undergo an proliferative phase of 480 hours before they mature for another 192 hours. BCR-ABL/ABL ratios are calculated based on the fraction of leukemic cells among this pool of maturing cells. The transition probabilities α and ω depend on the actual affinity (t) of the cell, on the fixed parameters ]AB and ]Ee and on the transition characteristics f 6 and f 2 , respectively: The transition characteristics f 6 and f 2 are crucial regulators for the description of the dynamic differences between healthy and leukaemic cells as well as the effect of tyrosine kinase inhibitors. Besides their cell type it depends also on the cell numbers of the target context. They are modelled by a general class of sigmoid functions: The parameters ν s , ν . , ν y and ν | determine the shape of the transition functions f 6/2 . Ñ _/` is a scaling factor of the corresponding compartment. ν s , ν . , ν y and ν | can be uniquely determined by In order to explain the competitive advantage of leukemic compared with healthy cells, we assume that the leukemic cells have an increased and unregulated proliferative activity. Technically, the transition characteristics f _ and f`, which describe the transit of cells between A and Ω, differ such that leukaemic cells have (I) an increased and cell number-independent probability per time step to be activated into cycle and (II) a slightly increased propensity to find an empty niche site, which is modelled by an increased probability to change to the niche-like signalling context A. Both assumptions are necessary to consistently explain characteristics of CML pathogenesis (Horn, Loeffler, & Roeder, 2008;Roeder et al., 2006). For a fully developed CML, this leads to a higher absolute number of leukemic cells compare to normal cells in the homeostatic situation. However, as the regulation of the normal cells is unaltered, those cells are limited in their proliferative activity and are therefore slowly outcompeted. This is an indirect effect caused by different regulating mechanism rather than a direct spatial competition. A conceptual sketch of the model is provided in Figure 2B of the main document.
Due to performance reasons, we precomputed simulations in certain parameter ranges, based on the results of previous studies. Technically, we varied the degradation rate rdeg, the rate of gradual TKI onset rtrans, the activation function f 2 345 for leukemic cells, the deactivation function f 6 345 for leukemic cells, and the initial leukemic tumour load LRATIOinit (see Suppl. Table S4). The simulation of all possible combinations lead to a parameter grid of 270,400 combinations = (f 2 345 , f 6 345 , r 89: , LRATIO ABAC , r CDEBF ), building the base of a look-up table. Given the full combinatorics and performing five simulation runs per design point to obtain a robust estimate of a typical BCR-ABL time course (which is a stochastic process) over 8 years of follow-up, we ended up with a total of 1,352,000 simulation runs. Suppl. Table S4: Parameter ranges as the basis of the lookup table. All possible combinations were presimulated in order to use efficiently for comparison with patients' data.

Figure S1
Supplementary Figure S1. Predictability of slope analysis for clinical outcomes. Cohorts are separated into a high and low risk group by maximizing the sum of sensitivity and specificity (ROC analysis) with respect to the intended outcome for dasatinib and imatinib individually. The determined cut points together with the results (i.e., p-values) of testing for equality of the cumulative incidence in high and low risk group are given for each treatment arm, respectively.   xiii Figure S4 Supplementary Figure S4. Comparison of dynamic modelling parameters. Subfigures show histograms of (A) the patient specific degradation rate rdeg , (B) the rate of gradual TKI onset rtrans , (C) the reduced activation function (fω CML ), (D) the rate by which leukemic stem cells re-enter into quiescence (fα CML ) and (E) the initial tumor load at therapy start LRATIOinit , and obtained for the most optimal parameter setting Q*min for each patient. Supplementary Figure S7. Prediction accuracy of estimated model parameters within the patient's confidence region omitting values with BCR-ABL ratio > 10%. False positives (FP) and false negatives (FN) rates for predictions of 5 year outcomes are given as a function of the available (shorter) observation period (n = 211). The slightly higher FP rate (compared to Fig. 4E and S6) results from the higher uncertainty (due to fewer measurements), while the overall trend of increasingly precise predictions with longer observation times remains unchanged.