Numerical analysis of critical parameter values for remission during imatinib treatment of chronic myelogenous leukemia

1.   Introduction

Chronic myelogenous leukemia (CML) is a form of blood cancer that stems from an overgrowth of myeloid cells in the bone marrow. Consisting of approximately 15% of leukemia cases ^[1], CML has been under heavy investigation in search of the best way to treat patients. Although allogenic stem cell or bone marrow transplantation is the only known curative treatment for CML ^[2], our work focuses on the treatment of the tyrosine kinase inhibitor imatinib mesylate, commercially known as Gleevec.

We categorize the success of the treatment plan based on the lowest total concentration of leukemia cells in the blood after the plan has been implemented. The three stages are hematologic, cytogenetic, and molecular remission, corresponding to having $10^6$, $10^8$, and $10^{10}$ leukemic cells, respectively. Assuming that the average person holds approximately 6 liters of blood, the concentration cutoffs for these three stages are shown in Table 1.

Under continuous imatinib treatment, most patients with CML attain some form of remission. Specifically, almost all patients achieve hematologic remission and 75% achieve cytogenetic remission ^[3]. However, CML is never completely eliminated from imatinib alone and will eventually return after stopping treatment ^[4,5]. Due to this, there are concerns about lingering leukemic cells that acquire imatinib resistance and weakening of the immune system as a result of prolonged usage ^[6,7,8,9].

Many mathematical models have been developed to study the dynamics of interaction between CML and an anti-leukemia immune response under various treatment strategies to mathematically realize clinical observations. In Bunimovich-Mendrazitsky and Shklyar ^[10], the authors sought to solve an optimal control problem on imatinib dosages combined with interferon-$\alpha$ (IFN) injections, centered around a system of two differential equations: one for all CML cells, and one for a population of T-cells representing immune response. Bunimovich-Mendrazitsky et al. ^[11] focused on a system of three ODEs, consisting of two CML growth stages and a T-cell population, with a Heaviside function to capture the timing of IFN injections. Numerical insights are gained from simulating treatment plans of imatinib with and without IFN.

In this work, we focus on the DDE model presented by Kim et al. ^[5], which is an extension of the ODE model presented by Michor et al. ^[4]. The authors defined a system of delay-differential equations (DDE) with eight leukemia sub-populations (four non-mutated and four with acquired resistance, each with respective growth and death rates) under continuous imatinib treatment; and an anti-leukemia immune response as a T-cell population, with terms that describe interactions between CML and T-cells.

The primary motivation for our work comes from Paquin et al. ^[12]. The authors conducted a sensitivity analysis on each model parameter in ^[5] under continuous imatinib with a strategic treatment interruption (STI) at some specified time (encoded as a reversion to unhindered CML growth rates for some specified duration). All parameters were simultaneously varied over a specified range using Latin hypercube sampling (LHS) and fed into the DDE system to simulate CML and immune dynamics, computing various correlation coefficients from each range.

Rather than focusing on simultaneous parameter changes, our work focuses on the effect of individual parameter variations toward improving the time to remission. This work aims to investigate the effectiveness of an STI implementation based on the range of possible values for model parameters that contribute toward curbing remission time under this setup. We compare these parameter ranges between numerical solutions to the model with and without an STI. We find the emergence of thresholds toward achieving remission, which are improved when introducing an interruption in imatinib treatment.

This paper is organized as follows: In Section 2, we verify existence, uniqueness, and positivity of solutions to systems of delay-differential equations. In Section 3, we introduce the delay-differential equations model in full, and how the model describes the assumed interactions between CML, imatinib, and the immune system. In Section 4, we investigate the presence of threshold values for various patient recovery metrics when individually varying parameter values while holding all others at a constant frame of reference. We compare thresholds both with and without an STI. In Section 5, we discuss future research and implications of these thresholds.

2.   Properties of a system of delay-differential equations

We begin with the one-dimensional case with a single delay $r > 0$; initial data $\phi: [s-r, s] \to \mathbb{R}$; and a given $s \in \mathbb{R}$. We are focused on the solution to the following initial-value problem:

As with systems of ordinary differential equations, preliminary properties (existence, uniqueness, positivity of solutions) for a general system of delay-differential equations are established (and proved in the literature as such). The theorems in this section come from ^[13].

The following theorem reads similarly to that for existence and uniqueness of ordinary initial-value problems:

Theorem 1. Let $f(t, x, y)$ and $f_x (t, x, y)$ be continuous on $\mathbb{R}^3$, $s \in \mathbb{R}$, and let $\phi : [s-r, s] \to \mathbb{R}$ be continuous. Then there exists $\sigma > s$ and a unique solution of the initial-value problem (2.1) on $[s-r, \sigma]$.

The goal is to extend the local existence of a solution from Theorem 1 to be defined for all $t \ge s$. A noncontinuable solution is that for (2.1) that cannot be extended to a larger interval. The following theorem shows that if this solution is not defined for all $t \ge s-r$, then it will grow to infinity as $t \rightarrow \sigma$, where $\sigma$ comes from Theorem 1.

Theorem 2. Let $f$ satisfy the hypothesis of Theorem 1, and let $x:[s-r, \sigma) \to \mathbb{R}$ be the noncontinuable solution of the initial-value problem (2.1). If $\sigma < \infty$, then

Using the following remark, we can apply Theorems 1 and 2 to our 5-dimensional system.

Remark 1. Theorems 1 and 2 extend immediately to the case that $x \in \mathbb{R}^n$ and $f: \mathbb{R} \times \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}^n$ with almost no change in the proof. It also extends to multiple discrete delays $r_0 < r_1 < \cdots < r_m$ where $f = f(t, y(t), y(t - r_0), y(t - r_1), \ldots y(t - r_m))$ with very little change.

Finally, to ensure that our system reflects a biological system, we use the following theorem to verify that non-negative initial conditions yield non-negative solutions.

Theorem 3. Suppose that $f: \mathbb{R} \times \mathbb{R}^n_+ \times \mathbb{R}^n_+ \to \mathbb{R}^n$ satisfies the hypotheses of Theorem 1, Remark 1; and

If the initial data in (2.1) satisfy $\phi \ge 0$, then the corresponding solution $x(t)$ of 2.1 satisfies $x(t) \ge 0$ for all $t \ge s$ where defined.

3.   A mathematical model for CML with imatinib treatment

We utilize a compartmental model presented by ^[5] to simulate the dynamics of CML, which extends a model from ^[4]. Both models begin by constructing the leukemia sub-populations that represent the CML cells at different life stages. The leukemia begins growing from stem cells to progenitor cells, then to differentiated cells, and finally terminally differentiated cells. We denote the concentration of cells in each sub-population at any given time $t$ with $y_0(t), y_1(t), y_2(t)$, and $y_3(t)$, respectively.

Within each sub-population, there are natural growth rates, death rates, and interactions from the immune response that move cells into and out of that stage. All sub-populations have growth rates into the next mature stage, with stem cells having an additional regeneration rate back into the same stage. The regeneration rate is $r_y$, and the growth rates from stem cells to progenitor, progenitor to differentiated, and differentiated to terminally differentiated are $a_y, b_y$, and $c_y$, respectively. The natural death rates at each stage are $d_0, d_1, d_2$, and $d_3$, respectively.

Under continuous imatinib treatment, the growth rates $a_y$ and $b_y$ are reduced by a factor of 100 and 750, respectively. To incorporate the possibility of drug resistance with prolonged use, both models introduce a mutation rate into the system. It is assumed that leukemic stem cells will morph into their drug-resistant counterparts at a rate of $u$ cells per day, which begins the construction of a new set of drug-resistant stem cells, progenitor cells, differentiated, and terminally differentiated cells, denoted as $z_0, z_1, z_2$, and $z_3$, respectively. Similarly to the non-mutated system, the corresponding growth rates are $a_z, b_z, c_z$ with stem cell regeneration rate $r_z$. We assume the death rates are identical to the non-mutated system.

Unlike in ^[4], the model presented in ^[5] considers the body's natural immune response toward CML, conveyed as a T-cell population and interaction terms with each leukemia sub-population. Denoting its population $T(t)$ at any time $t$, the changes in concentration are governed by a natural supply rate $s_T$, death rate $d_T$, those lost to stimulated division in the presence of CML, and those gained via stimulated division after a time delay of $\tau$ days per division. We illustrate the movements of all cell populations in Figure 1.

