Mathematical Modeling of Leukemia Chemotherapy in Bone Marrow

Acute Lymphoblastic Leukemia (ALL) accounts for the 80% of leukemias when coming down to pediatric ages. Survival of these patients has increased by a considerable amount in recent years. However, around 15-20% of treatments are unsuccessful. For this reason, it is definitely required to come up with new strategies to study and select which patients are at higher risk of relapse. Thus the importance to monitor the amount of leukemic cells to predict relapses in the first treatment phase. In this work we develop a mathematical model describing the behavior of ALL, examining the evolution of a leukemic clone when treatment is applied. In the study of this model it can be observed how the risk of relapse is connected with the response in the first treatment phase. This model is able to simulate cell dynamics without treatment, representing a virtual patient bone marrow behavior. Furthermore, several parameters are related to treatment dynamics, therefore proposing a basis for future works regarding childhood ALL survival improvement.


Introduction
Leukemia is a malignant disease originating in the bone marrow.Particularly, it arises from a disruption in hematopoiesis, the process in charge of blood cells production [1,2].Hematopoiesis is usually depicted as hierarchical tree, in which a hematopoietic stem cell (HSC) can differentiate into other cells from the lymphoid or myeloid line.Depending on the branch, platelet, red blood cells and lymphocytes are produced.Leukemia is not only distinguished by the linage but also depending on the maturation stage of the cancer cells.Acute Lymphoblastic Leukemia (ALL) is caused by cells with fast growth in the lymphoid branch, with special incidence in pediatric patients [3,4].
Each type and subtype of leukemia has an associated protocol that specifies the therapeutic recommendations, and that varies from country to country [5].In particular, SEHOP-PETHEMA-2013 protocol is used in Spain to treat pediatric patients with ALL.This protocol was designed by Spanish Society of Pediatric Hematology and Oncology along with Spanish Hematology Treatment Program and it is regularly reviewed [6,7].
There are different treatments and medications depending on the risk group assigned to a patient: standard, intermediate or high.In such protocols, patients go through successive phases: induction, consolidation, re-induction and maintenance.Each phase consists of the combination of chemotherapeutic agents and corticoids with varying schedules.As reflected in other works [8,9], the progressive improvement and modification of these protocols has increased survival rates to around 80 to 85%.Nonetheless, 15 to 20% of patients fail to achieve long term remission and therefore relapse [10].Apart from adding to and modifying these protocols, recent reviews have highlighted the need for new therapies and approaches that would tackle the disease differently [11,12].This also includes the improvement of risk assignment, the management of secondary effects and the study of personalized medicine, whose power lies not only in treatment but also in prevention [13].
In this sense, mathematical models have recently found a wide variety of applications in medicine [14,15].In particular, the discipline of mathematical oncology has been able to approach recent oncological issues [16,17].
To understand cell dynamics, several works have described biological processes related to cellular behaviour [18,19] and also related to hematopoiesis and lymphopoiesis [20][21][22].This has allowed the mathematical characterization of the hierarchical structure of blood cell lineages [23][24][25][26] and their development in bone marrow.Specifically, ALL has been considered in previous mathematical models [27][28][29] with an structural vision of leukemic cells similar to that of healthy cells.In order to fulfill their potential, these models should be validated against appropriate treatment data [30][31][32].Based on the current protocols to treat ALL patients [33,34], some mathematical models describe the drugs' effect regardless of the leukemia type [35,36], while other focus explicitly on the myeloid linage [37][38][39].Other mathematical models describe alternative strategies to chemotherapy, such as CAR T-cell therapy [40][41][42][43], which is successful in patients that do not respond to the standard treatment.
Inspired by the examples above, we set out in this study to build a mathematical model that describes the evolution of ALL and its response to treatment in the bone marrow.This type of leukemia has received less attention from the mathematical community, and can benefit just as much.The mathematical model would allow the derivation of theoretical results (like stability analysis) but also clinical information like the relative importance of parameters, the estimation of response times or the influence of dosage and schedule.It could eventually be used to fit individual patient data and obtain insights into the possibility of personalized therapy for this disease [5,28].
In this paper, we focus on the simulation of standard risk ALL pediatric patients in the first phase of the treatment, or induction phase, which lasts 28 days.In terms of drugs, this phase comprises the corticosteroid Prednisone, the chemotherapeutic agents Vincristine and Daunorubicine, and Asparaginase, each with their own dosages and infusion times.In addition, triple intrathecal is given to prevent central nervous system disease involvement.The modelling of this treatment stage requires an analysis of the bone marrow environment in addition to the assumptions related to the growth of leukemic cells.We study this scenario and simulate how treatment affects both leukemia and the lymphoid linage in bone marrow.Besides, we propose a classification of the patients depending on their leukemic cells development to predict which patients will respond to treatment correctly.
The work has been organised as follows: In Sec. 1 we describe an existing healthy lymphopoiesis mathematical model and extend it to include leukemia.This extended model will also include treatment administration.In Sec. 2 we show model simulations of leukemia growth and its evolution under therapy, obtaining realistic ranges for parameters and highlighting the most relevant ones in terms of response.In Sec. 3 we discuss the results obtained.

