## Abstract

A very high number of different types of blood cells must be generated daily through a process called haematopoiesis in order to meet the physiological requirements of the organism. All blood cells originate from a population of relatively few haematopoietic stem cells residing in the bone marrow, which give rise to specific progenitors through different lineages. Steady-state dynamics are governed by cell division and commitment rates as well as by population sizes, while feedback components guarantee the restoration of steady-state conditions. In this study, all parameters governing these processes were estimated in a computational model to describe the haematopoietic hierarchy in adult mice. The model consisted of ordinary differential equations and included negative feedback regulation. A combination of literature data, a novel divide et impera approach for steady-state calculations and stochastic optimization allowed one to reduce possible configurations of the system. The model was able to recapitulate the fundamental steady-state features of haematopoiesis and simulate the re-establishment of steady-state conditions after haemorrhage and bone marrow transplantation. This computational approach to the haematopoietic system is novel and provides insight into the dynamics and the nature of possible solutions, with potential applications in both fundamental and clinical research.

## 1. Introduction

All functionally distinct blood cell populations are generated through haematopoiesis arising from the rare haematopoietic stem cells (HSC) located in the bone marrow. HSC have the potential to self-renew and also to differentiate into any haematopoietic cell lineage, ensuring the supply of functionally differentiated mature blood cells [1,2]. Considerable progress has been made in delineating developmental pathways of several lymphoid and myeloid lineages and in identifying transcription factors that establish and maintain their fate (for review, see [3,4]), where HSC and early progenitors make an irrevocable choice between lymphoid and myeloid fates is still under debate (for review, see [5–7]). In the classical model, haematopoiesis is depicted as a hierarchical differentiation tree system where HSC give rise to multipotent progenitors of the two main lineages (myeloid and lymphoid) and then to precursors, committing to multiple or single pathways ultimately leading to the mature cells [8,9]. This unidirectional lineage specification process is thought to proceed by a series of irreversible decisions which commit the cell populations to increasingly differentiated states with decreasing potential for self-renewal, as well as a general and gradual increase in cell numbers [10,11].

The tight control and regulation of all these concurrent processes are fundamental for maintenance of homeostasis as the system responds to external perturbations such as haemorrhage or disease by supplying the necessary number of mature blood cells. The clinical importance of the haematopoietic system, as well as its privileged position as a model for stem cell biology, has inspired the development of mathematical models for more than three decades [12–14]. Since such a large system is not easily amenable to global modelling, most attempts have focused on specific aspects of the haematopoietic hierarchy in order to reduce complexity. A number of models have addressed the dynamical features of specific branches of the system, e.g. the oscillatory behaviour of HSC [15,16], regulation of erythropoiesis [17,18], T-cell production [19] or the effects of ageing on B-cell development [20]. Other models have focused on fundamental structural and functional properties of haematopoiesis, such as the regulation of self-renewal [21,22], the characterization of the stem cell population as the most sensitive to environmental signals out of all haematopoietic cell populations [11], the importance of signalling factors on the decision to self-renew or differentiate, depending, for example, on the number of fully differentiated cells [23,24], the role of feedback control mechanisms [25–27] and lineage specification principles [28]. More global attempts describe haematopoiesis from HSC to mature cells, but mostly implementing simplified architectures, either by not explicitly considering the different lineages and branching points [29], or by doing so at a course grained level. The latter has been particularly useful when modelling clinically relevant phenomena such as periodic haematological diseases [30–32] and response to chemotherapy [33,34]. As a whole, these models have provided insights into the steady-state properties of haematopoiesis, as well as the dynamics of response to perturbations, fundamental in maintaining homeostasis. However, these efforts have focused on details of specific features and populations, rather than a more fine-grained compartmental view of haematopoiesis as a whole.

