Development of the Analytical Model for the Assessment of the Efficiencies of Different Therapeutic Modalities in Leukemia

Different therapeutic regimes like chemotherapy, cytokine based immunotherapy and/or stem cell transplantation are generally be the suggested options for treatment of leukemia. However, the selection of the best therapeutic strategy for an individual patient remains uncertain till now. Several analytical models are proposed for each of the therapeutic strategies; however, no single analytical model is available which can make a comparative assessment regarding the long-term therapeutic efficacy among these strategies. This, in turn, may limit the clinical application of such analytical models. To address this issue we have developed an analytical model that can make a comparative assessment among these therapeutic regimes along with their different combinations. The model also has the flexibility to test their efficacies in different types of leukemia. Contrary to the previous models here we have incorporated several other important constraints like chemotherapeutic toxicity and chronic effects of HLA incompatibility that may limit therapeutic success. The simulation studies with the developed model indicate towards the immune potentiation as the determining factor to bring the leukemia free survival. It also reveals that the time optimization for immune boosting is another important crucial factor for providing leukemia free survival.

Keywords: Leukemia; Difference delay equation; Cachexia; Chemotherapy; HLA matching

Several mathematical models have been proposed to study the complex dynamical behavior of the hematopoiesis. These models have attempted to reveal the pathogenesis of different hematological diseases (HD) including hematological malignancies (HM). Amongst, the first work was targeted with acute myelogenous leukemias (AML) where systems of delay ordinary difference equations (DODE) are considered for modeling the hematopoietic dynamics in different compartments of hematopoietic lineages and the effect of chemotherapeutic intervention [1-3]. This modeling strategy has also been used in explaining the behaviour of periodic chronic myelogeneous leukemia and cyclical neutropenia diseases [4,5]. Similar modeling strategy was further extended for designing of long term targeted chemotherapeutic strategies (imatinib) in chronic myelogenous leukemias (CML) [6] and the effect of autologous vaccination procedure by T-cells [7]. Using a system of nonlinear delay differential equations, further work suggests that stability of the system is sensitive to the T cell related parameters like T cell stimulation and/or supply of T cell into the system [8].

However, models without delay have also been proposed by several workers [9-11]. These sorts of modeling approaches are helpful for the control and stability analysis of the system. Presently available analytical models for hematopoiesis consisting with multicompartmental models and the incorporation of ordinary differential delay equations makes it feasible to study the replication potential of the hematopoietic cells in several compartments as well as the chemotherapeutic efficiency, whereas other models without delay consideration are formulated for the simplicity of stability and control analysis.

For achieving better tolerance of the hematopoietic cell transplantation, hematopoietic stem cell (HSC) transplantation is suggested. Stem cell transplantation (SCT), an emerging therapeutic procedure in hematology involves a complex procedure,however, little has been done with the mathematical modeling, in spite of rich medical and clinical literatures are available in this field. In this regard works of DeConde et al, 2005 [12] can be mentioned where delay differential equations are used to explain the time evolution of six populations after allo-SCT in chronic myeloid leukemia and concluded that remission depends on the allo-reactive T cell immune function [7]. Previous modeling approaches were further modified to incorporate the drug resistant mutation in CML and suggest for combination therapy as a means of long-term remission in CML undergoing imatinib treatment [13]. Further works suggests that relapse is due to presence of dormant stem cells and/or symmetric division of stem cell acquire for the development of drug resistant phenotype [14,15].

However, allo-SCT transplantation is a complex therapeutic procedure as it may require intermittent intervention at different time points due to non-availability of ideal HLA matching in clinical scenarios. Hence, a flexible model with different input variables would be the ideal for the assessment of the potentiality of different therapeutic modalities for varying circumstances. In several instances of leukemia, partial remission is achieved by the application of myeloablative chemotherapy (MAC) and/or radiotherapy.Allogenic hematopoietic stem cell transplantation (allo-HSCT) is therefore suggested for achieving long-term leukemia free survival (LFS). Donor lymphocyte alloreactivity against the host hematopoietic system translates powerful graft-versus-leukemia (GVL) reactions [16]. However if the transplanted cells have the HLA incompatibility with the host then a severe graft-versushost diseases (GVHD) occurs which favors the opportunistic infection within the host. This has the significance of transplanted related mortality. Allo-HSC transplantation is regarded as the double edge sword as it protects the host from leukemia relapse and simultaneously the death of the host.

Though it is difficult to get, however complete remission through transplantation is possible with the ideal HLA matching with the syngenic (monogygotic twin of the patient) donor cells. Allo-HSC transplantation has two folds benefits – on one hand the conditioning regimen consisting of MAC destroys the hematopoietic process of patient which in turn reduce the tumor load in the patient, and following procedure of allo-HSC transplantation exerts a long lasting immunologic GVL by the donor derived T-lymphocytes which helps to remove the residual diseases of the patient after the MAC [17].

Now-a-days, allo-HSC transplantation (HSCT) is done with three sources – bone marrow (BM), peripheral blood (PBL) and umbilical cord blood (UCB). In recent time majority of cases HSCT is performed with PBL as the sources of HSC. For harvesting of HSC from the PBL it is needed to mobilize the HSC to the peripheral blood of the donor by the application of GMCSF or low dose cyclophosphamide. Allo-HSC of UCB origin is the better option compared to the adult HSC because of its greater tolerance within the host body due to low expression of the HLA on the stem cell surface [18].

Attempt has been made to study the post transplant immune responses in CML [2]; however, the management of the post transplant complications is still a challenge for different types of leukemias and yet to be determined. Though cancer vaccination strategy has been designed [7], however, in designing this model several parametric estimations like kinetic coefficient, probability of cancer cell death is not a very clinically feasible approach while one considers the individual leukemia cases. Therefore it is needed to develop a generalized model for different leukemias that could be more clinically compatible in the assessment of different therapeutic procedures including transplantation and its complications. Here we propose a simplified but clinically compatible analytical model for the assessment of the outcome of chemotherapeutic strategies, immune potentiation, HSCT and post-HSCT complications.

In present model we have considered a difference delay equation based three compartmental system model [19] having 7 cell types - the stem cell compartment having hematopoietic stem cell (S), progenitor compartment having three progenitor cells of three lineages namely erythroid (P1), leukocytoid (P2) and megakaryocytoid (P3) and matured cell compartment having mature cells of three lineages namely RBC (B1), WBC (B2) and platelet (THB) (B3). In each compartment three different category cells – normal cells (g), drug sensitive cells (s) and drug resistive cells (r) are considered. Hence in first compartment three cell types (S_g, S_s and S_r), in second compartment total nine cell types (P1_g, P2_g, P3_g, P1_s, P2_s, P3_s, P1_r, P2_r and P3_r) and in third compartment total nine cell types (B1_g, B2_g, B3_g, B1_s, B2_s, B3_s, B1_r, B2_r and B3_r) are considered. Therefore in present system model we have considered total 21 types of cells (Figure 1).

The equation (1) is used as a general form to represent the cellular dynamics of each cell type of any three compartment of the hematopoietic system at discrete time interval (k).

