Differential response to cytotoxic therapy explains treatment dynamics of AML patients: insights from a mathematical modelling approach

i. Patient Data

We used qPCR NPM1 time courses of n=275 patients from the AML Registry of the University Hospital Dresden Carl Gustav Carus and two phase 3 trials of the Study Alliance Leukemia (SAL), namely AML2003 (NCT00180102) and AML60+ (NCT00180167). All patients gave written informed consent to participate in the study in accordance with the Declaration of Helsinki. Bone marrow aspirates were taken at irregular time points, usually at diagnosis and during primary therapy, with additional measurements during follow-up and at relapse. NPM1 quantification was processed as described in [26]. As there were two different platforms used for PCR-quantification (i.e., Lightcycler480 and Taqman 7500, both with 5’-nuclease assays) a correction factor between the two values was estimated and applied to ensure comparability between measurements (for details see Supplementary Material). In total, 69 patients presented with a molecular relapse and 57 patients did not reach remission. The median age was 53 years (range: 20 - 79). The median number of NPM1 measurements per patients was 5 (range: 3 - 21). The medium number of therapy cycles received was 4 (range: 1 – 10). Measurements below the detection limit of 0.001 % NPM1/ABL (i.e. the lower limit of detection of the method) were handled as 0.001 %. All measurements after allogeneic stem cell transplantation were removed from the data set, as this treatment has a major effect on the disease dynamics and is not part of this analysis.

Patients were only included in our analysis of the accuracy of relapse time estimation if they did not reach remission at all, reached remission and relapsed within two years, or presented with later measurements confirming sustained remission during the entire two years.

We recently suggested different quantitative characteristics to describe a patients’ molecular time course (in [4]). These are: the elimination slope (α), the minimal NPM1/ABL level after primary treatment (induction + consolidation) within 9 months after treatment start (n), the maximum slope during relapse phase (β) and the time of molecular relapse (d) (for explanation see Figure 3A). We defined molecular relapse as the approximated time point when the NPM1/ABL ratio exceeds the relapse threshold of 1 % (following the suggestion in [26]). This time was approximated using a linear interpolation between the last point below the threshold and the first one above.

• Download figure
• Open in new tab

Figure 1:

A Schematic overview of the mathematical model describing the dynamics of leukaemia-initiating cells and healthy stem cells in the bone marrow of an AML patient. B Leukemic burden after therapy depending on the ratio of leukemic activation to healthy activation . Dashed line shows detection limit for clinically measured leukemic burden. For x the time course is shown in Figure 1C. C Time course of leukemic and healthy cell numbers relative to the numbers at diagnosis with a .

• Download figure
• Open in new tab

Figure 2:

A Model fit to the data (PatientID = 104) and 95%-confidence interval of an example patient suffering a relapse. Mean absolute error (MAE) = 0.6 log10% B Model fit to the data (PatientID =3137) and 95%-confidence interval of an example patient staying in long time remission. MAE = 0.04 log10%. The red shaded regions are the times of treatment. The red marks indicate the fitted parameter values for the leukemic activation and the leukemic proliferation (p_l) with their corresponding 95% confidence interval. C Time courses of cell numbers relative to the total cell numbers at diagnosis of model fit to patient from plot (A). D Time courses of cell numbers relative to the total cell numbers at diagnosis of model fit to patient from plot (B). H_Q - quiescent HSCs, L_Q - quiescent LSCs, H_A - active HSCs, L_A - active LSCs. E Mean absolute error (MAE) in log10% for all patient fits as a measure of goodness of fit. Black dashed line indicates threshold for quantitatively good fitting patients. Gray dashed line indicates threshold for qualitatively good fitting patients. For all patients with higher MAE the fit was poor.

• Download figure
• Open in new tab

Figure 3:

A Schematic overview of an NPM1 time course. α: elimination slope during primary therapy; n: minimal NPM1 level after primary treatment within 9 months after treatment start; d: time until molecular relapse; B Spearman correlation coefficients between the fitted parameters (i.e. leukemic activation , leukemic proliferation p_l and their ratio) and the patients molecular time course characteristics with adjusted p-values. C Relapse-free survival for intermediate and favorable ELN risk groups. P < 0.01 (logrank-test). HR = 1.86; D Relapse-free survival for high and low ratio of leukemic activation and leukemic proliferation as estimated by the mathematical model with a threshold between high and low of 18. P < 0.0001 (logrank-test). HR = 5.47. E Boxplot for the comparison of the molecular relapse times of patients with high (n=116) or low (n=10) ratio of leukemic activation and leukemic proliferation as estimated by the mathematical model with a threshold between high and low of 18. P < 0.0001 (U-test) F 1-year relapse probability (sigmoid function) as estimated by a logistic regression. Corresponding patient data (•) indicate the 1-year relapse status (relapse/no relapse) depending on the ratio of the two fitted parameters .

The whole dataset used for the present analyses is available as supplementary material (Patient-Data1.csv with all patient specific information; PatientData2.csv with NPM1/ABL measurements).

ii. The mathematical model

We developed a mathematical model describing the HSC dynamics in the bone marrow of AML patients (see Figure 1A). Since AML dynamics are largely driven by a pool of leukemia-initiating cells, we restrict our model to the stem cell dynamics (c.f. [17]). Similar to the general modelling concept applied for healthy haematopoiesis and chronic myeloid leukaemia [9, 22], we consider two stem cell states with different proliferative activity, i.e. quiescent and proliferating cells. Cells in these two states are assumed to differentially respond to chemotherapy. While quiescent cells, both leukemic (L_Q) and healthy (H_Q), are considered insensitive to S-phase specific drugs, activated leukemic (L_A) and healthy cells (H_A) are targeted by those substances. Assuming a finite (stem cell) niche capacity for both states (K_Q/A) [27] the transition between the states is modelled as a function of the number of cells in the target state. Additionally, assuming different transition rates between the states for leukemic and healthy cells, the following system of ordinary differential equation (ODEs) describes the change of the number of cells (H_Q/L_Q) in quiescent state Q over time (t):

The cell numbers (H_A/L_A) in active state A are additionally influenced by the proliferation of the active cells with rate p_h/l, and by their dependency on the total amount of cells in the current state, modelled by the carrying capacity K_A. Furthermore, active cells can leave the stem cell state with a differentiation rate d_h/l or they can be targeted by chemotherapy with a constant kill rate c during treatment. The cyclic chemotherapy administration is realized in the model by a fixed chemotherapy kill rate, which can be switched on and off. These assumptions lead to the following model description in the active state A:

The main readout of the model to be compared with the patients NPM1/ABL measurements is the leukemic burden as the percentage of leukemic cells in both states combined with respect to all cells in the states:

To compare model results with the data, all values below the detection limit of the qPCR method of 0.001 % were set to this limit. As it is common to have NPM1/ABL ratios of up to 1000 %, because of the higher abundance of the NPM1 transcript compared to the reference gene [28, 29] the percentage of leukemic cells from the model was increased by a factor of 100 to translate them into the NPM1/ABL ratio and make the model results comparable with the measurements.

ii. Model parameters

Most values for the model parameters were taken from the literature (see below) and treated as constant values for all patients. To account for patient-specificity, the aggressiveness of the leukemic clone and the individual chemosensitivity were modelled by the leukemic proliferation rate (p_l) and the leukemic activation rate , which were estimated per patient. The transition rate of leukemic cells into quiescence was kept constant for all patients and was set to 0.002 per day. It was obtained by fitting the model to all patients multiple times, each time with a different inactivation rate to find the globally best fitting rate. The other model parameters were set to plausible values based on available literature data:

• Chemotherapeutic kill rate (c) was kept constant to ensure identifiability of the fitted parameters. It was set to 0.99 per day to account for the high total leukemic cell kill, as reported for an in vivo mouse study [30].

