A quantitative study on the role of TKI combined with Wnt/$\beta$-catenin signaling and IFN-$\alpha$ in the treatment of CML through deterministic and stochastic approaches

TL;DR

This study develops deterministic and stochastic models to analyze the combined effects of TKI, Wnt/$\beta$-catenin signaling, and IFN-$\alpha$ in treating CML, aiming to optimize therapy and predict remission likelihood.

Contribution

It introduces the first comprehensive mathematical model integrating all three therapies and incorporates stochastic analysis for CML treatment dynamics.

Findings

Patient-specific models enable tailored therapy dosages.

Dual and triple therapies potentially improve outcomes and reduce side effects.

Likelihood of remission can be predicted based on stem cell levels.

Abstract

We propose deterministic and stochastic models for studying the pharmacokinetics of chronic myeloid leukemia (CML), upon administration of IFN-$\alpha$ (the traditional treatment for CML), TKI (the current frontline medication for CML) and Wnt/$\beta$-catenin signaling (the state-of-the art therapeutic breakthrough for CML). To the best of our knowledge, no mathematical model incorporating all these three therapeutic protocols are available in literature. Further, this work introduces a stochastic approach in the study of CML dynamics. The key contributions of this work are: (1) Determination of the patient condition, contingent upon the patient specific model parameters, which leads to prediction of the appropriate patient specific therapeutic dosage. (2) Addressing the question of how the dual therapy of TKI and Wnt/$\beta$-catenin signaling or triple combination of all three, offers potentially improved therapeutic responses, particularly in terms of reduced side effects of TKI or IFN-$\alpha$. (3) Prediction of the likelihood of CML extinction/remission based on the level of CML stem cells at detection.

Figures

Tables

Table 1. Table 1 : Possible state changes during Δ ​ t Δ 𝑡 \Delta t .

i	State change $\left({\left(\Delta X\right)}_i\right)$	Probability $\left(p_i\right)$	Description
1	${(1,0,0)}^{\top }$	$r{C}_s\left(1-\frac{C_s}{K}\right)\Delta t$	Source of a CML stem cell
2	${(-1,0,0)}^{\top }$	$d_1^{\prime }{C}_s\Delta t$	Neutralizing by Wnt or natural death of a CML stem cell
3	${(0,1,0)}^{\top }$	$\mu C_s\Delta t$	Production of a mature CML cell
4	${(0,-1,0)}^{\top }$	$\left(\mu T+d_2^{\prime }\right)C_m\Delta t$	Neutralizing by CTL/ TKI or natural death of a mature CML cell
5	${(0,0,1)}^{\top }$	$\alpha C_{m}T\Delta t$	Development of a CTL cell
6	${(0,0,-1)}^{\top }$	$d_3^{\prime }T\Delta t$	Stimulating by IFN-$\alpha$ along with natural death of a CTL cell

Table 2. Table 2 : Parameter values for numerical results

Parameter	Value [Reference]	Unit
$r$	$1$	day-1
$K$	$41.667$ [15]	cells ml-1
$\gamma$	$3.647\times {10}^{-6}$ [15]	ml cell-1 day-1
$\mu$	$10$	day-1
$d_1$	$0.0375$ [15]	day-1
$d_2$	$0.0375$ [15]	day-1
$d_3$	$0.5$ [31]	day-1
$\alpha$	$0.00041$ [31]	ml cell-1 day-1
${\beta }_1$	$0.05$	day-1
${\beta }_2$	$0.07$ [31]	day-1
${\beta }_3$	$0.065$ [31]	day-1

Table 3. Table 3 : Qualitative summary of the therapeutic results.

	${\beta }_2\ne 0$	${\beta }_1\ne 0,{\beta }_2\ne 0$	${\beta }_2\ne 0,{\beta }_3\ne 0$	${\beta }_1\ne 0,{\beta }_2\ne 0,{\beta }_3\ne 0$
$C_s$	Unchanged	Decline	Unchanged	Decline
$C_m$	Decline at high dosage	Decline at high dosage	Decline	Fast decline
$T$	Decline	Decline	Dosage dependent decline/increase	Dosage dependent decline/increase

Equations

$$\frac{d{C}_s(t)}{dt}$$

$$\frac{d{C}_m(t)}{dt}$$

$$\frac{dT(t)}{dt}$$

$$\frac{d{C}_s(t)}{dt}=r^{\prime }{C}_s(t)\left(1-\frac{C_s(t)}{K^{\prime }}\right)-d_1{C}_s(t)-{\beta }_1{C}_s(t),$$

$$\frac{dC_{s}}{dt}\bigg{|}_{C_{s}=0}=0,\leavevmode\nobreak\ \frac{dC_{m}}{dt}\bigg{|}_{C_{m}=0}=\mu C_{s}\leavevmode\nobreak\ \text{and }\frac{dT}{dt}\bigg{|}_{T=0}=0.$$

$$\frac{dX(t)}{dt}$$

$$\mathcal{D}=\left\{(C_s(t),C_m(t),T(t))\in {\mathbb{R}}_{+}^3\text{leavevmode}:\text{leavevmode}0\le C_s(t)\le K,\text{leavevmode}0\le C_m(t)\le \frac{sK}{d_{23}^{\prime }},\text{leavevmode}0\le T(t)\le \frac{\alpha sK}{\gamma d_{23}^{\prime }}\right\}.$$

$$J=\left[\begin{array}{ccc}r\left(1-\frac{2C_s}{K}\right)-d_1^{\prime }-rT& 0& -r{C}_s\\ \mu & -d_2^{\prime }-\gamma T& -\gamma C_m\\ 0& \alpha T& \alpha C_m-d_3^{\prime }\end{array}\right].$$

$$a_0$$

$$a_1$$

$$a_2$$

$$P^{\ast \ast }=\frac{\alpha \mu K(r-d_1^{\prime })-r{d}_2{d}_3^{\prime }}{r{d}_3^{\prime }}.$$

$$\dot{X}=S\Phi ,$$

$$S=\left[\begin{array}{ccccc}r-d_1^{\prime }& 0& 0& -\frac{r}{K}& 0\\ \mu & -d_2^{\prime }& 0& 0& -\gamma \\ 0& 0& -d_3^{\prime }& 0& \alpha \end{array}\right]\text{leavevmode}\text{and}\text{leavevmode}\Phi =\left[\begin{array}{c}{C}_s\\ C_m\\ T\\ C_s^2\\ C_{m}T\end{array}\right]$$

$$B=S\cdot \text{diag}(\sqrt{\Phi })=\left[\begin{array}{ccccc}(r-d_1^{\prime })\sqrt{C_s}& 0& 0& -\frac{r}{K}{C}_s& 0\\ \mu \sqrt{C_s}& -d_2^{\prime }\sqrt{C_m}& 0& 0& -\gamma \sqrt{C_{m}T}\\ 0& 0& -d_3^{\prime }\sqrt{T}& 0& \alpha \sqrt{C_{m}T}\end{array}\right].$$