In the above equations (1) – (3) Nx(k) represents the number of cell on k-th day. N represents the cell type and suffix x represents the category of that concerned cell type. N_xp(k-dk) is the precursor cell type from which the present cell type was differentiated dk time ago. Hence dk is the delay time required for cell maturation by the action of different cytokines, cellular signaling molecules and transcription factors. Each cell type grows exponentially with its own multiplication rate (m_Nx) and decays with its own apoptosis rate (a_Nx). Thus, for normal stem cell these variables are denoted with (m_Sg) and (a_Sg). Similarly for drug sensitive P1 type cell, these variables are represented by (m_P1s) and (a_P1s). These notational schemes are followed for other cell types and are represented in Table 1A-C. It is assumed that multiplication rate of each cell type at the matured cells compartment is zero.

For the cells which are differentiable such as stem cell (S), erythroid (P1), leukocytoid (P2) and megakaryocytoid (P3), they will have c1 = 1 in their system equation. For the system equation of RBC, WBC and platelet, c1 = 0 as they do not differentiate. When the effect of differentiation is absent, that is, when k is less than the delay time (dk) required for the formation of the cell type to the down-stream cell type of the same lineage, c2 = 0. Otherwise c2 will be a positive quantity. In the present model, during the stage of differentiation from stem cell to progenitor cell stage, value of c2 = 1/3 (as it is assumed that stem cell is distributed equally to three lineages after differentiation). However, in the system equation for the cells of the subsequent compartment, c2 =1 (as each progenitor cell type is differentiated into a single lineage only).

Again N_xdr is the differentiation rate of any cell type and N_xpdr is the differentiation rate of its concerned precursor cell type. For example, in the system equation of stem cells has only N_xdr and for sensitive category stem cell, N_xdr is represented as S_sdr. Similarly, for drug resistive category of P1 lineage cell type, N_xdr = P1_rdr and N_xpdr = S_xdr; for P3 lineage cell type of drug sensitive category, N_xdr = P3_sdr and N_xpdr = S_sdr and differentiation rates for other cell types are listed in Table 1A-C.

Now, if k<dk (delay time required for the formation of a progenitor cell from its previous compartment cell of the same lineage) the system equation of concerned differentiated cell types can be represented through the following generalized equation:

In equations (3) - (5), the value of f is different for different cell types. The above mentioned generalized equations may be used to develop the system model in form of a difference delay equation [19] which is mentioned below.

While considering the conversion from concerned cell type (i.e., N_x) to other two cell categories of same type concerned conversion rates represented by C_{N xy} and C_{N xz} are subtracted from the diagonal elements of A matrix. Again a_(i-2,j-1), a_(i-2,j), a_(i-1,j-2), a_(i-1,j), a_(i,j-2) and a_(i,j-1) elements in matrix A where i = j = 3, 6, 9, 12, 15, 18, 21 can be modified with the corresponding conversion rates to incorporate conversion from other two categories to concerned cell type. For example, putting i = j = 3 in (i,j), we get six matrix elements of A namely a_(1,2), a_(1,3), a_(2,1), a_(2,3), a_(3,1), a_(3,2); where a_(1,2) and a_(1,3) are the conversion rates from S_r and S_s to S_g i.e., C_Ssg and C_Srg respectively (Table 1A). For other category cells of stem cell type the other elements become operative similarly; for P1, i = j = 6 and likewise.

The number of different normal cells is maintained through physiological homeostatic mechanism. In the mature cell compartment if the number is decreased below a certain level then positive feedback system causes the respective progenitor cells to increase its number. Similarly, if the number increases, then it produces negative feedback. Functionally within the physiological system, if RBC counts decreases, then reticulocytes (autocrine feedback) and erythroblasts are stimulated by erythropoietin (EPO) to form mature RBC. Similarly, when platelet counts decreases, the normal count (number) is maintained by megakaryocytes and megakaryoblast. Likewise when WBC cell count (say T-lymphocyte) decreases, its' number is increased by IL-2 or other cytokines by stimulating the corresponding precursor cells or lymphoblasts. The negative feedback occurs within the leukocyte lineage when the number increases. For lymphoid lineage a feedback system is also considered from mature cell level upto the HSC level. Autocrine regulation is considered for blast cell level and HSC by the effect of cytokine; for example, regulation by EPO, IL-2, thrombopoietin and stem cell factor at erythroblast, lymphoblast, megakaryoblast and HSC level respectively. In the model, both autocrine and inter-compartmental feedback mechanisms are considered. Inter-compartmental feedback mechanism is being effective by modulating the differentiation rate and apoptosis rate (Figure 2).

It is considered that under normal condition, host immune system is capable of removing leukemic cells [20-22]. This has been incorporated by modifying equation (7a) into equation (7b):

where the functional effective (immunity related) killing of the malignant cells (kill_eff) by B2_g becomes effective depending upon the minimum cell number of B2_g (Min _B2g) and is operated by kill_im and CF. The assumed efficiency of immune related killing (kill_eff) by the host immune cells are in the order of B2_s>B2_r=P2_s>P2_r>S_s>S_r. It is assumed that initially there is total (functional) immune suppression i.e., kill_im = 0. Application of any immune-modulatory agent is required to make the immune related killing effective. It is also assumed that CF is a time varying variable that saturates up to 60% immunity level of the normal population after application of immune-modulatory agent for a certain period of time and reaching this level will be dependent upon the applied dose of that immune-modulatory agent.

About 50% of the cancer patients suffer from a progressive weight loss due to atropy of adipose tissue and skeletal muscle. This phenomenon is called cachexia. This is more often associated with preterminal cancer patients. Cachexia may be present in the early stages of tumor growth before any signs or symptoms of malignancy [23]. Recent time several tumor related soluble factors are identified those are responsible for cancer cachexia [24,25]. Though cancer cachexia is not a term frequently used in leukemia; however in recent time by ultrastructural analysis it has been shown that there is some morphological de-arrangement in the red cell membrane (phenotype effect) in different leukemias including chronic and myelodysplastic syndrome [26]. Such phenomenon is known as cancer cachexia and is considered in the model.

Leukmic blast cells may produce much more cachexia effect. It is considered that resistant variety of malignant cells will produce the most detrimental effect compared to other cell types. Two types of toxicity are considered here – phenotypic and genotypic. The gradations of phenotype effect are implemented by three factors MTox1, MTox2 and MTox3 and genotype effect implemented by factor Tox. Considered gradation of toxicity intensities are in the following order: MTox1 > MTox2 > MTox3 > Tox (Figure 3). Hence a_(i,j) of equation (7b) will be modified as:

Again these three factors are multiplied with three different factors to implement the severity of killing intensity in different compartments. It has been assumed that the order of killing intensity is decreasing from matured cell compartment to stem cell level. Hence MTox1, MTox2 and MTox3 are again multiplied with another three factors ch_MS (= 0.02) for stem cell, ch_MP (= 1.5) for progenitor cell compartment and ch_MB (= 5) for matured cell compartment. It is also assumed that ch_MP (= 0) when P2_r (k) + P2_s (k) <10,000.