In this work, we have addressed haematopoiesis from a global perspective and have implemented what, to our knowledge, is the most extensive model of the haematopoietic system, explicitly considering the main branching points and specific lineages from HSCs to the main mature cell type populations in adult mice. This was achieved through (i) exhaustive literature searches; (ii) an innovative constrained optimization strategy (i.e. divide et impera heuristic approach) that allowed the identification of a large number of possible configurations for the steady state; (iii) a series of selection and filtering criteria to determine the best configuration for the steady state; exploration of the parameter space defining the feedback activation by simulating (iv) 10 per cent haemorrhage and (v) bone marrow transplantation (BMT) after irradiation. Our model provided a detailed description of steady-state dynamics in the different lineages and was able to qualitatively capture the observed behaviour of the system under mild (10% haemorrhage) and heavy (BMT after irradiation) perturbations, thus providing insight into the relevance of feedback control mechanisms for the maintenance of homeostasis.

## 2. Methods: a mathematical model for the haematopoietic system

### 2.1. Architecture

Among the various representative structures of lineage branches in haematopoiesis [6,7], the classical haematopoietic hierarchic architecture [8,9] has been chosen as a scaffold for the dynamic model (figure 1). In this structure, long-term haematopoietic stem cells (LT-HSCs) reside at the top of the system, giving rise to short-term haematopoietic stem cells (ST-HSCs), which commit into multipotent progenitors (MPPs). These branch into common myeloid (CMP) and common lymphoid progenitors (CLP). The first commit into granulocyte–monocyte progenitors and precursors (GMPP) which differentiate into granulocytes and monocytes (GM). CMP also give rise to megakaryocyte–erythrocyte progenitors (MEP), which commit into platelets and erythrocytes, through platelet progenitors and precursors (PPP) and erythrocyte progenitors and precursors (EPP), respectively. On the other hand, CLP give rise to naive B lymphocytes, T lymphocytes and natural killer (NK) cells, through the corresponding progenitors and precursors, i.e. BPP, TPP and NKPP. The experimental evidence for the existence of two LT-HSC populations in mice has also been considered [35,36]. Specifically, the first population, dormant long-term haematopoietic stem cells (d-LT-HSC), is characterized by the highest self-renewal potential, but it is activated only in response to haematopoietic stress. The second, active long-term haematopoietic stem cells (a-LT-HSC), gives rise to progenitors and mature cells required for daily maintenance of haematopoiesis.

The final structure of the model is the result of intensive literature searches, compromising between the need for providing a global description of the haematopoiesis process and the limits related to compiling parameter information from a multitude of different experimental settings. This involved minimal and schematic lineages, neglecting some secondary branches and collapsing cascades of different populations of cells. Hence the names assigned to each population considered in the model should be considered as general descriptions. Moreover, the lineage of dendritic cells and naive B, T and NK cells has been omitted. This restriction has been imposed since it is beyond the scope of this study to model the full complexity of the immune system (e.g. distinction between naive, activated and memory cells; antigenic stimulation; presence of adjuvants; etc.), or characterize its steady state.

### 2.2. Population dynamics

Each cell population or compartment is assumed to have a homogeneous behaviour. The homeostasis of each compartment results from the balance between formation, either from commitment of upstream stem cells/progenitors or from cell division and loss, ascribable either to cell death or to commitment to downstream cells. The time evolution for a generic compartment *i* (any of the compartments in figure 1) is then given by
2.1where *N*_{i} (*t*) is the average number of cells in compartment *i* at time *t*; (day^{−1}) is the commitment rate of cells in compartment *i* into cells in compartment *j* and is assumed to depend on both time and number of cells in the receiving compartment, i.e. *N _{j}*; (day

^{−1}) is the net division rate (i.e. division rate corrected by death rate) for compartment

*i*and is assumed to be a function of both time and number of cells in the same compartment; is the commitment rate of cells in compartment

*l*into cells in compartment

*i*and is assumed to depend on both time and number of cells in the receiving compartment, i.e.

*N*.

_{i}In other words, the first term of equation (2.1) is the loss of cells in compartment *i* owing to the commitment in downstream compartments, while the second and third terms represent the formation of cells in compartment *i* ascribable either to net division or commitment from upstream compartments, respectively.

Equation (2.1) is valid for all but mature cells, i.e. platelets (P), erythrocytes (E), granulocytes and monocytes (GM). These cells are assumed to neither divide nor commit and to have a constant death rate. The corresponding equation for mature cells is then
2.2where ; (day^{−1}) is the death rate; the commitment rate (day^{−1}) from the respective progenitors (i.e. ) is assumed to depend on the number of mature cells *N _{i}* as well as on the time

