Bayesian modelling of response to therapy and drug-sensitivity in acute lymphoblastic leukemia

Acute lymphoblastic leukemia (ALL) accounts for around 25% of pediatric cancers (age $\leq$ 18 years) and is the most common cancer among children [18]. The cure rate for children with acute lymphoblastic leukemia (ALL) has improved greatly over the past few decades. Ten-year survival for childhood ALL has gone from around 11% in the first half of the 1960s to more than 90% by 2010 thanks to advances in understanding of the disease and clinical trials that have greatly improved outcomes [24, 17]. Research in the past has revealed cytogenetic characteristics of patients’ leukemic blast cells that are associated with prognosis [3]. For example, children with ETV6-RUNX1-positive ALL have a better outcome from chemotherapy than those who are negative for this fusion gene [11]. There are now many ALL subtypes defined by mutations, gene fusions as well as transcriptome with some directly informing new treatment strategies [14]. In addition to genomic determinants of prognosis, pharmacokinetic variability, some of which is influenced by pharmacogenomics, contributes to heterogeneity in treatment response and toxicity, both systemically and at the cellular level [28, 10].

The general treatment strategy for leukemia and many other hematologic malignancies is to first reduce cancer cell burden substantially, with the objective of inducing clinical remission. This initial phase of a chemotherapeutic regimen is called the induction phase. The induction phase often includes sequential treatment with a combination of cytotoxic and/or targeted-therapies. The strategy uses intermediate indicators of treatment response after a first course of the drugs and adapts subsequent drug doses and combinations accordingly. The goal is to achieve remission by the end of the induction phase, without which the patient’s prognosis is much poorer.

In ALL, the absence of minimal residual disease (MRD) has served as an early indicator of benefit from the anti-leukemia chemotherapy for decades [7]. Methodologies for MRD measurement detect and quantify residual tumor cells beyond the sensitivity level of cytomorphology. MRD assessment can be refined by the evaluation of additional genomic markers [9]. Nowadays the clinical impact of MRD is widely accepted, and MRD is considered the most important prognostic factor in the management of ALL. The prognostic significance of MRD can vary, depending on the timing in which it is measured. Early assessment of MRD during and at the end of the induction phase is often used for clinical decision making when treating patients with ALL [5].

Assays for MRD have a lower limit of detection. A patient is defined as MRD negative if the patient’s MRD levels fall below a detection threshold, typically less than 1/10000 cells. It is common in clinical practice to adopt a threshold of $>$ 0.01% to define MRD positivity. This value represents a limit of detection by flow cytometry, and retrospective analyses have shown that patients whose MRD ¿ 0.01% after induction therapy have a greater risk of relapse and poorer prognosis [4, 2, 15]. It is possible to achieve a routine sensitivity of 0.001% by PCR in clinical samples and similar sensitivity may be achieved by flow cytometry with specific B-cell ALL subtypes [26].

The study population for this work comprises patients who participated in three clinical trials at the St. Jude Children’s Research Hospital. The treatment regimens call for MRD assessments on day 15 and day 42 (end of induction) for determination of further treatment. Patients with MRD of 1% or higher on day 15 receive more intensive therapy for the remainder of the induction phase. Further intensification is reserved for patients with 5% of more leukemic cells on day 15. Patients with so-called standard risk ALL who have MRD of 0.01% or higher on day 42 are reclassified as high-risk ALL.

Additional information on the participants in the studies, collected at the beginning of the treatment, is available, such as age, gender, white blood cell count (WBC), and ALL subtype, as well as results of an ex vivo drug sensitivity screening to characterize the patients’ sensitivity profiles to anti-leukemic drugs.

The aim of this work is to investigate the relationship between patients’ leukemia subtype, leukemic cell drug sensitivity, and clinical benefit as measured by the early marker of treatment effect (MRD). To this end, we specify a joint model in a Bayesian framework for the MRD outcomes collected at days 15 and 42 able to account for the censoring and time dependence between the MRD measurements. Our statistical model combines multiple complexities and allows inference from the full Bayesian model. The complexities include the following. The model allows for the presence of leukemic cells below the assay’s lower limit of detection by treating MRD assessments as left censored if recorded as $<$0.01%. The model also considers MRD assessments on days 15 and 42 jointly via an autoregressive model. For day 42 MRD, inference is conditional on the presence of MRD on day 15, since no patient who is MRD negative on day 15 becomes MRD positive on day 42. This conditional autoregressive model allows for clearer interpretation of covariate effects. The model also includes Bayesian clustering as part of posterior fitting. The clustering of patients is based on ex vivo measurements of the sensitivity of each patient’s own leukemic cells to the anticancer drugs the patient receives during induction. The investigation also includes the genomic ALL subtype of each patient. A final component of the model is the use of the horseshoe prior [6] for selecting important covariates from the large number of features examined. This class of prior distributions for the regression coefficients provides a way to identify important covariates, thanks to the strong shrinkage effect of its heavy tails. In particular, the density function of this distribution presents a singularity at zero, leading to a more aggressive shrinkage of small coefficients towards zero than other standard prior distributions, and leaving important larger coefficients unaffected. A performance comparison between the horseshoe distribution and the spike-and-slab approach of [12] is offered in [6], showing consistency of the two posterior variable selection results. For further discussion on prior elicitation for the horseshoe distribution, see [20].