Conventionally, in clinical practice patients are treated with a myeloablative (MYL) chemotherapeutic regime. This has been incorporated in the model (Figure 4). However, before applying the drug strategy, cells are allowed to grow for certain period of time (say, 30 days). On the consecutive day myeloablative chemotherapy is applied. MYL is introduced to the system model as a specific subtractive term with a fixed percentage to the existing number of a specific cell population at discrete interval of time. The percentage of cell killing is dependent on the drug sensitivity of a particular cell type.

In the system model there is a provision of changing the drug dose, the number of cycles of drug application and interval between two successive drug applications. So in the system equations a_(i,j) of equation (7c) have been further modified as:

where, d_myl(k) is the amount of drug present within the system denoted either the amount of applied drug (drug_myldose) on the day of drug application or drug retention (drug_mylret) on the subsequent days as a fraction of the previous day's applied drug (it is assumed that in each corresponding day of drug application a certain amount of drug is cleared from the system);

Drug application is started on the day k_mylsd and n_myl is the number of cycles and t_myld is the interval of two subsequent drug applications in days and C_myl is the ON/OFF switch for drug application. Myl_Nx is the myeloablative drug sensitivity of the concerned cell N_x . The drug sensitivity Myl_Nx can be varied accordingly from different lineages and/or cell types (Table 2A).

Model for the supportive therapy to MYL: High dose of MYL therapy produces reduction in RBC and platelet level in the patient. In clinical practice, this has been tackled through intermittent transfusion of exogenous RBC and platelet. These sorts of therapeutic management are also incorporated in the model. So equation (7) is modified as

On the day of MYL application, minimal level (a set value) of RBC (RBC_llm) and platelets (Platelet_llm) are checked and if found lower than that level, then RBC and/or platelet are transfused. If both the values are sufficient then drug will be delivered. If either value is insufficient, then, as a precautionary measure, an upper lower level of the other (that is, if RBC_llm is less then Platelet_ulm or if Platelet_llm is less then RBC_ulm) will be checked and if both are found below, then transfusion decision of both cell types will be made accordingly. Drug is not delivered on the day of any type of transfusion is applied. And in that case MYL drug application will be delayed by Sft days from k_mylsd and after which subsequent days will be applied. This strategy prevents the unnecessary delay in subsequent drug application. The small delay period after transfusion is considered due to restoration of patient's physiological condition for tolerating the subsequent chemotherapy load; however, model has the flexibility to change or avoid this gap period by making Sft=0. Again, transfused cells are represented through following relations -

T_RBC and T_PLATELET represent amount of transfused RBC and platelets respectively. E1 and E2 may be fraction or whole number. Myl_TPlatelet and Myl_TRBC are the MYL drug sensitivities of T_RBC and T_PLATELET respectively. Application strategy of blood transfusion is demonstrated in Figure 5.

With a mismatched HLA between donor and host, the transplanted donor cells will also be recognized as foreign cell by the host immune system and vice versa. This produces Graft versus host disease (GVHD) (Figure 6). It is to be noted here that generally MYL therapy (as mentioned earlier) is applied as the conditioning regime before stem cell transplantation.

Transplanted cells' behavior within recipient's system: Transplanted cells' behavior can be represented by modifying the equation (8). In equation (8) each of the matrixes will be updated as follows. Hence,

where a_(i,j), b¹_(i,j), b_2(i,j), b³_(i,j), b⁴_(i,j) representing the elements of [28×28] matrixes A, B¹, B², B³, B⁴ respectively and RBC_{transfuse(i,j)}, PLATELET_{transfuse(i,j)} representing the elements of [28×1] matrixes RBC_transfuse and PLATELET_transfuse respectively. In equation (8a) the effect of cachexia on the donor cells are also included. In equation (8a), the effect of MYL drug on the donor cell is ignored as at the time of transplantation the quantity of residual drug in the system has been considered to be negligible. In HLA mismatched case, transplanted donor cells are killed due to GVHD reaction is incorporated into the system model by a factor kill_b2g→D. Transplantation with 100% HLA matching G = 0; otherwise G>0, that means with the increase in HLA mismatch level the value of G (Graft versus host diseases, GVHD) will increase. The killing intensity kill_mul is operative for GVHD effect that can be different for different cell types of donor.

Host cells' behavior towards transplanted cells: Model has been formulated with an assumption that donor stem cell (S_e) transplantation is made when chemotherapeutic drug level within the system is minimum say 10% of the last administered drug dose. However, this can be changed according to the choice of the investigator. Alloreactive lymphocytes (B2_e) are developed [along with other lineages i.e., erythrocytes (B1_e) and platelets (B3_e)] from transplanted stem cells lead to an immune response causing Graft-versus-tumor (GVT) (or GVL) [kill_b2e→hm] and/or GVHD [kill_b2e→hg] effect. So a_(i,j) represented in equation (8a) is modified further as:

In equation (8a-1) hm stands for host leukemic cells (S_s, S_r, P1_s, P1_r, P2_s, P2_r, P3_s, P3_r, B1_s, B1_r, B2_s, B2_r, B3_s, B3_r) and hg stands host normal cells (i.e., S_g, P1_g, P2_g, P3_g, B1_g, B2_g and B3_g). The GVHD reaction is operated into the system by modifying the kill factors kill_b2g→D and kill_b2e→hg. The killing intensity kill_mul is operative for GVHD and GVL effect that can be different for different cell types of host.

Nonlinear effect of GVHD: It is assumed that with the increase in the degree of HLA mismatch, effect of GVHD also increases simultaneously. The effect of GVHD (kill_b2e→hg) works on all the normal cells of the recipient. GVHD produces several aberrant unknown but cascade of biochemical reactions within the recipient that gradually increases with the progress of time.

It is assumed that the effect of GVHD will be operative when minimum cell numbers of B2_e is developed from transplanted HSC. As a result with the progress of time GVHD may be operative nonlinearly depending on the cell numbers of B2_e. Killing efficiency of B2_e to other hematopoietic cells of donor will be governed with a multiplying factor effect_procell which is governed though the following conditions.

In the above relation, effect_p1g, effect_p2g and effect_p3g have some constant values (Table 4). The equation (8a-1) will be modified further as follows:

In HLA mismatch the GVHD related killing intensity by host cells (kill_mu1) can be different for different cell types of donor cells. In the model it is assumed that host lymphocytes are unable to kill the malignant cells and at the time of transplantation the effects of cancer cachexia on donor as well as host stem cells are nil, as donor cells are transplanted in a condition when the tumor load is minimum. Different kill factors that have been considered in GVHD are represented in Table 4.

In our model we have considered the effect of chronic GVHD (cGVHD). Clinically the effect of cGVHD is assessed by different qualitative based clinical observations of changes in liver, gut, lungs [27]. Mathematically this can be represented through a marker (MR) that has been considered in the model as a time varying squared function of the number of leukocytes of donor origin. This function will vary according to the effect of different grades of mismatch (G) or GVHD. This is represented by the following relation, where Cal_GVHD is a time varying calibration factor.