*t*. The detailed dynamical equations for each compartment are included in electronic supplementary material, section Dynamical equations.

### 2.3 Model properties

#### 2.3.1. Feedback control

Communication between cell populations throughout the differentiation process is crucial in maintaining tight control over cell numbers according to physiological requirements, both in steady-state and upon perturbation. A number of secreted factors are known to constitute feedback control mechanisms between populations. Examples of such feedbacks are the regulation of erythrocyte production by erythropoietin [37,38], platelet production by thrombopoietin [39] or leucocyte production by granulocyte colony-stimulating factor [40]. We have assumed a simplified model of negative feedbacks, involving local interactions between compartments for the majority of the system. Upon perturbation, steady state in a given compartment will primarily be re-established by increasing the net division rate. After a critical threshold level is reached in terms of cell numbers, a second response is activated by increasing the commitment rate from the immediately upstream compartment. Piece-wise linear functions were used to formulate these feedbacks (figure 1*a*,*b*). See the electronic supplementary material, section Feedback control for details.

In order to assess the importance of communication between populations in maintaining homeostasis upon perturbation, we have implemented two alternative schemes: no feedbacks and no critical threshold between the first (increase in net division rate) and second (increase in commitment rate) responses. Specifically, in the first case net division and commitment rates are assumed constant, that is independent from the number of cells in the compartment they refer to. In the second case, as soon as the number of cells decrease in one compartment, its net division and commitment rates start to linearly increase simultaneously.

#### 2.3.2. Parameters

The dynamic model has in total 16 equations—one for each compartment. For the steady state, there are 50 parameters (in detail: 16 compartment sizes, 13 net division rates, 18 commitment rates and three death rates). For feedback activations, there are 44 additional parameters related to the step-wise linear functions (in detail: 13 maximum net division rates, 15 maximum commitment rates and 16 critical thresholds).The critical threshold *c*_{i} is one of the parameters describing the feedback control in our model. During perturbation, a compartment will try to re-establish its steady state by increasing first its net division rate. After the critical threshold is reached in terms of cell number, e.g. the current cell number is below a fraction *c _{i}* of its steady-state value, a second response is activated by increasing the commitment rate from the upstream compartment. The critical threshold has been set to 1 for the commitment rates in mature cells (i.e. P, E, GM) since these cells do not have any proliferative capacity. On the other hand, the same parameter has been fixed to 0.95 for the net division rate of d-LT-HSC given their properties of highest self-renewal potential and activation only in response to haematopoietic stress [35,36].

#### 2.3.3. Scaling of commitment rates

In order to support the observed general increase in cell numbers from top to bottom of the haematopoietic system [10], we have assumed increasing commitment rates for the vast majority of compartments going downstream in the system in steady state (see the electronic supplementary material, section Scaling of commitment rates for details). This scaling property is reasonable for all cell populations, but it is not necessarily fulfilled in the lymphoid branch. A possible reason is that while platelets, erythrocytes, granulocytes and monocytes do not proliferate, certain lymphocytes never need to become fully post-mitotic [41,42].

### 2.4. Steady state

Steady-state parameters retrieved from the literature for adult mice are summarized in table 1. In particular, the net division rate for d-LT-HSC () has been fixed to zero since we have assumed that the death rate equals the division rate given the negligible turnover rate of those cells in steady state (that is, one division every 171 days [45]) compared with the lifetime of a mouse.