$$B{B}^{\top }=\left[\begin{array}{ccc}(r-d_1^{\prime }{)}^2{C}_s+\frac{r^2}{K^2}{C}_s^2& \mu (r-d_1^{\prime })C_s& 0\\ \mu (r-d_1^{\prime })C_s& {\mu }^2{C}_s+{d_2^{\prime }}^2{C}_m+{\gamma }^2{C}_{m}T& -\alpha \gamma C_{m}T\\ 0& -\alpha \gamma C_{m}T& {d_3^{\prime }}^{2}T+{\alpha }^2{C}_{m}T\end{array}\right].$$

$$H=\left[\begin{array}{ccc}{h}_{11}& 0& 0\\ h_{21}& h_{22}& 0\\ 0& h_{32}& h_{33}\end{array}\right],$$

$$h_{11}=\sqrt{(r-d_1^{\prime }{)}^2{C}_s+\frac{r^2}{K^2}{C}_s^2},\text{leavevmode}{h}_{21}=\frac{\mu (r-d_1^{\prime })C_s}{h_{11}},\text{leavevmode}{h}_{22}=\sqrt{{\mu }^2{C}_s+{d_2^{\prime }}^2{C}_m+{\gamma }^2{C}_{m}T-h_{21}^2},$$

$$h_{32}=\frac{-\alpha \gamma C_{m}T}{h_{22}},\text{leavevmode}{h}_{33}=\sqrt{{d_3^{\prime }}^{2}T+{\alpha }^2{C}_{m}T-h_{32}^2}.$$

$$dX(t)=f(X(t))dt+H(X(t))dW(t),$$

$$d{C}_s(t)$$

$$d{C}_m(t)$$

$$dT(t)$$

$$t_p=\inf \left\{t\in [0,t_e):\text{at least one element of}\text{leavevmode}\{C_s(t),C_m(t),T(t)\}\notin \left(\frac{1}{p},p\right)\right\}\text{leavevmode}\text{for}\text{leavevmode}p\ge p_0$$

$$P\{t_p\le t_T\}\ge ϵ,\text{leavevmode}\forall \text{leavevmode}p\ge p_1.$$

$$V(X(t))=\Big{(}C_{s}(t)+1-\ln C_{s}(t)\Big{)}+\frac{{d}^{\prime}_{1}}{\mu}\Big{(}C_{m}(t)+1-\ln C_{m}(t)\Big{)}+\frac{\gamma{d}^{\prime}_{1}}{\alpha\mu}\Big{(}T(t)+1-\ln T(t)\Big{)}.$$

$$dV(X(t))$$

$$LV(X(t))$$

$$c_1=d_1^{\prime }-r+\frac{d_1^{\prime }{d}_2^{\prime }}{\mu }+\frac{\gamma d_1^{\prime }{d}_3^{\prime }}{\alpha \mu }\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{and}\text{leavevmode}\text{leavevmode}\text{leavevmode}{c}_2=2\left(r+\frac{r}{K}\right)+\frac{2\gamma d_1^{\prime }}{\mu }.$$

$$\left(r+\frac{r}{K}\right)C_s(t)+\frac{\gamma d_1^{\prime }}{\mu }T(t)\le c_{2}V(X(t)).$$

$$dV(X(t))$$

$${\int }_0^{t_p\wedge t_1}dV(X(t))$$

$$V(X(t_p\wedge t_1))$$

$$EV(X(t_p\wedge t_1))$$

Full text

Keywords: CML Stem Cells; Mature CML Cells; CTL Response; Stochastic Model

1 Introduction

The complex physiological process of hematopoiesis involving the differentiation and regeneration of hematopoietic stem cells (HSC) arising in the bone marrow is extremely regulated and critical for maintaining the blood supply in the body and ensuring a balance between the differentiated cell lines of red blood cells (erythrocytes), white blood cells (leukocytes) and platelets [1, 2, 3]. Leukemia which is characterized through elevated blood cells (generally white blood cells) count is triggered by mutated blood cells violating the highly regulated process of hematopoiesis. It can be classified as either acute or chronic (contingent on the progression speed of the disease) and is either myeloid or lymphocytic (contingent on the maturity type of the cells) [2, 4]. A vast majority of cases of chronic myeloid leukemia (CML), a myeloproliferative condition is attributed to a reciprocal translocation of one chromosome 9 and one chromosome 22, resulting in one chromosome 9 being longer and one chromosome 22 being shorter, than their respective normal counterparts [2, 4, 5]. The occurrence of resulting chromosomal abnormality, called Philadelphia chromosome (denoted Ph) leads to the production of the BCR-ABL fusion oncogene. The BCR-ABL oncogene encodes for tyrosine kinase (a protein) that triggers the white blood cells to proliferate at an abnormal rate as compared to normal white blood cells [2, 4, 6]. The progression of CML can typically be observed in three phases, namely, the chronic phase, the accelerated phase and finally the blast crisis, which is akin to acute leukemia [7].

Traditionally, the standard therapy for CML includes interferon alfa-2a (IFN-$\alpha$) whose curative effectiveness was highly limited [8, 9]. The enhanced life expectancy of CML patients, upon administration of IFN-$\alpha$ is believed to be achieved through activation of the effector immune cells [10], specifically the stimulation of the effector T-cells [11]. At the turn of this century, a new line of CML treatment by way of tyrose kinase inhibitor (TKI), such as imatinib, dasatinib and nilotinib brought about a great leap in terms of therapeutic promise [2]. Because the TKIs were very target specific, they targeted the CML cells by preventing the activation of proteins [2]. TKIs limit the proliferation of the leukemic cells and causes their apoptosis, with minimal invasive effect on the healthy cells [6]. Imatinib’s success in terms of remission and patient survival was reported to be 90% for a 10-year period in patients receiving continuous treatment [12]. The remission is determined on the basis of the ratio of BCR-ABL to a constant transcript, namely, BCR, ABL or GUS [2]. Despite achieving remission and significantly increased life expectancy through TKIs, a truly complete cure remains elusive, due to non-responsiveness of quiescent leukemic stem cells (LSCs) to TKIs [13, 14]. Remissions along with concurrent existence of low-levels of leukemic cell levels believed to be due to allogenic immune response to tumor, is pivotal for long term cure [15]. However, a combination therapy of imatinib and IFN-$\alpha$ provides improved response upon administration of IFN-$\alpha$ after remission induced by imatinib has taken place [8]. The toxicity and the consequent side effects resulting from multi-drug combination can be addressed by optimal treatment scheduling [16]. In addition to IFN-$\alpha$, allogenic bone marrow transplant [2, 17] acted as the predominant treatment for CML, predating the introduction of TKIs. While the mechanism of both IFN-$\alpha$ and allogenic bone marrow transplant relies on activation of immune response, the practical usability of the latter is extremely limited. The effectiveness of IFN-$\alpha$ involves mechanism such as apoptosis and activation of immune cells [18]. As already observed, the TKIs despite achieving great therapeutic success, is ineffective against quiescent LSCs, while the IFN-$\alpha$ is effective against them. The complementary effectiveness of both renders them as a good choice for combination therapy [2]. The TKIs drive the CML cells into remission, leaving it up to the IFN-$\alpha$ to act against the residual cells. Further, IFN-$\alpha$ is capable of driving the quiescent LSCs to a point where they can be targeted by the TKIs, which otherwise would not have been possible [19, 20]. The results from clinical trials after cessation of therapeutic treatment of imatinib also gives credence to the role of immune response in CML [2].