Due to immune suppression leukemic patients are susceptible to infections [28]. Such susceptibility becomes more pronounced after MYL therapy. Clinically this could be circumvented by the application of cytokines (CYT) like interferon, GM-CSF and/ or IL-2 [29-33]. Application of cytokines can also augment the number of immunocytes. This, in turn, may also kill the residual malignant cells. Such therapeutic strategy has incorporated in the model equation. This cytokine application scheme can also be operated single or in combination of other therapeutic strategy. Application strategy of MYL/Chemotherapy with transfusion followed by cytokine therapy is shown in Figure 7.

It is assumed that these cytokines will affect the multiplication rate of the progenitor cells of the leukocytic lineages, which in turn, also affect the differentiation process for the formation of mature immunocytes. To introduce this sorts of therapeutic procedure, the a_(i,j) of equation (8a-2) will be modified as follows:

k_CYTsd, n_CYT, t_CYTd, m_CYT and C_CYT indicate the starting day of cytokine (drug) application, the number of cycles and the interval of two subsequent cytokine drug applications (in days), cytokine induced effect (additive) on the multiplication rate of recipient's target cells and the ON/OFF switch for cytokine (drug) application. This modeling scheme can also be operative after HSC transplantation, if needed and then application of cytokine will affect the transplanted cells of donor origin through m_DCYT.

To keep the chronic GVHD under control immunosuppressive drugs (d_im) are generally suggested. Present model has the provision in the assessment of such therapeutic management. For modeling it has been considered that immunosuppressive drugs affect different leukocytic subpopulations of both host and transplanted cells i.e., P2_g, P2_r, P2_s, B2_g, B2_r, B2_s, P2_e and B2_e. Depending upon the severity of GVHD /matching, clinical decision regarding the immunosuppressant drug (d_im) strategy i.e., drug dose, interval of consecutive drug application and duration becomes important. Such issues are addressed in the developed model. Therefore a_(i,j) of equation (8a-3) will be modified as follows.

d_im(k) is the amount of immuno-suppressive drug present within the system denote either the amount of applied drug (immunosuppressive) on the day of drug application (drug_imdose)or drug retention (drug_imret) on the subsequent days as a fraction of the previous day's applied drug (it is assumed that in each corresponding day of drug application a certain amount of drug is cleared from the system). Drug application is started on the day k_imsd, n_im is the number of cycles, t_imd is the interval of two subsequent drug applications in days and C_im is the ON/OFF switch for drug application is introduced through above scheme. And d_imsensh and d_imsensD are two positive fractional numbers representing the immunosuppressive drug sensitivities of the host cells and donor cells respectively (Table 2C). Further the controlling efficacy of immunosuppressive drug will depend on degree of GVHD (G), minimum drug effective set day (GVHD_set), drug dose and drug type by affecting Cal_GVHD with a fractional multiplier (GVHD_res). The operational relationship is represented as follows:

A rigorous simulation exercise with equations (1-8) has been carried out in MATLAB 6.5. Simulation studies may help in revealing the rationale of different therapeutic strategies for HM and HD.

With the initial parametric values as mentioned in Table 1, extensive simulation exercises have been carried out. Simulation studies show that under freely growing condition, there is an exponential growth for the malignant cells while a gradual decay of normal cells (nonmalignant) in the long run (Figure 9) (Table 5). With the increase in malignant cell population, normal cells of the other lineages i.e., erythroid and megakaryocytoid lineages are gradually decreased. Simulation study with the increase in multiplication rate of S_g shows that nonmalignant stem cell population falls earlier than the previous condition due to generation of more number of malignant cells due to the conversion. Disappearance of normal stem has been found earlier with increase in multiplication rate of drug resistive leukoblast cell. These earlier fall in normal stem cell population is due to development of extra toxicity burden to the nonmalignant cell population. Free growth pattern of this simulation resembles the acute leukemia case.

Simulation is carried with an assumption that myeloablative (MYL) chemotherapeutic regime is applied on the 30^th day after the diagnosis of leukemia (day 1). Drug regime is considered as the amount of total drug = number of cycles × drug dose. Drug dose, cycle and interval between two successive drug application and MYL drug sensitivity for every cell are presented in Table 2. Efficiency of drug also depends on drug retention, in simulation it is considered as 80% of the previous day. In each case drug clearance rate per day has been also considered as 20% of the previous day (Table 2). With MYL drug application there is a gradual decay of different cell populations. After stoppage of drug regimen the malignant cells again start growing. In the peripheral blood, the effect of regaining of cells compared to the precursor level is delayed, however, the effect is more prominent in the stem and progenitor cells present in the bone marrow. This indicates that MYL chemotherapy is unable to bring the leukemic patients in remission (Figure 10) (Table 5). It is to be noted here that though RBC and platelet count reduced due to MYL therapy and after maintaining a reduced level for a longer period of time, both reduces further due to effect of cancer cachexia that has been incorporated into systems equations. Simulation is also carried out with varying number of dose, cycle and interval between two successive drug applications. With the increase in dose, though there is more chance of removal of malignant stem cells however, that may bring down the RBC and platelet population very low and with this level patient may not survive in reality [35-40] i.e., when there is no exogenous RBC and/or platelet transfusion, T_ex=0.

High dose Myeloablative chemotherapy application: System was simulated with different combinations of drug doses (Drug_myldose), cycle numbers (n_myl), intervals (t_myld) with the same initial conditions as depicted in Table 1 and followings are noted.

Combination I: Drug_myldose=0.7, n_myl=8, t_myld=15, T_ex=0 (Figure 10A)
Combination II: Drug_myldose=0.8, n_myl=8, t_myld=15, T_ex=0
Combination III: Drug_myldose=0.7, n_myl=9, t_myld=15, T_ex=0
Combination IV: Drug_myldose=0.7, n_myl=8, t_myld=20, T_ex=0

Observation with Combination II: When the system was simulated with increase in MYL dose, keeping other conditions unchanged, it was observed that nonmalignant stem cell survived for longer duration (15 days more) than the Combination I. Again after completion with Combination II the nonmalignant stem cell count took more time (136 days more) to reach normal level compared to Combination I. The system was found to continue with lower RBC level than the Combination I.

Observation with Combination III: With Combination III it was found that after completion of MYL therapy, the nonmalignant stem cell count took more time (34 days more) to reach normal level compared to Combination I however quicker by 102 days compared to Combination II. Again simulation with Combination III (keeping other conditions unchanged) it was found that nonmalignant stem cell survived for shorter duration (3 days less) compared to the simulation with Combination II, but 12 days more than with Combination I. However with this combination the system continued to stay with lower RBC cell count compared to Combination I though higher than Combination II.