The information retrieved from the literature reduced the number of unknown parameters from 50 to 32. Solving the 16 dynamic equations of the model in steady state (i.e. equations (2.1) and (2.2) equal to zero) in 32 unknowns clearly leads to a multitude of solutions. The exploration of the parameter space has been interpreted as an optimization problem and faced, in the first instance, with the simulated annealing algorithm. As the scaling property consists in constrains between unknown parameters, the parameter space continuously changes during a simulation, making this a very challenging problem for this optimization algorithm (see the electronic supplementary material, section Exploration of the steady-state solution space for details). This hurdle was overcome through the implementation of a divide et impera heuristic approach which consists in (see the electronic supplementary material, section Exploration of the steady-state solution space for details): (i) dividing the haematopoietic system in three modules, namely ROOT, MYELOID BRANCH and LYMPHOID BRANCH (figure 1), (ii) exploring the parameter space in each module by taking advantage of the scaling property. This step led to different possible configurations for each module. (iii) A merging algorithm was then used to connect the independently gained parameter configurations from the three modules (see the electronic supplementary material, section Exploration of the steady-state solution for details). The configuration with (i) the smallest value of merging error, defined as the sum of the squares of the differences between predicted values of the dynamical equations and zero, i.e. their steady-state value, (ii) respecting the scaling property and (iii) not exceeding the total number of BM cells (mature cells and TPP not included in the counting) has been selected as the best configuration of the steady state. The robustness of the best configuration for the steady state was tested (see the electronic supplementary material, section Robustness for details).

### 2.5. Perturbations: haemorrhage and bone marrow transplantation after irradiation

Once established the best configuration for the steady state (see the electronic supplementary material, figure S5), the parameter space defining feedback activation has been explored by simulation of the recovery from haemorrhage and BMT after irradiation with random values for these parameters. Only the configurations of parameters that allowed the recovery within a biologically reasonable time window have been considered. In detail, for the 13 parameters corresponding to maximum net division rates (), three possible options have been taken into account, i.e. 2 day^{−1}; 3 day^{−1}; 2 day^{−1} for d-LT-HSC, a-LT-HSC, ST-HSC, MPP, CMP, CLP and 3 day^{−1} for MEP, PPP, EPP, GMPP, BPP, TPP and NKPP (i.e. ). The 15 parameters for the maximum commitment rates (*a _{i}*

^{C}) have been uniformly randomly generated in the range [0,20] day

^{−1}with the Latin hypercube sampling technique [62] respecting the scaling property where applied, similarly to what has been done for the steady state. Finally, the 12 critical thresholds (

*c*) for a-LT-HSC, ST-HSC, MPP, CMP, MEP, PPP, EPP, GMPP, CLP, BPP, TPP and NKPP have been uniformly randomly generated in the range [

_{i}*x*,1], where

*x*has been fixed either to 0.5, 0.7 or 0.9, scaling the parameter from high to low values while moving upstream in the system (similarly to what was done in steady state for the commitment rates where the scaling property was valid) in order to capture the notion that, given their multipotent capacity, cells in upstream compartments should be as preserved as possible. The use of these different ranges allowed a wide coverage of the parameter space. All the possible combinations of these options have been applied to simulate a 10 per cent blood loss (in total 1.35 × 10

^{5}simulations; see the electronic supplementary material, table S1 for the initial conditions for the compartment sizes). Among these combinations the ones (7067) showing a recovery time (defined as the time elapsed between time 0 and the time when 98 per cent of the steady-state compartment size is reached) of two weeks [48] without the involvement of the d-LT-HSC have been selected to simulate the haematopoietic reconstruction by BMT after irradiation. In detail, two scenarios have been explored: the first one shows that the transplantation of only HSC and MPP leads to thrombocytopenia and anaemia (considered as loss of platelets and erythrocytes higher than 15% with respect to homeostasis), the second allows for the recovery of steady-state values provided by the transplantation of myeloid progenitors together with HSC and MPP (see the electronic supplementary material, table S1 for the initial conditions for the compartment sizes) [63]. Only 26 configurations among the 7067 valid for the 10 per cent haemorrhage simulation have been able to represent the second scenario.

All methods described above are summarized in the flowchart pictured in figure 2.

## 3. Results

### 3.1. Dynamical properties in steady state

The best solutions were used to characterize some of the steady-state features of the system. Taking into account the mean transit times between compartments, the inverse of the corresponding commitment rates versus number of cells in each compartment (figure 3), a number of dynamical principles can be observed:

— The mean transit time decreases with differentiation where the scale property has been applied (i.e. in all the hierarchy, except the lymphoid branch). Transition from a-LT-HSC into ST-HSC is by far the longest (over 30 days). A key implication of these findings is the necessity to co-transplant more rapidly differentiating precursor cells to allow sufficient recovery following myeloablation, which has experimental support [64];

