Research Journal of Recent Sciences ______ ______________________________ ______ ____ ___ ISSN 2277 - 2502 Vol. 1( 4 ), 1 4 - 21 , April (201 2 ) Res.J. Recent Sci. International Science Congress Association 14 Analysis o f G - CSF Treatment of C N using Fast Fourier Transform Balamuralitharan S. 1 and Rajasekaran S. 2 1 Department of Mathematics, Sri Ramanujar Engineering College, Chennai - 48, INDIA 2 Professor and COE, Department of Mathematics, BS Abdur Rahman University, Chennai - 48, INDIA Available online at: www.isca.in (Received 9 th February 2012 , revised 12 th February 2012 , accepted 18 th February 2012 ) Abstract In hi reerch er we ry o fin ou n inveie he F Fourier Trnfor oel n G - CSF Treatment of CN (Cyclicl Neuroeni)”, in eil. In hi nlyi of G - CSF reen of ‘Neuroeni’, we e ‘’ fro CN. They re rey c ollies. They are usually used to build an extended model of it. It produces the dynamics of circulating blood cells. They are found from the dogs with and without daily G - CSF therapy. It is a model which is very useful for collection of laboratory data. Th is mathematical model helps us to reproduce the large variation of data too. They occur from one dog to another. It has long term effects on the oscillations when the frequency of drug delivery is made. This model is also useful to account for the fe atures of untreated G - CSF. It is also useful for treatment of dogs with CN. Therefore this model is considered as an accomplished one. There is fitting parameters for 3 days and not for 4 dogs for estimation or evaluation. It is also essential and necessary to m odel the more samples for increase in Neutrophil amplification. The proposed interventions are practical. It may reduce the amount of G - CSF. It required potential maintenance. Sometimes, it may even improve the treatment effects too. This model gives us go od result in treatment. The changes would be practical and reduce the risk side as well as the cost of treatment in G - CSF. Key Words: FFT modeling, neutrophil, granulocyte colony - stimulating factor (G - CSF). Introduction All blood cells are derived from the hematopoietic stem cells (HSC), which are undifferentiated cells having a high proliferative potential 1. These multipotent stem cells can proliferate and mature to form all types of blood cells. Production in these cell lines is regulated by a variety of cytokines, including erythropoietin (EPO), which mediates the regulation of erythrocyte production, thrombopoietin (TPO), which regulates production of platelets, as well as granulocyte colony - stimulating factor (G - CSF), which regulates leu kocyte numbers 2. In a comprehensive mathematical model for the regulation of hematopoiesis was presented. This work was motivated by the existence of several hematological diseases that display a highly dynamic nature characterized by oscillations in one or more of the circulating cell lines. Examples of these are cyclical neutropenia, periodic chronic myelogenous leukemia, cyclical thrombocytopenia and periodic hemolytic anemia. In this chapter, we concentrate on cyclical neutropenia, a rare hematological disorder characterized by oscillations in the circulating neutrophil count. These levels fall from normal to barely detectable levels with a typical period of 19 to 21 days in humans, even though periods up to 40 days have been observed. These oscillation s in the neutrophil count are generally accompanied by oscillations with similar period in the platelets, lymphocytes and reticulocytes . Cyclical neutropenia also occurs in grey collies with periods on the order of 11 to 16 days . This animal model has prov ided extensive experimental data that has enriched our understanding of cyclical neutropenia 3. Though the gene modified responsible for canine cyclical neutropenia has been identified, the dynamic origin of the cycling is only partially understood. Because of its interesting dynamical nature, many mathematical models have been formulated to attempt to answer this question. While many have modeled cyclical neutropeni a as arising only from destabilization of neutrophil dynamics , the work of FFT and suggest that the origin of cyclical neutropenia lies in a destabilization of the combined HSC and neutrophil control system. The hypothesis that oscillations originate in the stem cells is supported by the observation that in cyclical neutropenia oscillations are also present in platelets and reticulocytes. Cyclical neutropenia in humans is often treated using granulocyte colony stimulating factor (G - CSF) , which is known to interfere with apoptosis 4 . In the authors showed that, depending on the starting date of the G - CSF treatment, the neutrophil count could either be stabilized or show large amplitude oscillations . Their model suggested that other G - CSF treatment schemes c ould be effective while using less G - CSF. However, this model included neither erythrocyte nor platelet dynamics even though clinical data indicates oscillations in those cell lines in cyclical neutropenia patients. Thus it is not known if the results woul d be consistent with observed platelet and reticulocyte data. Second, the simulations did not take into account the pharmacokinetics of G - CSF. In this chapter, we present a new model for assessing the effects of G - CSF treatment in cyclical neutropenia. To do this, we augment the comprehensive model of the hematopoietic system from by Research Journal of Recent Sciences ______ _ _ _______________________________ ______________ _ ________ ISSN 2277 - 2502 Vol. 1( 4 ), 1 4 - 21 , April (201 2 ) Res. J. Recent Sci. International Science Congress Association 15 coupling it with a two - compartment pharmacokinetic model that accounts for G - CSF kinetics 5. Material and Methods Periodic Auto - Immune Hemolytic Anemia: Le  (, α, β) e the oulion of recuror cell  ie , e α n re of roucion β. le V(R) and U(R) be the velocities of rate of production and maturation, which may depend on the hormone (EPO) concentration. If N (R) is the number of cells recruited into the rep roduce precursor population, and then the entry of new precursor cells into the structured model will satisfy the boundary condition, and boundary condition rate cells exchange Le he irh re for reroucin recuror cell e β n α represent the death rate through apoptosis. Let be the density of the distribution of maturity levels of the cells when released into the circulation blood, where represents the mean age of mature precursor cells and The disappearance rate function is given by: With these condition the age structured model for the population of precursor cells with �t 0, 0 < μ < μ 1 n 0 < μ < μ 2 satisfies: le  (, α, 0) e he oulion of ure non – reproducing cell  r n e α. Fro he iernce re funcion, he boundary condition for cells entering the mature population is given by where the maturity level μ 1 n μ 2 represents the maximum age for a cell reaching maturity. We assumed that destruction of RBC occurs by active removal of the old cells 6. Form a modeling point of view, this results in a moving boundary condition with the age of the ole RBC, ϴ() vryin in . The ounry coniion i hen given by where F i he fixe RBC reovl re. If γ(ϴ) i he eh re of mature cells , then the partial differential equation describing  (,ϴ)i iven y: The total population of mature cells function is given by: The differential equation for R (red blood cells production) is thus: Where k is the decay constant for the hormone and the rate of roucion β i iven y  onoone ecrein Hill funcion. The given equation is linear in R. IF (integrating factor) = = Cyclical Thrombocytopenia: Platelets are blood cells whose function is to take part in the clotting process, and thrombocytopenia denotes a reduced platelet (thrombocyte) count. In cyclical thrombocytopenia (CT), platelet counts oscillate generally from very low values (1 109 cells/L) to normal (150 − 450 109 platelets/L) or above normal levels (2000 109 cells/L). These oscillations have been observed with periods varying between 20 and 40 days . In addition, patients may exhibit a variety of clinica l symptoms indicative of impaired coagulation such as purpura, petechiae, epistaxis, gingival bleeding, menorrhagia, easy bruising, possibly premenstrually, and gastrointestinal bleeding . There are two proposed origins of cyclical thrombocytopenia. One is of auto - immune origin and most prevalent in females. The other is of amegakaryocytic origin, more common in males. Autoimmune cyclical thrombocytopenia is characterized by a shortened platelet lifespan at the time of decreasing platelet counts . This is co nsistent with normal to high levels of bone marrow megakaryocytes and with an increased destruction rate of circulating platelets . Autoimmune CT has also been postulated to be a rare form of idiopathic (immune) thrombocytopenic purpura (ITP) . Cyclical Neu tropenia: We presented a two variable delay differential equation (DDE) system that has negative feedback loops in both the peripheral loop and the stem cell loop and illustrates the four compartments of the model: the hematopoietic stem cell (HSC) compar tment (S), the neutrophil compartment (N), the erythrocyte compartment (R) and the platelets compartment (P). The HSCs are assumed to be self - renewing, and thus cells in the resting (S 0 ) phase can either enter the proliferative phase at rate K (S) or diffe rentiate into neutrophils (N) at rate F(N). As the neutrophil precursors differentiate, their numbers are amplified by a factor A, which accounts for both successive divisions and cell loss due to Research Journal of Recent Sciences ______ _ _ _______________________________ ______________ _ ________ ISSN 2277 - 2502 Vol. 1( 4 ), 1 4 - 21 , April (201 2 ) Res. J. Recent Sci. International Science Congress Association 16 apoptosis. It is also assumed that apoptosis occurs during the roliferive he  re γ n h ure neurohil ie  re α. A cn e een in f ig ure 1, the system is controlled by four negative feedback loops. The first one regulates the rate K (S) of reentry of HSCs to the proliferative cycle, and it operates wih  ely τ h ccoun for he ie reuire o rouce wo daughter cells from one mother cell. The equations for the two variables N, R, P and S can be derived from a time age maturation formulation, or written directly from consulting f ig ure 1. For the compartment N, the loss is the efflux to death N and the production of mature neutrophils is equ al to the influx F (N) S from the HSC compartment times the amplification A. Since one needs to take ino ccoun he rni ie τ, he ol roucion of ure neutrophils are AF(N (t - τ ))S( - τ), AF(R ( - τ ))S( - τ) n AF(P (t - τ ))S( - τ ) or equivalently AF(N )S, AF(R)S and AF(P)S (recall that N = N (t - τ ), R  R ( - τ ) n P  P ( - τ )). For the second variable, the loss from the compartment S is the flux reentering the proliferative phase, K (S)S, plus the efflux going into diff erentiation, F(N )S. The production of S is equal to the flux of cells reentering and surviving the proliferative phase, given by K (S)S , times the cell division factor 4. The dynamics of S is then described by The feed back functions F and K are monotone decreasing Hill functions. F controls the number of compartments while K(S) regulates the level of HSCs. It developed a comprehensive model that contains four compartments: the HSC (Q), the neutrophils (N), the erythrocytes (R) and the platelets (P). This model combines a number of compartmental models. The stem cell and neutrophil dynamics are based on the model and the erythrocyte and platelet compartment are simplified models. The circulat ing cells are coupled to each other via their common origin in stem cell compartment. Regulatory negative feedback loops determine how much differentiation from the stem cells each cell line will undergo. Since it takes several days to produce a mature cel l from a newly differentiated cell, time delays appear in the equations. The model consists of a set of four coupled delay differential equations. Periodic Chronic Myelogenous Leukemia: Leukemia is a cancer of the blood or bone marrow characterized by an abnormal proliferation of blood cells, usually leucocytes. Chronic myelogenous leukemia (CML) is distinguished from other leukemias by the presence of a genetic abnormality in blood cells , called the Philadelphia chromosome, which is a translocation between chromosomes 9 and 22 that leads to the formation of the BcrAbl fusion protein. This protein is thought to be responsible for the dysfunctional regulation of myelocyte growth and other f eatures of CML . DDE Models: Delay - differential equations (DDEs) are a large and important class of dynamical systems. They often arise in biological systems where time lags naturally occur . In particular, in hematology several processes are controlled th rough feedback loops and these feedbacks are generally operative only after a certain time, thus introducing a delay in the system feedback. The general form of a DDE for x(t) R n is where x is the delayed variable (x(t − )) and f is a functional operator in R R n C 1 . DDE with Constant Delays: Delay differential equations with constant delays take the form where the quantities i , i = 1, 2, ..n are positive constants. For simplicity, consid er the DDE with a single constant delay: To obtain a solution of e quation (26) for t � 0, one needs to specify a history function on [ − , 0]. Indeed, recall that for an ordinary differential equation (ODE) system with n variable s, one would only need to specify the initial values x(0) for each of the n state variables. In order to solve a DDE, one needs to specify not only the value at t = 0, but also all the past values of x(t) over the interval [ − , 0]. For example, let X( t) represent the circulating cell population of a certain type of blood cell, assume that is the random rate of loss of cells in the circulation and F is the flux of cells from the previous Research Journal of Recent Sciences ______ _ _ _______________________________ ______________ _ ________ ISSN 2277 - 2502 Vol. 1( 4 ), 1 4 - 21 , April (201 2 ) Res. J. Recent Sci. International Science Congress Association 17 compartment. Then, the dynamics of the number of circulating c ells will have the generic form DDE with Distributed Delays: A distribution of delays is then be more appropriate and the DDE becomes an integro - differential equation of the form The density G(u) of the distribution function is referred to as he memory function or the kernel and is normalized, i.e. This type of model can also be interpreted as allowing for a stochastic element in the duration of the delay . Also, we will see t hat for some densities G(u), e quation (27) can be equivalently viewed as a system of ordinary differential equations 7. ODE Models: Delay differential equations naturally arise in modeling biological systems. Consider the following DDE system with a distr ibuted delay: with the special choice of the density of the gamma distribution for the memory function where a is a positive number and p is a positive integer or zero. Note that the function G( u) has a maximum at u = p/a and that, as a and p increase, keeping p/a fixed, the kernel approaches a delta function and the distributed delay approaches the discrete time delay with = p/a. Moreover, it is clear that the following three properties are s atisfied: The central idea of the method is to replace the distributed delay by an extension of the set of variables. Define p + 1 new variable as Then, using the properties of G one can show that these new variables satisfy a sequence of linear ODEs Solving the following system is thus equivalent to solving the DDE problem (28), given that the new variables are given appropriate initial values: The linear chain trick could be useful for numerical computations since it reduces the problem to an ODE system, for which several numerical methods are available. However, this method cannot be used for all sorts of delays 8 . Age - Structured Models: Let x(t, a) be the cell density at ti me t and age a in a generic compartment. We assume that cells disappear (die) at a rate (t). We also assume that the cells in the compartment age with a velocity V (t) and that a cell enters a compartment at age a = 0 and exits this compartment at age a = . Therefore, the equation satisfied by x(t, a) is an time - age equation (advection, or reaction - convection, equation): The right hand side in this equation represents the rate at which cells in the age interval a to a + a disa ppear at time t. To represent the manner in which new cells enter the compartment, we define the boundary condition (B.C.) x(t, 0) = H(t). Finally, to fully represent the problem, we specify the initial condition (I.C.) x(0, a) = (a). Using the method of characteristics , we obtain the following delay differential equation: where X(t) is the total number of cells and T satisfies . Note that if is a constant, e quation (29) reduces to In addition, if the aging velocity is constant (V (t) = V ), we have that T satisfies =VT which implies that T = /V . Hence, if and V are constant, we obtain a delay differential equation with constant delay: Model Composition: The model has several negative feedback functions that regulate stem cell proliferation and differentiation into the three circ ulating cell types. In this model, as in the one presented, these are given by: Research Journal of Recent Sciences ______ _ _ _______________________________ ______________ _ ________ ISSN 2277 - 2502 Vol. 1( 4 ), 1 4 - 21 , April (201 2 ) Res. J. Recent Sci. International Science Congress Association 18 We must also specify an input function I(t) that represents the subcutaneous G - CSF injections. We assume that this input is brief in duration, and that the total amount of G - CSF added corresponds to the desired dosage, namely: Note that if is small, a Gaussian - like input approximates a Di rac - function, and we can write Therefore to simulate periodic injections, we let where H (t) denotes the Heaviside step function The day on which treatment is initiated is denoted by d, and the Heviie funcion ily urn he injecion on. The er ‘‘ o T’’ enure erioiciy, n we reuire h T so that the approximation to the integral remains valid. Finally, we en sure that holds by choosing the parameter a a = dosage 9. It remains only to describe how the G - CSF acts on the hematological portion of the model. Because we believe from previous modeling efforts that A N , S , and 1 are the parameters that need to change under G - CSF, we model G - CSF injections as causing fluctuations in these three parameters: The uercri ‘‘’’ n ‘‘u’’ reecively inice vlue corresponding to values that, in the model without the dynamics of G - CSF, match treated and untreated data respectively. The parameters m A , m g , and m t are slopes that specify how much A N , S , a nd 1 change in response to a given change in G - CSF concentration, G. is the average G - CSF concentration for each data set. These were computed using the G - CSF model alone, and using the average neutrophil levels in each data set. Results and Discussi on In each case, we found that the neutrophil amplification increases substantially under G - CSF treatment, as does the rate of stem cell apoptosis, and the differentiation into the neutrophil line 10. We therefore predict similar changes for the remaining dogs there is some redundancy in the model, in that increasing the neutrophil amplification and the differentiation into the neutrophil line from the stem cells has similar effects. This is not unexpec ted, since the primary effect of both changes is to raise neutrophil levels. Figure 2. shows the fit of the untreated and treated data for d ogs 100, 118 and 127. Note that the analysis matches the data reasonably well for the neutrophils as well as for the erythrocytes and platelets. This confirms that the new model, with the G - CSF coupled to the cell population dynamics, is capable of reproducing the data. The least squares differences between the simulations and the data were not significantly less than t he values reported. These simulations and data are for daily treatment. Figure 2 shows the data and analysis for the other four dogs ( d ogs 101, 113, 117 and 128), again with daily treatment. Recall that these were the estimated, not fitted, values for the treated parameters and note the quality of the fits. Thus, we are able to match observed data without automated parameter fitting based simply on an examination of the treated data and the parameter changes for d ogs 100, 118 and 127. For each dog, we perfo rmed simulations comparing daily treatment, treatment every other day, and every three days. We find that particularly for d ogs 100, 101, 118 and 127, changing the period of the treatment can significantly affect the nature of the oscillations. It shows th e results of treating d og 118 every other day, rather than every day. We have also explored the effects of changing the time at which the treatment is initiated. In most cases, this did not significantly change the long - term behavior. However, for d og 127 the amplitude of the oscillations was significantly reduced when the treatment was initiated in the latter half of the cycle. More specifically, measured from day 1, we find that smaller oscillations occur if treatment is initiated on day 8 or afterwards, or on days 2 or 5 (see f igure 2). When treatment was initiated on other days, larger oscillations in the model resulted. We were aware from our previous study of similar models that there is the possibility that two or more qualitatively different states c an be locally stable, and we have also found evidence for this in the present model. Namely, changing the treatment onset time from day 1 to day 8 for d og 127 caused the simulation to stabilize to two very different types of behavior. It should also be not ed that increasing the G - CSF dosage in the model sometimes helped to stabilize oscillations ( d og 127), but in several cases ( d ogs 100, 128 n 101)  oe incree fro 5 μ/k o  oe in he range 15 - 25 μ/k cue oe iulion o fil. In h analysis, the differentiation rate out of the stem cells was so high, and the apoptosis rate in the stem cells was so high, that the stem cell population was no longer able to maintain itself. For the other dogs, there was always a dosage that was suffici ently high to terminate the FFT analyze, but it was sometimes a factor of 10 higher than the actual dosage given. Research Journal of Recent Sciences ______ _ _ _______________________________ ______________ _ ________ ISSN 2277 - 2502 Vol. 1( 4 ), 1 4 - 21 , April (201 2 ) Res. J. Recent Sci. International Science Congress Association 19 Figure 2 Continuous neutrophil data and analyze for d ogs 100, 118, 127, 101, 113, 117 and 128. The left figure shows untreated data (red) and fit (normal area). The right figure shows treated data (green) and analyze for dogs under daily G - CSF treatment. Note that the model accounts for the different scaling in neutrophil counts. The calculations were obtained using parameters resulting from th e FFT analysis method. Neutrophil units are 108 cells - kg − 1 Research Journal of Recent Sciences ______ _ _ _______________________________ ______________ _ ________ ISSN 2277 - 2502 Vol. 1( 4 ), 1 4 - 21 , April (201 2 ) Res. J. Recent Sci. International Science Congress Association 20 We have developed a model of the hematopoietic system that includes a pharmacokinetic model of G - CSF dynamics in tissue and in circulation. The model is able to account for the features of untreated, and G - CSF - treated, data for dogs with cyclical neutropen ia. This is accomplished, starting with parameter fitting done , by fitting parameters for 3 dogs and thereby estimating, not fitting, parameters for 4 other dogs. One of the most intriguing observations resulting from the parameter fitting in this study is that to fit observed data for cyclical neutropenic dogs and human patients during G - CSF treatment it was necessary to assume that there was an increase in the rate of apoptosis in the stem cell compartment during G - CSF treatment, at the same time as the m ore expected increase in neutrophil amplification. Research Journal of Recent Sciences ______ _ _ _______________________________ ______________ _ ________ ISSN 2277 - 2502 Vol. 1( 4 ), 1 4 - 21 , April (201 2 ) Res. J. Recent Sci. International Science Congress Association 21 The study we report here about treatment schedules indicates that changing the period of the treatment from daily to every other day, and then to every third day, almost always significantly alters the nat ure of the oscillations. Since G - CSF is costly and may have undesirable side effects, it may be worth exploring this option further in humans. Furthermore, we found in one case (Dog 127) that changing the time of onset of treatment to the latter half of th e cycle (as measured by setting day 1 to be the day when the neutrophil level is minimal) results in much smaller amplitude oscillations in the treated FFT analysis. In the model, both of these interventions had more significant effects on the oscillations than did changing the G - CSF dosage. Indeed, increasing the dosage was not seen to be a viable option in our analysis, as it frequently led to the termination of the FFT rather than to the stabilization of oscillations 11 . The observed data are highly vari able from one dog to another, but the simulations can be individualized to account for this. Thi reen he oiiliy of uin rel ie”  for  iven dog to individualize model analysis and make predictions about the effects of different treatme nt schedules. Earlier modeling work also suggested that significantly different behavior would result from different G - CSF treatment schedules. Our model substantiates this, and quantifies the effects using realistic G - CSF dynamics and yielding analysis th at are directly comparable to observed data. Our central result is that in the model, changing the time of treatment initiation and or the period of treatment may result in equally good, or better, long - term outcomes and may require less G - CSF. These chang es would be practical to implement and, if less G - CSF were required, would reduce the risk of side effects as well as the cost of treatment 12 . Conclusion The model could be used to study different mechanisms of granulopoiesis and of G - CSF administration. In particular, we assumed G - CSF was acting on both the apoptosis rate and the amplification factor in the proliferative phase of neutrophil precursors. Since both effects are modeled separately, we could study in more detail the precise effects of each fac tor. As mentioned above, clinical data for alternative treatment schedules with G - CSF are needed to validate our results and make further model improvements. If such data were available, the individualized approach presented in the paper could be interesti ng to implement since it fits the model to data before and during treatment for a given subject. FFT simulations for this subject could then predict the possible outcomes of different treatment schemes. In conclusion, hematopoiesis and G - CSF effects are no t yet fully understood and several aspects are still being studied. Clinical findings and future studies will provide new insights and help to better understand the system. The models presented here would then have to be modified accordingly, as this is pa rt of the evolutionary process of mathematical modeling. References 1. Adimy M. and Crauste F. , Global stability of a partial differential equation with distributed delay due to cellular replication , Nonlin. Analysis., 54 , 1469 – 1491 (2003) 2. Balamuralitharan S. and Rajasekaran S., A Mathematical Age - structured Model on Aiha Using Delay Partial Differential Equations. Australian Journal of Basic and Applied Sciences, 5(4) , 9 - 15 (2011) 3. Balamuralitharan S. and Rajasekaran S., A Mathematical Age - structured Model on CN Using Delay Partial Differential Equations , Proceedings of the International Conference on Emerging Trends in Mathematics and Computer Applications , MEPCO Schlenk Engineering College, Sivakasi, India , 183 - 187 (2010) 4. Belair J., Mackey M.C., and Mahaffy J.M. , Age - structured and two - delay models for erythropoiesis , Math. Biosci., 128 , 317 – 346 (1995) 5. Bernard S., Belair J. and Mackey M. , Oscillations in cyclical neutropenia: New evidence based on mathematical modeling , J. Theor. Biol., 223 , 283 – 298 (2003) 6. Beuter A., Glass L., Mackey M.C. and Titcombe M.S. Nonlinear Dynamics in Physiology and Medicine , Springer - Verlag (2003) 7. Colijn C., Fowler A.C. and Mackey M.C. , High frequency spikes in long period blood ce ll oscillations , J. Math. Biol., 53 , 499 – 519 (2006) 8. Dale D.C. and Hammond W.P. , Cyclic neutropenia: A clinical review , Blood Rev., 2 , 178 – 185 (1988), 9. Foley C. and Mackey M. , Dynamic hematological disease , A review , J. Math. Biol., In press (2008) 10. Foley C. , Bernard S. and Mackey M. , Cost - effective G - CSF therapy strategies for cyclical neutropenia: Mathematical modelling based hypotheses , J. Theor. Biol., 238 , 754 – 763 (2006) 11. Haurie C., Dale D.C., Rudnicki R., and Mackey M.C. Modeling complex neutrophil dynamics in the grey collie , J. Theor. Biol., 204 , 504 – 519 (2000) 12. Mahaffy J., Belair J. and Mackey M. , Hematopoietic model with moving boundary condition and state dependent delay: Applications in erythropoiesis , J. Theor. Bio., 190 , 135 – 146 (1998)