Observation with Combination IV: Simulation with the increased intervals (Combination IV) it was found that nonmalignant stem cell survived for only 10 days more compared to Combination I; however, 5 days and 2 days lesser than with Combination II and Combination III respectively. Completion of MYL therapy with Combination IV nonmalignant stem cell count took 52 days more than Combination III to reach normal cell count level, however earlier by 54 days than Combination II but longer by 82 days than Combination I. And the system continued to stay with lower RBC cell count compared to Combination I though higher than Combination II or Combination III.

Low dose chemotherapeutic drug application: The system was simulated with the same initial conditions as depicted in Table 1 but with following combinations of drug doses (Drugmyldose), cycle numbers (nmyl) and intervals (t_myld).
Combination V: Drug_myldose=0.35, n_myl=8, t_myld=15, T_ex=0 (Figure 10B)
Combination VI: Drug_myldose=0.40, n_myl=8, t_myld=15, T_ex=0
Combination VII: Drug_myldose=0.35, n_myl=9, t_myld=15, T_ex=0
Combination VIII: Drug_myldose=0.35, n_myl=8, t_myld=20, T_ex=0

Observation with Combination VI: When the system was simulated with slightly increased lower dose of chemotherapeutic drug (Combination VI) then it was found that after completion of drug application, the nonmalignant stem cell count took almost same time as with Combination V (only 2 days more in later case) to reach normal cell count level. Again simulation with Combination VI shows that nonmalignant stem cell survived for longer duration (11 days more) than Combination V. It was found that with this combination, the system continued to stay with lower RBC cell count compared to Combination V.

Observation with Combination VII: When the system was simulated with increased cycle number (Combination VII) it was found that nonmalignant stem cell survived for 9 days more compared to Combination V, same as with combination VI; but after completion of Chemotherapy with Combination VII the nonmalignant stem cell count took same time as Combination V to reach normal cell count level though quicker than combination VI (only 2 days earlier). When simulated with Combination VII the system continued to stay with lower RBC cell count though higher than with combination VI.

Observation with Combination VIII: After completion of chemotherapeutic drug application with Combination VIII, the nonmalignant stem cell count took 5 days lesser than both the Combination V and VII to reach normal cell count level; however, 7 days earlier than Combination VI. Again simulation with Combination VIII shows that nonmalignant stem cell survived only 4 days more compared to Combination V, 7 days lesser compared to Combination VI and 5 days lesser compared to Combination VII. It was observed that the system continued to stay with lower RBC cell count compared to Combination V though higher than Combination VII.

In all the above cases though P2_r and P2_s were reducing initially with the application of drug but with the stoppage of drug they started growing exponentially. Similar observations have been found for malignant cells in mature cell compartment.

RBC and platelet transfusions are occasionally needed during high dose chemotherapy, as though the number of both RBC and platelet is regained after stoppage of MYL/chemo therapeutic drug, but it may reach so minimal level that a patient would not survive (as shown in Figure 10b). Hence during the application of MYL/chemotherapy drug exogenous RBC and platelet transfusions are needed as supportive therapy (i.e., T_ex=1). The effect of such supportive treatment procedure is included in model and studied through simulation. It is assumed that before application of each MYL/chemotherapy dose, the need for supportive therapy is being checked. On the day of MYL drug application if RBCs and/or platelets levels were found below RBC_llm and Platelet_llm respectively, packed cells of 1,00,000 of RBC and 50,000 of platelets are assumed to be transfused; however, for lower dose of chemotherapeutic drug application, packed cells of 8,00,000 of RBC and 8,000 of platelets are assumed to be transfused (Figure 11). Transfused RBC and platelets have the apoptosis rate which are denoted by a_{T_RBC} (= 0.07747) and a_{T_Platelet} (=0.017) respectively.

System was simulated in presence of RBC and platelet transfusions with Drug_myldose=0.7, n_myl=8 and t_myld=15. It has been observed that though the recovery of the nonmalignant stem cell population into the normal cell count level was delayed by 109 days but they reduced to zero by the toxicity generated by the malignant cells almost at the same time (only 4 days delayed in presence of transfusion compared to Combination I).

System was simulated in presence and absence of RBC and platelet transfusion with Combination V and observed that the recovery of the nonmalignant stem cell population into the normal level was same as the Combination V, but in presence of transfusion with Combination V the nonmalignant stem cell population remained in normal level 5 days more. However, RBC level after the completion of Combination V with transfusion showed slightly lesser than the Combination V alone. Though transfusion helped in maintaining the normal RBC and platelet level in the system but the gap period for transfusion allows the malignancy to increase that, in turn, generates toxicity burden and thereby decrease the nonmalignant cell population.

For the treatment of leukemia through immunotherapy, particularly, cytokine based therapies are suggested after MYL therapy. It is assumed that applied cytokine affects only P2_g cell types (through differentiation process) and increases B2_g cell population which in turn, help in removing leukemic cells from the system. Simulation results are shown in Table 5. With present parametric values simulations suggest that cytokine drug application after MYL are unsuccessful in removing leukemic cells from the system (Figure 12); however, cytokine is applied after low dose chemotherapeutic drug application then malignant cells can successfully be removed. This observation may corroborate the findings of O' Brien et al 1995 [29], where authors showed that residual leukemia after chemotherapy was successfully removed with interferon treatment.

Myeloablative (high) or low dose chemotherapeutic drug applications followed by cytokine drug applications were observed when the system was simulated with following combinations of cytokine dug application after High dose MYL (Drug_myldose=0.7, n_myl=8, t_myld=15, T_ex=1) and low dose chemotherapy (Drug_myldose=0.35, n_myl=8, t_myld=15, T_ex=1). The initial parametric values including the assumed effects of cytokine drug application are presented in Table 2. It was observed that cytokine drug dose (Drug_CYTdose), drug application cycle (n_CYT), interval between two successive drug applications (t_CYTd), minimum number of WBC (Min_B2g) required to start immunity related killing (i.e., to turn kill_im =1 ) in presence of cytokine drug and malignancy killing efficiency (kill_eff) play important role in terms of removal of malignancy and maintenance of normal cell population. Rigorous simulation exercises have been carried out to find out the optimum values for cytokine drug dose /application cycles/intervals/minimum WBC level to start immunity related killing (Min_B2g) and malignancy killing efficiency (kill_eff) by cytokine mediated killing.

Application of low dose chemotherapy (Drug_myldose=0.35, n_myl=8, t_myld=15, T_ex=1) followed by cytokine drug application: With parametric variations the followings observations are noted.

Observation I: The system was simulated with different drug (Drug_CYTdose) level when n_CYT =15, t_CYTd=25, Min_B2g =224 and kill_eff =80. The cytokine drug was started from day 150 i.e., 5 days after the application of last dose of chemotherapy. It was observed that drug resistive stem cell disappeared from the system on day 568, 661 and 706 when Drug_CYTdose were 0.8, 0.76 and 0.755 respectively. i.e., with application of increased drug level the removal of malignancy becomes faster. And in all the cases non-malignant cell population in all compartments were found either maintaining normal cell count level or growing to reach normal cell count level.