The manuscript is organized as follows. Section 2 describes the ALL data, while Section 3 introduces the modelling strategy. In Section 4 we present the analysis results. Section 5 concludes the paper with a discussion.

The motivating application consists of data for patients treated in three prospective pediatric ALL clinical trials at the St. Jude Children’s Research Hospital. These studies were part of the long-standing Total Therapy clinical trials program at the institution, initiated in 1962, comprising clinical trials that continue to build sequentially on each other [23]. These studies have shown the benefit of full-dose chemotherapy and treatment directed at the central nervous system for treating pediatric ALL, with greatly improved outcomes over the years, currently achieving event-free survival that exceeds 90% [24].

Our data set includes MRD measurements for childhood ALL patients in Total Therapy studies XV ($n=192$), XVI ($n=428$), and XVII ($n=168$). The treatment regimens for the studies include many of the same drugs during the six-week induction phase. There are several anti-leukemic agents that are common to the three studies; prednisone, vincristine, daunorubicin, PEG-asparaginase, methotrexate, cyclophosphamide, cytarabine, and mercaptopurine. MRD measurements are made on day 15 and at the end of the remission induction phase of treatment (day 42).

The detection limit for MRD is 0.01% (1 leukemia cell among 10,000 normal cells in the bone marrow), below which the exact values are not considered accurately measurable. Figure 3 displays the observed MRD values at day 15 and 42 (on the $\log_{10}$ scale). The values of MRD smaller than 0.01 are therefore censored and are jittered around 0.001 in the figure to improve presentation. The detection limit for MRD poses a modelling challenge, which we address in this work. The day 15 MRD assessment determines subsequent treatment during the remainder of the induction phase. Additionally, patients whose leukemic cells exhibit genomic variants for which targeted agents are available receive these chemotherapeutic agents (e.g. a tyrosine kinase inhibitor for patients whose ALL cells harbor a chromosomal translocation creating a BCR/ABL fusion). In our analysis, we focus on the drugs that all patients receive.

Additional information recorded at entry include age, gender, white blood cell count (WBC), ALL subtype, and treatment protocol number. Characteristics of the patients that are available across the three studies are summarised in Table 1. The median age of the children is between 4.9 years and 5.6 years, typical for childhood ALL. There are slightly more male than female patients. Baseline white blood cell counts are somewhat lower in Total XVI than in the other two studies.

The genomic subtype categories are shown in Figure 4. In the data set, 23 subtypes are observed, with the most common subtypes being hyperploid and ETV6-RUNX1, typical for childhood ALL. Section 1 includes information about the fusion product ETV6-RUNX1, which confers a favorable prognosis. Hyperdipoid ALL is also a common subtype that is also considered to have a more favorable prognosis [19]. We merged several subtypes that included fewer than 10 subjects each to improve degrees of freedom for this analysis. For example, we merged the Ph-like CRLF2 and the Ph-like non-CRLF2 subtypes [25]. The analysis included 12 distinct subtypes. The detected MRD levels vary for different genomic categories, as shown in Figure 1, displaying the measured MRD values on a $\log_{10}$-scale at each time point and across subtypes.

The data set includes estimates of patient-specific cancer cell sensitivities to various drugs in the combination chemotherapeutic regimen. The estimates arise from ex vivo dilution assays of patient samples. In essence, the patient’s leukemic cells are plated with each drug via dilution assays to determine sensitivity. For each drug concentration, a count of the number of surviving cells provides an estimate of the lethality of that concentration to the patient’s leukemia cells. From the roughly six concentrations per drug, an estimated LC₅₀, the concentration that led to 50% of the cells dying, is used as the measure of the sensitivity of each patient’s disease to the specific drug. We included five anti-leukemia drugs in these analyses, selected after excluding the compounds for which the missing rate is more than 40%. The compounds used for this analysis are asparaginase, prednisone, vincristine, 6TG and 6MP. The entire ALL pharmacotype dataset was previously published [16]. To highlight the relationship between the observed MRD values and the estimated LC₅₀s, in Figure 2 we present the Pearson correlations between these two variables for each of the five compounds and for each subtype category, averaged over 100 imputed $\log_{10}$ LC₅₀ data sets (see Section 3).