Material and methods
In this section we describe a mathematical model of healthy, homeostatic lymphopoiesis in the bone marrow.We then include the appearance of a leukemic cell and study its evolution.We finally approximate a treatment function following the current treatment scheme during induction phase and model its effect in the previous models.An overview of all such models is shown in Figure 1.

Healthy bone marrow model
Following the immunophenotypical characterization of B cells [21], previous mathematical models, have considered B lymphopoiesis as a sequential process dividing their maturation stages in Pro-B (early stage characterized by immunophenotypic markers CD10 + /CD45 -), Pre-B (intermediate stage with CD10 + /CD45 + ) and Transition cells (last stage with CD10 -/CD45 + ) [25,27].Mathematically, this can be regarded as three compartments respectively depending on each stage: C 1 (t), C 2 (t) and C 3 (t), along with their associated parameters related to proliferation and differentiation.This was translated in [26] into the compartmental where ρ i , for i = 1, 2, is the proliferation rate and α i , i = 1, 2, 3, is the transition rate related to each compartment i.
In this model the authors also took into account a signal in charge of regulating cell dynamics, s = s(t), with the following expression: where k is the inhibitory parameter and N the cells subpopulation in charge of the feedback signalling.According to the results in [26], we consider here N = 3 i=1 C i .This signalling was shown to affect predominately the proliferation rate due to biological reasons such as bounded growth or mathematical ones such as positivity of the solutions.Having reviewed this model, we expand it by including the appearance of a leukemic clone.

Modeling leukemia cells
Now, in addition to the three compartments C 1 , C 2 , C 3 , we consider a leukemic cells compartment, L. We assume that leukemic cells grow as a logistic function where the curve's maximum value is L max and ρ L is the logistic growth rate.Just like cells in a healthy hematopoiesis, a leukemic clone is characterized by a maturation stage [23,27].Therefore, we assume L to depend on the stage where the leukemic clone originates, i.e., from the Pro-B stage (ρ L = ρ 1 ) or from Pre-B stage (ρ L = ρ 2 ).Finally, γ L would denote the blood exit rate of leukemic cells.Here we do not consider clones coming from the transition stage due to the fact that ALL is produced by immature, proliferating cells.
We consider then the following equation to define the flux in the leukemic cells compartment: Finally, signalling s L in this compartment has a similar behavior to the healthy bone marrow compartments as in [26], and we assume it to affect the proliferation rate as in Model (1).Therefore, assuming leukemic cells should not be affected by their own signalling due to the evasion of growth suppressors [44].Nevertheless, healthy stages are influenced by cells in compartment L, due to the fact that the more cells in L, the less healthy cells develop due to the invasion of the bone marrow.Consequently, healthy signalling s in Model (1) should be modified to: In addition to these assumptions, we can also include a constant influx of healthy stem cells c 0 in our initial compartment C 1 [40], which results in the following model, including both leukemic and healthy B cells: with s as in Eq. ( 5) and s L as Eq. ( 4), where ρ L depends on the origin of the leukemic clone [24,27].
This model allows us to study different treatment regimes, to compare with actual results and to include and test what is known about the behaviour and effect of these drugs [5,34,45].

