Parameter identification of hematopoiesis mathematical model - periodic chronic myelogenous leukemia.

Periodic chronic myelogenous leukemia (PCML) is a dynamic hematopoietic disease which causes oscillations of circulating leukocytes, platelets and reticulocytes. Mathematical modeling is an invaluable tool to help in predicting hematopoiesis behavior. In this paper we modify the existing models based on improving the parameters of the model. Also more parameters are estimated regarding the proposed model. It is our major intention to construct a physiological model which can map major identified mechanisms of leukopoiesis to provide a deeper insight into this complex biological process. In the proposed model the leukocytes line has been modeled more precisely. In fact, precursor cells have been considered as two separate groups: proliferating precursor cells and non-proliferating precursor cells. As a result, more parameters have appeared in the model and identifying the new parameters has resulted in a better fit of clinical data and the data extracted from the model for both platelets and leukocytes. That is, the new model describes the leukocytes and platelets of the system in a way that is closer to clinical data, so the proposed model can be more useful for predicting the behavior of leukocytes and platelets for PCML disease. Compared with the previous works, it is shown that the new model has a better fit of the quantitative data on leukocytes and platelets.

Introduction

Leukemia is a kind of white blood cell cancer. It may be categorized based on the characteristics of the disease (acute or chronic). A patient with acute leukemia usually dies within several months, while a patient with chronic leukemia may live for years.

The diagnosis of chronic myelogenous leukemia (CML) is usually based on detection of the Philadelphia chromosome. This abnormality arises as the result of a problem in cell division within the bone marrow. It is a type of translocation, meaning that genetic material from two different chromosomes switches places. In this particular abnormality, chromosomes 9 and 22 are involved. CML arises in a bone marrow stem cell that is the precursor to all the types of blood cells (white blood cells, red blood cells and platelets). This hypothesis is derived from two pieces of evidence. First, in CML the Philadelphia chromosome is found in all of the hematopoietic lineages which come from the hematopoietic stem cell. Second, periodic oscillations are observed in three types of blood cells in one rare sort of CML, which leads to the mentioned hypothesis.

To provide a quantitative basis for prediction of behavior of stem cells and blood cells under hematopoietic diseases, several models have been developed. To our knowledge, Mackey, in 1978 [1] proposed the first model of hematopoietic stem cell dynamics, which has been developed by many authors, including Haurie and Mackey [2], who focused on some hematological diseases which originate from stem cells and specified differentiation and maturation processes in hematopoiesis, Bernard and Mackey [3], who concentrated on the production and oscillations of white blood cells, observed in patients with cyclic neutropenia, and Pujo-Menjouet and Mackey [4], who focused on applications of the model proposed in [3] to chronic myelogenous leukemia. Coljin and Mackey [5, 6] extracted a global model of the hematopoietic system that concentrates on either stem cell dynamics or the leukocyte, erythrocyte and platelet cell lines. Foley and Mackey [7] proposed a delay differential equation model for the regulation of neutrophil production under the effects of G-CSF. Models in the form of partial differential equations have been proposed by Scholz, Engel and Loeffler [8] that simulate the damage caused by cytotoxic drugs on the stem cells and leukocytes. The hematopoietic system has also been formulated as partial differential equations (PDE) by Ostby et al. [9] and ordinary differential equations (ODE) by Panetta et al. [10], Friberg et al. [11], Scholz et al. [8] and Shochat et al. [12].

Recently, a description of dynamics of chronic and acute myeloid leukemias was proposed by Adimy et al. [13, 14]. In [13], a compartmental model with discrete delays was suggested. The conditions for local stability of steady states and for Hopf bifurcation due to delay parameters were also given. In [14], a model with distributed delays was considered.

The goal of this paper is to identify some unknown parameters of the model and improve the discrepancy between the model and clinical data.

The model of the hematopoietic system

Differential equations (partial, ordinary or delay), stochastic processes, Boolean networks, etc. are some mathematical tools used in hematological modeling. The desired level of description of the model often determines the mathematical tools. Since delays arise in biological systems because of properties inherent to the different processes (time lag due to maturation, transmission of an impulse, etc.), delay differential equations are of great importance in describing biological systems such as the hematopoietic system.

Previously, leukocyte [15], erythrocyte [16] and platelet [17] models have been proposed separately. Mackey in [18] connected all these three distinct models. The model describes four distinct compartments: the hematopoietic stem cells (HSC) and the circulating leukocytes, platelets and erythrocytes. The stem cells are pluripotential as well as self-renewing, and are differentiable into the leukocytes, erythrocytes or platelets. Fig. 1 shows the full model of the system represented by Colijn and Macky [18]. The model considers the hematopoietic system as a set of compartments: the pluripotential, non-proliferating stem cells, leukocytes, erythrocytes and platelets are denoted by Q, N, R and P in unit of cells/kg, respectively. K (N), K (R) and K (P) represent the differentiation rate of HSC into these three lines in units of days–1, respectively. These differentiation rates are a monotone decreasing function of the numbers of circulating cells of the relevant type (negative feedback). More discussion of the forms of K (N), K (R) and K (P) will follow [18].

Fig. 1

A graphical representation of the full model of hematopoiesis including the HSC and three differentiated cell lines (leukocytes, erythrocytes and platelets) extracted from Mackey and Coljin 2005