Estimates of LC₅₀s for patients are not available for some medications if the number of available cells was not sufficient to test all medications. Therefore, we applied multiple imputation for the missing observations. We impute missing data via the R package mice [27].

The main goals of the statistical model are: (i) to characterize the MRD responses at both time points accounting for their time dependence and the censoring; (ii) to investigate the effect of baseline and in-trial covariate information on the MRD levels; (iii) to identify subgroups of patients presenting similar sensitivity to drug exposures and (iv) to investigate the relationship between the clustering structure and drug-subtype interactions.

Data are available for $N=788$ children, for whom both baseline and in-trial information are available. Covariates included in the analysis are age, gender, WBC at diagnosis, ALL subtype and therapy protocol, as well as drug sensitivity profiles obtained from ex vivo dilution assays. The regression model included eleven indicator variables for the twelve ALL subtypes. For each patient, at most one of these eleven parameters will equal 1. The baseline category is ETV6-RUNX1, the most prevalent subtype. If a patient is positive for ETV6-RUNX1, all eleven subtype variables will be 0.

The ALL data set presents two main statistical challenges: (i) censoring of the MRD levels due to the assay’s lower limit of detection (i.e., $<$ 0.01%); (ii) estimation of each patient’s drug sensitivity profile. As stated above, MRD values above the threshold 0.01% are observed. We chose to treat the lack of MRD as a censored observation instead, allowing for the fact that there may still be residual disease but at an undetectable level. This consideration is particularly important when introducing temporal dependence. In the likelihood, censoring can be accounted for using standard strategies from survival analysis and censored observations are imputed as part of the Markov chain Monte Carlo (MCMC) algorithm.

We account for the repeated measures aspect of the MRD assessments by including an autoregressive term (of order 1) to model dependence between the MRD at day 15 and day 42, only for those observations for which MRD at day 15 is above the detection threshold. Otherwise, the two time points are treated as conditionally independent. We opt for this strategy because in the dataset if MRD at day 15 is censored (i.e., $<$ 0.01%), then MRD at day 42 also falls below the detection threshold (see Figure 3). Introducing dependence among the observations censored at both times would bias the estimate of the temporal effect as well as of the covariate effects, since the statistical model includes covariate information via a linear regression term.

Another important feature of our modelling strategy is the inclusion of a patient-specific measure of the patient’s leukemic cells’ sensitivity to individual drugs in our model of MRD. Drug sensitivity profiles for each patient are estimated from ex vivo drug sensitivity assays. The model includes a finite mixture model for the profiles to allow the data to perhaps determine categories that relate to overall sensitivity or resistance to one or more of the anticancer drugs. A major advantage of the Bayesian framework is that it allows for joint estimation of the MRD model and the drug sensitivity profiles, enabling a probabilistically sound quantification and propagation of uncertainty.

We describe in detail the statistical strategy below.

Posterior inference under the proposed model is obtained via JAGS [22], interfaced with RStudio via the R package rjags [21]. We run the MCMC chain for 15000 iterations, discarding the first 5000 as burn-in and thinning every two.

Figure 7 displays the scatterplots of the posterior samples of the coefficients $\rho_{0}$, $\rho$ and of the variances $\sigma^{2(1)}$, $\sigma^{2(2)}$. The posterior distribution of the auto-regressive coefficient $\rho$ is concentrated around the value 1 (posterior median 0.98, posterior interquartile range 0.92-1.05), indicating that the responses at the two different times are strongly associated for those subjects with $\text{MRD}^{(1)}>0.01$. Furthermore, the residual variances $\sigma^{2(1)}$ and $\sigma^{2(2)}$ show little a posteriori correlation.

We now discuss the effect of the covariates $\bm{X}_{i}$ and $\bm{c}_{i}$ on the responses. The posterior distribution of the coefficients $\left(\bm{\beta}^{(t)},\bm{\gamma}^{(t)}\right)$ for $t=1,2$ are summarized in Figure 9, where we report their posterior means and 95% credible intervals. Each bar corresponds to one of the covariates included in the analysis. Most of the covariates are associated only with MRD measured at time $t=1$ (i.e., their 95% credible interval does not contain the value zero). Interestingly, we observe that age, $\log_{10}$ WBC, and protocol are all relevant factors at day 15, while gender is not. The ALL subtype of the patient is also a relevant factor. The subtypes’ positive posterior regression coefficients on day 15 show substantially higher MRD on day 15 for these subtypes, relative to the reference subtype ETV6-RUNX1, as expressed by many non-zero coefficients. This holds also for the drug sensitivity profiles (see the posterior distribution of variables Cluster 2 and Cluster 3), relative to Cluster 1, the reference cluster. We see that Cluster 2 is associated with lower MRD on day 15 that Cluster 1, while Cluster 3 patients tend to have higher day 15 MRD values than Cluster 1 patients. Cluster 2 has more ETV6-RUNX1 and hyperdiploid patients, which are associated with sensitivity to prednisone and asparaginase, consistent with previous publications.