• Proliferation rate of healthy cells (p_h) was set to once every 25 days. It was estimated that healthy hematopoietic stem cells divide once every 30 days [31]. But as the actual proliferation rate in the model depends on the niche space in the current state and hence this rate would only be reached for an empty state, the rate was set to this higher value.

• Stem cell differentiation rate (d_h/l) was set to the value of the proliferation rate without influence of the capacity (1/30 per day) for both leukemic and healthy cells.

• Activation rate for healthy cells was set to 0.01 per day and the inactivation rate to 0.2 per day, based on the functions for these values from the original model [9].

• Carrying capacities for each state (K_A/H) were based on the value of 105 cells as previously described [9], which is about a factor 10 smaller than actual stem cell numbers in a patient’s bone marrow.

The competitive advantage of leukemic cells compared to healthy cells originates from their increased proliferative potential [32, 33], which is represented in the model as an increased leukemic proliferation rate of up to once every 5 days.

When fitting the model to patient data, the leukemic activation rate and the leukemic proliferation rate (p_l) were estimated by minimizing the weighted log-likelihood, where negative NPM1/ABL measurements were double weighted to account for the possibility of these values being far below the detection limit. Differences between values were calculated on the log10 scale. For all patients, 100 random starting points were generated for fitting to eliminate the dependence on the start value and to increase the probability of reaching the global minimum. In order to avoid overfitting, we assessed the identifiability of the fitted parameters using the profile likelihood (see [34] for details). In the corresponding likelihood landscape in Supplementary Figure S1 it is shown that the identifiability is given for both free parameters.

For an additional, easier interpretable measure of goodness of fit the Mean Absolute Error (MAE) was estimated by computing the mean divergence of the model from the measurements for each patient fit on the log10 scale.

iv. Statistical analysis

General bivariate correlations were quantified using the Spearman rank correlation coefficient (ρ). P-values for multiple pairwise correlations were adjusted using Bonferroni correction. To specifically measure the agreement between relapse times estimated with the model and estimated from the clinical data we applied the Lin’s concordance correlation coefficient (ρ_c). The (Mann-Whitney) U-test was applied to compare rank distributions (i.e., shift in the medians) of two independent groups of patients. The Kaplan-Meier method was used for relapse-free survival analysis, using the time until recurrence/death. The logrank-test was applied to compare the survival of two independent groups. A cox regression was conducted for hazard ratio (HR) estimation. The threshold of the parameter ratio to discriminate between the two groups was determined by maximizing the HR with the additional requirement of at least 10 target events (i.e. relapse) in each group. Odds ratios were determined by logistic regression. Identifiability of the model parameters was assessed using the profile likelihood approach, as described in detail in [34]. The likelihood landscape is derived by estimating the profile likelihood for all parameter combination. The 95%-confidence intervals for the model parameters were derived following the likelihood-based confidence interval definition by Meeker and Escobar [35].

All analyses were done with MATLAB R2018b (The MathWorks, Inc., 1994-2018).

i. Reduction of leukemic burden can be modelled without assuming a different chemotherapeutic kill effect on normal and leukemic cells