The drug resistance of LSCs to TKIs leads to the consideration of alternative therapeutic protocols. The cells of the bone marrow microenvironment (BMM) regulate and protect LSCs from the effect of TKIs, resulting in the limited efficacy of TKIs against LSCs, which is responsible for disease persistence or relapse after discontinuation of the therapy [21, 22]. The Wnt which are a family of secreted glycoproteins plays an important role in the regulation of cell development, differentiation, proliferation and death [23]. The Wnt signaling can be modified by a Porcupine acyl transferase inhibitor. Such an inhibitor, WNT974, in combination with TKI regulates the Wnt/$\beta$-catenin signaling and potentially inhibits the proliferation of CML stem and progenitor cells [21, 22, 24]. The combination of Wnt/$\beta$-catenin signaling and TKI inhibits the BMM protection of LSCs from TKI targeted therapy and could enhance the therapeutic response in CML patients.

We discuss below some of the relevant mathematical models for CML dynamics as a precursor to the presentation of our proposed model. The discussion will mostly be centered around the dynamics of CML and the consequent immune response. The modeling of leukemia through a diffusion model was proposed in [25] and analyzed to study the spatial effect of leukemic cells. Accordingly, the populations of normal and leukemic cells were considered. The model is suggestive that there is eventual breakdown in the coexistence of normal and leukemic cells with eventual replacement of normal cells by the leukemic cells. The periodic CML (PCML) wherein oscillations are observed in the hematopoietic system is modeled in [3] with the key conclusion being that the levels of proliferating stem cells is critical in triggering the process of PCML. The dynamics of hematopoiesis was modeled by incorporating two different types of cells, namely, hematopoietic stem cells and differentiated cells for the normal as well as cancerous cells [1, 26, 27]. The global analysis suggested [1] that coexistence of normal and cancerous cells is not possible in the long run, with eventual convergence either to the safe equilibrium or the blast equilibrium. A model incorporating delay equations and three types of cell division along with clinical implication is dealt with in [26] wherein the evolution of different cell line populations is sought to be captured through the CML model studied. Cell divisions occurring during cell evolution in CML is considered in [27] along with the competition between normal and leukemic cells. A model involving the CML cancer cells and the naive T-cells as well as the CML specific effector T-cells was analyzed in [28]. Moore and Li [28] considered two types of immune cells, namely, the naive T-cells and the effector T-cells. The naive T-cell are activated only if they are CML-specific, whereas the effector T-cells are (apriori) CML-specific and can act immediately. In case a CML-specific naive T-cell is activated, it binds to a peptide major histocompatibility complex (MHC) pair and the presence of costimulators will eventually result in the proliferation of the T-cells. As week-long proliferation window is followed by differentiation of effector cells capable of producing immune response to the CML antigen [28]. A sampling based analysis led to the model prediction in which the Gompertzian growth rate as well as the natural death rate of CML cells (without altering the other model parameters) are the key parameters for the control of CML progression. The CML model due to Moore and Li [28] is revisited in [29] and a rigorous analysis of the global dynamics for the same is carried out. The authors established that the dynamics of the system is unstable around the tumor-free equilibrium and obtained the global stability conditions for the internal tumor equilibrium. A general model for CML to include normal, leukemic and resistant types of both stems cells and progenitor cells is considered in [30] and the global analysis for the disease free, healthy free and endemic equilibrium is discussed. Further, optimal control problems for several imatinib therapy scenarios and the suboptimal response in CML are studied. A mathematical model incorporating the combination treatment of imatinib to target the CML cells and IFN-$\alpha$ driven immunotherapy is discussed in [31]. The model involved the dynamics of the CML cells and the effector T-cells with a time varying delay of about a week for the therapeutic effect of IFN-$\alpha$ to kick in. The scheduling of the combination treatment of imatinib and IFN-$\alpha$ through an optimal control problem is presented in [8]. The model used involves the CML cells and the greater effector T-cells, akin to [31]. The key finding of the work was the duration of administration of imatinib and the cessation of treatment to be possible only vis-a-vis IFN-$\alpha$ therapeutic treatment. In Besse et al. [15], the authors reduced a model proposed in [32] (for which no theoretical analysis was provided) to one which accounts for cycling stem cells, mature leukemic cells and the immune cells. This simplification was driven by determination of the key characteristics of the CML dynamics reported in [32]. A bifurcation analysis based on the model system is suggestive of TKI treatment to be consistent with a treatment free low disease state being existent and stable.

The organization of the remainder of the paper is as follows. In Section 2, we present the deterministic model and its analysis. The formulation and analysis of the stochastic model are carried out in Section 3. This is followed by the numerical results and discussions in Section 4.

2 Deterministic Model

We propose the following model, involving the CML stem cells, the mature CML cells and the immune response in presence of combination therapy of Wnt/$\beta$-catenin signaling, TKIs and IFN-$\alpha$:

[TABLE]

Here $C_{s}(t)$ and $C_{m}(t)$ denote the concentrations of CML stem cells and mature CML cells, respectively, with $T(t)$ representing the number of effector T-cells cytotoxic to CML (CTL), all at time $t$. The first term in the equation for $C_{s}$ describes the logistic growth for the CML stem cells, with $r$ being the logistic growth rate and the constant $K$ being the carrying capacity, which represents the maximum possible concentration level of the CML stem cells that can be sustained. The natural death rate leading to the clearance of the CML stem cells is denoted by $d_{1}$. The inclusion of response to the targeted therapeutic protocol of Wnt/$\beta$-catenin signaling by a Porcupine acyl transferase inhibitor leading to the reduction in the level of CML stem cells is assumed to be taking place as a result of dosage $\beta_{1}$ of Wnt/$\beta$-catenin signaling. In the second equation which describes the dynamics of the mature CML cells, the production of the mature CML cells is assumed to happen as a result of differentiation of the CML stem cells at a rate $\mu$ accompanied by a natural cell death rate happening at a rate of $d_{2}$. The reduction of the mature CML cells is assumed to happen as a result of the influence of the CTL immune cells and also in response to the targeted treatment of TKIs (imatinib, nilotinib, or dasatinib). While the former is assumed to happen at a rate $\gamma$, the reduction of level of mature CML cells is assumed to take place as a result of TKI dosage of $\beta_{2}$. As CTL is an antigen-specific immune cell, we ignore its effect on CML stem cells. It is assumed that CTL does not have any effect on the CML stem cells, until they differentiate into mature ones. Finally, for the third equation, the production of the effector T-cells (CTLs) is driven by the immune response triggered by the presence of the mature CML cells, at a rate $\alpha$. The natural clearance of the CTLs is assumed to happen at a rate $d_{3}$. Further, the stimulation of CTL is also as a result of the therapeutic response to IFN-$\alpha$, whose administered dosage is assumed to be $\beta_{3}$. Note that, adding the decay term, $-\mu C_{s}$ (due to differentiation of the CML stem cells into mature ones), in the first equation of (2) results in the equation becoming

[TABLE]

