The gene regulatory network of HSCs
As a result of our manual search for information on the processes that regulate HSCs’ differentiation, we integrate a network of 200 nodes (Supplementary Material S1). We reduce this network by using a network reducing algorithm (41) into a 21 nodes GRN of HSCs that encompasses gene expression, signaling, and metabolism of HSCs (Supplementary Material S2 and Fig. 2). Concerning to gene expression of HSCs, the GRN of HSCs includes the following transcription factors: HIF-1, FOXO3a, p53, Runx1, Meis1, Mef2c, CEBPA GATA1, GATA2, PU.1, Ikzf1 and Gfi1. It has been reported that HIF-1 regulates the quiescence of HSCs as well as ROS homeostasis and metabolic adaptation to the hypoxic conditions of the hematopoietic niche (57). On the other hand, FOXO3a is important to sustain the regenerative properties of HSCs, since the loss of function of this protein decreases self-regeneration of these cell types (1,58). In this sense, FOXO3a is also implicated with p53 in the maintaining of HSCs’ quiescence as well as regulation of ROS intracellular levels (1,58). Selective repression of GATA1, GATA2 and PU.1 conditions the progression in the differentiation pathway of HSCs towards MEP and GMP cell types (59–62). Furthermore, CLP differentiation is controlled by Ikzf1 and Gfi-1 (63,64). Similarly, it has been observed that the differential expression of p53, Gfi1 (58,65), PU.1 and GATA1 (63,64) directly affects lymphopoiesis, either increasing or reducing its speed.
Experimental evidence shows that metabolic factors are essential to alter gene expression patterns that control HSCs differentiation (66,67). Indeed, these cells show a glycolysis-based metabolism whereas differentiation of such cells shifts towards an oxidative metabolism in order to fulfill cellular energetic demands (8). To be precise, it has been seen that AMPK (68), AKT and mTOR (69) exert critical functions in the regulation of the metabolic switch necessary to initiate the differentiation of HSCs. On the other hand, oxidative phosphorylation (OXPHOS) is one of the processes strongly regulated by the metabolic switch that controls the HSCs differentiation HSCs (70). This occurs because OXPHOS is a constant source of chemical energy in the form of ATP and ROS (71), which are suppressed by enzymes such as superoxide dismutase (SOD) (19) and catalases (72). Consequently, in this model, we represented as nodes the enzymes AKT, AMPK, mTOR and SOD, as well as OXPHOS and ROS in the form of H2O2 and the superoxide ion (O2−). The Antioxidants node represents the cellular ROS scavenger system. Collectively, the nodes that compose this GRN (Fig. 2) represent fundamental molecules that regulate the differentiation of HSCs.
The network conditions the differentiation dynamics of HSCs
Boolean formalisms were used to represent and computationally study all interactions described by the GRN of HSCs (Fig. 2 and Table 1). Next, the Boolean model was implemented in BoolNet (42) to study its dynamics. As a result of the above, the model was able to recover steady states equivalent to HSC, MEP, GMP, and CLP cell types (Fig. 3). The MEP, GMP and CLP attractors show mTOR, AKT, H2O2, and O2− activity as expected from active common progenitors (7,66,67). Additionally, MEP steady states show expression of the lineage marker GATA1 (44). GMP attractors are characterized by PU.1 expression (45) while CLP steady states exhibit Gfi1 activity (46). Concerning to HSCs, the model predicted that such cells will present high expression levels of GATA2 (43), FOXO3a, p53 and Meis1, as well as high activity of AMPK and the absence of ROS. Furthermore, this equilibrium can only be reached in the absence of oxygen, as has been previously reported (6). Therefore, the model is congruent with the reported knowledge about gene expression patterns and intracellular activity of HSC differentiation fates. To evaluate if the dynamic behavior observed in our model was independent of the mathematical formalism used it was necessary to transform the Boolean model into ordinary differential equations (ODEs) as is described in García-Gómez et al., (2017). After solving the continuous model, we obtained the equilibrium frequencies of each hematopoietic lineage (Supplementary Material S4), which were congruent with our previous outcomes, indicating that the steady states and the frequencies of common progenitors emerge as a result of restrictions imposed by the regulatory network and independently of the mathematical formalism utilized.
A systems-level approach to study hematopoiesis solves experimental discrepancies
The standard validation procedure of Boolean models is mutant simulations (41,73–75). This procedure evaluates the predictive ability of the model through a series of gain-of-function (GOF) and loss-of-function (LOF) tests. To this end, information was sought on the effect of mutations on the differentiation of HSCs. Next, we performed simulations to recreate the conditions described in the literature, and the effect of such gene alterations on the model was determined (see Methods). As a result of this procedure, we found that the LOF analysis (Table 2) and the GOF analysis (Table 3) showed a high level of consistency with respect to the experimental observations, thus proving that the model is able to recapitulate the mutant phenotypes observed experimentally.
Table 2
Outcomes of the Loss-of-function analysis of the different genes that participate in HSCs differentiation.
Node | Simulation | Experimental data | References | Congruency* |
FOXO3a | HSC basin of attraction is reduced | Compartment reduction of HSCs | (1) | Yes |
CEBPA | All attractors are recovered | GMP differentiation impaired | (156) | No |
Gfi1 | CLP attractors expansion | CLP decrease | (46) | Yes |
PU.1 | Loss of GMP attractor | Myeloid differentiation impairment | (157) | Yes |
GATA1 | MEP attractor is lost | Erythrocytes’ differentiation is blocked | (60) | Yes |
GATA2 | Only one HSC attractor is recovered | Loss of quiescence | (158) | Yes |
p53 | Neither of HSC attractors are recovered | Loss of quiescence | (58) | Yes |
Meis1 | Neither of HSC attractors are recovered | Loss of quiescence | (159) | Yes |
Runx1 | Loss of two MEP attractors | Erythrocytes’ differentiation is impaired | (160) | Yes |
HIF1 | Neither of HSC attractors are recovered | Loss of quiescence | (77) | Yes |
Ikzf1 | CLP attractors are not recovered | CLP differentiation is blocked | (161) | Yes |
Mef2c | Two GMP attractors are lost | Slight diminution of GMP numbers | (162) | Yes |
mTOR | HSC basin of attraction is augmented | Differentiation is blocked | (163) | Yes |
AMPK | HSC basin of attraction is diminished | Differentiation is enhanced | (9) | Yes |
AKT | Two GMP attractors are lost | Proliferation is inhibited | (9) | Yes |
H2O2 | HSC basin of attraction is augmented | Loss of quiescence | (164) | Yes |
O2- | Reduction in HSC basin of attraction | HIF-1 activity is impaired | (165) | Yes |
SOD | Reduction in HSC basin of attraction | Loss of quiescence | (1) | Yes |
Antioxidants | Neither of HSC attractors are recovered | Loss of quiescence | (1) | Yes |
OXPHOS | All attractors are recovered | Differentiation is blocked | (8) | No |
Oxygen | HSC basin of attraction is augmented | Quiescence is enhanced | (77) | Yes |
* This refers to the congruence that exists between the simulations vs. the experimental data. |
Table 3
Outcomes of Gain-of-function analysis
Node | Simulation | Experimental data | References | Congruency* |
FOXO3a | HSC basin of attraction is augmented | Enhanced quiescence and antioxidants activity | (166) | Yes |
CEBPA | Loss of HSC attractors | Myeloid differentiation is enhanced | (167) | Yes |
Gfi1 | Loss of HSC attractors | Differentiation is enhanced | (168) | Yes |
PU.1 | Only Myeloid attractors are recovered | Myeloid differentiation is enhanced | (169) | Yes |
GATA1 | Only MEP attractors are recovered | Erythrocyte differentiation enhanced | (23) | Yes |
GATA2 | Fewer MEP attractors recovered | MEP differentiation impaired | (170) | Yes |
p53 | HSC basin of attraction is augmented | Enhanced quiescence and antioxidants activity | (171) | Yes |
Meis1 | Fewer MEP attractors recovered | Erythrocyte differentiation impaired | (172) | Yes |
Runx1 | Fewer MEP attractors recovered | Erythrocytes’ differentiation impaired | (173) | Yes |
HIF1 | HSC basin of attraction is augmented | Enhanced quiescence and antioxidants activity | (77) | Yes |
Ikzf1 | Only CLP attractors recovered | CLP differentiation enhanced | (161) | Yes |
Mef2c | MEP attractors are recovered | GATA1 expression is reduced | (174) | No |
mTOR | HSC basin of attraction is reduced | Enhanced differentiation | (163) | Yes |
AMPK | HSC basin of attraction is augmented | Augmented quiescence | (9) | Yes |
AKT | Loss of HSC attractors | Loss of quiescence | (9) | Yes |
H2O2 | Loss of HSC attractors | Enhanced differentiation | (88) | Yes |
Antioxidants | HSC basin of attraction is augmented | Enhanced quiescence | (1) | Yes |
OXPHOS | HSC basin of attraction is reduced | Differentiation inhibited | (9) | Yes |
Oxygen | HSC attractors are lost | Loss of quiescence | (175) | Yes |
* This refers to the congruence that exists between the simulations vs. the experimental data. |
Interestingly, our results allow us get insights on controversial aspects about the role of some transcriptional factors in the differentiation of HSCs. Specifically, it has been reported that HIF-1 is not relevant to maintain the quiescence of HSCs (76). However, our data show that its absence (Table 2) abrogates the quiescent HSC phenotype. Therefore, our results are consistent with studies that highlight the importance of this transcriptional factor for quiescent HSCs (77). These data suggest that our model is useful for unraveling controversies about molecular aspects of HSC differentiation.
RNA-seq datasets prove the accuracy of GRN model
To further analyze de accuracy our Boolean model we evaluated its ability to reproduce the gene expression patterns observed in the cell types studied. To this end, we sought RNA-seq data on WT gene expression of the HSC, MEP, GMP, and CLP phenotypes. As a result of the search, we only found data for the HSC (accession number: GSE213596) (52), MEP (accession number: GSE112692) (53) and GMP (accession number: GSE88982) (54) phenotypes. Subsequently, the data were processed and discretized (see Methods) in order to compare them directly with the attractors produced by the network. It should be noted that, to determine if there was OXPHOS activity, we used the Atp6v0b gene as a marker, since it encodes a portion of the V0 domain of vacuolar ATPase (78), and it has been reported that its expression increases proportionally with OXPHOS activity (79). At the end of this procedure, we found that the gene expression patterns predicted by the Boolean model under normal conditions strongly coincide with the real patterns obtained from RNA-seq data (Fig. 4). Taken together, the mutant simulations and the reproduction of the RNA-seq data, demonstrates the effectiveness of our model as a tool to analyze the differentiation of HSCs.
Frequency of Common progenitors emerge from state transitions.
Under homeostatic conditions, hematopoiesis is a process that is tightly regulated, and cell populations maintain specific proportions over time (80). Hence, we aimed to understand how hematopoietic progenitor populations emerge within the niche of HSCs. To do this, we initially employed the size of the basins of attraction in the discrete model as an indicator for the magnitude of phenotypic populations. (Fig. 5a). The basin of attraction is defined as the state space of initial conditions that converge into the same stable state.
In qualitative terms, this approach predicts, from smallest to largest, the following order in the size of phenotypic populations: HSC, GMP, CLP and MEP. Quantitatively, the Boolean model shows that, less than 1% of initial states leads to HSC steady states. In contrast, MEP, GMP, and CLP represent 84.36%, 4.29%, and 11.17% respectively (Fig. 5a). Although promising, these results do not accurately reflect the quantitative nature of the experimental observations (80)
Studies on gene expression show that this process is strongly influenced by randomness (26,81–83), which made us to hypothesize that the proportions observed in vivo (Fig. 5d) might be originated by stochastic fluctuations. Therefore, we applied a stochastic approach to our model, by using software MaBoSS to generate Markov chains (50) (Supplementary Material S5). To obtain such expressions, stochastic fluctuations are applied to the logic rules of the Boolean model (Table 1) so that, at a given instant of time, the model may or may not give the result mathematically expected (50). As a result of this procedure, it was possible to determine the probabilities that the phenotypes have of transitioning between each other. These transition probabilities were used to calculate the stationary probabilities of the Markov chains, which is equivalent to directly calculating the population proportion corresponding to each phenotype (82,83). Since our model considered oxygen to be a signaling molecule from the microenvironment and is not regulated by any other element of the network, we performed this analysis both in the presence and absence of oxygen. In the presence of oxygen, the stochastic simulation recovered a frequency of 40.7% for MEP and 28.5%, 27.5 % and 3.3 % for CLP, GMP and HSC (Fig. 5b). These requencie are closer to the ones observed in vivo where MEP, GMP and CLP represent 35 %, 30% and 9% of CD34 + cells (80) (Fig. 5d). However, ontrary to observed in vivo, in our simulation CLP frequency is greater than GMP. Interestingly, in the absence of oxygen the GMP frequency (27.6 %) is greater than the CLP frequency (26.3 %) while HSC nd MEP show a frequency of 6.4 % and 39.7 respectively (Fig. 5c). Taken together ths data sugest that the population proportion of every cell type seen is defined not only by the constraints of the GRN (Fig. 5a) but also by state transitions between the common progenitors allowed by random genetic drifts (Fig. 5b) and by the hematopoietic stem cell niche´s microenvironment (Fig. 5c).
HSC differentiation is a stochastic process
To further investigate the impact of stochastic drift in hematopoietic stem cell differentiation we simulate the differentiation trajectory of a HSC cell into the common progenitors MEP, GMP and CLP. In order to do this, we initialize our simulation at the HSC stable state and allowed stochastic fluctuations of HSC active and inactive nodes both independently and combined (Fig. 6). For every simulation we track the HSC differentiation trajectory and calculate the entropy and transition entropy of the system (50). In the context of a GRN, maximum entropy is observed when all states have equal probabilities. Conversely, when a particular state has a probability of one, the entropy is reduced to zero (50) .Moreover, transition entropy measures the disorder of a system at single trajectory (50). When the transition entropy is higher than zero the system exhibits a non-deterministic behavior (50).
In the absence of stochastic fluctuations an HSC will not differentiate over time and both entropy and transition entropy are zero (Fig. 6a and Fig. 6e). However, when inactive nodes are allowed a 20% chance of stochastic fluctuation the HSC will differentiate into CLP, MEP and GMP progenitors (Fig. 6b). Transition entropy of this simulation is higher at be beginning converging to zero over time, indicating that the cell fate decision is not entirely deterministic (Fig. 6f).
Interestingly, HSC with 20% stochastic fluctuations of active nodes only differentiate into CLP and MEP progenitors and not GMP are produced (Fig. 6c). Entropy shows an increase during the early differentiation in accordance with recent studies that describe a peak of entropy at the onset of HSC differentiation (Fig. 6g) (84,85). However, stochastic fluctuation of HSC active nodes does not affect transition entropy. The lack of GMP progenitors from this simulation suggest that at least a basal expression of PU1 is required to reach the GMP cell fate, in agreement with PU1 the reported dose-dependent function (Mak et al., 2011).
When both, active and inactive nodes, are allowed to fluctuate stochastically the HSC state is able to differentiate into the MEP, CLP and GMP progenitors (Fig. 6d). Additionally, at the beginning of the simulation the transition entropy is higher than zero (Fig. 6h), thus suggesting that HSC differentiation is a stochastic process, and that basal genetic fluctuation are required for HSC to differentiate into all common progenitors.
Oxygen levels control differentiation of HSCs
Finally, we focused on evaluating the quantitative effect of variations of input signals such as oxygen on the dynamics of the GRN. In order to perform this, we used the stochastic model to explore the effect of gradual variations in oxygen levels. We run stochastic simulation using three different probability distributions of Oxygen as initial states, p[(Oxygen = 1)] = 0, p[(Oxygen = 1)] = 0.5 and, p[(Oxygen = 1)] = 1. Since in our network oxygen is no regulated by any other node, once we establish a probability distribution at the beginning of the simulation it remains constant throughout it.
When oxygen concentration was set to 0% (p[(Oxygen = 1)] = 0), we observed the same phenotypic frequencies mentioned before (Fig. 5c and Fig. 7a). This mimics the microenvironment found in the hematopoietic niche, where low levels of oxygen are necessary to maintain the HSCs undifferentiated and quiescent. However, when oxygen concentration was raised to 50% (p[(Oxygen = 1)] = 0.5), HSCs steady state showed a sharp decrease (Fig. 7b) that was further exacerbated at oxygen 100% (p[(Oxygen = 1)] = 1) (Fig. 7c). Additionally, at higher oxygen concentration our simulation recovers a higher MEP frequency.
These data are in agreement with the experimental evidence that shows that higher oxygen concentration is associated with loss of quiescence and an increment in the proportion of megakaryocytes (6,87,88). Interestingly, along with MEP, CLP frequency also increase at higher oxygen concentration. It has been shown that high oxygen tension is associated with T cell differentiation (89). However, whether oxygen regulates early lymphopoieses remains as an open question. Moreover, HSC and CLP reside within different bone marrow compartments, HSC reside in perivascular niches while CLP are localized at endosteal niches (90). Direct measurements of oxygen tension inside the bone barrow shows that endosteal niches exhibits higher oxygen tension than perivascular ones (91). Taken together, our results predict a potential link of oxygen tension and early lymphoid differentiation.