— The mean transition time from a-LT-HSC to any terminal mature compartment is never less than 45 days, with differentiation in the MYELOID faster than in the LYMPHOID BRANCH. This finding is compatible with the reconstitution kinetics of highly purified HSC [44];

— Terminal maturation in the erythroid branch involves a fast and dramatic increase in cell numbers in the latest stages, particularly, between MEP and the mature P and E (transiently through PPP and EPP, respectively). Again, this is supported experimentally by the recent observation that such later stages of erythroid development are characterized by dramatically enhanced cell cycling [65];

— Maturation in the LYMPHOID BRANCH is more time-demanding—the extreme case being the transition from CLP into TPP (around 2000 days). This value was calculated from parameters in the literature [43,50,51,53]. It is considerably higher than the average lifespan of a mouse (700–800 days). These results challenge the notion that the CLP would be a required and sufficient intermediate precursor of the T-cell lineage [66–69].

### 3.2. Perturbations allow for pruning of solution space and provide insight into feedback control

Steady-state conditions reflect the homeostatic behaviour of the system, but represent an incomplete picture of the process. Organisms are subject to different environmental perturbations and the production of blood cells must respond to these changes accordingly to maintain homeostasis. In order to reflect this structural feature of the system, feedbacks were implemented between a compartment size and its net division rate, as well as the commitment rate from the compartment immediately upstream (figure 1*a*,*b*). Upon perturbation on a given compartment, the net division rate increases up to a *plateau* level which then signals for the increase in the commitment rate. The commitment and net division rate maximum values, as well as the critical threshold for commitment were optimized, with the single constraint being the scaling of maximum commitment rates, from higher to lower values when moving upstream on the system (as was performed for steady-state conditions).

The initial optimization step yielded a high number of possible solutions satisfying the constraint. In order to prune the solution space, we considered two perturbation scenarios, reproducing experimental results:

(a) Simulation of 10 per cent blood loss (mild haemorrhage): only solutions that allowed for full recovery on a two-week period [48] were selected for subsequent analysis;

(b) BMT after irradiation: only solutions allowing for recovery of steady-state values without anaemia and thrombocytopenia, as provided by transplantation of myeloid progenitors together with HSC and MPP, were selected [63].

From the initial pool of 1.35 × 10^{5} solutions, 7067 were selected in step (a) and a final set of 26 was reached after step (b).

The solution space for the non-steady-state configurations was scrutinized in more detail (figure 4*a*). The effect of sequential selection of good solutions using the two perturbation scenarios was evident in the maximum commitment rate parameter (middle plot). The final 26 configurations (in blue) were located at lower values compared with the initial parameter sets (in black) and the 7067 configurations respecting the two-week recovery time for a 10 per cent blood loss (in red). In other words, to recapitulate biological behaviour under haemorrhage and BMT after irradiation, maximum commitment rates must be tightly controlled and at low levels, otherwise there is a risk of depletion of the primordial compartments. Analysing in detail the 26 final configurations (figure 4*b*), the critical threshold seems to be the more flexible parameter with respect to the maximum values for net division and commitment rates (bottom plot).

#### 3.2.1. Response to mild haemorrhage is confined to the myeloid branch