where $\displaystyle{r^{\prime}=r-\mu}$ and $K^{\prime}=K\left(1-\frac{\mu}{r}\right)$. However, the form of the first equation of (2) was motivated by the models in [15, 32]. From the biological perspective, all the model parameters have to be positive, and accordingly, the initial condition for the system (2) is $\left(C_{s}(0),C_{m}(0),T(0)\right)\in\mathbb{R}_{+}^{3}$. We now define $d_{1}^{\prime}=d_{1}+\beta_{1}$, $d_{2}^{\prime}=d_{2}+\beta_{2}$ and $d_{3}^{\prime}=d_{3}-\beta_{3}$ and assume that $d_{3}^{\prime}>0$. We now show the non-negativity and boundedness of the solution to the model system (2).

Theorem 1**.**

The solution to the system (2) with initial condition in $\mathbb{R}_{+}^{3}$ exists in $\mathbb{R}_{+}^{3}$ and are ultimately bounded.

Proof.

From the theory of functional differential equations [33], it follows that the system (2) with the non-negative initial condition admits a unique solution $(C_{s}(t),C_{m}(t),T(t))$. Further, from the system (2), we obtain,

[TABLE]

Therefore, it follows that the solution $\left(C_{s},C_{m},T\right)$ starting from $\mathbb{R}_{+}^{3}$ exists in $\mathbb{R}_{+}^{3}$.

Now, from the first equation of the model system (2), we obtain $\displaystyle{\frac{dC_{s}(t)}{dt}\leq rC_{s}\left(1-\frac{C_{s}}{K}\right)}$. Therefore, $\displaystyle{\limsup\limits_{t\rightarrow\infty}C_{s}\leq K}$. We now introduce a new variable defined as $\displaystyle{X(t)=C_{m}+\frac{\gamma}{\alpha}T}$. From the last two equations of (2), we obtain

[TABLE]

Therefore, $\displaystyle{\limsup\limits_{t\rightarrow\infty}X\leq\frac{sK}{d_{23}^{\prime}}}$. Hence, $\displaystyle{\limsup\limits_{t\rightarrow\infty}C_{m}\leq\frac{sK}{d_{23}^{\prime}}}$ and $\displaystyle{\limsup\limits_{t\rightarrow\infty}T\leq\frac{\alpha sK}{\gamma d_{23}^{\prime}}}$. Thus the system (2) with the non-negative initial condition has a unique, non-negative and ultimately bounded solution in the positively invariant set

[TABLE]

∎

We now present the stability analysis for the steady states of the deterministic model (2), which admits three steady states, namely,

$\displaystyle{\widetilde{E}=(\widetilde{C}_{s},\widetilde{C}_{m},\widetilde{T})=(0,0,0)}$. 2. 2.

$\displaystyle{\bar{E}=(\bar{C}_{s},\bar{C}_{m},\bar{T})=\left(\frac{K(r-d_{1}^{\prime})}{r},\frac{\mu K(r-d_{1}^{\prime})}{rd_{2}^{\prime}},0\right)}$. 3. 3.

$\displaystyle{E^{*}=(C_{s}^{*},C_{m}^{*},T^{*})=\left(\frac{K(r-d_{1}^{\prime})}{r},\frac{d_{3}^{\prime}}{\alpha},\frac{\alpha\mu K(r-d_{1}^{\prime})-rd_{2}^{\prime}d_{3}^{\prime}}{r\gamma d_{3}^{\prime}}\right)}$.

We now define a threshold parameter expressed as $\displaystyle{P^{*}=\frac{\alpha\mu K(r-d_{1}^{\prime})}{rd_{2}^{\prime}d_{3}^{\prime}}}$. The equilibrium point $\bar{E}$ exists if $r>d_{1}^{\prime}$ and $E^{*}$ exists if $r>d_{1}^{\prime}$ and $P^{*}>1$. The Jacobian for the model system is given by

[TABLE]

At $\widetilde{E}$, the eigenvalues of the Jacobian $J$ are $x=-d_{2}^{\prime}<0,\leavevmode\nobreak\ x=-d_{3}^{\prime}<0$ and $x=r-d_{1}^{\prime}$, which implies that $\widetilde{E}$ is stable when $r<d_{1}^{\prime}$ and unstable when $r>d_{1}^{\prime}$. For $\bar{E}$, we have $r>d_{1}^{\prime}$ and the eigenvalues of the Jacobian $J$ are $x=-d_{2}^{\prime}<0,\leavevmode\nobreak\ x=-r+d_{1}^{\prime}<0$ and $x=\frac{\alpha\mu k(r-d_{1}^{\prime})-rd_{2}^{\prime}d_{3}^{\prime}}{rd_{2}^{\prime}}<0$ which holds only if $P^{*}<1$. The eigenvalues of the Jacobian $J$ at $E^{*}$ are obtained by solving the characteristic equation given by $x^{3}+a_{2}x^{2}+a_{1}x+a_{0}=0$, where,

[TABLE]

Since $C_{s}^{*},C_{m}^{*},T^{*}>0$ if $P^{*}>1$, therefore $a_{0},a_{1},a_{2}>0$ and $a_{1}a_{2}-a_{0}=\frac{rC_{s}^{*}}{K}\left(d_{2}^{\prime}+\gamma T^{*}\right)\left(d_{2}^{\prime}+\gamma T^{*}+\frac{rC_{s}^{*}}{K}\right)>0$ if $P^{*}>1$. Thus, by the Routh-Hurwitz criteria, all the three eigenvalues of the Jacobian $J$ at $E^{*}$ are either negative or have negative real parts if $P_{*}>1$.

In summary, we can state the following theorem:

Theorem 2**.**

For the model system (2), the following results hold:

1. (a)

The disease free equilibrium $\widetilde{E}$ is stable when $r<d_{1}^{\prime}$ and unstable when $r>d_{1}^{\prime}$. 2. (b)

The disease but CTL free equilibrium $\bar{E}$ is locally asymptotically stable when $r>d_{1}^{\prime}$ with $P^{*}<1$ and unstable when $P^{*}>1$. 3. (c)

The disease and CTL coexistent equilibrium $E^{*}$ is locally asymptotically stable when $r>d_{1}^{\prime}$ with $P^{*}>1$.

The assumption $r>d_{1}^{\prime}$ leads to the disease free equilibrium $\widetilde{E}$ being unstable. However, in case of $r<d_{1}^{\prime}$, the disease free equilibrium $\widetilde{E}$ would be stable which is suggestive of remission or even cure as a result of the combined rate of natural death and Wnt/$\beta$-catenin signaling exceeding the growth rate of the CML stem cells. Further, in case a patient is in state $E^{*}$ (which corresponds to a more severe case than $\bar{E}$), an appropriate choice of dosage of $\beta_{2}$ can drive the system from $P^{*}>1$ to $P^{*}<1$, resulting in stability at $\bar{E}$. We now define

[TABLE]

The condition for the most severe disease equilibrium $E^{*}$ is stable if $P^{*}>1$ i.e., $\beta_{2}<P^{**}$. If $\beta_{2}$ is increased such that $\beta_{2}>P^{**}$ which is equivalent to $P^{*}<1$, then the system converges to the less severe disease equilibrium $\bar{E}$. Finally, if $\beta_{1}$ is pushed beyond $r-d_{1}$, then the CML may be sent to remission or cure.

3 Stochastic Model