Modeling treatment for leukemia
As explained previously, there are different treatment protocols that can be applied [28], depending on the drugs used, their doses, or even the time of infusion.We consider the Induction I'A phase, as described by the SEHOP-PETHEMA-2013 protocol, for a standard risk patient.This phase consists of four drugs and lasts 37 days.Triple intrathecal is not included in the study since it does not affect the lymphoid branch.Each drug has a schedule with corresponding doses and days of administration, depending on the patient's body surface, as indicated in Figure 2. To simplify, we consider treatment days as a scale which corresponds leukemia detection day with day +0 and treatment beginning with day +1.The most important days for monitoring the patient are treatment days +8, +15 and +33.In those days, bone marrow sample extractions are carried out and lymphocyte levels are studied.In day +8 a blood extraction is done to assess response to prednisone.According to SEHOP-PETHEMA-2013 protocol, there must be less than 10 6 leukemic cells per milliliters of blood in day +8.The Minimal Residual Disease (MRD) is studied in bone marrow, whose positivity is defined by the presence of 0.01% or more leukemic cells in bone marrow [46].Therefore, patient responds to the treatment if leukemic cells involve at most 0.01% in bone marrow in day +15, and there are no blasts in day +33.
Following studies on the behaviour of each drug [47][48][49][50][51], we consider the treatment µ j (t) = µ j as a function that describes the amount of drug j in the bone marrow at time t, being j ∈ J = {P, V, D, A}: Prednisone (µ P ), Vincristine (µ V ), Daunorubicin (µ D ) and Asparaginase (µ A ). Once the dose is administered, the drug has an exponential decrease associated with that medicine half-life.
We then define µ j : R + → R + as with λ j related to each drug j half-life τ j measured in days: Each drug j has a different dose q j in several days D j ⊂ N, as shown in Figure 2. We define the dose administered as: We finally consider the total treatment function µ, represented in Figure 3, as the weighted sum of all drugs: where δ j indicates the influence of drug j on the total effect of the treatment.2) and with parameters from Table 1.
Finally, we include a term in Model ( 6) that indicates cellular death by chemotherapy.We consider this does not affect transition cells, owing to the fact that they do not have any associated proliferation.The resulting model is shown below (11).

Parameters estimation
We consider literature data for healthy bone marrow properties [26] along with parameters related to leukemia [40].Drugs half lives are taken from pharmacokinetics and pharmacodynamics studies.All parameter values are included in Table 1.The estimations for dosage assume a 25 − 30 kg body mass for a child, which in turn implies a body surface area equivalent to 1 m 2 [52].
The study of all estimated parameters is presented in Supplementary Information document.Control parameters related to drug, δ j with j ∈ J, have been selected by searching for values that implies that the patient responds to treatment correctly, i.e. considering an amount of blasts less than 10 6 blasts/ml in blood in day +8 and MRD< 0.01% in bone marrow in day +15 [53].In preliminary simulations we find that if δ j q j < 10 0 , the dose administered is ineffective.We further consider that δ j q j > 10 1 involves the death of the patient due to drug toxicity.We then set all δ j q j to the same order of magnitude 10 0 , 10 1 , ∀j ∈ J.In this way we ensure Model (11) can represent a virtual patient whose behavior is similar to real patients according to the amount of cells in each bone marrow aspiration.

Leukemic and healthy cells dynamics without treatment are properly represented by the model.
We now consider no treatment is applied, Model (6), i.e., Model (11) with µ j = 0, ∀j ∈ J.This model combines the healthy dynamics of immature B-cells and how leukemic cells stability values as initial inputs for the model.We simulate the appearance of a leukemic cell (at t = 0), depending on the leukemic clone originates, either in a Pro-B (ρ L = ρ 1 ) or a Pre-B stage (ρ L = ρ 2 ).Results of solving in a range to [0, 300] days are shown in Figure 4, displaying the range of days with more information in each case.
We next show the simulations for each case.Figure 4 represents Model (6) cells evolution when ρ L = ρ 1 is considered for the first case, and secondly, model for ρ L = ρ 2 is simulated.In addition, for each time t measured in days, the bone marrow cells proportions are shown for the different maturation compartments.
As it can be observed in Figure 4, both graphics represent the total bone marrow invasion by leukemia.Depending on the stage from which the first leukemic cell comes, leukemic cells proliferate more or less quickly.The appearance of a leukemic clone is set to be detected when the leukemic population accounts for 80% of blasts from the total population of B lymphocytes, [54].When the first leukemic cell is originated in Pro-B stage, ρ L = ρ 1 , the disease is detected around day 150.However, if the original leukemic cell comes from the Pre-B stage, leukemia is appreciated from day 223.Both two cases imply a progressive decreasing in healthy cells due to the fact that the invasive ability of leukemic cells prevents the normal development of cells in the bone marrow.