Haemorrhage conditions were simulated by reducing by 10 per cent the steady-state cell number values of all terminal compartments in the MYELOID BRANCH (P, E and GM) and allowing the system to evolve until a new equilibrium was reached. This process was repeated for all 26 possible solutions in order to obtain the average behaviour for the system (figure 5). Full recovery (defined as the time required, after the perturbation (time zero), to reach the 98 per cent of steady-state compartment size) required the involvement of the immediately upstream compartments for all cell types (PPP, EPP and GMPP for P, E and GM; figure 5*a*). Recovery of E required the greatest participation of the progenitor cells, with approximately 4 per cent of EPP loss. Relative losses of 2 per cent and 1 per cent were observed for the PPP and GMPP compartments, respectively. Consistent with these observations, recovery times were longer for the E compartment (15 days) than the P (6 days) or GM (1 day) compartments (figure 5*b*,*c*). If we consider blood donation as an induced form of mild haemorrhage, these results are in accordance with data available in healthy adult male humans, where donated plasma is replaced after 2–3 days while red blood cells are restored at a slower rate, after around 36 days [70]. A more detailed description of the response could be obtained by analysing the changes in net division and commitment rates for the different compartments over time. Upon the initial perturbation, commitment rates from the progenitor into mature compartments increased to compensate the 10 per cent loss (figure 5*d*, left panel). In absolute terms, the commitment rate into the P compartment was the highest, followed by commitment into E and finally commitment into GM. However, and more significantly, the largest deviation from steady-state values was observed for the E compartment, whereas commitment rates for P and GM showed minor change. This could be explained by the significant difference between steady-state compartment sizes (see the electronic supplementary material, figure S5), with E being approximately one order of magnitude larger than P and four orders of magnitude larger than GM. Consequently, a 10 per cent decrease in compartment sizes required a more pronounced response from EPP than PPP or GMPP. This effect also propagated into the net division rates in these compartments (figure 5*d*, middle panel). EPP net division increased significantly as a response, whereas in PPP and GMPP changes were minimal. The difference between EPP and PPP net division changes could also be attributed to the steady-state commitment rates feeding from their common progenitor, MEP, with commitment into PPP almost twice as fast as commitment into EPP (figure 5*d*, right panel), allowing faster re-establishment of this compartment. Very importantly, commitment rates into these compartments remained constant at steady-state levels and all changes were limited to the lower levels of the MEYLOID BRANCH, with MEP and CMP being unaltered throughout the simulation.

Contrary to the 14-day recovery time for haemorrhage, this was not a selection criterion when screening for solutions, and could not have been predicted beforehand. It is an important feature of the system, which exhibited robustness within the context of a significant perturbation and contained its effects to only the most differentiated compartments. Compartments at the ROOT, where perturbations would have more serious consequences for the organism, were not perturbed at all.

We then simulated increasingly more pronounced haemorrhage scenarios: for 15 per cent blood loss the system still recovered within reasonable time (see the electronic supplementary material, figure S7); for 30 per cent, recovery was observed but only after 40 days (see the electronic supplementary material, figure S8); for 50 per cent the system did not recover within the 90-day simulation period (see the electronic supplementary material, figure S9).

#### 3.2.2. Co-transplantation of myeloid progenitor cells is necessary in order to reestablish homeostasis after irradiation

Bone marrow irradiation is an extreme perturbation which broadly wipes out the progenitor populations and constitutes, in conjunction with external BMT, a reboot to the haematopoietic system. The levels of erythrocytes and platelets may decrease to life-threatening levels, causing anaemia and thrombocytopenia. The process was simulated by reducing the levels of all compartments to zero except the HSC (d-LT-HSC, a-LT-HSC, ST-HSC), which were kept at 0.1 per cent of the steady-state value (simulating BMT), the MPP, which were kept at 1 per cent of the steady-state value and the mature compartments of the MYELOID BRANCH (P, E and GM), which were not affected (scenario 1). An alternative protocol was to include a fraction of myeloid progenitor cells, as described by Nakorn *et al.* [63] (scenario 2). In both scenarios, the system was evolved for all 26 possible solutions in order to obtain the average behaviour (figure 6). The transplantation of only HSC and MPP as simulated in scenario 1 did not allow for the sustained maintenance of E and P levels, leading to anaemia and thrombocytopenia within 10 days after transplantation (figure 6*a* and electronic supplementary material, figures S10–S11). The recovery of platelets within 30 days in this scenario, as emerging from the 95% confidence interval (calculated as: ), can be misleading, but the loss of platelets and erythrocytes within 10 days after transplantation is in any case incompatible with survival.

On the other hand, the beneficial effects of transplanting myeloid progenitors together with HSC were clearly observable, with levels of P and E suffering minor changes. Importantly, we tested different levels of myeloid progenitor transplantation and observed that the minimum number of progenitor cells that allow for recovery of steady-state values, in this context without anaemia and thrombocytopenia, is 20 per cent of the steady-state value. Considering scenario 2, there was barely any loss of E and P cells, although the recovery time for the latter compartment was longer than the first and still of approximately 7 days (figure 6*b*,*c*). The longer recovery of platelets was because of a faster and more dramatic loss (see the electronic supplementary material, figure S11B) related to their higher death rate (see the electronic supplementary material, figure S5) compared with erythrocytes.