Although the deterministic model analyzed in Section 2, provides useful insight on the dynamics of CML progression, it does not offer any description of the probability of extinction of the disease. Accordingly, the random fluctuations in the physiological factors, the complex mechanism of the in-host immune system and the random process of cell growth/death motivates the introduction of stochasticity in the model. Therefore, in order to introduce randomness and uncertainly in the model description, we derive a stochastic differential equation (SDE) model, based on the state changes during interaction among the populations of CML stem cells, mature CML cells and CTLs themselves (considered in formulation of the deterministic model (2)), and we formulate the SDE model following the approach in [34, 35].

Let $C_{s}$, $C_{m}$ and $T$ represent three random variables denoted in vector form as $X=(C_{s},C_{m},T)^{\top}$. The changes in these variables in a small time interval $\Delta t$ is given by $\Delta X=X(t+\Delta t)-X(t)$. Based on the mechanism of the deterministic model, the system (2) can be written as

[TABLE]

where

[TABLE]

represent the coefficient matrix and interaction strategy among the populations, respectively. Therefore the diffusion matrix is given by

[TABLE]

Hence

[TABLE]

In order to reduce the number of Wiener processes, we find an alternative diffusion matrix $H$ having less number of columns such that $BB^{\top}=HH^{\top}$. By the Cholesky decomposition method, we compute

[TABLE]

where

[TABLE]

Hence the Itô SDE model corresponding to the system (2) is

[TABLE]

where $W(t)=(W_{1}(t),W_{2}(t),W_{3}(t))^{\top}$ is a vector of three independent Wiener processes. Thus the SDE model can be explicitly expressed as follows

[TABLE]

where $d^{\prime}_{1}=d_{1}+\beta_{1},\leavevmode\nobreak\ d^{\prime}_{2}=d_{2}+\beta_{2}$ and $d^{\prime}_{3}=d_{3}-\beta_{3}$.

We will show that the solution of the stochastic system (3) is non-negative and bounded on $n$-th moment as well as stochastically bounded, in the sense of probability.

Theorem 3**.**

For any initial condition $\displaystyle{X_{0}=(C_{s}(0),C_{m}(0),T(0))\in\mathbb{R}^{3}_{+}}$, there exists a unique solution $\displaystyle{X(t)=}$ $\displaystyle{\left(C_{s}(t),C_{m}(t),T(t)\right)\in\mathbb{R}^{3}_{+}}$ to the stochastic system (3) for all $t\geq 0$ almost surely (a.s.) (i.e., with probability $1$), where $\displaystyle{\mathbb{R}^{3}_{+}=}$ $\displaystyle{\left\{(x_{1},x_{2},x_{3}):x_{i}>0,\leavevmode\nobreak\ i=1,2,3\right\}}$.

Proof.

It is easily seen that the coefficients of the system (3) satisfy the local Lipschitz condition [36, 37]. Therefore, for any initial condition $X_{0}\in\mathbb{R}^{3}_{+}$, there exists a unique local solution $X(t)\in\mathbb{R}^{3}_{+}$ to the system (3) for $\displaystyle{t\in[0,t_{e}),\leavevmode\nobreak\ t_{e}>0}$. In order to prove the solution is global, i.e., $\displaystyle{X(t)\in\mathbb{R}^{3}_{+}}$ a.s. for all $t\geq 0$, we have to show that $t_{e}=\infty$ a.s. We choose a sufficiently large number $p_{0}\geq 0$ such that each of $\displaystyle{C_{s}(0),C_{m}(0)\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ T(0)}$ lies in $\displaystyle{\left[\frac{1}{p_{0}},p_{0}\right]}$. We now define

[TABLE]

and $\displaystyle{\inf\phi=\infty}$, where $\phi$ is an empty set. We observe that $t_{p}$ is increasing as $p$ increases and denote $\displaystyle{t_{\infty}=\lim_{p\to\infty}t_{p}}$. This implies $\displaystyle{t_{\infty}\leq t_{e}}$ a.s. Therefore, in order to show that $t_{e}=\infty$ a.s., we prove that $t_{\infty}=\infty$ a.s. If possible, let $t_{\infty}<\infty$ a.s. Then there exists two constants $t_{T}>0$ and $\epsilon\in(0,1)$ such that $P\{t_{\infty}\leq t_{T}\}>\epsilon$. Hence there exists an integer $p_{1}\geq p_{0}$ such that

[TABLE]

We define a $\mathbb{C}^{2}$-function $V:\mathbb{R}^{3}_{+}\to\mathbb{R}_{+}$ by

[TABLE]

The function $V$ is positive definite, since $(x+1-\ln x)>0,\leavevmode\nobreak\ \forall\leavevmode\nobreak\ x>0$. Using Itô formula, we calculate the differential of $V$ along the solution trajectories of the stochastic system (3) as follows

[TABLE]

where

[TABLE]

Let

[TABLE]

Then using the relation $x\leq 2(x+1-\ln x),\leavevmode\nobreak\ \forall\leavevmode\nobreak\ x>0$, we obtain

[TABLE]

Hence

[TABLE]

where $c_{3}=\max\{c_{1},c_{2}\}$. By taking $t_{1}\leq t_{T}$, we obtain

[TABLE]

where $t_{p}\land t_{1}=\min\{t_{p},t_{1}\}$. This results in

[TABLE]

Taking expectation on both sides and using Fubini’s theorem as well as properties of Itô integral, we obtain

[TABLE]

Using the Gronwall’s inequality, we obtain

[TABLE]

where c_{4}=\Big{(}V(X_{0})+c_{3}t_{T}\Big{)}e^{c_{3}t_{T}}. Let $\Omega_{p}=\{t_{p}\leq t_{T}\}$ for $p\geq p_{1}$. Therefore, using (3.2), it follows that $P(\Omega_{p})\geq\epsilon$. Note that for every $\omega\in\Omega_{p}$, there is at least one of $C_{s}(t_{p},\omega),\leavevmode\nobreak\ C_{m}(t_{p},\omega)$ and $T(t_{p},\omega)$ which equal(s) either $p$ or $\displaystyle{\frac{1}{p}}$ and hence

[TABLE]

where

[TABLE]

Therefore (3.2) and (3.3) give

[TABLE]

where $1_{\Omega_{p}}$ is the indicator function of $\Omega_{p}$ and defined by

[TABLE]

Taking $p\to\infty$ leads to the contradiction $\infty>c_{4}=\infty$. Hence we can conclude that $t_{\infty}=\infty$ a.s. This completes the proof of Theorem 3. ∎

Theorem 4**.**

Let $\displaystyle{X(t)=\left(C_{s}(t),C_{m}(t),T(t)\right)\in\mathbb{R}^{3}_{+}}$ be a solution of the stochastic system (3) with an initial condition $\displaystyle{X_{0}=(C_{s}(0),C_{m}(0),T(0))\in\mathbb{R}^{3}_{+}}$. Then

[TABLE]

where $M_{i}(n),\leavevmode\nobreak\ i=1,2,3$ are finite positive quantities and $n\geq 1$.

Proof.

Let $U_{1}(C_{s}(t))=[C_{s}(t)]^{n},\leavevmode\nobreak\ n\geq 1$. Using Itô formula, we obtain from (3),

[TABLE]

Integrating both sides from [math] to $t$, we obtain

[TABLE]

Taking expectation on both sides and using Fubini’s theorem as well as properties of Itô integral, we obtain

[TABLE]

Differentiating both sides with respect to $t$, we obtain

[TABLE]