A portion of non-proliferating stem cells differentiates into leukocyte, erythrocyte and platelet cell lines at the rate of K (N), K (R) and K (P), respectively. The remaining portion proliferates at the rate of β (Q) (in units of days–1). In the proliferation compartment the cells divide taking a time of τ (units of days). So, in the equation describing Q, the number of cells passing the proliferation phase is twice the number of those entering this phase a time τ ago. But the cells traversing the proliferation compartment are not quite twice those entering, because some random losses occur at a rate of γ (in units of days–1) over the proliferating time. So, the number of non-proliferating cells entering per unit time from the proliferation compartment is 2e – β(Q(t – τ ). In this paper the notation X is used for X (t – τ) hereinafter. Thus, the delay differential equation describing Q is (1).

In each leukocyte, erythrocyte and platelet cell line there is a dimensionless amplification parameter A , A and A , respectively, which represents several stages of cell divisions.

We use the equations related to red blood cells and platelets without any changes from [18], but for the leukocyte line we take some ideas from the model represented by Foley and Mackey [19], which is mainly about modeling the effects of G-CSF on leukocytes and stem cells. We consider the proliferative and non-proliferative population of precursor cells of leukocytes separately. This separation leads to some parameters being added to the model. We assume that all these new parameters are fixed and find their best values so that the model has a better fit of the quantitative data on leukocytes and platelets. Fig. 2 provides an overview of the modified model structure.

Fig. 2

A graphical representation of the full model of hematopoiesis including the HSC and three differentiated cell lines (leukocytes, erythrocytes and platelets) modifying the model of leucocyte line

Data fitting based on Coljin and Mackey (2005)

Data fitting based on modified model

In the leukocyte line, the upper box and lower box are the compartment of proliferating precursor cells and non-proliferating precursor cells, respectively. The leukocytes which enter the circulating leukocyte compartment after amplification and maturation are denoted by N. They are randomly lost at a rate γ due to apoptosis. γ is the rate of random loss of non-proliferating precursor cells while γ represents the rate of random loss of proliferating precursor cells. We assumed that the parameters γ , A , K , γ and γ are independent of G-CSF. τ is non-proliferative neutrophil precursors maturation time and τ is proliferative neutrophil precursors maturation time. The equation governing the leukocytes is then with satisfying is the velocity that cells in the non-proliferating precursor compartment age with. An increase in V leads to decreased transit time τ.

Simulation and fitting

We extracted clinical data for platelets and leukocytes from [20] and tried to find the parameters of the leukocyte line such that the model fits real data more precisely. We have used B. Ermentrout's implementation of the Marquar-dt-Levenberg method, described by Press et al. [21]. Our method to do so is minimizing the function where y are the observed data points and y are the simulated data points. We use MATLAB software for all of our calculations. In this paper, hereinafter LS represents the least square value . Considering proliferating and non-proliferating precursor compartments in the leukocyte line and estimating the parameters γ , τ , τ and V displayed in Table 1 improves the LS value for both platelets and leukocytes data in comparison with Mackey's (2005) work (Table 2 and Table 3).

Table 1

Estimation of parameters of the model

Parameter name	Estimated value	Unit
γ NN	0.26	days–1
τ NN	3.27	days
γ NP	0.25	days–1
τ NP	1.2	days
v n	6	–
θ 2	0.31
A N	953690	–
A P	46.97 × 104	–

Table 2

Comparison of the least square value for leukocyte and platelet data between existing model and improved one

	Least square value for leukocytes	Least square value for platelets
Mackey's model	4.80 × 1010	1.93 × 1011
Our model	3.88 × 1010	1.69 × 1011

Table 3

Comparison of the difference of the mean values for leukocyte and platelet data between Mackey's work and our research

	Difference of mean values of model data and real data for leukocytes	Difference of mean values of model data and real data for platelets
Mackey's model	3.44 × 109	1.47 × 1010
Our model	0.1 × 109	1.15 × 1010

Results

Although the parameters γ , γ , τ and τ are just involved in leukocyte dynamics, a change in their value affects the count of both leukocytes and platelets.

High values of τ (proliferative neutrophil precursors maturation time) cause a decrease in the mean value of the leukocytes.

Variation in the parameter V has the least effect on the model in comparison with the other parameters.

The dynamics of the system are more sensitive to variations in τ in comparison with τ . Looking at equ. (2) more precisely supports this fact. As can be observed the parameter K depends on τ in the mentioned equation.

The lower the value of the parameter τ , the higher the mean value of the leukocytes is.

The higher the value of the parameter θ2 is, the more the amount of the peak and the mean value of N and P are.

The parameters τ and θ2 affect the shape of the model at its first and last oscillates.

In conclusions, this article identifies some parameters of the hematopoietic system so that the model has a better fit of the quantitative data on leukocytes and platelets in comparison with fitness done by the others. The effects of mentioned parameters on behavior of leukocytes and platelets are also investigated.

Discussion

The major intention of this paper is to construct a physiological model which can map major identified mechanisms of leukopoiesis to provide a deeper insight into this complex biological process. A new model of leukopoiesis has been proposed. In fact, precursor cells of the leukocyte line have been considered as two separate groups: proliferating precursor cells and non-proliferating precursor cells. As a result, more parameters have appeared in the model and identifying the new parameters has resulted in a better fit of clinical data and the data extracted from the model for both platelets and leukocytes. Using the new model can lead to better prediction of the behavior of leukocytes and platelets for PCML disease.