The interaction term $q_C p(C, T)y_i$, or $q_C p(C, T)z_i$, represents the concentration of the leukemia sub-population killed after a successful encounter with a T-cell. Expanding $p(C, T) = p_0 e^{-c_n C} kT$ helps to describe a successful interaction. Upon encountering a cancer cell with probability $p_0$, with a rate of interaction $k$ a T-cell is stimulated to divide, with it and its predecessors returning to the population after a time delay of $\tau$ days and the cancer cell dying with probability $q_C$. If $n$ cells are ready to divide, they bring in $2^n$ cells after $n\tau$ days. Finally, the model incorporates negative pressure, a biological mechanism that reduces successful T-cell interactions with increased cancer concentration, to govern the $p(C, T)$ term in each differential equation. In ^[5], the authors chose the exponential function $e^{-c_n C}$ to represent immune down-regulation, where $c_n$ represents the decay rate of immune responsivity.

3.1.   Existence, uniqueness, and positivity of CML system solutions

Using the theorems from Section 2, we verify existence, uniqueness, and positivity of solutions to the simplified model with no assumed mutation rate:

Proposition 1. The CML system, when $u = 0$, has a unique, positive solution defined for all $t \ge 0$.

Proof. Arrange the 5 non-mutated DDEs in the CML system as $x'(t) = f(t, x(t), x(t - n \tau))$. Let $x(t) = (y_0(t), y_1(t), y_2(t), y_3(t), T(t))$ denote the order in which cell populations are written in the solution. Expanding $p(C, T)$, each $f_i$ is defined as follows:

It is clear that each $f_i$ consists of continuously differentiable functions (note the state variables $y_0, y_1, y_2, y_3, T$ are assumed to be differentiable, and hence continuous). Setting our initial data $\phi$ to the initial conditions for the model, local existence and uniqueness are established by Theorem 1 (and Remark 1).

To verify positivity, suppose that $y_0 = y_1 = y_2 = y_3 = T = 0$. Then $f_1 = \cdots = f_4 = 0$, and $f_5 = s_T > 0$. Since the initial conditions are all positive, then by Theorem 3, the corresponding solution is non-negative for all $t$ where defined.

Finally, we verify that our solution is defined for all $t \ge 0$ with Theorem 2, i.e., that $\sigma = \infty$. We assess bounds on each cell population. Since $y_0(t) \ge 0$, then $y_0'(t) \le (r_y - d_0) y_0$. Using the variation of constants formula, we obtain

With this,

Similarly:

And:

Finally:

By Theorem 2, these estimates tell us that the solution $x(t)$ is defined for all $t \ge 0$. □

3.2.   Computational implementation

Due to the system's nonlinear nature, we numerically simulate the DDE model in MATLAB (Version 2022b) using dde23. A simulation under continuous imatinib treatment is shown in Figure 2. The observation of cancer concentration rebounding after approximately 24 months is the primary motivation for introducing STIs in the first place. To implement an STI in MATLAB, we momentarily revert growth rates $a_y$ and $b_y$ to pre-imatinib values (1.6 and 10, respectively) within a specified window, and then re-divide by 100 and 750 afterwards. To calculate times to remission, we indicate the first instance that the CML concentration falls below the corresponding threshold from Table 1. When incorporating an STI, we indicate the first instance after the interruption that CML falls back to a remission stage. If the setup is so that the concentration is already below a remission stage at the beginning of the interruption, then the time to that remission stage is defaulted at 300 days.

4.   Classification of critical thresholds for remission

Using the model described in the previous section, we investigate the presence of influential model parameters based on how they contribute to patient recovery. Specifically, the presence of critical threshold values can inform the effectiveness of a specific treatment plan. For this analysis, we focus on the time it takes for the system to achieve remission over the course of the simulation.