### 3.3. Feedback mechanisms are fundamental for homeostasis

We have assessed the importance of feedback mechanisms in maintaining homeostasis upon perturbations, by testing two alternatives to our minimal feedback model: no critical threshold between increase in net division and commitment rates; and no feedbacks altogether, with net division and commitment rates never increasing from steady-state values.

#### 3.3.1. The haemorrhage scenario

Here feedbacks allowed for the fast re-establishment of the E, P and GM compartments with minimal involvement of the respective progenitors, and without perturbations at the ROOT level (figure 7*a*,*b*). When no feedback responses were in place, changes were still limited to the mature compartments and their progenitors. However, since recovery depended on steady-state rates, re-establishment of the E compartment was achieved only after 72 days, incompatible with the physiological needs of the organism. Finally, when no critical threshold between feedback responses was enforced, we observed that the initial 10 per cent perturbation propagates throughout the system, displaying significant effects up to the MPP compartment (figure 7*a*). As a consequence of this significant disruption, homeostasis could not be reached, and steady-state values were not recovered within the 90-day period of the simulation, with the exception of the GMPP and GM compartments (figure 7*b*).

#### 3.3.2. The bone marrow transplantation after irradiation scenario

For this more extreme perturbation scenario, the feedback response also allowed for the recovery of the system, avoiding anaemia and thrombocytopenia (see the electronic supplementary material, figure S12A) and supporting the recovery of the system to steady-state levels within physiologically reasonable times (see the electronic supplementary material, figure S12B,C). Absence of feedbacks, on the other hand, led to anaemia and thrombocytopenia, since steady-state commitment and net division rates could not replenish the E and P compartments, depleted through cell death (see the electronic supplementary material, figure S12A). That was also the case for all other compartments in the system, since none of the 26 solutions allowed for recovery of steady-state values within 90 days (see the electronic supplementary material, figure S12B), evidencing that repopulation after transplantation was not possible without feedback control systems. Feedback responses without critical threshold also lead to anaemia and thrombocytopenia (see the electronic supplementary material, figure S12A), but in this case partial repopulation was observed. Nevertheless, this feedback scheme could not support the normal re-establishment of the myeloid branch, where (with the exception of GMPP and GM) less than 50 per cent of the 26 solutions allowed for full recovery (see the electronic supplementary material, figure S12B), and never within physiologically reasonable times (see the electronic supplementary material, figure S12C).

## 4. Discussion

The fundamental properties of differentiation, multipotentiality and self-renewal were described for HSC more than 40 years ago [71]. The considerable efforts on basic and clinical research that followed have made haematopoiesis a privileged system for the study of tissue-specific stem cell biology. Numerous mathematical models have been developed over the years, providing an understanding on the dynamical properties of haematopoiesis [12–34]. However, owing to the size and complexity of the system, these models have mostly focused on particular aspects and do not allow a more encompassing analysis over the whole system.

In this work, we have designed and implemented what is, to our knowledge, the most comprehensive global model of haematopoiesis, encompassing all the main branches and bifurcations. We adopted the classical model of the haematopoietic structure and compiled extensive experimental data from the literature in order to define parameters as accurately as possible. We also took into account a minimal set of assumptions regarding structural and functional principles of not only the haematopoietic, but also adult stem cell systems in general.

The model allowed for a characterization of steady-state properties of the system. Additionally, we assumed the existence of negative feedback loops enabling communication between the different compartments, essential in maintaining homeostasis of the system during perturbation events. The molecular mechanisms and functional implications underlying these feedbacks have been modelled in detail for the haematopoietic [25–27] and other stem cell systems [24]. Here, we assumed a common minimal model of local communication for the majority of the structure. Despite the need for further experimental validation, this feedback scheme provided robustness to the system, allowing recovery after perturbation and preserving the most undifferentiated compartments at the top of the hierarchy. These desirable properties were not observed when feedbacks were not present, or the critical threshold was not enforced (net division and commitment rates increased simultaneously). Furthermore, the final set of optimal feedback parameters provided some insight into the fundamental properties and limits of the response, namely the need for tight regulation and low values of maximum commitment rates in order to maintain homeostasis.

