Abstract
Background
Hematopoiesis is a complex process involving different cell types and feedback mechanisms mediated by cytokines. This complexity stimulated various models with different scopes and applications. A combination of complementary models promises to provide their mutual confirmation and to explain a broader range of scenarios. Here we propose a combination of an ordinary differential equation (ODE) model of human granulopoiesis and an agentbased model (ABM) of hematopoietic stem cell (HSC) organization. The first describes the dynamics of bone marrow cell stages and circulating cells under various perturbations such as GCSF treatment or chemotherapy. In contrast to the ODE model describing cell numbers, our ABM focuses on the organization of individual cells in the stem population.
Results
We combined the two models by replacing the HSC compartment of the ODE model by a difference equation formulation of the ABM. In this hybrid model, regulatory mechanisms and parameters of the original models were kept unchanged except for a few specific improvements: (i) Effect of chemotherapy was restricted to proliferating HSC and (ii) HSC regulation in the ODE model was replaced by the intrinsic regulation of the ABM. Model simulations of bleeding, chronic irradiation and stem cell transplantation revealed that the dynamics of hybrid and ODE model differ markedly in scenarios with stem cell damage. Despite these differences in response to stem cell damage, both models explain clinical data of leukocyte dynamics under four chemotherapy regimens.
Conclusions
ABM and ODE model proved to be compatible and were combined without altering the structure of both models. The new hybrid model introduces model improvements by considering the proliferative state of stem cells and enabling a cell cycledependent effect of chemotherapy. We demonstrated that it is able to explain and predict granulopoietic dynamics for a large variety of scenarios such as irradiation, bone marrow transplantation, chemotherapy and growth factor applications. Therefore, it promises to serve as a valuable tool for studies in a broader range of clinical applications, in particular where stem cell activation and proliferation are involved.
Background
Hematopoietic stem cells (HSC) have been in the focus of research since the beginning of last century [1]. Easy accessibility and handling, in combination with elegant experimental techniques like clonal assays [2,3] made the hematopoietic system the best studied mammalian stem cell system. As a consequence, the first models were designed in the 1960s [4,5]. The process of hematological homeostasis is characterized by a relative stability of the (small) stem cell pool and a massive amplification along the differentiation process, leading to a daily production of about 10^{11}10^{12} mature blood cells [6]. This observation led to the so called pedigree concept, which postulates that stem cells originate only from stem cells, i.e. either maintain the stem cell state or lose it irreversibly [7].
This concept represents a core assumption of most mathematical models for hematopoiesis that have been formulated complementary to experiments. Some do not explicitly model stem cells but include them as a source of cellular influx into the modeled differentiation stages of hematopoiesis [810]. Models that explicitly model the hematopoietic stem cell population mostly focus on the cell number of one [11,12], or more populations (such as a resting and proliferating cells [13]. Considering cell numbers these models implicitly ignore intercellular homogeneity. A few models do consider structured cell populations and introduce an additional cellular feature [14]. However, all these models share the concept of unidirectional cell flux towards differentiated states.
Following this concept, we also developed ordinary differential equations (ODE) based lineage models of human granulopoiesis, erythropoiesis and thrombopoiesis [1519]. All these models are supplied by the same stem cell model. They describe the dynamic regulation of HSC, proliferating and maturing progenitors, mature blood cells and cytokines of the hematopoietic system and aim at predicting the complex dynamics of hematopoiesis during combined chemotherapy and growth factor applications. A number of feedback loops control differentiation and amplification, e.g. via the probability of stem cell selfrenewal, amplification rates and maturation times of committed cells. Models of pharmacokinetics and dynamics of growth factor and chemotherapy applications were introduced recently, allowing precise predictions of clinical data in numerous scenarios [19].
The strictly hierarchical pedigree concept was challenged as experimental evidence for stem cell flexibility was found at the end of the last century. Cells from neural [20], skeletal [21,22] and vascular tissue [23] were shown to be capable of engraftment in irradiated hosts and to contribute subsequently to the production of mature blood cells. Most likely this flexibility is induced and controlled by the stem cell environment [24,25]. Our formerly developed agentbased model (ABM) of hematopoietic stem cell organization incorporates such a context dependent stem cell regulation by considering two stem cell growth environments (GE) [26]. In one of these two GE, which can be interpreted as a stem cell niche[24], the stem cells are quiescent and regenerate their ability to remain in the niche, while in the other GE they proliferate and lose this niche affinity. The broad range of successful applications of this modeling concept includes modeling of clonal competition [27,28], agedependent repopulation potential after irradiation [29] and treatment of chronic myeloid leukemia with Imatinib and interferonα [3032].
Certainly, the aim of systems biology and modeling is generation of consistent explanations for as many scenarios as possible. Therefore, the following natural question arises from the diversity of individual models and the different scenarios they explain: Can we assume their compatibility and even expect synergistic advantages from their combination? Desirably, the combined models should cover (many) complementary scenarios for maximal synergistic effects, but apply for the same scenarios for compatibility. Due to this premise, we pursue this question with an attempt of secondlevel modeling: a combination of our ODE model of human granulopoiesis that follows the hierarchical paradigm and our agentbased stem cell model. Both models were developed on the basis of different biological evidence, experiments and data and are well established. Thus our effort was to combine them without changing their assumptions, equations, parameters and application scenarios as far as possible.
If the models prove to be compatible, we can expect the mentioned synergistic effect for a set of reasons. The ODE model turned out to be highly sensitive to changes in its stem cell compartment (SCC) [17]. Exchanging its homogeneous stem cell population for a structured stem cell model enables modeling of phenomena caused by heterogeneity within the SCC, such as lineage priming of stem cells, clonal competition like leukemia development and treatment or cellcycle specific susceptibility to chemotherapy or radiation. Although ODE model simulations of granulocyte dynamics under chemotherapy agree well with clinical data, predicted reduction of stem cell numbers appears to be overestimated [33]. On the other hand, the ABM does neither consider effects of growth factors nor their pharmacokinetics. As a result, the ABM is not able to correctly account for shortterm effects in the reconstitution of blood cell number after chemotherapy which is often supported by growthfactor applications. Besides, ABM simulations for realistic cell numbers from stem cells to differentiated cells are computationally expensive.
Therefore, we here propose an integration of these model concepts guided by the motivation that systemsbiologic modeling is an iterative process. We substitute the ODE stem cell model by a difference equationbased formulation of the ABM describing continuous cell numbers [34]. By this approach, we expect to construct a more comprehensive model with a broader range of possible applications, e.g. in the context of planning of clinical trials.
The adaptations necessary to combine the models will be discussed in detail. We compare the new hybrid model with the former ODE model in a number of qualitative tests, e.g. with respect to chronic irradiation or bleeding. Finally, we used the hybrid model to simulate four different chemotherapy regimens comprising both, dose and time intensifications of therapy and growthfactor applications as well.
Methods
In the following, we briefly explain our ODE model for human granulopoiesis, our ABM stem cell model and its difference formulation. Building on this background we then establish our hybrid model.
ODE model for human granulopoiesis
We here summarize the central features of the ODE model of human granulopoiesis. The complete set of equations and parameters is given in the Appendix (A1.1 and A1.2). All model assumptions are comprehensively discussed in Scholz et al. [17]. The model describes the dynamics of concatenated cell compartments by a set of ODEs (Figure 1). Cell compartments represent cell numbers of morphologically distinguishable cell stages of granulopoiesis, namely HSC, colonyforming units, maturing blasts in the bone marrow and mature granulocytes in peripheral blood (see Table 1) [17]. The two growth factors granulocytemacrophage colony stimulating factor (GMCSF) and granulocyte colony stimulating factor (GCSF) are also explicitly modeled and regulate amplification and transit times in bone marrow compartments.
Figure 1. Schematic representation of the ODE model. Compartments representing cell types and growth factors are depicted as boxes (s. Table 1). Compartment MGB is substructured into G46 to model maturation with Γdistributed transit times. Cell fluxes between compartments are shown as black arrows, the effect of chemotherapy as gray arrows with associated kill rate k_{X} indicated by a label and feedbacks as colored arrows: intrinsic stem cell feedback (red), feedback from bonemarrow cells to stem cells (green) and feedback between later stages of granulopoiesis mediated by explicitly modeled growthfactors (blue). The colored arrows indicate the input for the feedback functions (s. Appendix A2.1) that dynamically control compartment parameters mentioned by the labels.
Table 1. List of all compartments in the ODE model and their biological equivalents
In all individual cell compartments, except for stem cells, a general balance equation reflects the dynamics of growth, differentiation and cell death due to chemotherapy: change in compartment size = influx × amplification  efflux  cell loss.
Influx is given by the efflux of the upstream compartment, while amplification and efflux of compartment X are computed from the amplificationrelated parameters proliferative fraction a_{X} and amplification A_{X} and the residencerelated parameters transit times T_{X} (or in compartment S probability of selfrenewal p). The stem cell compartment represents the root of the hematopoietic tree without influx. In this case, the first term of the balance equation is simply replaced by proliferation/amplification.
Amplification and differentiation is controlled by the four feedback loops sketched in Figure 1. Amplification and selfrenewal in the stem cell compartment are regulated by stem cell number C_{S} (Figure 1, red) and the number of granulopoietic bone marrow cells C_{G} (Figure 1, green). A decrease/increase in both cell numbers C_{S} and C_{G} results in an increase/decrease, respectively, in proliferation, while the effects of C_{S} and C_{G} on selfrenewal are contrary: low stem cell numbers C_{S} result in high selfrenewal with higher priority than low bone marrow cell number C_{G} causing increased differentiation. These two feedback functions are adopted from Wichmann and Loeffler (s. Appendix A1.1) [15].
The remaining two feedback loops (blue) control maturation and amplification in the whole bone marrow except for the SCC. They are mediated by the growth factors GCSF and GMCSF (s. Appendix A1.1) [3537]. GCSF acts in different modes: It increases proliferation, maturation and the release of bone marrow cells to the blood. GMCSF regulates amplification in CG. Endogenous production of both growth factors depends on cell numbers in the bone marrow and the blood. Subcutaneous and intravenous GCSF applications were included by a simple pharmacokinetic model. We model a delayed influx from the subcutaneous compartment by a twocompartment system. Degradation in the central compartment is assumed to happen unspecifically and specifically with the specific degradation being proportional to the number of circulating granulocytes. We assumed that applied GCSF acts in the same way as endogenous GCSF, i.e. novel GCSF pharmaceuticals such as pegylated GCSF with different pharmacodynamics were not considered here.
According to mouse data, chemotherapy was modeled as an acute transient depletion of bone marrow cell stages following a first order kinetics depending on the drugs applied and the cell compartments affected. (s. Appendix A1.1) [33].
All cell compartments are normalized with respect to the setting C_{S} = 1 in steadystate. Hence, all subsequent compartments are given in units of equilibrium stem cell number. Furthermore, we usually present relative compartment sizes Cx^{rel} = Cx / Cx^{nor}. The steadystate values Cx^{nor} for all compartments X are given in the Appendix (A1.3).
The agentbased stem cell model
In contrast to the representation of the stem cell population in the ODE model, our agentbased stem cell model represents each cell individually. A general scheme of the model is given in Figure 2. Each cell is characterized by the growth environment GE that accommodates the cell and a niche affinity a∈[a_{min}, a_{max} ] [26]. In each simulation step the cells can either remain in or switch between the two GEs. GE Α and Ω represent functionally distinct signaling contexts and may be interpreted as a niche and a nonniche environment, respectively. In the niche environment Α, cells are quiescent and regenerate niche affinity a, while in the nonniche GE Ω they proliferate and lose niche affinity. If a cell acquires an affinity lower than a_{min}, it loses its ability to switch back to Α, to regenerate affinity a and to selfrenew. Then it leaves the stem cell compartment and enters further differentiation and amplification stages. When a cell switches GEs, its affinity a remains constant. The evolution of each cell’s affinity a at time t is given by:
Figure 2. Scheme of the agentbased stem cell model and its coupling in the hybrid model. Each cell is characterized by its affiliation to one GE, niche affinity a < a_{max} and cell cycle position c. The two GEs represent functionally different environments: in GE Α affinity a increases with time and is limited by a_{max}. In GE Ω it decreases and cells proliferate. When a drops below the threshold a = a_{min}, the cells lose the ability to switch to Α and leave the stem cell compartment. In the hybrid model they enter the progenitor compartment CG. The effect of chemotherapy CX was modeled as cell loss in GE Ω without affecting GE A.
Here r > 1 and d > 1 are the regeneration and differentiation coefficients, respectively. The cells at a_{max} in GE Α conserve their state unless they switch GEs.
The transitions between the GEs happen randomly with transition intensities (probability per time step) that depend on cell numbers in the target GEs and the individual cell’s affinity. These dependences introduce an intrinsic regulation of stem cell proliferation and differentiation. For a cell of affinity a and a HSC population of N_{Α} cells in GE Α and N_{Ω} cells in GE Ω, they are given by:
where α and ω are the intensities of transitions Ω → Α and Α → Ω, respectively. The functions f_{α} and f_{ω} model a limited capacity of both GEs approaching zero for increasing cell numbers in the target environment. Their definition is given in the Appendix (A2.1).
Because at a_{max} the cells are quiescent and the probability of their transition to Ω, and with it activation, is minimal, these cells can be interpreted as dormant stem cells [38].
Proliferating cells in Ω are additionally characterized by their cell cycle position c. In each simulation step it is simply increased by the time step: c_{t+1} = c_{t} + Δt. The established parameter set for human HSC assumes 48 h for a complete cell cycle, which is divided into two subphases: an intermediate G_{1}phase of 32 h (0 h ≤ c < 32 h) and a mitotic phase of 16 h (32 h ≤ c ≤ 48 h), subsuming S/G_{2}/Mphases [32]. During mitotic phase, no transitions to GE Α are allowed. At the end of cell cycle, a cell divides into two daughter cells. They inherit the affinity of their mother cell and start cell cycle in G_{1}Phase (c = 0 h). If a cell switches to GE Α, it always becomes quiescent. If a cell switches to GE Ω, it enters the mitotic phase at c = 32.
A summary of all parameters of the stem cell model is given in the Appendix (A2.2).
Difference equation formulation
When considering realistic cell numbers, simulations of individual cells are computationally expensive. To replace the ODE SCC we thus decided to use a difference formulation of the ABM introduced in the last section which was established by Kim et al. [34]. This difference formulation is based on exactly the same equations and parameters as the ABM. It reduces computational effort by defining a discrete set of affinities and describing numbers of cells of those affinities instead of individual cells of arbitrary affinities [34]. Considering Eq. (1) this discrete set is constructed as follows:
1. If a cell in GE Ω starts with maximal affinity a_{max} = a_{0} = 1 and does not change GE, it will adopt the affinities a_{k} = d ^{−k} ; k∈ℕ_{0} and eventually leave GE Ω at a_{min}. Hence, defining , its niche affinities in GE Ω are fully described by the set with the associated affinities a_{k} = d ^{−k} . Analogously, all niche affinities in A can be mapped onto with and the affinities a_{k} = r^{−k}.
2. If a cell switches GEs, its affinity is conserved (Eq.1) and to set up the difference equations formulation two matching sequences of affinities are required in the two GEs. This is not the case in general, if d ≠ r. But if ln r / ln d is a rational number with ln r / ln d = v/w, the set of all possible affinities is given by with ρ = ln r / v = ln d / w. The established parameters of the human system are r = 1.1 and d = 1.05 with ln r / ln d ≈ 2. Setting ρ = 0.0488 results in the slightly modified coefficients r = e^{2ρ} = 1.10252 and d = e^{ρ} = 1.05001. Together with a_{min} = 0.002 from the human parameter set, this results in k_{max} = 127.
Hence, the difference formulation describes the number of cells Α_{k} and Ω_{k} in compartments that are associated with the affinities
While the quiescent cells in GE Α are fully characterized by their affinity a, in GE Ω the cellular state is given by the combination of affinity a and cell cycle position c. The number Ω_{k} of cells of affinity a_{k} in Ω is subdivided into the number of cells Ω_{k,c} of cell cycle position c. With c running from 0h to 48 h and a simulation time step of Δt = 1 h this leads to 49 steps, a cell cycle duration of τ = 49 h and a 128 × 49 matrix Ω_{k,c} of cellular states. Thus, the cell numbers N_{Α} and N_{Ω} that determine the transition characteristics are given by:
where A_{k} is the number of cells of affinity a_{k} in GE Α and Ω_{k,c} is the number of cells of affinity a_{k} and cell cycle position c in GE Ω. For each value of a_{k}, the transition probabilities α and ω are given by Eq. (2). For computation time, we approximated the binomially distributed number of cells that switch GEs with their expectation values ω(a_{k},N_{Ω}(t)) · A_{k}(t) and α(a_{k},N_{Α}(t)) · Ω_{k,c}(t) for c = 0…31. These settings finally result in the following set of difference equations:
Construction and implementation of the hybrid model
In the hybrid model, the difference formulation of the ABM introduced in the last section replaces the ODE SCC. Accordingly, stem cell number C_{S} is now given by N_{S} = N_{Α} + N_{Ω} (see Eqs. (4) and (5)). Naturally, the new SCC must comply with the coupling of SCC and other compartments in the ODE model regarding cell fluxes, chemotherapy effects and feedbacks (Figure 1).
Cells passing the threshold a_{min} lose the ability to return to GE Α and to selfrenew which means that they enter subsequent maturation stages. They constitute the efflux of the stem cell compartment and the influx into the committed progenitor compartment CG. The equilibrium efflux of the ABM stem cell compartment relative to equilibrium stem cell number is and, thus, considerably smaller than in the ODE model (). However, the stem cell uncertainty principle [39] makes is hard to identify stem cells and, thus, their total number precisely. On the other hand, the ODE model normalizes all cell numbers with regard to stem cell number, while the ABM uses the parameters Ñ_{Α} and Ñ_{Ω} for scaling stem cell number which have no influence on qualitative model behavior. Thus, we consider the discrepancy in output as a scaling problem and decided to match equilibrium efflux of the ABM with the one in the ODE model by using the scaling parameter k_{scaling}. Therefore, in the hybrid model the efflux from S and influx to CG given by:
For the interval of one ABM simulation step it is assumed to be constant.
In the ODE model, chemotherapy causes a general reduction of the homogeneous stem cell population. The difference equations model distinguishes between proliferative and nonproliferative cells. Hence, it is possible to model the effect of chemotherapeutic drugs by cellcycle specific toxicities [33]. In analogy to the ODE model, we model chemotherapy as a transient first order loss, but limited to proliferating cells in GE Ω. This is implemented by an intermediate state Ω′_{k,c} in each simulation step:
where Ψ_{CX} is the characteristic chemotherapy function of the regimen considered (s. Appendix A1.1). Subsequently, all other processes such as transition between GEs, cell divisions, differentiation and regeneration are computed from this intermediate state.
Finally, feedback loops 1 and 2 of the ODE model (Figure 1, red and green) remain to be integrated into the hybrid model. Loop 1 is replaced by the intrinsic regulation of the new SCC. To implement the effect of feedback loop 2 in the hybrid model, we considered two alternative mechanisms using the number of granulopoietic cells in the bone marrow as regulator. In analogy to an increased proliferative fraction, the first one is based on an activation of quiescent stem cells by regulation of the transition probability ω. The second increases activity in the SCC by a general acceleration of all processes. Both concepts were tested extensively with various feedback functions. It turned out that this feedback loop, which is of clear importance for the ODE model, can be dropped in the hybrid model without loss of accuracy for the scenarios considered (see Discussion).
We like to emphasize here that constructing the hybrid model we did not perform any parameter fittings, i.e. all parameter values derived from the ODE and ABM model were preserved. This also applies for parameters required only for certain model scenarios such as pharmacokinetic and –dynamic parameters of GCSF applications or chemotherapy toxicities.
Implementation and simulations
The hybrid model was implemented in Simulink 7.4 and Matlab 7.9.0.529 (R2009b). Using the Simulink interface of user defined Sfunctions we implemented the difference equation formulation of the stem cell model as a Matlab Sfunction. Simulations of the hybrid model were carried out on SUSE Linux Enterprise Server 11 (×86 64). All parameters for the ODE model, the difference formulation of the ABM and toxicities of drugs are given in A1.2, A2.2 and A3, respectively.
Results
A number of qualitative and quantitative simulations were performed to study the behavior of our hybrid model in comparison to the former ODE model. If not stated otherwise, we present normalized cell numbers C^{rel} in all figures.
Qualitative model behavior
Depletion of the granulocyte compartment
First, we considered an initial depletion of granulocytes setting and initializing all other compartments with their equilibrium number C^{nor}. We compared the two models regarding dynamics of granulocyte recovery. In both models, the depleted granulocyte compartment is quickly repopulated after about 8 h and equilibrates after a short overcompensation. Equilibrium is reached after about 1d. Since the loss in granulocyte numbers is signaled only via GCSF to PGB and MGB, the feedback to the SCC is not involved. Consequently, the results of hybrid and ODE model are virtually identical (not shown).
Depletion of the granulocyte and maturing granulopoietic precursor compartment
Next we simulated repopulation after a depletion of the two latest stages of granulopoiesis (Figure 3). After initially setting the cell numbers in GRA and all MGB subcompartments to 10^{15} of their equilibrium values and all remaining compartments to C^{nor}, the resulting repopulation dynamics are more complex compared to the previous scenario. In the ODE model, the loss in compartment MGB is transmitted to the SCC (Figure 1, green feedback loop) and results in both, an increased proliferative fraction and a decreased selfrenewal in S (s. Appendix A1.1). This induces an increased efflux and a reduction of S. In consequence, the ODE model quickly repopulates after about 3d, but the feedbacks cause damped oscillations for approximately 80 days. Without the stem cell activation via feedback 2, the hybrid model shows a very similar first response except for the stem cell compartment, but no oscillations occur.
Figure 3. Responses of ODE and hybrid model after depletion of compartments MGB and GRA. In both models the depleted compartments are repopulated within 3d. a) In the ODE model, the feedback to the SCC causes a transient decrease in stem cell number and damped oscillations. b) The hybrid model does not activate the SCC and the perturbation has vanished in all compartments after ~23 days without oscillations.
Simulations of bone marrow transplantation
Because the differences between the models lie in the two SCC, we addressed them directly by simulations mimicking bone marrow transplantation with GCSF support after myeloablative conditioning. For this purpose we initialized all bone marrow compartments with 1% and GRA with 50% of their equilibrium values. In all simulations GCSF was applied until recovery of granulocyte number is achieved. In the ODE model, GCSF was applied for 8 days. Under these conditions, the SCC regenerates 50% of its equilibrium value after 18 days and 100% after 61 days (Figure 4a). Due to the feedback mechanisms small oscillations in stem cell number occur. Repopulation of progenitor compartment CG is similar to S. The more mature compartments PGB and MGB recover faster than stem cells S and progenitors CG. Mature granulocytes GRA recover 25% of their equilibrium value after 4.6 days, 50% after 6 days and repopulate completely after 7 days. If GCSF application is extended to 20 days, recovery dynamics does not change in GRA, but clearly in the SCC despite of no direct effects of GCSF on S. Regeneration in the SCC is accelerated considerably, repopulating 50% after 13.5 days and 100% after only 19 days (data not shown).
Figure 4. Simulations of bone marrow transplantation with GCSF support after myeloablative conditioning. All bone marrow compartments were initialized with 1%, only GRA with 50% of its equilibrium value. We present relative cell numbers throughout. a) In the ODE model, S repopulates 50% after 18d and 100% after 60d. GRA repopulates 50% after 6d and 100% after 7d. b) Both GEs in the SCC of the hybrid model are reduced. Repopulation is slow and not completed after 100d. c) Stem cell reduction is limited to the proliferative GE Ω. S repopulates 98% in 5 days, PGB and GRA 100% after 9.5d. The other compartments repopulate within the next 40 days.
Analogously, we tested the performance of the hybrid model in this scenario (Figure 4b). In both GEs Α and Ω of the SCC the content was reduced to 1% of its equilibrium value. Here, recovery is very slow compared to the ODE model and not completed after 100d. Initially, only GRA and CG decrease clearly. Stem cell number S decreases only slightly, because the cells found at low affinities a in GE Ω rarely switch to Α, but rather leave the SCC. After these cells have left the SCC, the intrinsic regulation of the SCC first repopulates S before it generates progenitor cells again. Initially maximal amplification allows cell numbers in PGB, MGB and GRA to grow. Then the efflux of progenitor cells from S drops to only approximately 0.05% of its equilibrium value and causes a decline of all subsequent cell numbers. The nadir propagates through all compartments eventually reaching GRA, which starts to regenerate at d45.
If only Ω is reduced in the SCC, 82% of the stem cells remain in GE Α, repopulate GE Ω and restore the efflux of progenitors (Figure 4c). Reduced transit times in later departments and increased amplification lead to a repopulation of GRA after only 10 days.
Hence, fundamentally different behavior arose from those three scenarios. While recovery is fast in the ODE model, in the hybrid model it is not achieved within a reasonable time frame, if all stem cells are reduced. If only GE Ω is affected, maximal reduction of GRA is similar to the ODE model, but stem cell dynamics are largely different.
Preliminary adaptations of the hybrid model to this scenario were motivated by the effect of myeloablation on the bone marrow environment (s. Discussion). We altered the transition characteristics of the SCC that model limited capacities of both GEs. Without changing their values at the equilibrium cell numbers, the feedback was made more sensitive by increasing the sensitivity of the sigmoidal transition characteristics to cell number (s. Appendix A2.2). This resulted in much faster repopulation dynamics after reduction of both GEs (s. Appendix, Figure 9). Higher values of the transition characteristics at small cell numbers resulted in stem cell repopulation after 31 days and recovery of 25% of granulocytes after 36 days. In these simulations the fast repopulation dynamics of the stem cell compartment results in oscillations, too. However, we have not elaborated this approach for application to bone marrow transplants, because here we limit ourselves to the combination of the two models as a proof of principle.
Chronic irradiation
After single initial reductions we were interested in the behavior of the system when exposed to a continuous damage as in the case of chronic irradiation. In mice irradiated with a daily dose of up to 0.6 Gy/d this continuous damage results in a strong, dosedependent decline of CFUS numbers during the first days of irradiation. Then CFUS numbers were found to stabilize at a lower level (see [15] and references therein). Cell numbers of later cell stages also decrease in a dosedependent way, but until a dose of about 0.6 Gy/d their cell numbers are less sensitive compared to CFUS. Then they decline rapidly and the animal dies. These results suggest that until a certain threshold of damage, the system can compensate the decline of CFUS. We expect a similar behavior for human hematopoiesis and checked whether the hybrid model can reproduce this behavior, too.
Following Wichmann and Loeffler, chronic irradiation is modeled as a constant first order kill with rates k_{S} in S and k_{CG} = k_{PGB} = 0.66 k_{S} in CG and PGB [15]. In the hybrid model, we consider applications of the kill rate k_{S} to either both GEs or only to the proliferative GE Ω. For each scenario, we recorded the contents of cell compartments after 100d. In the ODE model, S and CG show the same dose dependence (Figure 5a). At small doses, their cell numbers are highly sensitive to small doseincrements. At higher doses the SCC responds to the stem cell loss by higher proliferation and selfrenewal and stabilizes at low stem cell numbers until the system dies for kill rates k_{S} > 0.025 h^{1}. At stem cell kill rates lower than 0.005 h^{1} the rapid loss of S and CG is compensated by increased proliferative fractions and amplifications resulting in only slightly reduced cell counts of circulating granulocytes (GRA). For doses of 0.005 h^{1} < k_{S} < 0.025 h^{1} the loss can no longer compensated and GRA starts to decline rapidly.
Figure 5. Simulations of chronic irradiation. Constant kill rates are applied to S, CG and PGB (with k_{CG} = k_{PGB} = 0.66 k_{S}). Cell numbers in S, CG, PGB and GRA are recorded after 100 days. a) In the ODE model reductions in S and CG are similar. PGB and GRA almost maintain cell numbers at low doses. At higher doses cell numbers in PGB and GRA decrease rapidly. b,c) In the hybrid model damage in S is smaller than in all other compartments. b) If both GEs Α and Ω are damaged, the system is sensitive to kill rates that are an order of magnitude smaller than for the ODE model. c) If damage is restricted to GE Ω, the SCC is more robust at kill rates comparable to those of the ODE model. All other compartments are more sensitive in comparison to stem cells.
In contrast, in the hybrid model reduction in S is smaller than in all other compartments and the greatest cell loss is always found in CG (Figure 5b,c). During the differentiation process the amplifying compartments compensate this loss to some extent, too. If both GEs are damaged, the system is more sensitive to the chronic damage and collapses at kill rates an order of magnitude smaller than the ODE model (Figure 5b). At a kill rate k_{S} ~ 10^{3} h^{1} proliferation in GE Ω fails to counterbalance the losses in both GEs, the SCC cannot stabilize anymore and dies. If only GE Ω is affected, S itself is very robust against irradiation, but the kill in GE Ω minimizes progenitor efflux. Yet at intermediate kill rates, this strong reduction of efflux cannot be compensated at later cell stages (Figure 5c).
Quantitative modeling  application to chemotherapy data
After studying the behavior of the hybrid model qualitatively in comparison to the ODE model, we now apply it quantitatively to clinical data provided by the NHLB trial of the German High Grade NonHodgkinLymphoma Study Group [40,41]. Since one of us (M. Loeffler) is the responsible biostatistician of this trial group, we have access to raw patient data for our modeling purposes. Leukocyte counts were determined in patients receiving one of four CHOPlike multicycle polychemotherapy regimens. The schedules differ in application of etoposide (doseintensification) and duration of therapy cycle (timeintensification). The timeintensified regimens require GCSF support. For comparison with the normalized simulation results leukocyte number is assumed to be proportional to granulocyte number with normal value of 7000 leukocytes per μl blood [16]. Details of therapies can be found in the Appendix, Table 5. The pooled patient data is available online as supplemental data (s. Appendix: Additional file 1).
Additional file 1. Pooled patient data of the NHLB trial of the German High Grade NonHodgkinLymphoma Study Group for the four considered chemotherapy regimen CHOP21, CHOP14, CHOEP21 and CHOEP14. For each of the four regimens, median and quartiles are supplied for each day for the duration of the protocol as given in Appendix, Table 5.
Format: XLSX Size: 32KB Download file
Chemotherapy regimen without GCSF
The hybrid model distinguishes between proliferative and quiescent stem cells. For the mitosisrelated effect of the applied drugs we restrict their toxic effect on proliferating stem cells. The effect of the three cytotoxic drugs applied in the CHOP regimen is modeled by a single set of toxicity parameters, one parameter for each compartment Ω, CG, PGB and MGB. Toxicity parameters were taken from [17]. In particular, we used the same parameter value for toxicity in Ω as for our former stem cell toxicity. The effect of prednisone is modeled as a prolongation of granulocyte halflife according to Bishop et al. [42] and Dale et al. [43]. The results of the simulation of CHOP21 in comparison to clinical data are shown in Figure 6a. After an initial increase of leukocytes due to prednisone application, the cell counts decrease to about 10% of their normal value until day 11, followed by a recovery phase until the start of the next cycle. Model results fit well to the clinical data in the sense that the model prediction is in the interquartile range of the clinical data for almost all time points. For the CHOEP regimen, an additional set of toxicity parameters representing the toxicity of the drug etoposide is used. Again, model predictions fit well to clinical data applying the same toxicity parameters as used in Scholz et al. [17] (Figure 6b).
Figure 6. Comparison of the predictions of the hybrid model with clinical data on leukocyte dynamics in peripheral blood during chemotherapy without GCSF. Clinical data of a) CHOP21 and b) CHOEP21 administration (Blue: median of patients, black: 25 and 75 percentiles) are compared with corresponding simulation results of the hybrid model (red). Simulation results fit well to clinical data in the sense that they lie in the interquartile range of data for almost all time points. Cell numbers are normalized with respect to the average WBC/leukocyte count value of healthy individuals (7000 cells/μl).
Chemotherapy regimen with GCSF
The CHOP14 and CHOEP14 regimen use the same schedule of drug administration, but reduce regeneration time. The damage caused by this intensification is counterbalanced by application of the growth factor GCSF between day 4 and day 13 of each cycle. Therefore we have used the same toxicity parameters as for a cycle duration of 21 days. GCSF applications are covered by the pharmacokinetic and dynamic model of GCSF provided by the ODE model and given in the Appendix A1.1 [17]. Again, for both scenarios a good agreement of hybrid model simulations and data was found for the original toxicity parameters and without other adaptations (Figure 7). Even though a completely different model of stem cell organization and regulation was used in the hybrid model, the resulting leukocyte dynamics under chemotherapy is consistent with both, the ODE model and the clinical data.
Figure 7. Comparison of the predictions of the hybrid model with clinical data on leukocyte dynamics in peripheral blood during chemotherapy with GCSF support. Clinical data of a) CHOP14 and b) CHOEP14 administration (Blue: median of patients, black: 25 and 75 percentiles) are compared with corresponding simulation results of the hybrid model (red). The effect of growth factor support is reflected by the peak approximately one day after starting the GCSF treatment at day 4 in each cycle. As for the regimens without GCSF, simulation results fit well to clinical data. Cell numbers are normalized with respect to the average WBC/leukocyte count value of healthy individuals (7000 cells/μl).
Comparison of stem cell dynamics during chemotherapy
To explain this phenomenon, we had a closer look at the behavior of both SCCs and their effluxes in the course of CHOP21 simulations (Figure 8). During CHOP21 application, total stem cell numbers in ODE and hybrid model are reduced (Figure 8a). But while the ODE SCC is almost eliminated, the SCC of the hybrid model maintains the quiescent population and, thus, more than 80% of its total number due to the limitation of cytotoxicity to proliferating cells. Although proliferative cells are equally reduced in both models, in the hybrid model the quiescent cells quickly repopulate the proliferative GE after the cytotoxic effect of drug administration has passed, and facilitate a nearly complete recovery in each cycle. In contrast, the SCC of the ODE model has to regrow from only 0.23% of its equilibrium population. Therefore it maximizes its proliferative fraction a_{S} and puts all cells to proliferation, but even so it only restores about 20% of its equilibrium population until the beginning of the next cycle.
Figure 8. Comparison of stem cell dynamics in ODE and hybrid model. a) Stem cell number in the ODE model (green) is heavily reduced and does not recover within one cycle. In the hybrid model, total stem cell number (red) is only reduced to ~80%. Depletion of proliferative stem cells (blue) is similar to the ODE model, but regeneration is fast. b) Effluxes from the ODE SCC (green) and the hybrid model (red) almost vanish after drug administration, but increase simultaneously around day 7.
Despite the huge differences in the dynamics of stem cell numbers after chemotherapy, the effluxes of both SCC show similar dynamics (Figure 8b). In the beginning of each cycle both almost vanish and around day 7 both increase rapidly to more than 50% of their equilibrium value. At this time efflux is almost restored in the hybrid model, because GE Ω recovered almost completely. In the ODE model the feedback loops regulating proliferative fraction a_{S} and probability of selfrenewal p allow an increase of efflux to about 50% of its equilibrium value although stem cell number C_{S} is still around 10% of its equilibrium value . In consequence, on the level of mature granulocytes, both models agree with each other and the clinical data as well, because the differences in efflux from the stem cell compartment are sufficiently small to be compensated by the feedback controlled amplification in subsequent bone marrow compartments.
Discussion
The number of attempts to model hematopoiesis or hematopoietic subprocesses is considerable (e.g. [812,1618,4449]). Usually established models are extended to cover new scenarios [45,46,4951]. However, despite the often complementary scenarios covered by yet existing models, little effort was made to unite them in order to show their compatibility/consistency, to seek confirmation or to address open questions. Here, we pursue such a crossvalidation of two established models for hematopoiesis developed in different contexts and on the basis of different experimental data by combining them into one comprehensive hybrid model for human granulopoiesis. The two models were chosen for their range of clinically relevant applications, established parameter sets, in particular for human hematopoiesis, and their complementary level of description. Rather than constructing a model from scratch, we rely on these complementary modeling works by adopting all their model assumptions, equations and parameter settings. Therefore only adaptations necessary to combine them were made. Due to the different contexts in which the models were developed, it was by far not clear whether such an attempt will be successful. Our approach can therefore be understood as a second level of modeling by combining established models in order to construct more comprehensive ones. To solve this task, we examined different hypotheses for model combination and performed a number of qualitative and quantitative benchmark simulations to compare the hybrid model with the established ODE model, to test whether the combined models are compatible and to check whether the new hybrid model has potential to explain a broader range of scenarios.
The ordinary differential equationbased lineage model for granulopoiesis developed by Scholz et al. [17] served as the backbone of our hybrid model. This model describes the dynamics of cell stages from stem cells to circulating granulocytes and the granulopoietic growth factors GCSF and GMCSF. It succeeded in modeling and prediction of a number of clinically relevant scenarios of chemotherapy and GCSF applications [17]. The model of chemotherapy is based on transient cell depletions of bone marrow compartments in dependence on drugs and drug doses. Alternative approaches to model chemotherapy are based on assumptions regarding increased apoptosis and impaired proliferation. We relied on our approach since our former modeling proved to be successful in explaining clinical data sets of different chemotherapies [16,17,19].
Because stem cells are represented as a homogenous population, modeling of processes influenced by stem cell heterogeneity such as cell cycle dependent drug toxicity is difficult in the ODE framework. This was motivation for us to combine it with a more detailed model of hematopoietic stem cells which was developed to understand e.g. clonal competition and immunechemotherapy effects in leukemia treatment [27,28,32]. Combination is realized by replacing the stem cell compartment of the ODE model by a difference equation formulation of the agentbased model of HSC organization [34].
Most of the hematopoietic models published so far do not provide a detailed view on the stem cell level, since they either do not model stem cell dynamics [810,14,48] or treat them similarly simple as the ODE model [11,12]. We chose the ABM for HSC introduced by Roeder and Loeffler [26] to replace the stem cell compartment of our ODE for two reasons: it represents the stem cell population in a highly structured way (two niches, affinity, cell cycle) and has been proven successful in a broad range of scenarios relevant for clinical applications such as Imatinibtreatment of chronic myeloid leukemia in humans [30,32]. Replacement is performed on the basis of a difference equations formulation of the ABM proposed by Kim et al. [34]. Because the difference formulation conserves the ideas and behavior of the ABM, combination of these different types of models (ABM, ODE) at different scales (single cell, cell populations) is also interesting from a general systemsbiological point of view. The stem cell models of Adimy [44] and Mackey [13] describe the stem cell population with less detail, but are continuous alternatives for substitution of the stem cell compartment of our ODE model. We have also experimented with ODE compartment models of proliferating and dormant stem cells in the framework of a murine model of haematopoiesis under chemotherapy with fairly good results (unpublished).
For similar computation times in both submodels, we chose to take advantage of two simplifications of the ABM: (i) a difference formulation, which does not deteriorate accuracy [34] and (ii) the approximation of the probabilistic transitions between GEs by their expectation values. The first simplification impedes the analysis of individual cell fates and clonal contributions, but separate matrices and parameter sets allow modeling of limited numbers of stem cell subpopulations, such as leukemic and normal cells. The second makes the model deterministic and results in continuous cell numbers, but eliminates stochastic variation among realizations that could be interpreted as variation among individuals. However, if desired, stochastic transitions can be easily reestablished by using binomially distributed random numbers.
A major challenge of our approach was to substitute the feedbackloop signaling from more mature bone marrow compartments back to stem cells. In the ODE model this and the internal stem cell feedback control proliferative fraction (a_{S}) and selfrenewal probability (p) in the SCC (s. Appendix A1.1). Because the ABM provides an intrinsic regulation of proliferation and stem cell maintenance, the stem cell feedback is replaced by this internal regulation. To translate the feedback of the bone marrow to stem cells into the hybrid model, we reinterpreted the stem cell parameters a_{S} and p of the ODE model in terms of the new stem cell model. In the ABM the proliferative fraction is given by the fraction N_{Ω}/( N_{Ω+}N_{Α}) of stem cells in GE Ω and is regulated by the transitions between the GEs. These in turn are controlled by the transition intensities α and ω. Accounting for experimental results on activation of dormant HSC after GCSF stimulation [38], we first tested a dependence of the probability ω for a transition Α → Ω on the number of granulopoietic cells in the bone marrow. However, simulations showed a long delay between stem cell activation and increase of progenitor efflux. This delay led to overcompensation at the level of granulocytes and, thus, a clear deviation from the clinical data. Alternatively we considered stem cell activation as an acceleration of all processes in the ABM SCC. The resulting shorter cell cycle and faster differentiation leads to an increase in cell divisions in GE Ω and also progenitor efflux from GE Ω. In combination with the regulation of transition intensity ω, it reduces the overcompensation mentioned above, but it does not eliminate it. Independent of this combination, it shifts the granulocyte nadir after chemotherapy by an earlier increase in production of progenitor cells. This causes additional deviation from the clinical data. Both, regulation of ω and acceleration, were implemented as continuous piecewise functions of the bone marrow cellularity G. In the sensitive range around equilibrium conditions we tested linear and quadratic dependences on bone marrow cellularity G. Outside this range we assumed constant minimal and maximal values. Various ranges and sensitivities of both feedback functions and their combinations were tested, but it turned out that the hybrid model agrees best with the quantitative data after chemotherapy if the native stem cell regulation of the ABM is used without an additional feedback. Therefore, we decided to perform all simulations only with the original stem cell regulation of the ABM.
Except for this elimination of a feedback and the restriction of the effect of chemotherapy to proliferating cells, we did not introduce any modifications of both submodels, in particular with respect to parameter settings. However, this substitution altered the number of model parameters. With the ODE stem cell compartment we eliminated eight parameters of the former ODE model including those of loop 2, but introduced 16 parameters with the ABM model. Because these parameters were neither fitted nor adapted, we added no additional degrees of freedom to achieve a better fit to the data. A fine tuning is possible, but contradicts the idea of combining the models as a proof of compatibility and confirmation. It might be necessary when applying the hybrid model to new and more complex scenarios not yet covered by both submodels, e.g. more complex chemotherapies, diseased granulopoiesis and bonemarrow transplantation (s. Section “Simulations of bone marrow transplantation” and A4), which is planned for the future.
In our qualitative simulations, we observed that after damaging mature compartments, the recovery of circulating cells is similar under the former ODE model and the new hybrid model. But at early cell stages, the ODE model shows damped oscillations which do not occur for the hybrid model. Lack of data prevents us to decide which behavior is in better agreement with biological reality.
Main differences between the hybrid and the ODE model arose after perturbations of the stem cell pool itself. In our ODE model simulations of myeloablative bone marrow transplants, the SCC repopulates completely within approximately 50 days and GRA recovers its equilibrium value after only 7 days. This is well below the lower boundary of the wide range of recovery times for absolute neutrophil count after bone marrow transplantation [52,53]. While repopulation in the ODE model is too fast, in the hybrid model it is far too slow. Only if damage in the SCC is restricted to proliferating cells, reconstitution is completed within reasonable times. However, the latter is not comparable to myeloablation since about 82% of the stem cells remain unaffected. Myeloablation very likely triggers a wide destruction of bone marrow structure and the niche environment as is known for irradiation in mice [5456]. Hence, our model does not cover this scenario so far. It would require adaptations of the model mechanisms that consider changes in environment and signaling context. One approach could be a stronger dependence of transitions on the capacity of the damaged bone marrow niche after conditioning and, consequently, a higher fraction of proliferative cells until the niche is restored. In preliminary simulations we used different transition characteristics with higher values for low cell numbers without changes for equilibrium numbers (s. Appendix A4, Figure 9). Consequently, cells are more likely to switch GEs, which can be interpreted as a reduced ability of each GE to keep HSC in a stable state. Driven by the instable quiescence the stem cells enabled a recovery of 25% of granulocyte number within 36 days. This is close to the upper boundary found in clinical studies [52,53] and demonstrates the potential of the hybrid model for applications to highdose chemotherapy and stem cell transplantation. The transition characteristics of the new stem cell compartment control the activation of stem cells and capacity of the two GEs. They also determine the ratio of dormant and active HSC and their rate of exchange. Hence, a detailed elaboration on changing growth environments and its application to highdose chemotherapy is possible, but goes beyond the scope of this article.
Similar arguments hold for our qualitative simulations of chronic irradiation. But here neither a relation of experimental irradiation dose and model toxicity parameters nor a clear matching of model compartments and CFUS exists, so our results can only be interpreted in a qualitative way. The compensation of stem cell loss by regulated amplification in more committed compartments observed for radiation experiments in mouse [15] is conserved in the hybrid model if one identifies CFUS with the compartment CG.
Despite the differences in stem cell dynamics for bone marrow transplantation and chronic irradiation, quantitative application of the hybrid model to clinical data of four scenarios of conventional chemotherapies with and without GCSF treatment was successful. Here, we restricted the damaging effect of chemotherapy to the proliferating stem cell subcompartment of the hybrid model and applied the stem cell kill rates estimated for the former ODE model. This biologically more plausible approach of cell cycle specific chemotherapeutic drugs in the hybrid model avoids the apparently unrealistic reduction that is predicted by the ODE model [17,33]. We chose the CHOPlike scenarios of the NHLtrial, because they can serve as benchmark scenarios for the modeling of general, more complex chemotherapies [17,19]. Agreement of hybrid model and data was good for all chemotherapy regimens without any adaptation of parameters and only the minute change of model structure that were necessary for the combination of the submodels.
The decision to drop the feedback of bone marrow cellularity to stem cells was based on the comparison with quantitative data of the CHOPlike chemotherapy regimens. These regimens intend to maximally affect tumors at the lowest toxicity for normal tissues. Hence, it is conceivable that under these therapies, the bone marrow environment remains widely intact. Our results therefore suggest that HSCs, in particular the dormant subpopulation, represent an important population that is highly protected and only activated after more severe damage than done by these four damageoptimized regimens. Activation of dormant HSCs only after heavy damage also explains why diseases involving HSCs require myeloablative therapies. For such strong damage, our preliminary model simulations suggest a reduced stability and activation of the dormant state due to this environmental destruction, which is reflected in transition characteristics that are more sensitive to bone marrow cellularity.
The reason for the agreement of the two models for the four CHOPlike regimens is the similarity of efflux dynamics of the two SCC. Both models agree on a proliferative fraction of about 15% of the stem cells at equilibrium, too, but assume different cycle times (8 h in the ODE model, 48 h in the ABM). Proliferative activity in the stem cell compartment is computed from proliferative fraction and cycle time. In both models, it covers a similar range, but follows opposite dynamics during chemotherapy. The ODE SCC activates all stem cells in order to regenerate the nearly eliminated SCC and sends them into proliferation (a_{S} = 1), but is unable to recover within one chemotherapy cycle. In contrast, the hybrid model predicts also a nearly elimination of proliferating HSC (a_{S} ~ 0) by the cytotoxic drugs, but quickly repopulates them from quiescent cells without activating more HSC than in equilibrium. Hence, the two model scenarios suggest contrary hidden dynamics, but both explain the clinical data equally well.
To evaluate which model better resembles biology in the situation at hand, measurements of stem cell numbers during chemotherapy would be required, which is impossible for humans. Even if biopsies would be available for analysis, the problem of identifying quiescent and proliferative HSC subpopulations and their cycle times remains. Generally, the question of the ratio of quiescent and active stem cells and their switching frequency is heavily discussed [38,57]. In the present situation, there are various estimates, but mainly for mouse HSC which seem to differ in many aspects from the human system [58,59]. Hence, for the better availability of experimental data, an adaptation of the model concept to the murine system is planned for the future. The advantage of the hybrid model to distinguish between the two populations of quiescent and proliferating stem cells becomes more powerful for murine hematopoiesis and related experimental possibilities. Here it provides a well suited tool to study the equilibrium of dormant/active HSC and its impact on hematopoiesis. Insights from mouse modeling could later improve our modeling in humans. This could result in a more universal model of human granulopoiesis, whose clinical applications could go beyond existing models. Additionally, the hybrid model will serve as the basis for a comprehensive model of complete hematopoiesis by integrating the model of stem cell lineage commitment and the other models of hematopoietic lineages developed by our group since they rely on the same stem cell model [15,17,18,47].
Conclusions
In our hybrid model we have successfully combined two complementary models for stem cell organization and granulopoiesis by replacing the former stem cell compartment of the lineage model and corresponding regulations. The hybrid model features a novel combination of both, a detailed representation of stem cells and the description of later cell stages of granulopoiesis including the effect of growth factor signaling. In particular, its representation of stem cell states allows the modeling of cell cycle dependent cytotoxic effects of chemotherapy. A number of qualitative and quantitative simulations of the hybrid model revealed very different behaviors of the former and the new stem cell model. Nevertheless, quantitative application of the hybrid model to four CHOPlike chemotherapy regimens showed that the hybrid model explains the clinically observed dynamics of leukocyte numbers equally well. Model validation would require more detailed experimental data on human HSC.
This novel combination paves the way for studying not yet addressed scenarios such as highdose chemotherapy, radiotherapy and diseased hematopoiesis. It serves as another step towards a comprehensive model of hematopoiesis comprising models of stem cell regulation, lineage commitment, bone marrow cell populations, growthfactors and chemotherapy actions.
Appendix
A1. Details of the ODE model
A1.1 Model equations
The regulatory Zfunction
This regulatory function is used to calculate amplification A in CG, PGB and subcompartment G6 in MGB as well as the transit time in MGB. It also determines the endogenous production of the growth factors.
This representation allows regulating the quantity between a minimum and a maximum in dependence on . In steadystate it is equal to Y^{nor}. The parameter b_{Y} allows choosing different degrees of sensitivity. Details can be found in [17].
Characteristic function of chemotherapy
Chemotherapy is modeled in all compartments as a transient depletion of cells following a first order kinetics. The rate of this kinetics is given as the product of a kill rate k_{X} (drug and compartment specific) and a characteristic function Ψ_{CX} modeling its time dependence:
where t_{i} are the time points of chemotherapy application. If drugs are applied at different schedules, the different chemotherapy functions are added.
Probability of selfrenewal in S
The probability of selfrenewal in the stem cell compartment S is controlled by numbers of stem and granulopoietic cells in the bone marrow:
The representation of p by tanh was used to regulate the quantity symmetrically between a minimum and a maximum and to guarantee maximum sensitivity in steadystate. See Wichmann/Loeffler for further details [34].
Proliferative fraction a_{X} in compartments S and CG
Also the proliferative fraction a_{X} in the compartments S and CG is regulated by the numbers of stem and granulopoietic cells in the bone marrow:
This regulation was also adopted from Wichmann/Loeffler p.65 [15] and Scholz et al. [17] respectively. In compartment X (=S or CG) the proliferative fraction a_{X} is regulated between a minimum and a maximum value, and . The quantity x represents the weighted logarithmic relative size of the regulating compartments (S and G). The normal value and the value for intensified stimulation corresponding to x = −ln2 are fixed parameters, too.
Amplification splitting in all proliferating compartments
According to Scholz et al. amplification is always divided into an amplification of the influx (“in”) and the efflux (“out”) respectively (see [17] for further details).
Stem cell compartment S
Granulopoietic progenitor cells CG
Proliferating granulopoietic precursor cells PGB
Maturing granulopoietic precursor cells MGB G46
To guarantee distributed transit times, the compartments G46 are further subdivided into N_{G4}, N_{G5} and N_{G6} subcompartments (for details see [17]). While the maturing time is regulated in all subcompartments, amplification is regulated only in G6 (A < 1) which reflects postmitotic apoptosis, while A = 1 in G4/G5. Below the equation for G4 are given. The corresponding equations for G5/G6 are obtained by replacing G4 with G5/G6 and PBG with G4/G5, respectively.
Granulocytes GRA
Ψ_{Pred} is the characteristic function of prednisone administration analogous to the characteristic function of chemotherapy (A2) but without first cycle effect (i.e. f_{fc} = 1).
GranulocyteMacrophage Colony Stimulating Factor GMCSF
Granulocyte Colony Stimulating Factor GCSF
The model includes two GCSF compartments: C_{GCSF} representing blood concentration and C_{SC} for subcutaneous GCSF administration which is subdivided into C_{SC_1} and C_{SC_2} to model delayed absorption.
Where t_{i} are the times when the infusions start and t_{inf} is the duration of GCSF application. Degradation is modeled by two independent processes, one dependent on GRA and another independent of GRA [17]:
A1.2 Model parameter
Here we provide the full set of parameters established by Scholz et.al. [17] and used in the simulations of the ODE model. The identical set was used without any adaptation in the simulations of the hybrid model.
A1.3 Equilibrium values of all compartments
Evaluation of the parameters given in Table 2 leads to the equilibrium numbers given in Table 3. Please note that these numbers are given in units of equilibrium stem cell number.
A2. Details of the difference model
A2.1 Model equations
Transition functions
The general class of sigmoidal functions used for calculation of transition characteristics and transition intensities is given by:
It is defined by a scaling parameter Ñ and four shape parameter ν_{i}, i = 1,…,4, which can be expressed more intuitively by f (0), f (Ñ_{Α} /2), f (Ñ_{Α}) and f (∞) by
with the auxiliary quantities
A2.2 Model parameter
The set of parameters used in the simulations of the hybrid model is given in Table 4. Please note again that this is exactly the same parameter set as used in the leukemia simulations of the ABM except for the small changes in differentiation coefficient and regeneration coefficient [32].
Table 4. Parameter of the stem cell compartment in the hybrid model
A3. Chemotherapy and toxicity parameter
A3.1 Schedule of the considered regimen
For quantitative simulations we have considered four established chemotherapy regimen that serve as benchmarks, because they comprise dose and timeintensification of therapy. Their schedules and cycle as applied in the NHLB trial [40,41] numbers are provided in Table 5.
Table 5. Details of the four considered CHOPlike chemotherapy regimens
A3.2 Toxicity parameter
Scholz et.al. established the toxicity parameters given in Table 6 for their simulations of the four CHOPlike chemotherapy which are summarized in Table 5[17]. We use exactly the same parameters in our hybrid model simulations, including the first cycle effect (f_{fc} = 1.3).
Table 6. Toxicity parameter for simulations of the CHOPlike regimens
A4 Preliminary simulation of highdose chemotherapy with bone marrow transplant
The damage induced by the CHOPlike chemotherapy regimen appears to be small enough to render the feedback from bone marrow cells to stem cells not only unnecessary, but even distorting. For the stronger damage conferred by myeloablative therapies the response of the hybrid model can be accelerated by an adaptation of the transition characteristics that makes them more sensitive to changing cell numbers (Figure 9). This adaptation does not influence equilibrium cell numbers, because values at equilibrium are kept constant.
Figure 9. Preliminary simulation of myeloablative bone marrow transplant with modified transition characteristics. While keeping the equilibrium values of the transition characteristics nearly unchanged, for smaller cell numbers N_{Α} and N_{Ω} the transition characteristics increase much stronger than for the established human parameter set. The more dynamic switching behavior results effectively in stem cell activation and much faster recovery. Stem cells repopulate completely after 31 days and GRA recovers 25% of its equilibrium value after 36 days. Again GCSF application was continued until the recovery of GRA.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
Model development: AK, IR, ML, MS, Design of simulation scenarios: AK, ML, MS, Model simulations and data analysis AK, Paper writing: AK, MS, Contribution to paper writing: IR, ML. All authors read and approved the final version of the manuscript.
Acknowledgements
AK was funded by a grand of the Federal Ministry of Education and Research of the Federal Republic of Germany ("'Haematosys"', BMBF / PTJ0315452A). M.S. was funded by LIFE  Leipzig Research Center for Civilization Diseases, University of Leipzig. LIFE is funded by means of the European Union, by the European Regional. Development Fund (ERDF) and by means of the Free State of Saxony within the framework of the excellence initiative.
The authors want to thank I. Glauche for helpful discussions. We also like to thank the German High Grade NonHodgkinLymphoma Study Group (chair M. Pfreundschuh) for the kind permission to use the data.
References

Maximow AA: Der Lymphozyt als gemeinsame Stammzelle der verschiedenen Blutelemente in der embryonalen Entwicklung und im postfetalen Leben der Säugetiere.

Ichikawa Y, Pluznik DH, Sachs L: In vitro control of the development of macrophage and granulocyte colonies.
Proc Natl Acad Sci USA 1966, 56:488495. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Till JE, Mc CE: A direct measurement of the radiation sensitivity of normal mouse bone marrow cells.
Radiat Res 1961, 14:213222. PubMed Abstract  Publisher Full Text

Wichmann HE: Computer modeling of Erythropoiesis. In Current Concepts in Erythropoiesis. Edited by Dunn CDR. Chichester: Wiley; 1983:99.

Lajtha LG, Oliver R, Gurney CW: Kinetic model of a bonemarrow stemcell population.
Br J Haematol 1962, 8:442460. PubMed Abstract  Publisher Full Text

Kaushansky K, Lichtman M, Beutler E, Kipps T, Prchal J, Seligsohn U: Williams Hematology. 8th edition. New York: McGrawHill Professional; 2010.

Loeffler M, Birke A, Winton D, Potten C: Somatic mutation, monoclonality and stochastic models of stem cell organization in the intestinal crypt.
J Theor Biol 1993, 160:471491. PubMed Abstract  Publisher Full Text

Crauste F, Demin I, Gandrillon O, Volpert V: Mathematical study of feedback control roles and relevance in stress erythropoiesis.
J Theor Biol 2010, 263:303316. PubMed Abstract  Publisher Full Text

Crauste F, PujoMenjouet L, Genieys S, Molina C, Gandrillon O: Adding selfrenewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis.
J Theor Biol 2008, 250:322338. PubMed Abstract  Publisher Full Text

Vainstein V, Ginosar Y, Shoham M, Ranmar DO, Ianovski A, Agur Z: The complex effect of granulocyte colonystimulating factor on human granulopoiesis analyzed by a new physiologicallybased mathematical model.
J Theor Biol 2005, 234:311327. PubMed Abstract  Publisher Full Text

Manesso E, Teles J, Bryder D, Peterson C: Dynamical modelling of haematopoiesis: an integrated view over the system in homeostasis and under perturbation.
J R Soc Interface 2013, 10:20120817. PubMed Abstract  Publisher Full Text

MarciniakCzochra A, Stiehl T, Ho AD, Jager W, Wagner W: Modeling of asymmetric cell division in hematopoietic stem cells–regulation of selfrenewal is essential for efficient repopulation.
Stem Cells Dev 2009, 18:377385. PubMed Abstract  Publisher Full Text

Mackey MC: Unified hypothesis for the origin of aplastic anemia and periodic hematopoiesis.
Blood 1978, 51:941956. PubMed Abstract  Publisher Full Text

Ostby I, Rusten LS, Kvalheim G, Grottum P: A mathematical model for reconstitution of granulopoiesis after high dose chemotherapy with autologous stem cell transplantation.
J Math Biol 2003, 47:101136. PubMed Abstract  Publisher Full Text

Wichmann H, Loeffler M: Mathematical Modeling of Cell Proliferation: Stem Cell Regulation in Hemopoiesis. BocaRaton, Florida: CRC Press; 1985.

Engel C, Scholz M, Loeffler M: A computational model of human granulopoiesis to simulate the hematotoxic effects of multicycle polychemotherapy.
Blood 2004, 104:23232331. PubMed Abstract  Publisher Full Text

Scholz M, Engel C, Loeffler M: Modelling human granulopoiesis under polychemotherapy with GCSF support.
J Math Biol 2005, 50:397439. PubMed Abstract  Publisher Full Text

Scholz M, Gross A, Loeffler M: A biomathematical model of human thrombopoiesis under chemotherapy.
J Theor Biol 2010, 264:287300. PubMed Abstract  Publisher Full Text

Scholz M, Schirm S, Wetzler M, Engel C, Loeffler M: Pharmacokinetic and dynamic modelling of GCSF derivatives in humans.
Theor Biol Med Model 2012, 9:32. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Bjornson CR, Rietze RL, Reynolds BA, Magli MC, Vescovi AL: Turning brain into blood: a hematopoietic fate adopted by adult neural stem cells in vivo.
Science 1999, 283:534537. PubMed Abstract  Publisher Full Text

Goodell MA, Jackson KA, Majka SM, Mi T, Wang H, Pocius J, Hartley CJ, Majesky MW, Entman ML, Michael LH, Hirschi KK: Stem cell plasticity in muscle and bone marrow.
Ann N Y Acad Sci 2001, 938:208218.
discussion 218–220
PubMed Abstract  Publisher Full Text 
Jackson KA, Mi T, Goodell MA: Hematopoietic potential of stem cells isolated from murine skeletal muscle.
Proc Natl Acad Sci USA 1999, 96:1448214486. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Montfort MJ, Olivares CR, Mulcahy JM, Fleming WH: Adult blood vessels restore host hematopoiesis following lethal irradiation.
Exp Hematol 2002, 30:950956. PubMed Abstract  Publisher Full Text

Moore KA, Lemischka IR: Stem cells and their niches.
Science 2006, 311:18801885. PubMed Abstract  Publisher Full Text

Schofield R: The relationship between the spleen colonyforming cell and the haemopoietic stem cell.
Blood Cells 1978, 4:725. PubMed Abstract

Roeder I, Loeffler M: A novel dynamic model of hematopoietic stem cell organization based on the concept of withintissue plasticity.
Exp Hematol 2002, 30:853861. PubMed Abstract  Publisher Full Text

Roeder I, Braesel K, Lorenz R, Loeffler M: Stem cell fate analysis revisited: interpretation of individual clone dynamics in the light of a new paradigm of stem cell organization.
J Biomed Biotechnol 2007, 2007:84656. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Roeder I, Horn K, Sieburg HB, Cho R, MullerSieburg C, Loeffler M: Characterization and quantification of clonal heterogeneity among hematopoietic stem cells: a modelbased approach.
Blood 2008, 112:48744883. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Glauche I, Thielecke L, Roeder I: Cellular aging leads to functional heterogeneity of hematopoietic stem cells: a modeling perspective.
Aging Cell 2011, 10:457465. PubMed Abstract  Publisher Full Text

Glauche I, Horn K, Horn M, Thielecke L, Essers MA, Trumpp A, Roeder I: Therapy of chronic myeloid leukaemia can benefit from the activation of stem cells: simulation studies of different treatment combinations.
Br J Cancer 2012, 106:17421752. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Horn M, Loeffler M, Roeder I: Mathematical modeling of genesis and treatment of chronic myeloid leukemia.
Cells Tissues Organs 2008, 188:236247. PubMed Abstract  Publisher Full Text

Roeder I, Horn M, Glauche I, Hochhaus A, Mueller MC, Loeffler M: Dynamic modeling of imatinibtreated chronic myeloid leukemia: functional insights and clinical implications.
Nat Med 2006, 12:11811184. PubMed Abstract  Publisher Full Text

Lohrmann HP, Schreml W: Cytotoxic drugs and the granulopoietic system. Berlin. New York: Springer; 1982.

Kim PS, Lee PP, Levy D: Modeling imatinibtreated chronic myelogenous leukemia: reducing the complexity of agentbased models.
Bull Math Biol 2008, 70:728744. PubMed Abstract  Publisher Full Text

Lenhoff S, Rosberg B, Olofsson T: Granulocyte interactions with GMCSF and GCSF secretion by endothelial cells and monocytes.
Eur Cytokine Netw 1999, 10:525532. PubMed Abstract  Publisher Full Text

Lord BI, Bronchud MH, Owens S, Chang J, Howell A, Souza L, Dexter TM: The kinetics of human granulopoiesis following treatment with granulocyte colonystimulating factor in vivo.
Proc Natl Acad Sci USA 1989, 86:94999503. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Schmitz S, Franke H, Brusis J, Wichmann HE: Quantification of the cell kinetic effects of GCSF using a model of human granulopoiesis.
Exp Hematol 1993, 21:755760. PubMed Abstract

Wilson A, Laurenti E, Oser G, van der Wath RC, BlancoBose W, Jaworski M, Offner S, Dunant CF, Eshkind L, Bockamp E, et al.: Hematopoietic stem cells reversibly switch from dormancy to selfrenewal during homeostasis and repair.
Cell 2008, 135:11181129. PubMed Abstract  Publisher Full Text

Potten CS, Loeffler M: Stem cells: attributes, cycles, spirals, pitfalls and uncertainties Lessons for and from the crypt.
Development 1990, 110:10011020. PubMed Abstract  Publisher Full Text

Pfreundschuh M, Trumper L, Kloess M, Schmits R, Feller AC, Rube C, Rudolph C, Reiser M, Hossfeld DK, Eimermacher H, et al.: Twoweekly or 3weekly CHOP chemotherapy with or without etoposide for the treatment of elderly patients with aggressive lymphomas: results of the NHLB2 trial of the DSHNHL.
Blood 2004, 104:634641. PubMed Abstract  Publisher Full Text

Pfreundschuh M, Trumper L, Kloess M, Schmits R, Feller AC, Rudolph C, Reiser M, Hossfeld DK, Metzner B, Hasenclever D, et al.: Twoweekly or 3weekly CHOP chemotherapy with or without etoposide for the treatment of young patients with goodprognosis (normal LDH) aggressive lymphomas: results of the NHLB1 trial of the DSHNHL.
Blood 2004, 104:626633. PubMed Abstract  Publisher Full Text

Bishop CR, Athens JW, Boggs DR, Warner HR, Cartwright GE, Wintrobe MM: Leukokinetic studies. 13. A nonsteadystate kinetic evaluation of the mechanism of cortisoneinduced granulocytosis.
J Clin Invest 1968, 47:249260. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Dale DC, Fauci AS, Wolff SM: Alternateday prednisone. Leukocyte kinetics and susceptibility to infections.
N Engl J Med 1974, 291:11541158. PubMed Abstract  Publisher Full Text

Adimy M, Crauste F, Ruan SG: A mathematical study of the hematopoiesis process with applications to chronic myelogenous leukemia.
Siam J Appl Math 2005, 65:13281352. Publisher Full Text

Brooks G, Provencher G, Lei J, Mackey MC: Neutrophil dynamics after chemotherapy and GCSF: the role of pharmacokinetics in shaping the response.
J Theor Biol 2012, 315:97109. PubMed Abstract  Publisher Full Text

Colijn C, Foley C, Mackey MC: GCSF treatment of canine cyclical neutropenia: a comprehensive mathematical model.
Exp Hematol 2007, 35:898907. PubMed Abstract  Publisher Full Text

Glauche I, Cross M, Loeffler M, Roeder I: Lineage specification of hematopoietic stem cells: mathematical modeling and biological implications.
Stem Cells 2007, 25:17911799. PubMed Abstract  Publisher Full Text

Ostby I, Kvalheim G, Rusten LS, Grottum P: Mathematical modeling of granulocyte reconstitution after highdose chemotherapy with stem cell support: effect of posttransplant GCSF treatment.
J Theor Biol 2004, 231:6983. PubMed Abstract  Publisher Full Text

Zhuge C, Lei J, Mackey MC: Neutrophil dynamics in response to chemotherapy and GCSF.
J Theor Biol 2012, 293:111120. PubMed Abstract  Publisher Full Text

Colijn C, Mackey MC: A mathematical model of hematopoiesis: II Cyclical neutropenia.
J Theor Biol 2005, 237:133146. PubMed Abstract  Publisher Full Text

Colijn C, Mackey MC: A mathematical model of hematopoiesis–I. Periodic chronic myelogenous leukemia.
J Theor Biol 2005, 237:117132. PubMed Abstract  Publisher Full Text

Berger C, Bertz H, Schmoor C, Behringer D, Potthoff K, Mertelsmann R, Finke J: Influence of recombinant human granulocyte colonystimulating factor (filgrastim) on hematopoietic recovery and outcome following allogeneic bone marrow transplantation (BMT) from volunteer unrelated donors.
Bone Marrow Transplant 1999, 23:983990. PubMed Abstract  Publisher Full Text

Ernst P, Bacigalupo A, Ringden O, Ruutu T, Kolb HJ, Lawrinson S, Skacel T: A phase 3, randomized, placebocontrolled trial of filgrastim in patients with haematological malignancies undergoing matchedrelated allogeneic bone marrow transplantation.
Arch Drug Inf 2008, 1:8996. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chute JP, Muramoto GG, Salter AB, Meadows SK, Rickman DW, Chen B, Himburg HA, Chao NJ: Transplantation of vascular endothelial cells mediates the hematopoietic recovery and survival of lethally irradiated mice.
Blood 2007, 109:23652372. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Fliedner TM, Bond VP, Cronkite EP: Structural, cytologic and autoradiographic (H3thymidine) changes in the bone marrow following total body irradiation.
Am J Pathol 1961, 38:599623. PubMed Abstract  PubMed Central Full Text

Shirota T, Tavassoli M: Alterations of bone marrow sinus endothelium induced by ionizing irradiation: implications in the homing of intravenously transplanted marrow cells.
Blood Cells 1992, 18:197214. PubMed Abstract

Takizawa H, Regoes RR, Boddupalli CS, Bonhoeffer S, Manz MG: Dynamic variation in cycling of hematopoietic stem cells in steady state and inflammation.
J Exp Med 2011, 208:273284. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Doulatov S, Notta F, Laurenti E, Dick JE: Hematopoiesis: a human perspective.
Cell Stem Cell 2012, 10:120136. PubMed Abstract  Publisher Full Text

Iwasaki H, Akashi K: Hematopoietic developmental pathways: on cellular basis.
Oncogene 2007, 26:66876696. PubMed Abstract  Publisher Full Text