Remission Stage | Hematologic | Cytogenetic | Molecular | Total Elimination |
Concentration, k/µL | <1.67 | <1.67⋅10−2 | <1.67⋅10−4 | <10−10 |
Chronic myelogenous leukemia (CML) is a cancer of the white blood cells that results from uncontrolled growth of myeloid cells in the bone marrow and the accumulation of these cells in the blood. The most common form of treatment for CML is imatinib, a tyrosine kinase inhibitor. Although imatinib is an effective treatment for CML and most patients treated with imatinib do attain some form of remission, imatinib does not completely eradicate all leukemia cells, and if treatment is stopped, all patients eventually relapse. Kim et al. constructed a system of delay differential equations to mathematically model the dynamics of anti-leukemia T-cell responses to CML during imatinib treatment, and demonstrated the usefulness of the mathematical model for studying novel treatment regimes to enhance imatinib therapy. Paquin et al. demonstrated numerically using this DDE model that strategic treatment interruptions (STIs) may have the potential to completely eradicate CML in certain cases. We conducted a comprehensive numerical study of the model parameters to identify the mathematical and numerical significance of the individual parameter values on the efficacy of imatinib treatment of CML. In particular, we analyzed the effects of the numerical values of the model parameters on the behavior of the system, revealing critical threshold values that impact the ability of imatinib treatment to achieve remission and/or elimination. We also showed that STIs provide improvements to these critical values, categorizing this change as it relates to parameters inherent to either CML growth or immune response.
Citation: Dana Paquin, Lizzy Gross, Avery Stewart, Giovani Thai. Numerical analysis of critical parameter values for remission during imatinib treatment of chronic myelogenous leukemia[J]. Mathematical Biosciences and Engineering, 2025, 22(6): 1551-1571. doi: 10.3934/mbe.2025057
[1] | Rafael Martínez-Fonseca, Cruz Vargas-De-León, Ramón Reyes-Carreto, Flaviano Godínez-Jaimes . Bayesian analysis of the effect of exosomes in a mouse xenograft model of chronic myeloid leukemia. Mathematical Biosciences and Engineering, 2023, 20(11): 19504-19526. doi: 10.3934/mbe.2023864 |
[2] | R. A. Everett, Y. Zhao, K. B. Flores, Yang Kuang . Data and implication based comparison of two chronic myeloid leukemia models. Mathematical Biosciences and Engineering, 2013, 10(5&6): 1501-1518. doi: 10.3934/mbe.2013.10.1501 |
[3] | Natalia L. Komarova . Mathematical modeling of cyclic treatments of chronic myeloid leukemia. Mathematical Biosciences and Engineering, 2011, 8(2): 289-306. doi: 10.3934/mbe.2011.8.289 |
[4] | Elena Fimmel, Yury S. Semenov, Alexander S. Bratus . On optimal and suboptimal treatment strategies for a mathematical model of leukemia. Mathematical Biosciences and Engineering, 2013, 10(1): 151-165. doi: 10.3934/mbe.2013.10.151 |
[5] | Tianqi Song, Chuncheng Wang, Boping Tian . Mathematical models for within-host competition of malaria parasites. Mathematical Biosciences and Engineering, 2019, 16(6): 6623-6653. doi: 10.3934/mbe.2019330 |
[6] | Pedro José Gutiérrez-Diez, Jose Russo . Design of personalized cancer treatments by use of optimal control problems: The case of chronic myeloid leukemia. Mathematical Biosciences and Engineering, 2020, 17(5): 4773-4800. doi: 10.3934/mbe.2020261 |
[7] | Katrine O. Bangsgaard, Morten Andersen, Vibe Skov, Lasse Kjær, Hans C. Hasselbalch, Johnny T. Ottesen . Dynamics of competing heterogeneous clones in blood cancers explains multiple observations - a mathematical modeling approach. Mathematical Biosciences and Engineering, 2020, 17(6): 7645-7670. doi: 10.3934/mbe.2020389 |
[8] | Olga Vasilyeva, Tamer Oraby, Frithjof Lutscher . Aggregation and environmental transmission in chronic wasting disease. Mathematical Biosciences and Engineering, 2015, 12(1): 209-231. doi: 10.3934/mbe.2015.12.209 |
[9] | Hongfan Lu, Yuting Ding, Silin Gong, Shishi Wang . Mathematical modeling and dynamic analysis of SIQR model with delay for pandemic COVID-19. Mathematical Biosciences and Engineering, 2021, 18(4): 3197-3214. doi: 10.3934/mbe.2021159 |
[10] | Dan Zhu, Qinfang Qian . Optimal switching time control of the hyperbaric oxygen therapy for a chronic wound. Mathematical Biosciences and Engineering, 2019, 16(6): 8290-8308. doi: 10.3934/mbe.2019419 |
Chronic myelogenous leukemia (CML) is a cancer of the white blood cells that results from uncontrolled growth of myeloid cells in the bone marrow and the accumulation of these cells in the blood. The most common form of treatment for CML is imatinib, a tyrosine kinase inhibitor. Although imatinib is an effective treatment for CML and most patients treated with imatinib do attain some form of remission, imatinib does not completely eradicate all leukemia cells, and if treatment is stopped, all patients eventually relapse. Kim et al. constructed a system of delay differential equations to mathematically model the dynamics of anti-leukemia T-cell responses to CML during imatinib treatment, and demonstrated the usefulness of the mathematical model for studying novel treatment regimes to enhance imatinib therapy. Paquin et al. demonstrated numerically using this DDE model that strategic treatment interruptions (STIs) may have the potential to completely eradicate CML in certain cases. We conducted a comprehensive numerical study of the model parameters to identify the mathematical and numerical significance of the individual parameter values on the efficacy of imatinib treatment of CML. In particular, we analyzed the effects of the numerical values of the model parameters on the behavior of the system, revealing critical threshold values that impact the ability of imatinib treatment to achieve remission and/or elimination. We also showed that STIs provide improvements to these critical values, categorizing this change as it relates to parameters inherent to either CML growth or immune response.
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 106, 108, and 1010 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.
Remission Stage | Hematologic | Cytogenetic | Molecular | Total Elimination |
Concentration, k/µL | <1.67 | <1.67⋅10−2 | <1.67⋅10−4 | <10−10 |
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-α (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.
We begin with the one-dimensional case with a single delay r>0; initial data ϕ:[s−r,s]→R; and a given s∈R. We are focused on the solution to the following initial-value problem:
x′(t)=f(t,x(t),x(t−r)),x(t)=ϕ(t), s−r≤t≤s. | (2.1) |
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 fx(t,x,y) be continuous on R3, s∈R, and let ϕ:[s−r,s]→R be continuous. Then there exists σ>s and a unique solution of the initial-value problem (2.1) on [s−r,σ].
The goal is to extend the local existence of a solution from Theorem 1 to be defined for all t≥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≥s−r, then it will grow to infinity as t→σ, where σ comes from Theorem 1.
Theorem 2. Let f satisfy the hypothesis of Theorem 1, and let x:[s−r,σ)→R be the noncontinuable solution of the initial-value problem (2.1). If σ<∞, then
limt→σ−|x(t)|=∞. |
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∈Rn and f:R×Rn×Rn→Rn with almost no change in the proof. It also extends to multiple discrete delays r0<r1<⋯<rm where f=f(t,y(t),y(t−r0),y(t−r1),…y(t−rm)) 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:R×Rn+×Rn+→Rn satisfies the hypotheses of Theorem 1, Remark 1; and
∀i,t,∀x,y∈Rn+:xi=0⟹fi(t,x,y)≥0. |
If the initial data in (2.1) satisfy ϕ≥0, then the corresponding solution x(t) of 2.1 satisfies x(t)≥0 for all t≥s where defined.
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 y0(t),y1(t),y2(t), and y3(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 ry, and the growth rates from stem cells to progenitor, progenitor to differentiated, and differentiated to terminally differentiated are ay,by, and cy, respectively. The natural death rates at each stage are d0,d1,d2, and d3, respectively.
Under continuous imatinib treatment, the growth rates ay and by 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 z0,z1,z2, and z3, respectively. Similarly to the non-mutated system, the corresponding growth rates are az,bz,cz with stem cell regeneration rate rz. 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 sT, death rate dT, those lost to stimulated division in the presence of CML, and those gained via stimulated division after a time delay of τ days per division. We illustrate the movements of all cell populations in Figure 1.
The interaction term qCp(C,T)yi, or qCp(C,T)zi, represents the concentration of the leukemia sub-population killed after a successful encounter with a T-cell. Expanding p(C,T)=p0e−cnCkT helps to describe a successful interaction. Upon encountering a cancer cell with probability p0, 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 τ days and the cancer cell dying with probability qC. If n cells are ready to divide, they bring in 2n cells after nτ 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−cnC to represent immune down-regulation, where cn represents the decay rate of immune responsivity.
dy0dt=(ry(1−u)−d0)y0−qCp(C,T)y0,dy1dt=ayy0−d1y1−qCp(C,T)y1,dy2dt=byy1−d2y2−qCp(C,T)y2,dy3dt=cyy2−d3y3−qCp(C,T)y3, |
dz0dt=(rz−d0)z0+ryuy0−qCp(C,T)z0,dz1dt=azz0−d1z1−qCp(C,T)z1,dz2dt=bzz1−d2z2−qCp(C,T)z2,dz3dt=czz2−d3z3−qCp(C,T)z3, |
dTdt=sT−dTT−p(C,T)C+2np(Cnτ,Tnτ)qTCnτ, |
p(C,T)=p0e−cnCkT,C(t)=∑yi(t)+∑zi(t),Cnτ=C(t−nτ), Tnτ=T(t−nτ). |
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≥0.
Proof. Arrange the 5 non-mutated DDEs in the CML system as x′(t)=f(t,x(t),x(t−nτ)). Let x(t)=(y0(t),y1(t),y2(t),y3(t),T(t)) denote the order in which cell populations are written in the solution. Expanding p(C,T), each fi is defined as follows:
f1=(ry−d0)y0−qCp0kexp[−cn∑yi]Ty0,f2=ayy0−d1y1−qCp0kexp[−cn∑yi]Ty1,f3=byy1−d2y2−qCp0kexp[−cn∑yi]Ty2,f4=cyy2−d3y3−qCp0kexp[−cn∑yi]Ty3,f5=sT−dTT−p0kexp[−cn∑yi]T∑yi+2nqTp0kexp[−cn∑yi,nτ]Tnτ∑yi,nτ. |
It is clear that each fi consists of continuously differentiable functions (note the state variables y0,y1,y2,y3,T are assumed to be differentiable, and hence continuous). Setting our initial data ϕ to the initial conditions for the model, local existence and uniqueness are established by Theorem 1 (and Remark 1).
To verify positivity, suppose that y0=y1=y2=y3=T=0. Then f1=⋯=f4=0, and f5=sT>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≥0 with Theorem 2, i.e., that σ=∞. We assess bounds on each cell population. Since y0(t)≥0, then y′0(t)≤(ry−d0)y0. Using the variation of constants formula, we obtain
y0(t)≤y0(0)exp[(ry−d0)t]. |
With this,
y′1(t)≤ayy0(t)≤ayy0(0)exp[(ry−d0)t]⟹y1(t)≤y1(0)+ayy0(0)[exp[(ry−d0)t]−1ry−d0]=ayy0(0)(1d1−1ry−d0+exp[(ry−d0)t]ry−d0). |
Similarly:
y′2(t)≤byy1(t)≤byy1(0)+aybyy0(0)[exp[(ry−d0)t]−1ry−d0]⟹y2(t)≤y2(0)+byy1(0)t−aybyy0(0)ry−d0t+aybyy0(0)exp[(ry−d0)t](ry−d0)2=aybyy0(0)[1d1d2+(1d1−1ry−d0)t+exp[(ry−d0)t](ry−d0)2]. |
And:
y′3(t)≤cyy2(t)≤aybycyy0(0)[1d1d2+(1d1−1ry−d0)t+exp[(ry−d0)t](ry−d0)2]⟹y3(t)≤y3(0)+aybycyy0(0)[1d1d2t+12(1d1−1ry−d0)t2+exp[(ry−d0)t](ry−d0)3]=aybycyy0(0)[1d1d2d3+1d1d2t+12(1d1−1ry−d0)t2+exp[(ry−d0)t](ry−d0)3]. |
Finally:
T′(t)≤sT+2np0qTkexp[−cn∑yi,nτ]Tnτ∑yi,nτ≤sT⟹T(t)≤T(0)+sTt. |
By Theorem 2, these estimates tell us that the solution x(t) is defined for all t≥0.
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 ay and by 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.
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 1012 leukemia cells [14]. The three standard types of remission are a 100-fold difference from each other, corresponding to 1010, 108, and 106 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 α to have a critical threshold value toward achieving remission stage R, denoted as αR, if αR acts as an upper or lower bound of possible α values contributing to achieving a certain remission stage.
● If for any αi>α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 α larger than αR), then αR is an upper critical threshold value.
● If for any αi<αR, the time to remission stage R is not defined, then αR is a lower critical threshold value.
We focus on the remission threshold values for molecular remission and total elimination, denoted as αM and α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.
Parameter Category | Characterization, Post-STI Implication | Examples |
Immune-specific | Anti-leukemia immune response. "Healthy" values should be large; remission thresholds act as lower bounds. Thresholds under STIs should decrease; patients can attain remission with weaker immune systems. | k, p0, qC, qT, n |
CML-specific | CML growth, compromised immune response. "Healthy" patient values should be small; remission thresholds act as upper bounds. Thresholds under STIs should increase; patients can attain remission with stronger CML strains. | d1, d2, dT, ry, y0(0) |
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).
Parameter | Description | Base Value | Range |
λ | Fractional adjustment constant | 0.75 | 0.5 to 1 |
d0 | Stem cell death rate | 0.003 λ/day | ±25% |
d1 | Progenitor cell death rate | 0.008 λ | ±25% |
d2 | Differentiated cell death rate | 0.05 λ | ±25% |
d3 | Terminal cell death rate | λ | ±25% |
ry | Stem cell regeneration rate | 0.008/day | ±25% |
ay | Stem cell growth rate | 1.6 (without imatinib) | ±25% |
1.6/100 (with imatinib) | |||
by | Progenitor cell growth rate | 10 (without imatinib) | ±25% |
10/750 (with imatinib) | |||
cy | Differentiated cell growth rate | 100 | ±25% |
rz | Imatinib-resistant mutation: stem cell regeneration rate | 0.023/day | ±25% |
az | Imatinib-resistant mutation: stem cell growth rate | 1.6 | Same as ay |
bz | Imatinib-resistant mutation: progenitor cell growth rate | 10 | Same as by |
cz | Imatinib-resistant mutation: differentiated cell growth rate | 100 | Same as cy |
k | Kinetic (mixing) coefficient | 1 (k/μL)−1/day | ±25% |
p0 | Probability of T-cell engaging cancer cell | 0.8 | ±25% |
qC | Probability of cancer cell death from T-cell encounter | 0.75 | ±25% |
qT | Probability of T-cell surviving encounter with cancer cell | 0.5 | ±25% |
τ | Duration of one T-cell division | 1 day | 0.5 to 1 |
n | Average number of T-cell divisions | 2.2 | 1 to 3 |
dT | Anti-leukemia T-cell death rate | 0.0022/day | 1×10−3 to 1×10−2 |
sT | Anti-leukemia T-cell supply rate | 8.8×10−7 (k/μL)/day | 1×10−6 to 1×10−5 |
cn | Decay rate of immune responsivity | 7 (k/μL)−1 | 0 to 10 |
y0(0) | Initial leukemia stem cell concentration | 2.4×10−6 k/μL | Patient data [4] |
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.
To illustrate our process, consider the parameter p0, 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 p0=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 p0 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 p0 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 dT. 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).
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 p0 and qC) 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 dT and ry) 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, y0(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 y0(0), which is 51% lower than its non-STI counterpart. (See Figure 7 for the exact percentages.)
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., d1 and y0(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 y0(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 y0(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,z0(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 p0, for instance, a molecular threshold α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.
The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.
This research was generously supported by the William and Linda Frost Fund in the Cal Poly Bailey College of Science and Mathematics.
The authors declare that there is no conflict of interest.
[1] |
E. Jabbour, H. Kantarjian, Chronic myeloid leukemia: 2020 update on diagnosis, therapy and monitoring, Am. J. Hematol., 95 (2020), 691–709. https://doi.org/10.1002/ajh.25792 doi: 10.1002/ajh.25792
![]() |
[2] |
J. Barrett, Allogeneic stem cell transplantation for chronic myeloid leukemia, Semin. Hematol., 40 (2003), 59–71. https://doi.org/10.1016/S0037-1963(03)70043-2 doi: 10.1016/S0037-1963(03)70043-2
![]() |
[3] |
J. Cortes, M. Talpaz, S. O'Brien, D. Jones, R. Luthra, J. Shan, et al., Molecular responses in patients with chronic myelogenous leukemia in chronic phase treated with imatinib mesylate, Clin. Cancer Res., 11 (2005), 3425–3432. https://doi.org/10.1158/1078-0432.CCR-04-2139 doi: 10.1158/1078-0432.CCR-04-2139
![]() |
[4] |
F. Michor, T. Hughes, Y. Iwasa, S. Branford, N. Shah, C. Sawyers, et al., Dynamics of chronic myeloid leukemia, Nature, 435 (2005), 1267–1270. https://doi.org/10.1038/nature03669 doi: 10.1038/nature03669
![]() |
[5] |
P. S. Kim, P. P. Lee, D. Levy, Dynamics and potential impact of the immune response to chronic myelogenous leukemia, PLOS Comput. Biol., 4 (2008). https://doi.org/10.1371/journal.pcbi.1000095 doi: 10.1371/journal.pcbi.1000095
![]() |
[6] |
N. Komarova, D. Wodarz, Drug resistance in cancer: Principles of emergence and prevention, Proc. Natl. Acad. Sci., 102 (2005), 9714–9719. https://doi.org/10.1073/pnas.0501870102 doi: 10.1073/pnas.0501870102
![]() |
[7] |
I. Roeder, I. Glauche, Pathogenesis, treatment effects, and resistance dynamics in chronic myeloid leukemia¨Cinsights from mathematical model analyses, J. Mol. Med., 86 (2008), 17–27. https://doi.org/10.1007/s00109-007-0241-y doi: 10.1007/s00109-007-0241-y
![]() |
[8] |
I. Roeder, M. Horn, I. Glauche, A. Hochhaus, M. Mueller, M. Loeffler, Dynamic modeling of imatinib-treated chronic myeloid leukemia: Functional insights and clinical implications, Nat. Med., 12 (2006), 1181–1184. https://doi.org/10.1038/nm1487 doi: 10.1038/nm1487
![]() |
[9] |
B. Werner, D. Lutz, T. H. Brümmendorf, A. Traulsen, S. Balabanov, Dynamics of resistance development to imatinib under increasing selection pressure: A combination of mathematical models and in vitro data, PLoS One, 6 (2011), e28955. https://doi.org/10.1371/journal.pone.0028955 doi: 10.1371/journal.pone.0028955
![]() |
[10] |
S. Bunimovich-Mendrazitsky, B. Shklyar, Optimization of combined leukemia therapy by finite-dimensional optimal control modeling, J. Optim. Theory Appl., 175 (2017), 218–235. https://doi.org/10.1007/s10957-017-1161-9 doi: 10.1007/s10957-017-1161-9
![]() |
[11] |
S. Bunimovich-Mendrazitsky, N. Kronik, V. Vainstein, Optimization of interferon-Alpha and imatinib combination therapy for chronic myeloid leukemia: A modeling approach, Adv. Theory Simul., 2 (2019), 1800081. https://doi.org/10.1002/adts.201800081 doi: 10.1002/adts.201800081
![]() |
[12] |
D. Paquin, P. S. Kim, P. Lee, D. Levy, Strategic treatment interruptions during imatinib treatment of chronic myelogenous leukemia, Bull. Math. Biol., 73 (2011), 1082–1100. https://doi.org/10.1007/s11538-010-9553-0 doi: 10.1007/s11538-010-9553-0
![]() |
[13] | H. Smith, An Introduction to Delay Differential Equations with Applications to the Life Sciences, Springer, 2011. |
[14] |
B. Lowenberg, Minimal residual disease in chronic myeloid leukemia, N. Engl. J. Med., 349 (2003), 1399–1401. https://doi.org/10.1056/NEJMp038130 doi: 10.1056/NEJMp038130
![]() |
[15] |
D. Paquin, D. Sacco, J. Shamshoian, An analysis of strategic treatment interruptions during imatinib treatment of chronic myelogenous leukemia with imatinib-resistant mutations, Math. Biosci., 262 (2015), 117–124. https://doi.org/10.1016/j.mbs.2015.01.011 doi: 10.1016/j.mbs.2015.01.011
![]() |
[16] |
S. Niculescu, P. Kim, K. Gu, P. Lee, D. Levy, Stability crossing boundaries of delay systems modeing immune dynamics in leukemia, Discrete Contin. Dyn. Syst. - Ser. B, 13 (2010), 129–156. https://doi.org/10.3934/dcdsb.2010.13.129 doi: 10.3934/dcdsb.2010.13.129
![]() |
[17] |
A. Besse, G. Clapp, S. Bernard, F. Nicolini, D. Levy, T. Lepoutre, Stability analysis of a model of interaction between the immune system and cancer cells in chronic myelogenous leukemia, Bull. Math. Biol., 80 (2018), 1084–1110. https://doi.org/10.1007/s11538-017-0272-7 doi: 10.1007/s11538-017-0272-7
![]() |
[18] |
B. Cahlon, D. Schmidt, On stability of systems of delay differential equations, J. Comput. Appl. Math., 117 (2000), 137–158. https://doi.org/10.1016/S0377-0427(99)00337-4 doi: 10.1016/S0377-0427(99)00337-4
![]() |
[19] | C. T. Baker, F. A. Rihan, Sensitivity Analysis of Parameters in Modelling with Delay-Differential Equations, Manchester Centre for Computational Mathematics, 1999. |
[20] |
F. A. Rihan, Sensitivity analysis for dynamic systems with time-lags, J. Comput. Appl. Math., 151 (2003), 445–462. https://doi.org/10.1016/S0377-0427(02)00659-3 doi: 10.1016/S0377-0427(02)00659-3
![]() |
[21] |
S. Bunimovich-Mendrazitsky, L. Shaikhet, Stability analysis of delayed tumor-antigen-activatedImmune response in combined BCG and IL-2immunotherapy of bladder cancer, Processes, 8 (2020), 1564. https://doi.org/10.3390/PR8121564 doi: 10.3390/PR8121564
![]() |
Remission Stage | Hematologic | Cytogenetic | Molecular | Total Elimination |
Concentration, k/µL | <1.67 | <1.67⋅10−2 | <1.67⋅10−4 | <10−10 |
Parameter Category | Characterization, Post-STI Implication | Examples |
Immune-specific | Anti-leukemia immune response. "Healthy" values should be large; remission thresholds act as lower bounds. Thresholds under STIs should decrease; patients can attain remission with weaker immune systems. | k, p0, qC, qT, n |
CML-specific | CML growth, compromised immune response. "Healthy" patient values should be small; remission thresholds act as upper bounds. Thresholds under STIs should increase; patients can attain remission with stronger CML strains. | d1, d2, dT, ry, y0(0) |
Parameter | Description | Base Value | Range |
λ | Fractional adjustment constant | 0.75 | 0.5 to 1 |
d0 | Stem cell death rate | 0.003 λ/day | ±25% |
d1 | Progenitor cell death rate | 0.008 λ | ±25% |
d2 | Differentiated cell death rate | 0.05 λ | ±25% |
d3 | Terminal cell death rate | λ | ±25% |
ry | Stem cell regeneration rate | 0.008/day | ±25% |
ay | Stem cell growth rate | 1.6 (without imatinib) | ±25% |
1.6/100 (with imatinib) | |||
by | Progenitor cell growth rate | 10 (without imatinib) | ±25% |
10/750 (with imatinib) | |||
cy | Differentiated cell growth rate | 100 | ±25% |
rz | Imatinib-resistant mutation: stem cell regeneration rate | 0.023/day | ±25% |
az | Imatinib-resistant mutation: stem cell growth rate | 1.6 | Same as ay |
bz | Imatinib-resistant mutation: progenitor cell growth rate | 10 | Same as by |
cz | Imatinib-resistant mutation: differentiated cell growth rate | 100 | Same as cy |
k | Kinetic (mixing) coefficient | 1 (k/μL)−1/day | ±25% |
p0 | Probability of T-cell engaging cancer cell | 0.8 | ±25% |
qC | Probability of cancer cell death from T-cell encounter | 0.75 | ±25% |
qT | Probability of T-cell surviving encounter with cancer cell | 0.5 | ±25% |
τ | Duration of one T-cell division | 1 day | 0.5 to 1 |
n | Average number of T-cell divisions | 2.2 | 1 to 3 |
dT | Anti-leukemia T-cell death rate | 0.0022/day | 1×10−3 to 1×10−2 |
sT | Anti-leukemia T-cell supply rate | 8.8×10−7 (k/μL)/day | 1×10−6 to 1×10−5 |
cn | Decay rate of immune responsivity | 7 (k/μL)−1 | 0 to 10 |
y0(0) | Initial leukemia stem cell concentration | 2.4×10−6 k/μL | Patient data [4] |
Remission Stage | Hematologic | Cytogenetic | Molecular | Total Elimination |
Concentration, k/µL | <1.67 | <1.67⋅10−2 | <1.67⋅10−4 | <10−10 |
Parameter Category | Characterization, Post-STI Implication | Examples |
Immune-specific | Anti-leukemia immune response. "Healthy" values should be large; remission thresholds act as lower bounds. Thresholds under STIs should decrease; patients can attain remission with weaker immune systems. | k, p0, qC, qT, n |
CML-specific | CML growth, compromised immune response. "Healthy" patient values should be small; remission thresholds act as upper bounds. Thresholds under STIs should increase; patients can attain remission with stronger CML strains. | d1, d2, dT, ry, y0(0) |
Parameter | Description | Base Value | Range |
λ | Fractional adjustment constant | 0.75 | 0.5 to 1 |
d0 | Stem cell death rate | 0.003 λ/day | ±25% |
d1 | Progenitor cell death rate | 0.008 λ | ±25% |
d2 | Differentiated cell death rate | 0.05 λ | ±25% |
d3 | Terminal cell death rate | λ | ±25% |
ry | Stem cell regeneration rate | 0.008/day | ±25% |
ay | Stem cell growth rate | 1.6 (without imatinib) | ±25% |
1.6/100 (with imatinib) | |||
by | Progenitor cell growth rate | 10 (without imatinib) | ±25% |
10/750 (with imatinib) | |||
cy | Differentiated cell growth rate | 100 | ±25% |
rz | Imatinib-resistant mutation: stem cell regeneration rate | 0.023/day | ±25% |
az | Imatinib-resistant mutation: stem cell growth rate | 1.6 | Same as ay |
bz | Imatinib-resistant mutation: progenitor cell growth rate | 10 | Same as by |
cz | Imatinib-resistant mutation: differentiated cell growth rate | 100 | Same as cy |
k | Kinetic (mixing) coefficient | 1 (k/μL)−1/day | ±25% |
p0 | Probability of T-cell engaging cancer cell | 0.8 | ±25% |
qC | Probability of cancer cell death from T-cell encounter | 0.75 | ±25% |
qT | Probability of T-cell surviving encounter with cancer cell | 0.5 | ±25% |
τ | Duration of one T-cell division | 1 day | 0.5 to 1 |
n | Average number of T-cell divisions | 2.2 | 1 to 3 |
dT | Anti-leukemia T-cell death rate | 0.0022/day | 1×10−3 to 1×10−2 |
sT | Anti-leukemia T-cell supply rate | 8.8×10−7 (k/μL)/day | 1×10−6 to 1×10−5 |
cn | Decay rate of immune responsivity | 7 (k/μL)−1 | 0 to 10 |
y0(0) | Initial leukemia stem cell concentration | 2.4×10−6 k/μL | Patient data [4] |