Prior to imatinib treatment, CML patients typically have approximately $10^{12}$ leukemia cells ^[14]. The three standard types of remission are a 100-fold difference from each other, corresponding to $10^{10}$, $10^{8}$, and $10^{6}$ leukemic cells. Using the general assumption that the average person has approximately 6 liters of blood, we obtain remission levels as concentrations, highlighted in Table 1. We are primarily focused on the time it takes for patients to achieve molecular remission or better (where total elimination equates to less than 1 CML cell in the blood).

We define a parameter $\alpha$ to have a critical threshold value toward achieving remission stage $R$, denoted as $\alpha_R$, if $\alpha_R$ acts as an upper or lower bound of possible $\alpha$ values contributing to achieving a certain remission stage.

● If for any $\alpha_i > \alpha_R$, the time for the CML system to achieve remission stage $R$ is not defined (meaning remission stage $R$ is not achievable for values of parameter $\alpha$ larger than $\alpha_R$), then $\alpha_R$ is an upper critical threshold value.

● If for any $\alpha_i < \alpha_R$, the time to remission stage $R$ is not defined, then $\alpha_R$ is a lower critical threshold value.

We focus on the remission threshold values for molecular remission and total elimination, denoted as $\alpha_M$ and $\alpha_E$, respectively.

Classifying parameters based on exuding either an upper or lower critical threshold value also has a biological interpretation. We call those with upper critical threshold values CML-specific parameters, as they represent characteristics regarding CML growth or a more compromised immune response. If these values are too large, the interactions of CML and immune responsivity overwhelm the T-cells, preventing the system from attaining remission. On the other hand, parameters with lower critical threshold values are called immune-specific, as they describe processes relating to the anti-leukemia immune response. If these values are too small, the success of T-cell interactions and stimulated divisions is hindered, preventing remission. Examples of such parameters in either category are listed in Table 2.

To analyze parameter thresholds, we start with a set of reference values coming from universal estimates, as well as a set of patient-specific numerical values; see ^[4], "Patient 4". Although the base state in our model corresponds to a specific patient, the range of numerical parameters we consider encompasses a broad spectrum of biologically meaningful values, ensuring that our analysis remains generalizable and applicable to diverse patient populations. The parameter ranges used in our numerical simulations follow key prior work in the mathematical modeling of CML and its interaction with the immune system, particularly the models presented in ^[5] and ^[12]. For each parameter, we use a range of values (either flat or percentage-based, coming from ^[12]) to vary over, while fixing all other parameter values at the base state. The parameter ranges that we consider in this study correspond to those analyzed using Latin hypercube sampling in these works. For each simulation, we numerically simulate the solution to the DDE system over the course of 3000 days. The model parameters, their base values, and their ranges of interest are summarized in Table 3. We compare thresholds over various patient recovery metrics to simulations with a single, 15-day STI from $t = 300$ to $t = 315$ days, as in the comprehensive sensitivity analysis in ^[12]. We focus here on the 15-day STI from day 300 to day 315 as implemented in ^[12], where it was shown through numerical simulations and a comprehensive sensitivity analysis using Latin hypercube sampling to be particularly effective in promoting immune system activation and achieving deep molecular remission or elimination in patients with CML. Figure 3 shows this simulated setup. In this case, the system (patient) is able to attain and maintain a low CML concentration (specifically, total elimination of CML; remission stages are summarized in Table 1).

We acknowledge that the patient tradeoff in implementing an STI is the initial spike in CML after interrupting imatinib. We aim to implement a treatment plan that minimizes both the maximum CML spike and the time out of any remission stage, post-STI; more work must be done to find this optimal timing, duration, and quantity, which Paquin et al. ^[12] moves toward.

We begin by exploring a particular parameter with a notable threshold value, and then expanding the process across all parameters that express thresholds to contributing to remission, comparing critical values to that of simulations that incorporate an STI. As will be seen, the threshold value to a certain remission stage (or complete elimination) may fall outside of a biologically feasible range of values for that parameter. This can illustrate the strength of an STI in bringing the system to a state that cannot be further altered by a more compromised immune system or stronger strain of CML. Additionally, while some of the threshold values identified in our numerical simulations lie outside the range of what is currently considered biologically feasible, these results still offer important insights for the development of future treatment strategies. In particular, identifying which parameters play a critical role in achieving remission—even if only under extreme or currently unrealistic values—can suggest new directions for biological investigation or therapeutic intervention. For instance, an immune-specific parameter threshold that is currently infeasible may highlight a biological process that, if better understood or more effectively targeted, could eventually be modulated to expand the feasible range. On the other hand, infeasible CML-specific parameter thresholds support the effectiveness of STIs in treating CML. Thus, our results can inform innovative approaches to treatment design, motivate experimental work to explore biological plasticity, and help shape future therapies that go beyond the limitations of current clinical protocols.