By Hölder’s inequality, it follows that

[TABLE]

Using the above inequality, one can obtain

[TABLE]

Letting $E[C^{n}_{s}(t)]\geq 1$ leads to

[TABLE]

Using the technique of solving Bernoulli differential equations and theory of differential inequalities, we obtain

[TABLE]

where

[TABLE]

Further, let $U_{2}(C_{m}(t))=[C_{m}(t)]^{n},\leavevmode\nobreak\ n\geq 1$. Using Itô formula, we obtain from (3),

[TABLE]

Taking expectation after integration from [math] to $t$, we obtain

[TABLE]

After differentiation with respect to $t$, the above equation becomes

[TABLE]

From inequality (3.4), there exists $L>0$ such that

[TABLE]

Therefore,

[TABLE]

Using Hölder’s inequality, we obtain

[TABLE]

Letting $E[C^{n}_{m}(t)]\geq 1$ leads to

[TABLE]

For $\displaystyle{n\leq 1+\frac{2}{\gamma}}$, we obtain

[TABLE]

Therefore

[TABLE]

where

[TABLE]

Finally, let $U_{3}(T(t))={[T(t)]}^{n},\leavevmode\nobreak\ n\geq 1$. Using Itô formula, we obtain from (3),

[TABLE]

Proceeding on similar lines, we obtain

[TABLE]

From inequality (3.7), there exists $\bar{L}>0$ such that

[TABLE]

Therefore,

[TABLE]

Using Hölder’s inequality, we obtain

[TABLE]

Assuming $d^{\prime}_{3}>\alpha\bar{L}$, it follows that

[TABLE]

where

[TABLE]

Hence the theorem follows. ∎

Theorem 5**.**

Let $\displaystyle{X(t)=\left(C_{s}(t),C_{m}(t),T(t)\right)\in\mathbb{R}^{3}_{+}}$ be a solution of the stochastic system (3) with an initial condition $\displaystyle{X_{0}=(C_{s}(0),C_{m}(0),T(0))\in\mathbb{R}^{3}_{+}}$. Then for every $\epsilon>0$, there is a constant $L$ for which

[TABLE]

where $||X(t)||=\sqrt{C^{2}_{s}(t)+C^{2}_{m}(t)+T^{2}(t)}$. That is, the solution $X(t)$ is stochastically bounded (or bounded in probability).

Proof.

From Theorem 4, it follows that

[TABLE]

Putting $n=2$, we obtain

[TABLE]

This implies that for $M>0$, there is a $t^{\prime}>0$ such that

[TABLE]

Since $E\left[||X(t)||^{2}\right]$ is continuous function of $t$, there exists $\bar{M}>0$ such that

[TABLE]

Assuming $B=\max\{{M},\leavevmode\nobreak\ \bar{M}\}$, we obtain

[TABLE]

Hence

[TABLE]

Using Tchebychev’s inequality, it follows for $L>0$ that

[TABLE]

which implies

[TABLE]

Taking $\epsilon=\frac{\bar{B}}{L^{2}}$, we obtain the following result

[TABLE]

∎

Theorem 6**.**

Suppose the following conditions are satisfied:

[TABLE]

where

[TABLE]

with $\eta=\min\{r^{2},\leavevmode\nobreak\ 2rK^{-1}\}$. Then the disease free equilibrium (trivial solution) of the stochastic system (3) is almost surely exponentially stable.

Proof.

Let

[TABLE]

Using Itô formula, we obtain

[TABLE]

where $\eta=\min\{r^{2},\leavevmode\nobreak\ 2rK^{-1}\}$. Note that the following expression

[TABLE]

can be expressed as $Y^{\top}PY$, where

[TABLE]

The matrix $P$ is negative definite under the following conditions being satisfied:

[TABLE]

where

[TABLE]

Let $\lambda_{\text{max}}$ be the largest (negative) eigenvalue of $P$. Then

[TABLE]

Therefore

[TABLE]

Note that

[TABLE]

Using inequality (3.12), we obtain

[TABLE]

Integrating and using the fact that

[TABLE]

we obtain

[TABLE]

i.e.,

[TABLE]

Hence $C_{s}(t)\to 0$, $C_{m}(t)\to 0$ and $T(t)\to 0$ a.s. as $t\to\infty.$ This completes the proof of Theorem 6. ∎

We now determine the probability of CML extinction in the disease persistent stage. In order to ascertain this probability, we identify the possible state changes during the time interval $\Delta t$ tabulated in Table 1.

We derive the offspring probability generating functions (PGFs) using the theory of multitype continuous-time branching process approximation [38, 39]:

The offspring PGF for $C_{s}$ (corresponding to $i=1,2,3$ in Table 1) is

[TABLE] 2. 2.

The offspring PGF for $C_{m}$ (corresponding to $i=4,5$ in Table 1) is

[TABLE]

Consequently, the expectation matrix for the PGFs is

[TABLE]

If $(q_{1},q_{2})\in(0,1)\times(0,1)$ is a fixed point, then

[TABLE]

This implies $\displaystyle{q_{1}=1,\leavevmode\nobreak\ \frac{d_{1}^{\prime}}{r}}$ and $\displaystyle{q_{2}=1}$. For the disease persistent stage, $r>d_{1}^{\prime}$, we consider $\displaystyle{q_{1}=\frac{d_{1}^{\prime}}{r}}$. Using Galton-Watson branching theory [40] , it follows that the probability of CML extinction is,

[TABLE]

4 Numerical Results

In this section, we dwell upon the numerical illustration of the results obtained for both the deterministic and the stochastic CML models. The computations were run using MATLAB® for a set of parameter values given in Table 2, except for the cases with some particular parameter values (which will be specified in the discussion). Accordingly, we take the initial condition to be $\displaystyle{\left(C_{s}(0),C_{m}(0),T(0)\right)=\left(20,500,1000\right)}$, for the deterministic as well as the stochastic model simulations, with the number of sample paths for the latter (generated using the Euler-Maruyama method [41]) being $5000$. In order to numerically show the stability of the equilibria, we choose, two different values of the parameter $\mu$, namely, $\mu=1$ and $\mu=10$, to generate various scenarios. The set of all other parameter values chosen from Table 2 with $\mu=1$ results in $P^{*}=0.3333(<1)$ which gives the existence of the equilibrium $\displaystyle{\bar{E}}=\left(38.0211,353.6850,0\right)$. Similarly, for $\mu=10$ we get $P^{*}=3.3335(>1)$ which gives the existence of the equilibrium $\displaystyle{E^{*}=\left(38.0211,1060.9756,68785.3432\right)}$.