In order to reproduce the molecular dynamics of AML we developed a mechanistic mathematical model of the competition of leukemic and healthy stem cells in the bone marrow. This modelling approach is in line with earlier works in which a minimal, two-compartment model of HSC organization was applied to describe various phenomena such as stem cell competition, aging and leukaemia [22, 9, 23, 36, 11, 37]. The model describes two functional states between which the cells can reversibly transit depending on the capacity in the target state. While in the “active state” the cells proliferate and can potentially differentiate, they stay inactive in the “quiescent state” (Figure 1A). Both healthy hematopoietic stem cells (HSCs) and leukaemia-initiating cells (LICs) can adopt each state. However, while their transition dynamics between the two states differ, leukemic and normal cells are competing for shared resources (such as limited niche space). The parameters for the healthy cells were taken from the literature as described in Materials and Methods. The routinely used chemotherapeutics D-arabinosyl cytosine (AraC) and Daunorubicin effect all dividing cells during S-phase equally [38]. Therefore, our first question to be answered was, whether it is possible to assume an equal cell kill on both active HSCs and LICs during chemotherapy and still achieve a leukemic clearance. Therefore, we defined all parameters for the leukemic cells to be the same as for the healthy ones, while only the leukemic activation rate was increased systematically. Figure 1B shows the leukemic burden at the end of 5 chemotherapy cycles depending on the assumed value of leukemic activation. The figure indicates that an elimination of the measurable leukemic burden (for a detection limit at 10-3 %) can be consistently explained even without differences in the chemo-induced kill rate if an elevated activation of leukemic cells by a factor of 50 is assumed compared to healthy cells. The exemplary time course in Figure 1C confirms the rapid decline of leukemic cells during therapy, while the healthy cells stay at a high level. This adheres to the general concept that leukemic cells are almost constantly proliferating without periods of extended quiescence and reflects the typical leukemic phenotype [39]. It shows further that this model is able to reproduce a broad bandwidth of treatment induced reduction of leukemic burden (reaching from 1 to 6 log-reduction), which is also confirmed in different AML patients during chemotherapy [29].

ii. A simple mechanistic model is able to mimic the molecular dynamics of most NPM1-mut AML patients

After having demonstrated the general ability of the model to mimic the leukemic reduction during therapy we analysed whether the model can also account for the interpatient-variability in treatment response and dynamic behaviour. For this analysis, 1567 NPM1/ABL measurements of 275 patients with a median number of 5 measurements per patient were available. We assessed the ability of the model to reflect a large variety of time courses, including remission (n=149), unresponsiveness to therapy (n=57) and relapse cases (n=69). Therefore, we adjusted the individual activation rate of leukemic cells , as well as the individual leukemic proliferation rate p_l for each patient’s time course of leukemic burden to find a patient specific parameter combination that optimally describes the individual disease dynamics, while all other parameters were kept constant. In Figure 2A and 2B there are two examples of fitted patient time courses along with the 95 % confidence interval for the optimal fit. These examples illustrate that the model reproduces the characteristic time course of AML patients presenting with either a molecular relapse or a sustained molecular remission. Although, both patients received very similar chemotherapies, it can be seen that for the relapsing patient the NPM1 burden is not eliminated after the therapy (Figure 2A), while for the other patient (Figure 2B) the value is below the detection limit. Hence, the slightly higher leukemic activation and the lower leukemic proliferation (p_l) are sufficient to explain the difference between relapse and sustained remission. Further, we observed, that a high leukemic activation rate is linked to a fast decline of leukemic burden during therapy, an effect which is illustrated clearly by the good responding patient in Figure 2B. Therefore, we conclude that this rate estimate can be interpreted as a measure for the patient’s individual chemosensitivity. Also, we observed, that a high leukemic proliferation rate is linked to fast leukemic regrowth, which can be seen in the poorly responding patient in Figure 2A. From this we conclude that the estimated leukemic proliferation rate provides a measure of the tumour intrinsic aggressiveness.

Looking in more detail how the different cell populations within our model respond to the treatment (Figure 2C and 2D), one observes that the difference between a good and a poor therapy response is within the leukemic cell population, with the good responder showing steeper decline during therapy as well as slower regrowth afterwards. In the poor responder the number of leukemic cells in the model also drops to a vanishingly low proportion of their initial numbers during treatment. Nonetheless, a relapse occurs within a few months, because of the fast regrowth. Furthermore, a decrease of the net proliferation rate during leukemic regrowth after chemotherapy can be observed (as the curves (L_A and L_Q) flatten after an exponential growth phase around 350 days after diagnosis). This observation agrees with a recent study by Akinduro et al., where the authors analysed the growth of AML cells after transplant in a mouse model [40]. The high number of quiescent HSCs is also in accordance with experimental data, as it was recently shown that in a leukemic mouse bone marrow HSCs were mostly in a quiescent state [41]. This mechanism is also referred to as “self-protection” by which HSCs avoid environmental effects in the quiescent state [42].