4.1.   Case study: $p_0$ and $d_T$

To illustrate our process, consider the parameter $p_0$, defined in Table 3 as the probability of a T-cell engaging a cancer cell. Sampling values 25% above and below the universal estimate in increments of 0.01, we obtain Figure 4(a). With continuous imatinib treatment where $p_0 = 0.8$, patients are at risk of a relapse after prolonged dosages, and although molecular remission is obtained for large encounter probabilities, total elimination of CML is never achieved under our base state.

For this investigation, there is a threshold value that $p_0$ must cross in order to positively contribute to decreasing the time to molecular remission. By implementing a single, 15-day STI, the system is able to not only remove the $p_0$ threshold for molecular remission, but also introduces a way to achieve total elimination of CML by varying this immune-specific parameter.

In the context of CML-specific parameters, an STI is also able to improve thresholds by allowing the system to take on more stress. Consider the T-cell death rate $d_T$. With a large enough value beyond the base state of 0.0022 T-cell deaths per day, Figure 5(a) highlights the system unable to achieve any type of remission beyond the cytogenetic stage. With an STI, however, threshold values emerge as upper bounds for molecular remission and total elimination, shown in Figure 5(b).

4.2.   Overall effectiveness and threshold comparison with strategic treatment interruption

We continue this analysis with all model parameters that initially possess a threshold. To organize our results, we can categorize all parameters of interest as immune-specific or CML-specific. In the context of patient recovery, immune-specific (such as $p_0$ and $q_C$) thresholds are preferably lower with an STI, implying that the patient does not have to be as healthy to attain remission under this treatment plan. Similarly, CML-specific thresholds (from parameters such as $d_T$ and $r_y$) are preferably higher with an STI, showing that the patient can still achieve remission with the treatment plan under worsening conditions.

In this section, we highlight the effectiveness of incorporating an STI into the CML system at achieving total elimination, as well as how these molecular remissions (as upper or lower bounds, depending on the type of parameter) compare with those from the system under continuous imatinib treatment. For both investigations, we focus on a subset of the model parameters in Table 2, which exhibit remission thresholds within the range specified in Paquin's sensitivity analysis.

We begin by showcasing how implementing this interruption improves thresholds for both immune-specific and CML-specific parameters. The pattern follows those from the case study discussed in Section 3.1. With this treatment interruption setup, both molecular and elimination thresholds from immune-specific parameters decrease relative to their respective base values, while thresholds for CML-specific parameters are much more impacted. (See the caption in Figure 6 for exact percentage changes.)

To compare both treatment plans, we measure the percent change from non-STI remission thresholds to the corresponding remission threshold in a system with a treatment interruption. Compared to the continuous treatment system, the STI system improves all but one parameter, $y_0(0)$, in the same trend as in Figure 6. All immune-specific molecular thresholds are reduced, allowing the system to achieve remission with compromised immune responses; and most of the CML-specific molecular thresholds are increased, allowing the patient described in the system to recover with more debilitating CML characteristics. The only outlier is $y_0(0)$, which is 51% lower than its non-STI counterpart. (See Figure 7 for the exact percentages.)

5.   Discussion and future work

The model in ^[5] has been used to more realistically capture the dynamics of CML by considering interactions between cancer cells and T-cells with time-delayed divisions. The authors in ^[12] use this model to conduct a sensitivity analysis over various parameter ranges and STI treatment plans. Using this model and motivated by these ranges, we gain insight into the critical behavior of the system by searching for where the system begins to transition from one state to another, rather than investigating where stability occurs.