We first examine the dynamics for the case $P^{*}<1$ corresponding to the case when $\bar{E}$ exists and is stable. The level of the CML stem cells increases rapidly, eventually stabilizing at the peak level of $C_{s}=38.0211$ shortly after the initial time (Figure 1(a)). As for the mature CML cells, the level gradually decreases and eventually stabilizes at the minimum level of $C_{m}=353.6850$ (Figure 1(b)). However, it takes more time to reach the corresponding stabilized level as compared to the case of CML stem cells. Further, the level of $T$ gradually decreases, eventually reaching zero level beyond which time there is cessation of the CTL immune response (Figure 1(c)). We conclude that the reduction in the level of the mature CML cells, consequently, results in gradual decline with eventual non-response of the CTLs. Consequently, the patient stabilizes at the less severe disease equilibrium $\bar{E}$. We now examine the dynamics for the case $P^{*}>1$ when $E^{*}$ exists and is stable. Similar to the case of $P^{*}<1$, the population level of the CML stem cells increases very fast and in the long run stabilizes at the peak level of $C_{s}=38.0211$ (Figure 2(a)). The level of the mature CML cells in this case increases very fast until the CTL immune response is stimulated to an efficient level. As a result of the affect of CTL immune response, the level of the mature CML cells keeps dropping. Consequently, due to less concentration of the mature CML cells, the generation of CTL immune response also decreases after elapse of some time and eventually the levels of the mature CML cells and the CTL immune response stabilize at the level of $C_{m}=1060.9756$ and $T=68785.3432$, respectively. These results are illustrated in Figure 2(b) and Figure 2(c), respectively. The simulations for the case of $P^{*}>1$ suggests that sufficient level of mature CML cells results in continued CTL immune response, eventually sending the patient to the most severe disease equilibrium.

We now assess the pharmacokinetic impact resulting from administration of Wnt/$\beta$-catenin signaling ($\beta_{1}$), TKI ($\beta_{2}$) and IFN-$\alpha$ ($\beta_{3}$) for several scenarios. Firstly, we consider the impact of monotherapy of TKI, followed by dual therapy of TKI with Wnt/$\beta$-catenin and IFN-$\alpha$, individually, with eventual consideration of all three. In order to present the effect of administration of TKI ($\beta_{2}$) in absence of Wnt/$\beta$-catenin signaling ($\beta_{1}$) and IFN-$\alpha$ ($\beta_{3}$), we ran the simulations with $\beta_{1}=0$ and $\beta_{3}=0$, for five different values of $\beta_{2}$, namely, [math], $0.05$, $0.1$, $0.5$ and $0.8$, having corresponding $P^{*}$ values of $8.7695$, $3.7584$, $2.3917$, $0.6118$ and $0.3927$. The simulations for CML stem cells, mature CML cells and the CTL immune response are presented in Figure 3. From Figure 3(a), it is observed that variation of $\beta_{2}$ has no impact at all on the dynamics of the CML stem cells. However, with an increase in the values of $\beta_{2}$, the values of the mature CML cells converge to ${C}_{m}^{*}$ for the first three values of $\beta_{2}$ (which corresponds to $P^{*}>1$), but converges to the lower value $\bar{C}_{m}$, with an increase for the other two values of $\beta_{2}$ (which corresponds to $P^{*}<1$) (Figure 3(b)). For the CTL immune response, one observes a decrease corresponding to increase in the values of $\beta_{2}$ (Figure 3(c)). For the last two values of $\beta_{2}$, the CTL immune response is eventually not sustained. From the simulations, we observe that the increase in $\beta_{2}$, upto a certain level does not affect the eventual stabilized level $E^{*}$ of the mature CML cells. Upon further increase of $\beta_{2}$, the impact is that the level of the mature CML cells decreases and eventually the progression stabilizes at a level $\bar{E}$. Thus, in this case, the administration of high dosage is essential for achieving desired therapeutic response in terms of reduction of mature CML cells, although there is no impact on CML stem cells at all. This result supports the observations made in [13, 14]

The impact of administration of dual therapy of Wnt/$\beta$-catenin signaling ($\beta_{1}$) and TKI ($\beta_{2}$), in absence of IFN-$\alpha$ ($\beta_{3}$), are presented in Figure 4. The set of $\left(\beta_{1},\beta_{2}\right)$ values chosen are $(0,0)$, $(0.5,0.05)$, $(0.1,0.2)$, $(0.1,0.5)$, $(0.3,0.5)$ and $(0.3,0.6)$, with the corresponding values of $P^{*}$ being $8.7695$, $1.8059$, $1.2408$, $0.5483$, $0.4211$ and $0.3551$. In this case, the values of CML stem cells show a decreasing trend resulting from the corresponding increase in the values of $\beta_{1}$ only (irrespective of $\beta_{2}$) (Figure 4(a)). For the values $(0.1,0.2)$, $(0.1,0.5)$ of $(\beta_{1},\beta_{2})$, the trajectories for CML stem cells converge along the same path as in the case where $\beta_{1}$ is left unchanged. A similar scenario is observed in case of $(0.3,0.5)$ and $(0.3,0.6)$. For the dynamics of mature CML cells, we observe that, for the case with $P^{*}>1$, there is an eventual convergence to a certain identical level of mature CML cells, while the convergent level gradually gets lowered with a gradual decrease in $P^{*}$ when $P^{*}<1$ (Figure 4(b)). The decrease in value of $P^{*}$ results in the corresponding decrease in the levels of CTL immune response and even converges to zero level for the cases where $P^{*}<1$ (Figure 4(c)). The simulation suggests that this dual therapy is not effective on mature CML cells in case of $P^{*}>1$. However, this therapy has an impact on CML stem cells, which was observed in [21, 22, 24].

We now observe the pattern of dual therapy of TKI ($\beta_{2}$) and IFN-$\alpha$ ($\beta_{3}$), in absence of Wnt/$\beta$-catenin signaling ($\beta_{1}$), in Figure 5. The set of $\left(\beta_{2},\beta_{3}\right)$ values chosen are $(0,0)$, $(0.1,0.1)$, $(0.3,0.1)$, $(0.3,0.2)$, $(0.8,0.05)$ and $(0.2,0.4)$, with the corresponding values of $P^{*}$ being $8.7695$, $2.9896$, $1.2179$, $1.6239$, $0.4363$ and $6.9233$. Interestingly, the absence of $\beta_{1}$ results in the dynamics of CML stem cells remaining unchanged irrespective of the various combinations of $\beta_{2}$ and $\beta_{3}$ being applicable (Figure 5(a)). From Figure 5(b), it is observed for the values $(0.1,0.1)$ and $(0.3,0.1)$ of $\left(\beta_{2},\beta_{3}\right)$ that, the concentration of mature CML cells does not change in terms of its stabilized level due to unchanged IFN-$\alpha$ dosage. Further, in case of the value $(0.8,0.05)$ of $(\beta_{2},\beta_{3})$, the level of mature CML cells reaches a low level, in spite of $\beta_{3}$ being very less, since higher value of $\beta_{2}$ results in $P^{*}<1$. For the value $(0.2,0.4)$, we observe that the mature CML cell stabilizes at the least level (Figure 5(b)), as CTL response increases to very high level in this case (Figure 5(c)). This shows that a low dosage of $\beta_{2}$ can be more effective in combination with a relatively high dosage of $\beta_{3}$. The simulations for the values $(0.3,0.1)$ and $(0.3,0.2)$ show that the level of mature CML cells gets reduced (Figure 5(b)), due to increase in dosage of IFN-$\alpha$, resulting in stimulation of CTL immune response (Figure 5(c)), a phenomenon which was not observed in case of monotherapy of TKI. Thus, the increase in IFN-$\alpha$ results in stimulation of CTL immune response, which is responsible for eventual decrease in the level of mature CML cells. However, this dual therapy does not have any effect on the CML stem cell population.