2.2.
The first blood extraction in day +8 of treatment allows to approximate the influence of prednisone (δ P ).
Treatment begins according to treatment schedule in SEHOP-PETHEMA-2013 protocol, after 80% of blasts are detected, producing a drop in both leukemic and healthy cells.
In order to find the optimal parameters to represent Model (11), we pay attention to the first blood extraction (day +8).To do so, we estimate relations between blasts in bone marrow and blood, owing to the fact that there is a correlation between them [55,56].1. Parameter values.L max is measured in number of cells and k in cell −1 while the rest of the presented parameters related to the proliferation and transition rates are measured in day −1 .On the other hand, medicine dosage q j is measured in mg/day and the influence of each drug is considered day/mg.Asparaginase is a particular case since is measured in U/day and its influence in day/U.Values of λ j are measured in day −1 due to the fact that they come from half-life values.Values used to initialize each case C i (0), i = 1, 2, 3, have been obtained from the steady values in healthy bone marrow models [26] along with our parameters.A single leukemic cell is supposed at the beginning, from which Model (6) starts, and the influx of cells c 0 is also obtained from Model (1) steady states.

Parameter
SEHOP-PETHEMA-2013 protocol defines a good response to prednisone if there is fewer than 10 6 blasts per milliliter of blood.We assume a pediatric patient with 25 − 30kg body mass, and due to the blood -body mass correlation [57,58], this assumption entails approximately 2.1 − 2.5 liters of blood.We consider a patient with 2.3l of blood to our study, and consequently, patients should have less than 2.3 × 10 9 blasts in blood.It implies that leukemic cells limit in the bone marrow should be less than 2.3 × 10 10 blasts at day +8 to obtain a good response to prednisone.We analyze all δ P possible values within our range in Table 1, taking into account leukemic cells evolution for the first eight days of treatment (see Figure 5).
As it is shown in Figure 5 (a), the amount of blasts in bone marrow in day +8 of treatment implies a direct consequence related to possible δ P values.Values lower than δ P0 imply that patients do not response to prednisone, while values higher than δ P0 entail a good prognosis for the patient due to the fact that in day +8 of treatment reach less than the limit of blasts assumed.
On the other hand, in Figure 5 (b) the development of leukemic cells according to Model ( 11) is depicted for five different values of δ P .Values which allow to reach blasts levels lower than 2.3 × 10 10 total blasts in day +8 corresponds to the green area in Figure 5 (a).

Prednisone and Vincristine influence values analysis implies a classification of patients
Let us how consider all parameters δ j in order to analyze their behaviour in Model (11), taking into account the number of blasts in bone marrow at day +15.
Based on the parameter sensitivity analysis (Supplementary Information), main influential parameter is δ P .To study their influence, we vary parameters δ P and δ V , related, respectively, to prednisone and vincristine, and check how many leukemic cells are there day +15 of treatment schedule in comparison to the healthy population (MRD).We fix parameters δ D = 2.5/30day/mg and δ A = 2.5/10 4 day/U given their minimal influence in the  11) is solved for 50 values of prednisone influence.For each δ P value, amount of blasts in day +8 is obtained and associated with a patient who responses (green points) or not (gray points) to prednisone.(b) Leukemic cells evolution depending on several δ P values.Green-dotted lines denote patients who response to prednisone, while gray-dashed lines represent bone marrow dynamics which do not reach any response to prednisone from data in day +8 (yellow vertical line).Horizontal solid red line represents the limit of blasts assumed by SEHOP-PETHEMA-2013 protocol to consider a good response of less than 2.3 × 10 10 blasts in bone marrow.sensitivity analysis.We focus on the variation of δ P values in comparison with δ V values to obtain the minimal value which implies the patient responds to treatment.We assume patient responds to treatment correctly if day +15 of treatment leukemic cells represent less than 0.01% of the whole lymphocyte population.
Values ranges considered for δ j , j ∈ J are those in Table 1.In Supplementary Information, we show a heat map which takes parameters in the widest range.That range of parameter for such heat map has been adapted to a more meaningful range.We show in Figure 6 a heat map which connects parameters variation with the blasts percentage in day +15.
This heat map leads to associate all values with two semiplanes which will be able to predict if the patient responds correctly depending on the values which match his data.Regarding Figure 6, we can approximate the influence of each drug in a patient and predict if the patient responds to the treatment correctly or not.Firstly, parameters related to the influence of drugs in cells death, δ j , can be approximated for a set of data from a specific patient in several days of evolution of the treatment.Thus, we classify these parameters according to Figure 6 in purple or white zone to predict which percentage of leukemic cells will be there in bone marrow day +15.There will be a good prognosis if leukemic compartment entails less than 0.01% (MDR) in bone marrow, white zone, while if those values are in the purple zone, the patient will not respond correctly to treatment.According to choose the minimal values which implies a good prognosis of the patient, δ P = 0.092day/mg and δ V = 2.11day/mg will be parameters used in the next results.