Observation II: The system was simulated with different kill_eff when Drug_CYTdose=1.2, n_CYT =15, t_CYTd=25 and Min_B2g = 224. The cytokine drug was started from day 160 i.e., 15 days after the application of last dose of chemotherapy. It was observed that drug resistive stem cell disappeared from the system on day 320, 341 and 387 when kill_eff were 100, 90 and 80 respectively. i.e., with application of increased killing efficiency the removal of malignancy becomes faster. And in all the cases, non-malignant cell populations in all level were found either maintaining normal cell count level or growing to reach normal cell count level.

Observation III: The system was simulated with different minimum WBC level (Min_B2g) to start immunity related killing when Drug_CYTdose=1, n_CYT =15, t_CYTd=25 and kill_eff =80. The cytokine drug was started from day 160 i.e., 15 days after the application of last chemo dose. It was observed that drug resistive stem cell disappeared from the system on day 393, 393 and 412 when (Min_B2g) was 223, 222 and 224 respectively; i.e., with the lowering the minimum effective WBC level ((Min_B2g) for starting the immunity related killing, malignancy removal becomes faster. And in every the cases and non-malignant cell population in all level were found either maintaining normal cell count level or growing to reach normal cell count level. But when simulated with (Min_B2g) =225 with other conditions unchanged malignancy found to increase and all the non-malignant cell populations reduced to zero (stem cell disappeared from the system reduced to zero on day 711). In this case the minimum Min_B2g has been found as 224. In this present parametric setting above this value the system will fail to control malignancy. In case of cytokine application (Drug_CYTdose = 1) after chemotherapeutic strategy with transfusion though there is a persistence of B1_r and B1_s cells having a number >2000 on day 1500 however with a decaying tendency.

Application of high dose MYL (Drug_myldose=0.7, n_myl=8, t_myld=15, T_ex=1) followed by cytokine drug application: With parametric variations the followings observations are noted.

The system was simulated with different drug (Drug_CYTdose) levels when n_CYT =15, t_CYTd=25, Min_B2g =224 and kill_eff =80. The cytokine drug was started from day 170 i.e., 15 days after the application of last MYL dose. It was observed that keeping all initial conditions as depicted in Table 1 unchanged with the increase in Drug_CYTdose level, nonmalignant cells remained in the system for longer time and recovery of nonmalignant stem cell to normal cell count level were faster (nonmalignant stem cell reached to normal cell count after completion of MYL on day 338^th, 304^th and 291^st when Drug_CYTdose level were assumed as 1, 4 and 6 respectively). For Drug_CYTdose =6 the malignancy were found to disappear from all compartment (resistive stem cell disappeared from the system on day 363) and normal cell populations were found either in normal cell count level or growing to reach normal cell count level. Though malignancy were found to be disappeared from all compartmental level with Drug_CYTdose = 6 but it may be unrealistic (as Drug_CYTdose×m_CYT =6×0.07=0.42 i.e., 42% extra growth in Leukoblast per day).

The failure to control malignancy with MYL followed by cytokine is due to the fact that high dose MYL bring down normal leukoblastic cells into such a lower level that may not be possible to boost to such high level so quickly to control malignancy. In present system model it has been assumed that with high dose MYL application when number of P2_g<30, the S_g>dt and P2_gdt are increased to 30 days and 20 days respectively than their normal value; whereas S_gdr and P2_gdr are decreased to 0.0165×{3-(m_Sg + a_Sg)} cells/day and 0.765×{1-(m_P2g + a_P2g)} cells/day respectively) with respect to their normal value and a_Sg gets increased (to 0.48 cells/day) with respect to its normal value.

For simulation it is assumed that HSC transplantation has been done 25 days (on 180^th day after the day of diagnosis) after the completion of last dose of myeloablative (MYL) drug regime (cycles 8, drug dose 0.7 and interval 15 days). The day for HSC transplantation is chosen when drug level within the system reaches below a desired level (say < 1% of the applied drug). It is to be noted here that in simulation MYL drug was applied on day 30^th, 45^th, 60^th, 75^th, 95^th, 115^th, 135^th, 155^th. The unequal gap between successive drug application were due to transfusion of RBC and Platelet on 90^th, 110^th, 130^th and 150^th days for low RBC and Platelet count. The assumed characteristics parameters of the transplanted HSC are presented in Table 3. As there is, no GVHD, so the maturated donor lymphocytes that are produced from the transplanted HSC will kill only the residual malignant cells with a factor kill _b2e→hm = 7.142 (Table 4). With the present simulation condition, it has been observed that around ~324 days after HSC transplantation all the malignant cells were totally disappeared from the recipient's system, no relapse was noticed even after 2000 days (Figure 13) (Table 5). Chimera (Mixture of donor and host cells) in different cellular level is observed. Simulation has been carried out with different killing efficiency of lymphocytes originated form transplanted cells.

Observation 1: Similar type of observations (as in Figure 13) were noted for kill_b2e→hm = 7.14, and 7.13 with other conditions unchanged.

Observation 2: With the other conditions unchanged simulation exercises with varying killing efficiencies show that kill_b2e→hm = 7.13 is the minimum required value for the complete removal of malignancy from the system. It signifies that the exogenous cell should have minimum killing efficiency for complete malignancy removal. Simulation with the value below 7.13 there is a decreasing trend of normal cell level due to toxicity generated by the increasing level of malignant cell. However, when the system was simulated with increased number of transplanted HSC (S_e = 12) having kill_b2e→hm = 7.125, the malignancy is removed from the system and normal cell population is sustained. It signifies that donor stem cells with lower killing efficiency can be compensated with increasing its number.

Observation 3: With other conditions unchanged when simulation is carried out with S_e = 12 and kill_b2e→hm = 7.125; however, transplantation is assumed to be on 190th day after the day of diagnosis, malignancy exist, in all level.

Observation 4: With other conditions unchanged donor stem cells with kill_b2e→hm = 7.125 and transplantation is assumed to be on 190^th day after the day of diagnosis with more number of HSC (S_e = 20) there is a complete removal of malignancy from the system. It signifies that the delay (in days) in transplantation may compensated (for removing leukemic cells) by increased donor cell number by multiple source of cord blood engraftment.

Therefore delay in transplantation can be compensated either with increased donor cell number or donor cells having property of generating lymphocytes of higher killing efficiency. Moreover, in the model, the application of MYL drug strategy (drug dose and cycle number and interval) can be changed/adjusted. This allows the investigator to check the condition of attaining a minimal tumor load suitable for transplantation.

HLA mismatch is manifested as graft versus host diseases (GVHD). Simulation is carried out considering chronic GVHD with two grades: high and low (Figure 14). The condition of HLA mismatch can also be studied through this model. The effects of HLA mismatch are incorporated in the model using several killing factors: exogenous lymphocytes developed from transplanted HSC reacted against the recipient (host) malignant cells (kill_b2e→hm) (same as for 100% HLA matching case) and host normal cells (kill_b2e→g). The effect of GVHD is manifested through development of GVHD against normal tissue host like liver, gut and lungs. And the intensity of such GVHD related tissue damage can be measured using a marker (MR), where MR= f(B2_e, G).