The quality of the individual fits was evaluated using the mean absolute error (MAE) of each model fit, which corresponds to the average residual between model fits and available data points on the log10 scale. We observed that fits with a MAE of up to 1 could quantitatively describe the patient’s course of disease. Hence, at least 227 of the 275 patients can be sufficiently described by the individual model fits (Figure 2E). Patients with an average error of up to 1.5 log could still be qualitatively described (n = 260).

Looking more closely at the 15 patients in which the model could not optimally fit to the time course data, it becomes clear that for most of them (10) the model as it is cannot account for the strong leukemic reduction found in the data (example patient in Supplementary Figure S3A). For 2 of the patients the very fast regrowth is not reproducible with the model and for 3 patients a chemo-resistance arises during treatment of a relapse, which is not yet captured by our current model (both can be seen in example patient in Supplementary Figure S3B).

Taken together, we could show that the proposed simple model is indeed able to reproduce most of the molecular behaviour of NPM1-mut AML patients, which can now be used to analyse the molecular characteristics of the course of disease.

iii. The model parameters are closely linked to time course characteristics

Based on the patient individual model fits we obtain a pair of model parameters (leukemic activation rate and leukemic proliferation rate p_l) for each specific patient time course (distribution of the parameters in the cohort can be found in Supplementary Figure S2). In order to investigate how these model parameters correlate with typical, clinically relevant characteristics of the NPM1/ABL dynamics, we perform a systematic correlation analysis. For a quantitative description of the relevant NPM1/ABL time course characteristics we use a parametrization based on the elimination slope (α), the minimal NPM1 level after primary treatment (induction + consolidation) within 9 months after treatment start (n), the maximum slope during relapse phase (β) and the time until molecular relapse (d_d) (see Figure 3A, [4]). Figure 3B provides a summary of pairwise comparisons between the individually fitted model parameters and the characteristics of treatment response. Generally, for all characteristics but the relapse slope β, a significant correlation with a correlation coefficient of up to 0.6 was found. Furthermore, the leukemic activation and leukemic proliferation (p_l) are conversely correlated with all described time course characteristics. That means, that later molecular relapse is linked to increased leukemic activation, but decreased leukemic proliferation. This supports the earlier conclusion, that the leukemic activation is a measure for chemosensitivity, whereas the leukemic proliferation provides a measure for the aggressiveness of the disease. Furthermore, this anti-correlation with the time course characteristics suggests that the ratio between the two parameters reflects a good measure for the overall severity of the patient’s individual course of disease.

From a clinical perspective it is most important to know whether and when a patient will relapse or not, and to estimate the chance for relapse-free survival. Using the approximated time point of molecular relapse as a surrogate for haematological relapse (including irresponsiveness), we studied whether a division into high and low ratio of estimated leukemic activation and leukemic proliferation (i.e. ; threshold at , determined by maximizing the difference between both groups, details in Materials and Methods) can improve the current ELN risk stratification (stratifying into high and low risk). Figure 3C and 3D support the notion that both variables are indeed correlated with relapse-free survival. However, the differences in the Kaplan-Meier plots indicate that this effect is considerably more pronounced using the ratio of leukemic activation and proliferation. Comparing the hazard ratio (HR) between the two stratifications further supports this impression (HRratio=5.5 (95% CI: 2.8 to 10.4); HRELN=1.9 (95% CI: 1.2; 2.8)). A clear difference in the time of relapse can be seen, when directly comparing the molecular relapse times of the high and the low parameter groups (Figure 3E). Patients with a high parameter ratio tended to relapse in median 15.3 months later. Moving away from the division into two risk groups, we apply a logistic regression approach to estimate to which extend the estimated ratio correlates with the probability to relapse within the first year after diagnosis (Figure 3F). The odds ratio of 1.17 [95% CI: 1.11 to 1.22] indicates, that an increase in the ratio by 5 units more than doubles the chance of experiencing a relapse within the first year. However, the nearly linear relationship between the ratio and the relapse risk indicates that the ratio alone cannot reliably predict relapse for a particular patient, although it can be useful in combination with other measures to improve risk prediction.