2.4.
A new model about bone marrow behavior in presence of treatment is proposed to study.
After proposing Model (11), we simulate the situation using data in Table 1 along with parameters associated to the δ j values study in previous sections.As proposed in Sec.1.2 it is possible to observe that leukemic cells have a higher proliferation than the rest of lymphocytes.Moreover, as said before, it is known that the blasts level in bone marrow when leukemia is detected is around 80%.We assume that treatment starts the next day and present results in Figure 7.
Therefore, according to the available data about different patients, we have been able to approximate the evolution of bone marrow cells in the presence of treatment as it is shown in Figure 7.It can be observed that the treatment affects both healthy and leukemic cells and that once the disease is eliminated, healthy cells compartments return to normal levels.Generally, in Figure 7, a leukemic bone marrow with treatment is represented.From usual cells compartments proportions, leukemia develops until leukemic cells exceed 80% of total occupancy in bone marrow.From those conditions, the therapy is applied for 37 days, and then, bone marrow normal levels are recovered.In a reduced range of days (Figure 7(b)) it can be appreciated how cells proportions are modified according to drugs administered each day.In fact, day +8 is the first aggressive day in which leukemic cells do not grow so fast, therefore, in following days this population continues its reduction.It is for that reason that the highest cells proportion in bone marrow for about 30 days of treatment is the transition compartment, C 3 , which are not influenced by treatment since they do not proliferate.Paying attention to the amount of cells (Figure 7(c)-(d)), it is possible to verify that all of the cells in bone marrow decrease gradually and finally, return to normal levels.