Simulation has been carried out in different conditions and results are presented as below.

Observation 1: Simulation study shows that in case of transplantation with HLA mismatch, removal of malignancy is found to be much faster than 100% HLA match case (Table 5). For example Sr reduced to zero on day 504^th when transplanted with 100% HLA matched donor cell; however, when transplanted with low and high mismatched donor, removal of S_r is on day 224^th and 217^th respectively.

Observation 2: As kill_b2e→hm increases while other factors unchanged normal cell population also reduces faster (Table 5) along with the faster removal of malignancy.

Observation 3: Keeping other conditions unchanged with delay in transplantation though normal cell survived for longer period but malignancy removal was found delayed relative to the previous cases.

Observation 4: Transplantation with high mis-match shows faster removal of malignancy faster than transplantation with lower mis-matched donor and it also reduces normal cell population earlier even more than free growth condition.

With the present parametric settings (Table 4) it is observed that with the increase in HLA mis-matching, GVHD appears much earlier (Figure 14). The effect is much more prominent on the mature cell level than the precursor cell level.

When allogenic hematopoietic stem cells (allo-HSC) with a certain degree of HLA mis-match are transplanted, patients need to be treated with an immunosuppressive drug for the management of GVHD. It is worthwhile to mention here that depending upon the severity of mismatch, determination of immunosuppressive (IMS) drug dose and drug schedule could be the important aspects for clinical cases. The developed model has the provision to test the efficacy of the IMS drug for the management of GVHD. This model has the flexibility of choosing IMS drug scheduling strategy.

In simulation different IMS drug strategies have been applied to different HLA mis-match cases as mentioned previously. In all the cases, allo-HSC transplantation has been performed on day 180 after completion of conventional MYL drug regime. Simulations have been carried out considering two different types of IMS drugs (i.e., IMS-1 and IMS-2). Both types of drugs have the same mechanism of action; however, IMS-2, signifies higher dose (double) than IMS-1 and once one drug is applied it continued throughout the simulation period. Strategy started 30 days after HSCT and IMS-2 drug strategy started 60 days after HSCT (Figure 15).

IMS drug dose, drug application interval, drug retention and the considered sensitivities of these drugs to different cells are depicted in Table 2C and results of such simulations are presented in Table 5. Simulation study indicates that IMS drug application in controlled level can keep GVHD under control to remove the malignancy from the system and thereby increase the sustaining period of normal cell population.

However, with the application of excessive IMS drug dose (keeping other conditions unchanged) malignancy may relapse in system.

So far our simulation exercises have shown the management strategies for chemo-unresponsive leukemia (worst) case. However, several leukemia cases respond to the conventional chemotherapy regimes. Simulation study with the increase in drug sensitivity ~10 times for the malignant cells only as mentioned in Table 2 and with the application of 8 cycles of high dose MYL in presence of transfusion can bring the considered worst case (as in Table 1) into remission; however, it does not bring to remission when 6 cycles of the same therapeutic regime is applied.

Simulation study with the developed model also carried out to show the features of chemo-responsive leukemia (remission) cases with change in some initial parametric values (Case II). Free growth condition of such case is initialized by setting m_Ss = 0.07, m_Sr = 0.075, a_Ss = 0.032, a_Sr =0.032, m_P2s = 0.0345, m_P2r =0.0345, a_P2s = 0.00359, a_P2r =0.00359, 1.3 times, 1.3 times and 1.5 times of S_rdr, P2_sdr, and P2_rdr of worst case respectively while rest of the values are kept unchanged as mentioned in Table 1 (Figure 16). Remission is possible with the application of 6 cycles of high dose of MYL with 5.5 times MYL drug sensitivity towards malignant cells as mentioned in Table 2. Simulation study with this case depicts that the number of malignant stem cells have slower growth rate than the worst case and the leukemic/malignant cells have greater drug sensitivity

Keeping all the parametric values and the MYL drug regime same as of the Case II remission case while kept back the parametric value of m_Sr and a_Sr as of the worst case fails to bring the system into long-term remission and simulation study depicts the presence of residual disease (Figure 17).

Previously simulation studies resemble acute leukemia cases. However, change in initial parametric values in the developed model can show the features of chronic leukemia cases. Chronic myeloid leukemia (CML) is considered as the differentiation defect with an arrest of cells at different stages of morphological maturation [40]. Free growth condition is initialized by changing S_rdr, P2_sdr and P2_rdr to 1.35 times, 1.55 times and 1.95 times respectively of worst case while rest of the values are kept unchanged as mentioned in Table 1 (Figure 18). In this case drug application is considered daily with drug dose 0.037 (6/1000 th of total high dose MYL drug) for 1500 days. It is to be noted here that in CML (chronic myeloid leukemia) cases are generally treated with the targeted drug imatinib and this drug is being applied daily [41-43]. Simulation is carried out with two different drug sensitivities for malignant leukemic cells – one simulation is done with the same sensitivity as mentioned in Table 2 (Case I) and another simulation is carried out with the 5.5 times of drug sensitivity for malignant cells (Case II). However in both the simulations it is also considered that such drug strategy has minimal effect on the normal cells of the hematopoietic system. Simulations studies depict that with the increase in drug sensitivity of malignant cells can able to bring the leukemic free survival in the long-run.

In reality medical treatments are dependent on two approaches – population based randomized clinical trials and qualitative judgments of the subjective experiences of clinicians. Hence treatment dealing procedure of an individual patient becomes empirical. However, physicalist approaches demand for mechanistic explanation. Mathematical models with dynamical nature may provide the way-out to an understanding of a disease process in a much more quantitative manner. It is also true that till now the qualitative judgments and subjective experiences cannot be denied while dealing with the individual clinical cases [44,45]. This fact is mostly denied as mere information and/or unable to comprehend its importance with the cultural training of physicalist school of thought [46]. Ironically majority of developed mathematical models are done with abstract assumptions of variables that are unable to hold any applicative significance to clinical practice for dealing of individual patients though the model behavior synchronizes some real data. Biological discipline has a long history of training program to comprehend the importance of biological variability and analyzing the natural phenomena. Though physiology and medicine both deal life processes, but due to technical and instrumental dependency, both follow reductionism; contrarily, biological analysis is viewed as and confined within the subjective experiences and as a result there is absence of time varying dynamical data (of individual cases) in medical literature. Systems Biology and its translation towards Systems Medicine seems to be essential to meet this gap. Towards this, presently, development of analytical model with rationality can be done [47,48]. Simulation studies present here may reflect the subjective experiences of clinicians and the developed model may be helpful in providing the mathematical rationality towards their domain knowledge.

