ARTICLE IN PRESS Journal of Theoretical Biology 237 (2005) 117–132 www.elsevier.com/locate/yjtbi A mathematical model of hematopoiesis—I. Periodic chronic myelogenous leukemia Caroline Colijna, Michael C. Mackeyb, a Department of Mathematics and Centre for Nonlinear Dynamics, McGill University, 3655 Promenade Sir William Osler, Montreal, Que., Canada H3G 1Y6 b Departments of Physiology, Physics and Mathematics and Centre for Nonlinear Dynamics, McGill University, 3655 Promenade Sir William Osler, Montreal, Que., Canada, H3G 1Y6 Received 16 December 2004; received in revised form 22 March 2005; accepted 30 March 2005 Available online 21 June 2005 Abstract Periodic chronic myelogenous leukemia (PCML) is an interesting dynamical disease of the hematopoietic system in which oscillating levels of circulating leukocytes, platelets and/or reticulocytes are observed. Typically all of these three differentiated cell types have the same oscillation period, but the relation of the oscillation mean and amplitude to the normal levels is variable. Given the appearance of the abnormal Philadelphia chromosome in all of the nucleated progeny of the hematopoietic stem cells (HSCs), the most parsimonious conclusion is that chronic myelogenous leukemia, and its periodic variant, arise from derangements partially involving the dynamics of the stem cells. Here, we have synthesized several previous mathematical models of HSC dynamics, and models for the regulation of neutrophils, platelets and erythrocytes into a comprehensive model for the regulation of the hematopoietic system. Based on estimates of parameters for a typical normal human, we have systematically explored the changes in some of these parameters necessary to account for the quantitative data on leukocyte, platelet and reticulocyte cycling in 11 patients with PCML. Our results indicate that the critical model parameter changes required to simulate the PCML patient data are an increase in the amplification in the leukocyte line, an increase in the differentiation rate from the stem cell compartment into the leukocyte line, and the rate of apoptosis in the stem cell compartment. Our model system is particularly sensitive to changes in stem cell apoptosis rates, suggesting that changes in the numbers of proliferating stem cells may be important in generating PCML. r 2005 Elsevier Ltd. All rights reserved. Keywords: Leukemia; Mathematical model 1. Introduction All blood cells arise from a common origin in the bone marrow, the hematopoietic stem cells (HSCs). These stem cells can differentiate into one of three major cell lines: the leukocytes, the platelets, and the erythrocytes. The exact details of how the numbers of circulating cells of each type are regulated remain Corresponding author. Tel.: +1 514 398 4336; fax: +1 514 398 7452. E-mail addresses: ccolijn@cnd.mcgill.ca (C. Colijn), mackey@cnd.mcgill.ca (M.C. Mackey). 0022-5193/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2005.03.033 somewhat obscure, though the broad outlines are clear. For the red blood cells, the cytokine erythropoietin (EPO) mediates a negative feedback loop that helps to regulate erythrocyte production (Adamson, 1974; Silva et al., 1996). A decrease in numbers of circulating erythrocytes leads to a decrease in tissue pO2 levels, which in turns triggers a production and release of EPO by renal macula cells. This elevation of EPO increases the net production of primitive erythroid precursors, and finally an increase in the numbers of circulating erythrocytes takes place some days later. In the production of white blood cells, granulocyte ARTICLE IN PRESS 118 C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 colony-stimulating factor (G-CSF) has been demonstrated to be very important and it too mediates a negative feedback loop (Price et al., 1996). As for the erythrocytes and leukocytes, the production of platelets involves a negative feedback loop, related to the production and maturation of megakaryocytes. These are platelet precursor cells that give rise to 1000–5000 platelets each (Beutler et al., 1995). Thrombopoietin (TPO), which is known to be involved in platelet production, mediates this negative feedback and TPO has been shown to have effects on other cell lines as well (Ratajczak et al., 1997; Tanimukai et al., 1997), which could imply that the three lines are not fully independent. The regulation of the HSC population is also of importance. Previous mathematical models of HSC dynamics have been based on a G 0 cell cycle model, and thus considered actively proliferating stem cells and quiescent, non-proliferating cells in the so-called G 0 phase (Mackey, 1978, 1979a, 2000; Fowler and Mackey, 2002; Pujo-Menjouet et al., 2001). There is typically a delay in mathematical models of the stem cells, reflecting the non-zero time that it takes the cells to complete the proliferative phase of the cell cycle which includes DNA synthesis, mitosis, and cytokinesis. However, the parameters determining the rate of proliferation in the stem cell compartment, which depends ultimately on the current number of stem cells, are not well-known (PujoMenjouet et al., 2001; Ogawa, 1993; Lemishka et al., 1986). Since stem cell oscillations are thought to drive oscillations in several periodic hematological diseases (Haurie et al., 1998), developing a better knowledge of stem cell dynamics is important. There is significant autonomy between the three cell lines, as evidenced by their distinct responses to varying circulating cell numbers in each line; i.e. to distinct demands for production. Several hematological diseases provide motivation for the mathematical study of hematopoietic regulation, as well as providing the opportunity to study early hematopoietic regulation and its effects on the relationships between the cell lines (Haurie et al., 1998). These disorders include (but are not limited to) cyclical neutropenia (CN), periodic chronic myelogenous leukemia (PCML), cyclical thrombocytopenia, and periodic hemolytic anemia. The first two involve oscillations in all circulating hematological cells with the same period. The regulatory mechanisms in early hematopoiesis are not as well understood as the peripheral hematopoietic regulatory elements. Thus, periodic hematological disorders in which oscillations occur in more than one cell line are of considerable interest since the most parsimonious explanation of oscillatory pathologies involving more than one cell type is that they arise in the HSC. Of particular interest here will be periodic CML, in which oscillations occur primarily in the leukocytes, but may also occur in the platelets and in some cases the reticulocytes (Fortin and Mackey, 1999). Leukemia is a disease of the blood and bone marrow, characterized by large and apparently uncontrolled accumulation of blood cells. Myelogenous leukemia involves the myelocytes, which are granular leukocyte precursors. In acute leukemia, the disease progresses rapidly, and produces large numbers of immature nonfunctioning white blood cells. In chronic forms of leukemia, the cells produced are initially more functional. CML is associated with a chromosomal abnormality, known as the Philadelphia chromosome, that occurs in all cell lineages in about 90% of cases. This abnormality consists of a transposition of parts of chromosomes 9 and 22, and results in the formation of the Bcr–Abl fusion protein (ODwyer et al., 2000). This protein is thought to be responsible for the dysfunctional regulation of myelocyte growth and other features of CML (Melo, 1996). (See Grignani (1985) for more detail about CML.) In periodic chronic myelogenous leukemia (PCML) the leukocyte count varies periodically, typically between values of 30 and 200
109 cells L1 . This is far above the normal value of 6
109 cells L1 . The variation occurs with a period in the range of 40–80 days, which is very long in comparison with the maturation and lifespan of the stem cells and the leukocytes. In addition, oscillations may occur in the platelets and occasionally also in the reticulocytes, and in these cases the platelet and reticulocyte periods are the same as the leukocyte periods (Henderson et al., 1996; Fortin and Mackey, 1999). It has been argued that this, in addition to the occurrence of the Philadelphia chromosome in all differentiated lineages, is indicative of the stem cell origins of PCML (Fortin and Mackey, 1999). In both PCML and CN, the hypothesis that oscillations originate in the stem cells is related to the fact that oscillations occur in different lines. However, in many earlier mathematical models, only one cell line, or one line coupled to the stem cells, is represented. In this context, it is not possible to examine the effects of a destabilization in one line or in the stem cell compartment on whole system. For example, while PujoMenjouet et al. (2001) explored how long-period oscillations (as seen in PCML) could arise within the context of a G 0 stem cell model, the stem cell model alone could not predict whether the leukocytes and platelets would oscillate at the levels observed in PCML. Similarly, Bernard et al. (2003) were able to duplicate various features of CN with an integrated mathematical model of the HSC and peripheral neutrophil control. However, since their model did not include platelet and erythrocyte regulation it is unknown if their simulated neutropenic conditions would be consistent with observed platelet and erythrocyte data in CN. ARTICLE IN PRESS 119 C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 This paper presents a model framework within which these questions can be addressed. We utilize a G 0 type model for the stem cell compartment, based on the work of Pujo-Menjouet et al. (2001), and couple it to a leukocyte model based on that of Bernard et al. (2003), as well as to two simplified models representing peripheral platelet and erythrocyte regulation. In Section 2.1, the model is described, including the various feedback and control functions that capture the essence of the control exercised by G-CSF, EPO and TPO. Section 2.2 deals with the parameter estimation for the model in the normal (non-pathophysiological) case. In Section 3, the published PCML data with which the model is compared are briefly described. This is followed in Section 4 with a description of the techniques used to fit the parameters in the model to the PCML data, and the results of the fitting approaches along with their interpretation. 2. A model of the hematopoietic system
τS QτS a = τS a=0 Q(t) s (t,a) γS κN κR AP AR AN κP τPM τNM τRM N(T) τPS τRS 2.1. Model formulation In previous work leukocyte (Hearn et al., 1998; Haurie and Mackey, 2000; Bernard et al., 2003), erythrocyte (Mackey, 1979b; Bélair et al., 1995; Mahaffy et al., 1998a) and platelet (Bélair and Mackey, 1987; Santillan et al., 2000) dynamics have been modelled separately, with the goal of building quantitative understanding of cellular production within the context of periodic hematological disorders. The goal of the present work is to link these models together, connecting models for the three distinct cell lines to a mathematical model of the stem cell population (Mackey, 1978). The present model thus has four distinct compartments representing the HSC and the circulating leukocytes, platelets and erythrocytes. The stem cells are pluripotential and self-renewing, and can differentiate into the leukocyte, erythrocyte or platelet lines. Alternatively, the stem cells may re-enter the proliferative phase of the stem cell compartment. The stem cell and leukocyte compartments are modelled using the stem cell model connected to a neutrophil population as in Bernard et al. (2003). The platelet and erythrocyte compartments are simplified approximations of earlier modelling efforts. The full model is illustrated in Fig. 1. The pluripotential, non-proliferating stem cells are denoted by Q (in units of 106 cells kg1 , see Fig. 1). The circulating leukocytes, erythrocytes and platelets are denoted N (units 109 cells kg1 ), R (units 1011 cells kg1 ) and P (units 1010 cells kg1 ), respectively. The rates of differentiation into these three lines are given by kN ðNÞ, kR ðRÞ and kP ðPÞ, respectively. These differentiation rates have units of days1 . The fact that these rates each depend on the numbers of circulating cells of the γN γP γR P(T) R(T) Fig. 1. A cartoon representation of the full model including the HSC and the three differentiated cell lines. See the text for full details as well as the notation. relevant type encapsulates the feedback between the circulating cell numbers and the production; the feedback is always negative so when the number of circulating mature cells of a given line falls, the relevant differentiation rate k has a corresponding compensatory increase. More discussion of the forms of kN ðNÞ, kR ðRÞ and kP ðPÞ will follow. With reference to Fig. 1, first consider the stem cell portion of the model. There are two ways that the stem cells can exit the non-proliferating compartment: they can either re-enter the proliferating phase at a rate bðQÞ (units of days1 ), or they can differentiate into any of the three cell lines, at rates kN ðNÞ, kR ðRÞ and kP ðPÞ for the leukocytes, erythrocytes and platelets respectively. After re-entering the proliferating phase, the cells divide, taking a time tS (units of days) to do so. The production of Q, where Q represents the density of non-proliferating stem cells, is given by the flux of stem cells entering the proliferative phase, multiplied by two to account for cell division, and delayed by the time tS that it takes them to multiply. However, the number of cells leaving the proliferating compartment is not simply twice the number that entered it a time tS ago, because during proliferation, there may be loss ARTICLE IN PRESS 120 C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 through cell death. We represent this loss as occurring at a rate gS (units of days1 ) over the time period tS . Thus, the number of non-proliferating stem cells entering per unit time (and augmenting Q) is 2egS tS bðQðt  tS ÞÞQðt  tS Þ. (In this paper, we always use the notation QtS  Qðt  tS Þ.) The loss from the non-proliferating compartment is given by the sum of the terms representing differentiation into the three cell lines, plus the loss due to the cells re-entering the proliferative phase. Thus, the differential delay equation governing Q is dQ ¼  bðQÞQ  ðkN þ kR þ kP ÞQ dt þ 2egS tS bðQðt  tS ÞÞQðt  tS Þ. ð1Þ There are three more differential delay equations in the model, each representing one cell line. These three models for peripheral control share some common features and display some differences. In each cell line there is an amplification stage representing many stages of cell division (typically 12–18, depending on the cell line). The dimensionless amplification parameters are AN , AR and AP . Following completion of amplification through cell division, the cells traverse a maturation compartment (duration in days denoted by tiM ), where i ¼ N; R; P) and then enter the circulation. While the leukocytes die at a constant rate gN (units of days1 ), the platelets and erythrocytes enter an ‘age-structured’ compartment where they persist over times tRS and tPS , respectively, and are primarily lost to senescence. However, there are also well documented random losses of platelets and erythrocytes represented by gR and gP (again units of days1 ). Specifically considering the leukocytes, after amplification and maturation they enter the circulating leukocyte compartment (denoted by N). They are randomly lost at a rate gN . The time for maturation is denoted by tN , and the equation governing the leukocytes is then, as in Bernard et al. (2003) dN ¼ gN N þ AN kN ðN tN ÞQtN . dt (2) The platelet and erythrocyte dynamics are somewhat more complex. Indeed, in previous work (Mahaffy et al., 1998b; Santillan et al., 2000), age-structured models were used to represent their dynamics in the bone marrow and the loss via senescence of circulating erythrocytes and platelets. In the present model we will replace the age distribution simply with discrete delays in the erythrocyte and platelet lines (Mackey, 1996). The equations for the erythrocyte and platelet dynamics are dR ¼  gR R þ AR fkR ðRtRM ÞQtRM dt  egR tRS kR ðRtRM þtRS ÞQtRM þtRS g ð3Þ and dP ¼  gP P þ AP fkP ðPtPM ÞQtPM dt  egP tPS kP ðPtPM þtPS ÞQtPM þtPS g. ð4Þ In Eq. (3), the term gR R reflects a random loss of erythrocytes at rate gR . The positive term in Eq. (3) represents the amplification of the cells entering the erythrocyte line from the stem cell compartment, minus the loss that accrues in the age-structured compartment (represented by the lower box in the erythrocyte line; see Fig. 1) over the time tRS . Eq. (4) for the platelet compartment has exactly the same structure as Eq. (3), and all the terms in Eq. (4) are analogous to those in Eq. (3). The lifetime tRS ¼ 120 days is very long, and consequently oscillations in the HSC or the leukocyte and platelet lines may not be transmitted visibly to the erythrocyte line due to the damping of the long timescale. However, oscillations may be visible in the erythrocyte precursors (the reticulocytes), and for this reason, we have divided the erythrocyte compartment into two sub-compartments for the purposes of comparing model simulations with existing reticulocyte data. Eq. (3) is thus replaced by two equations dL ¼  gR L þ AR fkR ðRtRM ÞQtRM dt  egR tret kR ðRtRM þtret ÞQtRM þtret g ð5Þ and dR ¼  gR R þ AR fegR tret kR ðRtRM þtret ÞQtRM þtret dt  egR tsum kR ðRtsum ÞQtsum g, ð6Þ L now represents the reticulocytes, and R still represents the erythrocytes, which are simply the ‘later’ component of the simplified age-structured box of Fig. 1. The delay time tsum is the total time spent after maturation, tsum ¼ 120 days. It remains to discuss the forms of the feedback functions kN , kP and kR , as well as the rate bðQÞ at which the stem cells return to the proliferative phase. In Bernard et al. (2003), the re-entry rate to the proliferative phase was given as bðQÞ ¼ k0 ys2 . ys2 þ Qs (7) The form of this equation is very common in enzyme kinetic relations (see the discussion in Bernard et al. (2003)), and offers a good match to existing data as well as providing convenient analytical properties. We will use the parameters given in Bernard et al. (2003) for this function (see Table 2). ARTICLE IN PRESS C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 For the leukocyte feedback function kN , we use the same function as was used in Bernard et al. (2003) kN ðNÞ ¼ f 0 yn1 . yn1 þ N n (8) 121 which is exactly analogous to Eq. (9). By the same reasoning, we use kR ðRÞ / E / 1 1 þ K r Rm (12) and write This function is derived from the assumption that the differentiation rate is proportional to the number of bound G-CSF receptors on the cells. While the concentration of G-CSF changes with time, for a given leukocyte concentration, the G-CSF concentration reaches equilibrium on a time-scale of hours, whereas the time-scale of interest in our study of periodic hematological disorders is on the order of weeks. We therefore make the additional assumption that the G-CSF concentration changes rapidly enough that it reaches equilibrium with respect to the number of circulating white blood cells. This is what allows us to write the rate of differentiation of the stem cells into the leukocyte line as a function of N rather than as a function of the G-CSF concentration; see Bernard et al. (2003) for a derivation of Eq. (8). Megakaryocytes are platelet precursors, each of which produces, on average, 1000–5000 platelets (Beutler et al., 1995). In Santillan et al. (2000) the authors assumed that the number of megakaryocytes of age 0, entering from the stem cell compartment, was simply proportional to the concentration of TPO. In our model, we preserve the approximation that the differentiation is proportional to the thrombopoiesis concentration, but assume, as in the case of G-CSF, that the TPO is in dynamic equilibrium with the platelet numbers, and hence that dT=dt ¼ 0, where T is the concentration of TPO. From Santillan et al. (2000), dT a ¼  kT, dt 1 þ KPr kR ðRÞ ¼ 1 1 þ K p Pr dQ ¼ bðQÞQ  ðkN þ kR þ kP ÞQ þ 2egS tS bðQtS ÞQtS , dt dN ¼ gN N þ AN kN ðN tN ÞQtN , dt dL ¼  gR L þ AR fkR ðRtRM ÞQtRM dt  egR tret kR ðRtRM þtret ÞQtRM þtret g, dR ¼  gR R þ AR fegR tret kR ðRtRM þtret ÞQtRM þtret dt  egR tsum kR ðRtsum ÞQtsum g, dP ¼  gP P þ AP fkP ðPtPM ÞQtPM dt  egP tPS kP ðPtPM þtPS ÞQtPM þtPS g, bðQÞ ¼ k0 k̄p . 1 þ K p Pr (11) The rate kR of differentiation into the erythrocyte line is modelled similarly. Indeed, in Mahaffy et al. (1998a), it was assumed that the number of new precursor cells entering the erythrocyte line was linearly proportional to the concentration of EPO. There, the authors develop the equation dE a ¼  kE, dt 1 þ KRm ys2 , ys2 þ Qs kN ðNÞ ¼ f 0 yn1 , yn1 þ N n kP ðPÞ ¼ k̄p , 1 þ K p Pr kR ðRÞ ¼ k̄r . 1 þ K r Rm (15) For reference, the delays in the system, and the time values composing them, are listed in Table 1. and can write kP ðPÞ ¼ ð14Þ where (9) (10) (13) In summary, the equations comprising the model are so we have kP ðPÞ / T / k̄r . 1 þ K r Rm Table 1 Delays in the model Delay Description tS tN tRM tret tRS tsum tPM tPS Stem cell proliferation time Leukocyte maturation time Erythrocyte maturation time Reticulocyte maturation time Erythrocyte aging time to senescence tRM þ tret þ tRS Platelet maturation time Platelet aging time to senescence ARTICLE IN PRESS 122 C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 2.2. Parameter estimation Throughout this discussion, the subscript indicates a steady-state value. The steady-state values are obtained (and related to other parameters) by setting the rates of change in Eq. (14) to zero: Q_ ¼ N_ ¼ L_ ¼ R_ ¼ P_  0. 2.2.1. The stem cell compartment The HSC parameters are the death rate gS in the proliferating phase, the time tS that cells spend in the proliferative phase, the parameters k0 , ys and s occurring in the Hill function bðQÞ, and the size of the non-proliferative HSC compartment Q at steady state, Q . The value of Q is derived from various data. Abkowitz et al. (1988) give a value of 8 stem cells per 105 nucleated bone marrow cells in cats, while in mice it is estimated that there are 1–50 stem cells per 105 nucleated bone marrow cells (Boggs et al., 1982; Micklem et al., 1987). Since there are approximately 1:4
1010 nucleated bone marrow cells per kg of body mass in mice, this would imply that there are between 1:4
105 and 7
106 HSC per kg body mass in mice, or 1:1
106 stem cells per kg body mass in cats. We have adopted the latter value so Q ¼ 1:1
106 stem cells per kg body mass. Mackey (2000) has estimated in the mouse that gS 2 ð0:069; 0:228Þ day1 , that the steady-state rate of re-entry from the non-proliferative phase to the proliferative phase is bðQ Þ 2 ð0:020; 0:053Þ day1 , tS 2 ð1:41; 4:25Þ days, and that the steady-state rate of differentiation into all cell lines from Q is between 0.010 and 0:024 day1 . The duration of the proliferative phase of the HSC cell cycle has been estimated as tS ¼ 2:8 days by Abkowitz et al. (1988) in the cat, and between 1.4 and 4.3 days by Mackey (2000) in the mouse. Bernard et al. (2003) give the rate of apoptosis in the stem cell compartment gS ¼ 0:01–0:2 days1 . In their discussion of long-period oscillation in a G 0 stem cell model, Pujo-Menjouet et al. (2001) give 0:18 days1 , and in Mackey (2000) the range for gS is given as 0.07–0:23 days1 . We began our simulations with gS ¼ 0:18. Abkowitz et al. (1988) estimated the mean proliferation rate of HSCs in mice to be bðQ Þ ¼ 0:06 days1 which, in combination with the other Hill function parameters given in Bernard et al. (2003), gives y2 ¼ 0:5 ð
106 Þ cell kg1 and k0 ¼ 3:0 days1 . These parameters are taken from a variety of experiments (Lyons, 1999; Oostendorp, 1995), in conjunction with fitting a visual inspection of the Bernard et al. (2003) model to experimental data. The parameter s may be interpreted as the number of ligand molecules that are required in the activation of a receptor site, which in the case of HSCs is not well understood. Bernard et al. (2003) estimate that the value of s is between 2 and 4. Based on our results and the discussion of oscillations in the stem cell model analysed by Pujo-Menjouet et al. (2001), we have used s ¼ 4. 2.2.2. The leukocyte compartment For the leukocyte control dynamics, we need to estimate the stationary state (mean) leukocyte count, N , as well as the loss rate gN , the leukocyte maturation time tN and the Hill function parameters f 0 , y1 and n. The loss rate of circulating leukocytes is gN ’ 2:4 days1 (Deubelbeiss et al., 1975). Novak and Necas (1994) give the daily leukocyte production as 2:4
109 cells kg1 day1 . This gives N ’ 0:5
109 cells kg1 day1 for humans. We have used a range of N ’ 0:35–0:7
109 cells kg1 day. Lebowitz and Rubinow (1969) estimate a maturation time for neutrophils of tN ¼ 3–10 days, following the analysis of Hearn et al. (1998) we take tN ¼ 3:5 days for the leukocyte compartment. The amplification parameter in the leukocyte line, AN , has been estimated by Mackey (2000) as 300,000 based on the data in Cheshier et al. (1999), and as 171,000 based on the data of Abkowitz et al. (2000). The previous model of Bernard et al. (2003) gives a range of 0–100,000 based on Novak and Necas (1994). We have taken a value of AN ¼ 75; 200, which is within this range. We will adjust the units of N so that the numerical value of AN we use in the steady-state is 75; the units of N are 109 cells kg1 body mass while the units of Q (the quiescent stem cells) are 106 cells kg1 . The Hill function parameters used in our model are similar to those given in Bernard et al. (2003), where the range of the parameter f 0 was given to be f 0 ¼ 0:4–1:5 days1 based on experiments in Sachs and Lotem (1994) and Ward et al. (1999) on the effects of G-CSF administration on cell production. Bernard et al. estimate y1 ¼ 0:36
108 cells kg1 from data given in Ward et al. (1999); Lotem and Sachs (1998); Akbarzadeh et al. (2002), and n ¼ 1. In our model, we compute kN ¼ kN ðN Þ ¼ gN N AN Q from Eq. (2) evaluated at steady-state conditions. Having found kN , we compute f 0 f 0 ¼ kN y1 þ N ¼ 0:4. y1 2.2.3. The erythrocyte compartment Mahaffy et al. (1998a) have already estimated most of the required parameters. The mean (steady state) erythrocyte number is 3:5
1011 cells kg1 body weight. The maturation time tRM is approximately 6 days, and the total lifespan, corresponding to the time spent in the lower compartment in our model, tsum , is 120 days. The reticulocyte lifespan tret is 2.8 days (Beutler et al., 1995). ARTICLE IN PRESS 123 C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 The loss rate from the erythrocyte compartment was estimated to be gR ¼ 0:001 days1 . Based on a comparison with the data presented in Erslev (1991) and Mahaffy et al. (1998a) estimated the Hill function parameters to be m ¼ 6:96 and K r ¼ 0:0382. We used the same approach as in the leukocyte compartment to estimate kR based on the steady-state values, and k̄r from kR kR ¼ gR R AR Q ð1  egR tRS Þ and k̄r ¼ kR ð1 þ K R Rm Þ. In these calculations we use AR ¼ 563; 000 as the amplification in the erythrocyte line. Beutler et al. (1995) estimated 10–12 divisions between the proerythroblast stage and the mature red blood cells, while Novak and Necas (1994) estimate 9 doubling stages between the stem cells and the proerythroblasts, for a total of 19–21 divisions. Our amplification factor is thus approximately 219 . 2.2.4. The platelet compartment For platelets, Santillan et al. (2000) have estimated many of the relevant parameters. Thus, the mean platelet count is P ¼ 2:5
108 cells mL1 blood. Using the fact that a 70 kg adult has 6 L of blood, we find P ¼ 2:14
1010 cells kg1 body mass. By fitting the TPO concentration versus the platelet count data given in Kuter (1996), Santillan et al. (2000) estimate the Hill function parameters to be K P ¼ 31:18 in units of P , and r ¼ 1:29. We convert K P to kg cell1 (as we are not going to work in units of P ) so that K P ¼ 11:66 ð
1010 cells kg1 Þ1 . However, due to the range of observed steady-state values of P, P 2 ð1:4; 4Þ
1010 cells kg1 body mass, we require some flexibility in P , which is implemented in our model by allowing K P to change: K P 2 ð3; 13Þ
1010 cells kg1 body mass. Finally, Santillan et al. (2000) estimated that gP ¼ 0:15, tPM ¼ 7 days and tPS ¼ 9:5 days. The estimation of kP in the platelet case has been done somewhat differently than in the erythrocyte and leukocyte compartments, because it is necessary to set kP to ensure that a non-zero steady state will exist in the model. In other words, if the efflux from the resting phase of the stem cell compartment is too large, the small influx to the resting-phase stem cell compartment from the proliferating phase will be insufficient and the stem cells numbers will decrease exponentially to zero. For this reason, we ensure that (with reference to Eq. (1)), ðkN þ kR þ kP þ b Þ þ 2egS tS b ¼ 0 (16) for the values of N , R and P found from the literature. k̄p is found in the same manner as k̄r , and we use a value of AP ¼ 282; 000 for the platelet amplification factor. Beutler et al. (1995) reports that approximately 35; 000 platelets are produced per mL of blood per day, which, translated into the context of our model, would mean that AP times the steady-state flux out of the stem cell compartment should be approximately this amount, or 2:57
109 cells kg1 day1 . The steadystate flux out of the stem cell compartment is kP Q . We can, therefore, solve for AP kP  2500. Since it turns out from Eq. (16) that kP is on the order of 102 , this is consistent with our amplification factor. A summary of the parameters used in the model is given in Table 2. Table 2 Parameters used for the steady state Parameter name Value used Stem cell compartment 1.1 Q gS 0.1 2.8 tS k0 3.0 0.5 y2 s 4 Unit Sources
106 cells kg1 Days1 Days Days1
106 cells kg1 (none) 1 1 1,2 1 1 1 Leukocyte compartment 3.55–7 N gN 2.4 3.5 tN AN 75,200 (75.2) f0 0.40 0.36 y1 n 1
108 cells kg1 Days1 Days (none) 2,3 1,4,5 1 1,3 Days1
108 cells kg1 (none) (calculated) 1 1 Erythrocyte compartment 3.5 R gR 0.001 6 tRM tsum 120 tret 2.8 AR 563,000 k̄r 1.1 0.0382 Kr m 6.96
1011 cells kg1 Days1 Days Days Days (none) Days1 ð
1011 cells kg1 Þm (none) 6 6 6 6 3 3,7 (calculated) 6 6 Platelet compartment P 2.14 gP 0.15 7 tPM tPS 9.5 AP 282,000 1.17 k̄p 11.66 Kp r 1.29
1010 cells kg1 Days1 Days Days (none) Days1 ð
1010 cells kg1 Þr (none) 8 8 8 8 3 (calculated) 8 8 Sources: 1 ¼ Bernard et al. (2003), 2 ¼Abkowitz et al. (1988), 3 ¼Beutler et al. (1995), 4 ¼Deubelbeiss et al. (1975), 5 ¼Haurie and Mackey (2000), 6 ¼Mahaffy et al. (1998b), 7 ¼Novak and Necas (1994), 8 ¼Santillan et al. (2000). ARTICLE IN PRESS 124 C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 3. Data We have analysed data taken from published studies of periodic CML. These data were previously analysed for significant periodicity using Lomb periodogram (Lomb, 1976) techniques by Fortin and Mackey (1999). Each primary study presented time series of patient leukocyte counts, and some also provided platelet and reticulocyte data. In all but two of the studies the patient was untreated. Fortin and Mackey (1999) digitized the published data, converting the result to PostScript format using PhotoShop 4.0. The times and values of each data point were then digitally extracted. Upon Lomb periodogram analysis, Fortin and Mackey (1999) found 9 sets of data showing significant periodicity in the leukocytes and platelets (at significance levels pp0:05, in addition to several cases in which the platelets were reported to be oscillating, but the platelet data were not published. We have analysed the 9 data sets in which leukocytes and platelets were oscillating; in two of these, oscillations in the reticulocytes were observed as well. We have also analysed two data sets for which only leukocyte data were available. These 11 data sets are the data which we are interested in simulating with our model. The original studies in which the data are published are: Kennedy (1970), Rodriguez and Lutcher (1976), Gatti et al. (1973), Vodopick et al. (1972), Yamauchi and Ide (1992) and Umemura et al. (1986), Morley et al. (1967), Iizuka et al. (1984), Chikkappa et al. (1976) and Delobel et al. (1973). In the study by Kennedy (1970), the patient was undergoing constant hydroxyurea treatment; because the dosage was constant it is not believed that the treatment initiated the cyclicity. The patient described in Rodriguez and Lutcher (1976) was receiving bulsafan treatment; however, there was evidence of cycling before treatment began. None of the other patients were undergoing treatment. Only leukocyte data are given in Delobel et al. (1973) and Gatti et al. (1973), while platelet data are given in all the other cases; only Chikkappa et al. (1976) and Iizuka et al. (1984) give reticulocyte data. There is considerable variation not only in the periods, which range from 40 to 80 days, but also in the leukocyte means, for which the range spans a factor of 6, between 20
109 and 120
109 cells L1 . 4. Simulation and fitting There are several possibilities for the simulation of PCML in the context of this model. It has been suggested, for example, that CML originates in the HSC compartment, both because the Philadelphia chromosome is found in is all of the hematopoietic lineages, as mentioned in Section 1, and because in PCML all cell lines oscillate with the same period. It is of interest to discover to what extent the observed relationship between the cell lines is duplicated in the model. The periods in PCML are also of interest, being between 40 and 80 days. This is significantly longer than the periods in other oscillating hematopoietic diseases like CN, and is two orders of magnitude longer than the duration of the proliferating phase of the cell cycle in the stem cells. A variety of approaches to simulating PCML, based on existing notions of the mechanisms of the disease, were initially tried. For example, based on the hypothesis that oscillations are caused by changes in the dynamics of the HSC, we changed the apoptosis rate gs in the stem cell compartment to determine whether this might induce oscillations like those seen in PCML. However, with the parameters as estimated in Section 2.2, the oscillations generated had periods in the 20 day range, which is well below the entire observed range in PCML. In addition, the leukocytes did not oscillate significantly above the steady state in any of these simulations in contrast to the clinical data. In PCML there is a large proliferation of leukocytes, so we simulated this excessive proliferation by increasing the parameter AN from its steady-state value of approximately 75 to much higher values of 1000–5000. However, this also did not result in leukocyte numbers at the observed levels, and in general did not lead to a bifurcation in the stability of the steady state, and sustained oscillations. Other attempts to simulate PCML based on intuition were not successful in duplicating the observed behavior of the system. However, it must be recognized that the parameter space for this model is large, and intuitive explorations are capable of describing the behavior in only a tiny subset of the parameter space. Thus, some notion of where in the parameter space to begin the simulations is essential. In Pujo-Menjouet et al. (2001), some exploration of long-period oscillations in a G 0 model of HSCs has been carried out. In fact, Eq. (1) reduces to the equation used in PujoMenjouet et al. (2001) in the case that all the differentiation rates k are constant. Pujo-Menjouet et al. (2001) examined how to create long-period oscillations in the model, both at smaller values of s (the exponent in the Hill function bðQÞ), and as s ! 1. Their results (for s ¼ 4) provided a starting point at which our model generated long-period oscillations. It is also necessary to specify the initial history of the variables, since there are delays in the system. In our simulations, the histories over times t 2 ½240; 120 were set to the healthy, default values. Then, while t 2 ½120; 0, the system was integrated, allowing transient behavior to disappear; in most cases, since of course the parameters were significantly different from the steady-state (healthy) parameters, oscillations ARTICLE IN PRESS C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 occurred during this time. Simulations were compared to observed data starting at t ¼ 0. Thus, the histories immediately prior to comparison with the data are not constant functions, but vary depending on the parameter choices for the particular simulation. We use and compare two distinct approaches to the fitting of the model to the observed data. In the first approach, we have used B. Ermentrout’s implementation of the Marquardt–Levenberg method, described in Press et al. (1992), which is included with the software xppaut. The Marquardt–Levenberg method minimizes the function w2 ¼ N X ðys  yi Þ2 , s2i i¼1 (17) where yi are the observed data points, ys are the simulated data points, si are the estimated standard deviations of the noise in the observed data points, and N is the total number of points available. However, we do not have access to values of si for each data point, and have had to assume that si ¼ s ¼ constant (18) over each data set. In our case, xppaut sets s2 to be the standard deviation of the data for the leukocytes, platelets and reticulocytes separately, to provide scaling information to the algorithm. Leukocytes Data S¼ N X ðys  yi Þ2 for each data type, i.e. for the leukocytes, platelets and reticulocytes separately. The primary advantage of this method of fitting is that it quite rapidly converges to a nearby minimum in w2 , using a variable step size in the parameter adjustment to aid convergence. However, it has the disadvantage that in a parameter space in which there are many local ‘energy’ (w2 ) minima, the algorithm will converge to a local minimum and will not be able to escape that minimum even if there is a global minimum elsewhere in parameter space. Because of the problem of multiple local minima, it was necessary to make every effort to begin the Marquardt–Levenberg convergence algorithm at a point in parameter space from which it was likely that a ‘good’ minimum would be reached. This process was carried out by hand, using our putative understanding of the possible biological processes involved in PCML as well as experience with generating oscillations in the model. For example, the parameter AN was set higher than the steady state, while tS , k0 and gS were set within estimated ranges to generate long-period oscillations. The other parameters were initiated at their steady-state values. Full Simulation 25 20 20 20 15 15 15 10 10 10 5 5 5 0 200 Data 0 0 200 Full Simulation 25 0 20 20 15 15 15 10 10 10 5 5 5 0 200 0 Sampled Simulation 0 200 Sampled Simulation 25 20 0 (19) i¼1 25 25 Platelets In this case, minimizing w2 is simply the same as minimizing the sum of squares 25 0 125 0 0 200 0 200 Fig. 2. An example of the results of the Marquardt–Levenberg fitting to the Vodopick et al. (1972) data, Case 2. Leukocytes are in units of 109 kg1 , while platelets are in units of 1010 kg1 . In the left-hand panels are the observed data; the centre panels shows the full simulation and the right-hand panel plots show the data and the simulation sub-sampled at the observed data times reported by Vodopick et al. (1972). ARTICLE IN PRESS 126 C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 Leukocytes Data 20 20 20 10 10 10 0 Platelets Sampled Simulation Full Simulation 0 200 400 0 0 0 200 400 15 15 15 10 10 10 5 5 5 0 0 0 200 400 200 400 0 200 400 0 200 400 0 0 200 400 0.5 0.5 Reticulocytes 0.5 0 0 0 200 400 0 0 200 400 0 Fig. 3. An example of the results of the Marquardt–Levenberg fitting to the Iizuka et al. (1984) data. Leukocytes are in units of 109 kg1 , while platelets are in units of 1010 kg1 . In the left-hand panels are the observed data; the centre panels shows the full simulation and the right-hand panel plots show the data and the simulation sub-sampled at the observed data times reported by Iizuka et al. (1984). Table 3 Parameter estimates for the 11 patients of this study based on the Marquardt–Levenberg method Default Iizuka Chikkappa Delobel Gatti Kennedy Morley Rodriquez Umemura Vodopick (1) Vodopick (2) Yamauchi AN f0 AP kP gS k0 75.00 0.40 28.00 1.17 0.18 3.00 167.7 277.5 283.6 169.0 190.1 291.7 335.9 153.3 236.7 327.1 138.8 (1.23) (1.19) (2.87) (0.25) (0.31) (10.08) (19.52) (33.11) (0.04) (0.04) (64.16) 5.81 5.42 6.70 5.81 3.50 30.09 3.54 1.39 5.74 9.61 1.83 (1.17) (0.65) (0.99) (0.10) (0.13) (9.41) (9.09) (16.55) (0.01) (0.11) (109.56) 39.17 59.44 28.14 28.21 13.96 30.14 1.57 16.52 45.56 65.50 17.94 (0.00) (1.06) (0.00) (0.00) (51.80) (0.59) (854.00) (48.31) (0.63) (3.38) (47.39) 0.45 0.55 0.44 0.44 0.32 0.56 0.06 1.46 0.83 0.93 0.92 (0.00) (0.51) (0.00) (0.00) (24.71) (0.26) (280.39) (6.02) (0.20) (0.67) (29.40) 0.18 0.18 0.18 0.19 0.19 0.13 0.18 0.19 0.18 0.18 0.19 (0.02) (0.02) (0.02) (0.02) (0.02) (0.02) (0.02) (0.01) (0.02) (0.02) (0.01) 3.01 3.30 3.01 3.04 2.30 4.17 1.54 2.47 2.93 2.51 2.03 (0.99) (1.57) (0.52) (0.00) (0.00) (0.03) (7.83) (10.98) (0.04) (0.41) (14.56) The top line of the table gives the default values from Table 2 for comparison. Error estimates, in percent, are shown in brackets. Sample plots showing the data and the fit are given in Figs. 2 and 3. As shown in the plots, there are longperiod oscillations forming an envelope around much higher-frequency oscillations, particularly in the leukocytes. This was a typical feature of most of our simulations, and has also occurred in previous models based on the G 0 description of the stem cells (PujoMenjouet et al., 2001). It is therefore not an artifact of particular parameter sets chosen by the fitting algorithms. It should be noted that the data are quite sparsely sampled, over significant time periods. For example, in the Vodopick et al. (1972) data, the first few points are at t ¼ 10; 40; 54 . . . days, where near the end of the data set there are points at every 2 or 3 days. The sparse sampling means that over wide time ranges, it is not possible for the fitting algorithm to distinguish between simulations which contain higher-frequency oscillations and those which do not. In addition, there are segments of the observed data which do seem to show rapid ARTICLE IN PRESS C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 changes; see for example Iizuka et al. (1984) and Gatti et al. (1973). The parameters found using the Marquardt–Levenberg algorithm are given in Table 3. In this table, the values in parentheses are the errors in the estimates, given as a percentage of the corresponding parameter value. These estimates are based on the covariance matrices computed by the fit procedure. In these results, several trends of importance emerge. First, the values of the amplification parameter AN in the leukocyte line are all well above the steady-state value of 75. The rate of differentiation into the leukocyte line, f 0 , is also significantly greater than the steady-state value, implying that in PCML not only is there more amplification of the cellular precursors entering the leukocyte line, but that there are increased numbers of stem cells entering the leukocyte pathway. Both of these conclusions drawn from the model fitting procedure are consistent with clinical observations, giving us added confidence in the fitting results. It was found in simulations that changing the parameter y2 , part of the Hill function regulating the rate of proliferation of the quiescent stem cells, does not significantly improve the fit. y2 has not been included in the table, nor in subsequent fitting procedures. However, k0 , which scales the overall proliferation rate of the stem cells, has varied between 1.5 and 4 in these results, consistent with the (Bernard et al., 2003) range of k0 of 2.0–10.0. The parameter estimates for gS with the Marquardt –Levenberg method require some discussion, as in some cases the system was so sensitive to small changes in gS that either convergence was apparently impossible, or given convergence the estimated error was so small that it seemed unreasonable. Where gS could not be included in the fitting procedure, it was left fixed at a value near the default value, with the requirement that the fitting procedure gave reasonable convergence in the other parameters. It is for this reason that the cited errors in Table 3 for gS are very small. Also for this reason, it is not meaningful to draw conclusions about the relationship between the results for gS and the apoptosis in the stem cell compartment in PCML. This limitation does not apply, however, to the second fitting approach, discussed below. In the second parameter estimation scheme, we implemented a simulated annealing procedure. Simulated annealing is an approach to optimization which simulates physical cooling of metals. Originally derived from computer simulation of physical equilibrium (Metropolis et al., 1953), the algorithm explores the parameter space looking for a minimum in an ‘energy function’, which in this case is the sum of squares as given in Eq. (19). However, unlike other optimization procedures, in simulated annealing the current parameter set may be altered to one which gives a higher energy, when one yielding a lower energy cannot be 127 found. In this way, the system is perturbed out of local minima. The cooling, or annealing, is done by gradually decreasing the likelihood that such an upward move will be accepted. Simulated annealing is known to be highly advantageous when the energy function contains many local minima. Our energy function is given by vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uM  s uX ðL  Li Þ2 ðPs  Pi Þ2  i þ i E¼t ¯2 L P¯2 i¼1  s ðN̄  N̄Þ2 ðP̄s  P̄Þ2 þd þ Þ 2 2 N̄ P̄  ðvarðN s Þ  varðNÞÞ2 ðvarðPs Þ  varðPÞÞ2 þ . þd varðNÞ2 varðPÞ2 ð20Þ Here, N and P refer to leukocytes and platelets, while the superscript s indicates simulation, and the lack of a superscript indicates observed data. The bars indicate that the mean has been taken, and ‘var’ refers to the variance. M is the number of data points available. In the cases where platelet data are not available, the terms involving platelets are simply left out. Where reticulocyte data are available, the energy function is vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uM  s uX ðN  Li Þ2 ðPs  Pi Þ2  ðLs  Li Þ2 i E¼t þ i 2 þ i 2 2 N̄ P̄ L̄ i¼1  s ðN̄  N̄Þ2 ðP̄s  P̄Þ2 ðL̄s  L̄Þ2 þ Þþ Þ þd 2 2 2 N̄ P̄ L̄  ðvarðN s Þ  varðNÞÞ2 ðvarðPs Þ  varðPÞÞ2 þd þ varðNÞ2 varðPÞ2 þ ðvarðLs Þ  varðLÞÞ2 , varðLÞ2 ð21Þ where the notation is as it is in Eq. (20), and L refers to the reticulocytes. This energy function is simpler than it may appear. The first and dominant term is the usual sum of squares, with the leukocytes and platelets each normalized by their means to prevent the fit being dominated by one of them. (In other words, the presence of N̄, P̄ and L̄ in Eqs. (20) and (21) provides scaling information.) The square root deforms this function monotonically, making the energy landscape less steep. Such monotonic deformations are common in simulated annealing approaches. The remaining terms involving the means and variances of the leukocytes and platelets are scaled by a small parameter d; these terms encourage oscillatory solutions (recall that all of our data have been previously tested for periodicity using the Lomb algorithm (Fortin and Mackey, 1999)). Without these terms, non-oscillatory solutions would frequently emerge from the annealing process, because a constant ARTICLE IN PRESS 128 C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 solution can have a lower sum-of-squares energy function than an out-of-phase but oscillatory solution, or an oscillatory solution at a slightly different frequency Adding these terms to the energy function successfully reduced the likelihood that the simulations would converge to a steady-state solution. The details of our simulated annealing procedure are as follows. We use the Metropolis acceptance rule (Metropolis et al., 1953) for determining whether a proposed perturbation of the parameters will be accepted, namely that it is accepted if the move will decrease the energy. If none of the moves tested decrease the energy, then an uphill move is accepted with probability P ¼ eDE=T , where T is the ‘‘temperature’’. Thus, as T decreases, the probability that an uphill move will be accepted decreases. In determining how and when to change the temperature, we have used a geometric cooling schedule, in which the temperature is decreased according to the rule T n ¼ T 0 an , (22) where a 2 ð0:99520:999Þ. This is done after the system undergoes a random walk of a given length, where the moves in the random walk are generated by perturbations of the relevant parameters and selected by the Metropolis algorithm. The initial temperature T 0 is chosen using the criterion that a given fraction (usually half) of the attempted uphill moves will be accepted in the initial random walk. The framework for our program is due to Salamon et al. (2002), which can be consulted for a more detailed discussion of simulated annealing. We have applied this approach to the PCML data, comparing the simulations for leukocytes and platelets, as well as reticulocytes where applicable. The results of the parameter estimates, along with errors in percentages, are given in Table 4. In terms of the qualitative features of the simulations, the high-frequency oscilla- tions observed in the Marquardt–Levenberg fits appeared again in the simulated annealing fits, which is not surprising as the system has no way of selecting for the removal of these high-frequency oscillations with the sparsely sampled data available. Figs. 4 and 5 show examples of the simulated annealing output, both with fine sampling (left-hand plots) and with sampling at the same times as the observed data (right-hand plots). Note that the highfrequency oscillations are not visible at the lower sampling rate corresponding to the rate at which the clinical data were collected. As in the results for the Marquardt–Levenberg procedure, the most prominent result is the high overall rate of amplification in the leukocyte line. In some cases these are an order of magnitude above the default value of 75, though in those cases the errors are typically quite high as well. Here, error estimation was done by finding how far from the estimated parameter one had to go for the energy function to increase by 10 percent. This value is based on the observed similarity of results within this range of energy values. In the simulated annealing results, the values of f 0 range widely above the default values as they do for the Marquardt–Levenberg results. In both simulations there is a wide range of values of the parameter AP , which is necessary if the simulations are to match the platelet data. The rate of apoptosis in the stem cell compartment, gS , is usually found to be around the default range of 0.16–0.2, but occasionally was altered by the annealing procedure to a much lower value. This occurred in the fits to the data of Chikkappa et al. (1976) and Yamauchi and Ide (1992), and it is interesting to note that in these cases f 0 was closer to its default value of 0.4. Increasing f 0 and AN and perhaps also decreasing gS seem to be ways that the model can produce long-period oscillations at the leukocyte levels seen in PCML. Physiologically, this is Table 4 Parameter estimates for the 11 patients of this study based on the simulated annealing method Default Iizuka Chikkappa Delobel Gatti Kennedy Morley Rodriquez Umemura Vodopick (1) Vodopick (2) Yamauchi AN f0 AP kP gS k0 75.00 0.40 28.00 1.17 0.18 3.00 423.8 314.4 187.6 251.3 124.5 763.0 388.6 68.8 316.3 359.9 532.4 (4.51) (39.35) (4.74) (6.87) (9.32) (2.60) (2.59) (0.81) (0.68) (19.6) (42.0) 2.22 12.74 19.91 8.44 1.98 20.48 5.57 3.76 2.75 9.07 0.53 (12.5) (6.97) (0.77) (3.60) (4.89) (2.03) (12.6) (0.5) (2.60) (8.30) (27.6) 24.30 4.21 34.16 37.57 5.57 26.88 8.48 13.85 55.97 19.07 8.94 (10.00) (12.39) (16.85) (7.39) (12.90) (19.9) (13.8) (2.60) (1.15) (6.21) (4.92) 1.07 4.55 0.83 0.53 3.53 0.75 0.05 1.74 1.58 5.06 1.93 (10.4) (20.15) (4.99) (5.62) (11.34) (2.27) (28.1) (1.3) (5.00) (13.2) (5.6) 0.17 0.12 0.17 0.15 0.13 0.18 0.18 0.17 0.16 0.14 0.04 (2.42) (5.39) (0.23) (5.00) (3.20) (0.08) (7.38) (0.04) (2.50) (5.08) (19.38) 3.09 2.81 2.89 1.02 0.98 3.77 1.63 2.32 3.88 1.30 1.05 The top line of the table gives the default values from Table 2 for comparison. Error estimates, in percent, are shown in brackets. (7.81) (13.9) (0.47) (8.16) (7.98) (4.95) (12.2) (0.2) (1.9) (8.4) (14.3) ARTICLE IN PRESS 129 C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 Leukocytes Full Simulation Data 25 25 20 20 20 15 15 15 10 10 10 5 5 5 0 0 0 200 0 0 Data 25 Platelets Sampled Simulation 25 200 0 Sampled Simulation Full Simulation 25 25 20 20 20 15 15 15 10 10 10 5 5 5 0 0 0 200 200 0 0 200 0 200 Fig. 4. An example of the results of the simulated annealing fitting to the Vodopick et al. (1972) data, Case 2. Leukocytes are in units of 109 kg1 , while platelets are in units of 1010 kg1 . In the left-hand panels are the observed data; the centre panels shows the full simulation and the right-hand panel plots show the data and the simulation sub-sampled at the observed data times reported by Vodopick et al. (1972). Leukocytes Data 20 20 20 10 10 10 0 0 0 Platelets Sampled Simulation Full Simulation 200 400 0 0 200 400 15 15 15 10 10 10 5 5 5 0 0 0 200 400 200 400 0 200 400 0 200 400 0 0 200 400 0.5 0.5 Reticulocytes 0.5 0 0 0 0 200 400 0 0 200 400 Fig. 5. An example of the results of the simulated annealing fitting to the Iizuka et al. (1984) patient data. Leukocytes are in units of 109 kg1 , while platelets are in units of 1010 kg1 and reticulocytes are in units of 1011 kg1 . On the left hand side are the observed clinical data; the centre panels show the full simulation and the right-hand plots show the data and the simulation sub-sampled at the observed data times reported by Iizuka et al. (1984). ARTICLE IN PRESS 130 C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 consistent with excess proliferation in the leukocyte line as well at the stem cell level. Because this model does not explicitly represent each stage of leukocyte development, perhaps a decrease in gS is as close as the system can come to representing excess proliferation of stem cells that are differentiating into the leukocyte line. Our simulations also allow us to comment on the likely origin of the oscillations observed in PCML. Using the parameter sets presented in Tables 3 and 4, we ran simulations in a model system containing only the stem cells and leukocytes; in other words, the coupling to the platelet and erythrocyte compartments was removed. In most cases, this ‘‘decoupled’’ system still oscillates, implying that the parameter changes in the stem cell and leukocyte compartments are sufficient to destabilize the steady state and lead to oscillations. In the remaining (non-oscillating) cases, the decoupled system was close, in parameter space, to one that did oscillate. (These nearby parameter sets were found either by examining the results of the simulated annealing procedure that generated the parameter sets of Table 4, or by running the annealing procedure with the decoupled model beginning from the parameter sets in Table 4.) In addition, when the leukocyte compartment is then removed, leaving only the stem cell compartment in the model, the oscillations disappear, and the stem cell population decreases to zero. Thus, we feel that we can assert with some certainty that the oscillatory nature of PCML is due to a bifurcation in the combined dynamics of the regulation of the HSC and the neutrophil control. The nature of this bifurcation will be studied in a future, more mathematically oriented, communication. Overall, the quality of the simulated annealing fits is notably better than that of the Marquardt–Levenberg fits, almost certainly because less prior knowledge (of where to begin the fit) is required, and the simulated annealing approach is less likely to converge to a local minimum. In order to compare the relative quality of the fits from the two procedures, we computed a w2 -like function for all the fits; this function was the same as the energy function defined in Eqs. (20) and (21), except that the terms involving the means and variances of the simulation and data (those terms that were scaled by d in the energy function) were omitted. Thus, we computed vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uM  s uX ðN  N i Þ2 ðPs  Pi Þ2 ðLs  Li Þ2  i w¼t (23) þ i 2 þ i 2 2 N̄ P̄ L̄ i¼1 (leaving out the platelet and reticulocyte terms when data were not available). This w function is not the same as the usual w2 ; notably, the means of the data are in our function in place of the standard deviations usually used in computing w2 . However, since there are indications that fluctuations of cell numbers are not normally distributed even in normal circumstances (Perazzo et al., 2000), there is no significant advantage in using the usual w2 function for comparison. There is an advantage in using the mean values, as we have done, for this ensures that the platelets and the neutrophils are ‘weighted’ equally. Using the function given in Eq. (23), we find that the simulated annealing fits have lower w values by an order of magnitude, on average, than the Marquardt–Levenberg fit results. This quantitatively confirms the observation that the simulations generated by the simulated annealing fit parameters ‘look better’ when compared to the data. 5. Discussion and conclusions PCML is a strange cycling disease of the hematopoietic system in which oscillating levels of circulating leukocytes, platelets and/or reticulocytes are observed (Fortin and Mackey, 1999; Haurie et al., 1998). Typically it is observed that all of these three differentiated cell types have the same oscillation period, but that the relation of the oscillation mean and amplitude to the normal levels is variable. Given the common appearance of the abnormal Philadelphia chromosome in all of the nucleated progeny of the HSCs, the simplest conclusion is that CML, and its periodic variant, arise from derangements partially involving the dynamics of the stem cells. Here, we have synthesized several previous mathematical models of hematopoietic stem cell dynamics (Mackey, 1978, 1979a, 2000; Pujo-Menjouet et al., 2001), and models for the regulation of neutrophils (Hearn et al., 1998; Haurie and Mackey, 2000; Bernard et al., 2003), platelets (Bélair and Mackey, 1987; Santillan et al., 2000) and erythrocytes (Mackey, 1979b; Bélair et al., 1995; Mahaffy et al., 1998a) into a single model for the regulation of the hematopoietic system. Based on estimates of parameters for a typical normal human, we have systematically explored the changes in some of these parameters necessary to account for the quantitative data on leukocyte, platelet and reticulocyte cycling in 11 patients with PCML. We have used two different fitting procedures (the Marquardt–Levenberg procedure as well as simulated annealing), and there is an interesting and encouraging concordance between the results of both. In the most general terms, both lead to qualitative and quantitative agreement with the published data on PCML in reproducing the period, amplitudes and mean values of the oscillating cell types as well as the relative phase differences between them. It is significant that the present model is capable of duplicating the overall features of the coupled oscillations in several cell lines. The platelet oscillations in the fitted simulations, particularly with the simulated annealing method, ARTICLE IN PRESS C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 matched the observed data fairly well, including the relative phase. Where reticulocyte data were available, the oscillations were also matched. This indicates that the model is capable of generating the observed coupling between the cell lines. Note that there is often a high-frequency oscillation on top of these long time periods, but there is nothing in the available data that would rule out the possibility of high-frequency oscillations in cell numbers. However, it is likely that if the model included distributions of maturation delay times, rather than single discrete delays, some smoothing of this high-frequency output might occur. Our results indicate that the oscillatory nature of PCML is probably generated through a bifurcation in the dynamics of the coupled HSC compartment and the regulation of differentiated leukocytes. Based on the simulations, the critical model parameter changes required to simulate the periodic chronic myelogenous leukemia patient data are the amplification in the leukocyte line (AN ), the differentiation rate from the stem cell compartment into the leukocyte line (f 0 ), and the rate of apoptosis in the stem cell compartment (gS ). Our model system is also particularly sensitive to changes in gS , suggesting that changes in the numbers of proliferating stem cells may be important in generating PCML. However, there was no indication that the rate of re-entry from the quiescent stem cell compartment to the proliferating phase stem cells was involved. The platelet parameters seem to be of secondary importance, and changes in these seem to merely adjust levels of mean platelet levels and amplitude of the excursions from the mean. Acknowledgements This work was supported by MITACS (Canada) and the Natural Sciences and Engineering Research Council of Canada. References Abkowitz, J., Golinelli, D., Harrison, D., Guttorp, P., 2000. The in vivo kinetics of murine hemopoietic stem cells. Blood 96, 3399–3405. Abkowitz, J.L., Holly, R.D., Hammond, W.P., 1988. Cyclic hematopoiesis in dogs: studies of erythroid burst forming cells confirm an early stem cell defect. Exp. Hematol. 16, 941–945. Adamson, J., 1974. The relationship of erythropoietin and iron metabolism to red blood cell production in humans. Semin. Oncol. 21, 9–15. Akbarzadeh, S., Golinelli, D., harrison, D., Guttorp, P., 2002. Tyrosine residues of the granulocyte colony-stimulating factor receptor transmit proliferation and differentiation signals in murine bone marrow cells. Blood 99, 879–887. Bélair, J., Mackey, M., 1987. A model for the regulation of mammalian platelet. Ann. NY Acad. Sci. 504, 280–282. 131 Bélair, J., Mackey, M., Mahaffy, J., 1995. Age-structured and twodelay models for erythropoiesis. Math. Biosci. 128, 317–346. Bernard, S., Belair, J., Mackey, M., 2003. Oscillations in cyclical neutropenia: new evidence based on mathematical modeling. J. Theor. Biol. 223. Beutler, E., Lichtman, M.A., Coller, B.S., Kipps, T., 1995. Williams Hematology. McGraw-Hill, New York. Boggs, D., Boggs, S., Saxe, D., Gress, L., Canfield, D., 1982. Hematopoietic stem cells with high proliferative potential: assay of their concentration in marrow by the frequency and duration of cure of W/Wv mice. J. Clin. Invest. 70, 242–253. Cheshier, S., Morrison, S., Liao, X., Weissman, I., 1999. In vivo proliferation and cell cycle kinetics of long term self-renewing haematopoietic stem cells. Proc. Natl Acad. Sci. USA 96, 3120–3125. Chikkappa, G., Borner, A.H.B., Chanana, A.D., Cronkite, E.P., Ohl, S., Pavelec, M., Robertso, J.S., 1976. Periodic oscillation of blood leukocytes, platelets, and reticulocytes in a patient with chronic myelocytic leukemia. Blood 47, 1023–1030. Delobel, J., Charbord, P., Passa, P., Bernard, J., 1973. Evolution cyclique spontanée de la leucocytose dans un cas de leucémie myéloide chronique. Nouv. Rev. franc Hémat. 13, 221–228. Deubelbeiss, K.A., Dancey, J.T., Harker, L.A., Finch, C.A., 1975. Neutrophil kinetics in the dog. J. Clin. Invest. 55, 833–839. Erslev, A., 1991. Erythropoietin titers in health and disease. Semin. Hematol., 779–791. Fortin, P., Mackey, M., 1999. Periodic chronic myelogenous leukemia: Spectral analysis of blood cell counts and etiological implications. Br. J. Haematol. 104, 336–245. Fowler, A., Mackey, M., 2002. Relaxation oscillations in a class of delay differential equations. SIAM J. Appl. Math. 63, 299–323. Gatti, R.A., Robinson, W.A., Deinard, A.S., Nesbit, M., McCullough, J.J., Ballow, M., Good, R.A., 1973. Cyclic leukocytosis in chronic myelogenous leukemia. Blood 41, 771–782. Grignani, F., 1985. Chronic myelogenous leukemia. Crit. Rev. Oncol.Hematol. 4, 31–66. Haurie, C., Mackey, M., 2000. Modeling complex neutrophil dynamics in the grey collie. J. Theor. Biol. 204, 504–519. Haurie, C., Mackey, M.C., Dale, D.C., 1998. Cyclical neutropenia and other periodic hematological diseases: a review of mechanisms and mathematical models. Blood 92, 2629–2640. Hearn, T., Haurie, C., Mackey, M., 1998. Cyclical neutropenia and the peripherial control of white blood cell production. J. Theor. Biol 192, 167–181. Henderson, E.S., Lister, T.A., Greaves, M.F. (eds.), 1996. Leukemia. Saunders, London. Iizuka, Y., Horikoshi, A., Sekiya, S., Sawada, U., Ohshima, T., Amaki, I., 1984. Periodic fluctuation of leukocytes, platelets and reticulocytes in a case of chronic myelocytic leukemia: The relation between leukocyte counts, cfu-c colony formation, csa and cia. Acta Haematol. Japan 47, 71–79. Kennedy, B.J., 1970. Cyclic leukocyte oscillations in chronic myelogenous leukemia during hydroxyurea therapy. Blood 35, 751–760. Kuter, D., 1996. The physiology of platelet production. Stem Cells 14, 88–101. Lebowitz, J., Rubinow, S., 1969. Grain count distributions in labeled cell populations. J. Theor. Biol. 23, 99–123. Lemishka, I.R., Raulet, D.H., Mulligan, R.C., 1986. Developmental potential and dynamic behavior of hemopoietic stem cells. Cell 45, 917. Lomb, N.R., 1976. Least-squares frequency analysis of unequally spaced data. Astrophys. Space Sci. 39, 447–462. Lotem, J., Sachs, L., 1998. In vivo control of differentiation of myeloid leukemic cells by recombinant granulocyte–macrophage colonystimulating factor and interleukin 3. Blood 71, 375–382. Lyons, A.B., 1999. Divided we stand: tracing cell proliferation with carboxy uorescein diacetate syccinimidyl ester. Immunol. Cell Biol. 77. ARTICLE IN PRESS 132 C. Colijn, M.C. Mackey / Journal of Theoretical Biology 237 (2005) 117–132 Mackey, M.C., 1978. A unified hypothesis for the origin of aplastic anemia and periodic haematopoiesis. Blood 51, 941–956. Mackey, M.C., 1979a. Dynamic haematological disorders of stem cell origin. In: Vassileva-Popova, J.G., Jensen, E.V. (Eds.), Biophysical and Biochemical Information Transfer in Recognition. Plenum Publishing Corp, New York, pp. 373–409. Mackey, M.C., 1979b. Periodic auto-immune hemolytic anemia: an induced dynamical disease. Bull. Math. Biol. 41, 829–834. Mackey, M.C., 1996. Mathematical models of hematopoietic cell replication and control. In: Othmer, H., Adler, F., Lewis, M., Dallon, J. (Eds.), The Art of Mathematical Modeling: Case Studies in Ecology, Physiology and Biofluids. Prentice-Hall, New York, pp. 149–178. Mackey, M.C., 2000. Cell kinetic status of haematopoietic stem cells. Cell Prolif. 34, 71–83. Mahaffy, J., Belair, J., Mackey, M., 1998a. Hematopoietic model with moving boundary condition and state dependent delay: applications in erythropoiesis. J. Theor. Biol. 190, 135–146. Mahaffy, J.M., Bélair, J., Mackey, M., 1998b. Hematopoietic model with moving boundary condition and state dependent delay. J. Theor. Biol. 190, 135–146. Melo, J., 1996. The diversity of bcr–abl fusion proteins and their relationship to leukemia phenotype. Blood 88, 2375. Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A., Teller, E., 1953. Equation of state calculations by fast computing machines. J. Chem. Phys. 21, 1087–1092. Micklem, H., Lennon, J., Ansell, J., Gray., R.A., 1987. Numbers and dispersion of repopulating hematopoietic cell clones in radiation chimeras as functions of injected cell dose. Exp. Hematol. 15(3), 251–257. Morley, A.A., Baikie, A., Galton, D., 1967. Cyclic leukocytosis as evidence for retention of normal homeostatic control in chronic granulocytic leukaemia. Lancet, 1320–1322. Novak, J.P., Necas, E., 1994. Proliferation differentiation pathways of murine haematopoiesis: correlation of lineage fluxes. Cell. Prolif. 27, 597–633. ODwyer, M., Druker, B.J., Mauro, M., Talpaz, M., Resta, D., Peng, B., Buchdunger, E., Ford, J., Reese, S.F., Capdeville, R., Sawyers, C.L., 2000. Sti571: a tyrosine kinase inhibitor for the treatment of cml. Ann. Oncol. 11, 155. Ogawa, M., 1993. Differentiation and proliferation of hematopoietic stem cells. Blood 81, 2844–2853. Oostendorp, R.A., Audet, J., Eaves, C.J., 1995. High-resolution tracking of cell division suggests similar cell cycle kinetics of hematopoietic stem cells stimulated in vitro and in vivo. Blood, 855–862. Perazzo, C., Fernandez, E., Chialvo, P., 2000. Large scale-invariant fluctuations in normal blood cell counts: a sign of criticality? Fractals 8, 279–283. Press, W., Teukolsky, S., Vetterling, W., Flannery, B., 1992. Numerical Recipes in C, second ed. Cambridge University, Cambridge. Price, T.H., Chatta, G.S., Dale, D.C., 1996. Effect of recombinant granulocyte colony stimulating factor on neutrophil kinetics in normal young and elderly humans. Blood 88, 335–340. Pujo-Menjouet, L., Bernard, S., Mackey, M., 2001. Long period oscillations in a g0 model of hematopoietic stem cells. J. Theor. Biol. 209. Ratajczak, M., Ratajczak, J., Marlicz, W., Pletcher, C., Machalinski, B., Moore, J., Hung, H., Gewirtz, A., 1997. Recombinant human thrombopoietin (TPO) stimulates erythropoiesis by inhibiting erythroid progenitor cell apoptosis. Br. J. Haematol. 98, 8–17. Rodriguez, A.R., Lutcher, C.L., 1976. Marked cyclic leukocytosis leukopenia in chronic myelogenous leukemia. Am. J. Med. 60, 1041–1047. Sachs, L., Lotem, J., 1994. The network of hematopoietic cytokines. Proc. Soc. Exp. Biol. Med. 206, 170–175. Salamon, P., Sibani, P., Frost, R., 2002. Facts, Conjectures, and Improvements for Simulated Annealing. Society for Industrial and Applied Mathematics, Philadelphia. Santillan, M., Mahaffy, J., Belair, J., Mackey, M., 2000. Regulation of platelet production: the normal response to perturbation and cyclical platelet disease. J. Theor. Biol. 206, 585–603. Silva, M., Grillot, D., Benito, A., Richard, C., Nunez, G., FernandezLuna, J., 1996. Erythropoietin can promote erythroid progenitor survival by repressing apoptosis through bcl-1 and bcl-2. Blood 88, 1576–1582. Tanimukai, S., Kimura, T., Sakabe, H., Ohmizono, Y., Kato, T., Miyazaki, H., Yamagishi, H., Sonoda, Y., 1997. Recombinant human c-Mpl ligand (thrombopoietin) not only acts on megakaryocyte progenitors, but also on erythroid and multipotential progenitors in vitro. Exp. Hematol. 25, 1025–1033. Umemura, T., Hirata, J., Kaneko, S., Nishimura, J., Motomura, S., Kozuru, M., Ibayashi, H., 1986. Periodical appearance of erythropoietin-independent erythropoiesis in chronic myelogenous leukemia with cyclic oscillation. Acta Haematol. 76, 230–234. Vodopick, H., Rupp, E.M., Edwards, C.L., Goswitt, G.A., Beauchamp, J., 1972. Spontaneous cyclic leukocytosis and thrombocytosis in chronic granulocytic leukemia. New Eng. J. Med. 286, 284–290. Ward, A.C., van Aesch, Y.M., Gits, J., Schelen, A.M., de Koning, J.P., van Leeuwen, D., Freedman, M.H., Touw, I.P., 1999. Novel point mutation in the extracellular domain of the granulocyte colonystimulating factor (G-CSF) receptor and in a case of severe congenital neutropenia hyporesponsive to G-CSF treatment. J. Exp. Med. 190, 497–508. Yamauchi, K., Ide, A., 1992. Spontaneous remission with cyclic leukocytosis in chronic myelogenous leukemia. Acta Haematol. 88, 136–138.