Discussion
Leukemia is the cancer with the highest incidence in pediatric age, particularly, Acute Lymphoblastic Leukemia and 20% of these diagnoses do not respond correctly to treatment, and as a consequence, there is a relapse.This problem is the main reason of our study and the principal question to answer is why the treatment does not work for those patients.Our intention is to find patterns in bone marrow along with medicines applied, so that, it is possible to predict relapses and improve current treatments.We approximate treatment influence in cell death to analyze the behavior of each bone marrow based on flow cytometry data throughout treatment time.Modeling current treatment in leukemia B leads to advances towards optimizing the application of medicines.After studying leukemic cells behavior it is possible to analyze and advance in this area.Available clinical data are important to determinate parameters values according to different bone marrow samples from the same patient.
First of all, we have reviewed a healthy bone marrow model in Sec.1.1.That study allows us to include the assumption of the appearance of a leukemic cell.It is found that leukemic cells have a strong development due to their properties of malignant disease.The leukemic cell develops as it has been described in Model (6), and it could be more or less accelerated depending on the inherited proliferation rate for ρ L , as shown in Figure 4.In fact, available data [30,54] of hematologic patients reinforce that the disease is diagnosed when there are about 80% leukemic cells in bone marrow.As it is mentioned in Sec.2.1, six steady states are obtained by studying Model (6).The unique stable steady state is the one that represents the total occupancy of the bone marrow by leukemic cells.The rest of solutions are unstable.There are several possible explanations for this result.On one hand, both two solutions which have negative compartments, are biologically senseless results.In spite of the fact that there are two more solutions which are non negative, they lack biological meaning.Moreover, there is a solution which reproduces healthy bone marrow behavior without leukemia, Sec.1.1.Therefore, Model (6) replicates reality and its study allows to adjust parameters related to leukemic population growth.In the case in which leukemia appears in the Pro-B stage, ρ L = ρ 1 for Model (6), the amount of leukemic cells is the same as Pro-B cells around 127 days; it is equal to Transition cells in the day 132; and Pre-B cells are reached by leukemic cells in 135 days.Simulations for model which assume the appearance of the leukemic cell in Pre-B stage, ρ L = ρ 2 , have a similar behavior, but the process is slower.In general, both cases have the same behavior.Every simulation in the presence of the leukemic clone allow for deducing that healthy and leukemic cells coexist for a period of time which depends on the case, and then, leukemic cells take up all the bone marrow capacity, therefore, there is no space for healthy cells.According to the Model (6) and its simulations, it is important to take into account that cells development is supposed without any treatment which is responsible for reducing the leukemic cells proportion.Hence, it explains why there are that suddenly changes three months later.
From the information collected in the SEHOP-PETHEMA-2013 protocol, we model the treatment behavior in the bone marrow when it is administered.Without loss of generality, we consider a leukemic cell originated in Pro-B compartment.From the first leukemic cell appearance, the proliferation of that population takes about 150 days to be diagnosed with the result that when it is detected, the disease reproduces more quickly.That is why treatment begins to be administered in the following days, thus, there are no data related to the same patient in different times without treatment.We have focused on a standard risk patient in Induction I'A which is the first phase of the treatment and it is a determining phase in the therapy.If there are signs which do not correspond to proper treatment response, the protocol for this patient is modified and associated to other risk level.Therefore, we consider a therapy based on Prednisone, Vincristine, Daunorubicine and Asparaginase, administered according to the schedule shown in Figure 2. Taking into account the total amount of drugs in the body during the treatment, we assume an exponential decay of each drug from its application in order to different references about timing and dosage drugs effects [28,59,60].In addition, we include the term δ j as the influence of each drug in cells death, to obtain the function µ which imply leukemic and healthy cells death by treatment depending on their proliferation.
When a sample from a patient is studied as a possible diagnosis of ALL is owing to the patient has previously developed some of the hallmarks of cancer [44].Generally, bone marrow first extractions indicate that leukemia takes up around 80% of bone marrow.Therefore, protocol is immediately activated.That treatment approximation implies an improvement about the perception related to drugs actions in the body.Depending on each patient, the influence of each drug will be more or less aggressive, hence, the importance of the personalized medicine [13].It is for that reason that a detailed study about drugs influences must be contemplated to study those behaviors.The aim is to provide each new patient with the best treatment.A patient is considered with good prognosis if some conditions are reached.On the treatment day +8, blood sample should show less than 10 6 blasts per milliliters of blood, on the day +15, leukemic cells must entail values smaller than 0.01% and the day +33 blasts in bone marrow must not be observed, considering that this is the moment of the transplant.
A theoretical model of the bone marrow dynamics, Model (11), is presented to study, based on previous works along with new assumptions about leukemic cells behaviour and treatment.In order to study parameters related to prednisone, we have focused on the first eight days of treatment.In those days only prednisone is administered, therefore, it is possible to approximate optimal δ P value which implies the good prognosis condition in day +8.From range of δ P possible values in Table 1, Model ( 11) is solved for 50 values in that range, obtaining amount of blasts in blood in day +8 related to each one.It is obtained a value δ P0 ≈ 0.09 from which, values imply blasts levels smaller than the limit of amount of cells, consequently, we obtain a range of δ P values in which patient response to treatment.This result imply two important aspects to consider.On one hand, a patient whose sample is extracted in day +8 could be associated to the corresponding δ P and predict the evolution of the disease.On the other hand, if we analyze blood extraction in other day, +5 for instance, we could be able to associated δ P value before day +8 and predict results in that day.Thus, it would be possible to change treatment protocol on time.
Apart from the study of the day +8 conditions, we analyze good response conditions in day +15.According to sensitivity study, prednisone and vincristine are the most significant parameters and hence, we propose a study about influence parameters related to them.On day +15, MRD should be less than 0.01% to obtain good results in patients, therefore, we search for values which lead to those leukemic cells levels when Model ( 11) is solved for those values.Leukemic compartment percentages in bone marrow on day +15 are obtained for each case and presented in Figure 6 in a heat map.As it is observed, there are two semiplanes corresponding to patients who response and patients who do not response correctly to treatment.With the purpose of minimizing both δ P and δ V , we take average minimal values which imply a good response.What is interesting in this result is the fact that the δ P value obtained by this method is the same value found by the study of the blasts levels in day +8.This finding reinforces that studies must be based on fitting influence parameters to obtain more predictive information of each patient.
Finally, simulations supported by previous results in this work are presented.Previously, parameter estimation is performed and the bone marrow behavior is studied assuming a leukemic cell appearance.The disease is diagnosed when it entails an important part of cells in bone marrow and it is treated according to protocol.As it is observed in the Figure 7 and generally, treatment application leads to patients recovery and the stabilization of healthy cells levels from fifty days after treatment finish.Treatment is applied when leukemia is diagnosed, that it is supposed in day 150 from the first leukemic cell appearance since there are about 80% of blasts in bone marrow.During eight days prednisone is applied and it is possible to observe a drop in all leukemic and healthy cells but in day +8 is more significant because vincristine and daunorubicin are added to the therapy.From day +12 in which asparaginase is administered, leukemic cells values are insignificant in comparison to healthy cells.The fact that this situation occurs agree with parameters choice which guarantee that MRD is less than 0.01% on day +15.Graphs interpretation is facilitated by the comparison between amount of cells and proportions.In those simulations, it should be highlighted that leukemic cells do not appear the following 150 days.For that reason, it is considered a partial patient recovery, not completely due to the logistic growth of leukemia.Moreover, the disease always tend to appear again and that idea explains why treatment should continue with other phases which drop the possibility of any rest of leukemic cells develops.
Being limited to available data, this study lacks accuracy related to drug influence estimated parameters.On one hand, an exhaustive study of each patient would required to be able to model the leukemia growth without treatment.To do so, two bone marrow must be extracted, to observe how cells develop and obtain the parameter related to leukemic proliferation, ρ L , in each patient.This option is unfeasible since two bone marrow extractions are rare and unnecessary, as drugs are to be administered to the child as soon as possible.Hence, we propose a future study to relate blood blasts and search for a possible correlation.Another limitation of our analysis would be the simplicity of the treatment function.Drugs evolution inside the body has a complicated and more detailed process which has been reduced and it produces a lack of information.Absorption time or action mechanism are some of the conditions to be taken into account in order to get closer to real data.Moreover, only the Induction A phase has been formulated instead of the full treatment.We have focused on this phase of treatment since effects of the next stages (Induction B and Consolidation) are not visible enough.In addition, Induction A is the most interesting phase to study owing to its critical influence throughout treatment.
Notwithstanding with these limitations, the study suggests that it is possible to approximate contemplated parameters related to treatment provided that parameters values for leukemia are corroborated.Additionally, different compartment levels in patients bone marrows with ALL corresponds to values obtained from this previous study in which a virtual patient is simulated.These findings suggest several courses of action for improving techniques used by protocols.There is a definite need for personalizing medicine.We suggest, therefore, that studies consider mechanisms of action of drugs along with an exhaustive monitoring to prove how they inhibit healthy and leukemic cells proliferation as well as how they destroy cells and the environment.We propose a clinical study about model with those conditions.Comparisons between models including new assumptions will be able to offer some approximations to specific models capable of obtaining realistic parameters.Several patients data are essential to validate the results presented.The effectiveness of that resource will allow to predict relapses in patients who have been diagnosed with ALL in a short period of time after treatment begins.The challenge now is to model full treatment, which takes into account available data and suitable parameters which lead to new research to improve current protocols and, consequently, survival expectancy in childhood ALL.
In conclusion, in this paper, a mathematical model of the appearance of leukemic cells in a bone marrow has been presented based on a mathematical model of cells development in a healthy bone marrow.Therefore, we have studied leukemic cells behavior to characterise a model describing the lymphopoiesis process when a cancer cell appears.Furthermore, we have simulated the model proposed and obtained results which coincide with biological basis of cancer.Afterwards, treatment application is included according to current protocol, meaning healthy and leukemic cells death.Dynamics in bone marrow are modeled and studied based on parameters estimation to simulate a virtual patient who responds to treatment.This model provides a basis of future action protocols of treatments of ALL.