Finally, we analyze the clinical implications of the combination therapy of Wnt/$\beta$-catenin signaling ($\beta_{1}$), TKI ($\beta_{2}$) and IFN-$\alpha$ ($\beta_{3}$), as illustrated in Figure 6. The sets of $\left(\beta_{1},\beta_{2},\beta_{3}\right)$ values are taken to be $(0,0,0)$, $(0.1,0.2,0)$, $(0,0.2,0.1)$, $(0.1,0.2,0.1)$, $(0.6,0.5,0.4)$ and $(0.3,0.4,0.1)$. Accordingly, the corresponding values of $P^{*}$ are $8.7695$, $1.2408$, $1.7308$, $1.5509$, $1.1521$ and $0.6467$. The pattern of eventual convergence of the CML stem cells consistently decreases with the corresponding increase in the value of Wnt/$\beta$-catenin signaling (Figure 6(a)). The sets $(0.1,0.2,0)$ and $(0.1,0.2,0.1)$ of $\left(\beta_{1},\beta_{2},\beta_{3}\right)$ show that a small addition of IFN-$\alpha$ results in a sufficient therapeutic effectiveness so as to achieve reduction in the levels of mature CML cells (Figure 6(b)), with an increase of CTL immune response (Figure 6(c)). The combinations of $(0,0.2,0.1)$ and $(0.1,0.2,0.1)$ shows that the change in Wnt/$\beta$-catenin signaling (without changing TKI or IFN-$\alpha$) does not have any effect on the stabilized level of mature CML cells, but the convergence to the stable level is achieved at different times for different values of Wnt/$\beta$-catenin signaling (Figure 6(b)). The administration of the targeted therapy of Wnt/$\beta$-catenin results in decreasing level of CML stem cells (Figure 6(a)) and consequently a decline in the level of CTL immune response is observed (Figure 6(c)). Thus the combination dosage $(0.1,0.2,0.1)$ has a significant impact on the level of both CML stem cells and mature CML cells (observed as compared to the dosage $(0,0,0)$). A choice of high combination dosage of $(0.6,0.5,0.4)$ results in an adequate therapeutic response in terms of reduction in the concentration of CML stem cells as well as mature CML cells. Lastly, in case of $(0.3,0.4,0.1)$, a relatively high amount of Wnt/$\beta$-catenin signaling and TKI in combination with less IFN-$\alpha$ corresponding to $P^{*}<1$ results in decrease being observed for both the CML stem cells as well the mature CML cells. Consequently, the CTL immune response eventually goes to zero. Thus, a suitable combination of Wnt/$\beta$-catenin signaling, TKI and IFN-$\alpha$ can be very effective in the suppression and reduction of both the types of CML cells as well as in the stimulation of CTL response.

We qualitatively summarize the numerical results depicting the dynamics corresponding to the effect of therapeutic protocols, in Table 3.

We now perform the numerical simulation for the stochastic model (3) and present the histograms, in order to understand the probability distribution, at $t=50$ (chosen for illustrative purpose) based on $5000$ independent simulations, for the CML stem cells, mature CML cells and CTL immune response. An illustrative sample path along with the mean sample path, mean sample path minus standard deviation and mean sample path plus standard deviation is presented in Figures 7(a), 7(b) and 7(c), for the CML stem cells, mature CML cells and CTL immune response, respectively. The mean ($\mu_{X_{i}}$) and standard deviation ($\sigma_{X_{i}}$) at $t=50$, resulting from the $5000$ simulations are

[TABLE]

and

[TABLE]

respectively. The standard deviation for the case of CML stem cells is relatively less than that of mature CML cells and CTL immune response, with respect to the population size. Moreover, the stochastic mean, $\left(\mu_{C_{s}},\mu_{C_{m}},\mu_{T}\right)$ is very close to its deterministic equilibrium value $E^{*}$. The stochastic fluctuation has a lesser impact on the dynamics of CML stem cells and it has a greater impact on CTL population.

Further, we present the probability histogram plots in Figure 8, and observe that the probability distribution of CML stem cells and mature CML cells fit the normal curve. The mean and variance sets used to plot the normal distribution in Figures 8(a), 8(b) and 8(c) are $\left(\mu_{C_{s}},\sigma_{C_{s}}^{2}\right)$, $\left(\mu_{C_{m}},\sigma_{C_{m}}^{2}\right)$ and $\left(\mu_{T},\sigma_{T}^{2}\right)$, respectively.

Recall that the probability of CML extinction is dependent on both Wnt/$\beta$-catenin signaling ($\beta_{1}$) as well as the initial level of the CML stem cells ($C_{s}(0)$) which is depicted in Figure 9. Figure 9 illustrates the probability of disease extinction depending on various dosages of Wnt/$\beta$-catenin signaling as well as various initial levels of the CML stem cells. It is observed from Figure 9(a) that the probability of CML extinction is inversely proportional to the initial value of CML stem cells and it becomes very less (near to zero) if the initial CML stem cells exceeds $2$. Figure 9(b) suggests that the probability of disease extinction can be improved by increasing the dosage of Wnt/$\beta$-catenin signaling. It can be seen that the probability of CML extinction is greater when the initial level of the CML stem cells is less and Wnt/$\beta$-catenin signaling is high. The tabulated values of the parameters (Table 2) gives $P^{*}=3.26051$ and $q_{1}=0.1075$. Then the probability of CML extinction is $\displaystyle{q=\left(0.1075\right)^{C_{s}(0)}}$. As an illustration, we consider $C_{s}(0)=20$. Then $\displaystyle{q=\left(0.1075\right)^{C_{s}(0)}\approx 4.2478\times 10^{-20}}$, which shows that the probability of cure is almost non-existent and the prognosis is very poor. If we want the probability of CML extinction to be greater than a desired value (say $q_{e})$, then the dosage of the required Wnt/$\beta$-catenin signaling is as follows,

[TABLE]

5 Conclusion

This work focuses on the modeling and quantitative analysis for the disease progression dynamics of CML in consideration with CTL immune response. The pharmacokinetics in presence of the traditional therapeutic protocol of IFN-$\alpha$ and TKI is extended to include Wnt/$\beta$-catenin signaling. The deterministic model is proposed and analyzed. The stability of the severe and less severe equilibria (in case of the natural progression rate of CML stem cells exceeding the sum of the rate of natural clearance and the efficacy of TKI dosage) is presented with subject to a condition on a threshold parameter $P^{*}$. However, the deterministic model does not provide any insight into the likelihood of remission. This is addressed by the introduction of stochastic model which encapsulates this critical aspect of long term prognosis of CML remission and the probability of extinction of the disease. The well-posedness of both the deterministic model and stochastic model (in probability sense) is proved by showing the existence, uniqueness, boundedness and non-negativity of the solution. The stochastic mean solution is very close to its deterministic severe equilibrium value $E^{*}$. The stochastic noise is relatively more impactful on the population of mature CML cells and CTL cells. The stochastic model also showed the almost surely exponential stability of the disease free equilibrium under certain conditions. Further, the model predicts improved therapeutic response upon administration of Wnt/$\beta$-catenin signaling in addition to the IFN-$\alpha$ and TKI, particularly due to the resulting reduced side effects of the latter two. The results obtained demonstrates that the combination therapy is most successful in terms of the reduction of CML stem cells and mature CML cells with the stimulation of CTL immune response.