iv. Clinically relevant characteristics can be estimated using the mathematical model

In the clinical context the time of molecular relapse is a highly relevant measure. For this reason, we compared the molecular relapse time approximated from the data (d_d, definition can be found in Material and Methods) with the estimated molecular relapse time predicted by the model (d_m). Figure 4A shows the corresponding scatter plot for all patients for which a relapse time within 2 years could be approximated (n=175). As expected, there is a good accordance for most patients, while only for 13 out of the 175 patients (7%) the estimated molecular relapse time from the model diverges more than half a year from the molecular relapse time approximated from the data. The mean divergence for all patients is 1.9 months. However, the larger differences for some patients raises the question why the model is incapable of mimicking their course of disease.

• Download figure
• Open in new tab

Figure 4:

A Scatter plot for the comparison of the molecular relapse time approximated from the patient data d_d and the molecular relapse time estimated from the fitted model d_m for each patient. ρ_c = 0.90. Perfect accordance is indicated with dashed red line. Gray dashed line indicates divergence by half a year. Points with divergence larger than half a year are marked blue. B Scatter plot for the comparison of the molecular relapse time estimated from the fitted model to the first 9 months d_m of the patients and the molecular relapse time approximated from the data of these patients d_d. ρ_c = 0.37. C Example fit (PatientID = 3751) and 95% confidence interval for a patient where the model is not able to capture the molecular course of disease. Dashed line shows remission/relapse threshold. d_d is 0 days (no remission reached), d_m is 396 days. D Example fit (PatientID = 3621) and 95% confidence interval for a patient where the sparseness of data points makes it impossible to reliably approximate the molecular relapse time. Dashed line shows remission/relapse threshold. d_d is 0 days (no remission reached), d_m is 326 days.

Studying the 13 critical patients more closely it appears that the model is insufficient to describe at least 6 of them for technical reasons, i.e. very fast regrowth, very high chemosensitivity or stable tumour levels, that are not depicted in the model (see example in Figure 4C). To explain a stable tumour level over a longer time (as in the example patient), an additional intrinsic leukemic suppression (such as an immunological component) would be necessary in the model. For the other 7 patients with poor accordance of the relapse times, it is the sparsity of the available measurements that limits a reliable approximation of the relapse times from the data. For the example patient in Figure 4D the measurements suggest that a remission could not be reached, but it is conceivable that the tumour burden decreases further during the time of therapy. Therefore, it is very likely that a temporary remission occurs at the end of therapy, as suggested by the model fit. But just from the available data both conclusions cannot be excluded. Therefore, it is possible that the estimated relapse time from the fit is closer to the true relapse time than the approximation from the data. Hence, we conclude that the model, although not able to perfectly describe every single course of disease, might indeed harbour additional information about the true course between available measurements.

v. Predicting the relapse time based on the early time course is not reliable

We could show that the ratio is correlated to the approximated time of molecular relapse and suggested that the model can help to estimate the particular time point of molecular relapse occurrence. Moving a step further, we raise the question whether an early (and clinically even more relevant) prediction of the relapse behaviour can be made with comparable precision. From a clinical perspective, 9 months after treatment start is a time point when a decision about future treatment options for a particular patient is still relevant. We chose this time point as a cut-off until which the models were fitted, given there are at least 3 data points within this time period (n=89). Figure 4B compares the relapse times predicted by the model based on these truncated time courses with the approximated relapse time from the complete data. It appears that reliable model predictions cannot be achieved for the majority of cases. In fact, for more than 40% of the patients the divergence is more than 6 months and the mean difference over all patients is 6.6 months. This visualization indicates that a reliable prediction of relapse occurrence cannot be achieved if only sparse measurements during the initial treatment response during of the first 9 months of treatment are considered. However, it cannot be excluded from our current study that additional diagnostic parameters together with close meshed short term time course data convey more information for better predictions.