Figure 1 .
Figure 1.Diagram of cell development along with the appearance of leukemia and its treatment.a) Immature healthy B cells from the bone marrow can develop into more mature cells as follows: In an early maturation stage, a Pro-B cell can either proliferate in its own compartment with a rate ρ 1 or progress to the next compartment (Pre-B) with rate α 1 .Analogously, a Pre-B cell can remain in a Pre-B stage with a proliferation rate ρ 2 and develop to the transition compartment with rate α 2 .Finally we consider that transition cells do not proliferate in their own compartment and they arrive to the blood flux with rate α 3 .Growth rate is regulated by the feedback signal s. b) Leukemic cells could come from a Pro-B stage or Pre-B stage, inheriting the proliferation rate ρ L of the related compartment.A regulatory signal s L controls leukemic cell development.Furthermore, for the case of ALL we assume that leukemic cells do not differentiate into a next stage and therefore γ L indicates the blood exit rate for leukemic cells into blood.c) When drug is administered, both healthy and leukemic cells are affected by the treatment, with cells dying at a rate proportional to their proliferation rate.

Figure 3 .
Figure 3.Total treatment influence in bone marrow.Percentage of total drug dose for Induction phase, following SEHOP-PETHEMA-2013 (Figure2) and with parameters from Table1.