Haemorrhage and BMT after irradiation were used as perturbation scenarios at a second stage of our modelling process. Although fundamental in pruning the solution space, the study of these scenarios did not in itself influence the model design, and it is thus not trivial that solutions reproducing biological behaviour could be reached for both scenarios. Also, selection of solutions was based solely on the ability of the system to recover, and not on particular routes for that recovery. The best solutions provided us with a detailed description of the dynamical response to perturbations over the system and allowed us to estimate times of recovery and maximum loss for all compartments. Interestingly, for a mild perturbation such as 10 per cent blood loss, perturbations were only observed on the lower compartments of the myeloid branch, reinforcing the notion of system robustness. For BMT after irradiation, we confirmed the need for co-transplantation of myeloid progenitor cells together with HSC and MPP in order to assure recovery of steady-state values. Even with limitations in our model that does not consider factors, such as contaminations, saturation of niches or graft-versus-host effects, we could also predict the minimum fraction of progenitor cells able to confer this property.

Using relevant clinical scenarios, such as haemorrhage and BMT after irradiation, we were able to reduce the possible solution set from 7067 to 26. Although this is a dramatic reduction, we recognize that having 26 possible configurations is still not ideal. Nevertheless, taking into account the complexity of the problem, these results can be considered encouraging and provide a very promising global framework for the dynamical properties of haematopoiesis.

Another fundamental question in this study was the structure of the haematopoietic system itself. The definition of different compartments is not straightforward and depends on operative and functional criteria that frequently change within the community, as new experimental methods are developed. Although the resolution of these methods tends to increase, leading to the description of more sub-compartments, we have purposefully addressed the problem within what we believe to be a reasonable level of detail. We have taken into account major compartments which in some instances could be further subdivided, but would raise more problems from the modelling perspective. In recent years, the notion that HSC make an irrevocable choice between myeloid and lymphoid lineages early on has been increasingly questioned, and a number of alternatives to the classical model of haematopoiesis have been suggested (for review, see [6,7]). In 2005, Adolfsson *et al.* [72] proposed a revised version of the classical haematopoietic hierarchy where pluripotent HSC develop into lymphoid primed multipotent progenitors (LMPP) upon loss of megakaryocyte and erythrocyte potential, and subsequently generate CLP upon loss of granulocyte–monocyte potential. We have performed an exploratory comparative analysis by implementing the corresponding changes in our model architecture. The most significant difference was the fact that mean transit times for most lineages were significantly smaller in the LMPP model, with differences between 5 and 15 days (data not shown). Although preliminary, these results support the view that our global model can provide a ground for more detailed comparisons between architectures in the future.

We have presented here a comprehensive global model of the haematopoietic system, compiling extensive experimental data and implementing minimal assumptions regarding fundamental properties of stem cell systems. Our model was able to describe and predict fundamental features of the system both in steady-state and under perturbation scenarios, recapitulating observed behaviours for haemorrhage and BMT after irradiation. In doing so, while at the same time being unconstrained by any particular branch or property of the system at the design level, we believe that this model constitutes a general framework on top of which new models can be built, addressing more specific questions and providing new insight regarding the dynamical principles of haematopoiesis and their application to fundamental and clinical research.

## Acknowledgements

This work was in part funded by the Swedish Foundation for Strategic Research Senior Individual grant (A3 06:215) and the Swedish Research Council (Vr 621-2008-3074). J.T. is thankful to the PhD Program in Computational Biology at Instituto Gulbenkian de Ciencia, Oeiras, Portugal and was financially supported by Fundacao para a Ciencia e Tecnologia (SFRH/BD/33208/2007). We are indebted to Ellen Rothenberg at Caltech for valuable insights on properties of the B- and T-cell branches. We also thank our colleagues in the Department of Astronomy and Theoretical Physics at Lund University for their useful comments.

- Received October 5, 2012.
- Accepted November 27, 2012.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.