Previously DODE is being applied by several workers for mathematical modeling of both chronic and acute leukemia. Hence it is imperative to state that the same mathematical strategy can be applicable for both chronic and acute leukemia with varying in parametric values for different variables. Such mathematical strategies are also used for the understanding of the efficacy of different therapeutic strategies like myeloablative chemotherapy, tyrosine kinase inhibitor and/or T-cell based therapy, but different workers have used different variables in the development of the model with a focus on different therapeutic strategies. Hence each model is being apart from other. However, in clinic there is a need to evaluate and/ or predict the response of a particular individual patient to a particular therapeutic strategy among the different available therapeutic options. Hence there is an immense need for the development of an analytical model that can make a comparative assessment among the different available therapeutic options. To address this issue we have developed our model using DODE with an incorporation of feedback at different systems level. Our model can be aptly be fitted with both acute and chronic leukemia cases by setting the initial parametric values of the variables. Similarly setting of different drug parameters can be used for the assessment of the efficacy of different drugs that are being used currently in conventional clinical practice. The model has the flexibility in chemotherapeutic drug application (dose, duration and cycles), so the model can also be suited for the assessment and prediction regarding the dynamical features of different patients with different patho-physiological states.

The present model includes the complex nonlinear behavior of the different time varying system elements of the hematopoietic system and their inter-relationships. The complex inter-relationships have been reflected when observed with minor parametric variations. In freely growing condition it has been found that minor increase in multiplication rate of any normal cell type failed to increase its population rather it decreases. It is due to the consideration of mutability of normal cells to the malignant phenotypes. That in turn generates toxicity cachexia to bring down the normal cell population. This may be major cause of infection and bleeding in leukemia.

With the inclusion of such complex nature in the system model, the management of MYL drug application becomes more difficult for the maintenance of normal cell count of the normal cell population in one hand and removal of malignancy on other hand. With increase in MYL dose amount the malignant cell count can be found to decrease. However, simulation shows that excessive drug application for removal of malignancy may collapse the normal cell population due to the drug related toxicity which in turn, enhances the proneness to infection and bleeding. This reflects real clinical scenario. Contrarily, residual disease may exist due to application of lower dose of chemotherapy.

So application of MYL drug dose should be controlled judiciously with a tradeoff between maintenance of normal cell count for longer period of time versus elimination of malignant cell populations from system. This can only possible if there is an increase in drug clearance rate and/or drug sensitivity for malignant calls.

To circumvent the MYL drug related toxicity, RBC and/or platelet transfusion is the regular practice in clinic. This has been included in our system model. Simulation study shows that application of transfusion strategy may make a delay in the consecutive MYL drug application; which in turn, may allow the malignant cells to grow. This again becomes detrimental for the patient due to toxicity generation by the malignant cells.

Our simulation study shows failure of MYL therapy with or without transfusion. This signifies that we have considered a worst clinical case. Simulation has been carried out further to reveal the expected management procedure for drug resistive worst leukemia case. Considering this worst case, further simulation shows that the application of MYL drug followed by cytokine is also unable to bring the desired success under the same leukemic condition; however the application of lower doses of chemotherapeutic drug (preparatory regime) followed by cytokine strategy may bring the desired therapeutic success. This is due to prior application of MYL before cytokine application is unable to potentiate the immune-competent lymphocyte due to removal of cells of the existing leukocytic lineage. This observation corroborates the earlier observation which indicates that the leukemia free survival is dependent on the existing functional T-cell [7,12].

The considered cell count can be routinely assessed by hematological investigation, different immunological parameters like killing efficiency of lymphocytes can be assessed with different immunological investigations, and the considered chronic GVHD can be assessed by different qualitative clinical investigations.

The present model has the provision to test the efficacy of transplantation under different conditions in terms of HLA matched condition, transplantation day, delay in transplantation, malignancy killing efficiency by the lymphocyte of donor origin and number of transplanted cells. The considered cell number for transplantation with respect to volume can be equated in terms of body weight. With 100% HLA matching, simulation with donor cell having lower killing efficiency no improvement in terms of malignancy removal was noticed though there was no GVHD effect but GVL is failed. However increase in transplanted donor stem cell number (in reality this can be done with multiple donor engraftment) keeping other conditions unchanged, are able to remove the malignancy by developing GVL to threshold level. Again when simulated with delay in transplantation day there is a failure in controlling leukemia. Since during delay time leukemic cells reached to such level that develop toxicity which in turn, kills donor cell population. So determination of minimal tumor load on the day of transplantation is another major consideration. With this model, determination of such criteria can also be assessed. Hence, simulations suggest that this failure may be compensated by increasing initial transplanted cell number or donor cell having higher killing efficiency. This can be achieved by increasing the lymphocyte population of donor origin by applying cytokines. Model has such provision to incorporate such effect by m_DCYT.

With increase in HLA mis-match keeping other factors unchanged removal of malignancy also becomes faster. Though transplantation with high GVHD are able to remove malignancy faster than transplantation with lower mis-matched donor but it also reduces normal cell population faster even more than free growth condition. To tackle this, immunosuppressive drug application is suggested. Model has the provision in the assessment of such therapeutic strategy. Even the model flexibility is also made in a manner to test the second choice of immunosuppressive drug application, if first one fails. Our simulation again suggests that judicious choice of immunosuppressive drug application may eliminate both malignancy and GVHD from the system, as malignancy may predominate in the system due to application of excessive amount of immunosuppressive drug. Excessive application of immunosuppressive drug may enhance the susceptibility of the patient to infection and bleeding. This is due to removal of functional donor lymphocytes by the excessive amount of immunosuppressive drug. However such efficacy of immunosuppressive drug will also depend further on the degree of HLA mis-match.

It is worthwhile to mention here that with the availability of individual patient's dynamical data, fine adjustment of simulation studies with this model could be possible. It is to be mentioned here that our developed model is flexible in nature. In our model all the parametric values of every variable can be changed and can be adjusted accordingly to the individual clinical cases as well as choice of the investigator. The major objectivity of the present work is to provide a mathematical/analytical platform to make a comparative assessment among the different available therapeutic options [47,48]; so that, a clinician can make a prior assessment regarding the efficacy of a particular therapeutic strategy before applying it to an individual patient. This model can also be aptly fitted to test the efficacy of different therapeutic strategies on different individual cases with different types of leukemia. Moreover, model can be applied for prognostic assessment in individual cases. This quantitative assessment tool will help to transcendent the empirically based clinical practice in future.

Recent time several therapeutic modalities for the treatment of leukemia are being suggested by the scientific community; however, due to individual variation in response to a particular therapy it is difficult to rationalize which therapeutic modality would be the choice of the treatment for a particular individual patient. The developed model can help in the apriori assessment regarding the efficacies of different therapeutic modalities. Simulation studies with this developed model indicate that under the same leukemic condition, when conventional chemotherapeutic regime is failed to provide long-term leukemia free survival immune-boosting mechanism by any means like cytokine application or allo-HSCT could be the determining factor for long-term leukemia free survival. Our simulation study also provides the insight that the time of immune stimulation is another determining factor for the success of leukemia therapy. The developed model is flexible in nature. Clinicians can change any parametric value according to the need of the individual clinical cases and make an in-silico assessment regarding the suitability of a particular therapeutic regime depending on the individual leukemia cases before applying it in reality.

Authors acknowledge the eminent founder members of Society for Systems Biology & Translational Research for their valuable comments to enrich this work.