In the context of patient recovery, this translates to finding when molecular remission or total elimination is achievable. Depending on the choice of model parameters, which can be categorized as inherently describing CML or the body's anti-leukemia immune response, critical threshold values emerge as either upper or lower bounds (respectively) for these stages of remission. The chosen treatment plan of continuous treatment with a single, 15-day treatment interruption from 300 days to 315 days is shown to improve thresholds of each parameter's non-STI counterpart, either by increasing upper bounds or decreasing lower bounds. This confirms the effectiveness of this treatment plan, providing implications for patients to recover with either a stronger strain of CML or a more compromised immune system. As a remark, some thresholds clearly lie outside a biologically feasible range; although we were able to approximate critical values, this does not necessarily imply that CML will be able to obtain these radical transformations. Some upper bounds drastically change (i.e., $d_1$ and $y_0(0)$) and may never be achievable, establishing the STI implementation as beneficial to the patient in the long term (with the drawback of enduring rising CML, and associated symptoms, with an initial period of time off medication).

Limitations to this analysis are the assumptions to start with our choice of base values and treatments for the model setup. For instance, our base state was chosen due to the corresponding patient having the most data (T-cell concentrations) recorded across their visits with researchers. Ideally, a universal estimate or more well-established, biologically feasible ranges for patient-dependent parameters such as $n$ and $y_0(0)$ are used for model simulation instead of patient data. From there, we can compare simulated T-cell concentrations with real-world observations. This ambiguity in model parameter reference points may have led to plots and tables with results that do not align with intuition, i.e., CML-specific parameters yielding decreasing remission times with larger values. Further, this may reveal that $y_0(0)$ does not behave well as a parameter, serving better as a fixed initial condition; future study may consider an alternate cancer metric, such as total CML concentration $C(t)$.

There are many avenues for future investigation, starting with increasing the complexity of the model and treatment plan. Our results stem from a setup that assumes no emergence of imatinib-resistant CML cells ($u = 0, z_0(0) = 0$). Preliminary investigations have shown that the effect of STIs are greatly altered when including more terms into $C(t)$, requiring a different quantity and duration. Even in the non-mutated setup, our choice of STI poses a natural question: What is the best treatment plan? So long as the CML concentration in Figure 3 never exceeds a concentration of 25 k/µL, there are many ways to experiment with the quantity and duration of STIs: multiple interruptions (see ^[15]), earlier start times, later cutoffs, etc.

Furthermore, sampling over these ranges with finer increments reveals intermediate thresholds. In the plots for $p_0$, for instance, a molecular threshold $\alpha_M$ occurs at a value of 0.77 with continuous treatment, but a system with that value will take upwards of 7 years to achieve molecular remission, during which drug-resistant strains of CML can emerge. The true (or medically realistic) threshold ought to be 0.78, where time to remission drops under 2 years. Since these values are very close, the percent changes in Figures 6 and 7 are left relatively unchanged. However, this brings up the question of classifying critical sub-values, two for each parameter, and the use-case for each. The validation of our methods against other clinical datasets may help in generalizing remission threshold trends observed with the current base state, and assist in refining parameter estimates, especially for patient-specific attributes.

Finally, there is insight to be gained from studying the analytical behavior of this DDE system in hopes to understand the behavior of a closed-form solution similar to that of ^[4]. Our current work focuses on finding a bridge between stability from a clinical perspective (the ability to achieve remission) and a dynamical systems perspective (asymptotic stability of the system at a critical point). In this work, we have proven general results on existence, uniqueness, and positivity of solutions of the delay differential equations model of the interaction between CML and the immune response, following the theory of ^[13]. In future work, we will follow motivation from Niculescu et al. ^[16] and Besse et al. ^[17] in solving for critical points of the CML DDE system, finding eigenvalues from corresponding Jacobian matrices, and classifying stability. Comparing results for a simplified, non-delayed version of our model may yield insight to the original system of DDEs. The model at the center of our work is highly nonlinear due to the exponential functions, so simplifying the DDE's to appear as the system from ^[16] will assist in deriving this information. See ^{[18,19,20,21]} for other formal sensitivity and stability analyses of immune dynamics for leukemia and other cancers.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was generously supported by the William and Linda Frost Fund in the Cal Poly Bailey College of Science and Mathematics.

Conflict of interest

The authors declare that there is no conflict of interest.