Figure 4 .
Figure 4. Dynamics of cell compartments proportion in the bone marrow in the presence of the leukemic clone.Initial data and parameters from Table 1.(a) Case ρ L = ρ 1 = 0.6931 day −1 .Leukemic cells, L, (solid purple line) have a logistic growth and reach 80% around day 150.At the same time, Pro-B cells (C 1 , solid blue line), Pre-B cells (C 2 , dash-dotted blue line) and Transition cells (C 3 , dotted blue line) decrease.(b) Case ρ L = ρ 2 = 0.4621 day −1 .Dynamics is similar to the first case but leukemia is detected around the day 223.

Figure 5 .
Figure 5. Variation of δ P values to analyze results in day +8 of treatment.As presented in Table 1, δ P ∈ [1/60, 1/6] = [δ Pmin , δ Pmax ].Five values of δ P are highlighted: the midpoint between the maximum and minimum values δ P0 ≈ 0.09, δ P1 ≈ 0.05 the middle between δ Pmin and δ P0 , and δ P2 ≈ 0.13 the center of the interval [δ P0 , δ Pmax ] .(a) Model (11) is solved for 50 values of prednisone influence.For each δ P value, amount of blasts in day +8 is obtained and associated with a patient who responses (green points) or not (gray points) to prednisone.(b) Leukemic cells evolution depending on several δ P values.Green-dotted lines denote patients who response to prednisone, while gray-dashed lines represent bone marrow dynamics which do not reach any response to prednisone from data in day +8 (yellow vertical line).Horizontal solid red line represents the limit of blasts assumed by SEHOP-PETHEMA-2013 protocol to consider a good response of less than 2.3 × 10 10 blasts in bone marrow.

Figure 6 .
Figure 6.Blasts percentage in bone marrow day +15 of treatment.We consider 21 values up to 0.167day/mg for Prednisone and 21 values up to 4.22day/mg for Vincristine.For each values combination, Model (11) is solved for parameters in Table 1 with δ D = 2.5/30day/mg, δ A = 2.5/10 4 day/U .The value in each box corresponds to the percentage of leukemic cells in day +15 depending on values δ P and δ V .

Figure 7 .
Figure 7. Dynamics of cell compartments in the bone marrow in the presence of the leukemic clone along with therapy.Initial data, parameters related to leukemia and parameters about treatment are those presented in Table 1 except for δ P = 0.092day/mg, δ V = 2.11day/mg, δ D = 2.5/30day/mg, δ A = 2.5/10 4 day/U and ρ L = ρ 1 = 0.6931 day −1 case which we consider without loss of generality.(a) Cells proportion evolution in bone marrow from the first leukemic cell appearance along with the treatment application.(b) Bone marrow dynamics centered between day 150 and day 180, corresponding to +0 and +30 from therapy, respectively.(c) Compartment values measured in amount of cells in ten first days of treatment.It is observed that day +8 is the most aggressive of them.(d) Treatment finishes in day +37.The extinction of leukemic cells implies healthy cells back to normal.