Recall that inference for day 42 depends on the presence of MRD on day 15. We see that Protocol T15 is associated with a higher mean MRD value than the other two studies. T16 and T17 occurred later and included some newer targeted anticancer drugs. We also see an improvement in the effect for the hyperdiploid subtype on day 42, relative to the reference genomic subtype. We also present posterior inference on the censored MRD values, which are imputed within the MCMC algorithm. In Figure 10 we show the posterior means for the censored observations, as well as the observed MRD values. Note that the bottom left quadrant corresponds to MRD values censored at both time points, while the bottom right panel to observations censored only at day 42. No observation lies in the upper-left quadrant, indicating that none of the subjects in the study experiences a relapse between days 15 and 42.

Clustering of the subjects is achieved via model (3), in which the sensitivity profiles are modelled via a three-component mixture distribution. The posterior estimate of the partition, obtained by minimizing Binder’s loss function as implemented in the R package salso [8], presents three clusters of sizes 271, 267 and 250, respectively. We provide an interpretation of the cluster features by investigating the average values of $\log_{10}$ LC₅₀s within each cluster, as well as the proportion of subjects characterised by each subtype. The results are summarised in the radar plots in Figure 11. We observe differences among clusters in the values of the $\log_{10}$ LC₅₀ relative to the drugs prednisone and asparaginase, and in the proportion of subjects with subtype ETV6-RUNX1. Cluster 2 is associated with the lowest MRD on day 15 (best prognosis). This cluster has more ETV6-RUNX1 and hyperdiploid patients, which are associated with sensitivity to prednisone and asparaginase, consistent with previous publications (e.g., [13, 1]).

Minimal residual disease is used in clinical practice to assess treatment response in patients with ALL. MRD is typically regarded as an independent prognostic factor and is used in clinical trials for risk assignment and to guide clinical therapeutic choices in the management of the treatment of the cancer. In this work we combine information from MRD measurements taken at two time points with patient-specific drug sensitivity estimates and leukemia genomics, with the goal of better understanding treatment response and the role of potential prognostic or predictive (i.e., relevant for treatment selection) factors. Our results highlight that clustering of individuals in terms of sensitivity to drugs that are part of the induction phase of therapy is mainly driven by sensitivity to two anticancer agents, prednisone and asparaginase. The analysis also highlights the role of genomic subtypes as important determinants of response to induction therapy as measured by MRD.

From a methodological point of view, the Bayesian framework allows us to jointly estimate the parameters that relate MRD at different time points during the course of therapy to patient-specific covariates, including simultaneously inferred drug sensitivity profiles. Moreover, it is straightforward to account for censoring. A key component of our modelling strategy is that we model MRD at day 42 conditionally on the event that MRD on day 15 is detected. This allows for extra flexibility and accounts for the fact that in our data set remission at day 15 implies remission at day 42.

A drawback of our modelling strategy is pre-imputation of missing LC₅₀ values, which leads to underestimation of uncertainty. We opted for this strategy because of the substantial missing rate. In theory, it is straightforward to impute missing observations within the MCMC algorithm. Future work includes a further investigation of the relationship between subtypes and drug sensitivity. In the data, estimation of LC₅₀ occasionally led to estimates below and above the range of concentrations for the patient. Our model did not account for the censoring in estimation of LC₅₀s for affected drugs with these patients’ samples.

Another limitation is that the analysis was confined to a subset of ALL subtypes, pooling subtypes with fewer than ten patients. We did consider hierarchical modeling but still pooling. Future work will consider all ALL subtypes as separate categories, including those with small numbers. A robust multi-level model may be useful to account for the uncertainty of the rarer subtypes. Such an extension would be useful, as we can imagine that more and more “subtypes” with very few patients will be defined in the future.

Future research will consider the incorporation of a greater number of serial MRD measures in the model. Such a strategy could include measures after remission induction that could be part of monitoring the disease in peripheral blood rather than bone marrow and help fine tune treatment decision rules.